Main-Repo/SourceCode/AIReconstruction/AIReconstructionCppCore/FusionHelper.cpp

#include "FusionHelper.h"


/// <summary>
/// 构造函数
/// </summary>
FusionHelper::FusionHelper()
{
	_fixedImage = VolumeImageType::New();
	_movingImage = VolumeImageType::New();
	_fixedImageImpporter = ImportFilterType::New();
	_movingImageImporter = ImportFilterType::New();
	_spacing = 0;
	_colorType = Gray8;
}

/// <summary>
/// 析构函数
/// </summary>
FusionHelper::~FusionHelper()
{

}

/// <summary>
/// 读入基础体数据
/// </summary>
/// <param name="baseVolumeInfo"></param>
bool FusionHelper::LoadBaseVolumeData(UniformVolumeDataInfo baseVolumeInfo)
{
	try
	{
		_spacing = baseVolumeInfo.spacing;
		_colorType = baseVolumeInfo.colorType;
		_fixedImageImpporter = CreateITKImageFromImportVolumeDatas(baseVolumeInfo,true);
		_fixedImage = _fixedImageImpporter->GetOutput();
		_fixedImage->Update();

#if ENABLE_OBSERVE
		WriteMhaImage(_fixedImage, "fixedImage");
#endif
		return true;
	}
	catch (const std::exception& ex)
	{
		ErrorMsg::SetErrorMsg(FusionError, { ex.what() });
		return false;
	}
}

