在上一节:【3D 图像分割】基于 Pytorch 的 VNet 3D 图像分割6(数据预处理) 中,我们已经得到了与mhd
图像同seriesUID
名称的mask nrrd
数据文件了,可以说是一一对应了。
并且,mask
的文件,还根据结节被多少人同时标注,区分成了4个文件夹,分别是标注了一、二、三、四次,一共就4
个医生参与标注。
再加上官方已经给整理好的肺实质分割的文件,我们就获得了以下这些数据:
-
ct
图像数据; - 肺实质分割数据;
- 包含结节位置的
mask
数据。
那本篇做什么呢?看下面的导言。
一、导言
上述得到的这些,就满足了我们的需求了,都是一一对应的,无论是后续的数据预处理,还是拿过来用于训练,都非常的方便。
但是呢,对于原始的ct
数据,他在Z
轴上的层厚是不同的,这点可以在dicom
文件里面看到,也可以在mhd
文件的查询到关于层厚的信息。在这点上,不同的序列,差异是非常大的。表现在一个3维数组的结节上面,在这个维度上就是被压扁,和拉长的样子。
在x
和y
方向,其实也是存在spacing
的差异的,但是这种差异没有像z
轴那么夸张的,这里可以选择处理和不处理均可(有些论文进行了处理,有些没有。默认都是512x512
大小,resample
后会变小)。
至此,本篇的目的就很明确了,是要做下面几件事:
- 对原始图像进行肺实质提取,将肺区外的部分进行裁剪,或者改为固定像素值;
- 对图像和结节
mask
进行resample
操作,本篇是zyx
均进行resample
为1mm
。
二、具体实施
怎么做的部分,我们分三部分:
- 肺实质裁剪
-
image
和nodule mask
进行resample
操作 - 获取结节中心点坐标和半径
下面就一一展开
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
:和名称一样,进行肺实质的crop
和resample
操作; -
generateBBoxes_label
:将处理完毕的结节mask
,得到结节中心的坐标和半径。
2.2、肺实质裁剪
这小块的步骤,大概如下:
- 首先,就是数据读取,这部分的详细介绍,可以参考我之前的这篇文章:【医学影像数据处理】nii、npz、npy、dcm、mhd 的数据格式互转,及多目标分割处理汇总
- 其次,就是将
hu
值,转化为0-255
的值,也就是函数HU2uint8()
,对于这部分,可以参考hu
值是如何转为0-255
的可视化部分的介绍:【医学影像数据处理】 Dicom 文件格式处理汇总 - 另外,就是将肺区
mask
作用到图像上,肺实质外采用pad valud
补充 - 最后,将处理好的
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。异或的本质是没有进位的加法。
用真值表来看一下异或的运算逻辑:
看起来非常直观并且简单明了 但是如果我们将布尔值之间的异或换成数字之间的异或会发生什么???让我们来试一下:
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
阶段使用到。
可能会有这样的疑问:为什么不直接image
和lung mask
相乘,得到一个分割肺实质后留下来的image
呢?反而需要采用凸包优化的方式,多此一举呢?
答:在lung mask
里面,肺实质的分割是有误差的。也就是肺实质的分割是沿着肺区边缘的,但是某些结节的位置,恰好在肺区的边界上,且密度很大。那么mask
就会呈现一个内凹的一个状态。如果采用上面的方法,这样结节就被抠除了。采用凸包优化,就可以利用稍微扩展肺实质边缘,达到将更多肺区留下来的效果。
但是,对于肺结核等等大病灶的疾病,采用上述取出肺实质的方法就不行。主要是因为肺结核的病种范围比较大,尽管采用了凸包优化,最终还是会切除很大一块肺区位置,这样肺区就不完整了,有些得不偿失。
下面是skimage.morphology.convex_hull_image
官方给出的实例,如下:点击直达
作用到我们项目里面,切割后的样子如下:
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
前后的变化,如下所示:
发现这个scipy
库还真是好用,后续找时间全面的补充下这块的知识。
2.4、存储到本地
这部分就比较的简单了,主要就是说下数组存储的一些新的:
-
npy
文件存储一些简单的数组,比如下文的spacing
、坐标等等; -
nrrd
文件存储多维数组,比如下面的image
和mask
数组图像,大小是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
的数据处理基本上就完成了,也完成了我们最早希望得到的内容:文章来源:https://www.toymoban.com/news/detail-728604.html
-
images
和mask
数组,文件名一一对应; -
resample
操作到1mm
; - 肺实质外的部分丢弃;
6 和 7 这两个篇章,都是对前面几个章节数据部分的补充,你参考这两篇进行数据处理也行,参考其他的数据处理也行,最终得到的数据形式,只要是一样的就行。文章来源地址https://www.toymoban.com/news/detail-728604.html
到了这里,关于【3D 图像分割】基于 Pytorch 的 VNet 3D 图像分割7(数据预处理)的文章就介绍完了。如果您还想了解更多内容,请在右上角搜索TOY模板网以前的文章或继续浏览下面的相关文章,希望大家以后多多支持TOY模板网!