使用 SimpleElastix 手动注册
Manual registration with SimpleElastix
我正在使用 SimpleElastix(https://simpleelastix.github.io/) for the registration (Affine) of the two 2D images (see attached) 。为此,我正在使用此代码:
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 提供了很好的配准结果。
有什么想法吗?
谢谢
应用转换的最好和最安全的方法是使用 sitk.TransformixImageFilter()
,但我假设您有理由采用不同的方法。有了这个...
第一个问题:你必须考虑旋转中心。
总矩阵执行以下操作:
- 将中心转换为原点
- 应用矩阵
T
你有
- 将结果翻译回来,像这样
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 )
如果您真的必须自己进行转换,则您的代码中有两处需要更改:
- 轴相对于 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
我正在使用 SimpleElastix(https://simpleelastix.github.io/) for the registration (Affine) of the two 2D images (see attached)
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 的结果进行比较(见附图)
有什么想法吗?
谢谢
应用转换的最好和最安全的方法是使用 sitk.TransformixImageFilter()
,但我假设您有理由采用不同的方法。有了这个...
第一个问题:你必须考虑旋转中心。 总矩阵执行以下操作:
- 将中心转换为原点
- 应用矩阵
T
你有 - 将结果翻译回来,像这样
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 )
如果您真的必须自己进行转换,则您的代码中有两处需要更改:
- 轴相对于 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