<div dir="ltr">hi tim <div><br></div><div>there are 2 significant issues with this question, neither of which have to do with code</div><div><br></div><div>1) demons deformation fields may not be invertible - this is why diffeomorphic demons was developed</div>
<div><br></div><div>2) your approach to estimating an inverse is only accurate if deformations are sub-pixel i.e. motion is less than 1 voxel</div><div><br></div><div>one method for estimating an inverse is illustrated here:</div>
<div>
<p class="">./ITKv4/Modules/Filtering/DisplacementField/test/itkInvertDisplacementFieldImageFilterTest.cxx</p></div></div><div class="gmail_extra"><br clear="all"><div><div><br></div>brian<br><div><br></div><div><br></div>
</div>
<br><br><div class="gmail_quote">On Thu, Jan 2, 2014 at 4:32 PM, Tim Bhatnagar <span dir="ltr"><<a href="mailto:tim.bhatnagar@gmail.com" target="_blank">tim.bhatnagar@gmail.com</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
<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: <a href="tel:%28604%29%20675-8845" value="+16046758845" target="_blank">(604) 675-8845</a><br>Fax: <a href="tel:%28604%29%20675-8820" value="+16046758820" target="_blank">(604) 675-8820</a><br>
Web: <a href="http://oibg.mech.ubc.ca" target="_blank">oibg.mech.ubc.ca</a><br>
</div></div></div></div>
<br>_____________________________________<br>
Powered by <a href="http://www.kitware.com" target="_blank">www.kitware.com</a><br>
<br>
Visit other Kitware open-source projects at<br>
<a href="http://www.kitware.com/opensource/opensource.html" target="_blank">http://www.kitware.com/opensource/opensource.html</a><br>
<br>
Kitware offers ITK Training Courses, for more information visit:<br>
<a href="http://www.kitware.com/products/protraining.php" target="_blank">http://www.kitware.com/products/protraining.php</a><br>
<br>
Please keep messages on-topic and check the ITK FAQ at:<br>
<a href="http://www.itk.org/Wiki/ITK_FAQ" target="_blank">http://www.itk.org/Wiki/ITK_FAQ</a><br>
<br>
Follow this link to subscribe/unsubscribe:<br>
<a href="http://www.itk.org/mailman/listinfo/insight-users" target="_blank">http://www.itk.org/mailman/listinfo/insight-users</a><br>
<br></blockquote></div><br></div>