<div dir="ltr">Emma,<div><br></div><div>It is likely that your settings for the Grid of the BSplineDeformableTransform are not fitting well to the FixedImage.</div><div><br></div><div>You may want to use the helper class:</div>
<div><a href="http://www.itk.org/Doxygen/html/classitk_1_1BSplineTransformInitializer.html">http://www.itk.org/Doxygen/html/classitk_1_1BSplineTransformInitializer.html</a><br></div><div><br></div><div>that will take care of doing this properly, so you don&#39;t have to adjust those parameters by hand.</div>
<div><br></div><div>More details are available in the article in the Insight Journal</div><div><a href="http://www.insight-journal.org/browse/publication/216">http://www.insight-journal.org/browse/publication/216</a><br></div>
<div><br></div><div><br></div><div>   Thanks</div><div><br></div><div>        Luis</div></div><div class="gmail_extra"><br><br><div class="gmail_quote">On Wed, Dec 11, 2013 at 5:26 PM, Emma Saunders <span dir="ltr">&lt;<a href="mailto:emmasaunders123@gmail.com" target="_blank">emmasaunders123@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"><div dir="ltr">Hi everyone<div><br></div><div>I appear to be getting a step artefact at the boundaries of my registered image, using a multi-resolution Bspline approach.  The distance of the artefact co-indices with the origin of the original image.  I can&#39;t seem to find why this occurs.  I set the Bspline grid wrt the original fixed image so unsure why this is happening.</div>

<div><br></div><div>I would really appreciate any advice:</div><div><br></div><div>Below is the code if anyone is interested:</div><div><br></div><div> Many Thanks</div><div><br></div>
<div>Emma</div><div><br></div><div><br></div><div><br></div><div><br></div><div><br></div><div><br></div><div><div>#include &quot;itkMedianImageFilter.h&quot;</div><div>#include &quot;itkDiscreteGaussianImageFilter.h&quot;</div>

<div>#include &quot;itkImageRegistrationMethod.h&quot;</div><div>#include &quot;itkTimeProbesCollectorBase.h&quot;</div><div>#include &quot;itkSpatialObjectToImageFilter.h&quot;</div><div>#include &quot;itkEllipseSpatialObject.h&quot;</div>

<div>#include &quot;itkRegularStepGradientDescentOptimizer.h&quot; </div><div>#include &quot;itkImageFileReader.h&quot;</div><div>#include &quot;itkImageFileWriter.h&quot;</div><div>#include &quot;itkRegularStepGradientDescentOptimizer.h&quot; </div>

<div>#include &quot;itkAffineTransform.h&quot;</div><div>#include &quot;itkBSplineDeformableTransform.h&quot; //version 3</div><div>#include &quot;itkMattesMutualInformationImageToImageMetric.h&quot;</div><div>#include &quot;itkTimeProbesCollectorBase.h&quot;</div>

<div>#include &quot;itkMemoryProbesCollectorBase.h&quot;</div><div>#include &lt;itkHistogramMatchingImageFilter.h&gt;</div><div>#include &quot;itkBSplineResampleImageFunction.h&quot;</div><div>#include &quot;itkLBFGSBOptimizer.h&quot;</div>

<div>#include &quot;itkMattesMutualInformationImageToImageMetric.h&quot;</div><div>#include &quot;itkLBFGSOptimizer.h&quot;</div><div>#include &quot;itkImageFileWriter.h&quot;</div><div>#include &quot;itkResampleImageFilter.h&quot;</div>

<div>#include &quot;itkCastImageFilter.h&quot;</div><div>#include &quot;itkSquaredDifferenceImageFilter.h&quot;</div><div>#include &quot;itkBSplineResampleImageFunction.h&quot;</div><div>#include &quot;itkIdentityTransform.h&quot;</div>

<div>#include &quot;itkBSplineDecompositionImageFilter.h&quot;</div><div>#include &quot;itkCommand.h&quot;</div><div>#include &lt;fstream&gt;</div><div>#include &quot;itkCommand.h&quot;</div><div><br></div><div>unsigned int Counter = 0;</div>

<div>unsigned int Countold = 0;</div><div><br></div><div><br></div><div>//This is for the affine registration</div><div>class CommandIterationUpdate : public itk::Command</div><div>{</div><div>public:</div><div>  typedef  CommandIterationUpdate   Self;</div>

<div>  typedef  itk::Command             Superclass;</div><div>  typedef itk::SmartPointer&lt;Self&gt;   Pointer;</div><div>  itkNewMacro( Self );</div><div>protected:</div><div>  CommandIterationUpdate() {};</div><div>public:</div>

<div>  typedef itk::RegularStepGradientDescentOptimizer OptimizerType;</div><div>  typedef   const OptimizerType *                  OptimizerPointer;</div><div><br></div><div>  void Execute(itk::Object *caller, const itk::EventObject &amp; event)</div>

<div>    {</div><div>    Execute( (const itk::Object *)caller, event);</div><div>    }</div><div><br></div><div>  void Execute(const itk::Object * object, const itk::EventObject &amp; event)</div><div>    {</div><div>    OptimizerPointer optimizer =</div>

<div>                      dynamic_cast&lt; OptimizerPointer &gt;( object );</div><div>    if( ! itk::IterationEvent().CheckEvent( &amp;event ) )</div><div>      {</div><div>      return;</div><div>      }</div><div>     std::cout &lt;&lt; optimizer-&gt;GetValue() &lt;&lt; &quot;   &quot; &lt;&lt;std::endl;</div>

<div><br></div><div>    }</div><div>};</div><div><br></div><div> </div><div>//This is for B spline</div><div>   class BCommandIterationUpdate : public itk::Command</div><div>  {</div><div>  public:</div><div>    typedef  BCommandIterationUpdate                     Self;</div>

<div>    typedef  itk::Command                               Superclass;</div><div>    typedef  itk::SmartPointer&lt;BCommandIterationUpdate&gt;  Pointer;</div><div>    itkNewMacro( BCommandIterationUpdate );</div><div>  protected:</div>

<div>    BCommandIterationUpdate() {};</div><div>    </div><div>  public:</div><div>  </div><div>  typedef itk::LBFGSBOptimizer <span style="white-space:pre-wrap">                                        </span>   OptimizerType;</div><div>  typedef   const OptimizerType *                  OptimizerPointer;</div>

<div><br></div><div><br></div><div>  public:</div><div><br></div><div>    void Execute(itk::Object *caller, const itk::EventObject &amp; event)</div><div>      {</div><div>        Execute( (const itk::Object *)caller, event);</div>

<div>      }</div><div><br></div><div>    void Execute(const itk::Object * object, const itk::EventObject &amp; event)</div><div>      {</div><div><span style="white-space:pre-wrap">                </span>  </div><div><span style="white-space:pre-wrap">        </span>const OptimizerPointer <span style="white-space:pre-wrap">        </span> optimizer =</div>

<div>                      dynamic_cast&lt; OptimizerPointer &gt;( object );</div><div><span style="white-space:pre-wrap">                </span>  </div><div><span style="white-space:pre-wrap">                </span> </div><div>        if( !(itk::IterationEvent().CheckEvent( &amp;event )) )</div>

<div>          {</div><div>          return;</div><div>          }</div><div>          </div><div>        Countold=Counter;</div><div>        Counter = optimizer-&gt;GetCurrentIteration();</div><div>        </div><div>        if (Counter!=Countold)</div>

<div>        {</div><div> </div><div>        std::cout  &lt;&lt;  optimizer-&gt;GetValue() &lt;&lt; &quot; &quot; &lt;&lt;std::endl;</div><div><span style="white-space:pre-wrap">        </span>}</div><div>      }</div><div>  };</div>

