19 | 如何可视化CT影像数据，2d和3d图像

机器学习之禅

发布于 2022-07-11 15:50:29

1.4K0

发布于 2022-07-11 15:50:29

这几天研究的跟PyTorch没啥关系，跟深度学习也没啥关系，但是跟做项目关系很大，那就是怎么把CT影像可视化。要完美复现代码，需要安装的扩展包simpleITK，ipyvolume，diskcache，cassandra-driver

安装扩展包推荐使用conda install package-you-need，比如说安装ipyvolume，就用

conda install ipyvolume

但是conda经常非常慢，卡在solving environment阶段，或者是找不到合适的包，这个时候可以用pip来安装。在disk.py文件里，可能代码比较老了，有个引用的包失效了，这里需要修改一下，不然这一句会报错，from diskcache.core import BytesType, MODE_BINARY, BytesIO，修改后的代码如下

import gzipfrom diskcache import FanoutCache, Diskfrom cassandra.cqltypes import BytesTypefrom diskcache import FanoutCache, Disk,corefrom diskcache.core import iofrom io import BytesIOfrom diskcache.core import MODE_BINARYfrom util.logconf import logging
log = logging.getLogger(__name__)# log.setLevel(logging.WARN)log.setLevel(logging.INFO)# log.setLevel(logging.DEBUG)

这里构建一个展示图像的py文件，名字是vis.py，首先是引用信息

import matplotlib
matplotlib.use('nbagg')  #代码在远程执行，而在本地绘图，在我们本机上貌似没有明显效果import numpy as npimport matplotlib.pyplot as plt  #使用matplotlib绘图from testproject.dsets import Ct, LunaDataset

然后是一个全局变量

clim=(-1000.0, 300) #用来设置展示图像的颜色上下限

定义两个方法，第一个方法是获取正样本

def findPositiveSamples(start_ndx=0, limit=100):  #获取样本数据
    ds = LunaDataset() #从我们之前构建的LunaDataset实例化
    positiveSample_list = []
    for sample_tup in ds.candidateInfo_list: #从小圆点集合读取数据
        if sample_tup.isNodule_bool: #如果是结节
            print(len(positiveSample_list), sample_tup) #输出样本集现在的长度，并展示样本
            positiveSample_list.append(sample_tup) #把样本数据加入样本集

        if len(positiveSample_list) >= limit: #如果长度超出我们的限制跳出，这里开头限制是100
            break

    return positiveSample_list

第二个方法是绘图，主要的意思我已经都写在代码注释里了

