<div dir="ltr"><div><div>Hello all,<br><br></div>Would anyone have a WORKING c++/source file that successfully generates/acquires the inverse displacement field from a Demon's non-rigid registration field?<br><br></div>
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...<br>
<br>any help is much appreciated.. my c++ code is below:<br><br>#include "itkInverseDisplacementFieldImageFilter.h"<br>#include "itkImageFileWriter.h"<br>#include "itkImageFileReader.h"<br>#include "itkFilterWatcher.h"<br>
#include "itkVectorScaleImageFilter.h"<br>#include "itkImage.h"<br>#include "itkArray.h"<br><br>int main( int argc, char * argv[] )<br>{<br><br> if( argc < 4 )<br> {<br> std::cerr << "Missing Parameters " << std::endl;<br>
std::cerr << "Usage: " << argv[0];<br> std::cerr << "inputDemonsField subSampFactor outputInvertedField" << std::endl;<br> return EXIT_FAILURE;<br> }<br><br> const unsigned int Dimension = 3;<br>
typedef float VectorComponentType;<br><br> typedef itk::Vector< VectorComponentType, Dimension > VectorType;<br><br> typedef itk::Image< VectorType, Dimension > DisplacementFieldType;<br><br> typedef itk::Array< VectorType > ArrayType;<br>
<br> typedef itk::InverseDisplacementFieldImageFilter<<br> DisplacementFieldType,<br> DisplacementFieldType<br> > FilterType;<br><br> FilterType::Pointer filter = FilterType::New();<br><br> FilterWatcher watcher(filter);<br>
<br> // Read a DisplacementField<br> typedef itk::ImageFileReader< DisplacementFieldType > ReaderType;<br><br> ReaderType::Pointer reader = ReaderType::New();<br><br> reader->SetFileName( argv[1] );<br><br>
reader->Update();<br><br><br> // Obtaining an input displacement field<br> DisplacementFieldType::Pointer field = reader->GetOutput();<br><br> DisplacementFieldType::RegionType inputRegion;<br> DisplacementFieldType::SizeType inputSize;<br>
DisplacementFieldType::DirectionType inputDirection;<br> DisplacementFieldType::PointType inputOrigin;<br> DisplacementFieldType::SpacingType inputSpacing;<br><br> inputRegion = field->GetLargestPossibleRegion();<br>
inputSize = field->GetLargestPossibleRegion().GetSize();<br> inputDirection = field->GetDirection();<br> inputOrigin = field->GetOrigin();<br> inputSpacing = field->GetSpacing();<br> inputRegion.SetIndex(field->GetLargestPossibleRegion().GetIndex());<br>
<br> DisplacementFieldType::Pointer outField = DisplacementFieldType::New();<br><br> <br><br><br> itk::ImageRegionIteratorWithIndex< DisplacementFieldType > it( field, inputRegion );<br><br> filter->SetOutputSpacing( field->GetSpacing() );<br>
filter->SetOutputOrigin( field->GetOrigin() );<br> filter->SetSize( inputSize );<br><br> filter->SetInput( field );<br> filter->SetSubsamplingFactor( atof(argv[2]) );<br><br> try<br> {<br> filter->UpdateLargestPossibleRegion();<br>
}<br> catch( itk::ExceptionObject & excp )<br> {<br> std::cerr << "Exception thrown " << std::endl;<br> std::cerr << excp << std::endl;<br> }<br><br> // Make template to ensure proper indexing of new def field<br>
<br> DisplacementFieldType::Pointer outField = filter->GetOutput();<br> DisplacementFieldType::RegionType outputRegion;<br><br> outputRegion.SetIndex( inputRegion.GetIndex() );<br><br> // Write an image for regression testing<br>
typedef itk::ImageFileWriter< DisplacementFieldType > WriterType;<br><br> WriterType::Pointer writer = WriterType::New();<br><br> outField->Update();<br> writer->SetInput (outField );<br> writer->SetFileName( argv[3] );<br>
<br> try<br> {<br> writer->Update();<br> }<br> catch( itk::ExceptionObject & excp )<br> {<br> std::cerr << "Exception thrown by writer" << std::endl;<br> std::cerr << excp << std::endl;<br>
return EXIT_FAILURE;<br> }<br><br><br><br><br> // Now, test for loop invariant (acts as filter validation)<br> // f^-1(f(p1) + p1 ) - f(p1) = 0<br> it.GoToBegin();<br> while( !it.IsAtEnd() )<br> {<br><br> DisplacementFieldType::PointType p1;<br>
field->TransformIndexToPhysicalPoint(it.GetIndex(), p1);<br><br> DisplacementFieldType::PixelType fp1 = it.Get();<br><br> DisplacementFieldType::PointType p2;<br> p2[0] = p1[0] + fp1[0];<br> p2[1] = p1[1] + fp1[1];<br>
p2[2] = p1[2] + fp1[2];<br><br> DisplacementFieldType::IndexType id2;<br> outField->TransformPhysicalPointToIndex(p2,id2);<br> //filter->GetOutput()->TransformPhysicalPointToIndex(p2,id2);<br> std::cerr << "hurrah" << id2 << std::endl;<br>
DisplacementFieldType::PixelType fp2 = outField->GetPixel(id2);<br> //DisplacementFieldType::PixelType fp2 = filter->GetOutput()->GetPixel(id2);<br> std::cerr << "hurrah2" << std::endl;<br>
<br><br> if(vcl_abs(fp2[0] + fp1[0]) > 100.000<br> || vcl_abs(fp2[1] + fp1[1]) > 100.000<br> || vcl_abs(fp2[2] + fp1[2]) > 100.000)<br> {<br> std::cerr<<"Loop invariant not satisfied for index "<<it.GetIndex()<<" : f^-1(f(p1) + p1 ) + f(p1) = 0"<< std::endl;<br>
std::cerr<<"f(p1) = "<<fp1<<std::endl;<br> std::cerr<<"f^-1(f(p1) + p1 ) = "<<fp2<<std::endl;<br> std::cerr<<"diff: "<<fp1[0]+fp2[0]<<", "<<fp1[1]+fp2[1]<<", "<<fp1[2]+fp2[2]<<std::endl;<br>
// return EXIT_FAILURE;<br> }<br> ++it;<br> }<br><br> return EXIT_SUCCESS;<br><br>}<br clear="all"><div><div><div><br>-- <br>Tim Bhatnagar<br>PhD Candidate<br>Orthopaedic Injury Biomechanics Group<br>Department of Mechanical Engineering<br>
University of British Columbia<br><br>Rm 5000 - 818 West 10th Ave.<br>Vancouver, BC<br>Canada<br>V5Z 1M9<br><br>Ph: (604) 675-8845<br>Fax: (604) 675-8820<br>Web: <a href="http://oibg.mech.ubc.ca" target="_blank">oibg.mech.ubc.ca</a><br>
</div></div></div></div>