介紹如何在 Python 中使用 ITK 函式庫的 v4 影像對準框架處理二維影像的對準，並提供 hello world 範例程式碼。
ITKv4 影像對準框架
ITKv4 影像對準框架跟上一版比較起來更為彈性，主要差異就是對準過程都在新的虛擬影像（virtual image）中進行，固定影像（fixed image）與調動影像（moving image）都透過轉換（transforms）與內插（interpolators）在虛擬影像空間以 metric 衡量配適程度，然後將結果交給 optimizer 更新轉換參數，重複這個過程直到結果收斂。
在這張 ITKv4 影像對準框架圖中，虛線的項目代表資料物件（data objects），實線的項目代表處理物件（process objects）。
設定影像類型
在 ITKv4 的架構下要進行影像對準時，首先要決定影像的類型，包含像素的類型以及影像的維度。
import itk
import matplotlib.pyplot as plt
# 影像像素資料類型
PixelType = itk.F
# 實體座標系統用的資料類型
ScalarType = itk.D
# 影像維度
Dimension = 2
# 影像類型
FixedImageType = itk.Image[PixelType, Dimension]
MovingImageType = itk.Image[PixelType, Dimension]
讀取影像
根據影像類型建立影像 reader，並且讀取固定影像與調動影像。
# 建立影像 Reader
FixedReaderType = itk.ImageFileReader[FixedImageType]
fixedImageReader = FixedReaderType.New()
fixedImageReader.SetFileName("BrainProtonDensitySliceBorder20.mhd")
fixedImageReader.Update()
# 讀取 Fixed 影像
fixedImage = fixedImageReader.GetOutput()
# 建立影像 Reader
MovingReaderType = itk.ImageFileReader[MovingImageType]
movingImageReader = MovingReaderType.New()
movingImageReader.SetFileName("BrainProtonDensitySliceShifted13x17y.mhd")
movingImageReader.Update()
# 讀取 Moving 影像
movingImage = movingImageReader.GetOutput()
這裡影像 reader 所指定的影像類型是指讀取影像之後所輸出的影像類型，如果這個類型跟影像檔案中實際的資料不符，影像 reader 會自動處理影像類型的轉換。
若要查看影像內部的詳細資訊，可以直接用 print 輸出影像物件：
# 查看影像詳細資訊
print(fixedImage)
Image (0x25a89c0)
RTTI typeinfo: itk::Image<float, 2u>
Reference Count: 2
Modified Time: 246
Debug: Off
Object Name:
Observers:
none
Source: (0x2960380)
Source output name: Primary
Release Data: Off
Data Released: False
Global Release Data: Off
PipelineMTime: 59
UpdateMTime: 247
RealTimeStamp: 0 seconds
LargestPossibleRegion:
Dimension: 2
Index: [0, 0]
Size: [221, 257]
BufferedRegion:
Dimension: 2
Index: [0, 0]
Size: [221, 257]
RequestedRegion:
Dimension: 2
Index: [0, 0]
Size: [221, 257]
Spacing: [1, 1]
Origin: [0, 0]
Direction:
1 0
0 1
IndexToPointMatrix:
1 0
0 1
PointToIndexMatrix:
1 0
0 1
Inverse Direction:
1 0
0 1
PixelContainer:
ImportImageContainer (0x3dd22d0)
RTTI typeinfo: itk::ImportImageContainer<unsigned long, float>
Reference Count: 1
Modified Time: 244
Debug: Off
Object Name:
Observers:
none
Pointer: 0x3fed990
Container manages memory: true
Size: 56797
Capacity: 56797使用 matplotlib 的 imshow 顯示兩張影像的內容。
# 分開顯示兩張影像
fig, axs = plt.subplots(1, 2)
axs[0].imshow(itk.GetArrayViewFromImage(fixedImage), cmap='gray')
axs[0].set_title('Fixed Image')
axs[1].imshow(itk.GetArrayViewFromImage(movingImage), cmap='gray')
axs[1].set_title('Moving Image')
plt.show()
影像對準
進行影像對準時，需要定義好 transform（影像轉換，例如線性或非線性轉換）、optimizer（尋找最佳解的方法）、metric（衡量兩張影像距離的方法）、interpolator（計算內差的方法）的類型，接著將各元件串接起來，進行影像對準。
# 影像轉換類型
TransformType = itk.TranslationTransform[ScalarType, Dimension]
# 建立初始轉換
initialTransform = TransformType.New()
# 建立 Optimizer 物件
OptimizerType = itk.RegularStepGradientDescentOptimizerv4[ScalarType]
optimizer = OptimizerType.New()
# 建立 Metric 物件
MetricType = itk.MeanSquaresImageToImageMetricv4[FixedImageType, MovingImageType]
metric = MetricType.New()
# 影像對準物件
RegistrationType = itk.ImageRegistrationMethodv4[FixedImageType, MovingImageType]
registration = RegistrationType.New()
# 設定 Metric
registration.SetMetric(metric)
# 設定 Optimizer
registration.SetOptimizer(optimizer)
# 設定初始轉換
registration.SetInitialTransform(initialTransform)
# 設定輸入影像
registration.SetFixedImage(fixedImage)
registration.SetMovingImage(movingImage)
# Interpolator 類型
FixedLinearInterpolatorType = itk.LinearInterpolateImageFunction[FixedImageType, ScalarType]
MovingLinearInterpolatorType = itk.LinearInterpolateImageFunction[FixedImageType, ScalarType]
# 建立 Interpolator 物件
fixedInterpolator = FixedLinearInterpolatorType.New()
movingInterpolator = MovingLinearInterpolatorType.New()
# 設定 Interpolator
metric.SetFixedInterpolator(fixedInterpolator)
metric.SetMovingInterpolator(movingInterpolator)
設定調動影像的初始轉換，先取得預設參數，檢查參數數量。
# Moving 影像所使用的初始轉換
movingInitialTransform = TransformType.New()
initialParameters = movingInitialTransform.GetParameters()
# 查看參數數量
print(initialParameters.GetNumberOfElements())
2
設定初始轉換的參數，此處兩個參數分別代表 X 與 Y 軸方向的位移。
# 設定初始轉換參數
initialParameters[0] = 0
initialParameters[1] = 0
movingInitialTransform.SetParameters(initialParameters)
registration.SetMovingInitialTransform(movingInitialTransform)
固定影像所使用的初始轉換通常都是 identity（維持不變），由於預設值就是 identity，所以亦可省略不寫：
# Fixed 影像所使用的初始轉換
identityTransform = TransformType.New()
identityTransform.SetIdentity()
registration.SetFixedInitialTransform(identityTransform)
optimizer 的 learning rate 是非常重要的參數，其指定 optimizer 每次的步長（step length），不適當的 learning rate 會讓影像對準無法收斂。
learning rate 如果設定太大，會造成 optimizer 不穩定，若設定太小則會讓收斂過程非常緩慢，簡易的設定方式是一開始先設定為較小的數值（例如 1.0 到 5.0 左右），等到其他參數都調整好、整體都可以穩定收斂之後，再慢慢升高 learning rate 直到開始出現收斂不穩定的狀況為止，最佳的 learning rate 就是讓疊代次數最少，但依然維持收斂的穩定性。
learning rate 在實際運算時會跟 metric 的梯度（gradient）相乘產生步長，所以步長也會跟 metric 的數值成正比，在 metric 數值尺度比較大的狀況下，要維持相同的 optimizer 行為，就需要較小的 learning rate。
# 設定 Learning Rate
optimizer.SetLearningRate(4)
每當導數（derivative）方向發生劇烈改變時，通常就是代表剛越過局部極值點（local extrema），此時 optimizer 就會降低 learning，而降低 learning rate 的比例就是由 relaxation factor 參數來控制，這個參數值介於 0 到 1 之間，預設值是 0.5。
# 設定 Relaxation Factor
optimizer.SetRelaxationFactor(0.5)
隨著影像對準的過程，learning rate 會漸漸減小，當其減小到非常小的值的時候，就代表已經收斂了，這個判斷收斂的門檻值可以透過 minimum step length 參數來設定，預設值為 1e-4。
# 設定 Minimum Step Length
optimizer.SetMinimumStepLength(0.001)
在疊代的過程也有可能會出現無法收斂的狀況，因此會設定疊代次數的上限值。
# 設定最高疊代次數
optimizer.SetNumberOfIterations(200)
ITKv4 支援多層級的影像對準（multi-level registration），在不同階段中採用不同解析度的虛擬空間（virtual space），並搭配不同的平滑度（smoothness）的固定影像與調動影像，而這些條件必須在影像對準計算開始之前就設定好，否則就只能採用預設值進行影像對準。
這裡我們採用最簡單的單一層級對準方式，亦不使用任何的降解析度（shrinking）或平滑化（smoothing）。
# 設定 Levels 相關參數
registration.SetNumberOfLevels(1)
registration.SetSmoothingSigmasPerLevel([0])
registration.SetShrinkFactorsPerLevel([1])
所有的組件與參數都設定好之後，就可以實際進行影像對準的計算了。
# 進行影像對準
registration.Update()
計算完成之後，取得影像對準結果並輸出幾個重要的參數值。
# 取得影像對準結果
transform = registration.GetTransform()
finalParameters = transform.GetParameters()
translationAlongX = finalParameters.GetElement(0)
translationAlongY = finalParameters.GetElement(1)
numberOfIterations = optimizer.GetCurrentIteration()
bestMetricValue = optimizer.GetValue()
# 輸出結果
print("Registration result: ")
print(" Translation X = " + str(translationAlongX))
print(" Translation Y = " + str(translationAlongY))
print(" Iterations = " + str(numberOfIterations))
print(" Metric value = " + str(bestMetricValue))
Registration result: Translation X = 13.001531179313957 Translation Y = 16.99948401485064 Iterations = 21 Metric value = 0.0012263178458750189
TranslationTransform 轉換的參數意義非常單純，兩個參數就是代表 X 軸與 Y 軸的位移。這裡的 Iterations 就是實際計算所使用的疊代次數，這個數字如果很大，可能代表 learning rate 設定太小，會造成計算時間拉長，或是結果尚未收斂。
轉換影像
在影像對準計算完成之後，只是得到了轉換函數的各個參數值，如果要得到調動影像轉換後的結果，必須把轉換函數套用至調動影像上進行轉換。
由於此處我們有指定調動影像的初始轉換，因此這裡必須以 CompositeTransform 將初始轉換與最終計算出來的轉換結合，再套用至調動影像上，才能正確將調動影像轉換至固定影像的空間中。
# 建立對準結果的轉換
CompositeTransformType = itk.CompositeTransform[ScalarType, Dimension]
outputCompositeTransform = CompositeTransformType.New()
outputCompositeTransform.AddTransform(movingInitialTransform)
outputCompositeTransform.AddTransform(registration.GetModifiableTransform())
轉換影像的工作必須使用 resampler 搭配轉換函數來處理：
# 建立輸出對準結果用的 Resampler
ResampleFilterType = itk.ResampleImageFilter
resampler = ResampleFilterType[MovingImageType, FixedImageType].New()
resampler.SetInput(movingImage)
resampler.SetTransform(outputCompositeTransform)
resampler.SetUseReferenceImage(True)
resampler.SetReferenceImage(fixedImage)
resampler.SetDefaultPixelValue(100)
顯示影像對準的結果：
# 分開顯示兩張影像
fig, axs = plt.subplots(1, 2)
axs[0].imshow(itk.GetArrayViewFromImage(fixedImage), cmap='gray')
axs[0].set_title('Fixed Image')
axs[1].imshow(itk.GetArrayViewFromImage(resampler.GetOutput()), cmap='gray')
axs[1].set_title('Result Image')
plt.show()
