LUNA16结节可视化及世界（CT)坐标与体素（图像）坐标的转换

 这个是关于DCM的医学坐标系的一些背景知识，本文将在此基础上做些补充。

DICOM医学图像读取涉及到的医学坐标体系_sunyao_123的博客-CSDN博客_dicom坐标系

关于坐标名称，看了许多博文，有像素坐标，真实坐标，体素坐标，图像坐标，世界坐标搞得一头雾水。本文中CT坐标指的是设备拍摄时的坐标，单位为mm。图像坐标指的是转化为numpy后的坐标，单位为单位1。自己的理解世界坐标应该是CT坐标，体素坐标，像素坐标，真实坐标及图像坐标应该都是图像坐标。

一 、 肺部图像的转化

首先，CT图像的坐标原点是左上角，和将CT图像转为numpy格式的坐标原点一致，方便处理。

其次，CT图像的像素间距与像素原点每个case都不同具体和拍摄设备参数相关，会将这些参数存储到mha文件中。以下图为例，此case中的CT像素间隔和原点坐标分别为

Spacing[x ,y, z] = [0.73242199, 0.73242199, 1.25] 

Origin[x ,y, z] = [-163.199997, -145.5, -380.] 单位都是mm。而像素个数为512*512*305也就是转成numpy后，数组的shape为[305, 512, 512]。

我的理解是在图像坐标中图像的像素间隔为单位1，CT坐标下像素间隔为spacing，真实的图像大小应为转为numpy的shape乘上spacing。

new_size = data.shape * numpySpacing[::-1]  data.shape的输出顺序为（z, y, x)故相乘时需要将numpySpacing也反序。new_size为CT图像的真实大小单位为mm。

在可视化肺结节时我们使用的是图像坐标，故需要将肺结节的坐标也转化为图像坐标。

itkImage = sitk.ReadImage('D:/BaiduNetdiskDownload/LUNA16/subset9/1.3.6.1.4.1.14519.5.2.1.6279.6001.153985109349433321657655488650.mhd')
    numpyOrigin = np.array(list(itkImage.GetOrigin()))
    numpySpacing = np.array(list(itkImage.GetSpacing()))
    print(numpyOrigin,numpySpacing)
    data = sitk.GetArrayFromImage(itkImage)
    print(data.shape)
    new_size = data.shape * numpySpacing[::-1]
    print(new_size)
# 输出结果
[-163.199997 -145.5      -380.      ] [0.73242199 0.73242199 1.25      ]
(305, 512, 512)
[381.25       375.00006104 375.00006104]

二、肺结节的坐标转化

 以上述case为例，肺结节的坐标信息如下。分别是ID,x,y,z,半径。单位是mm。

 转化的公式：图像坐标 = CT坐标系下的（结节位置 - 原点位置） / 像素间隔

def worldToVoxelCoord(worldCoord, origin, spacing):
    # 图像坐标系的结节位置 = CT坐标系下的（结节位置 - 原点位置） / 体素间隔
    stretchedVoxelCoord = np.absolute(worldCoord - origin)
    voxelCoord = stretchedVoxelCoord / spacing
    return voxelCoord

 三、可视化

显示代码转载自： LUNA16数据集（二）肺结节可视化 - 走看看

def show_nodules(ct_scan, nodules, radius=20, pad=5, max_show_num=3):  # radius是正方形边长一半，pad是边的宽度,max_show_num最大展示数

    show_index = []
    for idx in range(nodules.shape[0]):  # lable是一个nx4维的数组，n是肺结节数目，4代表x,y,z,以及直径
        if idx < max_show_num:
            if abs(nodules[idx, 0]) + abs(nodules[idx, 1]) + abs(nodules[idx, 2]) + abs(nodules[idx, 3]) == 0: continue

            x, y, z = int(nodules[idx, 0]), int(nodules[idx, 1]), int(nodules[idx, 2])

            data = ct_scan[z]

            # 注意 y代表纵轴，x代表横轴
            data[max(0, y - radius):min(data.shape[0], y + radius),
            max(0, x - radius - pad):max(0, x - radius)] = 3000  # 竖线
            data[max(0, y - radius):min(data.shape[0], y + radius),
            min(data.shape[1], x + radius):min(data.shape[1], x + radius + pad)] = 3000  # 竖线
            data[max(0, y - radius - pad):max(0, y - radius),
            max(0, x - radius):min(data.shape[1], x + radius)] = 3000  # 横线
            data[min(data.shape[0], y + radius):min(data.shape[0], y + radius + pad),
            max(0, x - radius):min(data.shape[1], x + radius)] = 3000  # 横线

            if z in show_index:  # 检查是否有结节在同一张切片，如果有，只显示一张
                continue
            show_index.append(z)
            plt.figure
            plt.imshow(data, cmap='gray')

    plt.show()

数据准备代码： 

itkImage = sitk.ReadImage('D:/BaiduNetdiskDownload/LUNA16/subset9/1.3.6.1.4.1.14519.5.2.1.6279.6001.153985109349433321657655488650.mhd')
    numpyOrigin = np.array(list(itkImage.GetOrigin()))
    numpySpacing = np.array(list(itkImage.GetSpacing()))
    print(numpyOrigin,numpySpacing)
    data = sitk.GetArrayFromImage(itkImage)
    print(data.shape)
    new_size = data.shape * numpySpacing[::-1]
    print(new_size)
    annos = np.array(pandas.read_csv('D:/BaiduNetdiskDownload/LUNA16/CSVFILES/annotations.csv'))
    # 图像相应的结节信息
    annos = annos[198,1:]
    # CT坐标转为图像坐标
    annos1 = worldToVoxelCoord(annos[:3],numpyOrigin,numpySpacing)
    # 结节半径也转为图像坐标
    annos = np.concatenate([annos1,[annos[3]/numpySpacing[1]]])
    annos = np.expand_dims(annos,axis=0)
    print(annos.shape)
    show_nodules(data,annos)

结果显示：