<html><head></head><body style="word-wrap: break-word; -webkit-nbsp-mode: space; -webkit-line-break: after-white-space; "><div>Hi,</div><div><br></div><div>I'm trying to use a simple metric like MeanSquaresImageToImageMetric with a BsplineInterpolatorImageFunction and RegularStepGradientDescendetOptimizer to solve a particular registration problem but the optimizer can achieve the convergence (any convergence) and the result destroy the image. I used ITK-4.2.0 and ITK-4.1 with the same problem.</div><div><br></div><div><br></div><div>In my code and in the past I've used the optimizer of LBFGSBoptimizer but I've never got a stable convergence and change this for the RegularStepGradient (any help with LBFGSoptimizer is accepted too).</div><div><br></div><div>If anyone have an example how to use MeansSquares metric, Bspline interpolation and RegularStepGradient with the successfully result please send me the code.</div><div><br></div><div><br></div><div>I use the code with following parameters and two image of 200x200:</div><div><br></div><div><div>--numberTrheads 1</div><div>--numberLevels 5</div><div>--space 16</div><div>--metricType SSD</div><div>--optimizerType RSGD</div><div>--interpolatorType BSI</div><div>--maxIt 500</div><div>--relaxFact 0.8</div><div>--steplength 0.001 2</div><div><br></div></div><div><br></div><div><br></div><div>An extract of my code consern to the RESGOptimizer, SSDMetric and BsplineInterpolator is:</div><div><br></div><div><br></div><div><div>typedef itk::SingleValuedNonLinearOptimizer<span class="Apple-tab-span" style="white-space:pre"> </span>OptimizerType;//Abstract Optimizer</div><div>typedef itk::RegularStepGradientDescentOptimizer <span class="Apple-tab-span" style="white-space:pre"> </span>RSGDOptimizerType;</div><div><span class="Apple-tab-span" style="white-space:pre"> </span></div><div><br></div><div>typedef itk::ImageToImageMetric<FixedImageType,MovingImageType > MetricType;//Abstract Metric</div><div>typedef itk::MeanSquaresImageToImageMetric<FixedImageType,MovingImageType > SSDMetricType;</div><div><br></div><div>typedef itk::InterpolateImageFunction< MovingImageType,double> InterpolatorType;</div><div>typedef itk::BSplineInterpolateImageFunction<MovingImageType,double, double><span class="Apple-tab-span" style="white-space:pre"> </span> BsplineInterpolatorType;</div><div><br></div><div>typedef itk::MultiResolutionImageRegistrationMethod<FixedImageType,MovingImageType > MultiRegistrationType;</div><div><br></div><div><br></div><div>typedef itk::MultiResolutionPyramidImageFilter<FixedImageType,MovingImageType > FixedImagePyramidType;</div><div>typedef itk::MultiResolutionPyramidImageFilter<FixedImageType,MovingImageType > MovingImagePyramidType;</div><div><br></div><div>typedef itk::ResampleImageFilter<FixedImageType,MovingImageType > ResampleFilterType;</div><div><br></div><div><br></div><div><br></div><div>typename MetricType::Pointer metric = NULL;</div><div>typename InterpolatorType::Pointer interpolator = NULL;</div><div>typename OptimizerType::Pointer optimizer = NULL;</div><div><br></div><div>optimizer = RSGDOptimizerType::New();</div><div>metric = SSDMetricType::New();</div><div>interpolator = BsplineInterpolatorType::New();</div><div><br></div><div><br></div><div>typedef BsplineInterpolatorType * InterpolatorPointer;</div><div>InterpolatorPointer int_aux = dynamic_cast< InterpolatorPointer >(interpolator.GetPointer());</div><div>int_aux->SetSplineOrder(3);</div><div><br></div><div><br></div><div><br></div><div>typename TransformType::Pointer transform = TransformType::New();</div><div>typename MultiRegistrationType::Pointer multiregistration = MultiRegistrationType::New();</div><div><br></div><div>multiregistration->SetNumberOfThreads(this->m_NumberOfThreads);</div><div>typename FixedImagePyramidType::Pointer fixedImagePyramid = FixedImagePyramidType::New();</div><div>typename MovingImagePyramidType::Pointer movingImagePyramid = MovingImagePyramidType::New();</div><div><br></div><div>typename FixedImageType::RegionType fixedRegion = fixedImage->GetBufferedRegion();</div><div><br></div><div>multiregistration->SetMetric( metrics[i] );</div><div>multiregistration->SetOptimizer( optimizer );</div><div>multiregistration->SetInterpolator( interpolator );</div><div>multiregistration->SetTransform( transform );</div><div><br></div><div>// Setup the registration</div><div>multiregistration->SetFixedImagePyramid( fixedImagePyramid );</div><div>multiregistration->SetMovingImagePyramid( movingImagePyramid );</div><div>multiregistration->SetFixedImage( fixedImage );</div><div>multiregistration->SetMovingImage( movingImage);</div><div>multiregistration->SetFixedImageRegion( fixedRegion );</div><div><br></div><div><br></div><div><br></div><div>// Here we define the parameters of the BSplineDeformableTransform grid. We</div><div>// arbitrarily decide to use a grid with $5 \times 5$ nodes within the image.</div><div>// The reader should note that the BSpline computation requires a</div><div>// finite support region ( 1 grid node at the lower borders and 2</div><div>// grid nodes at upper borders). Therefore in this example, we set</div><div>// the grid size to be $8 \times 8$ and place the grid origin such that</div><div>// grid node (1,1) coincides with the first pixel in the fixed image.</div><div><br></div><div>unsigned int numberOfGridNodesInOneDimension = 1 + fixedImage->GetLargestPossibleRegion().GetSize()[0] / space;</div><div><br></div><div><br></div><div><br></div><div>#if ITK_VERSION_MAJOR < 4</div><div><br></div><div><br></div><div><span class="Apple-tab-span" style="white-space:pre"> </span>typename TransformType::RegionType bsplineRegion;</div><div><span class="Apple-tab-span" style="white-space:pre"> </span>typename TransformType::RegionType::SizeType gridSizeOnImage;</div><div><span class="Apple-tab-span" style="white-space:pre"> </span>typename TransformType::RegionType::SizeType gridBorderSize;</div><div><span class="Apple-tab-span" style="white-space:pre"> </span>typename TransformType::RegionType::SizeType totalGridSize;</div><div><br></div><div><span class="Apple-tab-span" style="white-space:pre"> </span>gridSizeOnImage.Fill( numberOfGridNodesInOneDimension );</div><div><span class="Apple-tab-span" style="white-space:pre"> </span>gridBorderSize.Fill( SplineOrder ); // Border for spline order = 3 ( 1 lower, 2 upper )</div><div><span class="Apple-tab-span" style="white-space:pre"> </span>totalGridSize = gridSizeOnImage + gridBorderSize;</div><div><br></div><div><span class="Apple-tab-span" style="white-space:pre"> </span>bsplineRegion.SetSize( totalGridSize );</div><div><br></div><div><br></div><div><span class="Apple-tab-span" style="white-space:pre"> </span>typename TransformType::SpacingType spacing = fixedImage->GetSpacing();</div><div><br></div><div><br></div><div><span class="Apple-tab-span" style="white-space:pre"> </span>typename TransformType::OriginType origin = fixedImage->GetOrigin();</div><div><br></div><div><span class="Apple-tab-span" style="white-space:pre"> </span>typename FixedImageType::SizeType fixedImageSize = fixedRegion.GetSize();</div><div><br></div><div><span class="Apple-tab-span" style="white-space:pre"> </span>for(unsigned int r=0; r<Dimension; r++)</div><div><span class="Apple-tab-span" style="white-space:pre"> </span>{</div><div><span class="Apple-tab-span" style="white-space:pre"> </span>spacing[r] *= static_cast<double>(fixedImageSize[r] - 1) /</div><div><span class="Apple-tab-span" style="white-space:pre"> </span> static_cast<double>(gridSizeOnImage[r] - 1);</div><div><span class="Apple-tab-span" style="white-space:pre"> </span>}</div><div><br></div><div><span class="Apple-tab-span" style="white-space:pre"> </span>typename FixedImageType::DirectionType gridDirection = fixedImage->GetDirection();</div><div><span class="Apple-tab-span" style="white-space:pre"> </span>typename TransformType::SpacingType gridOriginOffset = gridDirection * spacing;</div><div><br></div><div><span class="Apple-tab-span" style="white-space:pre"> </span>typename TransformType::OriginType gridOrigin = origin - gridOriginOffset;</div><div><br></div><div><span class="Apple-tab-span" style="white-space:pre"> </span>transform->SetGridSpacing( spacing );</div><div><span class="Apple-tab-span" style="white-space:pre"> </span>transform->SetGridOrigin( gridOrigin );</div><div><span class="Apple-tab-span" style="white-space:pre"> </span>transform->SetGridRegion( bsplineRegion );</div><div><span class="Apple-tab-span" style="white-space:pre"> </span>transform->SetGridDirection( gridDirection );</div><div><br></div><div>#else</div><div><span class="Apple-tab-span" style="white-space:pre"> </span>typename TransformType::PhysicalDimensionsType fixedPhysicalDimensions;</div><div><span class="Apple-tab-span" style="white-space:pre"> </span>typename TransformType::MeshSizeType meshSize;</div><div><span class="Apple-tab-span" style="white-space:pre"> </span>typename TransformType::OriginType fixedOrigin;</div><div><br></div><div><span class="Apple-tab-span" style="white-space:pre"> </span>for( unsigned int k=0; k< Dimension; k++ )</div><div><span class="Apple-tab-span" style="white-space:pre"> </span>{</div><div><span class="Apple-tab-span" style="white-space:pre"> </span>fixedOrigin = fixedImage->GetOrigin()[k];</div><div><span class="Apple-tab-span" style="white-space:pre"> </span>fixedPhysicalDimensions[k] = fixedImage->GetSpacing()[k] *</div><div><span class="Apple-tab-span" style="white-space:pre"> </span> static_cast<double>(</div><div><span class="Apple-tab-span" style="white-space:pre"> </span> fixedImage->GetLargestPossibleRegion().GetSize()[k] - 1 );</div><div><span class="Apple-tab-span" style="white-space:pre"> </span>}</div><div><span class="Apple-tab-span" style="white-space:pre"> </span>meshSize.Fill( numberOfGridNodesInOneDimension - SplineOrder );</div><div><br></div><div><span class="Apple-tab-span" style="white-space:pre"> </span>transform->SetTransformDomainOrigin( fixedOrigin );</div><div><span class="Apple-tab-span" style="white-space:pre"> </span>transform->SetTransformDomainPhysicalDimensions( fixedPhysicalDimensions );</div><div><span class="Apple-tab-span" style="white-space:pre"> </span>transform->SetTransformDomainMeshSize( meshSize );</div><div><span class="Apple-tab-span" style="white-space:pre"> </span>transform->SetTransformDomainDirection( fixedImage->GetDirection() );</div><div><br></div><div>#endif</div><div><br></div><div><span class="Apple-tab-span" style="white-space:pre"> </span>const unsigned int numberOfParameters = transform->GetNumberOfParameters();</div><div><br></div><div><span class="Apple-tab-span" style="white-space:pre"> </span>typename TransformType::ParametersType parameters( numberOfParameters );</div><div><br></div><div><span class="Apple-tab-span" style="white-space:pre"> </span>parameters.Fill( 0.0 );</div><div><br></div><div><span class="Apple-tab-span" style="white-space:pre"> </span>transform->SetParameters( parameters );</div><div><br></div><div><span class="Apple-tab-span" style="white-space:pre"> </span>// We now pass the parameters of the current transform as the initial</div><div><span class="Apple-tab-span" style="white-space:pre"> </span>// parameters to be used when the registration process starts.</div><div><br></div><div><br></div><div>multiregistration->SetInitialTransformParameters( transform->GetParameters() );</div><div><br></div><div><br></div><div>multiregistration->SetNumberOfLevels( numberLevels );</div><div><br></div><div><br></div><div>// Set the parameters of the RegularStepGradientDescentOptimizer Optimizer.</div><div>typedef RSGDOptimizerType * OptimizerPointer;</div><div><br></div><div>OptimizerPointer opt = dynamic_cast< OptimizerPointer >(optimizer.GetPointer());</div><div>opt->SetNumberOfIterations( maxIt );</div><div>opt->SetRelaxationFactor( relaxFact );</div><div>opt->SetMaximumStepLength( steplength[1] );</div><div>opt->SetMinimumStepLength( steplength[0] );</div><div><br></div><div><br></div><div><br></div><div>std::cout << "Initial Parameters = " << std::endl;</div><div>std::cout << transform->GetParameters() << std::endl;</div><div><br></div><div><br></div><div>CommandIterationUpdate::Pointer observer = CommandIterationUpdate::New();</div><div>optimizer->AddObserver( itk::IterationEvent(), observer );</div><div><br></div><div><br></div><div>std::string metric_name = metric->GetNameOfClass();</div><div><br></div><div>std::cout << std::endl << "Starting " <<metric_name<<" Registration" <<std::endl;</div><div>multiregistration->StartRegistration();</div><div>std::cout << "Optimizer stop condition = "</div><div><span class="Apple-tab-span" style="white-space:pre"> </span><< multiregistration->GetOptimizer()->GetStopConditionDescription()</div><div><span class="Apple-tab-span" style="white-space:pre"> </span><< std::endl;</div><div><br></div><div>std::string prefix = "SSD";</div><div><br></div><div>// the parameters type are the same for both optimizer.</div><div>typename RSGDOptimizerType::ParametersType finalParameters;</div><div>finalParameters = multiregistration->GetLastTransformParameters();</div><div>std::cout << "Last Transform Parameters" << std::endl;</div><div><span class="Apple-tab-span" style="white-space:pre"> </span>std::cout << finalParameters << std::endl;</div><div><br></div><div>unsigned int numberOfIterations = 0;</div><div>double bestValue = 0;</div><div><br></div><div>typedef RSGDOptimizerType * OptimizerPointer;</div><div><br></div><div><br></div><div>OptimizerPointer opt = dynamic_cast< OptimizerPointer >(optimizer.GetPointer());</div><div>numberOfIterations = opt->GetCurrentIteration();</div><div>bestValue = opt->GetValue();</div><div><br></div><div>// Print out results</div><div>//</div><div>std::cout << "Result: " << std::endl;</div><div>std::cout << " Iterations = " << numberOfIterations << std::endl;</div><div>std::cout << " Metric value = " << bestValue << std::endl;</div><div><br></div><div>transform->SetParameters( finalParameters );</div><div><br></div><div><br></div><div>typename ResampleFilterType::Pointer resample = ResampleFilterType::New();</div><div><br></div><div>resample->SetTransform( transform );</div><div>resample->SetInput( movingImage );</div><div>resample->SetSize( fixedImage->GetLargestPossibleRegion().GetSize() );</div><div>resample->SetOutputOrigin( fixedImage->GetOrigin() );</div><div>resample->SetOutputSpacing( fixedImage->GetSpacing() );</div><div>resample->SetOutputDirection( fixedImage->GetDirection() );</div><div>resample->SetDefaultPixelValue( 100 );</div><div>resample->Update();</div><div><br></div><div>typedef itk::CastImageFilter<FixedImageType,OutputImageType > CastFilterType;</div><div><br></div><div><br></div><div>typename OutputWriterType::Pointer writer = OutputWriterType::New();</div><div>typename CastFilterType::Pointer caster = CastFilterType::New();</div><div><br></div><div><br></div><div><br></div><div>if(salidaResample)</div><div>{</div><div>writer->SetFileName( (prefix+"_"+filenameResample).c_str() );</div><div><br></div><div>caster->SetInput( resample->GetOutput() );</div><div>writer->SetInput( caster->GetOutput() );</div><div><br></div><div><span class="Apple-tab-span" style="white-space:pre"> </span>try</div><div><span class="Apple-tab-span" style="white-space:pre"> </span>{</div><div><span class="Apple-tab-span" style="white-space:pre"> </span>writer->Update();</div><div><span class="Apple-tab-span" style="white-space:pre"> </span>}</div><div><span class="Apple-tab-span" style="white-space:pre"> </span>catch( itk::ExceptionObject & err )</div><div><span class="Apple-tab-span" style="white-space:pre"> </span>{</div><div><span class="Apple-tab-span" style="white-space:pre"> </span>std::cerr << "ExceptionObject caught !" << std::endl;</div><div><span class="Apple-tab-span" style="white-space:pre"> </span>std::cerr << err << std::endl;</div><div><span class="Apple-tab-span" style="white-space:pre"> </span>return EXIT_FAILURE;</div><div><span class="Apple-tab-span" style="white-space:pre"> </span>}</div><div>}</div></div><div><br></div><div><br></div><div><br></div><div><br></div><div><br></div><div><br></div><div>Please any help will be useful</div><div apple-content-edited="true"><div style="word-wrap: break-word; -webkit-nbsp-mode: space; -webkit-line-break: after-white-space; "><div style="word-wrap: break-word; -webkit-nbsp-mode: space; -webkit-line-break: after-white-space; "><div style="word-wrap: break-word; -webkit-nbsp-mode: space; -webkit-line-break: after-white-space; "><div style="word-wrap: break-word; -webkit-nbsp-mode: space; -webkit-line-break: after-white-space; "><div style="word-wrap: break-word; -webkit-nbsp-mode: space; -webkit-line-break: after-white-space; "><div style="word-wrap: break-word; -webkit-nbsp-mode: space; -webkit-line-break: after-white-space; "><div style="word-wrap: break-word; -webkit-nbsp-mode: space; -webkit-line-break: after-white-space; "><div style="word-wrap: break-word; -webkit-nbsp-mode: space; -webkit-line-break: after-white-space; "><div style="word-wrap: break-word; -webkit-nbsp-mode: space; -webkit-line-break: after-white-space; "><div style="word-wrap: break-word; -webkit-nbsp-mode: space; -webkit-line-break: after-white-space; "><div><div style="color: rgb(0, 0, 0); font-family: Helvetica; font-size: medium; font-weight: normal; font-style: normal; "><div style="color: rgb(0, 0, 0); font-family: Helvetica; font-size: medium; font-weight: normal; font-style: normal; ">Thanks,</div><div style="color: rgb(0, 0, 0); font-family: Helvetica; font-size: medium; font-weight: normal; font-style: normal; ">__________________________________<br><font class="Apple-style-span" color="#b7bfc5">|</font> Ariel Hernán Curiale Ph.D Student<br><font class="Apple-style-span" color="#b7bfc5">|</font> ETSI Telecomunicación<br><font class="Apple-style-span" color="#b7bfc5">|</font> Universidad de Valladolid<br><font class="Apple-style-span" color="#b7bfc5">|</font> Campus Miguel Delibes<br><font class="Apple-style-span" color="#b7bfc5">|</font> 47011 Valladolid, Spain<br><font class="Apple-style-span" color="#b7bfc5">|</font> Phone: 983-423000 ext. 5590</div><div><font class="Apple-style-span" color="#b7b7b7" style="font-family: Helvetica; font-size: medium; font-weight: normal; font-style: normal; ">|</font> Web: <a href="http://www.curiale.com.ar/"><font class="Apple-style-span" color="#084ebe">www.curiale.com.ar</font></a><br><font class="Apple-style-span" color="#b7b7b7" style="font-family: Helvetica; font-size: medium; font-weight: normal; font-style: normal; ">|</font>_________________________________</div></div></div></div></div></div></div></div></div></div></div></div></div>
</div>
<br></body></html>