基于B 样条变换的配准 itk 实现

#include "itkImageRegistrationMethodv4.h"
#include "itkMattesMutualInformationImageToImageMetricv4.h"

#include "itkTimeProbesCollectorBase.h"
#include "itkMemoryProbesCollectorBase.h"

#include "itkBSplineTransform.h"
#include "itkLBFGSBOptimizerv4.h"


#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"

#include "itkResampleImageFilter.h"
#include "itkCastImageFilter.h"
#include "itkSquaredDifferenceImageFilter.h"

//  The following section of code implements a Command observer
//  used to monitor the evolution of the registration process.
//
#include "itkCommand.h"
#include "itkNiftiImageIO.h"


std::string Rigid = "Rigid.nii.gz";
std::string Affine = "Affine.nii.gz";
std::string NonRigid = "NonRigid.nii.gz";

class CommandIterationUpdate : public itk::Command
{
public:
	using Self = CommandIterationUpdate;
	using Superclass = itk::Command;
	using Pointer = itk::SmartPointer<Self>;
	itkNewMacro(Self);

protected:
	CommandIterationUpdate() = default;

public:
	using OptimizerType = itk::LBFGSBOptimizerv4;
	using OptimizerPointer = const OptimizerType *;

	void
		Execute(itk::Object * caller, const itk::EventObject & event) override
	{
		Execute((const itk::Object *)caller, event);
	}

	void
		Execute(const itk::Object * object, const itk::EventObject & event) override
	{
		auto optimizer = static_cast<OptimizerPointer>(object);
		if (!(itk::IterationEvent().CheckEvent(&event)))
		{
			return;
		}
		std::cout << optimizer->GetCurrentIteration() << "   ";
		std::cout << optimizer->GetCurrentMetricValue() << "   ";
		std::cout << optimizer->GetInfinityNormOfProjectedGradient() << std::endl;
	}
};

