m_UseSequentialSampling error while performing registration of two 3D images of same modality

I am trying to execute the code from following page for my registration
I removed all this exception handling code …as I am not much aware of it. and tried executing. the error comes while executing registration->update() and it takes to error at m_UseSequentialSampling. How to resolve this? please help me. I have uploaded screenshots also.


I tried with the data available in ITK

  • brainweb1e1a10f20.mha
  • brainweb1e1a10f20Rot10Tx15.mha
    I get this error when regestration is performing. How can I rectify this? what does this error mean?

(edited After sometime)
Now i got to know that, for my input images, if i define dimension as 2 instead of 3, My code runs without any error and i am getting registration output image for one slice.
Why is it happening? do i need to declare anything else when playing with 3D images??

You did not post the code you are using. And DeformableRegistration7 is already 3D.

And the problem was line 387, the resize call which allocates memory. VisualStudio usually points to the next line instead of the one which causes the crash.

please find the following code

#include "itkImageRegistrationMethod.h"
#include "itkMeanSquaresImageToImageMetric.h"
#include "itkTimeProbesCollectorBase.h"
#include "itkMemoryProbesCollectorBase.h"
#include "itkBSplineTransform.h"
#include "itkLBFGSBOptimizer.h"
#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"
#include "itkResampleImageFilter.h"
#include "itkCastImageFilter.h"
#include "itkSquaredDifferenceImageFilter.h"