/// <summary>
/// 将已有的体数据与另一个融合在一起
/// </summary>
/// <param name="anotherVolumeInfo"></param>
bool FusionHelper::FuseWithAnotherVolumeData( UniformVolumeDataInfo  anotherVolumeInfo)
{
	try
	{
		// 判断spacing是否一致
		if (anotherVolumeInfo.spacing != _spacing)
		{
			char strMsgBuff[32];
			std::snprintf(strMsgBuff, 32, "unexpected spacing value ( %f ) ", anotherVolumeInfo.spacing);
			ErrorMsg::SetErrorMsg(FusionError, { strMsgBuff });
			return false;
		}
		// 创建移动图像
		_movingImageImporter = CreateITKImageFromImportVolumeDatas(anotherVolumeInfo);
		_movingImage = _movingImageImporter->GetOutput();
		_movingImage->Update();

#if ENABLE_OBSERVE
		WriteMhaImage(_movingImage, "movingImage");
#endif
		//// 第一阶段配准，平移变换，归一化互相关测度，将两个体数据在Z轴对齐
		auto parameters3d = Get3DTransformationsWithCorrelationMetrics();
		const double translation3dX = parameters3d[0];
		const double translation3dY = parameters3d[1];
		const double translation3dZ = parameters3d[2];
		//const double translation3dX = -0.28574895842683723;
		//const double translation3dY = 0.27410058291910605;
		//const double translation3dZ = 0.055393289309434761;

		////第二阶段配准。根据第一阶段得到的Z轴偏移量，找到左右两侧体数据的Z轴中心位置附近的断面2D图像
		////2D平移变换 作为变换在X-O-Y平面内配准
		////采用梯度差值测度 Gradient Difference，该测度在旧的ITK配准框架里
		auto parameter2d = Get2DTransformationsWithGradientDifferenceMetric(translation3dZ);
		const double translation2dX = parameter2d[0];
		const double translation2dY = parameter2d[1];
		//const double translation2dX = 0.38425032458639652;
		//const double translation2dY = -0.0059800695074633745;

		// ToDo 观察拼接缝隙处

		// 对3D体数据进行重采样
		// 根据配准计算得到的变换矩阵，对移动图像进行重采样，此时两个体数据在世界坐标对齐
		// 注：配准过程是在世界坐标系下进行的，每进行一次变换，图像中的每个像素都按照矩阵M在世界坐标系中做一次重采样
		// 一般情况下，是直接将移动图像根据计算得到的变换矩阵进行变换，然后将其与固定图像拼在一起即可
		// 但考虑到一般图像的边缘处成像效果都不是很好（和皮肤接触的不够好，可能有黑色的区域），
		// 如果直接将移动图像进行平移后拼接，可能会在接缝处有黑色区域，影响视觉效果
		// 因此，这里将原本应该全部由移动图像进行的平移，分一部分给固定图像，将两个图像分别移动一部分，这样接缝处不会有黑色区域，整体视觉效果会比较好
		// 如果resampleFixedImage为false，则全部由移动图像完成平移\
		// 设置固定图像的平移初值（固定不动，全0）
		Transform3DType::OutputVectorType translation3dFixed;
		translation3dFixed[1] = 0.0;
		translation3dFixed[2] = 0.0;
		// 设置移动图像的平移初值
		Transform3DType::OutputVectorType translation3dMoving;
		// 考虑到扫查时，y方向都是贴紧皮肤的，可以认为Y方向偏移量为0
		translation3dMoving[1] = translation2dY;
		translation3dMoving[2] = translation3dZ;
		// 一般情况下，
		// 如果translation2dX>0，直接将移动图像X方向裁掉abs(translation2dX)，然后将其与固定图像拼在一起即可
		// 如果translation2dX<0, 直接将固定图像X方向裁掉abs(translation2dX)，然后将其与固定图像拼在一起即可
		// 但考虑到一般图像的边缘处成像效果都不是很好（和皮肤接触的不够好，可能有黑色的区域），
		// 如果直接将移动图像进行平移后拼接，可能会在接缝处有黑色区域，影响视觉效果
		// 因此，这里将原本应该全部由移动图像进行的平移，分一部分给固定图像，将两个图像分别移动一部分，这样接缝处不会有黑色区域，整体视觉效果会比较好
		// 如果要移动的距离太小（小于transThreshold），则不再将平移量分给两个图像分别进行，由一个图像完成即可
		bool resampleFixedImage = true;
		double transThreshold = 3 * _spacing;
		double absTranslation2dX = std::abs(translation2dX);
		if(translation2dX <0)
		{
			if(absTranslation2dX >= transThreshold)
			{
				translation3dFixed[0] = absTranslation2dX * 0.7;
				translation3dMoving[0] = absTranslation2dX * 0.3;
			}
			else
			{
				translation3dFixed[0] = absTranslation2dX;
				translation3dMoving[0] = 0;
			}
		}
		if(translation2dX ==0)
		{
			translation3dFixed[0] = 0;
			translation3dMoving[0] = 0;
			resampleFixedImage = false;
		}
		if(translation2dX > 0)
		{
			if(absTranslation2dX >= transThreshold)
			{
				translation3dFixed[0] = absTranslation2dX * 0.3;
				translation3dMoving[0] = absTranslation2dX * 0.7;
			}
			else
			{
				translation3dFixed[0] = 0;
				translation3dMoving[0] = absTranslation2dX;
				resampleFixedImage = false;
			}
		}
		// 分别对两幅图像进行平移变换
		// 平移固定图像
		VolumeImageType::Pointer translatedFixedImage;
		if(resampleFixedImage)
		{
			translatedFixedImage = Translate3DImage(_fixedImage, translation3dFixed);
		}
		else
		{
			translatedFixedImage = _fixedImage;
		}
		// 平移移动图像
		VolumeImageType::Pointer translatedMovingImage = Translate3DImage(_movingImage, translation3dMoving);

#if ENABLE_OBSERVE
		WriteMhaImage(translatedFixedImage, "translatedFixImage");
		WriteMhaImage(translatedMovingImage, "translatedMovingImage");
#endif
		// 将两幅图拼接在一起，拼好后的图像直接替换fixedImage（以便继续下一次拼接）
		// 由于固定图像在创建的时候是翻转了的，因此这里得到的拼接图像也应该是翻转的
		// 因此是movingImage在左，将movingImage翻转，fixedImage在右，fixedImage不翻转
		_fixedImageImpporter = ConcatTwoVolumeImageAlongX(translatedMovingImage, true,translatedFixedImage,false);
		_fixedImage = _fixedImageImpporter->GetOutput();
		_fixedImage->Update();

#if ENABLE_OBSERVE
		WriteMhaImage(_fixedImage, "fusedImage");
#endif
		return true;
	}
	catch (const std::exception& ex)
	{
		ErrorMsg::SetErrorMsg(FusionError, { ex.what() });
		return false;
	}
}

/// <summary>
/// 得到融合后的体数据
/// </summary>
/// <param name="fusedVolumeInfo"></param>
bool FusionHelper::GetFusedVolumeData(UniformVolumeDataInfo& fusedVolumeInfo)
{
	try
	{
		// 将固定图像（即融合后的图像）沿着x轴翻转回来
		using FlipFilterType = itk::FlipImageFilter<VolumeImageType>;
		FlipFilterType::FlipAxesArrayType flipArray;
		flipArray[0] = 1;
		flipArray[1] = 0;
		flipArray[2] = 0;
		FlipFilterType::Pointer flipFilter = FlipFilterType::New();
		flipFilter->SetFlipAxes(flipArray);
		flipFilter->SetInput(_fixedImage);
		VolumeImageType::Pointer flippedImageFixed = flipFilter->GetOutput();
		flippedImageFixed->Update();

		// 将得到的图像进行转换
		using CastFilterType = itk::CastImageFilter<VolumeImageType, OutputImageType>;
		CastFilterType::Pointer caster = CastFilterType::New();
		caster->SetInput(flippedImageFixed);
		OutputImageType::Pointer outputImage= caster->GetOutput();
		outputImage->Update();

		// 将融合后的图像复制到fusedVolumeInfo里
		VolumeImageType::RegionType imgRegion = flippedImageFixed->GetLargestPossibleRegion();
		VolumeImageType::SizeType imgSize = imgRegion.GetSize();
		int x = imgSize[0];
		int y = imgSize[1];
		int z = imgSize[2];
		OutputPixelType* imageData = outputImage->GetBufferPointer();
		int pixelNum = x * y * z;
		// ToDo 暂未将图像类型转换成不同的ColorType，只支持8位灰度图
		memcpy(fusedVolumeInfo.dataBuffer, imageData, pixelNum);
		fusedVolumeInfo.x = x;
		fusedVolumeInfo.y = y;
		fusedVolumeInfo.z = z;
		fusedVolumeInfo.spacing = _spacing;
		fusedVolumeInfo.colorType = _colorType;
		return true;
	}
	catch (const std::exception& ex)
	{
		ErrorMsg::SetErrorMsg(FusionError, { ex.what() });
		return false;
	}
}