#第二个方法是绘制图像，输入的是一个CT图像uid，以及块索引#**kwargs 允许你将不定长度的键值对, 作为参数传递给一个函数。如果你想要在一个函数里处理带名字的参数, 你应该使用**kwargs。def showCandidate(series_uid, batch_ndx=None, **kwargs):
    ds = LunaDataset(series_uid=series_uid, **kwargs)
    pos_list = [i for i, x in enumerate(ds.candidateInfo_list) if x.isNodule_bool]#处理块索引信息，如果没有输入，就取正样本list的第一个作为索引；如果正样本list是空的，就设为0
    if batch_ndx is None:
        if pos_list:
            batch_ndx = pos_list[0]
        else:
            print("Warning: no positive samples found; using first negative sample.")
            batch_ndx = 0

    ct = Ct(series_uid) #根据uid获取CT影像数据
    ct_t, pos_t, series_uid, center_irc = ds[batch_ndx]
    ct_a = ct_t[0].numpy()

    fig = plt.figure(figsize=(30, 50))  #设置图像尺寸#这个group信息是用来取切片位置，这里取了9个切片，你可以试试把里面的值换一下
    group_list = [
        [9, 11, 13],
        [15, 16, 17],
        [19, 21, 23],
    ]#第一个子图 #add_subplot的参数(3,4,9)意思是绘制在一个三行四列的子图的第9个位置。
    subplot = fig.add_subplot(len(group_list) + 2, 3, 1) #这里len(group_list)是3，也就是add_subplot(5,3,1)，绘制在5行三列的图像第1个位置
    subplot.set_title('index {}'.format(int(center_irc[0])), fontsize=30)  #设置子图像的标题    for label in (subplot.get_xticklabels() + subplot.get_yticklabels()): #设置子图像的轴信息
        label.set_fontsize(20)
    plt.imshow(ct.hu_a[int(center_irc[0])], clim=clim, cmap='gray') #取整个CT数据的index信息展示，对应我们的俯视图

    subplot = fig.add_subplot(len(group_list) + 2, 3, 2)
    subplot.set_title('row {}'.format(int(center_irc[1])), fontsize=30)
    for label in (subplot.get_xticklabels() + subplot.get_yticklabels()):
        label.set_fontsize(20)
    plt.imshow(ct.hu_a[:,int(center_irc[1])], clim=clim, cmap='gray')#取行信息展示，对应我们的正视图
    plt.gca().invert_yaxis()

    subplot = fig.add_subplot(len(group_list) + 2, 3, 3)
    subplot.set_title('col {}'.format(int(center_irc[2])), fontsize=30)
    for label in (subplot.get_xticklabels() + subplot.get_yticklabels()):
        label.set_fontsize(20)
    plt. imshow(ct.hu_a[:,:,int(center_irc[2])], clim=clim, cmap='gray')#取列信息展示，对应我们的侧视图
    plt.gca().invert_yaxis()

    subplot = fig.add_subplot(len(group_list) + 2, 3, 4)
    subplot.set_title('index {}'.format(int(center_irc[0])), fontsize=30)
    for label in (subplot.get_xticklabels() + subplot.get_yticklabels()):
        label.set_fontsize(20)
    plt.imshow(ct_a[ct_a.shape[0]//2], clim=clim, cmap='gray')#这个是取候选小块的信息，除2是放大？

    subplot = fig.add_subplot(len(group_list) + 2, 3, 5)
    subplot.set_title('row {}'.format(int(center_irc[1])), fontsize=30)
    for label in (subplot.get_xticklabels() + subplot.get_yticklabels()):
        label.set_fontsize(20)
    plt.imshow(ct_a[:,ct_a.shape[1]//2], clim=clim, cmap='gray')#候选小块的行信息
    plt.gca().invert_yaxis()

    subplot = fig.add_subplot(len(group_list) + 2, 3, 6)
    subplot.set_title('col {}'.format(int(center_irc[2])), fontsize=30)
    for label in (subplot.get_xticklabels() + subplot.get_yticklabels()):
        label.set_fontsize(20)
    plt.imshow(ct_a[:,:,ct_a.shape[2]//2], clim=clim, cmap='gray') #候选小块的列信息
    plt.gca().invert_yaxis()#这里是绘制group的9张切片图
    for row, index_list in enumerate(group_list): 
        for col, index in enumerate(index_list):
            subplot = fig.add_subplot(len(group_list) + 2, 3, row * 3 + col + 7)
            subplot.set_title('slice {}'.format(index), fontsize=30)
            for label in (subplot.get_xticklabels() + subplot.get_yticklabels()):
                label.set_fontsize(20)
            plt.imshow(ct_a[index], clim=clim, cmap='gray')


    print(series_uid, batch_ndx, bool(pos_t[0]), pos_list)

接下来看一下效果。打开一个Jupyter notebook。

%matplotlib inline  #用于指定matplotlib绘制的图像嵌入到notebook中from testproject.vis import findPositiveSamples, showCandidate #引入两个方法positiveSample_list = findPositiveSamples() #获取正样本

这个时候的输出很长，输出信息是100条正样本

然后随便找一条数据，比如说第11条，取出它的uid，皇天不负有心人，可以看到很成功的输出了图像。图像分五行，跟我们之前的设定是一样的。

series_uid = positiveSample_list[11][2]
showCandidate(series_uid)

我们再看一下第18个，展示效果很不错。

series_uid = positiveSample_list[18][2]
showCandidate(series_uid)

除了这种平面可视化，因为我们的数据是立体数据，这里还有一种方式可以把它展示成一个3D动态效果图。这里就需要用到我们上面提到的ipyvolume扩展包。ipyvolume是一个支持在网页渲染WebGL图像的包，可以让我们在Jupyter Notebook中渲染3D图像。详细的项目情况可以参考这里GitHub - widgetti/ipyvolume: 3d plotting for Python in the Jupyter notebook based on IPython widgets using WebGL

先用ipyvolume画一个立方体看看。这是一个两层嵌套的立方体

import numpy as np
import ipyvolume as ipvV = np.zeros((128,128,128)) # our 3d array# outer boxV[30:-30,30:-30,30:-30] = 0.75V[35:-35,35:-35,35:-35] = 0.0# inner boxV[50:-50,50:-50,50:-50] = 0.25V[55:-55,55:-55,55:-55] = 0.0ipv.quickvolshow(V, level=[0.25, 0.75], opacity=0.03, level_width=0.1, data_min=0, data_max=1)

点击执行之后可能看不到图像，你不要担心，这个事情我也研究明白了，WebGL需要开启渲染加速才能使用，而且有些浏览器是不支持的，我这里用的是chrome浏览器，需要配置一下才能启动。在chrome执行下面的步骤，

在浏览器地址栏输入 chrome://flags
找到 Override software rendering 列表选择 Enabled
点击 Relaunch Now 按钮，重启浏览器然后再执行代码就可以看到图像了，图像上方和下方有一些参数可以调节，下面的图像会发生实时的变化。

图像看起来挺丑的，但是好处是我们可以用鼠标旋转它。

接下来绘制一下CT数据

ct = getCt(series_uid)
ipv.quickvolshow(ct.hu_a, level=[0.25, 0.5, 0.9], opacity=0.1, level_width=0.1, data_min=-1000, data_max=1000)

这次得到的是一个绿色的柱体，如果你不仔细看都没办法发现这是一个人体，这个快速展示的效果真的很不好看，尤其是加上这个让人难受的绿色。看起来就像一个长满青苔的木桩。

为了让它能更好的展示我们的数据，这里还需要一些复杂的操作，这里又用到一个包 scipy.ndimage.morphology，使用这个包对我们的数据进行了形态学变换，你可以理解成进行了一些基础的图像处理，通过这些处理能够把不同密度的组织数据提取出来，如果你对这个包感兴趣可以挨个查一查都是什么意思。

import scipy.ndimage.morphologydef build2dLungMask(ct, mask_ndx, threshold_gcc = 0.7):
    dense_mask = ct.hu_a[mask_ndx] > threshold_gcc
    denoise_mask = scipy.ndimage.morphology.binary_closing(dense_mask, iterations=2)
    tissue_mask = scipy.ndimage.morphology.binary_opening(denoise_mask, iterations=10)
    body_mask = scipy.ndimage.morphology.binary_fill_holes(tissue_mask)
    air_mask = scipy.ndimage.morphology.binary_fill_holes(body_mask & ~tissue_mask)

    lung_mask = scipy.ndimage.morphology.binary_dilation(air_mask, iterations=2)

    return air_mask, lung_mask, dense_mask, denoise_mask, tissue_mask, body_maskdef build3dLungMask(ct):
    air_mask, lung_mask, dense_mask, denoise_mask, tissue_mask, body_mask = mask_list = \        [np.zeros_like(ct.hu_a, dtype=np.bool) for _ in range(6)]

    for mask_ndx in range(ct.hu_a.shape[0]):
        for i, mask_ary in enumerate(build2dLungMask(ct, mask_ndx)):
            mask_list[i][mask_ndx] = mask_ary    return air_mask, lung_mask, dense_mask, denoise_mask, tissue_mask, body_mask

接着我们把CT数据拿过来，这里我们只绘制骨头和肺部

from testproject.dsets import getCt
ct = getCt(series_uid)air_mask, lung_mask, dense_mask, denoise_mask, tissue_mask, body_mask = build3dLungMask(ct)bones = ct.hu_a * (ct.hu_a > 1.5)lungs = ct.hu_a * air_mask
ipv.figure()ipv.pylab.volshow(bones + lungs, level=[0.17, 0.17, 0.23], data_min=100, data_max=900)ipv.show()

初始还是一个俯视图，

我们给它转到正面，可以比较清楚的看到骨头，但是感觉看不到肺部啊，马马虎虎，但是总体比上面的绿色苔藓好多了。

关于数据可视化估计还有很多其他的方法和包，在做项目的过程中使用可视化方法可以有效帮助我们检测数据的错误，因为可以直接用肉眼看嘛，有问题一下就看出来了，比如说加载出来的图像不符合我们的认知，或者压根加载不出图像等等情况。再一个，图像有助于帮助看不懂代码的人了解算法效果，比如说产品运营，拿来做项目汇报也是不错的。

这两天终于把可视化的代码捣鼓清楚了，明天继续看书。

本文参与腾讯云自媒体同步曝光计划，分享自微信公众号。

原始发表：2022-06-11，如有侵权请联系 cloudcommunity@tencent.com 删除

webgl