Python ITK 3D 剛性影像對準 Rigid Registration 置中初始化教學與範例

介紹如何使用 VersorRigid3DTransform 剛性影像轉換搭配 CenteredTransformInitializer 置中初始化對影像進行對準。

影像對準

首先定義各種資料與物件的型別。

import itk
import matplotlib.pyplot as plt

# 定義資料型別
FixedImageType   = itk.Image[itk.F, 3]
MovingImageType  = itk.Image[itk.F, 3]
TransformType    = itk.VersorRigid3DTransform[itk.D]
OptimizerType    = itk.RegularStepGradientDescentOptimizerv4[itk.D]
RegistrationType = itk.ImageRegistrationMethodv4[FixedImageType,
                                                 MovingImageType]
MetricType       = itk.MeanSquaresImageToImageMetricv4[FixedImageType,
                                                       MovingImageType]

建立影像 reader，並指定影像來源檔案：

# 建立影像 Reader
fixedImageReader  = itk.ImageFileReader[FixedImageType].New()
movingImageReader = itk.ImageFileReader[MovingImageType].New()

# 指定影像來源檔案
fixedImageReader.SetFileName("brainweb1e1a10f20.mha")
movingImageReader.SetFileName("brainweb1e1a10f20Rot10Tx15.mha")

這裡所使用的原始影像可以從 Kitware Data 的網站上取得。

建立 registration、optimizer、metric 物件：

# 建立 Registration
registration = RegistrationType.New()

# 指定 Registration 的影像來源
registration.SetFixedImage(fixedImageReader.GetOutput())
registration.SetMovingImage(movingImageReader.GetOutput())

# 設定 Optimizer
optimizer = OptimizerType.New()
registration.SetOptimizer(optimizer)

# 設定 Metric
metric = MetricType.New()
registration.SetMetric(metric)

建立並設定初始轉換：

# 建立並設定初始轉換
initialTransform = TransformType.New()

# 使用 CenteredTransformInitializer 初始化影像位置
TransformInitializerType = itk.CenteredTransformInitializer[
    TransformType, FixedImageType, MovingImageType];
transformInitializer = TransformInitializerType.New()
transformInitializer.SetTransform(initialTransform)
transformInitializer.SetFixedImage(fixedImageReader.GetOutput())
transformInitializer.SetMovingImage(movingImageReader.GetOutput())
transformInitializer.MomentsOn()
transformInitializer.InitializeTransform()

# 自行指定影像初始旋轉角度
axis = [0, 0, 1]
angle = 0.05  # 弧度（radian），0 ~ 2*pi
rotation = initialTransform.GetVersor()
# rotation.SetIdentity()
# rotation.SetRotationAroundX(angle)
rotation.Set(axis, angle)
initialTransform.SetRotation(rotation)

# 設定初始轉換
registration.SetInitialTransform(initialTransform)

這裡我們將所有的轉換都放在同一個 initialTransform 初始轉換中，不再另外個別設定 fixed 與 moving 影像的初始轉換。

設定各參數尺度，旋轉的單位是弧度（radian），平移的單位就要看影像的 voxel size 而定：

# 設定各參數尺度
numOfParam = initialTransform.GetNumberOfParameters()
optimizerScales = itk.OptimizerParameters.D(numOfParam)
translationScale = 1.0 / 1000.0
optimizerScales.SetElement(0, 1.0)  # 旋轉
optimizerScales.SetElement(1, 1.0)  # 旋轉
optimizerScales.SetElement(2, 1.0)  # 旋轉
optimizerScales.SetElement(3, translationScale) # X 軸平移
optimizerScales.SetElement(4, translationScale) # Y 軸平移
optimizerScales.SetElement(5, translationScale) # Z 軸平移
optimizer.SetScales(optimizerScales)

設定其他的參數：

# 設定疊代次數上限
optimizer.SetNumberOfIterations(200)

# 設定 Learning Rate
optimizer.SetLearningRate(0.2)

# 設定最小步長
optimizer.SetMinimumStepLength(0.001)

