Hi Luis,<div><br></div><div>I am still struggling with this, do you have any ideas?</div><div><br></div><div>Thanks,</div><div><br></div><div>Steve<br><br><div class="gmail_quote">On Thu, Mar 18, 2010 at 10:37 PM, Steve Lancey <span dir="ltr">&lt;<a href="mailto:steve.lancey@gmail.com">steve.lancey@gmail.com</a>&gt;</span> wrote:<br>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex;">Hi Luis,<div><br></div><div>thanks for your answer:</div><div><br></div><div>1)</div><div class="im"><div>&gt;&gt;In the code that you posted, we see that you are<br>
&gt;&gt;reading the deformation field from a file and passing<br>&gt;&gt;it to the WrapImageFilter:<br>
&gt;&gt;    warper-&gt;SetDeformationField( fieldReader-&gt;GetOutput() );<br>&gt;&gt;Could you please describe how you generated that file ?</div></div><div>I am using the following code:</div><div><br></div><div>typedef itk::fem::FEMRegistrationFilter&lt;ImageType,ImageType&gt; RegistrationType;</div>

<div>...</div><div>...</div><div>...</div><div>  RegistrationType::Pointer registrationFilter = RegistrationType::New(); </div><div>...</div><div>...</div><div>...</div><div><div>  registrationFilter-&gt;WriteWarpedImage(</div>

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