【3D 图像分割】基于 Pytorch 的 VNet 3D 图像分割7(数据预处理)

这篇具有很好参考价值的文章主要介绍了【3D 图像分割】基于 Pytorch 的 VNet 3D 图像分割7(数据预处理)。希望对大家有所帮助。如果存在错误或未考虑完全的地方,请大家不吝赐教,您也可以点击"举报违法"按钮提交疑问。

在上一节:【3D 图像分割】基于 Pytorch 的 VNet 3D 图像分割6(数据预处理) 中,我们已经得到了与mhd图像同seriesUID名称的mask nrrd数据文件了,可以说是一一对应了。

并且,mask的文件,还根据结节被多少人同时标注,区分成了4个文件夹,分别是标注了一、二、三、四次,一共就4个医生参与标注。

再加上官方已经给整理好的肺实质分割的文件,我们就获得了以下这些数据:

  1. ct图像数据;
  2. 肺实质分割数据;
  3. 包含结节位置的mask数据。

那本篇做什么呢?看下面的导言。

一、导言

上述得到的这些,就满足了我们的需求了,都是一一对应的,无论是后续的数据预处理,还是拿过来用于训练,都非常的方便。

但是呢,对于原始的ct数据,他在Z轴上的层厚是不同的,这点可以在dicom文件里面看到,也可以在mhd文件的查询到关于层厚的信息。在这点上,不同的序列,差异是非常大的。表现在一个3维数组的结节上面,在这个维度上就是被压扁,和拉长的样子。

xy方向,其实也是存在spacing的差异的,但是这种差异没有像z轴那么夸张的,这里可以选择处理和不处理均可(有些论文进行了处理,有些没有。默认都是512x512大小,resample后会变小)。

至此,本篇的目的就很明确了,是要做下面几件事:

  1. 对原始图像进行肺实质提取,将肺区外的部分进行裁剪,或者改为固定像素值;
  2. 对图像和结节mask进行resample操作,本篇是zyx均进行resample1mm

二、具体实施

怎么做的部分,我们分三部分:

  1. 肺实质裁剪
  2. imagenodule mask进行resample操作
  3. 获取结节中心点坐标和半径

下面就一一展开

2.1、主函数部分

由于这部分数据量比较多,所以在主函数部分采用了多进程的模式,加快处理速度。需要读进来的数据也就是前面篇章已经处理好的,这里都可以直接使用。

下面就是主函数

import sys
import numpy as np
import scipy.ndimage
from skimage import measure, morphology
import SimpleITK as sitk
from multiprocessing import Pool
import os
import nrrd

###############################################################################
# 将标记的mask,和ct原图,加入左右肺区分割的图像,生成去除noise的,剩下肺区的ct和结节mask
###############################################################################
def main():
    n_consensus = 4
    do_resample = True
    img_dir = './LUNA16/image_combined'
    lung_mask_dir = './LUNA16/seg-lungs-LUNA16'
    nod_mask_dir = os.path.join('./LUNA16/nodule_masks', str(n_consensus))

    save_dir = os.path.join('./LUNA16/preprocessed', str(n_consensus))
    os.makedirs(save_dir, exist_ok=True)

    params_lists = []
    # 多进程处理
    for nrrd_name in os.listdir(nod_mask_dir):
        #                         seg-lungs-LUNA16, masks_test/3, seg-lungs-LUNA16, preprocessed_test/3, True
        pid = nrrd_name[:-5]
        params_lists.append([pid, lung_mask_dir, nod_mask_dir, img_dir, save_dir, do_resample])

    pool = Pool(processes=4)
    pool.map(cropResample_process, params_lists)

    pool.close()
    pool.join()

    pool = Pool(processes=4)
    pool.map(generateBBoxes_label, params_lists)

    pool.close()
    pool.join()

if __name__ == '__main__':
    main()

有两个部分,

  • cropResample_process:和名称一样,进行肺实质的cropresample操作;
  • generateBBoxes_label:将处理完毕的结节mask,得到结节中心的坐标和半径。

2.2、肺实质裁剪

