Hi all,<br><br>By referring to &quot;itkInverseDisplacementFieldImageFilterTest.cxx&quot;, I tried the following code piece for estimating the inverse of a BSpline transform through using itk::InverseDisplacementFieldImageFilter.  <br>
<br>However, the code crashes on the line <br>&quot;vFilter-&gt;UpdateLargestPossibleRegion();&quot;<br><br>Please point me the reason, thanks!<br><br><br>//====== Start of the code piece ===========================================<br>
void InvertBSpline(const ImageType::Pointer vMovImg, const ImageType::Pointer vFixImg, const BSplineDeformableTransform::Pointer vTransform)<br>{<br>    typedef itk::Vector&lt;float,3&gt;        VectorType;<br>    typedef itk::Image&lt;VectorType, 3&gt;   DisplacementFieldType;<br>
    DisplacementFieldType::Pointer vDisplacementField = DisplacementFieldType::New();<br>    vDisplacementField-&gt;CopyInformation(vFixImg);<br>    vDisplacementField-&gt;Allocate();<br>    typedef itk::ImageRegionIteratorWithIndex&lt;ItkRawType&gt; FixIteratorType;<br>
    typedef itk::ImageRegionIteratorWithIndex&lt;DisplacementFieldType&gt; DisIteratorType;<br>    FixIteratorType itFix(vFixImg, vFixImg-&gt;GetRequestedRegion());<br>    DisIteratorType itDis(vDisplacementField, vDisplacementField-&gt;GetRequestedRegion());<br>
    for (itFix.GoToBegin(), itDis.GoToBegin(); !itFix.IsAtEnd(); ++itFix, ++itDis)<br>    {<br>        ImageType::IndexType idx = itFix.GetIndex();<br>        itk::Point&lt;double,3&gt; vFixPoint, vMovPoint;<br>        vFixImg-&gt;TransformIndexToPhysicalPoint(idx, vFixPoint);<br>
        vMovPoint = vTransform-&gt;TransformPoint(vFixPoint);<br>        VectorType vTmp;<br>        vTmp[0] = vMovPoint[0] - vFixPoint[0];<br>        vTmp[1] = vMovPoint[1] - vFixPoint[1];<br>        vTmp[2] = vMovPoint[2] - vFixPoint[2];<br>
        itDis.Set(vTmp);<br>    }<br><br>    //- Create inverse displacement filter<br>    //<br>    typedef itk::InverseDisplacementFieldImageFilter&lt;<br>                                    DisplacementFieldType,<br>                                    DisplacementFieldType&gt;  FilterType;<br>
    FilterType::Pointer vFilter = FilterType::New();<br>    vFilter-&gt;SetOutputSpacing(vMovImg-&gt;GetSpacing());<br>    vFilter-&gt;SetOutputOrigin(vMovImg-&gt;GetOrigin());<br>    vFilter-&gt;SetSize(vMovImg-&gt;GetLargestPossibleRegion().GetSize());<br>
    vFilter-&gt;SetInput(vDisplacementField);<br>    vFilter-&gt;SetSubsamplingFactor(16);<br><br>    try<br>    {<br>        vFilter-&gt;UpdateLargestPossibleRegion();<br>    }<br>    catch( itk::ExceptionObject &amp; excp )<br>
    {<br>        std::cerr &lt;&lt; &quot;Exception thrown &quot; &lt;&lt; std::endl;<br>        std::cerr &lt;&lt; excp &lt;&lt; std::endl;<br>    }<br>}<br><br>//====== End of the code piece =============================================<br>
<br>