# 紀錄並傳回最佳解
optimizer.SetReturnBestParametersAndValue(True)

設定回呼函數：

from IPython.display import clear_output

# 初始事件回呼函數
def cb_registration_start():
    global cb_metric_values
    global cb_current_iteration_number
    cb_metric_values = []
    cb_current_iteration_number = -1

# 結束事件回呼函數
def cb_registration_end():
    global cb_metric_values
    global cb_current_iteration_number
    del cb_metric_values
    del cb_current_iteration_number
    plt.close()

# 疊代事件回呼函數
def cb_registration_iteration():
    global cb_metric_values
    global cb_current_iteration_number

    # 忽略重覆產生的疊代事件
    if optimizer.GetCurrentIteration() == cb_current_iteration_number:
        return

    cb_current_iteration_number = optimizer.GetCurrentIteration()
    cb_metric_values.append(optimizer.GetValue())

    # 清空輸出
    clear_output(wait=True)

    # 繪製 Metric 圖形
    plt.plot(cb_metric_values)
    plt.xlabel('Iteration Number', fontsize=12)
    plt.ylabel('Metric Value', fontsize=12)
    plt.show()

# 加上 Observer，設定事件與回呼函數的對應關係
optimizer.AddObserver(itk.StartEvent(), cb_registration_start)
optimizer.AddObserver(itk.EndEvent(), cb_registration_end)
optimizer.AddObserver(itk.IterationEvent(), cb_registration_iteration)

設定 levels 相關參數：

# 設定 Levels 相關參數
registration.SetNumberOfLevels(1)
registration.SetSmoothingSigmasPerLevel([0])
registration.SetShrinkFactorsPerLevel([1])

進行影像對準，並取得影像對準結果：

# 進行影像對準
registration.Update()

# 取得影像對準結果
transform = registration.GetTransform()
finalParameters = transform.GetParameters()
versorX = finalParameters[0]
versorY = finalParameters[1]
versorZ = finalParameters[2]
finalTranslationX = finalParameters[3]
finalTranslationY = finalParameters[4]
finalTranslationZ = finalParameters[5]
numberOfIterations = optimizer.GetCurrentIteration()
bestMetricValue = optimizer.GetValue()

print("Registration result: ")
print(" versor X      = " + str(versorX))
print(" versor Y      = " + str(versorY))
print(" versor Z      = " + str(versorZ))
print(" Translation X = " + str(finalTranslationX))
print(" Translation Y = " + str(finalTranslationY))
print(" Translation Z = " + str(finalTranslationZ))
print(" Iterations    = " + str(numberOfIterations))
print(" Metric value  = " + str(bestMetricValue))

Registration result: 
 versor X      = -0.0004894967032320457
 versor Y      = 4.237733030565726e-05
 versor Z      = -0.08724172378868676
 Translation X = 2.647231867780532
 Translation Y = -17.466151639587682
 Translation Z = -0.0026474097082672884
 Iterations    = 21
 Metric value  = 43.983994481735564

影像轉換

根據影像對準所得到的參數，建立一個轉換：

# 根據參數建立轉換
finalTransform = TransformType.New()
finalTransform.SetFixedParameters(registration.GetOutput().Get().GetFixedParameters())
finalTransform.SetParameters(finalParameters)
matrix = finalTransform.GetMatrix()
offset = finalTransform.GetOffset()
print("Matrix:\n", itk.GetArrayFromMatrix(matrix))
print("Offset:\n", offset)

Matrix:
 [[ 9.84777760e-01  1.73818110e-01  1.69840568e-04]
 [-1.73818193e-01  9.84777284e-01  9.67866412e-04]
 [ 9.77576960e-07 -9.82654697e-04  9.99999517e-01]]
Offset:
 itkVectorD3 ([-15.0303, -0.0642312, 0.104991])

建立輸出對準結果用的 resampler，搭配轉換產生影像對準的結果：

