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

[brian avants]

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>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->TransformPhysicalPointToIndex(p2,id2);
>     //filter->GetOutput()->TransformPhysicalPointToIndex(p2,id2);
>     std::cerr << "hurrah" << id2 << std::endl;
>     DisplacementFieldType::PixelType fp2 = outField->GetPixel(id2);
>     //DisplacementFieldType::PixelType fp2 =
> filter->GetOutput()->GetPixel(id2);
>     std::cerr << "hurrah2" << std::endl;
>
>
>     if(vcl_abs(fp2[0] + fp1[0]) > 100.000
>        || vcl_abs(fp2[1] + fp1[1]) > 100.000
>        || vcl_abs(fp2[2] + fp1[2]) > 100.000)
>       {
>       std::cerr<<"Loop invariant not satisfied for index
> "<<it.GetIndex()<<" : f^-1(f(p1) + p1 ) + f(p1)  = 0"<<   std::endl;
>       std::cerr<<"f(p1) = "<<fp1<<std::endl;
>       std::cerr<<"f^-1(f(p1) + p1 ) = "<<fp2<<std::endl;
>       std::cerr<<"diff: "<<fp1[0]+fp2[0]<<", "<<fp1[1]+fp2[1]<<",
> "<<fp1[2]+fp2[2]<<std::endl;
>     //  return EXIT_FAILURE;
>       }
>     ++it;
>     }
>
>   return EXIT_SUCCESS;
>
> }
>
> --
> 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
>
> _____________________________________
> 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
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://www.itk.org/pipermail/insight-users/attachments/20140103/62e4228d/attachment.htm>