/// <summary>
/// 从输入的VolumeData里创建ITKImage
/// </summary>
/// <param name="volumeInfo"></param>
/// <param name="image"></param>
/// <param name="flip"></param>
/// <returns></returns>
ImportFilterType::Pointer FusionHelper::CreateITKImageFromImportVolumeDatas(UniformVolumeDataInfo volumeInfo, bool flip)
{
	uint8_t* dataBuffer = volumeInfo.dataBuffer;
	int x = volumeInfo.x;
	int y = volumeInfo.y;
	int z = volumeInfo.z;
	// region的size
	SizeType size;
	size[0] = x;
	size[1] = y;
	size[2] = z;
	// region的origin
	IndexType start;
	start.Fill(0);
	// region
	RegionType region;
	region.SetIndex(start);
	region.SetSize(size);
	ImportFilterType::Pointer importFilter = ImportFilterType::New();
	importFilter->SetRegion(region);
	// origin
	OriginType origin[3] = { 0.0,0.0,0.0 };
	importFilter->SetOrigin(origin);
	// spacing
	SpacingType spacing[3];
	spacing[0] = volumeInfo.spacing;
	spacing[1] = volumeInfo.spacing;
	spacing[2] = volumeInfo.spacing;
	importFilter->SetSpacing(spacing);
	// 创建自己的Buffer（因为要把输入的volumeData里的数据转换一下数据类型）
	int numPixels = x * y * z;
	// 注意：itk采用SmartPointer（引用计数为0时，会自动调用Delete），
	// 如果itkOwnTheBuffer设为true，则由itk来管理所占用的空间，
	// 因为itkOwnTheBuffer为true，所以这里new出来的可以不用delete
	PixelType* localDataBuffer = new PixelType[numPixels];
	int bytesPerPixel = ImageHelper::GetBytesPerPixel(volumeInfo.colorType);
	// ToDo 这里没有根据ColorType将彩色图像转换成所需的灰度值
	if(flip)
	{
		for(int ni=0; ni<z*y; ni++)
		{
			for(int nx=0; nx<x;nx++)
			{
				int index = ni * x + x - nx;
				auto pixelVal = dataBuffer[index * bytesPerPixel];
				localDataBuffer[ni * x + nx] = static_cast<PixelType>(pixelVal);
			}
		}
	}
	else
	{
		for (int ni = 0; ni < numPixels; ni++)
		{
			auto pixelVal = dataBuffer[ni * bytesPerPixel];
			localDataBuffer[ni] = static_cast<PixelType>(pixelVal);
		}
	}
	importFilter->SetImportPointer(localDataBuffer, numPixels, ItkOwnTheBuffer);
	// 返回图像
	return importFilter;
}

