[Insight-users] c++ code for getting the inverse displacement field form Demon's registration

[Tim Bhatnagar]

Thanks Brian,

I think I did not do enough research into the stipulations on using each
type of deformable registration method. I'll be looking into the
diffeomorphic demons method
 To see if I can obtain more useable results.

Thanks again,

Tim

On Friday, January 3, 2014, brian avants wrote:

> hi tim
>
> there are 2 significant issues with this question, neither of which have
> to do with code
>
> 1) demons deformation fields may not be invertible - this is why
> diffeomorphic demons was developed
>
> 2) your approach to estimating an inverse is only accurate if deformations
> are sub-pixel i.e. motion is less than 1 voxel
>
> one method for estimating an inverse is illustrated here:
>
>
> ./ITKv4/Modules/Filtering/DisplacementField/test/itkInvertDisplacementFieldImageFilterTest.cxx
>
>
> brian
>
>
>
>
> On Thu, Jan 2, 2014 at 4:32 PM, Tim Bhatnagar <tim.bhatnagar at gmail.com<javascript:_e({}, 'cvml', 'tim.bhatnagar at gmail.com');>
> > wrote:
>
>> Hello all,
>>
>> Would anyone have a WORKING c++/source file that successfully
>> generates/acquires the inverse displacement field from a Demon's non-rigid
>> registration field?
>>
>> Many users have been helping me in achieving this, though I'm worried
>> that I'm losing some of the finer points in my code implementation.. this
>> is what I have so far, and it compiles, but the program always crashes at
>> the same point, after a negative index is created for my new field...
>>
>> any help is much appreciated.. my c++ code is below:
>>
>> #include "itkInverseDisplacementFieldImageFilter.h"
>> #include "itkImageFileWriter.h"
>> #include "itkImageFileReader.h"
>> #include "itkFilterWatcher.h"
>> #include "itkVectorScaleImageFilter.h"
>> #include "itkImage.h"
>> #include "itkArray.h"
>>
>> int main( int argc, char * argv[] )
>> {
>>
>>   if( argc < 4 )
>>     {
>>     std::cerr << "Missing Parameters " << std::endl;
>>     std::cerr << "Usage: " << argv[0];
>>     std::cerr << "inputDemonsField subSampFactor outputInvertedField" <<
>> std::endl;
>>     return EXIT_FAILURE;
>>     }
>>
>>   const     unsigned int Dimension = 3;
>>   typedef   float VectorComponentType;
>>
>>   typedef   itk::Vector< VectorComponentType, Dimension > VectorType;
>>
>>   typedef itk::Image< VectorType,  Dimension > DisplacementFieldType;
>>
>>   typedef itk::Array< VectorType > ArrayType;
>>
>>   typedef itk::InverseDisplacementFieldImageFilter<
>>     DisplacementFieldType,
>>     DisplacementFieldType
>>     >  FilterType;
>>
>>   FilterType::Pointer filter = FilterType::New();
>>
>>   FilterWatcher watcher(filter);
>>
>>   // Read a DisplacementField
>>   typedef itk::ImageFileReader<  DisplacementFieldType  > ReaderType;
>>
>>   ReaderType::Pointer reader = ReaderType::New();
>>
>>   reader->SetFileName( argv[1] );
>>
>>   reader->Update();
>>
>>
>>   // Obtaining an input displacement field
>>   DisplacementFieldType::Pointer field = reader->GetOutput();
>>
>>   DisplacementFieldType::RegionType inputRegion;
>>   DisplacementFieldType::SizeType inputSize;
>>   DisplacementFieldType::DirectionType inputDirection;
>>   DisplacementFieldType::PointType inputOrigin;
>>   DisplacementFieldType::SpacingType inputSpacing;
>>
>>   inputRegion = field->GetLargestPossibleRegion();
>>   inputSize = field->GetLargestPossibleRegion().GetSize();
>>   inputDirection =  field->GetDirection();
>>   inputOrigin =  field->GetOrigin();
>>   inputSpacing =  field->GetSpacing();
>>   inputRegion.SetIndex(field->GetLargestPossibleRegion().GetIndex());
>>
>>   DisplacementFieldType::Pointer outField = DisplacementFieldType::New();
>>
>>
>>
>>
>>    itk::ImageRegionIteratorWithIndex< DisplacementFieldType > it( field,
>> inputRegion );
>>
>>    filter->SetOutputSpacing( field->GetSpacing() );
>>    filter->SetOutputOrigin( field->GetOrigin() );
>>    filter->SetSize( inputSize );
>>
>>   filter->SetInput( field );
>>   filter->SetSubsamplingFactor( atof(argv[2]) );
>>
>>   try
>>     {
>>     filter->UpdateLargestPossibleRegion();
>>     }
>>   catch( itk::ExceptionObject & excp )
>>     {
>>     std::cerr << "Exception thrown " << std::endl;
>>     std::cerr << excp << std::endl;
>>     }
>>
>>   // Make template to ensure proper indexing of new def field
>>
>>   DisplacementFieldType::Pointer outField = filter->GetOutput();
>>   DisplacementFieldType::RegionType outputRegion;
>>
>>    outputRegion.SetIndex( inputRegion.GetIndex() );
>>
>>   // Write an image for regression testing
>>   typedef itk::ImageFileWriter<  DisplacementFieldType  > WriterType;
>>
>>   WriterType::Pointer writer = WriterType::New();
>>
>>   outField->Update();
>>   writer->SetInput (outField );
>>   writer->SetFileName( argv[3] );
>>
>>   try
>>     {
>>     writer->Update();
>>     }
>>   catch( itk::ExceptionObject & excp )
>>     {
>>     std::cerr << "Exception thrown by writer" << std::endl;
>>     std::cerr << excp << std::endl;
>>     return EXIT_FAILURE;
>>     }
>>
>>
>>
>>
>>   // Now, test for loop invariant (acts as filter validation)
>>   // f^-1(f(p1) + p1 ) - f(p1)  = 0
>>   it.GoToBegin();
>>   while( !it.IsAtEnd() )
>>     {
>>
>>     DisplacementFieldType::PointType p1;
>>     field->TransformIndexToPhysicalPoint(it.GetIndex(), p1);
>>
>>     DisplacementFieldType::PixelType fp1 = it.Get();
>>
>>     DisplacementFieldType::PointType p2;
>>     p2[0] = p1[0] + fp1[0];
>>     p2[1] = p1[1] + fp1[1];
>>     p2[2] = p1[2] + fp1[2];
>>
>>     DisplacementFieldType::IndexType id2;
>>     outField->Transform
>> _____________________________________
>> Powered by www.kitware.com
>>
>> Visit other Kitware open-source projects at
>> http://www.kitware.com/opensource/opensource.html
>>
>> Kitware offers ITK Training Courses, for more information visit:
>> http://www.kitware.com/products/protraining.php
>>
>> Please keep messages on-topic and check the ITK FAQ at:
>> http://www.itk.org/Wiki/ITK_FAQ
>>
>> Follow this link to subscribe/unsubscribe:
>> http://www.itk.org/mailman/listinfo/insight-users
>>
>>
>

-- 
Tim Bhatnagar
PhD Candidate
Orthopaedic Injury Biomechanics Group
Department of Mechanical Engineering
University of British Columbia

Rm 5000 - 818 West 10th Ave.
Vancouver, BC
Canada
V5Z 1M9

Ph: (604) 675-8845
Fax: (604) 675-8820
Web: oibg.mech.ubc.ca
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://www.itk.org/pipermail/insight-users/attachments/20140105/1212c6ba/attachment.htm>