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

[Tim Bhatnagar]

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
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://www.itk.org/pipermail/insight-users/attachments/20140102/e364ed55/attachment.htm>