这小块的步骤,大概如下:

  1. 首先,就是数据读取,这部分的详细介绍,可以参考我之前的这篇文章:【医学影像数据处理】nii、npz、npy、dcm、mhd 的数据格式互转,及多目标分割处理汇总
  2. 其次,就是将hu值,转化为0-255的值,也就是函数HU2uint8(),对于这部分,可以参考hu值是如何转为0-255的可视化部分的介绍:【医学影像数据处理】 Dicom 文件格式处理汇总
  3. 另外,就是将肺区mask作用到图像上,肺实质外采用pad valud补充
  4. 最后,将处理好的image、mask和相关参数存储到本地

代码如下,就该说明的部分都进行注释,相信能轻易看懂。

def load_itk_image(filename):
    """Return img array and [z,y,x]-ordered origin and spacing
    """
    itkimage = sitk.ReadImage(filename)
    numpyImage = sitk.GetArrayFromImage(itkimage)

    numpyOrigin = np.array(list(reversed(itkimage.GetOrigin())))
    numpySpacing = np.array(list(reversed(itkimage.GetSpacing())))

    return numpyImage, numpyOrigin, numpySpacing

def HU2uint8(image, HU_min=-1200.0, HU_max=600.0, HU_nan=-2000.0):
    """
    Convert HU unit into uint8 values. First bound HU values by predfined min
    and max, and then normalize
    image: 3D numpy array of raw HU values from CT series in [z, y, x] order.
    HU_min: float, min HU value.
    HU_max: float, max HU value.
    HU_nan: float, value for nan in the raw CT image.
    """
    image_new = np.array(image)
    image_new[np.isnan(image_new)] = HU_nan

    # normalize to [0, 1]
    image_new = (image_new - HU_min) / (HU_max - HU_min)
    image_new = np.clip(image_new, 0, 1)
    image_new = (image_new * 255).astype('uint8')

    return image_new

def convex_hull_dilate(binary_mask, dilate_factor=1.5, iterations=10):
    """
    Replace each slice with convex hull of it then dilate. Convex hulls used
    only if it does not increase area by dilate_factor. This applies mainly to
    the inferior slices because inferior surface of lungs is concave.
    binary_mask: 3D binary numpy array with the same shape of the image,
        that only region of interest is True. One side of the lung in this
        specifical case.
    dilate_factor: float, factor of increased area after dilation
    iterations: int, number of iterations for dilation
    return: 3D binary numpy array with the same shape of the image,
        that only region of interest is True. Each binary mask is ROI of one
        side of the lung.
    """
    binary_mask_dilated = np.array(binary_mask)
    for i in range(binary_mask.shape[0]):
        slice_binary = binary_mask[i]

        if np.sum(slice_binary) > 0:
            slice_convex = morphology.convex_hull_image(slice_binary)

            if np.sum(slice_convex) <= dilate_factor * np.sum(slice_binary):
                binary_mask_dilated[i] = slice_convex

    struct = scipy.ndimage.generate_binary_structure(3, 1)
    binary_mask_dilated = scipy.ndimage.binary_dilation(
        binary_mask_dilated, structure=struct, iterations=10)

    return binary_mask_dilated

def apply_lung_mask(image, binary_mask1, binary_mask2, pad_value=170,
                    bone_thred=210, remove_bone=False):
    """
    Apply the binary mask of each lung to the image. Regions out of interest
    are replaced with pad_value.
    image: 3D uint8 numpy array with the same shape of the image.
    binary_mask1: 3D binary numpy array with the same shape of the image,
        that only one side of lung is True.
    binary_mask2: 3D binary numpy array with the same shape of the image,
        that only the other side of lung is True.
    pad_value: int, uint8 value for padding image regions that is not
        interested.
    bone_thred: int, uint8 threahold value for determine parts of image is
        bone.
    return: D uint8 numpy array with the same shape of the image after
        applying the lung mask.
    """
    binary_mask = binary_mask1 + binary_mask2
    binary_mask1_dilated = convex_hull_dilate(binary_mask1)
    binary_mask2_dilated = convex_hull_dilate(binary_mask2)
    binary_mask_dilated = binary_mask1_dilated + binary_mask2_dilated
    binary_mask_extra = binary_mask_dilated ^ binary_mask

    # replace image values outside binary_mask_dilated as pad value
    image_new = image * binary_mask_dilated + pad_value * (1 - binary_mask_dilated).astype('uint8')

    # set bones in extra mask to 170 (ie convert HU > 482 to HU 0;
    # water).
    if remove_bone:
        image_new[image_new * binary_mask_extra > bone_thred] = pad_value

    return image_new

