[Insight-users] Deformation field is not applied correctly when using WarpImageFilter

[Steve Lancey]

Hi all,

I have been using itk for quite some time, mainly the FEM package. I run FEM
registration on 3D Volumes with around 100 iteration. As a result of the
registration procedure, the warped image file and the multi component
deformation field are stored.
When visualizing the warped image + deformation field in Paraview it appears
to be correct.

The problem I am facing is the following:

I want to manually warp the moving image to obtain the registered image
(moving + deformation field = registered image)
I am using the WarpImageFilter, more or less the same code that is used in
the FEMRegistrationFilter. I am unable to reproduce the same result, it
seems the image has been barely 'warped' when visualizing the result. There
is a big difference between the warped image obtained through the
FEMRegistrationFilter and by manually warping it.

Is my assumption correct that the Deformation field returned from the
FEMRegistrationFilter holds the total deformation and not only the
deformation of the last iteration.

This is the code I am working with (it is not optimized)


#if defined(_MSC_VER)
#pragma warning ( disable : 4786 )
#endif

#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"
#include "itkMinimumMaximumImageFilter.h"
#include "itkShiftScaleImageFilter.h"
#include "itkImageRegionIterator.h"


#include "itkCastImageFilter.h"
#include "itkWarpImageFilter.h"
#include "itkLinearInterpolateImageFunction.h"

#include "itkSquaredDifferenceImageFilter.h"
#include "itkCheckerBoardImageFilter.h"
#include "itkTimeProbe.h"

#include "itkImageRegionIteratorWithIndex.h"
#include "itkPointSet.h"
#include <fstream>

#include "itkCastImageFilter.h"

int main( int argc, char *argv[] )
{
    if( argc < 4 )
    {
        std::cerr << "Missing Parameters " << std::endl;
        std::cerr << "Usage: " << argv[0];
        std::cerr << " fixedImagefile movingImageFile deformField.vtk ";
        std::cerr << " outputImageFile " << std::endl;
        return 1;
    }


    const unsigned int Dimension = 3;
    typedef float PixelType;


    typedef itk::Image< PixelType, Dimension >  FixedImageType;
    typedef itk::ImageFileReader< FixedImageType > FixedImageReaderType;
    FixedImageReaderType::Pointer fixedImageReader =
FixedImageReaderType::New();
    fixedImageReader->SetFileName( argv[1] );


    typedef itk::Image< PixelType, Dimension >  MovingImageType;
    typedef itk::ImageFileReader< MovingImageType > MovingImageReaderType;
    MovingImageReaderType::Pointer movingImageReader =
MovingImageReaderType::New();
    movingImageReader->SetFileName( argv[2] );

    typedef itk::Vector< float, Dimension >    VectorPixelType;
    typedef itk::Image<  VectorPixelType, Dimension > DeformationFieldType;

    typedef itk::ImageFileReader< DeformationFieldType >  FieldReaderType;
    FieldReaderType::Pointer fieldReader = FieldReaderType::New();

    fieldReader->SetFileName( argv[3] );


   try {

    fieldReader->Update();
movingImageReader->Update();
fixedImageReader->Update();
   }

   catch ( itk::ExceptionObject & err )
   {
    std::cerr << "ExceptionObject caught !" << std::endl;
    std::cerr << err << std::endl;
    return -1;
    }



    typedef itk::WarpImageFilter<FixedImageType,
MovingImageType,DeformationFieldType  >     WarperType;
    WarperType::Pointer warper = WarperType::New();
 typedef WarperType::CoordRepType WarperCoordRepType;
typedef itk::LinearInterpolateImageFunction<MovingImageType,
WarperCoordRepType  >  InterpolatorType;

    InterpolatorType::Pointer interpolator = InterpolatorType::New();
    FixedImageType::Pointer fixedImage = fixedImageReader->GetOutput();
FixedImageType::Pointer warpedImage;

  warper->SetInput( movingImageReader->GetOutput() );
      warper->SetDeformationField( fieldReader->GetOutput() );
      warper->SetInterpolator( interpolator );
      warper->SetOutputOrigin( fixedImage->GetOrigin() );
      warper->SetOutputSpacing( fixedImage->GetSpacing() );
      warper->SetOutputDirection( fixedImage->GetDirection() );
  MovingImageType::PixelType padValue = 0;
      warper->SetEdgePaddingValue( padValue );
     warper->Update();
 warpedImage = warper->GetOutput();

    /*typedef itk::ImageFileWriter< MovingImageType >  WriterType;

    WriterType::Pointer writer =  WriterType::New();
    writer->SetFileName( argv[4] );
    writer->SetInput( warper->GetOutput()   );
    writer->Update();*/



  typedef itk::MinimumMaximumImageFilter<MovingImageType> MinMaxFilterType;
  MinMaxFilterType::Pointer minMaxFilter = MinMaxFilterType::New();
  minMaxFilter->SetInput( warpedImage );
  minMaxFilter->Update();
  float min = minMaxFilter->GetMinimum();
  double shift = -1.0 * static_cast<double>( min );
  double scale = static_cast<double>( minMaxFilter->GetMaximum() );
  scale += shift;
  scale = 255.0 / scale;
  typedef itk::ShiftScaleImageFilter<MovingImageType, MovingImageType>
FilterType;
  FilterType::Pointer filter = FilterType::New();
  filter->SetInput( warpedImage );
  filter->SetShift( shift );
  filter->SetScale( scale );
  filter->Update();


  typedef itk::Image<PixelType,Dimension>                    WriteImageType;
  typedef itk::CastImageFilter<MovingImageType,WriteImageType> CasterType1;
  CasterType1::Pointer caster1 = CasterType1::New();
  caster1->SetInput(filter->GetOutput());
  caster1->Update();
  typedef itk::ImageFileWriter<WriteImageType> WriterType;
  WriterType::Pointer writer = WriterType::New();
  writer->SetFileName(argv[4]);
  writer->SetInput(caster1->GetOutput() );
  writer->Write();
    return 0;
}
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://www.itk.org/pipermail/insight-users/attachments/20100310/289a3471/attachment.htm>