Блог

Я использую 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 обеспечивает хороший результат регистрации.

Есть идеи?

Спасибо

Лучший и безопасный способ применить преобразование — использовать 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