def cropResample_process(params):
    #    seg-lungs-LUNA16, masks_test/3, seg-lungs-LUNA16, preprocessed_test/3, True
    pid, lung_mask_dir, nod_mask_dir, img_dir, save_dir, do_resample = params

    print('Preprocessing %s...' % (pid))

    img_org, origin, spacing = load_itk_image(os.path.join(img_dir, '%s.mhd' % (pid)))
    lung_mask, _, _ = load_itk_image(os.path.join(lung_mask_dir, '%s.mhd' % (pid)))
    nodule_mask, _ = nrrd.read(os.path.join(nod_mask_dir, '%s.nrrd' % (pid)))

    # 4-右肺   3-左肺   5-气管
    binary_mask_r = lung_mask == 4
    binary_mask_l = lung_mask == 3
    binary_mask = binary_mask_r + binary_mask_l

    img_org = HU2uint8(img_org)
    img_lungRL = apply_lung_mask(img_org, binary_mask_r, binary_mask_l)

有一个点前面从没有说明过,那就是官方提供的lung mask数组,在这里简要的记录下:

  • 数字3,表示左肺
  • 数字4,表示右肺
  • 数字5,表示气管

还是第一次看到这个按位异或运算符(^),简单的学习了下:

按位异或运算符(^)用于将两个操作数的每个对应位进行逻辑异或操作。如果两个对应位的值相同,则结果为0,否则为1。异或的本质是没有进位的加法。

用真值表来看一下异或的运算逻辑:
【3D 图像分割】基于 Pytorch 的 VNet 3D 图像分割7(数据预处理),人工智能(AI)医学影像,3d,pytorch,人工智能,LUNA16,LIDC- IDRI
看起来非常直观并且简单明了 但是如果我们将布尔值之间的异或换成数字之间的异或会发生什么???让我们来试一下:

0 ^ 0
0 ^ 1
1 ^ 0
1 ^ 1

>>>
0
1
1
0

结果告诉我们: 数字相同异或值为0 ;数字不相同异或值为1。

让我们再试试0,1除外的数字:

5 ^ 3
>>>6

卧槽?为什么答案会是6,不应该是1么??这背后发生了什么?

异或是基于二进制基础上按位异或的结果 5 ^ 3 的过程 其实是将5和3分别转换为二进制:

5 = 0101(b)

3 = 0011(b)

按位异或:

0^0 ->0

1^0 ->1

0^1 ->1

1^1 ->0

排起来就是0110(b) 转换为十进制:6。更多异或内容,可参考这里:浅谈Python逻辑运算符 异或xor


dilate膨胀后的binary mask和原始的binary mask求异或运算,对应位的值相同,结果为0,否则为1。那么,得到的结果也就是膨胀出来的那部分,就是bone,这部分在去除bone阶段使用到。

可能会有这样的疑问:为什么不直接imagelung mask相乘,得到一个分割肺实质后留下来的image呢?反而需要采用凸包优化的方式,多此一举呢?

:在lung mask里面,肺实质的分割是有误差的。也就是肺实质的分割是沿着肺区边缘的,但是某些结节的位置,恰好在肺区的边界上,且密度很大。那么mask就会呈现一个内凹的一个状态。如果采用上面的方法,这样结节就被抠除了。采用凸包优化,就可以利用稍微扩展肺实质边缘,达到将更多肺区留下来的效果。

但是,对于肺结核等等大病灶的疾病,采用上述取出肺实质的方法就不行。主要是因为肺结核的病种范围比较大,尽管采用了凸包优化,最终还是会切除很大一块肺区位置,这样肺区就不完整了,有些得不偿失。

