【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 = 4do_resample = Trueimg_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 pid in os.listdir(nod_mask_dir):#                         seg-lungs-LUNA16, masks_test/3, seg-lungs-LUNA16, preprocessed_test/3, Trueparams_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, numpySpacingdef 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 minand max, and then normalizeimage: 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_newdef convex_hull_dilate(binary_mask, dilate_factor=1.5, iterations=10):"""Replace each slice with convex hull of it then dilate. Convex hulls usedonly if it does not increase area by dilate_factor. This applies mainly tothe 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 thisspecifical case.dilate_factor: float, factor of increased area after dilationiterations: int, number of iterations for dilationreturn: 3D binary numpy array with the same shape of the image,that only region of interest is True. Each binary mask is ROI of oneside 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_convexstruct = scipy.ndimage.generate_binary_structure(3, 1)binary_mask_dilated = scipy.ndimage.binary_dilation(binary_mask_dilated, structure=struct, iterations=10)return binary_mask_dilateddef 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 interestare 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 notinterested.bone_thred: int, uint8 threahold value for determine parts of image isbone.return: D uint8 numpy array with the same shape of the image afterapplying the lung mask."""binary_mask = binary_mask1 + binary_mask2binary_mask1_dilated = convex_hull_dilate(binary_mask1)binary_mask2_dilated = convex_hull_dilate(binary_mask2)binary_mask_dilated = binary_mask1_dilated + binary_mask2_dilatedbinary_mask_extra = binary_mask_dilated ^ binary_mask# replace image values outside binary_mask_dilated as pad valueimage_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_valuereturn image_newdef cropResample_process(params):#    seg-lungs-LUNA16, masks_test/3, seg-lungs-LUNA16, preprocessed_test/3, Truepid, lung_mask_dir, nod_mask_dir, img_dir, save_dir, do_resample = paramsprint('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 == 4binary_mask_l = lung_mask == 3binary_mask = binary_mask_r + binary_mask_limg_org = HU2uint8(img_org)img_lungRL = apply_lung_mask(img_org, binary_mask_r, binary_mask_l)

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

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

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

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

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

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

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

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

下面是skimage.morphology.convex_hull_image官方给出的实例,如下:点击直达
0作用到我们项目里面,切割后的样子如下:

2

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. 1x1x1image: 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 into1x1x1 mm.order: int, order for resample function scipy.ndimage.interpolation.zoomreturn: 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_shaperesize_factor = new_shape / image.shapeimage_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())):# 一个结节,一个结节的resamplemask = (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 pltfig = 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前后的变化,如下所示:
1

发现这个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 = paramsmasks, _ = 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 这两个篇章,都是对前面几个章节数据部分的补充,你参考这两篇进行数据处理也行,参考其他的数据处理也行,最终得到的数据形式,只要是一样的就行。

本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如若转载,请注明出处:http://www.mzph.cn/news/124770.shtml

如若内容造成侵权/违法违规/事实不符,请联系多彩编程网进行投诉反馈email:809451989@qq.com,一经查实,立即删除!

相关文章

设计模式_状态模式

状态模式 介绍 设计模式定义案例问题堆积在哪里解决办法状态模式一个对象 状态可以发生改变 不同的状态又有不同的行为逻辑游戏角色 加载不同的技能 每个技能有不同的&#xff1a;攻击逻辑 攻击范围 动作等等1 状态很多 2 每个状态有自己的属性和逻辑每种状态单独写一个类 角色…

git本地搭建服务器[Vmware虚拟机访问window的git服务器]

先按照https://zhuanlan.zhihu.com/p/494988089说明下载好Gitblit然后复制到tomcat的webapps目录下,如下: 双击"startup.bat"启动tomcat: 然后访问"http://127.0.0.1:8080/gitblit/"即可看到git的界面: 说明git服务器已经能够成功运行了! Vmware虚拟机…

k8s基本操作命令

目录 1、//查看资源对象简写 2、//查看集群信息 3、//配置kubectl自动补全 4、//node节点查看日志 5、//查看 master 节点状态 6、//查看命令空间 7、//查看default命名空间的所有资源 8、//创建命名空间app 9、//删除命名空间app 10、//在命名空间kube-public 创建…

Java电商平台 - API 接口设计之 token、timestamp、sign 具体架构与实现|电商API接口接入

一&#xff1a;token 简介 Token&#xff1a;访问令牌access token, 用于接口中, 用于标识接口调用者的身份、凭证&#xff0c;减少用户名和密码的传输次数。一般情况下客户端(接口调用方)需要先向服务器端申请一个接口调用的账号&#xff0c;服务器会给出一个appId和一个key, …

anaconda+tensorflow安装完整步骤【亲测可用】

anacondatensorflow安装完整步骤 anaconda安装tensorflow1.安装anaconda2.下载windows版本进行下载并安装3.打开Anaconda Prompt4. 安装tensorflow PyCharm下载与安装1.官网下载pycharm社区版2.PyCharm环境配置3.测试 anaconda安装tensorflow 1.安装anaconda 官网下载anacond…

常用排序算法的理解

1.插入排序 插入排序的思想是将一个记录插入到已经排好序的有序表中&#xff0c;从而形成一个新的、记录数加1的有序表。在其实现过程使用双层循环&#xff0c;外层循环是进行插入的次数&#xff08;也可以理解为比较的轮数&#xff09;&#xff0c;内层循环是当前记录查找插入…

centos部署tomcat

