VoxelMorph-无监督医学图像配准-代码讲解

VoxelMorph官方Github地址：https://github.com/voxelmorph/voxelmorph，本文按照官方的Tutorial提供的路线进行讲解。

原文：Visit the VoxelMorph tutorial to learn about VoxelMorph and Learning-based Registration. Here’s an additional small tutorial on warping annotations together with images, and another on template (atlas) construction with VoxelMorph.

主要分为：

• Additional small tutorial 根据注释对图像进行变换
• VoxelMorph tutorial VoxelMorph教程
• Template (atlas) construction 模版搭建教程
  本文Github地址：https://github.com/MaybeRichard/VoxelMorph-explain

第一部分：Additional small tutorial 根据注释对图像进行变换

环境及背景介绍：

本部分的官方代码地址：https://colab.research.google.com/drive/1V0CutSIfmtgDJg1XIkEnGteJuw0u7qT-#scrollTo=h1KXYz-Nauwn
这一部分主要介绍的是如何使用vxm库里的方法对图像进行变换，代码中的方法是随机生成一个矩阵，然后根据该矩阵对图像进行仿射变换。
环境要求：tensorflow2.4,VoxelMorph
完整代码如下：

# 安装和导包
!pip install voxelmorph
import voxelmorph as vxm
import tensorflow as tf
import numpy as np
import matplotlib.pyplot as plt
 
# 对输入图像进行适当的预处理
pad_amt = 10
(x_train,_),_ = tf.keras.datasets.mnist.load_data()
# float64占用的内存是float32的两倍，是float16的4倍；比如对于CIFAR10数据集，如果采用float64来表示，需要60000323238/1024**3=1.4G，光把数据集调入内存就需要1.4G；如果采用float32，只需要0.7G，如果采用float16，只需要0.35G左右；占用内存的多少，会对系统运行效率有严重影响；（因此数据集文件都是采用uint8来存在数据，保持文件最小）
im = x_train[0,...].astype('float')/255
# np.pad(需要填充的array，((上，下),(左，右))，mode=constant...),这一步是为了增加边缘，累死padding，作用是防止后面的平移导致其超出范围
im = np.pad(im,((pad_amt,pad_amt),(pad_amt,pad_amt)))

# 手工创建变换矩阵
aff = np.eye(3) # 创建主对角矩阵
aff[:2,:2]+=np.random.randn(2,2)*0.1 # 在上半部分的2*2区域加入随机噪声
aff[:2, 2] = np.random.uniform(-10, 10, (2, )) # 均匀分布，(low,high,size) aff[:2, 2]的尺寸是(2,)
aff_inv = np.linalg.inv(aff)

margin=10
nb_annotations = 5
annotations = [np.random.uniform(margin,f-margin,nb_annotations) for f in im.shape] # 创建两个注释，(48,48)表示两个
annotations = np.stack(annotations,1)

# np.newaxis 的功能是增加新的维度，但是要注意 np.newaxis 放的位置不同，产生的矩阵形状也不同。放在第一个，给行上增加维度，放在最后一个，给列上增加维度
im_keras = im[np.newaxis,...,np.newaxis]
aff_keras = aff[np.newaxis,:2,:]
annotations_keras = annotations[np.newaxis,...]

# 进行仿射变换
im_warped = vxm.layers.SpatialTransformer()([im_keras, aff_keras])
im_warped = im_warped[0, ..., 0]

# get dense field of inverse affine
field_inv = vxm.utils.affine_to_dense_shift(aff_inv[:-1,:], im.shape, shift_center=True)[np.newaxis, ...]

# warp annotations
data = [tf.convert_to_tensor(f, dtype=tf.float32) for f in [annotations_keras, field_inv]]
annotations_warped = vxm.utils.point_spatial_transformer(data)[0, ...].numpy()

# 结果可视化
plt.figure()
# note that x and y need to be flipped due to xy indexing in matplotlib
plt.subplot(1, 2, 1)
plt.imshow(im, cmap='gray')
plt.plot(*[annotations[:, f] for f in [1, 0]], 'o')  
plt.axis('off')

plt.subplot(1, 2, 2)
plt.imshow(im_warped, cmap='gray')
plt.plot(*[annotations_warped[:, f] for f in [1, 0]], 'o')
plt.axis('off');

代码分析与讲解：

1. 库的导入：

!pip install voxelmorph
import voxelmorph as vxm
import tensorflow as tf
import numpy as np
import matplotlib.pyplot as plt

2. 图像的输入与预处理：

教程中使用的是mnist的数据集，数据集的预处理步骤包括:

• 对原始图像进行边缘填充
• 灰度归一化（具体作用参考：https://blog.csdn.net/qq_41383956/article/details/88593538）

# 加载mnist数据集
# 其标准输出应为： (x_train, y_train), (x_test, y_test)，但是只需要x_train数据展示，所以其他的丢掉
(x_train,_),_ = tf.keras.datasets.mnist.load_data()
# 灰度归一化，从0-255压缩到0-1
im = x_train[0,...].astype('float')/255
# 边缘填充
# 这一步的目的是，在后面对图像进行变换时，原本的Mnist数据集的28*28在变换后，
# 数字可能会移出图像区域，所以扩大原始数据的大小，也就是空白部分，方便展示变换的效果。
# pad_amt设置为10，及补充的区域为10个pixel
pad_amt = 10
# np.pad(需要填充的array，((上，下),(左，右))，mode=constant...),这一步是为了增加边缘，可以理解为padding
# 原始数据28*28，填补大小，上下左右各10，处理后数据48*48
im = np.pad(im,((pad_amt,pad_amt),(pad_amt,pad_amt)))

数据处理前后效果：

3. 手动创建变换矩阵：

# 手动生成仿射变换矩阵，方便后面affine操作

# 创建主对角矩阵
aff = np.eye(3) 
 # 在左上半部分的2*2区域加入随机噪声
aff[:2,:2]+=np.random.randn(2,2)*0.1
# 前两行的第三列的内容使用(-10,10)之间的均匀随机采样数字来替换
#  np.random.uniform(low,high,size)，使用(2,)的原因是aff[:2,2]数组就是一个两行一列的值
aff[:2, 2] = np.random.uniform(-10, 10, (2, ))
# 对上面计算后的矩阵求逆
aff_inv = np.linalg.inv(aff)

# 手动生成annotation变换矩阵，方便后面warp操作
margin=10
nb_annotations = 5
 # 创建一个列表，其中包含两个annotations，每个中包含nb_annotations个随机数字，范围在(margin,f-margin)之间
annotations = [np.random.uniform(margin,f-margin,nb_annotations) for f in im.shape]
# np.stack的简单用法在我的notion中有说明：
# https://sandy-property-d5e.notion.site/np-stack-48a69e31be084aa98cd15ce7d093c2ec
annotations = np.stack(annotations,1)

处理后的数据分别为：
aff_inv:
[ 1.01 − 3.137 − 7.602 5.857 8.561 3.138 0 0 1 ] \begin{bmatrix} 1.01& -3.137 & -7.602\\ 5.857& 8.561& 3.138\\ 0& 0 &1 \end{bmatrix}