[Insight-users] Problems with the convergence using Regular step descendent and Bsplines
Ariel Hernán Curiale
curiale at gmail.com
Sat Sep 29 08:14:49 EDT 2012
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
> |_________________________________
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://www.itk.org/pipermail/insight-users/attachments/20120929/56151fde/attachment.htm>
More information about the Insight-users
mailing list