/// <summary>
/// 计算fixedImage和movingImage之间的平移
/// 用基于互相关的测度，用梯度下降的优化器
/// </summary>
/// <returns></returns>
Transform3DType::ParametersType FusionHelper::Get3DTransformationsWithCorrelationMetrics()
{
	using Optimizer3DType = itk::RegularStepGradientDescentOptimizerv4<double>;
	using Metric3DType = itk::CorrelationImageToImageMetricv4<VolumeImageType, VolumeImageType>;
	using Registration3DType = itk::ImageRegistrationMethodv4<VolumeImageType, VolumeImageType, Transform3DType>;

	// 设置regist3d
	Optimizer3DType::Pointer optimizer3d = Optimizer3DType::New();
	Metric3DType::Pointer metric3d = Metric3DType::New();
	Registration3DType::Pointer regist3d = Registration3DType::New();
	regist3d->SetOptimizer(optimizer3d);
	regist3d->SetMetric(metric3d);

	// 设置初始值（xyz两个方向的初始值都设为0）
	// 读入两幅图的尺寸
	Transform3DType::Pointer initTrans3d = Transform3DType::New();
	Optimizer3DType::ParametersType initialPara3d(initTrans3d->GetNumberOfParameters());
	initialPara3d[0] = 0.0;
	initialPara3d[1] = 0.0;
	initialPara3d[2] = 0.0;
	initTrans3d->SetParameters(initialPara3d);
	regist3d->SetMovingInitialTransform(initTrans3d);

	// 设置固定图像和移动图像
	regist3d->SetFixedImage(_fixedImage);
	regist3d->SetMovingImage(_movingImage);
	regist3d->SetObjectName("TranslationRegistration");

	// 设置收缩和平滑参数
	constexpr unsigned int numberOfLevels = 1;
	Registration3DType::ShrinkFactorsArrayType shrinkFactorsPerLevel;
	shrinkFactorsPerLevel.SetSize(numberOfLevels);
	shrinkFactorsPerLevel[0] = 8;
	Registration3DType::SmoothingSigmasArrayType smoothingSigmasPerLevel1;
	smoothingSigmasPerLevel1.SetSize(numberOfLevels);
	smoothingSigmasPerLevel1[0] = 2.0;
	regist3d->SetNumberOfLevels(numberOfLevels);
	regist3d->SetShrinkFactorsPerLevel(shrinkFactorsPerLevel);
	regist3d->SetSmoothingSigmasPerLevel(smoothingSigmasPerLevel1);

	// 设置学习率和补偿
	optimizer3d->SetLearningRate(20 * _spacing);
	optimizer3d->SetMinimumStepLength(0.5 * _spacing);
	//每次导数的方向突然改变时，优化器假定已经通过了一个局部极值，
	//并通过一个松弛因子(Relaxation Factor)来减小步长。取值范围 0~1，默认0.5
	optimizer3d->SetRelaxationFactor(0.9);
	optimizer3d->SetNumberOfIterations(100);
	optimizer3d->SetReturnBestParametersAndValue(true);

#if ENABLE_OBSERVE
	CommandIterationUpdate1::Pointer observer1 = CommandIterationUpdate1::New();
	optimizer3d->AddObserver(itk::IterationEvent(), observer1);
	using TranslationCommandType = RegistrationInterfaceCommand1<Registration3DType>;
	TranslationCommandType::Pointer command1 = TranslationCommandType::New();
	regist3d->AddObserver(itk::MultiResolutionIterationEvent(), command1);
#endif

	regist3d->Update();

#if ENABLE_OBSERVE
	std::cout << std::endl;
	std::cout
		<< "regist3d stop condition: "
		<< regist3d->GetOptimizer()->GetStopConditionDescription()
		<< std::endl;
#endif

	Transform3DType::ConstPointer transTransform = regist3d->GetTransform();
	Transform3DType::ParametersType parameters3d = transTransform->GetParameters();

#if ENABLE_OBSERVE
	const unsigned int numberOfIterations1 = optimizer3d->GetCurrentIteration();
	const double bestValue1 = optimizer3d->GetValue();
	const double  translation3d_x = parameters3d[0];
	const double  translation3d_y = parameters3d[1];
	const double  translation3d_z = parameters3d[2];
	std::cout << "3D transformation Result = " << std::endl;
	std::cout << " Iterations    = " << numberOfIterations1 << std::endl;
	std::cout << " Metric value  = " << bestValue1 << std::endl;
	std::cout << " Translation X = " << translation3d_x << std::endl;
	std::cout << " Translation Y = " << translation3d_y << std::endl;
	std::cout << " Translation Z = " << translation3d_z << std::endl;
#endif

	return parameters3d;
}

