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

[Luis Ibanez]

Hi Steve,

Thanks for your detailed description of the problem.


In the code that you posted, we see that you are
reading the deformation field from a file and passing
it to the WrapImageFilter:

     warper->SetDeformationField( fieldReader->GetOutput() );


Could you please describe how you generated that file ?


Did you save in there the vector image returned by the
FEMRegistrationFilter via the method

             FieldType* GetDeformationField()   ?




Also the warped image that you are using as reference
for comparison... did you generated this image by
calling the method :

  warpedImage =
    femRegistration->WarpImage( movingImage )   ?



Please let us know,


     Thanks


          Luis


--------------------------------------------------------------------------
On Wed, Mar 10, 2010 at 11:38 PM, Steve Lancey <steve.lancey at gmail.com> wrote:
> 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;
> }
>
>
>
> _____________________________________
> 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.html
>
> 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
>
>