[ITK Community] [Insight-users] c++ code for getting the inverse displacement field form Demon's registration
brian avants
stnava at gmail.com
Fri Jan 3 09:10:19 EST 2014
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://public.kitware.com/pipermail/community/attachments/20140103/62e4228d/attachment.html>
-------------- next part --------------
_____________________________________
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
More information about the Community
mailing list