使用 Scipy 对具有非均匀分辨率的 3D 图像进行仿射变换的行为

Question

我希望在不同分辨率的图像上应用affine transformation，在同质坐标中定义，但是当一个斧头的分辨率与其他斧头的分辨率不同时，我遇到了一个问题。

通常，由于只有仿射的平移部分取决于分辨率，所以我通过分辨率对平移部分进行归一化，并使用 scipy.ndimage.affine_transform 在图像上应用相应的仿射。

如果图像的分辨率对于所有轴都相同，则它可以完美运行，您可以在下面看到相同变换（缩放+平移，或旋转+平移，请参见下面的代码）的图像应用于图像，它的下采样版本（意味着较低的分辨率）。图像（几乎）完美匹配，据我所知，体素值的差异主要是由插值错误引起的。

但您可以看到，下采样变换后的图像和变换后（下采样用于比较）图像之间的形状重叠

Scale affine transformation applied on the same image, at two different (uniform) resolutions

Rotation affine transformation applied on the same image, at two different (uniform) resolutions

不幸的是，如果一个图像轴的分辨率与另一个不同（参见下面的代码），它适用于具有空非对角项（如平移或缩放）的仿射变换，但变换的结果给出一个完全错误的结果。

Rotation affine transformation applied on the same image, at two different (non-uniform) resolutions

在这里您可以看到代码的最小工作示例：

---

import numpy as np
import nibabel as nib
from scipy.ndimage import zoom
from scipy.ndimage import affine_transform
import matplotlib.pyplot as plt

################################
#### LOAD ANY 3D IMAGE HERE ####
################################
#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% TO BE DEFINED BY USER
orig_img = np.zeros((100, 100, 100))
orig_img[25:75, 25:75, 25:75] = 1

ndim = orig_img.ndim

################################
##### DEFINE RESOLUTIONS #######

#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% TO BE DEFINED BY USER
# Comment/uncomment to choose the resolutions (in mm) of the images
# ORIG_RESOLUTION = [1., 1., 1.]
# TARGET_RESOLUTION = [2., 2., 2.]

ORIG_RESOLUTION = [1., 0.5, 1.]
TARGET_RESOLUTION = [2., 2., 2.]

#####################################
##### DEFINE AFFINE TRANSFORM #######

affine_scale_translation = np.array([[1.5, 0.0, 0.0, -25.],
                                     [0.0, 0.8, 0.0, 0. ],
                                     [0.0, 0.0, 1.0, 0.  ],
                                     [0.0, 0.0, 0.0, 1.0]])
a = np.sqrt(2)/2.
affine_rotation_translation = np.array([[a  , a  , 0.0, -25.],
                                     [-a , a  , 0.0, 50. ],
                                     [0.0, 0.0, 1.0, 0.0 ],
                                     [0.0, 0.0, 0.0, 1.0]])

# #%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% TO BE DEFINED BY USER
# Comment/uncomment to choose the transformation to be applied 

# affine_tf, name_affine = affine_scale_translation, "Tf scale"
affine_tf, name_affine = affine_rotation_translation, "Tf rotation"

######################################################
######## DOWNSAMPLE IMAGE TO LOWER RESOLUTION ########
######################################################

downsample_img = zoom(orig_img,
                      zoom=np.array(ORIG_RESOLUTION)/np.array(TARGET_RESOLUTION),
                      prefilter=False, order=1)

##############################################################################
######## APPLY AFFINE TRANSFORMATION TO ORIGINAL AND DOWNSAMPLE IMAGE ########
##############################################################################

affine_st_full_res, affine_st_low_res = affine_tf.copy(), affine_tf.copy()

# Inverse transform as affine_transform apply the tf from the target space to the original space
affine_st_full_res, affine_st_low_res = np.linalg.inv(affine_st_full_res), np.linalg.inv(affine_st_low_res)

# Normalize translation part (normally expressed in millimeters) for the resolution
affine_st_full_res[:ndim, ndim] = affine_st_full_res[:ndim, ndim] / ORIG_RESOLUTION
affine_st_low_res[:ndim, ndim] = affine_st_low_res[:ndim, ndim] / TARGET_RESOLUTION


# Apply transforms on images of different resolutions

orig_tf_img = affine_transform(orig_img, affine_st_full_res, prefilter=False, order=1)
downsample_tf_img = affine_transform(downsample_img, affine_st_low_res, prefilter=False, order=1)


# Downsample result at full resolution to be compared to result on downsample image

downsample_orig_tf_img = zoom(orig_tf_img, zoom=np.array(
                              ORIG_RESOLUTION)/np.array(TARGET_RESOLUTION),
                              prefilter=False, order=1)

# print(orig_img.shape)
# print(downsample_img.shape)
# print(orig_tf_img.shape)
# print(downsample_orig_tf_img.shape)

###############################
######## VISUALISATION ########
###############################
# We'll visualize in 2D the slice at the middle of the z (third) axis of the image, in both resolution
mid_z_slice_full, mid_z_slice_low = orig_img.shape[2]//2, downsample_img.shape[2]//2

fig, ((ax1, ax2, ax3), (ax4, ax5, ax6)) = plt.subplots(nrows=2, ncols=3)

ax1.imshow(orig_img[:, :, mid_z_slice_full], cmap='gray')
ax1.axis('off')
ax1.set_title('1/ Origin image, at full res: {}'.format(ORIG_RESOLUTION))

ax2.imshow(downsample_img[:, :, mid_z_slice_low], cmap='gray')
ax2.axis('off')
ax2.set_title('2/ Downsampled image, at low res: {}'.format(TARGET_RESOLUTION))

ax3.imshow(downsample_tf_img[:, :, mid_z_slice_low], cmap='gray')
ax3.axis('off')
ax3.set_title('3/ Transformed downsampled image')

ax4.imshow(orig_tf_img[:, :, mid_z_slice_full], cmap='gray')
ax4.axis('off')
ax4.set_title('4/ Transformed original image')

ax5.imshow(downsample_orig_tf_img[:, :, mid_z_slice_low], cmap='gray')
ax5.axis('off')
ax5.set_title('5/ Downsampled transformed image')


error = ax6.imshow(np.abs(downsample_tf_img[:, :, mid_z_slice_low] -\
                          downsample_orig_tf_img[:, :, mid_z_slice_low]), cmap='hot')
ax6.axis('off')
ax6.set_title('Error map between 3/ and 5/')

fig.colorbar(error)

plt.suptitle('Result for {} applied on {} and {} resolution'.format(name_affine, ORIG_RESOLUTION, TARGET_RESOLUTION))
plt.tight_layout()

plt.show()

---

Answer 1

下采样基本上改变了图像的基向量。这使得下采样就像翻译一样。

但旋转和平移不可交换（不可互换）。

https://gamedev.stackexchange.com/questions/16719/what-is-the-correct-order-to-multiply-scale-rotation-and-translation-matrices-f

因此，在进行下采样后，您必须将图像的内容转换回原始位置。这可以通过使用ORIG_RESOLUTION 组成仿射变换矩阵并在实际变换之前应用它来完成。 或者使用矩阵乘法（例如使用np.dot）来修改实际的转换。