[Insight-users] Problems with the convergence using Regular step descendent and Bsplines
Demian Wassermann
demian at bwh.harvard.edu
Sun Sep 30 09:36:51 EDT 2012
Hello Ariel,
I don’t know if I’m stating something very evident to you but keep in mind that the registration objective functions usually have lots of local minima of which many will not yield the result you want (like destroying the image). One solution to this is to register in different steps using transforms of increasing degrees of freedom to approach the initial scenario of a complex transform with a simpler one.
I would advice performing a multi-step registration, first registering with rigid and/or affine transforms and then using a multi-scale approach for the deformable registration.
HTH
Demian
--
Demian Wassermann, PhD
demian at bwh.harvard.edu
LMI / PNL / SPL Labs
Harvard Medical School
Brigham and Women's Hospital
1249 Boylston, Boston, MA, USA
http://lmi.bwh.harvard.edu/~wassermann/
On Sep 29, 2012, at 9:14 AM, Ariel Hernán Curiale <curiale at gmail.com> wrote:
> Hi,
> Where I wrote "the optimizer can achieve ..." I meant "the optimizer can't achieve..."
>
> Sorry
> __________________________________
> | Ariel Hernán Curiale Ph.D Student
> | ETSI Telecomunicación
> | Universidad de Valladolid
> | Campus Miguel Delibes
> | 47011 Valladolid, Spain
> | Phone: 983-423000 ext. 5590
> | Web: www.curiale.com.ar
> |_________________________________
>
> El 29/09/2012, a las 13:54, Ariel Hernán Curiale escribió:
>
>> Hi,
>>
>> 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.
>>
>>
>> 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).
>>
>> If anyone have an example how to use MeansSquares metric, Bspline interpolation and RegularStepGradient with the successfully result please send me the code.
>>
>>
>> I use the code with following parameters and two image of 200x200:
>>
>> --numberTrheads 1
>> --numberLevels 5
>> --space 16
>> --metricType SSD
>> --optimizerType RSGD
>> --interpolatorType BSI
>> --maxIt 500
>> --relaxFact 0.8
>> --steplength 0.001 2
>>
>>
>>
>> An extract of my code consern to the RESGOptimizer, SSDMetric and BsplineInterpolator is:
>>
>>
>> typedef itk::SingleValuedNonLinearOptimizer OptimizerType;//Abstract Optimizer
>> typedef itk::RegularStepGradientDescentOptimizer RSGDOptimizerType;
>>
>>
>> typedef itk::ImageToImageMetric<FixedImageType,MovingImageType > MetricType;//Abstract Metric
>> typedef itk::MeanSquaresImageToImageMetric<FixedImageType,MovingImageType > SSDMetricType;
>>
>> typedef itk::InterpolateImageFunction< MovingImageType,double> InterpolatorType;
>> typedef itk::BSplineInterpolateImageFunction<MovingImageType,double, double> BsplineInterpolatorType;
>>
>> typedef itk::MultiResolutionImageRegistrationMethod<FixedImageType,MovingImageType > MultiRegistrationType;
>>
>>
>> typedef itk::MultiResolutionPyramidImageFilter<FixedImageType,MovingImageType > FixedImagePyramidType;
>> typedef itk::MultiResolutionPyramidImageFilter<FixedImageType,MovingImageType > MovingImagePyramidType;
>>
>> typedef itk::ResampleImageFilter<FixedImageType,MovingImageType > ResampleFilterType;
>>
>>
>>
>> typename MetricType::Pointer metric = NULL;
>> typename InterpolatorType::Pointer interpolator = NULL;
>> typename OptimizerType::Pointer optimizer = NULL;
>>
>> optimizer = RSGDOptimizerType::New();
>> metric = SSDMetricType::New();
>> interpolator = BsplineInterpolatorType::New();
>>
>>
>> typedef BsplineInterpolatorType * InterpolatorPointer;
>> InterpolatorPointer int_aux = dynamic_cast< InterpolatorPointer >(interpolator.GetPointer());
>> int_aux->SetSplineOrder(3);
>>
>>
>>
>> typename TransformType::Pointer transform = TransformType::New();
>> typename MultiRegistrationType::Pointer multiregistration = MultiRegistrationType::New();
>>
>> multiregistration->SetNumberOfThreads(this->m_NumberOfThreads);
>> typename FixedImagePyramidType::Pointer fixedImagePyramid = FixedImagePyramidType::New();
>> typename MovingImagePyramidType::Pointer movingImagePyramid = MovingImagePyramidType::New();
>>
>> typename FixedImageType::RegionType fixedRegion = fixedImage->GetBufferedRegion();
>>
>> multiregistration->SetMetric( metrics[i] );
>> multiregistration->SetOptimizer( optimizer );
>> multiregistration->SetInterpolator( interpolator );
>> multiregistration->SetTransform( transform );
>>
>> // Setup the registration
>> multiregistration->SetFixedImagePyramid( fixedImagePyramid );
>> multiregistration->SetMovingImagePyramid( movingImagePyramid );
>> multiregistration->SetFixedImage( fixedImage );
>> multiregistration->SetMovingImage( movingImage);
>> multiregistration->SetFixedImageRegion( fixedRegion );
>>
>>
>>
>> // Here we define the parameters of the BSplineDeformableTransform grid. We
>> // arbitrarily decide to use a grid with $5 \times 5$ nodes within the image.
>> // The reader should note that the BSpline computation requires a
>> // finite support region ( 1 grid node at the lower borders and 2
>> // grid nodes at upper borders). Therefore in this example, we set
>> // the grid size to be $8 \times 8$ and place the grid origin such that
>> // grid node (1,1) coincides with the first pixel in the fixed image.
>>
>> unsigned int numberOfGridNodesInOneDimension = 1 + fixedImage->GetLargestPossibleRegion().GetSize()[0] / space;
>>
>>
>>
>> #if ITK_VERSION_MAJOR < 4
>>
>>
>> typename TransformType::RegionType bsplineRegion;
>> typename TransformType::RegionType::SizeType gridSizeOnImage;
>> typename TransformType::RegionType::SizeType gridBorderSize;
>> typename TransformType::RegionType::SizeType totalGridSize;
>>
>> gridSizeOnImage.Fill( numberOfGridNodesInOneDimension );
>> gridBorderSize.Fill( SplineOrder ); // Border for spline order = 3 ( 1 lower, 2 upper )
>> totalGridSize = gridSizeOnImage + gridBorderSize;
>>
>> bsplineRegion.SetSize( totalGridSize );
>>
>>
>> typename TransformType::SpacingType spacing = fixedImage->GetSpacing();
>>
>>
>> typename TransformType::OriginType origin = fixedImage->GetOrigin();
>>
>> typename FixedImageType::SizeType fixedImageSize = fixedRegion.GetSize();
>>
>> for(unsigned int r=0; r<Dimension; r++)
>> {
>> spacing[r] *= static_cast<double>(fixedImageSize[r] - 1) /
>> static_cast<double>(gridSizeOnImage[r] - 1);
>> }
>>
>> typename FixedImageType::DirectionType gridDirection = fixedImage->GetDirection();
>> typename TransformType::SpacingType gridOriginOffset = gridDirection * spacing;
>>
>> typename TransformType::OriginType gridOrigin = origin - gridOriginOffset;
>>
>> transform->SetGridSpacing( spacing );
>> transform->SetGridOrigin( gridOrigin );
>> transform->SetGridRegion( bsplineRegion );
>> transform->SetGridDirection( gridDirection );
>>
>> #else
>> typename TransformType::PhysicalDimensionsType fixedPhysicalDimensions;
>> typename TransformType::MeshSizeType meshSize;
>> typename TransformType::OriginType fixedOrigin;
>>
>> for( unsigned int k=0; k< Dimension; k++ )
>> {
>> fixedOrigin = fixedImage->GetOrigin()[k];
>> fixedPhysicalDimensions[k] = fixedImage->GetSpacing()[k] *
>> static_cast<double>(
>> fixedImage->GetLargestPossibleRegion().GetSize()[k] - 1 );
>> }
>> meshSize.Fill( numberOfGridNodesInOneDimension - SplineOrder );
>>
>> transform->SetTransformDomainOrigin( fixedOrigin );
>> transform->SetTransformDomainPhysicalDimensions( fixedPhysicalDimensions );
>> transform->SetTransformDomainMeshSize( meshSize );
>> transform->SetTransformDomainDirection( fixedImage->GetDirection() );
>>
>> #endif
>>
>> const unsigned int numberOfParameters = transform->GetNumberOfParameters();
>>
>> typename TransformType::ParametersType parameters( numberOfParameters );
>>
>> parameters.Fill( 0.0 );
>>
>> transform->SetParameters( parameters );
>>
>> // We now pass the parameters of the current transform as the initial
>> // parameters to be used when the registration process starts.
>>
>>
>> multiregistration->SetInitialTransformParameters( transform->GetParameters() );
>>
>>
>> multiregistration->SetNumberOfLevels( numberLevels );
>>
>>
>> // Set the parameters of the RegularStepGradientDescentOptimizer Optimizer.
>> typedef RSGDOptimizerType * OptimizerPointer;
>>
>> OptimizerPointer opt = dynamic_cast< OptimizerPointer >(optimizer.GetPointer());
>> opt->SetNumberOfIterations( maxIt );
>> opt->SetRelaxationFactor( relaxFact );
>> opt->SetMaximumStepLength( steplength[1] );
>> opt->SetMinimumStepLength( steplength[0] );
>>
>>
>>
>> std::cout << "Initial Parameters = " << std::endl;
>> std::cout << transform->GetParameters() << std::endl;
>>
>>
>> CommandIterationUpdate::Pointer observer = CommandIterationUpdate::New();
>> optimizer->AddObserver( itk::IterationEvent(), observer );
>>
>>
>> std::string metric_name = metric->GetNameOfClass();
>>
>> std::cout << std::endl << "Starting " <<metric_name<<" Registration" <<std::endl;
>> multiregistration->StartRegistration();
>> std::cout << "Optimizer stop condition = "
>> << multiregistration->GetOptimizer()->GetStopConditionDescription()
>> << std::endl;
>>
>> std::string prefix = "SSD";
>>
>> // the parameters type are the same for both optimizer.
>> typename RSGDOptimizerType::ParametersType finalParameters;
>> finalParameters = multiregistration->GetLastTransformParameters();
>> std::cout << "Last Transform Parameters" << std::endl;
>> std::cout << finalParameters << std::endl;
>>
>> unsigned int numberOfIterations = 0;
>> double bestValue = 0;
>>
>> typedef RSGDOptimizerType * OptimizerPointer;
>>
>>
>> OptimizerPointer opt = dynamic_cast< OptimizerPointer >(optimizer.GetPointer());
>> numberOfIterations = opt->GetCurrentIteration();
>> bestValue = opt->GetValue();
>>
>> // Print out results
>> //
>> std::cout << "Result: " << std::endl;
>> std::cout << " Iterations = " << numberOfIterations << std::endl;
>> std::cout << " Metric value = " << bestValue << std::endl;
>>
>> transform->SetParameters( finalParameters );
>>
>>
>> typename ResampleFilterType::Pointer resample = ResampleFilterType::New();
>>
>> resample->SetTransform( transform );
>> resample->SetInput( movingImage );
>> resample->SetSize( fixedImage->GetLargestPossibleRegion().GetSize() );
>> resample->SetOutputOrigin( fixedImage->GetOrigin() );
>> resample->SetOutputSpacing( fixedImage->GetSpacing() );
>> resample->SetOutputDirection( fixedImage->GetDirection() );
>> resample->SetDefaultPixelValue( 100 );
>> resample->Update();
>>
>> typedef itk::CastImageFilter<FixedImageType,OutputImageType > CastFilterType;
>>
>>
>> typename OutputWriterType::Pointer writer = OutputWriterType::New();
>> typename CastFilterType::Pointer caster = CastFilterType::New();
>>
>>
>>
>> if(salidaResample)
>> {
>> writer->SetFileName( (prefix+"_"+filenameResample).c_str() );
>>
>> caster->SetInput( resample->GetOutput() );
>> writer->SetInput( caster->GetOutput() );
>>
>> try
>> {
>> writer->Update();
>> }
>> catch( itk::ExceptionObject & err )
>> {
>> std::cerr << "ExceptionObject caught !" << std::endl;
>> std::cerr << err << std::endl;
>> return EXIT_FAILURE;
>> }
>> }
>>
>>
>>
>>
>>
>>
>> Please any help will be useful
>> Thanks,
>> __________________________________
>> | Ariel Hernán Curiale Ph.D Student
>> | ETSI Telecomunicación
>> | Universidad de Valladolid
>> | Campus Miguel Delibes
>> | 47011 Valladolid, Spain
>> | Phone: 983-423000 ext. 5590
>> | Web: www.curiale.com.ar
>> |_________________________________
>>
>
> _____________________________________
> Powered by www.kitware.com
>
> Visit other Kitware open-source projects at
> http://www.kitware.com/opensource/opensource.html
>
> Kitware offers ITK Training Courses, for more information visit:
> http://www.kitware.com/products/protraining.php
>
> Please keep messages on-topic and check the ITK FAQ at:
> http://www.itk.org/Wiki/ITK_FAQ
>
> Follow this link to subscribe/unsubscribe:
> http://www.itk.org/mailman/listinfo/insight-users
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://www.itk.org/pipermail/insight-users/attachments/20120930/f16fb9e7/attachment.htm>
More information about the Insight-users
mailing list