ITK配准框架示例及代码解析（Python)

初学医学图像配准，把ITK里自带的配准示例看了几遍，由于没有注释，所以自己整理了一下。
一、配准方法
采用刚性变换的方法进行全局变换
使用基于均方差的准则来量化变换后图像之间的相似性
优化器选择使用梯度下降法来根据相似性测度中得到的经过量化的测度优化变换参数

二、配准步骤
1、定义待配准图像类型: 维数, 像素类型，并输入待配准的两幅图像，参考图 Fixed Image,浮动图 Moving Image；
2、定义图像转换类型，把参考图像空间映射到待配准图像空间，这里采用的是刚性变换；
3、定义优化器类型及设置优化参数，这里采用梯度下降法，设置初始振幅为4、最小步长为0.001、松弛因子为0.5、最大迭代次数为200；
4、相似性测度模块采用均方差的准则；
5、定义配准方法类型并实例化已经定义的相关类型；
6、初始化配准类型的参数，并将参数传递给配准过程；
7、相似度函数输入到优化模块中进行最优化计算得到最终变换参数，这个过程一般通过
迭代来实现，即重复2~4步直到取得最大值；
8、得到最优化后的各项参数并输出；
9、通过最终得到的参数重采样待配准图像，得到配准结果；
10、比较配准前后和配准过程后的差异。并且输出最后的配准结果以及配准前后的差异结果。
配准过程是一个优化问题，配准过程每进行一次迭代，得到一测度值，将该测度值与
我们所设定的值进行比较，如果达到预期的效果则停止迭代，得到最终配准结果。当然，
迭代可能无限制进行，所以我们还需要设置一迭代上限。

三、配准结果

fixed image

moving image

output image

differenceImageBefore

differenceImageAfter

四、代码

import itk

fixedImageFile = 'D:/peizhunDATA/fixed.mhd'
movingImageFile = 'D:/peizhunDATA/moving.mhd'
outputImageFile = 'D:/peizhunDATA/output.png'
differenceImageAfterFile = 'D:/peizhunDATA/differenceImageAfter.png'
differenceImageBeforeFile = 'D:/peizhunDATA/differenceImageBefore.png'


# 定义待配准图像类型: 维数, 像素类型：
PixelType = itk.ctype('float')
fixedImage = itk.imread(fixedImageFile, PixelType)
movingImage = itk.imread(movingImageFile, PixelType)
Dimension = fixedImage.GetImageDimension()
FixedImageType = itk.Image[PixelType, Dimension]
MovingImageType = itk.Image[PixelType, Dimension]


# 定义图像转换类型，把参考图像空间映射到待配准图像空间的转换：
TransformType = itk.TranslationTransform[itk.D, Dimension]  # 转换类型:刚性（移动）变换
initialTransform = TransformType.New()


# 定义优化器的类型及优化参数
optimizer = itk.RegularStepGradientDescentOptimizerv4.New(      # 梯度下降法的一种
    LearningRate=4,                # 初始振幅
    MinimumStepLength=0.001,       # 最小步长
    RelaxationFactor=0.5,          # 松弛因子
    NumberOfIterations=200)        # 最大迭代次数（以防止无限次迭代）


# 定义相似性测度类型
metric = itk.MeanSquaresImageToImageMetricv4[              # 均方根距离
    FixedImageType,
    MovingImageType].New()


# 定义配准方法类型并实例化已经定义的相关类型
registration = itk.ImageRegistrationMethodv4.New(FixedImage=fixedImage,
                                                 MovingImage=movingImage,
                                                 Metric=metric,
                                                 Optimizer=optimizer,
                                                 InitialTransform=initialTransform)

"""    
应该注意: 配准过程是一个优化的过程, 优化过程每迭代一次, 就与相似性测度进行一次比较，
当相似性测度值达到我们设定的预期值,或者达到设定的迭代上限时就停止迭代,得到最终结果. 
"""

