Are there any instances about the registration of 3D medical image by landmark?

Beginner Questions
,
1

When I set LandmarkBasedTransformInitializer instead of CenteredTransformInitializer, it seems a wrong result. Any suggestions?

2

Hi,

Here is a LandmarkBasedTransformInitializer example.

3

Many thanks, but I think 3D landmark is quite difference from 2D landmark.

4

ITK handles both 2D and 3D.

5

The following is my code, is that correct? I can’t get the right output from this registration.

using TransformType = itk::VersorRigid3DTransform<double>;

using OptimizerType = itk::RegularStepGradientDescentOptimizerv4<double>;
using MetricType = itk::MattesMutualInformationImageToImageMetricv4<InternalImageType, InternalImageType>;

using RegistrationType = itk::
	ImageRegistrationMethodv4<InternalImageType, InternalImageType, TransformType>;

auto metric = MetricType::New();
auto optimizer = OptimizerType::New();
auto registration = RegistrationType::New();

registration->SetMetric(metric);
unsigned int  iBins = 24;
metric->SetNumberOfHistogramBins(iBins);
metric->SetUseFixedImageGradientFilter(false);
metric->SetUseMovingImageGradientFilter(false);

using LandmarkBasedTransformInitializerType =
	itk::LandmarkBasedTransformInitializer<TransformType, InternalImageType, InternalImageType>;

LandmarkBasedTransformInitializerType::Pointer landmarkBasedTransformInitializer =
	LandmarkBasedTransformInitializerType::New();
//  Create source and target landmarks.
using LandmarkContainerType = LandmarkBasedTransformInitializerType::LandmarkPointContainer;
using LandmarkPointType = LandmarkBasedTransformInitializerType::LandmarkPointType;

LandmarkContainerType fixedLandmarks;
LandmarkContainerType movingLandmarks;

LandmarkPointType fixedPoint;
LandmarkPointType movingPoint;

for (size_t i = 0; i < fixedpoints.size(); i++)
{
	fixedPoint[0] = fixedpoints[i].x;
	fixedPoint[1] = fixedpoints[i].y;
	fixedPoint[2] = fixedpoints[i].z;
	fixedLandmarks.push_back(fixedPoint);
}

for (size_t i = 0; i < movedpoints.size(); i++)
{
	movingPoint[0] = movedpoints[i].x;
	movingPoint[1] = movedpoints[i].y;
	movingPoint[2] = movedpoints[i].z;
	movingLandmarks.push_back(movingPoint);
}

landmarkBasedTransformInitializer->SetFixedLandmarks(fixedLandmarks);
landmarkBasedTransformInitializer->SetMovingLandmarks(movingLandmarks);

auto transform = TransformType::New();
transform->SetIdentity();
using VersorType = TransformType::VersorType;
using VectorType = VersorType::VectorType;
VersorType rotation;
VectorType axis;
axis[0] = 0.0;
axis[1] = 0.0;
axis[2] = 1.0;
constexpr double angle = 0;
rotation.Set(axis, angle);
transform->SetRotation(rotation);

landmarkBasedTransformInitializer->SetTransform(transform);
landmarkBasedTransformInitializer->SetReferenceImage(casterfix->GetOutput());
landmarkBasedTransformInitializer->InitializeTransform();

registration->SetOptimizer(optimizer);
registration->SetFixedImage(casterfix->GetOutput());
registration->SetMovingImage(castermove->GetOutput());
registration->SetInitialTransform(transform);

using OptimizerScalesType = OptimizerType::ScalesType;
OptimizerScalesType optimizerScales(transform->GetNumberOfParameters());
const double translationScale = 1.0 / 1000.0;
optimizerScales[0] = 1.0;
optimizerScales[1] = 1.0;
optimizerScales[2] = 1.0;
optimizerScales[3] = translationScale;
optimizerScales[4] = translationScale;
optimizerScales[5] = translationScale;
optimizer->SetScales(optimizerScales);
optimizer->SetNumberOfIterations(100);
optimizer->SetLearningRate(0.5);
optimizer->SetMinimumStepLength(0.001);
optimizer->SetReturnBestParametersAndValue(true);

optimizer->SetRelaxationFactor(0.8);

auto observer = CommandIterationUpdatev4::New();
optimizer->AddObserver(itk::IterationEvent(), observer);

constexpr unsigned int numberOfLevels = 3;

RegistrationType::ShrinkFactorsArrayType shrinkFactorsPerLevel;
shrinkFactorsPerLevel.SetSize(3);
shrinkFactorsPerLevel[0] = 4;
shrinkFactorsPerLevel[1] = 2;
shrinkFactorsPerLevel[2] = 1;

RegistrationType::SmoothingSigmasArrayType smoothingSigmasPerLevel;
smoothingSigmasPerLevel.SetSize(3);
smoothingSigmasPerLevel[0] = 2;
smoothingSigmasPerLevel[0] = 1;
smoothingSigmasPerLevel[0] = 1;

registration->SetNumberOfLevels(numberOfLevels);
registration->SetSmoothingSigmasPerLevel(smoothingSigmasPerLevel);
registration->SetShrinkFactorsPerLevel(shrinkFactorsPerLevel);

RegistrationType::MetricSamplingStrategyEnum samplingStrategy =
	RegistrationType::MetricSamplingStrategyEnum::RANDOM;
double samplingPercentage = 0.2;

registration->SetMetricSamplingStrategy(samplingStrategy);
registration->SetMetricSamplingPercentage(samplingPercentage);
registration->MetricSamplingReinitializeSeed(121213);

try
{
	registration->Update();
	std::cout << "Optimizer stop condition: "
		<< registration->GetOptimizer()->GetStopConditionDescription()
		<< std::endl;
}
catch (const itk::ExceptionObject & err)
{
	std::cerr << "ExceptionObject caught !" << std::endl;
	std::cerr << err << std::endl;
}
6

There are extra calls that are not needed, but I do not see anything standing out that would cause an issue. Was the landmark-based initialization evaluated outside the broader registration?

7

I just don’t know whether this is the right process of registration based on landmark or not.

8

Registering based on landmarks should be as simple as this:

The code in C++ should be only a little bit more verbose, but otherwise the same.