# 建立輸出對準結果用的 Resampler
ResampleFilterType = itk.ResampleImageFilter[MovingImageType, FixedImageType]
resampler = ResampleFilterType.New()
resampler.SetInput(movingImageReader.GetOutput())
resampler.SetTransform(finalTransform)
resampler.SetUseReferenceImage(True)
resampler.SetReferenceImage(fixedImageReader.GetOutput())
resampler.SetDefaultPixelValue(100)

# 顯示各切面影像
fixedImageArray = itk.GetArrayViewFromImage(fixedImageReader.GetOutput())
resultImageArray = itk.GetArrayViewFromImage(resampler.GetOutput())

fig, axs = plt.subplots(2, 3, figsize=(10, 7))
axs[0,0].imshow(fixedImageArray[fixedImageArray.shape[0]//2,:,:], cmap='gray')
axs[0,0].set_title('Fixed Image(X)')
axs[0,1].imshow(fixedImageArray[:,fixedImageArray.shape[1]//2,:], cmap='gray')
axs[0,1].set_title('Fixed Image(Y)')
axs[0,2].imshow(fixedImageArray[:,:,fixedImageArray.shape[2]//2], cmap='gray')
axs[0,2].set_title('Fixed Image(Z)')
axs[1,0].imshow(resultImageArray[resultImageArray.shape[0]//2,:,:], cmap='gray')
axs[1,0].set_title('Result Image(X)')
axs[1,1].imshow(resultImageArray[:,resultImageArray.shape[1]//2,:], cmap='gray')
axs[1,1].set_title('Result Image(Y)')
axs[1,2].imshow(resultImageArray[:,:,resultImageArray.shape[2]//2], cmap='gray')
axs[1,2].set_title('Result Image(Z)')
plt.show()

將 fixed 影像跟對準結果的影像合併為 RGB 影像：

# 轉換影像數值
OutputPixelType = itk.UC
OutputImageType = itk.Image[OutputPixelType, 3]

fixedIntensityRescaler = itk.RescaleIntensityImageFilter[FixedImageType, OutputImageType].New()
fixedIntensityRescaler.SetInput(fixedImageReader.GetOutput())
fixedIntensityRescaler.SetOutputMinimum(itk.NumericTraits[OutputPixelType].min())
fixedIntensityRescaler.SetOutputMaximum(itk.NumericTraits[OutputPixelType].max())

resultIntensityRescaler = itk.RescaleIntensityImageFilter[FixedImageType, OutputImageType].New()
resultIntensityRescaler.SetInput(resampler.GetOutput())
resultIntensityRescaler.SetOutputMinimum(itk.NumericTraits[OutputPixelType].min())
resultIntensityRescaler.SetOutputMaximum(itk.NumericTraits[OutputPixelType].max())

# 合併為 RGB 影像
ComponentType = itk.UC
RGBPixelType = itk.RGBPixel[ComponentType]
RGBImageType = itk.Image[RGBPixelType, 3]

composeFilter = itk.ComposeImageFilter[OutputImageType, RGBImageType].New()
composeFilter.SetInput(0, fixedIntensityRescaler.GetOutput())
composeFilter.SetInput(1, resultIntensityRescaler.GetOutput())
composeFilter.SetInput(2, fixedIntensityRescaler.GetOutput())

# 顯示各切面影像
composeImageArray = itk.GetArrayViewFromImage(composeFilter.GetOutput())

fig, axs = plt.subplots(1, 3, figsize=(10, 5))
axs[0].imshow(composeImageArray[composeImageArray.shape[0]//2,:,:], cmap='gray')
axs[0].set_title('Result Image(X)')
axs[1].imshow(composeImageArray[:,composeImageArray.shape[1]//2,:], cmap='gray')
axs[1].set_title('Result Image(Y)')
axs[2].imshow(composeImageArray[:,:,composeImageArray.shape[2]//2], cmap='gray')
axs[2].set_title('Result Image(Z)')
plt.show()

參考資料：ImageRegistration5.py、ImageRegistration5.cxx、ImageRegistration6.cxx、ImageRegistration7.cxx、ImageRegistration8.cxx