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')