使用SimpleElastix手动注册

问题描述

我正在使用SimpleElastix（https://simpleelastix.github.io/）来注册（仿射）两个2D图像（参见附件）

。为此，我正在使用以下代码：

import SimpleITK as sitk 
elastixImageFilter = sitk.ElastixImageFilter()
elastixImageFilter.SetFixedImage(sitk.ReadImage("fixed_image.nii"))
elastixImageFilter.SetMovingImage(sitk.ReadImage("float_image.nii"))
elastixImageFilter.SetParameterMap(sitk.GetDefaultParameterMap("affine"))
resultimage=elastixImageFilter.Execute()
sitk.WriteImage(resultimage,"registred_affine.nii")

执行后者后，我获得包含转换矩阵的以下TransformParameters0.txt：

(Transform "AffineTransform")
(NumberOfParameters 6)
(TransformParameters 0.820320 0.144798 -0.144657 0.820386 -13.106613 -11.900934)
(InitialTransformParametersFileName "NoInitialTransform")
(UseBinaryFormatForTransformationParameters "false")
(HowToCombineTransforms "Compose")

// Image specific
(FixedImageDimension 2)
(MovingImageDimension 2)
(FixedInternalImagePixelType "float")
(MovingInternalImagePixelType "float")
(Size 221 257)
(Index 0 0)
(Spacing 1.0000000000 1.0000000000)
(Origin 0.0000000000 0.0000000000)
(Direction 1.0000000000 0.0000000000 0.0000000000 1.0000000000)
(UseDirectionCosines "true")

// AdvancedAffineTransform specific
(CenterOfRotationPoint 110.0000000000 128.0000000000)

// ResampleInterpolator specific
(ResampleInterpolator "FinalBSplineInterpolator")
(FinalBSplineInterpolationorder 3)

// Resampler specific
(Resampler "DefaultResampler")
(DefaultPixelValue 0.000000)
(ResultimageFormat "nii")
(ResultimagePixelType "float")
(CompressResultimage "false")

我的目标是使用这种矩阵变换来注册浮动图像并获得类似于SimpleElastix获得的注册图像。为此，我使用了这个小脚本：

import SimpleITK as sitk
import numpy as np

T= np.array([[0.82,0.144,-13.1],[-0.144,0.82,-11.9],[0,1]] ) #matrix transformation
 
img_moved_orig = plt.imread('moved.png')
img_fixed_orig = plt.imread('fixed.png')

img_transformed = np.zeros((img_moved_orig.shape[0],img_moved_orig.shape[1])) 
for i  in range(img_moved_orig.shape[0]): 
    for j in range(img_moved_orig.shape[1]): 
        pixel_data = img_moved_orig[i,j] 
        input_coords = np.array([i,j,1]) 
        i_out,j_out,_ = T @ input_coords 
        img_transformed[int(i_out),int(j_out)] = pixel_data

我获得了该注册图像，并将其与SimpleElastix的结果进行比较（请参见所附图像）

。我们可以观察到缩放比例尚未操作，转换存在问题。我想知道我是否错过了转换矩阵中的某些内容，因为SimpleElastix提供了良好的注册结果。

有什么想法吗？

谢谢

解决方法

应用转换的最佳和最安全的方法是使用sitk.TransformixImageFilter()，但是我想您有理由以其他方式进行转换。有了这种方式...

第一个问题：您必须考虑旋转中心。总矩阵执行以下操作：

将中心转换为原点
应用您拥有的矩阵T
像这样将结果翻译回

T = np.array([[0.82,0.144,-13.1],[-0.144,0.82,-11.9],[0,1]] )

center = np.array([[1,110],1,128],1]] )
center_inverse = np.array([[1,-110],-128],1]] )

total_matrix = center @ T @ center_inverse

我强烈建议您使用scikit-image进行转换。

from skimage.transform import AffineTransform
from skimage.transform import warp

total_affine = AffineTransform( matrix=total_matrix )
img_moving_transformed = warp( img_moved_orig,total_affine )

如果您真的必须自己进行转换，则需要在代码中进行两件事更改：

轴相对于elastix期望发生了翻转
从固定坐标到移动坐标的转换

img_transformed = np.zeros((img_moved_orig.shape[0],img_moved_orig.shape[1])) 
for i in range(img_moved_orig.shape[0]): 
    for j in range(img_moved_orig.shape[1]): 
        
        # j is the first dimension for the elastix transform
        j_xfm,i_xfm,_ = total_matrix @ np.array([j,i,1]) 

        pixel_data = 0
        # notice this annoying check we have to do that skimage handles for us
        if( j_xfm >= 0 and j_xfm < img_moved_orig.shape[1] and i_xfm >=0 and i_xfm < img_moved_orig.shape[0] ):
            # transformed coordinates index the moving image
            pixel_data = img_moved_orig[int(i_xfm),int(j_xfm),0] # "nearest-neighbor" interpolation

        # "loop" indices index the output space
        img_transformed[i,j] = pixel_data

elastix-itk image-processing medical-imaging python