下面是skimage.morphology.convex_hull_image官方给出的实例,如下:点击直达
【3D 图像分割】基于 Pytorch 的 VNet 3D 图像分割7(数据预处理),人工智能(AI)医学影像,3d,pytorch,人工智能,LUNA16,LIDC- IDRI作用到我们项目里面,切割后的样子如下:

【3D 图像分割】基于 Pytorch 的 VNet 3D 图像分割7(数据预处理),人工智能(AI)医学影像,3d,pytorch,人工智能,LUNA16,LIDC- IDRI

2.3、resample操作

本篇对resample的操作,在zyx的各个维度上,就雨露均沾,通通调整到1mm的状态,这样得到的一个像素大小,表示的也就是物理大小,不会引起任何一个维度上变形的情况。

代码如下所示:

def resample(image, spacing, new_spacing=[1.0, 1.0, 1.0], order=1):
    """
    Resample image from the original spacing to new_spacing, e.g. 1x1x1
    image: 3D numpy array of raw HU values from CT series in [z, y, x] order.
    spacing: float * 3, raw CT spacing in [z, y, x] order.
    new_spacing: float * 3, new spacing used for resample, typically 1x1x1,
        which means standardizing the raw CT with different spacing all into
        1x1x1 mm.
    order: int, order for resample function scipy.ndimage.interpolation.zoom
    return: 3D binary numpy array with the same shape of the image after,
        resampling. The actual resampling spacing is also returned.
    """
    # shape can only be int, so has to be rounded.
    new_shape = np.round(image.shape * spacing / new_spacing)

    # the actual spacing to resample.
    resample_spacing = spacing * image.shape / new_shape

    resize_factor = new_shape / image.shape

    image_new = scipy.ndimage.zoom(image, resize_factor, mode='nearest', order=order)

    return (image_new, resample_spacing)

if do_resample:
    print('Resampling...')
    img_lungRL, resampled_spacing = resample(img_lungRL, spacing, order=3)
    seg_nod_mask = np.zeros(img_lungRL.shape, dtype=np.uint8)
    for i in range(int(nodule_mask.max())):
        # 一个结节,一个结节的resample
        mask = (nodule_mask == (i + 1)).astype(np.uint8)
        mask, _ = resample(mask, spacing, order=3)
        seg_nod_mask[mask > 0.5] = i + 1

其中在resample函数里面,使用到了scipy.ndimage.zoom操作,直接将原始数据,zoom到新的shape

scipy.ndimage.zoom(input, zoom, output=None, order=3, mode='constant', cval=0.0, prefilter=True, *, grid_mode=False)[source]

函数中:

  • input:The input array
  • zoom:The zoom factor along the axes

下面是一段官方案例,展示了zoom前后的变化,可以参考:点击链接直达

from scipy import ndimage, datasets
import matplotlib.pyplot as plt

fig = plt.figure()
ax1 = fig.add_subplot(121)  # left side
ax2 = fig.add_subplot(122)  # right side
ascent = datasets.ascent()
result = ndimage.zoom(ascent, 3.0)
ax1.imshow(ascent, vmin=0, vmax=255)
ax2.imshow(result, vmin=0, vmax=255)
plt.show()

zoom前后的变化,如下所示:
【3D 图像分割】基于 Pytorch 的 VNet 3D 图像分割7(数据预处理),人工智能(AI)医学影像,3d,pytorch,人工智能,LUNA16,LIDC- IDRI

发现这个scipy库还真是好用,后续找时间全面的补充下这块的知识。

2.4、存储到本地

这部分就比较的简单了,主要就是说下数组存储的一些新的:

  • npy文件存储一些简单的数组,比如下文的spacing、坐标等等;
  • nrrd文件存储多维数组,比如下面的imagemask数组图像,大小是240x320x320大小的;
    以前喜欢用nii作为存储文件,现在发现不太好用,nrrd也可以存储数组,还能存储header头。

下面是代码:

lung_box = get_lung_box(binary_mask, img_lungRL.shape)  # 获取肺区分割的外轮廓

z_min, z_max = lung_box[0]
y_min, y_max = lung_box[1]
x_min, x_max = lung_box[2]