<div>  </div><div>  </div><div><br></div><div>int main( int argc, char *argv[] )</div><div>{</div><div>  if( argc &lt; 4 )</div><div>    {</div><div><br></div><div>    std::cerr &lt;&lt; &quot;Missing Parameters &quot; &lt;&lt; std::endl;</div>

<div>    std::cerr &lt;&lt; &quot;Usage: &quot; &lt;&lt; argv[0];</div><div>    std::cerr &lt;&lt; &quot;   fixedImageFile  movingImageFile &quot; &lt;&lt; std::endl;</div><div>    std::cerr &lt;&lt; &quot;   outputImagefile&quot; &lt;&lt; std::endl;</div>

<div>     std::cerr &lt;&lt; &quot;   outputDeformationField &quot; &lt;&lt; std::endl;</div><div>     </div><div>    return EXIT_FAILURE;</div><div>   </div><div>    }</div><div><br></div><div> // READ AND CREATE THE IMAGES AND DISPLACEMENT FIELD</div>

<div> </div><div>  const    unsigned int    Dimension = 3;</div><div>  typedef float        PixelType;</div><div>  typedef itk::Vector&lt; double, 3 &gt;           VectorPixelType;</div><div>   </div><div>  typedef itk::Image&lt; PixelType, Dimension &gt;  FixedImageType;</div>

<div>  typedef itk::Image&lt; PixelType, Dimension &gt;  MovingImageType;</div><div>  typedef itk::Image&lt; PixelType, Dimension &gt;  OutputImageType;</div><div>  typedef itk::Image&lt;  VectorPixelType, 3 &gt; DisplacementFieldType;</div>

<div><br></div><div>  typedef itk::ImageFileReader&lt; FixedImageType  &gt; FixedImageReaderType;</div><div>  typedef itk::ImageFileReader&lt; MovingImageType &gt; MovingImageReaderType;</div><div>  typedef itk::ImageFileWriter&lt; OutputImageType &gt;  ImageWriterType;</div>

<div>  typedef itk::ImageFileWriter&lt; DisplacementFieldType &gt;  FieldWriterType;</div><div><br></div><div>  typedef itk::ImageFileWriter&lt; DisplacementFieldType &gt;  FieldWriterType;</div><div>  FieldWriterType::Pointer FieldWriter = FieldWriterType::New();</div>

<div><br></div><div>  FixedImageReaderType::Pointer  FixedImageReader  = FixedImageReaderType::New();</div><div>  MovingImageReaderType::Pointer MovingImageReader = MovingImageReaderType::New();</div><div>  ImageWriterType::Pointer ImageWriter = ImageWriterType::New();</div>

<div>  FixedImageReader-&gt;SetFileName(  argv[1] );</div><div>  MovingImageReader-&gt;SetFileName( argv[2] );</div><div>  ImageWriter-&gt;SetFileName( argv[3] );</div><div>  FieldWriter-&gt;SetFileName( argv[4] );</div>
<div>
  </div><div>  </div><div>  </div><div>  FixedImageType::Pointer FixedImage = FixedImageReader-&gt;GetOutput();</div><div>  MovingImageType::Pointer MovedImage = MovingImageReader-&gt;GetOutput();</div><div>  OutputImageType::Pointer OutImage = OutputImageType::New();</div>

<div>  DisplacementFieldType::Pointer DeformField = DisplacementFieldType::New();</div><div>  std::cout &lt;&lt; &quot;First Fixed Image origin is &quot; &lt;&lt; FixedImage-&gt;GetOrigin() &lt;&lt; std::endl ;</div><div>

  </div><div>  FixedImage-&gt;Update();</div><div>  MovedImage-&gt;Update();</div><div>  </div><div>  </div><div><br></div><div>  </div><div>  //Median and Gaussian Filter</div><div>  </div><div> typedef itk::MedianImageFilter&lt;FixedImageType, FixedImageType &gt; FilterType;</div>

<div><br></div><div> </div><div>  // Create and setup a median filter</div><div>  FilterType::Pointer medianFilterFixed = FilterType::New();</div><div>  FilterType::InputSizeType radius;</div><div>  radius.Fill(1);</div>
<div>
  medianFilterFixed-&gt;SetRadius(radius);</div><div>  medianFilterFixed-&gt;SetInput( FixedImage );</div><div>  medianFilterFixed-&gt;Update();</div><div> </div><div>   FilterType::Pointer medianFilterMoved = FilterType::New();</div>

<div>  medianFilterMoved-&gt;SetRadius(radius);</div><div>  medianFilterMoved-&gt;SetInput( MovedImage );</div><div>  medianFilterMoved-&gt;Update();</div><div> </div><div>  typedef itk::DiscreteGaussianImageFilter&lt;</div>

<div>                                      FixedImageType,</div><div>                                      FixedImageType</div><div>                                                    &gt; GaussianFilterType;</div><div>                                                    </div>

<div>              </div><div>  GaussianFilterType::Pointer fixedSmoother  = GaussianFilterType::New();</div><div>  GaussianFilterType::Pointer movingSmoother = GaussianFilterType::New();</div><div><br></div><div>  fixedSmoother-&gt;SetVariance(1.0 );</div>

<div>  movingSmoother-&gt;SetVariance(1.0 );</div><div><br></div><div>  fixedSmoother-&gt;SetInput(  medianFilterFixed-&gt;GetOutput() );</div><div>  movingSmoother-&gt;SetInput ( medianFilterMoved-&gt;GetOutput() );</div>

<div><br></div><div>  fixedSmoother-&gt;Update();</div><div>  movingSmoother-&gt;Update();</div><div><br></div><div>   std::cout&lt;&lt; &quot;The fixed origin is &quot; &lt;&lt; fixedSmoother-&gt;GetOutput()-&gt;GetOrigin() &lt;&lt;std::endl;</div>

<div>   std::cout &lt;&lt; &quot;Fixed Image origin is &quot; &lt;&lt; FixedImage-&gt;GetOrigin() &lt;&lt; std::endl ;</div><div><br></div><div><br></div><div>//Define the transformations to be used</div><div><br></div><div>

//PERFROM AN INITIAL AFFINE</div><div><br></div><div>  typedef itk::AffineTransform&lt;</div><div>                                  double,</div><div>                                  Dimension  &gt;     AffineTransformType;</div>

<div>                                  </div><div>                                  </div><div>//This will be followed by Bspline              </div><div>                                 </div><div>                                 </div>

<div>  const unsigned int SpaceDimension = Dimension;</div><div>  const unsigned int SplineOrder = 3;</div><div>  typedef double CoordinateRepType;</div><div><br></div><div>  typedef itk::BSplineDeformableTransform&lt;</div>

<div>                            CoordinateRepType,</div><div>                            SpaceDimension,</div><div>                            SplineOrder &gt;     DeformableTransformType; //version 3</div><div><br></div>

<div>             </div><div> //Now for the Optimizer, Metric and Interpolatror  and Registration                               </div><div> //which are connected to the registration</div><div>   </div><div> //typedef itk::LBFGSOptimizer       bOptimizerType;</div>

<div>  typedef itk::LBFGSBOptimizer       bOptimizerType;</div><div>   bOptimizerType::Pointer      boptimizer     = bOptimizerType::New();</div><div><br></div><div> typedef itk::RegularStepGradientDescentOptimizer       OptimizerType;</div>

<div>  </div><div>  typedef itk::MattesMutualInformationImageToImageMetric&lt;</div><div>                                    FixedImageType,</div><div>                                    MovingImageType &gt;    MetricType;</div>

<div>                                   </div><div>  </div><div>  typedef itk:: LinearInterpolateImageFunction&lt;</div><div>                                   MovingImageType,</div><div>                                   double          &gt;    InterpolatorType;</div>

<div>                                   </div><div>  typedef itk::ImageRegistrationMethod&lt;</div><div>                                    FixedImageType,</div><div>                                    MovingImageType &gt;    RegistrationType;</div>

<div><br></div><div>  MetricType::Pointer         metric        = MetricType::New();</div><div>  OptimizerType::Pointer      optimizer     = OptimizerType::New();</div><div>  InterpolatorType::Pointer   interpolator  = InterpolatorType::New();</div>

<div>  RegistrationType::Pointer   registration  = RegistrationType::New();</div><div><br></div><div>  // Setup the metric parameters</div><div>  </div><div>  metric-&gt;ReinitializeSeed( 76926294 );</div><div>  metric-&gt;SetUseCachingOfBSplineWeights(1);</div>

<div><br></div><div>    const unsigned long numberOfImagePixels = </div><div>    FixedImage-&gt;GetLargestPossibleRegion().GetNumberOfPixels();</div><div><br></div><div>      const unsigned long numberOfSpatialSamplesAffine =</div>

<div>    static_cast&lt;unsigned long&gt;( numberOfImagePixels * 1 );</div><div><br></div><div>  metric-&gt;SetNumberOfSpatialSamples( numberOfSpatialSamplesAffine );</div><div>  metric-&gt;SetNumberOfHistogramBins( 40 ); // Reduce the number of bins if a deformable registration fails. If the number of bins is too large, the estimated PDFs will be a field of impulses and will inhibit reliable registration estimation.</div>

<div><br></div><div>  registration-&gt;SetMetric(        metric        );</div><div>  registration-&gt;SetOptimizer(     optimizer     );</div><div>  registration-&gt;SetInterpolator(  interpolator  );</div><div><br></div>

<div>//create the Affine transform</div><div><br></div><div>  AffineTransformType::Pointer  Affinetransform = AffineTransformType::New();</div><div>  registration-&gt;SetTransform( Affinetransform );  </div><div> </div><div>

//input the images</div><div><br></div><div>  registration-&gt;SetFixedImage(fixedSmoother-&gt;GetOutput());</div><div>  registration-&gt;SetMovingImage(movingSmoother-&gt;GetOutput());</div><div>  </div><div>//input the region</div>

<div><br></div><div>  registration-&gt;SetFixedImageRegion(FixedImage-&gt;GetBufferedRegion() );</div><div><br></div><div>//set the initial parameters</div><div>  typedef RegistrationType::ParametersType AffineParametersType;</div>

<div>  AffineParametersType initialParameters( Affinetransform-&gt;GetNumberOfParameters() );</div><div>  </div><div>  Affinetransform-&gt;SetIdentity();</div><div>  registration-&gt;SetInitialTransformParameters(</div><div>

                             Affinetransform-&gt;GetParameters() );</div><div>  </div><div><br></div><div>//set scale values for optimizer to search sensibly and set up optimizer</div><div> </div><div>  typedef OptimizerType::ScalesType       OptimizerScalesType;</div>

<div>  OptimizerScalesType optimizerScales( Affinetransform-&gt;GetNumberOfParameters() );</div><div><br></div><div>  FixedImageType:: SpacingType fixedspacing = FixedImage-&gt;GetSpacing();</div><div>  FixedImageType::RegionType region=FixedImage-&gt;GetLargestPossibleRegion(); </div>

<div>  const unsigned int numberOfPixels = region.GetNumberOfPixels();</div><div>  FixedImageType:: SizeType fixedsize =region.GetSize();</div><div><br></div><div><br></div><div>  double translationScalea = 1.0 /( 10.0 * fixedspacing[0]  *  fixedsize[0]  );</div>

<div>  double translationScaleb = 1.0 /( 10.0 * fixedspacing[1]  *  fixedsize[1]  );</div><div>  double translationScalec = 1.0 /( 10.0 * fixedspacing[2]  *  fixedsize[2]  );</div><div><br></div><div><br></div><div>//for 3D use below</div>

<div>  optimizerScales[0] =  1.0;</div><div>  optimizerScales[1] =  1.0;</div><div>  optimizerScales[2] =  1.0;</div><div>  optimizerScales[3] =  1.0;</div><div>  optimizerScales[4] =  1.0;</div><div>  optimizerScales[5] =  1.0;</div>

<div>  optimizerScales[6] =  1.0;</div><div>  optimizerScales[7] =  1.0;</div><div>  optimizerScales[8] =  1.0;</div><div>  optimizerScales[9]  =  translationScalea;</div><div>  optimizerScales[10] =  translationScaleb;</div>

<div>  optimizerScales[11] =  translationScalec;</div><div><br></div><div>  //step length was 0.1 larger step length quicker converge but may get erreaic values in the metric as it decreases</div><div>  unsigned int maxNumberOfIterations = 60; </div>

<div>  optimizer-&gt;SetMaximumStepLength( 0.05 );</div><div>  optimizer-&gt;SetMinimumStepLength( 0.001 ); //was 0.00001</div><div>  optimizer-&gt;SetNumberOfIterations( maxNumberOfIterations );</div><div>  optimizer-&gt;MaximizeOff();</div>

<div><br></div><div><br></div><div>//Perform Affine Registration put observer on the Affine optimizer</div><div><br></div><div>  CommandIterationUpdate::Pointer observera = CommandIterationUpdate::New();</div><div>  optimizer-&gt;AddObserver( itk::IterationEvent(), observera );</div>

<div><br></div><div>  try</div><div>    {</div><div>    registration-&gt;Update();</div><div>    std::cout &lt;&lt; &quot; &quot; &lt;&lt;std::endl;</div><div>    std::cout &lt;&lt; &quot;Optimizer stop condition: &quot;</div>

<div>           &lt;&lt; registration-&gt;GetOptimizer()-&gt;GetStopConditionDescription()</div><div>           &lt;&lt; std::endl;</div><div>    }</div><div>  catch( itk::ExceptionObject &amp; err )</div><div>    {</div>

<div>    std::cerr &lt;&lt; &quot;ExceptionObject caught !&quot; &lt;&lt; std::endl;</div><div>    std::cerr &lt;&lt; err &lt;&lt; std::endl;</div><div>    return EXIT_FAILURE;</div><div>    }</div><div><br></div><div>  std::cout &lt;&lt; &quot;Affine Registration completed&quot; &lt;&lt; std::endl;</div>

<div>  std::cout &lt;&lt; std::endl;</div><div>  optimizer-&gt;RemoveAllObservers();</div><div>  optimizer-&gt;RemoveAllObservers();</div><div><br></div><div>//Obtain the last Affine transform which will be used as a bulk transform in the Bspline</div>

<div><br></div><div>  Affinetransform-&gt;SetParameters(registration-&gt;GetLastTransformParameters() );</div><div><br></div><div>   </div><div>//////////////////////////////////////////////////////////////////////////////////////////////</div>

<div>///NOW FOR THE Bspline Registration</div><div>//  We construct two transform objects, each one will be configured for a resolution level.</div><div>//  Notice than in this multi-resolution scheme we are not modifying the</div>

<div>//  resolution of the image, but rather the flexibility of the deformable</div><div>//  transform itself.</div><div><br></div><div>//SO must define a transform and a grid for each level</div><div><br></div><div>  DeformableTransformType::Pointer  BtransformLow = DeformableTransformType::New();</div>

<div>  typedef DeformableTransformType::RegionType RegionType;</div><div>  RegionType bsplineRegionlow;</div><div>  RegionType::SizeType   gridBorderSize;</div><div>  RegionType::SizeType   gridSizeOnImagelow;</div><div>
  RegionType::SizeType   totalGridSizelow;</div>
<div>  </div><div>  </div><div>//.////////////////////////////////////////////////////////////////////// </div><div>///GRID SIZE FOR LOW RESOLUTION  </div><div>/////////////////////////////////////////////////////////////////////////  </div>

<div> </div><div>  gridSizeOnImagelow[0]=12;  //lat coronal, sup-inf //want 5 2 7 10 4 15</div><div>  gridSizeOnImagelow[1]=4;    // want</div><div>  gridSizeOnImagelow[2]=12;</div><div>  gridBorderSize.Fill( 3 );    // Border for spline order = 3 ( 1 lower, 2 upper )</div>

<div>  totalGridSizelow = gridSizeOnImagelow + gridBorderSize;</div><div>  bsplineRegionlow.SetSize( totalGridSizelow );</div><div>  typedef DeformableTransformType::SpacingType SpacingTypelow;</div><div>  SpacingTypelow spacinglow = FixedImage-&gt;GetSpacing();</div>

<div>  typedef DeformableTransformType::OriginType OriginTypelow;</div><div>  OriginTypelow originlow = FixedImage-&gt;GetOrigin();</div><div>  FixedImageType::RegionType FixedRegionlow = FixedImage-&gt;GetBufferedRegion();</div>

<div>  FixedImageType::SpacingType   fixedSpacinglow    = FixedImage-&gt;GetSpacing();</div><div>  FixedImageType::SizeType FixedImageSizelow = FixedRegionlow.GetSize();</div><div><br></div><div>  for(unsigned int r=0; r&lt;Dimension; r++)</div>

<div>    {</div><div>    </div><div>    spacinglow[r] = floor( fixedSpacinglow[r] * (FixedImageSizelow[r] - 1) / (gridSizeOnImagelow[r] - 1) ); //new way</div><div>    </div><div>    }</div><div>    </div><div>  FixedImageType::DirectionType gridDirection = FixedImage-&gt;GetDirection();</div>

<div>  SpacingTypelow gridOriginOffsetlow = gridDirection * spacinglow;</div><div>  OriginTypelow gridOriginlow = originlow - gridOriginOffsetlow;</div><div>  BtransformLow-&gt;SetGridSpacing( spacinglow );</div><div>  BtransformLow-&gt;SetGridOrigin( gridOriginlow );</div>

<div>  BtransformLow-&gt;SetGridRegion( bsplineRegionlow ); //i.e region of gridsize</div><div>  BtransformLow-&gt;SetGridDirection( gridDirection );</div><div>  std::cout &lt;&lt; &quot;grid region for low Bspline is &quot; &lt;&lt; BtransformLow-&gt;GetGridRegion() &lt;&lt;std::endl;</div>

<div><br></div><div><br></div><div>//////////////////////////////////////////////////////////////////////////////</div><div>///GRID FOR MEDIUM</div><div>////////////////////////////////////////////////////////////////////////////////</div>

<div><br></div><div>  DeformableTransformType::Pointer  BtransformMed= DeformableTransformType::New();</div><div>  RegionType bsplineRegionMed;</div><div>  RegionType::SizeType   gridSizeOnImageMed;</div><div>  RegionType::SizeType   totalGridSizeMed;</div>

<div> </div><div>//.////////////////////////////////////////////////////////////////////// </div><div>///GRID SIZE FOR MEDIUM RESOLUTION  </div><div>/////////////////////////////////////////////////////////////////////////  </div>

<div>  gridSizeOnImageMed[0]=24;  //lat coronal, sup-inf</div><div>  gridSizeOnImageMed[1]=8;    // want</div><div>  gridSizeOnImageMed[2]=24;</div><div>  totalGridSizeMed = gridSizeOnImageMed + gridBorderSize;</div><div>

  bsplineRegionMed.SetSize( totalGridSizeMed );</div><div> typedef DeformableTransformType::SpacingType SpacingTypeMed;</div><div>  SpacingTypeMed spacingMed = FixedImage-&gt;GetSpacing();</div><div>  typedef DeformableTransformType::OriginType OriginTypeMed;</div>

<div>  OriginTypeMed originMed = FixedImage-&gt;GetOrigin();</div><div>  FixedImageType::RegionType FixedRegionMed = FixedImage-&gt;GetBufferedRegion();</div><div>  FixedImageType::SpacingType   fixedSpacingMed    = FixedImage-&gt;GetSpacing();</div>

<div>  FixedImageType::SizeType FixedImageSizeMed = FixedRegionMed.GetSize();</div><div> </div><div>      for(unsigned int r=0; r&lt;Dimension; r++)</div><div>    {</div><div>    </div><div>    spacingMed[r] = floor( fixedSpacingMed[r] * (FixedImageSizeMed[r] - 1) / (gridSizeOnImageMed[r] - 1) ); //new way</div>

<div>    </div><div>    }</div><div>   </div><div>   SpacingTypeMed gridOriginOffsetMed = gridDirection * spacingMed; </div><div>   OriginTypeMed gridOriginMed = originMed - gridOriginOffsetMed;</div><div> </div><div>  BtransformMed-&gt;SetGridSpacing( spacingMed );</div>

<div>  BtransformMed-&gt;SetGridOrigin( gridOriginMed );</div><div>  BtransformMed-&gt;SetGridRegion( bsplineRegionMed ); //i.e region of gridsize</div><div>  BtransformMed-&gt;SetGridDirection( gridDirection );</div><div>

  std::cout &lt;&lt; &quot;grid region for Med Bspline is &quot; &lt;&lt; BtransformMed-&gt;GetGridRegion() &lt;&lt;std::endl;</div><div><br></div><div>//////////////////////////////////////////////////////////////////////////////</div>

<div>///GRID FOR HIGH</div><div>////////////////////////////////////////////////////////////////////////////////</div><div><br></div><div>//Now the Higher Transform</div><div><br></div><div>  DeformableTransformType::Pointer  BtransformHigh = DeformableTransformType::New();</div>

<div>  RegionType bsplineRegionhigh;</div><div>  RegionType::SizeType   gridSizeOnImagehigh;</div><div>  RegionType::SizeType   totalGridSizehigh;</div><div> </div><div>//.////////////////////////////////////////////////////////////////////// </div>

<div>///GRID SIZE FOR HIGH RESOLUTION  </div><div>/////////////////////////////////////////////////////////////////////////  </div><div>  gridSizeOnImagehigh[0]=50;  //lat coronal, sup-inf</div><div>  gridSizeOnImagehigh[1]=16;    // want</div>

<div>  gridSizeOnImagehigh[2]=50;</div><div>  totalGridSizehigh = gridSizeOnImagehigh + gridBorderSize;</div><div>  bsplineRegionhigh.SetSize( totalGridSizehigh );</div><div> typedef DeformableTransformType::SpacingType SpacingTypehigh;</div>

<div>  SpacingTypehigh spacinghigh = FixedImage-&gt;GetSpacing();</div><div>  typedef DeformableTransformType::OriginType OriginTypehigh;</div><div>  OriginTypehigh originhigh = FixedImage-&gt;GetOrigin();</div><div>  FixedImageType::RegionType FixedRegionhigh = FixedImage-&gt;GetBufferedRegion();</div>

<div>  FixedImageType::SpacingType   fixedSpacinghigh    = FixedImage-&gt;GetSpacing();</div><div>  FixedImageType::SizeType FixedImageSizehigh = FixedRegionhigh.GetSize();</div><div> </div><div>      for(unsigned int r=0; r&lt;Dimension; r++)</div>

<div>    {</div><div>    </div><div>    spacinghigh[r] = floor( fixedSpacinghigh[r] * (FixedImageSizehigh[r] - 1) / (gridSizeOnImagehigh[r] - 1) ); //new way</div><div>    </div><div>    }</div><div>   </div><div>   SpacingTypehigh gridOriginOffsethigh = gridDirection * spacinghigh; </div>

<div>   OriginTypehigh gridOriginhigh = originhigh - gridOriginOffsethigh;</div><div> </div><div>  BtransformHigh-&gt;SetGridSpacing( spacinghigh );</div><div>  BtransformHigh-&gt;SetGridOrigin( gridOriginhigh );</div><div>

  BtransformHigh-&gt;SetGridRegion( bsplineRegionhigh ); //i.e region of gridsize</div><div>  BtransformHigh-&gt;SetGridDirection( gridDirection );</div><div>  </div><div>  std::cout &lt;&lt; &quot;grid region for high Bspline is &quot; &lt;&lt; BtransformHigh-&gt;GetGridRegion() &lt;&lt;std::endl;</div>

<div>  </div><div>///////////////////////////////////////////////////////////////////////</div><div>///METRIC NUMBER OF SAMPLES</div><div>//////////////////////////////////////////////////////////////////////</div><div><br>

</div><div><br></div><div>//Low</div><div>  const unsigned int numberOfBSplineParametersLow =</div><div>               BtransformLow-&gt;GetNumberOfParameters();</div><div>  </div><div>  typedef DeformableTransformType::ParametersType     ParametersTypeLow;</div>

<div>  ParametersTypeLow parametersLow( numberOfBSplineParametersLow );</div><div>  </div><div>//Medium</div><div><br></div><div>  const unsigned int numberOfBSplineParametersMed =</div><div>               BtransformMed-&gt;GetNumberOfParameters();</div>

<div>  </div><div>  typedef DeformableTransformType::ParametersType     ParametersTypeMed;</div><div>  ParametersTypeMed parametersMed( numberOfBSplineParametersMed );</div><div><br></div><div>  const unsigned int numberOfBSplineParametersHigh =</div>

<div>               BtransformHigh-&gt;GetNumberOfParameters();</div><div>  </div><div>  typedef DeformableTransformType::ParametersType     ParametersTypeHigh;</div><div>  ParametersTypeHigh parametersHigh( numberOfBSplineParametersHigh );</div>

<div> </div><div><br></div><div>     const unsigned long numberOfSamplesBsplineHigh =</div><div>       static_cast&lt;unsigned long&gt;(</div><div>        vcl_sqrt( static_cast&lt;double&gt;( numberOfBSplineParametersHigh ) *</div>

<div>                static_cast&lt;double&gt;( numberOfPixels ) ) );</div><div><br></div><div><br></div><div><br></div><div>///////////////////////////////////////////////////////////////////////</div><div>///LOW PARAMETERS</div>

<div>//////////////////////////////////////////////////////////////////////</div><div><br></div><div>  typedef DeformableTransformType::ParametersType     ParametersTypelow;</div><div><br></div><div>  const unsigned int numberOfBSplineParameterslow =</div>

<div>               BtransformLow-&gt;GetNumberOfParameters();</div><div><br></div><div>  ParametersTypelow parameterslow( numberOfBSplineParameterslow );</div><div>  parameterslow.Fill( 0.0 );</div><div>  BtransformLow-&gt;SetBulkTransform( Affinetransform );</div>

<div> </div><div> //////////////////////////////////////////////////////////////////////////////</div><div>///OPTIMIZER FOR LOW RESOLUTION</div><div>////////////////////////////////////////////////////////////////////////////////</div>

<div><br></div><div> // if get IFLAG= -1. LINE SEARCH FAILED</div><div> // increase step lentgth but keep other parameters the same</div><div><br></div><div>/////////////////////////////////////////////////////////////////////</div>

<div>// SEE HERE FOR OPTIMIZER ADVICE</div><div>//<a href="http://www.itk.org/pipermail/insight-users/2011-March/040286.html" target="_blank">http://www.itk.org/pipermail/insight-users/2011-March/040286.html</a></div><div>
<br></div><div>
///FOR the LBFGSBOptimizer got LOW</div><div><br></div><div>  bOptimizerType::BoundSelectionType boundSelectlow( BtransformLow-&gt;GetNumberOfParameters() );</div><div>  bOptimizerType::BoundValueType upperBoundlow( BtransformLow-&gt;GetNumberOfParameters() );</div>

<div>  bOptimizerType::BoundValueType lowerBoundlow( BtransformLow-&gt;GetNumberOfParameters() );</div><div>  boundSelectlow.Fill( 0 );</div><div>  upperBoundlow.Fill( 0.0 );</div><div>  lowerBoundlow.Fill( 0.0 );</div><div>

  boptimizer-&gt;SetBoundSelection( boundSelectlow );</div><div>  boptimizer-&gt;SetUpperBound( upperBoundlow );</div><div>  boptimizer-&gt;SetLowerBound( lowerBoundlow );</div><div>  boptimizer-&gt;SetCostFunctionConvergenceFactor( 1.e1 );</div>

<div>  boptimizer-&gt;SetProjectedGradientTolerance( 1e-6 );  </div><div>  boptimizer-&gt;SetMaximumNumberOfIterations( 200 ); //was 200</div><div>  boptimizer-&gt;SetMaximumNumberOfEvaluations( 300 ); //</div><div>  boptimizer-&gt;SetMaximumNumberOfCorrections( 5 );</div>

<div>///////////////////////////////////////////////////////////////</div><div><br></div><div>  metric-&gt;SetNumberOfSpatialSamples( numberOfPixels*1 ); //</div><div>  RegistrationType::Pointer   bregistration  = RegistrationType::New();</div>

<div>  bregistration-&gt;SetNumberOfThreads(1);</div><div>  bregistration-&gt;SetInitialTransformParameters( BtransformLow-&gt;GetParameters() );</div><div>  bregistration-&gt;SetTransform(BtransformLow);</div><div>  bregistration-&gt;SetFixedImageRegion(FixedImage-&gt;GetBufferedRegion() );</div>

<div>  bregistration-&gt;SetMetric(        metric        );</div><div>  bregistration-&gt;SetOptimizer(     boptimizer     );</div><div>  bregistration-&gt;SetInterpolator(  interpolator  );</div><div>  bregistration-&gt;SetFixedImage(fixedSmoother-&gt;GetOutput());</div>

<div>  bregistration-&gt;SetMovingImage(movingSmoother-&gt;GetOutput());</div><div>  bregistration-&gt;SetNumberOfThreads(1);</div><div><br></div><div><br></div><div>//////////////////////////////////////////////////////////////////////</div>

<div>///REGISTRATION FOR LOW</div><div>///////////////////////////////////////////////////////////////////////</div><div><br></div><div>  BCommandIterationUpdate::Pointer observerl = BCommandIterationUpdate::New();</div>
<div>
  boptimizer-&gt;AddObserver( itk::IterationEvent(), observerl );</div><div>  try</div><div>    {</div><div><br></div><div><span style="white-space:pre-wrap">        </span></div><div>    bregistration-&gt;Update();</div><div>
    std::cout &lt;&lt; &quot; &quot; &lt;&lt;std::endl;</div><div>    std::cout &lt;&lt; &quot;Optimizer stop condition: &quot;</div><div>           &lt;&lt; bregistration-&gt;GetOptimizer()-&gt;GetStopConditionDescription()</div>

<div>           &lt;&lt; std::endl;</div><div><br></div><div>    }</div><div>  catch( itk::ExceptionObject &amp; err )</div><div>    {</div><div>    std::cerr &lt;&lt; &quot;ExceptionObject caught !&quot; &lt;&lt; std::endl;</div>

<div>    std::cerr &lt;&lt; err &lt;&lt; std::endl;</div><div>    return EXIT_FAILURE;</div><div>    }</div><div><br></div><div> boptimizer-&gt;RemoveAllObservers();</div><div><br></div><div><br></div><div><br></div><div>

/////////////////////////////////////////////////////////////////////</div><div> ///MEDIUM PARAMETERS </div><div>////////////////////////////////////////////////////////////////////  </div><div>  </div><div>  typedef DeformableTransformType::ParametersType     ParametersTypeMed;</div>

<div>  parametersMed.Fill( 0.0 ); </div><div><br></div><div>//  Now we need to initialize the BSpline coefficients of the higher resolution</div><div>//  transform. This is done by first computing the actual deformation field </div>

<div>//  at the higher resolution from the lower resolution BSpline coefficients. </div><div>//  Then a BSpline decomposition is done to obtain the BSpline coefficient of </div><div>//  the higher resolution transform.</div>

<div> </div><div>  unsigned int counterMed = 0;</div><div><br></div><div>  for ( unsigned int k = 0; k &lt; SpaceDimension; k++ )</div><div>    {</div><div>    typedef DeformableTransformType::ImageType ParametersImageTypeMed;</div>

<div>    typedef itk::ResampleImageFilter&lt;ParametersImageTypeMed,ParametersImageTypeMed&gt; ResamplerTypeMed;</div><div>    ResamplerTypeMed::Pointer upsamplerMed = ResamplerTypeMed::New();</div><div><br></div><div>    typedef itk::BSplineResampleImageFunction&lt;ParametersImageTypeMed,double&gt; FunctionTypeMed;</div>

<div>    FunctionTypeMed::Pointer functionMed = FunctionTypeMed::New();</div><div><br></div><div>    typedef itk::IdentityTransform&lt;double,SpaceDimension&gt; IdentityTransformTypeMed;</div><div>    IdentityTransformTypeMed::Pointer identityMed = IdentityTransformTypeMed::New();</div>

<div><br></div><div>    upsamplerMed-&gt;SetInput( BtransformLow-&gt;GetCoefficientImages()[k] );</div><div>    upsamplerMed-&gt;SetInterpolator( functionMed );</div><div>    upsamplerMed-&gt;SetTransform( identityMed );</div>

<div>    upsamplerMed-&gt;SetSize( BtransformMed-&gt;GetGridRegion().GetSize() );</div><div>    upsamplerMed-&gt;SetOutputSpacing( BtransformMed-&gt;GetGridSpacing() );</div><div>    upsamplerMed-&gt;SetOutputOrigin( BtransformMed-&gt;GetGridOrigin() );</div>

<div>    upsamplerMed-&gt;SetOutputDirection( FixedImage-&gt;GetDirection() );</div><div>    upsamplerMed-&gt;Update();</div><div><br></div><div>    typedef itk::BSplineDecompositionImageFilter&lt;ParametersImageTypeMed,ParametersImageTypeMed&gt;</div>

<div>      DecompositionTypeMed;</div><div>    DecompositionTypeMed::Pointer decompositionMed = DecompositionTypeMed::New();</div><div><br></div><div>    decompositionMed-&gt;SetSplineOrder( SplineOrder );</div><div>    decompositionMed-&gt;SetInput( upsamplerMed-&gt;GetOutput() );</div>

<div>    decompositionMed-&gt;Update();</div><div><br></div><div>    ParametersImageTypeMed::Pointer newCoefficientsMed = decompositionMed-&gt;GetOutput();</div><div><br></div><div>    // copy the coefficients into the parameter array</div>

<div>    typedef itk::ImageRegionIterator&lt;ParametersImageTypeMed&gt; IteratorMed;</div><div>    IteratorMed it( newCoefficientsMed, BtransformMed-&gt;GetGridRegion() );</div><div>    while ( !it.IsAtEnd() )</div><div>
      {</div>
<div>      parametersMed[ counterMed++ ] = it.Get();</div><div>      ++it;</div><div>      }</div><div><br></div><div>    }</div><div>  </div><div>  BtransformMed-&gt;SetParameters( parametersMed );</div><div><br></div><div>

////////////////////////////////////////////////////////////////</div><div><br></div><div>///FOR the LBFGSBOptimizer MEDIUM</div><div><br></div><div>  bOptimizerType::BoundSelectionType boundSelectMed( BtransformMed-&gt;GetNumberOfParameters() );</div>

<div>  bOptimizerType::BoundValueType upperBoundMed( BtransformMed-&gt;GetNumberOfParameters() );</div><div>  bOptimizerType::BoundValueType lowerBoundMed( BtransformMed-&gt;GetNumberOfParameters() );</div><div>  boundSelectMed.Fill( 0 );</div>

<div>  upperBoundMed.Fill( 0.0 );</div><div>  lowerBoundMed.Fill( 0.0 );</div><div>  boptimizer-&gt;SetBoundSelection( boundSelectMed );</div><div>  boptimizer-&gt;SetUpperBound( upperBoundMed );</div><div>  boptimizer-&gt;SetLowerBound( lowerBoundMed );</div>

<div>  boptimizer-&gt;SetCostFunctionConvergenceFactor( 1.e12 ); //was 1.e7</div><div> //set/Get the CostFunctionConvergenceFactor. Algorithm terminates</div><div> //  when the reduction in cost function is less than factor * epsmcj</div>

<div> // where epsmch is the machine precision.</div><div> // Typical values for factor: 1e+12 for low accuracy;</div><div> // 1e+7 for moderate accuracy and 1e+1 for extremely high accuracy.</div><div>  boptimizer-&gt;SetProjectedGradientTolerance( 1e-6 );</div>

<div>  //Gradient tolerance determines a lower threshold which cuases the optimization to stop if its</div><div>//value reached, recommended values for gradient tolerance are 1e-6 to 1e-10, higher values might result in premature elimination of the optimization process</div>

<div>  boptimizer-&gt;SetMaximumNumberOfIterations( 200); //was 200</div><div>  boptimizer-&gt;SetMaximumNumberOfEvaluations( 300 );</div><div>  boptimizer-&gt;SetMaximumNumberOfCorrections( 5 );</div><div>  </div><div>  //M Number of correction- number of function/gradient pairs used to build  model</div>

<div>  //This parameter is passed to the function which is used to create optimizer object.</div><div>  // Increase of number of corrections decreases number of function evaluations </div><div>  //and iterations needed to converge, but increases computational overhead associated with iteration.</div>

<div>  // On well-conditioned problems can be as small as 3-10. </div><div>  //If function changes rapidly in some directions and slowly in other ones, </div><div>  //then you can try increasing M in order to increase convergenc</div>

<div>  </div><div>///////////////////////////////////////////////////////////////////////</div><div>///REGISTRATION FOR MEDIUM</div><div>/////////////////////////////////////////////////////////////////////</div><div>  //For time probe</div>

<div>  </div><div>   metric-&gt;SetNumberOfSpatialSamples( numberOfPixels*1 ); //</div><div>  bregistration-&gt;SetInitialTransformParameters( BtransformMed-&gt;GetParameters() );</div><div>  bregistration-&gt;SetTransform( BtransformMed );</div>

<div>  </div><div>  itk::TimeProbesCollectorBase chronometerm;</div><div>  itk::MemoryProbesCollectorBase memorymeterm;</div><div><br></div><div>  BCommandIterationUpdate::Pointer observerm = BCommandIterationUpdate::New();</div>

<div>  boptimizer-&gt;AddObserver( itk::IterationEvent(), observerm );</div><div><br></div><div><br></div><div>  try </div><div>    { </div><div>    </div><div>    bregistration-&gt;Update(); </div><div>    </div><div>    </div>

<div>    } </div><div>  catch( itk::ExceptionObject &amp; err ) </div><div>    { </div><div>    std::cerr &lt;&lt; &quot;ExceptionObject caught !&quot; &lt;&lt; std::endl; </div><div>    std::cerr &lt;&lt; err &lt;&lt; std::endl; </div>

<div>    return EXIT_FAILURE;</div><div>    } </div><div><br></div><div>  boptimizer-&gt;RemoveAllObservers();</div><div>  </div><div> BtransformMed-&gt;SetParameters( bregistration-&gt;GetLastTransformParameters() );</div>

<div><br></div><div>//NOW HIGH RESOLUTION</div><div>/////////////////////////////////////////////////////////////////////</div><div> ///HIGH PARAMETERS </div><div>////////////////////////////////////////////////////////////////////  </div>

<div>  </div><div>  parametersHigh.Fill( 0.0 ); </div><div><br></div><div>  unsigned int counter = 0;</div><div><br></div><div>  for ( unsigned int k = 0; k &lt; SpaceDimension; k++ )</div><div>    {</div><div>    typedef DeformableTransformType::ImageType ParametersImageType;</div>

<div>    typedef itk::ResampleImageFilter&lt;ParametersImageType,ParametersImageType&gt; ResamplerType;</div><div>    ResamplerType::Pointer upsampler = ResamplerType::New();</div><div><br></div><div>    typedef itk::BSplineResampleImageFunction&lt;ParametersImageType,double&gt; FunctionType;</div>

<div>    FunctionType::Pointer function = FunctionType::New();</div><div><br></div><div>    typedef itk::IdentityTransform&lt;double,SpaceDimension&gt; IdentityTransformType;</div><div>    IdentityTransformType::Pointer identity = IdentityTransformType::New();</div>

<div><br></div><div>    upsampler-&gt;SetInput( BtransformMed-&gt;GetCoefficientImages()[k] );</div><div>    upsampler-&gt;SetInterpolator( function );</div><div>    upsampler-&gt;SetTransform( identity );</div><div>    upsampler-&gt;SetSize( BtransformHigh-&gt;GetGridRegion().GetSize() );</div>

<div>    upsampler-&gt;SetOutputSpacing( BtransformHigh-&gt;GetGridSpacing() );</div><div>    upsampler-&gt;SetOutputOrigin( BtransformHigh-&gt;GetGridOrigin() );</div><div>    upsampler-&gt;SetOutputDirection( FixedImage-&gt;GetDirection() );</div>

<div>    upsampler-&gt;Update();</div><div><br></div><div><br></div><div>    typedef itk::BSplineDecompositionImageFilter&lt;ParametersImageType,ParametersImageType&gt;</div><div>      DecompositionType;</div><div>    DecompositionType::Pointer decomposition = DecompositionType::New();</div>

<div><br></div><div>    decomposition-&gt;SetSplineOrder( SplineOrder );</div><div>    decomposition-&gt;SetInput( upsampler-&gt;GetOutput() );</div><div>    decomposition-&gt;Update();</div><div><br></div><div>    ParametersImageType::Pointer newCoefficients = decomposition-&gt;GetOutput();</div>

<div><br></div><div>    // copy the coefficients into the parameter array</div><div>    typedef itk::ImageRegionIterator&lt;ParametersImageType&gt; Iterator;</div><div>    Iterator it( newCoefficients, BtransformHigh-&gt;GetGridRegion() );</div>

<div>    while ( !it.IsAtEnd() )</div><div>      {</div><div>      parametersHigh[ counter++ ] = it.Get();</div><div>      ++it;</div><div>      }</div><div><br></div><div>    }</div><div>  </div><div>  BtransformHigh-&gt;SetParameters( parametersHigh );</div>

<div><br></div><div><br></div><div>/////////////////////////////////////////////////////////////////////////////////////////</div><div><br></div><div>///FOR the LBFGSBOptimizer HIGH</div><div><br></div><div>  bOptimizerType::BoundSelectionType boundSelecthigh( BtransformHigh-&gt;GetNumberOfParameters() );</div>

<div>  bOptimizerType::BoundValueType upperBoundhigh( BtransformHigh-&gt;GetNumberOfParameters() );</div><div>  bOptimizerType::BoundValueType lowerBoundhigh( BtransformHigh-&gt;GetNumberOfParameters() );</div><div>  boundSelecthigh.Fill( 0 );</div>

<div>  upperBoundhigh.Fill( 0.0 );</div><div>  lowerBoundhigh.Fill( 0.0 );</div><div>  boptimizer-&gt;SetBoundSelection( boundSelecthigh );</div><div>  boptimizer-&gt;SetUpperBound( upperBoundhigh );</div><div>  boptimizer-&gt;SetLowerBound( lowerBoundhigh );</div>

<div>  boptimizer-&gt;SetCostFunctionConvergenceFactor( 1.e12 );</div><div>  boptimizer-&gt;SetProjectedGradientTolerance( 1e-6 );//was 1e-6</div><div>  boptimizer-&gt;SetMaximumNumberOfIterations( 1000 ); //was 200</div>

<div>  boptimizer-&gt;SetMaximumNumberOfEvaluations( 300 );</div><div>  boptimizer-&gt;SetMaximumNumberOfCorrections( 5 );</div><div>  </div><div>  </div><div>  </div><div>///////////////////////////////////////////////////////////////////////</div>

<div>///REGISTRATION FOR HIGH</div><div>/////////////////////////////////////////////////////////////////////</div><div>  //For time probe</div><div>  </div><div>   metric-&gt;SetNumberOfSpatialSamples( numberOfPixels*1 ); //</div>

<div>  bregistration-&gt;SetInitialTransformParameters( BtransformHigh-&gt;GetParameters() );</div><div>  bregistration-&gt;SetTransform( BtransformHigh );</div><div>  </div><div>  itk::TimeProbesCollectorBase chronometerh;</div>

<div>  itk::MemoryProbesCollectorBase memorymeterh;</div><div><br></div><div>  BCommandIterationUpdate::Pointer observerh = BCommandIterationUpdate::New();</div><div>  boptimizer-&gt;AddObserver( itk::IterationEvent(), observerh );</div>

<div>  </div><div> try </div><div>    { </div><div> </div><div>    bregistration-&gt;Update(); </div><div>    </div><div>   std::cout &lt;&lt; &quot; &quot; &lt;&lt;std::endl;</div><div>    std::cout &lt;&lt; &quot;Optimizer stop condition: &quot;</div>

<div>           &lt;&lt; bregistration-&gt;GetOptimizer()-&gt;GetStopConditionDescription()</div><div>           &lt;&lt; std::endl;</div><div>    } </div><div>  catch( itk::ExceptionObject &amp; err ) </div><div>    { </div>

<div>    std::cerr &lt;&lt; &quot;ExceptionObject caught !&quot; &lt;&lt; std::endl; </div><div>    std::cerr &lt;&lt; err &lt;&lt; std::endl; </div><div>    return EXIT_FAILURE;</div><div>    } </div><div><br></div><div>

  // Finally we use the last transform parameters in order to resample the image.</div><div>  //</div><div> </div><div>  boptimizer-&gt;RemoveAllObservers();</div><div> </div><div> BtransformHigh-&gt;SetParameters( bregistration-&gt;GetLastTransformParameters() );</div>

<div><br></div><div><br></div><div>  typedef itk::ResampleImageFilter&lt;</div><div>                            MovingImageType,</div><div>                            FixedImageType &gt;    ResampleFilterType;</div><div>
<br>
</div><div>  ResampleFilterType::Pointer resample = ResampleFilterType::New();</div><div><br></div><div><br></div><div>  resample-&gt;SetTransform( BtransformHigh );</div><div>  resample-&gt;SetInput( MovedImage );</div>
<div>
  resample-&gt;SetSize(    FixedImage-&gt;GetLargestPossibleRegion().GetSize() );</div><div>  resample-&gt;SetOutputOrigin(  FixedImage-&gt;GetOrigin() );</div><div>  resample-&gt;SetOutputSpacing( FixedImage-&gt;GetSpacing() );</div>

<div>  resample-&gt;SetOutputDirection( FixedImage-&gt;GetDirection() );</div><div>  resample-&gt;SetDefaultPixelValue( 0 );</div><div>  resample-&gt;Update();</div><div><br></div><div>  ImageWriter-&gt;SetInput(resample-&gt;GetOutput());</div>

<div><br></div><div>    try</div><div>      {</div><div>      ImageWriter-&gt;Update();</div><div>      }</div><div>    catch( itk::ExceptionObject &amp; excp )</div><div>      {</div><div>      std::cerr &lt;&lt; &quot;Exception thrown &quot; &lt;&lt; std::endl;</div>

<div>      std::cerr &lt;&lt; excp &lt;&lt; std::endl;</div><div>      return EXIT_FAILURE;</div><div>      }</div><div><br></div><div>    </div><div>  FixedImageType::RegionType FixedRegion = FixedImage-&gt;GetBufferedRegion();</div>

<div>  DeformField-&gt;SetRegions( FixedRegion );</div><div>  DeformField-&gt;SetOrigin( FixedImage-&gt;GetOrigin() );</div><div>  DeformField-&gt;SetSpacing( FixedImage-&gt;GetSpacing() );</div><div>  FixedImageType:: DirectionType Fixeddirection =FixedImage-&gt;GetDirection();</div>

<div>  DeformField-&gt;SetDirection(Fixeddirection);</div><div>  DeformField-&gt;Allocate();</div><div>  DeformField-&gt;Update();</div><div>  </div><div>//WRITE THE DEFORMATION FIELD</div><div><br></div><div>    typedef itk::ImageRegionIterator&lt; DisplacementFieldType &gt; FieldIterator;</div>

<div>    FieldIterator fi( DeformField, FixedRegion );</div><div><br></div><div>    fi.GoToBegin();</div><div><br></div><div>    DeformableTransformType::InputPointType  fixedPointa;</div><div>    DeformableTransformType::OutputPointType movingPointa;</div>

<div>    DisplacementFieldType::IndexType indexa;</div><div><br></div><div>    VectorPixelType displacement;</div><div><br></div><div>    while( ! fi.IsAtEnd() )</div><div>      {</div><div>      indexa = fi.GetIndex();</div>

<div>      DeformField-&gt;TransformIndexToPhysicalPoint( indexa, fixedPointa );</div><div>      movingPointa = BtransformHigh-&gt;TransformPoint( fixedPointa );</div><div>      displacement = movingPointa - fixedPointa;</div>

<div>      fi.Set( displacement );</div><div>      ++fi;</div><div>      }</div><div><br></div><div>    FieldWriter-&gt;SetInput( DeformField );</div><div>    //std::cout &lt;&lt; &quot;Writing deformation field ...&quot;;</div>

<div><br></div><div>    try</div><div>      {</div><div>      FieldWriter-&gt;Update();</div><div>      }</div><div>    catch( itk::ExceptionObject &amp; excp )</div><div>      {</div><div>      std::cerr &lt;&lt; &quot;Exception thrown &quot; &lt;&lt; std::endl;</div>

<div>      std::cerr &lt;&lt; excp &lt;&lt; std::endl;</div><div>      return EXIT_FAILURE;</div><div>      }</div><div><br></div><div>    </div><div>  return EXIT_SUCCESS;</div><div>}</div><div> </div></div><div><br>
</div><div><br></div></div>
<br>_____________________________________<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.php" target="_blank">http://www.kitware.com/products/protraining.php</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></blockquote></div><br></div>