[Insight-users] Deformation field is not applied correctly when using WarpImageFilter
Steve Lancey
steve.lancey at gmail.com
Sat Mar 27 19:48:07 EDT 2010
Hi Luis,
I am still struggling with this, do you have any ideas?
Thanks,
Steve
On Thu, Mar 18, 2010 at 10:37 PM, Steve Lancey <steve.lancey at gmail.com>wrote:
> 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/20100327/a089f40d/attachment-0001.htm>
More information about the Insight-users
mailing list