# 初始化配准类型的参数，并将参数传递给配准过程
movingInitialTransform = TransformType.New()
initialParameters = movingInitialTransform.GetParameters()  # 得到浮动图像的初始变换参数
# 用于变换的参数队列由沿着每一维的方向的平移值构成。设置浮动图像初始变换参数值到零以便将变换初始成恒等变换。
initialParameters[0] = 0    # 沿x方向的偏移值设置为0
initialParameters[1] = 0    # 沿y方向的偏移值设置为0
movingInitialTransform.SetParameters(initialParameters)  # 将设置后的参数队列设为浮动图像新的初始变换参数
registration.SetMovingInitialTransform(movingInitialTransform)  # 设置完成后传递给配准过程
identityTransform = TransformType.New()
identityTransform.SetIdentity()
registration.SetFixedInitialTransform(identityTransform)  # 将参考图像的初始变换参数传递给配准过程
registration.SetNumberOfLevels(1)
registration.SetSmoothingSigmasPerLevel([0])
registration.SetShrinkFactorsPerLevel([1])


# 更新配准，使每次迭代生效
registration.Update()


# 得到最后的X,Y方向的变换参数
transform = registration.GetTransform()
finalParameters = transform.GetParameters()
translationAlongX = finalParameters.GetElement(0)
translationAlongY = finalParameters.GetElement(1)


# 一共迭代的次数
numberOfIterations = optimizer.GetCurrentIteration()


# 迭代完成后的相似性测度值
bestValue = optimizer.GetValue()


# 输出迭代完成后的各项参数
print("Result = ")
print(" Translation X = " + str(translationAlongX))
print(" Translation Y = " + str(translationAlongY))
print(" Iterations    = " + str(numberOfIterations))
print(" Metric value  = " + str(bestValue))


# 定义复合变换
CompositeTransformType = itk.CompositeTransform[itk.D, Dimension]
outputCompositeTransform = CompositeTransformType.New()
outputCompositeTransform.AddTransform(movingInitialTransform)
outputCompositeTransform.AddTransform(registration.GetModifiableTransform())


# 设置重采样过滤器的相应参数
resampler = itk.ResampleImageFilter.New(Input=movingImage,
                                        Transform=outputCompositeTransform,
                                        UseReferenceImage=True,
                                        ReferenceImage=fixedImage)
resampler.SetDefaultPixelValue(100)   # 设置默认灰度值,用来 "突出" 显示映射后在浮动图像之外的区域


# 定义输出图像和像素的类型
OutputPixelType = itk.ctype('unsigned char')
OutputImageType = itk.Image[OutputPixelType, Dimension]


# 重采样过滤器的输出
caster = itk.CastImageFilter[FixedImageType,
                             OutputImageType].New(Input=resampler)


# 输出配准后图像
writer = itk.ImageFileWriter.New(Input=caster, FileName=outputImageFile)
writer.SetFileName(outputImageFile)
writer.Update()


# 通过itk.SubtractImageFilter比较"经过变换的浮动图像"与"固定图像"之间的差异
difference = itk.SubtractImageFilter.New(Input1=fixedImage,
                                         Input2=resampler)


# 由于两幅图像之间的差异可能对应非常小的亮度值，使用itk.RescaleIntensityImageFilter
# 重新调节亮度值,使其更加明显
intensityRescaler = itk.RescaleIntensityImageFilter[FixedImageType,
                                                    OutputImageType].New(
    Input=difference,
    OutputMinimum=itk.NumericTraits[OutputPixelType].min(),
    OutputMaximum=itk.NumericTraits[OutputPixelType].max())


# 输出配准后的差异比较图
resampler.SetDefaultPixelValue(1)
writer.SetInput(intensityRescaler.GetOutput())
writer.SetFileName(differenceImageAfterFile)
writer.Update()


# 输出配准前的差异比较图
resampler.SetTransform(identityTransform)
writer.SetFileName(differenceImageBeforeFile)
writer.Update()

ITK配准示例网址