[Insight-users] Problems with the convergence using Regular step descendent and Bsplines

[Ariel Hernán Curiale]

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/19a7a994/attachment-0001.htm>