/// <summary>
/// 计算fixedImage和movingImage之间的平移
/// 用基于梯度差值的测度，用梯度下降的优化器
/// </summary>
/// <returns></returns>
Transform2DType::ParametersType FusionHelper::Get2DTransformationsWithGradientDifferenceMetric(double translationZ)
{
	// 第一步：抽取出z方向中心位置处的一张图（取_fixedImage和_moxingImage重叠部分的中间值）
	int pixelTransZ = std::abs(std::ceil(translationZ / _spacing));
	// 先得到固定图像和移动图像的z
	VolumeImageType::RegionType fixedImgRegion = _fixedImage->GetLargestPossibleRegion();
	VolumeImageType::SizeType fixedImgSize = fixedImgRegion.GetSize();
	int fixedImgZ = fixedImgSize[2];
	VolumeImageType::RegionType movingImgRegion = _movingImage->GetLargestPossibleRegion();
	VolumeImageType::SizeType movingImgSize = movingImgRegion.GetSize();
	int movingImgZ = movingImgSize[2];
	// 如果固定图像和移动图像之间没有重合区域，则直接报错
	int minImgZ = std::min(movingImgZ, fixedImgZ);
	if(pixelTransZ > minImgZ)
	{
		char strMsgBuff[32];
		std::snprintf(strMsgBuff, 32, "translationZ ( %f ) out of range.", translationZ);
		throw std::invalid_argument(strMsgBuff);
	}
	int indexCenter = (pixelTransZ + minImgZ) /2;
	int indexInFix, indexInMoving;
	if(translationZ <0)
	{
		indexInFix = indexCenter;
		indexInMoving = indexCenter - pixelTransZ;
	}
	else
	{
		indexInFix = indexCenter + pixelTransZ;
		indexInMoving = indexCenter;
	}
	// 提取出固定图像和移动图像中所需的那一帧
	auto sliceFixed = ExtractOneFrame(_fixedImage,indexInFix);
	auto sliceMoving = ExtractOneFrame(_movingImage,indexInMoving);
	sliceFixed->Update();
	sliceMoving->Update();

#if ENABLE_OBSERVE
	WriteBmpImage(sliceFixed, "sliceInCenterFixed");
	WriteBmpImage(sliceMoving, "sliceInCenterMoving");
#endif

	//// 第二步：把两幅图分别flip (为了让我们感兴趣的重叠区域在图像的左侧)
	// 注：因为改变策略，直接让固定图像翻转了，因此这里不用再翻转
	//using FlipFilterType = itk::FlipImageFilter<SliceImageType>;
	//// 只翻转0轴，即x轴
	//FlipFilterType::FlipAxesArrayType flipArray;
	//flipArray[0] = 1;
	//flipArray[1] = 0;
	//// 翻转固定图像
	//FlipFilterType::Pointer filterFixed = FlipFilterType::New();
	//filterFixed->SetFlipAxes(flipArray);
	//filterFixed->SetInput(sliceFixed);
	//SliceImageType::Pointer flipImageFixed = filterFixed->GetOutput();
	//flipImageFixed->Update();
	//// 翻转移动图像
	//FlipFilterType::Pointer filterMoving = FlipFilterType::New();
	//filterMoving->SetFlipAxes(flipArray);
	//filterMoving->SetInput(sliceMoving);
	//SliceImageType::Pointer flipImageMoving = filterMoving->GetOutput();
	//flipImageMoving->Update();

	// 第三步：对图像做高斯平滑
	using GaussianFilterType = itk::DiscreteGaussianImageFilter<SliceImageType, SliceImageType>;
	double sigma = 0.25;
	double variance = sigma * sigma;
	// 最大的核尺寸，默认是32
	const unsigned int maxKernelWidth = 32;
	// 平滑固定图像
	GaussianFilterType::Pointer gaussianFilterFixed = GaussianFilterType::New();
	gaussianFilterFixed->SetInput(sliceFixed);
	gaussianFilterFixed->SetVariance(variance);
	gaussianFilterFixed->SetMaximumKernelWidth(maxKernelWidth);
	gaussianFilterFixed->Update();
	SliceImageType::Pointer gaussianImageFixed = gaussianFilterFixed->GetOutput();
	gaussianImageFixed->Update();
	// 平滑移动图像
	GaussianFilterType::Pointer gaussianFilterMoving = GaussianFilterType::New();
	gaussianFilterMoving->SetInput(sliceMoving);
	gaussianFilterMoving->SetVariance(variance);
	gaussianFilterMoving->SetMaximumKernelWidth(maxKernelWidth);
	gaussianFilterMoving->Update();
	SliceImageType::Pointer gaussianImageMoving = gaussianFilterMoving->GetOutput();
	gaussianImageMoving->Update();

#if ENABLE_OBSERVE
	unsigned int dim = gaussianFilterFixed->GetFilterDimensionality();
	itk::FixedArray<double, 2> sigmaArr = gaussianFilterFixed->GetSigmaArray();
	double sigmaval = gaussianFilterFixed->GetSigma();
	itk::FixedArray<double, 2> errors = gaussianFilterFixed->GetMaximumError();
	int kernelWidth = gaussianFilterFixed->GetMaximumKernelWidth();
	std::cout << " dim = " << dim << " sigmaArr = " << sigmaArr << " sigma1 = " << sigmaval
		<< " errors = " << errors << " kernelWidth = " << kernelWidth << std::endl;

	WriteBmpImage(gaussianImageFixed, "gaussianImageFixed");
	WriteBmpImage(gaussianImageMoving, "gaussianImageMoving");
#endif

	// 第四步 设置ROI
	using RoiFilterType = itk::RegionOfInterestImageFilter<SliceImageType, SliceImageType>;
	float offsetX = 0.2f;
	int offsetXInPixel = offsetX * movingImgSize[0];
	// 给固定图像设置ROI
	RoiFilterType::Pointer roiFilterFixed = RoiFilterType::New();
	SliceImageType::RegionType roiFixed;
	SliceImageType::SizeType roiSizeFixed;
	roiSizeFixed[0] = fixedImgSize[0];
	roiSizeFixed[1] = fixedImgSize[1];
	roiSizeFixed[0] -= offsetXInPixel;
	SliceImageType::IndexType roiStartFixed;
	roiStartFixed.Fill(0);
	roiFixed.SetIndex(roiStartFixed);
	roiFixed.SetSize(roiSizeFixed);
	roiFilterFixed->SetRegionOfInterest(roiFixed);
	roiFilterFixed->SetInput(gaussianImageFixed);
	SliceImageType::Pointer roiImageFixed = roiFilterFixed->GetOutput();
	roiImageFixed->Update();
	// 给移动图像设置ROI
	RoiFilterType::Pointer roiFilterMoving = RoiFilterType::New();
	SliceImageType::RegionType roiMoving;
	SliceImageType::SizeType roiSizeMoving;
	roiSizeMoving[0] = movingImgSize[0];
	roiSizeMoving[1] = movingImgSize[1];
	roiSizeMoving[0] -= offsetXInPixel;
	SliceImageType::IndexType roiStartMoving;
	roiStartMoving.Fill(0);
	roiMoving.SetIndex(roiStartMoving);
	roiMoving.SetSize(roiSizeMoving);
	roiFilterMoving->SetRegionOfInterest(roiMoving);
	roiFilterMoving->SetInput(gaussianImageMoving);
	SliceImageType::Pointer roiImageMoving = roiFilterMoving->GetOutput();
	roiImageMoving->Update();

	// 第五步，将固定图像和移动图像进行2D配准
	using Optimizer2DType = itk::RegularStepGradientDescentOptimizer;
	using Registration2DType = itk::ImageRegistrationMethod<SliceImageType, SliceImageType>;
	using Metric2DType = itk::GradientDifferenceImageToImageMetric<SliceImageType, SliceImageType>;
	using InterpolatorType = itk::LinearInterpolateImageFunction<SliceImageType, double>;
	// 配准相关参数初始化
	Transform2DType::Pointer initTrans2d = Transform2DType::New();
	Optimizer2DType::Pointer optimizer2d = Optimizer2DType::New();
	InterpolatorType::Pointer interpolator2d = InterpolatorType::New();
	Registration2DType::Pointer regist2d = Registration2DType::New();
	Metric2DType::Pointer metric2d = Metric2DType::New();
	// 设置配准参数
	regist2d->SetOptimizer(optimizer2d);
	regist2d->SetTransform(initTrans2d);
	regist2d->SetInterpolator(interpolator2d);
	metric2d->SetDerivativeDelta(0.5);
	regist2d->SetMetric(metric2d);
	// 设置移动和固定图像
	regist2d->SetFixedImage(roiImageFixed);
	regist2d->SetMovingImage(roiImageMoving);
	regist2d->SetFixedImageRegion(gaussianImageFixed->GetBufferedRegion());
	// 设置初始值
	Registration2DType::ParametersType initialPara2d(initTrans2d->GetNumberOfParameters());
	initialPara2d[0] = 0.0;
	initialPara2d[1] = 0.0;
	regist2d->SetInitialTransformParameters(initialPara2d);
	// 设置步长
	// 根据spacing来设置步长？（Danzel原来代码里写的1mm是0.1个像素，所以步长设为了1.5, 最小步长是0.1）
	optimizer2d->SetMaximumStepLength(15 * _spacing);
	optimizer2d->SetMinimumStepLength(0.5 * _spacing);
	optimizer2d->SetNumberOfIterations(100);
	optimizer2d->SetGradientMagnitudeTolerance(1e-40);
	optimizer2d->MaximizeOn();

#if ENABLE_OBSERVE
	CommandIterationUpdate2::Pointer observer = CommandIterationUpdate2::New();
	optimizer2d->AddObserver(itk::IterationEvent(), observer);
#endif

	regist2d->Update();

#if ENABLE_OBSERVE
	std::cout << "Optimizer stop condition: "
		<< regist2d->GetOptimizer()->GetStopConditionDescription()
		<< std::endl;
#endif

	const Transform2DType::ParametersType parameters2d = regist2d->GetLastTransformParameters();

#if ENABLE_OBSERVE
	const unsigned int numberOfIterations2 = optimizer2d->GetCurrentIteration();
	const double bestValue2 = optimizer2d->GetValue();
	double translation2d_x = parameters2d[0];
	double translation2d_y = parameters2d[1];
	// Print out results
	std::cout << "在 X-O-Y 平面的配准 Result = " << std::endl;
	std::cout << " Iterations      = " << numberOfIterations2 << std::endl;
	std::cout << " Metric value    = " << bestValue2 << std::endl;
	std::cout << " Translation X   = " << translation2d_x << std::endl;
	std::cout << " Translation Y   = " << translation2d_y << std::endl;
#endif

	return parameters2d;
}

