Hi all,<div><br></div><div>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.</div>
<div>When visualizing the warped image + deformation field in Paraview it appears to be correct.</div><div><br></div><div>The problem I am facing is the following:</div><div><br></div><div>I want to manually warp the moving image to obtain the registered image (moving + deformation field = registered image)</div>
<div>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.</div>
<div><br></div><div>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.</div><div><br></div><div>This is the code I am working with (it is not optimized)</div>
<div><br></div><div><br></div><div><div>#if defined(_MSC_VER)</div><div>#pragma warning ( disable : 4786 )</div><div>#endif</div><div><br></div><div>#include "itkImageFileReader.h" </div><div>#include "itkImageFileWriter.h" </div>
<div>#include "itkMinimumMaximumImageFilter.h"</div><div>#include "itkShiftScaleImageFilter.h"</div><div>#include "itkImageRegionIterator.h"</div><div><br></div><div><br></div><div>#include "itkCastImageFilter.h"</div>
<div>#include "itkWarpImageFilter.h"</div><div>#include "itkLinearInterpolateImageFunction.h"</div><div><br></div><div>#include "itkSquaredDifferenceImageFilter.h"</div><div>#include "itkCheckerBoardImageFilter.h"</div>
<div>#include "itkTimeProbe.h" </div><div><br></div><div>#include "itkImageRegionIteratorWithIndex.h"</div><div>#include "itkPointSet.h"</div><div>#include <fstream></div><div><br>
</div><div>#include "itkCastImageFilter.h"</div><div><br></div><div>int main( int argc, char *argv[] )</div><div>{</div><div> if( argc < 4 )</div><div> {</div><div> std::cerr << "Missing Parameters " << std::endl;</div>
<div> std::cerr << "Usage: " << argv[0];</div><div> std::cerr << " fixedImagefile movingImageFile deformField.vtk ";</div><div> std::cerr << " outputImageFile " << std::endl;</div>
<div> return 1;</div><div> }</div><div><br></div><div> </div><div> const unsigned int Dimension = 3;</div><div> typedef float PixelType;</div><div><br></div><div><br></div><div> typedef itk::Image< PixelType, Dimension > FixedImageType;</div>
<div> typedef itk::ImageFileReader< FixedImageType > FixedImageReaderType;</div><div> FixedImageReaderType::Pointer fixedImageReader = FixedImageReaderType::New();</div><div> fixedImageReader->SetFileName( argv[1] );</div>
<div><br></div><div> </div><div> typedef itk::Image< PixelType, Dimension > MovingImageType;</div><div> typedef itk::ImageFileReader< MovingImageType > MovingImageReaderType;</div><div> MovingImageReaderType::Pointer movingImageReader = MovingImageReaderType::New();</div>
<div> movingImageReader->SetFileName( argv[2] );</div><div><br></div><div> typedef itk::Vector< float, Dimension > VectorPixelType;</div><div> typedef itk::Image< VectorPixelType, Dimension > DeformationFieldType;</div>
<div><br></div><div> typedef itk::ImageFileReader< DeformationFieldType > FieldReaderType;</div><div> FieldReaderType::Pointer fieldReader = FieldReaderType::New();</div><div><br></div><div> fieldReader->SetFileName( argv[3] );</div>
<div><br></div><div> </div><div> try {</div><div><br></div><div> fieldReader->Update();</div><div><span class="Apple-tab-span" style="white-space:pre"> </span>movingImageReader->Update();</div><div><span class="Apple-tab-span" style="white-space:pre"> </span>fixedImageReader->Update();</div>
<div> }</div><div><br></div><div> catch ( itk::ExceptionObject & err ) </div><div> {</div><div> std::cerr << "ExceptionObject caught !" << std::endl; </div><div> std::cerr << err << std::endl; </div>
<div> return -1;</div><div> } </div><div><br></div><div><br></div><div> </div><div> typedef itk::WarpImageFilter<FixedImageType, MovingImageType,DeformationFieldType > WarperType;</div><div> WarperType::Pointer warper = WarperType::New();</div>
<div><span class="Apple-tab-span" style="white-space:pre"> </span></div><div><span class="Apple-tab-span" style="white-space:pre"> </span>typedef WarperType::CoordRepType WarperCoordRepType;</div><div><span class="Apple-tab-span" style="white-space:pre"> </span>typedef itk::LinearInterpolateImageFunction<MovingImageType, WarperCoordRepType > InterpolatorType;</div>
<div><br></div><div> InterpolatorType::Pointer interpolator = InterpolatorType::New();</div><div> FixedImageType::Pointer fixedImage = fixedImageReader->GetOutput();</div><div><span class="Apple-tab-span" style="white-space:pre"> </span>FixedImageType::Pointer warpedImage;</div>
<div><br></div><div><span class="Apple-tab-span" style="white-space:pre"> </span> warper->SetInput( movingImageReader->GetOutput() );</div><div> warper->SetDeformationField( fieldReader->GetOutput() );</div>
<div> warper->SetInterpolator( interpolator );</div><div> warper->SetOutputOrigin( fixedImage->GetOrigin() );</div><div> warper->SetOutputSpacing( fixedImage->GetSpacing() );</div><div> warper->SetOutputDirection( fixedImage->GetDirection() );</div>
<div><span class="Apple-tab-span" style="white-space:pre"> </span> MovingImageType::PixelType padValue = 0;</div><div> warper->SetEdgePaddingValue( padValue );</div><div><span class="Apple-tab-span" style="white-space:pre"> </span></div>
<div> warper->Update();</div><div><span class="Apple-tab-span" style="white-space:pre"> </span> warpedImage = warper->GetOutput();</div><div><br></div><div> /*typedef itk::ImageFileWriter< MovingImageType > WriterType;</div>
<div> </div><div> WriterType::Pointer writer = WriterType::New();</div><div> writer->SetFileName( argv[4] );</div><div> writer->SetInput( warper->GetOutput() );</div><div> writer->Update();*/</div>
<div><br></div><div><br></div><div><br></div><div> typedef itk::MinimumMaximumImageFilter<MovingImageType> MinMaxFilterType;</div><div> MinMaxFilterType::Pointer minMaxFilter = MinMaxFilterType::New();</div><div> minMaxFilter->SetInput( warpedImage );</div>
<div> minMaxFilter->Update();</div><div> float min = minMaxFilter->GetMinimum();</div><div> double shift = -1.0 * static_cast<double>( min );</div><div> double scale = static_cast<double>( minMaxFilter->GetMaximum() );</div>
<div> scale += shift;</div><div> scale = 255.0 / scale;</div><div> typedef itk::ShiftScaleImageFilter<MovingImageType, MovingImageType> FilterType;</div><div> FilterType::Pointer filter = FilterType::New();</div>
<div> filter->SetInput( warpedImage );</div><div> filter->SetShift( shift );</div><div> filter->SetScale( scale );</div><div> filter->Update();</div><div><br></div><div><br></div><div> typedef itk::Image<PixelType,Dimension> WriteImageType;</div>
<div> typedef itk::CastImageFilter<MovingImageType,WriteImageType> CasterType1;</div><div> CasterType1::Pointer caster1 = CasterType1::New();</div><div> caster1->SetInput(filter->GetOutput());</div><div> caster1->Update();</div>
<div> typedef itk::ImageFileWriter<WriteImageType> WriterType;</div><div> WriterType::Pointer writer = WriterType::New();</div><div> writer->SetFileName(argv[4]);</div><div> writer->SetInput(caster1->GetOutput() );</div>
<div> writer->Write();</div><div><span class="Apple-tab-span" style="white-space:pre"> </span></div><div> return 0;</div><div>}</div></div><div><br></div><div><br></div><div><br></div>