介紹如何在 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()
