CT片居然可以这么玩：用头部CT断层扫描片复原三维头像

太过爱你忘了你带给我的痛 2022-09-02 03:41 331阅读 0赞

CT是现代医学影像的主力设备，寻常百姓并不陌生。通常，一张CT片由多张连续断层扫描的图像组成。在医生眼中，CT片展示了人体器官的形态和性质，是判断病人健康状况的重要依据。对于普通人而言，CT片则像天书，几无大用。不过呢，总有不甘寂寞喜欢折腾的程序员，把看似无用的CT片玩出了新花样——用头部CT断层扫描片成功复原出了逼真的三维头像。

CT片是如何变成三维头像的呢？说白了，其实很简单，只需要3步：

1.  准备一套头部CT断层扫描片。怎么准备？要么去网上找（我就是从网上下载的），要么自己扫描胶片。
2.  按照断层顺序逐张读入图片，并垂直堆叠成一个numpy数组。
3.  导入3D工具库WxGL，调用capsule函数，绘制三维头像。

我手头这套头部CT断层扫描片共有109张，是侧卧位的水平断层扫描，从一侧开始至另一侧结束。下图从左至右分别是第11张、第21张、第55张、第89张、第99张。不难看出，中间的一张面积最大，是鼻梁正中的断层图片，鼻尖朝向图片上方，越往两侧的图片，面积越小。

![在这里插入图片描述][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3h1Zml2ZQ_size_16_color_FFFFFF_t_70_pic_center]  
       有了CT断层扫描片，接下来就是打开文件读取数据了。做这一步之前需要先导入pillow模块和numpy模块。假定CT文件保存在D:\\XufiveGithub\\wxgl\\res\\headCT这个路径下，每个文件名都以head+序号的方式命名，下面是读取数据的代码。

>>> from PIL import Image
    >>> import numpy as np
    >>> base_dir = r'D:\XufiveGithub\wxgl\res\headCT'
    >>> data = list()
    >>> for i in range(109):
    	im = np.array(Image.open(r'%s\head%d.png'%(base_dir, i)))
    	data.append(im)
    
    >>> data = np.stack(data, axis=0)

其实，读取数据的代码可以写得更简洁，使用列表推导式，只用一行就可以了。

>>> data = np.stack([np.array(Image.open(r'%s\head%d.png'%(base_dir, i))) for i in range(109)], axis=0)
    >>> data.shape
    (109, 256, 256, 4)

从data.shape可以看出，总共有109个断层，每个断层分辨率是256x256像素，每个像素有RGBA四个通道。比较每个通道的数据，就会发现，每个断层的RGBA四个通道是完全相同的，这是胶片扫描的特点决定的。

三维重建的基本原理是滤除完全透明的像素，用具有相同透明度且位置相邻的像素构建surface，因此数组data中只有透明通道（A通道）是我们需要的。

data = data[...,3] # 提取透明通道
    >>> data.shape
    (109, 256, 256)

现在data变成一个三维数组了。数组的0轴，即第1张到第109张CT片的方向，对应着头像的左右；数组的1轴，即每张CT片自上而下的方向，对应着头像的前后；数组的2轴，即每张CT片从左至右的方向，对应着头像的上下。3D工具库WxGL习惯将数据数组的0轴视为高度方向，不做调整的话，绘制出的头像将是水平侧卧的，就像CT片展示的那样。为了让头像立起来，需要将头像的上下方对应的数组2轴翻转到0轴之后，再上下翻转。

>>> data = np.rollaxis(data, 2, 0) # 2轴翻转至0轴，让头像立起来（倒立）
    >>> data = np.flipud(data) # 上下翻转，头像正立
    >>> data.shape
    (256, 109, 256)

至此，数据处理算是完成了，不过，这些数据还没有空间定位，只是保持了相对的位置关系，我们需要给数据限定一个空间位置。由于CT片没有给出头部尺寸和分层厚度，只能按照通常的头像比例，给出头像在xyz轴上的空间范围。

x = np.linspace(-1, 1, data.shape[2]) # 头像前后范围：[-1, 1]
    y = np.linspace(-0.65, 0.65, data.shape[1]) # 头像左右范围：[-0.65, 0.65]
    z = np.linspace(-1, 1, data.shape[0]) # 头像高度范围：[-1, 1]

另外，还要滤除几乎完全透明的那些像素。我们暂且以最大值的1/8作为阈值，调整该阈值将会得到不同的重建效果。

level = data.max()/8 # 忽略低于此阈值的数据

万事俱备，只欠东风。是时候让wxgl一显身手了。

plt.capsule(data, x=x, y=y, z=z, level=level, color='#EEE9D4')
    plt.show()

不要惊讶，只需要两行代码，一个逼真的三维头像呈现在眼前。这一切，仅仅依赖一套CT断层扫描片和WxGL这个3D工具库。

![在这里插入图片描述][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3h1Zml2ZQ_size_16_color_FFFFFF_t_70_pic_center 1]  
       刚刚发布的0.7.3版3D工具库WxGL，新增了几个有趣的功能，其中之一就是刚刚用到的capsule函数。capsule意为胶囊，顾名思义，capsule函数可以从三维数据中找出等值面，并尝试构建出囊性结构（如果存在的话）。WxGL的安装非常简单，直接使用pip命令即可。如果在安装和使用WxGL的过程中遇到问题，请在文后留言，也可以直接到[我的GitHub][GitHub]反馈问题。

> pip install wxgl

最后，附上用头部CT断层扫描片复原三维头像的完整代码。

import numpy as np
    from PIL import Image
    from wxgl import wxplot as plt
    
    base = r'D:\XufiveGithub\wxgl\res\headCT'
    data = np.stack([np.array(Image.open(r'%s\head%d.png'%(base, i))) for i in range(109)], axis=0)
    
    data = data[...,3] # 提取透明通道
    data = np.rollaxis(data, 2, 0) # 2轴翻转至0轴，让头像立起来（倒立）
    data = np.flipud(data) # 上下翻转，头像正立
    
    x = np.linspace(-1, 1, data.shape[2]) # 头像前后范围：[-1, 1]
    y = np.linspace(-0.65, 0.65, data.shape[1]) # 头像左右范围：[-0.65, 0.65]
    z = np.linspace(-1, 1, data.shape[0]) # 头像高度范围：[-1, 1]
    level = data.max()/8 # 忽略低于此阈值的数据
    
    plt.capsule(data, x=x, y=y, z=z, level=level, color='#EEE9D4')
    plt.show(rotate=-1)

下图直观显示了不同的阈值对于三维重建结果的影响。

![在这里插入图片描述][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3h1Zml2ZQ_size_16_color_FFFFFF_t_70_pic_center 2]

[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3h1Zml2ZQ_size_16_color_FFFFFF_t_70_pic_center]: /images/20220829/978a7cba8a3d4053b26a054963fec379.png
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3h1Zml2ZQ_size_16_color_FFFFFF_t_70_pic_center 1]: /images/20220829/2c8a1254c6aa485e8eab019ffd723e0f.png
[GitHub]: https://github.com/xufive/wxgl
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3h1Zml2ZQ_size_16_color_FFFFFF_t_70_pic_center 2]: /images/20220829/19a23dc0fcca462eb59a4ea3eb80c8ac.png