Hi Luis,<div><br></div><div>thanks for your answer:</div><div><br></div><div>1)</div><div>>>In the code that you posted, we see that you are<br>>>reading the deformation field from a file and passing<br>>>it to the WrapImageFilter:<br>
>> warper->SetDeformationField( fieldReader->GetOutput() );<br>>>Could you please describe how you generated that file ?</div><div>I am using the following code:</div><div><br></div><div>typedef itk::fem::FEMRegistrationFilter<ImageType,ImageType> RegistrationType;</div>
<div>...</div><div>...</div><div>...</div><div> RegistrationType::Pointer registrationFilter = RegistrationType::New(); </div><div>...</div><div>...</div><div>...</div><div><div> registrationFilter->WriteWarpedImage(</div>
<div> (registrationFilter->GetResultsFileName()).c_str());</div></div><div><br>2)</div><div>>>Did you save in there the vector image returned by the<br>>>FEMRegistrationFilter via the method<br>>> FieldType* GetDeformationField() ?<br>
</div><div>...</div><div>...</div><div>...</div><div>registrationFilter->WriteDisplacementFieldMultiComponent();</div><div><br></div><div>Please let me know if you need additional information.</div><div><br></div><div>
Steve</div><div><br></div><div><br></div><div><br></div><div><br><div class="gmail_quote">On Tue, Mar 16, 2010 at 6:08 PM, Luis Ibanez <span dir="ltr"><<a href="mailto:luis.ibanez@kitware.com">luis.ibanez@kitware.com</a>></span> wrote:<br>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex;">Hi Steve,<br>
<br>
Thanks for your detailed description of the problem.<br>
<br>
<br>
In the code that you posted, we see that you are<br>
reading the deformation field from a file and passing<br>
it to the WrapImageFilter:<br>
<br>
warper->SetDeformationField( fieldReader->GetOutput() );<br>
<br>
<br>
Could you please describe how you generated that file ?<br>
<br>
<br>
Did you save in there the vector image returned by the<br>
FEMRegistrationFilter via the method<br>
<br>
FieldType* GetDeformationField() ?<br>
<br>
<br>
<br>
<br>
Also the warped image that you are using as reference<br>
for comparison... did you generated this image by<br>
calling the method :<br>
<br>
warpedImage =<br>
femRegistration->WarpImage( movingImage ) ?<br>
<br>
<br>
<br>
Please let us know,<br>
<br>
<br>
Thanks<br>
<br>
<br>
Luis<br>
<br>
<br>
--------------------------------------------------------------------------<br>
<div><div></div><div class="h5">On Wed, Mar 10, 2010 at 11:38 PM, Steve Lancey <<a href="mailto:steve.lancey@gmail.com">steve.lancey@gmail.com</a>> wrote:<br>
> Hi all,<br>
> I have been using itk for quite some time, mainly the FEM package. I run FEM<br>
> registration on 3D Volumes with around 100 iteration. As a result of the<br>
> registration procedure, the warped image file and the multi component<br>
> deformation field are stored.<br>
> When visualizing the warped image + deformation field in Paraview it appears<br>
> to be correct.<br>
> The problem I am facing is the following:<br>
> I want to manually warp the moving image to obtain the registered image<br>
> (moving + deformation field = registered image)<br>
> I am using the WarpImageFilter, more or less the same code that is used in<br>
> the FEMRegistrationFilter. I am unable to reproduce the same result, it<br>
> seems the image has been barely 'warped' when visualizing the result. There<br>
> is a big difference between the warped image obtained through the<br>
> FEMRegistrationFilter and by manually warping it.<br>
> Is my assumption correct that the Deformation field returned from the<br>
> FEMRegistrationFilter holds the total deformation and not only the<br>
> deformation of the last iteration.<br>
> This is the code I am working with (it is not optimized)<br>
><br>
> #if defined(_MSC_VER)<br>
> #pragma warning ( disable : 4786 )<br>
> #endif<br>
> #include "itkImageFileReader.h"<br>
> #include "itkImageFileWriter.h"<br>
> #include "itkMinimumMaximumImageFilter.h"<br>
> #include "itkShiftScaleImageFilter.h"<br>
> #include "itkImageRegionIterator.h"<br>
><br>
> #include "itkCastImageFilter.h"<br>
> #include "itkWarpImageFilter.h"<br>
> #include "itkLinearInterpolateImageFunction.h"<br>
> #include "itkSquaredDifferenceImageFilter.h"<br>
> #include "itkCheckerBoardImageFilter.h"<br>
> #include "itkTimeProbe.h"<br>
> #include "itkImageRegionIteratorWithIndex.h"<br>
> #include "itkPointSet.h"<br>
> #include <fstream><br>
> #include "itkCastImageFilter.h"<br>
> int main( int argc, char *argv[] )<br>
> {<br>
> if( argc < 4 )<br>
> {<br>
> std::cerr << "Missing Parameters " << std::endl;<br>
> std::cerr << "Usage: " << argv[0];<br>
> std::cerr << " fixedImagefile movingImageFile deformField.vtk ";<br>
> std::cerr << " outputImageFile " << std::endl;<br>
> return 1;<br>
> }<br>
><br>
> const unsigned int Dimension = 3;<br>
> typedef float PixelType;<br>
><br>
> typedef itk::Image< PixelType, Dimension > FixedImageType;<br>
> typedef itk::ImageFileReader< FixedImageType > FixedImageReaderType;<br>
> FixedImageReaderType::Pointer fixedImageReader =<br>
> FixedImageReaderType::New();<br>
> fixedImageReader->SetFileName( argv[1] );<br>
><br>
> typedef itk::Image< PixelType, Dimension > MovingImageType;<br>
> typedef itk::ImageFileReader< MovingImageType > MovingImageReaderType;<br>
> MovingImageReaderType::Pointer movingImageReader =<br>
> MovingImageReaderType::New();<br>
> movingImageReader->SetFileName( argv[2] );<br>
> typedef itk::Vector< float, Dimension > VectorPixelType;<br>
> typedef itk::Image< VectorPixelType, Dimension > DeformationFieldType;<br>
> typedef itk::ImageFileReader< DeformationFieldType > FieldReaderType;<br>
> FieldReaderType::Pointer fieldReader = FieldReaderType::New();<br>
> fieldReader->SetFileName( argv[3] );<br>
><br>
> try {<br>
> fieldReader->Update();<br>
> movingImageReader->Update();<br>
> fixedImageReader->Update();<br>
> }<br>
> catch ( itk::ExceptionObject & err )<br>
> {<br>
> std::cerr << "ExceptionObject caught !" << std::endl;<br>
> std::cerr << err << std::endl;<br>
> return -1;<br>
> }<br>
><br>
><br>
> typedef itk::WarpImageFilter<FixedImageType,<br>
> MovingImageType,DeformationFieldType > WarperType;<br>
> WarperType::Pointer warper = WarperType::New();<br>
> typedef WarperType::CoordRepType WarperCoordRepType;<br>
> typedef itk::LinearInterpolateImageFunction<MovingImageType,<br>
> WarperCoordRepType > InterpolatorType;<br>
> InterpolatorType::Pointer interpolator = InterpolatorType::New();<br>
> FixedImageType::Pointer fixedImage = fixedImageReader->GetOutput();<br>
> FixedImageType::Pointer warpedImage;<br>
> warper->SetInput( movingImageReader->GetOutput() );<br>
> warper->SetDeformationField( fieldReader->GetOutput() );<br>
> warper->SetInterpolator( interpolator );<br>
> warper->SetOutputOrigin( fixedImage->GetOrigin() );<br>
> warper->SetOutputSpacing( fixedImage->GetSpacing() );<br>
> warper->SetOutputDirection( fixedImage->GetDirection() );<br>
> MovingImageType::PixelType padValue = 0;<br>
> warper->SetEdgePaddingValue( padValue );<br>
> warper->Update();<br>
> warpedImage = warper->GetOutput();<br>
> /*typedef itk::ImageFileWriter< MovingImageType > WriterType;<br>
><br>
> WriterType::Pointer writer = WriterType::New();<br>
> writer->SetFileName( argv[4] );<br>
> writer->SetInput( warper->GetOutput() );<br>
> writer->Update();*/<br>
><br>
><br>
> typedef itk::MinimumMaximumImageFilter<MovingImageType> MinMaxFilterType;<br>
> MinMaxFilterType::Pointer minMaxFilter = MinMaxFilterType::New();<br>
> minMaxFilter->SetInput( warpedImage );<br>
> minMaxFilter->Update();<br>
> float min = minMaxFilter->GetMinimum();<br>
> double shift = -1.0 * static_cast<double>( min );<br>
> double scale = static_cast<double>( minMaxFilter->GetMaximum() );<br>
> scale += shift;<br>
> scale = 255.0 / scale;<br>
> typedef itk::ShiftScaleImageFilter<MovingImageType, MovingImageType><br>
> FilterType;<br>
> FilterType::Pointer filter = FilterType::New();<br>
> filter->SetInput( warpedImage );<br>
> filter->SetShift( shift );<br>
> filter->SetScale( scale );<br>
> filter->Update();<br>
><br>
> typedef itk::Image<PixelType,Dimension> WriteImageType;<br>
> typedef itk::CastImageFilter<MovingImageType,WriteImageType> CasterType1;<br>
> CasterType1::Pointer caster1 = CasterType1::New();<br>
> caster1->SetInput(filter->GetOutput());<br>
> caster1->Update();<br>
> typedef itk::ImageFileWriter<WriteImageType> WriterType;<br>
> WriterType::Pointer writer = WriterType::New();<br>
> writer->SetFileName(argv[4]);<br>
> writer->SetInput(caster1->GetOutput() );<br>
> writer->Write();<br>
> return 0;<br>
> }<br>
><br>
><br>
><br>
</div></div>> _____________________________________<br>
> Powered by <a href="http://www.kitware.com" target="_blank">www.kitware.com</a><br>
><br>
> Visit other Kitware open-source projects at<br>
> <a href="http://www.kitware.com/opensource/opensource.html" target="_blank">http://www.kitware.com/opensource/opensource.html</a><br>
><br>
> Kitware offers ITK Training Courses, for more information visit:<br>
> <a href="http://www.kitware.com/products/protraining.html" target="_blank">http://www.kitware.com/products/protraining.html</a><br>
><br>
> Please keep messages on-topic and check the ITK FAQ at:<br>
> <a href="http://www.itk.org/Wiki/ITK_FAQ" target="_blank">http://www.itk.org/Wiki/ITK_FAQ</a><br>
><br>
> Follow this link to subscribe/unsubscribe:<br>
> <a href="http://www.itk.org/mailman/listinfo/insight-users" target="_blank">http://www.itk.org/mailman/listinfo/insight-users</a><br>
><br>
><br>
</blockquote></div><br></div>