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

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

首先来介绍一下形变场，一个大小为[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 np


def grid2contour(grid, title):
    '''
    grid--image_grid used to show deform field
    type: numpy ndarray, shape： (h, w, 2), value range：(-1, 1)
    '''
    assert grid.ndim == 3
    x = 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 line
    Z1 = Z1[::-1]  # vertical flip
    Z2 = grid[:, :, 1] + 2

    plt.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 ticks
    plt.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_norm
    grid2contour(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_arr
    grid2contour(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 plt

import 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 grid
        vectors = [torch.arange(0, s) for s in size]
        grids = torch.meshgrid(vectors)
        grid = torch.stack(grids)  # y, x, z
        grid = torch.unsqueeze(grid, 0)  # add batch
        grid = grid.type(torch.FloatTensor)
        self.register_buffer('grid', grid)

        self.mode = mode

    def 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 + flow

        shape = flow.shape[2:]

        # Need to normalize grid values to [-1, 1] for resampler
        for 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,1
            new_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）可以改变网格线的密集程度。

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