# 裁剪操作,去除肺区外的
img_lungRL = img_lungRL[z_min:z_max, y_min:y_max, x_min:x_max]
if do_resample:
    seg_nod_mask = seg_nod_mask[z_min:z_max, y_min:y_max, x_min:x_max]
else:
    seg_nod_mask = nodule_mask[z_min:z_max, y_min:y_max, x_min:x_max]

np.save(os.path.join(save_dir, '%s_origin.npy' % (pid)), origin)  # origin (3,) 记录三维图像origin坐标信息
if do_resample:
    np.save(os.path.join(save_dir, '%s_spacing.npy' % (pid)), resampled_spacing)  # 记录了resample前后x\y\z三个维度的scale系数
np.save(os.path.join(save_dir, '%s_ebox_origin.npy' % (pid)), np.array((z_min, y_min, x_min)))

nrrd.write(os.path.join(save_dir, '%s_clean.nrrd' % (pid)), img_lungRL)  # 去除掉非肺区后的CT图像
nrrd.write(os.path.join(save_dir, '%s_mask.nrrd' % (pid)), seg_nod_mask)  # 去除掉非肺区后的结节MASK图像

2.5、获取结节中心点坐标和半径

这里获取标记结节的中心点坐标和半径,目的还是为了在裁剪patch等操作时候,能够直接从已经获得的结节里面拿取,直接进行crop操作。

这块的步骤和前面get_lung_box差不多,唯一的区别在于保存下来的是中心点,而不是上面的最大、最小边界坐标。

代码如下:

def generateBBoxes_label(params):
    pid, lung_mask_dir, nod_mask_dir, img_dir, save_dir, do_resample = params
    masks, _ = nrrd.read(os.path.join(save_dir, '%s_mask.nrrd' % (pid)))

    bboxes = []
    instance_nums = [num for num in np.unique(masks) if num]
    for i in instance_nums:
        mask = (masks == i).astype(np.uint8)
        zz, yy, xx = np.where(mask)
        d = max(zz.max() - zz.min() + 1, yy.max() - yy.min() + 1, xx.max() - xx.min() + 1)
        bboxes.append(np.array([(zz.max() + zz.min()) / 2., (yy.max() + yy.min()) / 2., (xx.max() + xx.min()) / 2., d]))

    bboxes = np.array(bboxes)
    if not len(bboxes):
        print('%s does not have any nodules!!!' % (pid))

    print('Finished masks to bboxes %s' % (pid))
    np.save(os.path.join(save_dir, '%s_bboxes.npy' % (pid)), bboxes)

三、总结

到这里,本篇内容,结合上一篇的内容,我们对Luna16的数据处理基本上就完成了,也完成了我们最早希望得到的内容:

  1. imagesmask数组,文件名一一对应;
  2. resample操作到1mm
  3. 肺实质外的部分丢弃;

6 和 7 这两个篇章,都是对前面几个章节数据部分的补充,你参考这两篇进行数据处理也行,参考其他的数据处理也行,最终得到的数据形式,只要是一样的就行。文章来源地址https://www.toymoban.com/news/detail-728604.html

到了这里,关于【3D 图像分割】基于 Pytorch 的 VNet 3D 图像分割7(数据预处理)的文章就介绍完了。如果您还想了解更多内容,请在右上角搜索TOY模板网以前的文章或继续浏览下面的相关文章,希望大家以后多多支持TOY模板网!

本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如若转载,请注明出处: 如若内容造成侵权/违法违规/事实不符,请点击违法举报进行投诉反馈,一经查实,立即删除!

领支付宝红包 赞助服务器费用

