SimpleITK 三维图像分析

1、去除3D 小连通域
　　在一些计算机视觉任务中，需要对模型的输出做一些后处理以优化视觉效果，连通域就是一种常见的后处理方式。尤其对于分割任务，有时的输出mask会存在一些假阳（小的无用轮廓），通过3D连通域找出面积较小的独立轮廓并去除可以有效地提升视觉效果。

　　二维图像连通域一般包括 4连通、8连通。对于三维数据一般包括6连通、18连通和26联通。
下面的代码只保留最大3D连通域。

 1 # -*- coding : UTF-8 -*-
 2 # @file   : prob2label.py
 3 # @Time   : 2021-10-19 9:35
 4 # @Author : wmz
 5 
 6 import os
 7 import SimpleITK as sitk
 8 from glob import glob
 9 import numpy as np
10 
11 
12 def getFiles(path, suffix):
13     return [os.path.join(root, file) for root, dirs, files in os.walk(path) for file in files if file.endswith(suffix)]
14 
15 
16 def connected_domain_2(image, mask=True):
17     cca = sitk.ConnectedComponentImageFilter()
18     cca.SetFullyConnected(True)
19     _input = sitk.GetImageFromArray(image.astype(np.uint8))
20     output_ex = cca.Execute(_input)
21     stats = sitk.LabelShapeStatisticsImageFilter()
22     stats.Execute(output_ex)
23     num_label = cca.GetObjectCount()
24     num_list = [i for i in range(1, num_label+1)]
25     area_list = []
26     for l in range(1, num_label +1):
27         area_list.append(stats.GetNumberOfPixels(l))
28     num_list_sorted = sorted(num_list, key=lambda x: area_list[x-1])[::-1]
29     largest_area = area_list[num_list_sorted[0] - 1]
30     final_label_list = [num_list_sorted[0]]
31 
32     # for idx, i in enumerate(num_list_sorted[1:]):  # 大于第一个的十分之一的都保留，注释掉之后只保留最大连通域
33     #     if area_list[i-1] >= (largest_area//10):
34     #         final_label_list.append(i)
35     #     else:
36     #         break
37     output = sitk.GetArrayFromImage(output_ex)
38 
39     for one_label in num_list:
40         if one_label in final_label_list:
41             continue
42         x, y, z, w, h, d = stats.GetBoundingBox(one_label)
43         one_mask = (output[z: z + d, y: y + h, x: x + w] != one_label)
44         output[z: z + d, y: y + h, x: x + w] *= one_mask
45 
46     if mask:
47         output = (output > 0).astype(np.uint8)
48     else:
49         output = ((output > 0)*255.).astype(np.uint8)
50     return output
51 
52 
53 def save_prob2label(prob_dir, save_labeldir):
54     # all_prob_seg = glob(os.path.join(prob_dir, "*.nrrd"))
55     all_prob_seg = getFiles(prob_dir, ".nrrd")
56     for index, file in enumerate(all_prob_seg):
57         print("processing", index + 1, '/', len(all_prob_seg), file)
58         label_file = file.replace(prob_dir, save_labeldir).replace(".nrrd", ".nii.gz")
59         prob_img = sitk.ReadImage(file)
60         prob_arr = sitk.GetArrayFromImage(prob_img)
61         label_arr = (prob_arr > Dice_value) * 1
62         label_arr = connected_domain_2(label_arr)
63         label_img = sitk.GetImageFromArray(label_arr)
64         label_img.SetOrigin(prob_img.GetOrigin())
65         label_img.SetDirection(prob_img.GetDirection())
66         dst_dir = label_file.rsplit('\\', 1)[0]
67         if not os.path.exists(dst_dir):
68             os.makedirs(dst_dir)
69         sitk.WriteImage(label_img, label_file)
70 
71 
72 if __name__ == '__main__':
73 
74     prob_nrrd_dir = r'C:\Users\wmz\Desktop\input'
75     save_label_dir = r'C:\Users\wmz\Desktop\test'
76     Dice_value = 0.5
77     save_prob2label(prob_nrrd_dir, save_label_dir)

