Hi, I’m getting an Error when I try to increase the number of control points in the Image Registration using B-Splines. The error I got is:
terminate called after throwing an instance of ‘itk::ExceptionObject’
what(): /home/felippe/SimpleITK-build/ITK/Modules/Core/Common/src/itkMultiThreader.cxx:399:
itk::ERROR: MultiThreader(0x4f58f40): Exception occurred during SingleMethodExecute
std::bad_alloc
The code I’m using is:
// This one header will include all SimpleITK filters and external
// objects.
#include <SimpleITK.h>
#include <iostream>
#include <stdlib.h>
#include <iomanip>
namespace sitk = itk::simple;
// use sitk's output operator for std::vector etc..
using sitk::operator<<;
class IterationUpdate
: public sitk::Command
{
public:
IterationUpdate( const sitk::ImageRegistrationMethod &m, const sitk::BSplineTransform &tx)
: m_Method(m),
m_BSplineTransform(tx)
{}
// Override method from parent which is called at for the requested event
virtual void Execute( )
{
if (m_Method.GetOptimizerIteration() == 0)
{
// The BSpline is resized before the first optimizer
// iteration is completed per level. Print the transform object
// to show the adapted BSpline transform.
std::cout << m_BSplineTransform.ToString() << std::endl;
}
// stash the stream state
std::ios state(NULL);
state.copyfmt(std::cout);
std::cout << std::fixed << std::setfill(' ') << std::setprecision( 5 );
std::cout << std::setw(3) << m_Method.GetOptimizerIteration();
std::cout << " = " << std::setw(10) << m_Method.GetMetricValue() << std::endl;
std::cout.copyfmt(state);
}
private:
const sitk::ImageRegistrationMethod &m_Method;
const sitk::BSplineTransform &m_BSplineTransform;
};
class MultiResolutionIterationUpdate
: public sitk::Command
{
public:
MultiResolutionIterationUpdate( const sitk::ImageRegistrationMethod &m)
: m_Method(m)
{}
// Override method from parent which is called at for the requested event
virtual void Execute( )
{
// The sitkMultiResolutionIterationEvent occurs before the
// resolution of the transform. This event is used here to print
// the status of the optimizer from the previous registration level.
if (m_Method.GetCurrentLevel() > 0)
{
std::cout << "Optimizer stop condition: " << m_Method.GetOptimizerStopConditionDescription() << std::endl;
std::cout << " Iteration: " << m_Method.GetOptimizerIteration() << std::endl;
std::cout << " Metric value: " << m_Method.GetMetricValue() << std::endl;
}
std::cout << "--------- Resolution Changing ---------" << std::endl;
}
private:
const sitk::ImageRegistrationMethod &m_Method;
};
int main()
{
//******************************************************
//Reading the Fixed Image
sitk::Image fixed = sitk::ReadImage( "/home/felippe/Área de Trabalho/Felippe/Mestrado/REGISTRO/Deformation_Registering_MATLAB/Volumes/Teste1/MNI152_T1_0.5mm.nii", sitk::sitkFloat32 );
//fixed = sitk::Normalize( fixed );
//Printing the dimensions of the fixed image
std::vector<unsigned int> fixed_dims = fixed.GetSize();
std::cout << "Fixed Image Dimensions: ";
for (auto i = fixed_dims.begin(); i != fixed_dims.end(); ++i)
std::cout << *i << ' ';
std::cout << std::endl << "Fixed Image Pixel Type: " << fixed.GetPixelIDTypeAsString() << std::endl;
std::cout << std::endl;
//Reading the Moving Image
sitk::Image moving = sitk::ReadImage( "/home/felippe/Área de Trabalho/Felippe/Mestrado/C_plus_plus/Codigos/build-Registration_ITK_CMAKE-Desktop_Qt_5_12_3_GCC_64bit-Default/mri_vbm6_transformed_affine.mha", sitk::sitkFloat32 );
//moving = sitk::Normalize( moving );
//Printing the dimensions of the moving image
std::vector<unsigned int> moving_dims = moving.GetSize();
std::cout << "Moving Image Dimensions: ";
for (auto i = moving_dims.begin(); i != moving_dims.end(); ++i)
std::cout << *i << ' ';
std::cout << std::endl << "Moving Image Pixel Type: " << moving.GetPixelIDTypeAsString() << std::endl;
std::cout << std::endl;
std::cout << std::endl;
//****************************************************************************
std::cout << "Spacing in the Fixed Image: ";
std::vector<double> fixed_spacing = fixed.GetSpacing();
for (auto i = fixed_spacing.begin(); i != fixed_spacing.end(); ++i)
std::cout << *i << ' ';
std::cout << std::endl << std::endl;
std::vector<unsigned int> transformDomainMeshSize(fixed.GetDimension(), 4);
sitk::BSplineTransform tx = sitk::BSplineTransformInitializer(fixed, transformDomainMeshSize);
std::cout << "Transform Domain Mesh Size: ";
for (auto i = transformDomainMeshSize.begin(); i != transformDomainMeshSize.end(); ++i)
std::cout << *i << ' ';
std::cout << std::endl << std::endl;
std::cout << "Initial Number of Parameters:" << tx.GetNumberOfParameters() << std::endl;
sitk::ImageRegistrationMethod R;
unsigned int number_of_histogram_bins = 32;
std::vector<double> samplingPercentage(3);
samplingPercentage[0] = 0.01;
samplingPercentage[1] = 0.01;
samplingPercentage[2] = 0.01;
R.SetMetricAsMattesMutualInformation(number_of_histogram_bins);
R.SetMetricSamplingStrategy(R.RANDOM);
R.SetMetricSamplingPercentagePerLevel(samplingPercentage);
R.MetricUseFixedImageGradientFilterOn();
const double learningRate = 5.0;
const unsigned int numberOfIterations = 100u;
const double convergenceMinimumValue = 1e-4;
const unsigned int convergenceWindowSize = 5;
R.SetOptimizerAsGradientDescentLineSearch( learningRate,
numberOfIterations,
convergenceMinimumValue,
convergenceWindowSize);
//R.SetOptimizerAsLBFGSB(convergenceMinimumValue, numberOfIterations);
R.SetInterpolator(sitk::sitkLinear);
const unsigned int numberOfLevels = 3;
std::vector<unsigned int> scaleFactors(numberOfLevels);
scaleFactors[0] = 1;
scaleFactors[1] = 2;
scaleFactors[2] = 5;
const bool inPlace = true;
R.SetInitialTransformAsBSpline(tx,
inPlace,
scaleFactors);
std::vector<unsigned int> shrinkFactors( numberOfLevels );
shrinkFactors[0] = 4;
shrinkFactors[1] = 2;
shrinkFactors[2] = 1;
R.SetShrinkFactorsPerLevel( shrinkFactors );
std::vector<double> smoothingSigmas( numberOfLevels );
smoothingSigmas[0] = 4.0;
smoothingSigmas[1] = 2.0;
smoothingSigmas[2] = 1.0;
R.SetSmoothingSigmasPerLevel( smoothingSigmas );
IterationUpdate cmd1(R, tx);
R.AddCommand( sitk::sitkIterationEvent, cmd1);
MultiResolutionIterationUpdate cmd2(R);
R.AddCommand( sitk::sitkMultiResolutionIterationEvent, cmd2);
std::cout << "Initializing the Registration!!" << std::endl << std::endl;
sitk::Transform outTx = R.Execute( fixed, moving );
std::cout << "-------" << std::endl;
std::cout << outTx.ToString() << std::endl;
std::cout << "Optimizer stop condition: " << R.GetOptimizerStopConditionDescription() << std::endl;
std::cout << " Iteration: " << R.GetOptimizerIteration() << std::endl;
std::cout << " Metric value: " << R.GetMetricValue() << std::endl;
//Saving the transformation
sitk::WriteTransform(outTx, "simple_transform.tfm");
//Applying the transformation
sitk::Image transformed_image = sitk::Resample(moving, fixed, outTx, sitk::sitkLinear, 0, moving.GetPixelID());
sitk::WriteImage(transformed_image, "mri_vbm6_transformed_Affine_BS_444.mha");
return 0;
}
When I set transformDomainMeshSize to 2 I can get sucess in the registration, but when I set 4 I get this error which I think is related with memory leak. So how can I solve this problem?
Thanks.