Java Downloads | Oracle 上面是下载网址 Tomcat是由Apache开发的一个Servlet容器&#xff0c;实现了对Servlet和JSP的支持&#xff0c;并提供了作为Web服务器的一些特有功能&#xff0c;如Tomcat管理和控制平台&#xff0c;安全域管理和Tomcat阀 简单来说&#xff1a;Tomcat…

Django之登录注册

最近在准备上线一个网站&#xff08;基于django的编程技术学习与外包服务网站&#xff09;&#xff0c;所以会将自己的在做这个项目的过程中遇到的模块业务以及所涉及到的部分技术记录在CSDN平台里&#xff0c;一是希望可以帮到有需要的同学&#xff0c;二十以供自己后续回顾学…

RabbitMQ生产者的可靠性

目录 MQ使用时会出现的问题 生产者的可靠性 1、生产者重连 2、生产者确认 3、数据持久化 交换机持久化 队列持久化 消息持久化 LazyQueue懒加载 MQ使用时会出现的问题 发送消息时丢失&#xff1a; 生产者发送消息时连接MQ失败生产者发送消息到达MQ后未找到Exchange生…

Jetson Xavier NX FFmpeg支持硬件编解码

最近在用Jetson Xavier NX板子做视频处理&#xff0c;但是CPU进行视频编解码&#xff0c;效率比较地下。 于是便考虑用硬解码来对视频进行处理。 通过jtop查看&#xff0c;发现板子是支持 NVENC硬件编解码的。 1、下载源码 因为需要对ffmpeg进行打补丁修改&#xff0c;因此需…

【c++|opencv】一、基础操作---2.图像信息获取

every blog every motto: You can do more than you think. https://blog.csdn.net/weixin_39190382?typeblog 0. 前言 图像信息获取&#xff0c;roi 1. 图像信息获取 // 获取图像信息#include <iostream> #include <opencv2/opencv.hpp>using namespace cv; …

JVM与Java体系结构

目录 一、Java虚拟机整体架构祥图 二、Java代码执行过程详图 三、汇编语言、机器语言、高级语言关系 四、JVM的架构模型 五、JVM的生命周期 &#xff08;一&#xff09;虚拟机的启动 &#xff08;二&#xff09;虚拟机的执行 &#xff08;三&#xff09;虚拟机的退出 …

使用示例和应用程序全面了解高效数据管理的Golang MySQL数据库

Golang&#xff0c;也被称为Go&#xff0c;已经成为构建强大高性能应用程序的首选语言。在处理MySQL数据库时&#xff0c;Golang提供了一系列强大的库&#xff0c;简化了数据库交互并提高了效率。在本文中&#xff0c;我们将深入探讨一些最流行的Golang MySQL数据库库&#xff…

数据库管理-第113期 Oracle Exadata 04-硬件选择(20231020)

数据库管理-第113期 Oracle Exadata 04-硬件选择&#xff08;2023010290&#xff09; 本周没写文章&#xff0c;主要是因为到上海参加了Oracle CAB/PAB会议&#xff0c;这个放在后面再讲&#xff0c;本期讲一讲Exadata&#xff0c;尤其是存储节点的硬件选择及其对应的一些通用…

mac安装并使用wireshark

mac安装并使用wireshark 1 介绍 我们在日常开发过程中&#xff0c;遇到了棘手的问题时&#xff0c;免不了查看具体网络请求情况&#xff0c;这个时候就需要用到抓包工具。比较著名的抓包工具就属&#xff1a;wireshark、fildder。我这里主要介绍wireshark。 2 安装 以mac安装为…

C#,数值计算——分类与推理Svmpolykernel的计算方法与源程序

1 文本格式 using System; namespace Legalsoft.Truffer { public class Svmpolykernel : Svmgenkernel { public int n { get; set; } public double a { get; set; } public double b { get; set; } public double d { get; set; …

多线程---synchronized特性+原理

文章目录 synchronized特性synchronized原理锁升级/锁膨胀锁消除锁粗化 synchronized特性 互斥 当某个线程执行到某个对象的synchronized中时&#xff0c;其他线程如果也执行到同一个对象的synchronized就会阻塞等待。 进入synchronized修饰的代码块相当于加锁 退出synchronize…

IDEA MyBatisX插件介绍

一、前言 前几年写代码的时候&#xff0c;要一键生成DAO、XML、Entity基础代码会采用第三方工具&#xff0c;比如mybatis-generator-gui等&#xff0c;现在IDEA或Eclipse都有对应的插件&#xff0c;像IDEA中MyBatisX就是一个比较好用的插件。 二、MyBatisX安装配置使用 MyBa…

GO学习之 通道(nil Channel妙用)

GO系列 1、GO学习之Hello World 2、GO学习之入门语法 3、GO学习之切片操作 4、GO学习之 Map 操作 5、GO学习之 结构体 操作 6、GO学习之 通道(Channel) 7、GO学习之 多线程(goroutine) 8、GO学习之 函数(Function) 9、GO学习之 接口(Interface) 10、GO学习之 网络通信(Net/Htt…

Mac电脑Android Studio和VS Code配置Flutter开发环境(图文超详细)

一、安装Android Studio 官网地址&#xff1a; https://developer.android.google.cn/ 历史版本下载地址&#xff1a; https://developer.android.com/studio/archive?hlzh-cn 二、安装Xcode 到App Store下载安装最新版本&#xff0c;如果MacOS更新不到13.0以上就无法安装…