int main()
{

    const    unsigned int    ImageDimension = 3;
typedef  float           PixelType;
typedef itk::Image< PixelType, ImageDimension >  FixedImageType;
typedef itk::Image< PixelType, ImageDimension >  MovingImageType;
const unsigned int SpaceDimension = ImageDimension;
const unsigned int SplineOrder = 3;
typedef double CoordinateRepType;
typedef itk::BSplineTransform<
	CoordinateRepType,
	SpaceDimension,
	SplineOrder >     TransformType;

typedef itk::LBFGSBOptimizer       OptimizerType;
typedef itk::MeanSquaresImageToImageMetric<
	FixedImageType,
	MovingImageType >    MetricType;
typedef itk::LinearInterpolateImageFunction<
	MovingImageType,
	double          >    InterpolatorType;
typedef itk::ImageRegistrationMethod<
	FixedImageType,
	MovingImageType >    RegistrationType;
MetricType::Pointer         metric = MetricType::New();
OptimizerType::Pointer      optimizer = OptimizerType::New();
InterpolatorType::Pointer   interpolator = InterpolatorType::New();
RegistrationType::Pointer   registration = RegistrationType::New();
registration->SetMetric(metric);
registration->SetOptimizer(optimizer);
registration->SetInterpolator(interpolator);

TransformType::Pointer  transform = TransformType::New();
registration->SetTransform(transform);

typedef itk::ImageFileReader< FixedImageType  > FixedImageReaderType;
typedef itk::ImageFileReader< MovingImageType > MovingImageReaderType;
FixedImageReaderType::Pointer  fixedImageReader = FixedImageReaderType::New();
MovingImageReaderType::Pointer movingImageReader = MovingImageReaderType::New();
fixedImageReader->SetFileName("C:\\VISUAL STUDIO\\projects\\bspline3dregistration1\\bin\\Debug\\brainweb1e1a10f20.mhd");
movingImageReader->SetFileName("C:\\VISUAL STUDIO\\projects\\bspline3dregistration1\\bin\\Debug\\brainweb1e1a10f20Rot10Tx15.mhd");
FixedImageType::ConstPointer fixedImage = fixedImageReader->GetOutput();
registration->SetFixedImage(fixedImage);
registration->SetMovingImage(movingImageReader->GetOutput());
fixedImageReader->Update();
FixedImageType::RegionType fixedRegion = fixedImage->GetBufferedRegion();
registration->SetFixedImageRegion(fixedRegion);

typedef TransformType::RegionType RegionType;
unsigned int numberOfGridNodes = 8;
TransformType::PhysicalDimensionsType   fixedPhysicalDimensions;
TransformType::MeshSizeType             meshSize;
TransformType::OriginType               fixedOrigin;
for (unsigned int i = 0; i < SpaceDimension; i++)
{
	fixedOrigin[i] = fixedImage->GetOrigin()[i];
	fixedPhysicalDimensions[i] = fixedImage->GetSpacing()[i] *
		static_cast<double>(
			fixedImage->GetLargestPossibleRegion().GetSize()[i] - 1);
}
meshSize.Fill(numberOfGridNodes - SplineOrder);
transform->SetTransformDomainOrigin(fixedOrigin);
transform->SetTransformDomainPhysicalDimensions(
	fixedPhysicalDimensions);
transform->SetTransformDomainMeshSize(meshSize);
transform->SetTransformDomainDirection(fixedImage->GetDirection());
typedef TransformType::ParametersType     ParametersType;
const unsigned int numberOfParameters =
	transform->GetNumberOfParameters();
ParametersType parameters(numberOfParameters);
parameters.Fill(0.0);
transform->SetParameters(parameters);

registration->SetInitialTransformParameters(transform->GetParameters());

std::cout << "Intial Parameters = " << std::endl;
std::cout << transform->GetParameters() << std::endl;

OptimizerType::BoundSelectionType boundSelect(transform->GetNumberOfParameters());
OptimizerType::BoundValueType upperBound(transform->GetNumberOfParameters());
OptimizerType::BoundValueType lowerBound(transform->GetNumberOfParameters());
boundSelect.Fill(0);
upperBound.Fill(0.0);
lowerBound.Fill(0.0);
optimizer->SetBoundSelection(boundSelect);
optimizer->SetUpperBound(upperBound);
optimizer->SetLowerBound(lowerBound);
optimizer->SetCostFunctionConvergenceFactor(1e+12);
optimizer->SetProjectedGradientTolerance(1.0);
optimizer->SetMaximumNumberOfIterations(500);
optimizer->SetMaximumNumberOfEvaluations(500);
optimizer->SetMaximumNumberOfCorrections(5);

itk::TimeProbesCollectorBase chronometer;
itk::MemoryProbesCollectorBase memorymeter;
std::cout << std::endl << "Starting Registration" << std::endl;
try
{
	memorymeter.Start("Registration");
	chronometer.Start("Registration");
	registration->Update();
	chronometer.Stop("Registration");
	memorymeter.Stop("Registration");
	std::cout << "Optimizer stop condition = "
		<< registration->GetOptimizer()->GetStopConditionDescription()
		<< std::endl;
}
catch (itk::ExceptionObject & err)
{
	std::cerr << "ExceptionObject caught !" << std::endl;
	std::cerr << err << std::endl;
	return EXIT_FAILURE;
}
OptimizerType::ParametersType finalParameters =
	registration->GetLastTransformParameters();
std::cout << "Last Transform Parameters" << std::endl;
std::cout << finalParameters << std::endl;
// Report the time taken by the registration
chronometer.Report(std::cout);
memorymeter.Report(std::cout);
// Software Guide : BeginCodeSnippet
transform->SetParameters(finalParameters);
// Software Guide : EndCodeSnippet
typedef itk::ResampleImageFilter<
	MovingImageType,
	FixedImageType >    ResampleFilterType;
ResampleFilterType::Pointer resample = ResampleFilterType::New();
resample->SetTransform(transform);
resample->SetInput(movingImageReader->GetOutput());
resample->SetSize(fixedImage->GetLargestPossibleRegion().GetSize());
resample->SetOutputOrigin(fixedImage->GetOrigin());
resample->SetOutputSpacing(fixedImage->GetSpacing());
resample->SetOutputDirection(fixedImage->GetDirection());
resample->SetDefaultPixelValue(100);
typedef  unsigned char  OutputPixelType;
typedef itk::Image< OutputPixelType, ImageDimension > OutputImageType;
typedef itk::CastImageFilter<
	FixedImageType,
	OutputImageType > CastFilterType;
typedef itk::ImageFileWriter< OutputImageType >  WriterType;
WriterType::Pointer      writer = WriterType::New();
CastFilterType::Pointer  caster = CastFilterType::New();
writer->SetFileName("C:\\VISUAL STUDIO\\projects\\bspline3dregistration1\\bin\\Debug\\outputimage1.mhd");
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;
}
 }

std::bad_alloc means that your computer does not have enough memory for what you are trying to do. Try closing other programs, increasing virtual memory size, and changing your pixel type from float to unsigned char (brainweb1e1a10f20.mha has UCHAR type, so you don’t need float for that input).

I tried as you said. But no use. My system has 16 GB RAM. I tried giving maximum allocation. But no use. Any other idea?

Try setting metric sampling percentage lower (e.g. 0.05, or 0.01). You can also try reducing numberOfGridNodes from 8 to 4. You should also use memory profiler to narrow down the problem.

1 Like