/// <summary>
/// 取其中的一帧图
/// </summary>
/// <param name="image"></param>
/// <param name="indexZ"></param>
/// <returns></returns>
SliceImageType::Pointer FusionHelper::ExtractOneFrame(VolumeImageType::Pointer image,int indexZ)
{
	using FilterType = itk::ExtractImageFilter<VolumeImageType, SliceImageType>;
	FilterType::Pointer filter = FilterType::New();
	filter->SetDirectionCollapseToSubmatrix();
	filter->InPlaceOn();
	RegionType imgRegion = image->GetLargestPossibleRegion();
	SizeType size = imgRegion.GetSize();
	size[2] = 0;
	IndexType start;
	start[0] = 0;
	start[1] = 0;
	start[2] = indexZ;
	RegionType desiredRegion;
	desiredRegion.SetSize(size);
	desiredRegion.SetIndex(start);
	filter->SetExtractionRegion(desiredRegion);
	filter->SetInput(image);
	SliceImageType::Pointer sliceImage = filter->GetOutput();
	sliceImage->Update();
	return sliceImage;
}

/// <summary>
/// 将体数据写成mha格式的文件到指定路径
/// </summary>
/// <param name="image"></param>
/// <param name="fileName"></param>
void FusionHelper::WriteMhaImage(VolumeImageType::Pointer imagePointer, std::string fileName)
{
	// 创建目录
	if(access(_outputFolderName.c_str(),0) == -1)
	{
		mkdir(_outputFolderName.c_str());
	}
	std::string subFolder = _outputFolderName +"\\Mha";
	if (access(subFolder.c_str(), 0) == -1)
	{
		mkdir(subFolder.c_str());
	}
	using ImageIOType = itk::MetaImageIO;
	ImageIOType::Pointer metaIO = itk::MetaImageIO::New();
	using WriteImageType = itk::Image<unsigned char, 3>;
	using CastSliceFilterType = itk::CastImageFilter<VolumeImageType, WriteImageType>;
	using WriterType = itk::ImageFileWriter<WriteImageType>;
	WriterType::Pointer writer = WriterType::New();
	CastSliceFilterType::Pointer cast = CastSliceFilterType::New();
	cast->SetInput(imagePointer);
	writer->SetImageIO(metaIO);
	writer->SetInput(cast->GetOutput());
	writer->SetFileName(subFolder+"\\"+fileName+".mha");
	writer->Update();
}

