医学图像配准之形变场可视化（绘制形变场）

本文讲述医学图像配准中形变场的可视化，包括一种直接通过工具查看的方法和两种手工绘制的方法。

首先来介绍一下形变场，一个大小为[W,H]的二维图像对应的形变场的大小是[W,H,2]，其中第三个维度的大小为2，分别表示在x轴和y轴方向的位移。同理，一个大小为[D,W,H]的三维图像对应的形变场的大小是[D,W,H,3]，其中第三个维度的大小为3，分别表示在x轴、y轴和z轴方向的位移。下图是一个二维脑部图像配准后得到的形变场。

如一开始所说，下面介绍3中方式来将形变场可视化。

1. 使用ITK-Snap可视化形变场

这种方式是最便捷也是我认为效果最好的，但是只适用于三维的形变场。首先将形变场保存为".nii"图像，用ITK-Snap工具打开，首先在下图红框的位置悬停一段时间，出现下图中的图标后点击，选择“Multi-Component Display”下面的“Grid”，就可以看到下图所示的效果。

2. 毛毛大神的方法

具体思路请看毛毛大神的文章plt.contour 绘制图像形变场（Deformation Field）。其实具体思路我没有完全看懂，大约是先产生一个规则的网格，然后加上形变场后展示出来。但是如果网格是[W,H]大小的，但是形变场是[W,H,2]大小的，所以只能选择一个维度相加，感觉有点奇怪。代码如下：

import matplotlib.pyplot as plt
import SimpleITK as sitk
import numpy as npdef grid2contour(grid, title):'''grid--image_grid used to show deform fieldtype: numpy ndarray, shape： (h, w, 2), value range：(-1, 1)'''assert grid.ndim == 3x = np.arange(-1, 1, 2.0 / grid.shape[1])y = np.arange(-1, 1, 2.0 / grid.shape[0])X, Y = np.meshgrid(x, y)Z1 = grid[:, :, 0] + 2  # remove the dashed lineZ1 = Z1[::-1]  # vertical flipZ2 = grid[:, :, 1] + 2plt.figure()plt.contour(X, Y, Z1, 15, levels=50, colors='k') #改变levels的值，可以改变形变场的密集程度plt.contour(X, Y, Z2, 15, levels=50, colors='k')plt.xticks(()), plt.yticks(())  # remove x, y ticksplt.title(title)plt.show()def show_grid():img = sitk.ReadImage("./2D1.nii")img_arr = sitk.GetArrayFromImage(img)img_shape = img_arr.shape# 起点、终点、步长（可为小数）x = np.arange(-1, 1, 2 / img_shape[1])y = np.arange(-1, 1, 2 / img_shape[0])X, Y = np.meshgrid(x, y)regular_grid = np.stack((X, Y), axis=2)grid2contour(regular_grid, "regular_grid")rand_field = np.random.rand(*img_shape[:2], 2)  # 参数前加*是以元组形式导入rand_field_norm = rand_field.copy()rand_field_norm[:, :, 0] = rand_field_norm[:, :, 0] * 2 / img_shape[1]rand_field_norm[:, :, 1] = rand_field_norm[:, :, 1] * 2 / img_shape[0]sampling_grid = regular_grid + rand_field_normgrid2contour(sampling_grid, "sampling_grid")img_arr[..., 0] = img_arr[..., 0] * 2 / img_shape[1]img_arr[..., 1] = img_arr[..., 1] * 2 / img_shape[0]img_grid = regular_grid + img_arrgrid2contour(img_grid, "img_grid")if __name__ == "__main__":show_grid()print("end")

在代码的第20行，改变levels的值，可以改变形变场中网格的密集程度。

3. 先生成规则网格，再用空间变换网络进行变形

这种方式是配准群里的一位朋友提出的，自己进行实现。思路是先生成规则网格并保存为图片，再用空间变换网络进行变形。代码如下：

