Ручная регистрация с SimpleElastix

Question

Я использую SimpleElastix (https://simpleelastix.github.io/) для регистрации (Affine) двух 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 дает хороший результат регистрации.

Есть идеи?

Спасибо

Answer 1

Самый лучший и безопасный способ применить преобразование - использовать 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. оси перевернуты относительно ожиданий эластикса
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