void  main()
{

	constexpr unsigned int ImageDimension = 3;
	using PixelType = float;

	using FixedImageType = itk::Image<PixelType, ImageDimension>;
	using MovingImageType = itk::Image<PixelType, ImageDimension>;
	using ImageIOType = itk::NiftiImageIO;

	ImageIOType::Pointer gNiftiImageIO = ImageIOType::New();
	ImageIOType::Pointer gNiftiImageIO1 = ImageIOType::New();
	ImageIOType::Pointer gNiftiImageIO2= ImageIOType::New();

	using FixedImageReaderType = itk::ImageFileReader<FixedImageType>;
	FixedImageReaderType::Pointer fixedImageReader =
		FixedImageReaderType::New();
	fixedImageReader->SetFileName(Rigid.data());
	fixedImageReader->SetImageIO(gNiftiImageIO);
	fixedImageReader->Update();
	FixedImageType::ConstPointer fixedImage = fixedImageReader->GetOutput();
	FixedImageType::RegionType   fixedRegion = fixedImage->GetBufferedRegion();

	using MovingImageReaderType = itk::ImageFileReader<MovingImageType>;
	MovingImageReaderType::Pointer movingImageReader =
		MovingImageReaderType::New();
	movingImageReader->SetFileName(Affine.data());
	movingImageReader->SetImageIO(gNiftiImageIO1);
	movingImageReader->Update();
	MovingImageType::ConstPointer movingImage = movingImageReader->GetOutput();

	// Software Guide : BeginCodeSnippet
	const unsigned int     SpaceDimension = ImageDimension;
	constexpr unsigned int SplineOrder = 3;
	using CoordinateRepType = double;

	using TransformType =
		itk::BSplineTransform<CoordinateRepType, SpaceDimension, SplineOrder>;


	using RegistrationType =
		itk::ImageRegistrationMethodv4<FixedImageType, MovingImageType>;
	RegistrationType::Pointer registration = RegistrationType::New();

	
	TransformType::Pointer transform = TransformType::New();
	

	unsigned int numberOfGridNodesInOneDimension = 7;



	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(numberOfGridNodesInOneDimension - SplineOrder);

	transform->SetTransformDomainOrigin(fixedOrigin);
	transform->SetTransformDomainPhysicalDimensions(fixedPhysicalDimensions);
	transform->SetTransformDomainMeshSize(meshSize);
	transform->SetTransformDomainDirection(fixedImage->GetDirection());

	registration->SetInitialTransform(transform);
	registration->InPlaceOn();

	using ParametersType = TransformType::ParametersType;

	const unsigned int numberOfParameters = transform->GetNumberOfParameters();

	ParametersType parameters(numberOfParameters);

	parameters.Fill(0.0);

	transform->SetParameters(parameters);
	//  Software Guide : EndCodeSnippet

	using MetricType =
		itk::MattesMutualInformationImageToImageMetricv4<FixedImageType,
		MovingImageType>;
	MetricType::Pointer metric = MetricType::New();
	metric->SetNumberOfHistogramBins(32);
	metric->SetUseMovingImageGradientFilter(false);
	metric->SetUseFixedImageGradientFilter(false);
	metric->SetUseSampledPointSet(false);

	using OptimizerType = itk::LBFGSBOptimizerv4;
	OptimizerType::Pointer optimizer = OptimizerType::New();

	// Software Guide : BeginCodeSnippet
	const unsigned int numParameters = transform->GetNumberOfParameters();
	OptimizerType::BoundSelectionType boundSelect(numParameters);
	OptimizerType::BoundValueType     upperBound(numParameters);
	OptimizerType::BoundValueType     lowerBound(numParameters);

	boundSelect.Fill(0);
	upperBound.Fill(0.0);
	lowerBound.Fill(0.0);

	optimizer->SetBoundSelection(boundSelect);
	optimizer->SetUpperBound(upperBound);
	optimizer->SetLowerBound(lowerBound);

	optimizer->SetCostFunctionConvergenceFactor(1.e7);
	optimizer->SetGradientConvergenceTolerance(1e-35);
	optimizer->SetNumberOfIterations(200);
	optimizer->SetMaximumNumberOfFunctionEvaluations(200);
	optimizer->SetMaximumNumberOfCorrections(7);
	// Software Guide : EndCodeSnippet

	// Create the Command observer and register it with the optimizer.
	//
	CommandIterationUpdate::Pointer observer = CommandIterationUpdate::New();
	optimizer->AddObserver(itk::IterationEvent(), observer);

	// One level registration is performed using the shrink factor 1 and
	// smoothing sigma 1
	//
	constexpr unsigned int numberOfLevels = 1;

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

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

	registration->SetFixedImage(fixedImage);
	registration->SetMovingImage(movingImage);
	registration->SetMetric(metric);
	registration->SetOptimizer(optimizer);
	registration->SetNumberOfLevels(numberOfLevels);
	registration->SetSmoothingSigmasPerLevel(smoothingSigmasPerLevel);
	registration->SetShrinkFactorsPerLevel(shrinkFactorsPerLevel);


	// Add time and memory probes
	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");

		const OptimizerType::ConstPointer outputOptimizer =
			dynamic_cast<const OptimizerType *>(registration->GetOptimizer());
		if (outputOptimizer.IsNotNull())
		{
			std::cout << "Optimizer stop condition = "
				<< outputOptimizer->GetStopConditionDescription()
				<< std::endl;
		}
		else
		{
			std::cerr << "Output optimizer is null." << std::endl;
			return ;
		}
	}
	catch (const itk::ExceptionObject & err)
	{
		std::cerr << "ExceptionObject caught !" << std::endl;

		return ;
	}

	// While the registration filter is run, it updates the output transform
	// parameters with the final registration parameters
	OptimizerType::ParametersType finalParameters = transform->GetParameters();

	// Report the time and memory taken by the registration
	chronometer.Report(std::cout);
	memorymeter.Report(std::cout);

	using ResampleFilterType =
		itk::ResampleImageFilter<MovingImageType, FixedImageType>;

	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());

	// This value is set to zero in order to make easier to perform
	// regression testing in this example. However, for didactic
	// exercise it will be better to set it to a medium gray value
	// such as 100 or 128.
	resample->SetDefaultPixelValue(0);

	using OutputPixelType = unsigned char;

	using OutputImageType = itk::Image<OutputPixelType, ImageDimension>;

	using CastFilterType =
		itk::CastImageFilter<FixedImageType, OutputImageType>;

	using WriterType = itk::ImageFileWriter<OutputImageType>;


	WriterType::Pointer     writer = WriterType::New();
	CastFilterType::Pointer caster = CastFilterType::New();


	writer->SetFileName(NonRigid.data());
	writer->SetImageIO(gNiftiImageIO2);

	caster->SetInput(resample->GetOutput());
	writer->SetInput(caster->GetOutput());


	try
	{
		writer->Update();
	}
	catch (const itk::ExceptionObject & err)
	{
		std::cerr << "ExceptionObject caught !" << std::endl;
		
		return ;
	}


	return ;
}