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

[Steve Lancey]

Hi Luis,

thanks for your answer:

1)
>>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 ?
I am using the following code:

typedef itk::fem::FEMRegistrationFilter<ImageType,ImageType>
RegistrationType;
...
...
...
  RegistrationType::Pointer registrationFilter = RegistrationType::New();
...
...
...
  registrationFilter->WriteWarpedImage(
        (registrationFilter->GetResultsFileName()).c_str());

2)
>>Did you save in there the vector image returned by the
>>FEMRegistrationFilter via the method
>>            FieldType* GetDeformationField()   ?
...
...
...
registrationFilter->WriteDisplacementFieldMultiComponent();

Please let me know if you need additional information.

Steve




On Tue, Mar 16, 2010 at 6:08 PM, Luis Ibanez <luis.ibanez at kitware.com>wrote:

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