2、【医学图像处理】之腹部骨骼提取（SimpleITK）
1.内容
步骤：
1.读取Dicom序列
2.设置固定阈值为100，把骨骼和心脏及主动脉都分割出来
3.形态学开运算+最大连通域提取,粗略的心脏和主动脉图像
4.将step1的结果与step2的结果相减,得到骨骼部分
5.最大连通域提取，去除小连接
6.将得到的图像与原始图像进行逻辑与操作

 1 import SimpleITK as sitk
 2 
 3 # 最大连通域提取
 4 def GetLargestConnectedCompont(binarysitk_image):
 5     cc = sitk.ConnectedComponent(binarysitk_image)
 6     stats = sitk.LabelIntensityStatisticsImageFilter()
 7     stats.SetGlobalDefaultNumberOfThreads(8)
 8     stats.Execute(cc, binarysitk_image)
 9     maxlabel = 0
10     maxsize = 0
11     for l in stats.GetLabels():
12         size = stats.GetPhysicalSize(l)
13         if maxsize < size:
14             maxlabel = l
15             maxsize = size
16     labelmaskimage = sitk.GetArrayFromImage(cc)
17     outmask = labelmaskimage.copy()
18     outmask[labelmaskimage == maxlabel] = 255
19     outmask[labelmaskimage != maxlabel] = 0
20     outmask_sitk = sitk.GetImageFromArray(outmask)
21     outmask_sitk.SetDirection(binarysitk_image.GetDirection())
22     outmask_sitk.SetSpacing(binarysitk_image.GetSpacing())
23     outmask_sitk.SetOrigin(binarysitk_image.GetOrigin())
24     return outmask_sitk
25 
26 # 逻辑与操作
27 def GetMaskImage(sitk_src, sitk_mask, replacevalue=0):
28     array_src = sitk.GetArrayFromImage(sitk_src)
29     array_mask = sitk.GetArrayFromImage(sitk_mask)
30     array_out = array_src.copy()
31     array_out[array_mask == 0] = replacevalue
32     outmask_sitk = sitk.GetImageFromArray(array_out)
33     outmask_sitk.SetDirection(sitk_src.GetDirection())
34     outmask_sitk.SetSpacing(sitk_src.GetSpacing())
35     outmask_sitk.SetOrigin(sitk_src.GetOrigin())
36     return outmask_sitk
37 
38 
39 # 读取Dicom序列
40 pathDicom = 'D:/PyCharm 2019.3.3/data/LIDC_nodul'
41 reader = sitk.ImageSeriesReader()
42 filenamesDICOM = reader.GetGDCMSeriesFileNames(pathDicom)
43 reader.SetFileNames(filenamesDICOM)
44 sitk_src = reader.Execute()
45 
46 # step1.设置固定阈值为100，把骨骼和心脏及主动脉都分割出来
47 sitk_seg = sitk.BinaryThreshold(sitk_src, lowerThreshold=100, upperThreshold=3000, insideValue=255, outsideValue=0)
48 sitk.WriteImage(sitk_seg, 'step1.mha')
49 
50 # step2.形态学开运算+最大连通域提取,粗略的心脏和主动脉图像
51 sitk_open = sitk.BinaryMorphologicalOpening(sitk_seg != 0, 2)
52 sitk_open = GetLargestConnectedCompont(sitk_open)
53 sitk.WriteImage(sitk_open, 'step2.mha')
54 
55 # step3.再将step1的结果与step2的结果相减,得到骨骼部分
56 array_open = sitk.GetArrayFromImage(sitk_open)
57 array_seg = sitk.GetArrayFromImage(sitk_seg)
58 array_mask = array_seg - array_open
59 sitk_mask = sitk.GetImageFromArray(array_mask)
60 sitk_mask.SetDirection(sitk_seg.GetDirection())
61 sitk_mask.SetSpacing(sitk_seg.GetSpacing())
62 sitk_mask.SetOrigin(sitk_seg.GetOrigin())
63 sitk.WriteImage(sitk_mask, 'step3.mha')
64 
65 # step4.最大连通域提取，去除小连接
66 skeleton_mask = GetLargestConnectedCompont(sitk_mask)
67 sitk.WriteImage(skeleton_mask, 'step4.mha')
68 
69 # step5.将得到的图像与原始图像进行逻辑与操作
70 sitk_skeleton = GetMaskImage(sitk_src, skeleton_mask, replacevalue=-1500)
71 sitk.WriteImage(sitk_skeleton, 'step5.mha')