0

我正在使用 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, 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 提供了良好的配准结果。

有任何想法吗 ？

谢谢

4

1 回答 1

0

应用转换的最佳和最安全的方法是使用sitk.TransformixImageFilter()，但我认为您有理由以不同的方式进行转换。有了这个...

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

  1. 将中心转换为原点
  2. T应用您拥有的矩阵
  3. 将结果翻译回来，像这样
T = np.array([[0.82, 0.144, -13.1], [-0.144, 0.82, -11.9], [0, 0, 1]] )

center = np.array([[1, 0, 110], [0, 1, 128], [0, 0, 1]] )
center_inverse = np.array([[1, 0, -110], [0, 1, -128], [0, 0, 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 )

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

  1. 轴相对于 elastix 期望翻转
  2. 变换是从固定坐标到移动坐标
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
于 2020-11-03T03:01:45.620 回答