import numpy as np
import SimpleITK as sitk
import matplotlib.pyplot as pltimport torch
import torch.nn as nn
import torch.nn.functional as F# 空间转换网络
class SpatialTransformer(nn.Module):# 1.生成网格grid；2.new_grid=grid+flow，即旧网格加上一个位移；3.将网格规范化到[-1,1]；4.根据新网格对原图进行采样def __init__(self, size, mode='bilinear'):"""Instiatiate the block:param size: size of input to the spatial transformer block:param mode: method of interpolation for grid_sampler"""super(SpatialTransformer, self).__init__()# Create sampling gridvectors = [torch.arange(0, s) for s in size]grids = torch.meshgrid(vectors)grid = torch.stack(grids)  # y, x, zgrid = torch.unsqueeze(grid, 0)  # add batchgrid = grid.type(torch.FloatTensor)self.register_buffer('grid', grid)self.mode = modedef forward(self, src, flow):"""Push the src and flow through the spatial transform block:param src: the original moving image:param flow: the output from the U-Net"""new_locs = self.grid + flowshape = flow.shape[2:]# Need to normalize grid values to [-1, 1] for resamplerfor i in range(len(shape)):new_locs[:, i, ...] = 2 * (new_locs[:, i, ...] / (shape[i] - 1) - 0.5)if len(shape) == 2:new_locs = new_locs.permute(0, 2, 3, 1)  # 维度置换，变为0,2,3,1new_locs = new_locs[..., [1, 0]]elif len(shape) == 3:new_locs = new_locs.permute(0, 2, 3, 4, 1)new_locs = new_locs[..., [2, 1, 0]]return F.grid_sample(src, new_locs, mode=self.mode)# 生成网格图片
def create_grid(size, path):num1, num2 = (size[0] + 10) // 10, (size[1] + 10) // 10  # 改变除数（10），即可改变网格的密度x, y = np.meshgrid(np.linspace(-2, 2, num1), np.linspace(-2, 2, num2))plt.figure(figsize=((size[0] + 10) / 100.0, (size[1] + 10) / 100.0))  # 指定图像大小plt.plot(x, y, color="black")plt.plot(x.transpose(), y.transpose(), color="black")plt.axis('off')  # 不显示坐标轴# 去除白色边框plt.gca().xaxis.set_major_locator(plt.NullLocator())plt.gca().yaxis.set_major_locator(plt.NullLocator())plt.subplots_adjust(top=1, bottom=0, left=0, right=1, hspace=0, wspace=0)plt.margins(0, 0)plt.savefig(path)  # 保存图像# plt.show()if __name__ == "__main__":out_path = r"C:\Users\zuzhiang\Desktop\new_img.jpg"  # 图片保存路径# 读入形变场phi = sitk.ReadImage("./2D1.nii")  # [324,303,2]phi_arr = torch.from_numpy(sitk.GetArrayFromImage(phi)).float()phi_shape = phi_arr.shape# 产生网格图片create_grid(phi_shape, out_path)img = sitk.GetArrayFromImage(sitk.ReadImage(out_path))[..., 0]img = np.squeeze(img)[np.newaxis, np.newaxis, :phi_shape[0], :phi_shape[1]]# 用STN根据形变场对网格图片进行变形STN = SpatialTransformer(phi_shape[:2])phi_arr = phi_arr.permute(2, 0, 1)[np.newaxis, ...]warp = STN(torch.from_numpy(img).float(), phi_arr)# 保存图片warp_img = sitk.GetImageFromArray(warp[0, 0, ...].numpy().astype(np.uint8))sitk.WriteImage(warp_img, out_path)print("end")

需要一提的是，由于不能精确的决定网格图片的大小，所以在生成的时候多生成了10个像素，然后在读取的时候再进行裁剪。改变上述代码中第57行的除数（10）可以改变网格线的密集程度。

感觉这种方式更优雅一点，但是效果和第二种是差不多的，第二种和第三种方法的效果如下图所示：