相关文章

  • 关于图像分割的预处理 transform

    目录 1. 介绍 2. 关于分割中的 resize 问题 3. 分割的 transform 3.1 随机缩放 RandomResize 3.2 随机水平翻转 RandomHorizontalFlip 3.3 随机竖直翻转 RandomVerticalFlip 3.4 随机裁剪 RandomCrop 3.5 ToTensor 3.6 normalization 3.7 Compose 3.8 中心裁剪 3.9 Resize 缩放 4. 预处理结果可视化 图像分割的预处理不像

    2024年02月04日
    浏览(50)
  • 如何基于香橙派AIpro对视频/图像数据进行预处理

    本文分享自华为云社区《如何基于香橙派AIpro对视频/图像数据进行预处理》,作者: 昇腾CANN。 受网络结构和训练方式等因素的影响,绝大多数神经网络模型对输入数据都有格式上的限制。在计算机视觉领域,这个限制大多体现在图像的尺寸、色域、归一化参数等。如果源图

    2024年04月22日
    浏览(52)
  • Pytorch学习笔记(3):图像的预处理(transforms)

      目录  一、torchvision:计算机视觉工具包  二、transforms的运行机制 (1)torchvision.transforms:常用的图像预处理方法 (2)transforms运行原理   三、数据标准化 transforms.Normalize() 四、数据增强  4.1 transforms—数据裁剪 (1)transforms.CentorCrop (2)transforms.RandomCrop (3)RandomResiz

    2023年04月13日
    浏览(46)
  • 【PyTorch】第八节:数据的预处理

    作者 🕵️‍♂️:让机器理解语言か 专栏 🎇:PyTorch 描述 🎨:PyTorch 是一个基于 Torch 的 Python 开源机器学习库。 寄语 💓:🐾没有白走的路,每一步都算数!🐾  ​   torchvision.transforms  是一个包含了常用的图像变化方法的工具包,该工具包主要用于图像预处理、数据增

    2023年04月23日
    浏览(58)
  • 深度学习中基于python的预处理和图像扩增方法

    容易出现的报错: 错误原因通常为保存的路径不正确: 应改为: 即第一个参数应该写到文件的名称,而不能只写到文件夹就停止。 灰度图片和黑白图片有些相似,但并不完全相同。 灰度图片是指每个像素点的颜色由灰度值来表示,通常使用8位无符号整数(0-255)表示。灰

    2024年02月08日
    浏览(45)
  • pytorch入门2--数据预处理、线性代数的矩阵实现、求导

    数据预处理是指将原始数据读取进来使得能用机器学习的方法进行处理。 首先介绍csv文件: CSV 代表逗号分隔值(comma-separated values),CSV 文件就是使用逗号分隔数据的文本文件。 一个 CSV 文件包含一行或多行数据,每一行数据代表一个记录。每个记录包含一个或多个数值,

    2024年02月04日
    浏览(43)
  • 【Computer Vision】图像数据预处理详解

    活动地址:[CSDN21天学习挑战赛](https://marketing.csdn.net/p/bdabfb52c5d56532133df2adc1a728fd) 作者简介 :在校大学生一枚,华为云享专家,阿里云星级博主,腾云先锋(TDP)成员,云曦智划项目总负责人,全国高等学校计算机教学与产业实践资源建设专家委员会(TIPCC)志愿者,以及编程

    2024年02月06日
    浏览(51)
  • 比赛准备笔记 --- TensotFlow、软件调试、数据预处理(图像,csv数据)

    TensorFlow是由Google团队开发的一个开源深度学习框架,完全基于Python语言设计。它的初衷是以最简单的方式实现机器学习和深度学习的概念,结合了计算代数的优化技术,使计算许多数学表达式变得简单。 优势: 强大的计算能力,支持多种硬件和分布式计算 灵活的数据流模型

    2024年02月06日
    浏览(55)
  • Pytorch和Tensoflow对比学习第三周--Day 19-20: 数据加载和预处理

    这两天的学习重点是掌握在PyTorch和TensorFlow中进行数据加载和预处理的方法。正确的数据处理是训练有效模型的关键步骤。 数据加载和预处理: 学习了如何使用PyTorch的DataLoader和Dataset类以及TensorFlow的数据API来加载和预处理数据。 理解了数据标准化、转换和批处理的重要性。

    2024年01月20日
    浏览(58)

觉得文章有用就打赏一下文章作者

支付宝扫一扫打赏

博客赞助

微信扫一扫打赏

请作者喝杯咖啡吧~博客赞助

支付宝扫一扫领取红包,优惠每天领

二维码1

领取红包

二维码2

领红包