/// <summary>
/// 将切面图像写成bmp格式的文件到指定路径
/// </summary>
/// <param name="image"></param>
/// <param name="fileName"></param>
void FusionHelper::WriteBmpImage(SliceImageType* image, std::string fileName)
{
	// 创建目录
	if (access(_outputFolderName.c_str(), 0) == -1)
	{
		mkdir(_outputFolderName.c_str());
	}
	std::string subFolder = _outputFolderName + "\\Bmp";
	if (access(subFolder.c_str(), 0) == -1)
	{
		mkdir(subFolder.c_str());
	}
	using WriteImageType = itk::Image<unsigned char, 2>;
	using CastSliceFilterType = itk::CastImageFilter<SliceImageType, WriteImageType>;
	using WriterType = itk::ImageFileWriter<WriteImageType>;
	WriterType::Pointer writer = WriterType::New();
	CastSliceFilterType::Pointer cast = CastSliceFilterType::New();
	cast->SetInput(image);
	writer->SetImageIO(itk::BMPImageIO::New());
	writer->SetInput(cast->GetOutput());
	writer->SetFileName(subFolder+"\\" + fileName + ".bmp");
	writer->Update();
}

/// <summary>
/// 移动一幅图像
/// </summary>
/// <param name="translation"></param>
/// <param name="image"></param>
/// <returns></returns>
VolumeImageType::Pointer FusionHelper::Translate3DImage(VolumeImageType::Pointer image,Transform3DType::OutputVectorType translation)
{
	// 平移固定图像
	using ResampleFilterType = itk::ResampleImageFilter<VolumeImageType, VolumeImageType>;
	OutputImageType* outputImage = nullptr;
	Transform3DType::Pointer transform = Transform3DType::New();
	transform->SetOffset(translation);
	ResampleFilterType::Pointer resampler = ResampleFilterType::New();
	resampler->SetTransform(transform);
	resampler->SetInput(image);
	VolumeImageType::RegionType imgRegion = image->GetLargestPossibleRegion();
	VolumeImageType::SizeType imgSize = imgRegion.GetSize();
	resampler->SetSize(imgSize);
	auto origin = image->GetOrigin();
	resampler->SetOutputOrigin(origin);
	auto spacing = image->GetSpacing();
	resampler->SetOutputSpacing(spacing);
	auto direction = image->GetDirection();
	resampler->SetOutputDirection(direction);
	resampler->SetDefaultPixelValue(0);
	auto resultImage = resampler->GetOutput();
	resultImage->Update();
	return resultImage;
}

