[Insight-users] Deformation field is not applied correctly when using WarpImageFilter
Luis Ibanez
luis.ibanez at kitware.com
Mon Mar 29 18:07:28 EDT 2010
Hi Steve,
It looks like a bug, but I haven't had a chance to
attempt replicating it.
Could you please file a bug report in the bug
tracker and attach to it the source code that you
are running and that is showing this behavior ?
Thanks
Luis
----------------------------------------------------------------
On Sat, Mar 27, 2010 at 7:48 PM, Steve Lancey <steve.lancey at gmail.com> wrote:
> 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
>>> >
>>> >
>>
>
>
More information about the Insight-users
mailing list