-
Notifications
You must be signed in to change notification settings - Fork 40
/
lung_segmentation.py
98 lines (71 loc) · 3.31 KB
/
lung_segmentation.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
import numpy as np
import skimage.measure
import skimage.segmentation
import skimage.morphology
import skimage.filters
import scipy.ndimage
import utils_plots
def segment_HU_scan(x):
mask = np.asarray(x < -350, dtype='int32')
for iz in xrange(mask.shape[0]):
skimage.segmentation.clear_border(mask[iz], in_place=True)
skimage.morphology.binary_opening(mask[iz], selem=skimage.morphology.disk(5), out=mask[iz])
if np.sum(mask[iz]):
mask[iz] = skimage.morphology.convex_hull_image(mask[iz])
return mask
def segment_HU_scan_frederic(x, threshold=-350):
mask = np.copy(x)
binary_part = mask > threshold
selem1 = skimage.morphology.disk(8)
selem2 = skimage.morphology.disk(2)
selem3 = skimage.morphology.disk(13)
for iz in xrange(mask.shape[0]):
# fill the body part
filled = scipy.ndimage.binary_fill_holes(binary_part[iz]) # fill body
filled_borders_mask = skimage.morphology.binary_erosion(filled, selem1)
mask[iz] *= filled_borders_mask
mask[iz] = skimage.morphology.closing(mask[iz], selem2)
mask[iz] = skimage.morphology.erosion(mask[iz], selem3)
mask[iz] = mask[iz] < threshold
return mask
def segment_HU_scan_elias(x, threshold=-350, pid='test'):
mask = np.copy(x)
binary_part = mask > threshold
selem1 = skimage.morphology.disk(8)
selem2 = skimage.morphology.disk(2)
selem3 = skimage.morphology.disk(9)
for iz in xrange(mask.shape[0]):
# fill the body part
filled = scipy.ndimage.binary_fill_holes(binary_part[iz]) # fill body
filled_borders_mask = skimage.morphology.binary_erosion(filled, selem1)
mask[iz] *= filled_borders_mask
mask[iz] = skimage.morphology.closing(mask[iz], selem2)
mask[iz] = skimage.morphology.erosion(mask[iz], selem3)
mask[iz] = mask[iz] < threshold
#discard disconnected regions
for iz in range(mask.shape[0]/2, mask.shape[0]-1):
overlap = mask[iz] * mask[iz+1]
label_image = skimage.measure.label(mask[iz+1])
for idx, region in enumerate(skimage.measure.regionprops(label_image)):
total_overlap = 0
for coordinates in region.coords:
total_overlap += overlap[coordinates[0], coordinates[1]]
if total_overlap < 2:
print 'region', idx, 'in slice z=', iz-1, 'has no overlap and will be discarded'
for coordinates in region.coords:
mask[iz+1, coordinates[0], coordinates[1]] = 0
for iz in range(mask.shape[0]/2,0,-1 ):
overlap = mask[iz] * mask[iz-1]
label_image = skimage.measure.label(mask[iz-1])
for idx, region in enumerate(skimage.measure.regionprops(label_image)):
total_overlap = 0
for coordinates in region.coords:
total_overlap += overlap[coordinates[0], coordinates[1]]
if total_overlap < 2:
print 'region', idx, 'in slice z=', iz-1, 'has no overlap and will be discarded'
for coordinates in region.coords:
mask[iz-1, coordinates[0], coordinates[1]] = 0
#utils_plots.plot_all_slices(x, mask, pid, './plots/')
return mask
if __name__ == "__main__":
print 'main function to test lung segmentation'