/// <summary>
/// 沿着x轴方向将两个Itk三维图像直接拼到一起
/// </summary>
/// <param name="oneImage"></param>
/// <param name="anotherImage"></param>
/// <returns></returns>
ImportFilterType::Pointer FusionHelper::ConcatTwoVolumeImageAlongX(VolumeImageType::Pointer oneImage, bool flipOneImage,
	VolumeImageType::Pointer anotherImage, bool flipAnotherImage)
{
	// 获取两幅图像的数据指针
	PixelType* imageDataOne = oneImage->GetBufferPointer();
	PixelType* imageDataAnother = anotherImage->GetBufferPointer();
	// 获取两幅图像的尺寸信息
	VolumeImageType::RegionType imgRegionOne = oneImage->GetLargestPossibleRegion();
	VolumeImageType::SizeType imgSizeOne = imgRegionOne.GetSize();
	int sizeXOne = imgSizeOne[0];
	int sizeYOne = imgSizeOne[1];
	int sizeZOne = imgSizeOne[2];
	VolumeImageType::RegionType imgRegionAnother = anotherImage->GetLargestPossibleRegion();
	VolumeImageType::SizeType imgSizeAnother = imgRegionAnother.GetSize();
	int sizeXAnother = imgSizeAnother[0];
	int sizeYAnother = imgSizeAnother[1];
	int sizeZAnother = imgSizeAnother[2];
	// 确定输出图像的尺寸
	int sizeXConcat = sizeXOne + sizeXAnother;
	int sizeYConcat = std::max(sizeYOne, sizeYAnother);
	int sizeZConcat = std::max(sizeZOne, sizeZAnother);
	// 创建所需空间
	int numPixels = sizeXConcat * sizeYConcat * sizeZConcat;
	PixelType* imageDataConcat = new PixelType[numPixels];
	memset(imageDataConcat, 0, sizeof(PixelType) * numPixels);
	int indexDst = 0;
	int indexSrc = 0;
	// 拷贝第一幅图到指定位置(只能拷在图像范围内的像素值)
	for (int nz = 0; nz < sizeZOne; nz++)
	{
		for (int ny = 0; ny < sizeYOne; ny++)
		{
			for (int nx = 0; nx < sizeXOne; nx++)
			{
				if(flipOneImage)
				{
					indexSrc = nz * sizeYOne * sizeXOne + ny * sizeXOne + sizeXOne - nx;
				}
				else
				{
					indexSrc = nz * sizeYOne * sizeXOne + ny * sizeXOne + nx;
				}
				indexDst = nz * sizeYConcat * sizeXConcat + ny * sizeXConcat + nx;
				imageDataConcat[indexDst] = imageDataOne[indexSrc];
			}
		}
	}
	// 拷贝第二幅图到指定位置（只能拷在图像范围内的像素值）将其直接贴在第一幅图的右侧
	for (int nz = 0; nz < sizeZAnother; nz++)
	{
		for (int ny = 0; ny < sizeYAnother; ny++)
		{
			for (int nx = 0; nx < sizeXAnother; nx++)
			{
				if(flipAnotherImage)
				{
					indexSrc = nz * sizeYAnother * sizeXAnother + ny * sizeXAnother + sizeXAnother - nx;
				}
				else
				{
					indexSrc = nz * sizeYAnother * sizeXAnother + ny * sizeXAnother + nx;
				}
				indexDst = nz * sizeYConcat * sizeXConcat + ny * sizeXConcat + nx + sizeXOne;
				imageDataConcat[indexDst] = imageDataAnother[indexSrc];
			}
		}
	}
	// 从已知dataBuffer里创建图像
	ImportFilterType::Pointer importFilter = ImportFilterType::New();
	// region的size
	SizeType size;
	size[0] = sizeXConcat;
	size[1] = sizeYConcat;
	size[2] = sizeZConcat;
	// region的origin
	IndexType start;
	start.Fill(0);
	// region
	RegionType region;
	region.SetIndex(start);
	region.SetSize(size);
	importFilter->SetRegion(region);
	// origin
	OriginType origin[3] = { 0.0,0.0,0.0 };
	importFilter->SetOrigin(origin);
	// spacing
	SpacingType spacing[3];
	spacing[0] = _spacing;
	spacing[1] = _spacing;
	spacing[2] = _spacing;
	importFilter->SetSpacing(spacing);
	// 图像区指针
	importFilter->SetImportPointer(imageDataConcat, numPixels, ItkOwnTheBuffer);
	// 返回图像
	return importFilter;
}