def seperate_lungs_3d(image):
#Creation of the markers as shown above:
marker_internal, marker_external, marker_watershed = generate_markers_3d(image)
#Creation of the Sobel-Gradient
sobel_filtered_dx = ndimage.sobel(image, axis=2)
sobel_filtered_dy = ndimage.sobel(image, axis=1)
sobel_gradient = np.hypot(sobel_filtered_dx, sobel_filtered_dy)
sobel_gradient *= 255.0 / np.max(sobel_gradient)
#Watershed algorithm
watershed = morphology.watershed(sobel_gradient, marker_watershed)
#Reducing the image created by the Watershed algorithm to its outline
outline = ndimage.morphological_gradient(watershed, size=(1,3,3))
outline = outline.astype(bool)
#Performing Black-Tophat Morphology for reinclusion
#Creation of the disk-kernel and increasing its size a bit
blackhat_struct = [[0, 0, 1, 1, 1, 0, 0],
[0, 1, 1, 1, 1, 1, 0],
[1, 1, 1, 1, 1, 1, 1],
[1, 1, 1, 1, 1, 1, 1],
[1, 1, 1, 1, 1, 1, 1],
[0, 1, 1, 1, 1, 1, 0],
[0, 0, 1, 1, 1, 0, 0]]
blackhat_struct = ndimage.iterate_structure(blackhat_struct, 8)
blackhat_struct = blackhat_struct[np.newaxis,:,:]
#Perform the Black-Hat
outline += ndimage.black_tophat(outline, structure=blackhat_struct) # very long time
#Use the internal marker and the Outline that was just created to generate the lungfilter
lungfilter = np.bitwise_or(marker_internal, outline)
#Close holes in the lungfilter
#fill_holes is not used here, since in some slices the heart would be reincluded by accident
##structure = np.ones((BINARY_CLOSING_SIZE,BINARY_CLOSING_SIZE)) # 5 is not enough, 7 is
structure = morphology.disk(BINARY_CLOSING_SIZE) # better , 5 seems sufficient, we use 7 for safety/just in case
structure = structure[np.newaxis,:,:]
lungfilter = ndimage.morphology.binary_closing(lungfilter, structure=structure, iterations=3) #, iterations=3) # was structure=np.ones((5,5))
### NOTE if no iterattions, i.e. default 1 we get holes within lungs for the disk(5) and perhaps more
#Apply the lungfilter (note the filtered areas being assigned -2000 HU)
segmented = np.where(lungfilter == 1, image, -2000*np.ones(marker_internal.shape))
return segmented, lungfilter, outline, watershed, sobel_gradient, marker_internal, marker_external, marker_watershed
python类disk()的实例源码
def seperate_lungs(image):
#Creation of the markers as shown above:
marker_internal, marker_external, marker_watershed = generate_markers(image)
#Creation of the Sobel-Gradient
sobel_filtered_dx = ndimage.sobel(image, 1)
sobel_filtered_dy = ndimage.sobel(image, 0)
sobel_gradient = np.hypot(sobel_filtered_dx, sobel_filtered_dy)
sobel_gradient *= 255.0 / np.max(sobel_gradient)
#Watershed algorithm
watershed = morphology.watershed(sobel_gradient, marker_watershed)
#Reducing the image created by the Watershed algorithm to its outline
outline = ndimage.morphological_gradient(watershed, size=(3,3))
outline = outline.astype(bool)
#Performing Black-Tophat Morphology for reinclusion
#Creation of the disk-kernel and increasing its size a bit
blackhat_struct = [[0, 0, 1, 1, 1, 0, 0],
[0, 1, 1, 1, 1, 1, 0],
[1, 1, 1, 1, 1, 1, 1],
[1, 1, 1, 1, 1, 1, 1],
[1, 1, 1, 1, 1, 1, 1],
[0, 1, 1, 1, 1, 1, 0],
[0, 0, 1, 1, 1, 0, 0]]
blackhat_struct = ndimage.iterate_structure(blackhat_struct, 8)
#Perform the Black-Hat
outline += ndimage.black_tophat(outline, structure=blackhat_struct)
#Use the internal marker and the Outline that was just created to generate the lungfilter
lungfilter = np.bitwise_or(marker_internal, outline)
#Close holes in the lungfilter
#fill_holes is not used here, since in some slices the heart would be reincluded by accident
##structure = np.ones((BINARY_CLOSING_SIZE,BINARY_CLOSING_SIZE)) # 5 is not enough, 7 is
structure = morphology.disk(BINARY_CLOSING_SIZE) # better , 5 seems sufficient, we use 7 for safety/just in case
lungfilter = ndimage.morphology.binary_closing(lungfilter, structure=structure, iterations=3) #, iterations=3) # was structure=np.ones((5,5))
### NOTE if no iterattions, i.e. default 1 we get holes within lungs for the disk(5) and perhaps more
#Apply the lungfilter (note the filtered areas being assigned -2000 HU)
segmented = np.where(lungfilter == 1, image, -2000*np.ones((512, 512))) ### was -2000
return segmented, lungfilter, outline, watershed, sobel_gradient, marker_internal, marker_external, marker_watershed
def seperate_lungs_3d(image):
#Creation of the markers as shown above:
marker_internal, marker_external, marker_watershed = generate_markers_3d(image)
#Creation of the Sobel-Gradient
sobel_filtered_dx = ndimage.sobel(image, axis=2)
sobel_filtered_dy = ndimage.sobel(image, axis=1)
sobel_gradient = np.hypot(sobel_filtered_dx, sobel_filtered_dy)
sobel_gradient *= 255.0 / np.max(sobel_gradient)
#Watershed algorithm
watershed = morphology.watershed(sobel_gradient, marker_watershed)
#Reducing the image created by the Watershed algorithm to its outline
outline = ndimage.morphological_gradient(watershed, size=(1,3,3))
outline = outline.astype(bool)
#Performing Black-Tophat Morphology for reinclusion
#Creation of the disk-kernel and increasing its size a bit
blackhat_struct = [[0, 0, 1, 1, 1, 0, 0],
[0, 1, 1, 1, 1, 1, 0],
[1, 1, 1, 1, 1, 1, 1],
[1, 1, 1, 1, 1, 1, 1],
[1, 1, 1, 1, 1, 1, 1],
[0, 1, 1, 1, 1, 1, 0],
[0, 0, 1, 1, 1, 0, 0]]
blackhat_struct = ndimage.iterate_structure(blackhat_struct, 8)
blackhat_struct = blackhat_struct[np.newaxis,:,:]
#Perform the Black-Hat
outline += ndimage.black_tophat(outline, structure=blackhat_struct) # very long time
#Use the internal marker and the Outline that was just created to generate the lungfilter
lungfilter = np.bitwise_or(marker_internal, outline)
#Close holes in the lungfilter
#fill_holes is not used here, since in some slices the heart would be reincluded by accident
##structure = np.ones((BINARY_CLOSING_SIZE,BINARY_CLOSING_SIZE)) # 5 is not enough, 7 is
structure = morphology.disk(BINARY_CLOSING_SIZE) # better , 5 seems sufficient, we use 7 for safety/just in case
structure = structure[np.newaxis,:,:]
lungfilter = ndimage.morphology.binary_closing(lungfilter, structure=structure, iterations=3) #, iterations=3) # was structure=np.ones((5,5))
### NOTE if no iterattions, i.e. default 1 we get holes within lungs for the disk(5) and perhaps more
#Apply the lungfilter (note the filtered areas being assigned -2000 HU)
segmented = np.where(lungfilter == 1, image, -2000*np.ones(marker_internal.shape))
return segmented, lungfilter, outline, watershed, sobel_gradient, marker_internal, marker_external, marker_watershed
def segment_lung_mask(image, speedup=4):
def largest_label_volume(im, bg=-1):
vals, counts = np.unique(im, return_counts=True)
counts = counts[vals != bg]
vals = vals[vals != bg]
if len(counts) > 0:
return vals[np.argmax(counts)]
else:
return None
if speedup>1:
smallImage = transform.downscale_local_mean(image,(1,speedup,speedup));
else:
smallImage = image;
# not actually binary, but 1 and 2.
# 0 is treated as background, which we do not want
binary_image = np.array((smallImage < -320) & (smallImage>-1400), dtype=np.int8)
#return binary_image;
for i, axial_slice in enumerate(binary_image):
axial_slice = 1-axial_slice
labeling = measure.label(axial_slice)
l_max = largest_label_volume(labeling, bg=0)
if l_max is not None: #This slice contains some lung
binary_image[i][(labeling!=l_max)] = 1
# Remove other air pockets insided body
labels = measure.label(binary_image, background=0)
m = labels.shape[0]//2;
check_layers = labels[m-12:m+20:4,:,:];
l_max = largest_label_volume(check_layers, bg=0)
while l_max is not None: # There are air pockets
idx = np.where(check_layers==l_max);
ii = np.vstack(idx[1:]).flatten();
if np.max(ii)>labels.shape[1]-24/speedup or np.min(ii)<24/speedup:
binary_image[labels==l_max] = 0;
labels = measure.label(binary_image, background=0)
m = labels.shape[0]//2;
check_layers = labels[m-12:m+20:4,:,:];
l_max = largest_label_volume(check_layers, bg=0)
else:
binary_image[labels != l_max] = 0
break
if speedup<=1:
return binary_image
else:
res = np.zeros(image.shape,dtype=np.uint8);
for i,x in enumerate(binary_image):
orig = np.copy(x);
x = binary_dilation(x,disk(5))
x = binary_erosion(x,disk(5))
x = np.logical_or(x,orig)
y = transform.resize(x*1.0,image.shape[1:3]);
res[i][y>0.5]=1
return res;
def unet_candidates():
cands = glob.glob("../data/predictions_epoch9_23_all/*.png")
#df = pd.DataFrame(columns=['seriesuid','coordX','coordY','coordZ','class'])
data = []
imname = ""
origin = []
spacing = []
nrimages = 0
for name in tqdm(cands):
#image = imread(name)
image_t = imread(name)
image_t = image_t.transpose()
#Thresholding
image_t[image_t<THRESHOLD] = 0
image_t[image_t>0] = 1
#erosion
selem = morphology.disk(1)
image_eroded = image_t
image_eroded = morphology.binary_erosion(image_t,selem=selem)
label_im, nb_labels = ndimage.label(image_eroded)
imname3 = os.path.split(name)[1].replace('.png','')
splitted = imname3.split("slice")
slice = splitted[1]
imname2 = splitted[0][:-1]
centers = []
for i in xrange(1,nb_labels+1):
blob_i = np.where(label_im==i,1,0)
mass = center_of_mass(blob_i)
centers.append([mass[1],mass[0]])
if imname2 != imname:
if os.path.isfile("../data/1_1_1mm_512_x_512_annotation_masks/spacings/{0}.pickle".format(imname2)):
with open("../data/1_1_1mm_512_x_512_annotation_masks/spacings/{0}.pickle".format(imname2), 'rb') as handle:
dic = pickle.load(handle)
origin = dic["origin"]
spacing = dic["spacing"]
imname = imname2
nrimages +=1
for center in centers:
coords = voxel_2_world([int(slice),center[1]+(512-324)*0.5,center[0]+(512-324)*0.5],origin,spacing)
data.append([imname2,coords[2],coords[1],coords[0],'?'])
#if nrimages == 5:
# break
df = pd.DataFrame(data,columns=CANDIDATES_COLUMNS)
save_candidates("../data/candidates_unet_final_23.csv",df)
def estimate_background(self, method, sigma_clip=True, sigma_level=3.0):
"""
This function ...
:param method:
:param sigma_clip:
:param sigma_level:
:return:
"""
# Make a distinction between the "PTS" way of estimating the background and all other methods.
# For the PTS way, in case sigma clipping is disabled, this means it is disabled only for the 'polynomial fitting step'
# of the background estimation, so provide two distinct masks to the interpolated() method: the clipped mask for
# the 'background noise' estimation and the non-clipped mask for the 'polynomial fitting step'.
if method == "pts":
if sigma_clip:
try:
mask = statistics.sigma_clip_mask(self.cutout, sigma_level=sigma_level, mask=self.mask)
except TypeError:
#plotting.plot_box(self.cutout)
#plotting.plot_mask(self.mask)
#print("xsize", self.cutout.xsize, self.cutout.ysize)
radius = int(round(0.25 * self.cutout.xsize))
#print("radius", 0.25*self.cutout.xsize, radius)
disk = morphology.disk(radius, dtype=bool)
mask = Mask.empty_like(self.cutout)
x_min = int(round(0.5 * (self.cutout.xsize - disk.shape[1])))
y_min = int(round(0.5 * (self.cutout.ysize - disk.shape[0])))
#plotting.plot_mask(mask)
mask[y_min:y_min+disk.shape[0], x_min:x_min+disk.shape[1]] = disk
#plotting.plot_mask(mask)
no_clip_mask = None
else:
mask = statistics.sigma_clip_mask(self.cutout, sigma_level=sigma_level, mask=self.mask)
no_clip_mask = self.mask
else:
# Perform sigma-clipping on the background if requested
if sigma_clip: mask = statistics.sigma_clip_mask(self.cutout, sigma_level=sigma_level, mask=self.mask)
else: mask = self.mask
no_clip_mask = None
if self.contamination is not None:
mask = mask + self.contamination
if no_clip_mask is not None: no_clip_mask = no_clip_mask + self.contamination
# Perform the interpolation
self.background = self.cutout.interpolated(mask, method, no_clip_mask=no_clip_mask, plot=self.special)
if self.special: self.plot(title="background estimated")
# -----------------------------------------------------------------
def estimate_background(self, method, sigma_clip=True, sigma_level=3.0):
"""
This function ...
:param method:
:param sigma_clip:
:param sigma_level:
:return:
"""
# Make a distinction between the "PTS" way of estimating the background and all other methods.
# For the PTS way, in case sigma clipping is disabled, this means it is disabled only for the 'polynomial fitting step'
# of the background estimation, so provide two distinct masks to the interpolated() method: the clipped mask for
# the 'background noise' estimation and the non-clipped mask for the 'polynomial fitting step'.
if method == "pts":
if sigma_clip:
try:
mask = statistics.sigma_clip_mask(self.cutout, sigma_level=sigma_level, mask=self.mask)
except TypeError:
#plotting.plot_box(self.cutout)
#plotting.plot_mask(self.mask)
#print("xsize", self.cutout.xsize, self.cutout.ysize)
radius = int(round(0.25 * self.cutout.xsize))
#print("radius", 0.25*self.cutout.xsize, radius)
disk = morphology.disk(radius, dtype=bool)
mask = Mask.empty_like(self.cutout)
x_min = int(round(0.5 * (self.cutout.xsize - disk.shape[1])))
y_min = int(round(0.5 * (self.cutout.ysize - disk.shape[0])))
#plotting.plot_mask(mask)
mask[y_min:y_min+disk.shape[0], x_min:x_min+disk.shape[1]] = disk
#plotting.plot_mask(mask)
no_clip_mask = None
else:
mask = statistics.sigma_clip_mask(self.cutout, sigma_level=sigma_level, mask=self.mask)
no_clip_mask = self.mask
else:
# Perform sigma-clipping on the background if requested
if sigma_clip: mask = statistics.sigma_clip_mask(self.cutout, sigma_level=sigma_level, mask=self.mask)
else: mask = self.mask
no_clip_mask = None
if self.contamination is not None:
mask = mask + self.contamination
if no_clip_mask is not None: no_clip_mask = no_clip_mask + self.contamination
# Perform the interpolation
self.background = self.cutout.interpolated(mask, method, no_clip_mask=no_clip_mask, plot=self.special)
if self.special: self.plot(title="background estimated")
# -----------------------------------------------------------------
def get_masks(im):
'''
Step 1: Convert into a binary image.
'''
print('step1')
binary = im < 604
# plt.imshow(binary,cmap=plt.cm.gray)
# plt.show()
'''
Step 2: Remove the blobs connected to the border of the image.
'''
print('step2')
cleared = clear_border(binary)
# plt.imshow(cleared,cmap=plt.cm.gray)
# plt.show()
'''
Step 3: Label the image.
'''
print('step3')
label_image = label(cleared)
# plt.imshow(label_image,cmap=plt.cm.gray)
# plt.show()
'''
Step 4: Keep the labels with 2 largest areas.
'''
print('step4')
areas = [r.area for r in regionprops(label_image)]
areas.sort()
if len(areas) > 2:
for region in regionprops(label_image):
if region.area < 10 and region.area > 3:
print(region.centroid,region.area)
# print(region.area)
centroid = region.centroid
plot_im(im,centroid)
# label_image[int(centroid[0]),int(centroid[1])] = 1000
# for coordinates in region.coords:
# label_image[coordinates[0], coordinates[1]] = 0
# binary = label_image > 999
# plt.imshow(binary,cmap=plt.cm.gray)
# plt.show()
'''
Step 5: Erosion operation with a disk of radius 2. This operation is
seperate the lung nodules attached to the blood vessels.
'''
# print('step5')
# selem = disk(2)
# binary = binary_erosion(binary, selem)
# plt.imshow(binary,cmap=plt.cm.gray)
# plt.show()
def img_tesseract_detect(c_rect, im):
# ????minAreaRect??????-90~0??????????????????
# ???????????????????????????????????????
pts = c_rect.reshape(4, 2)
rect = np.zeros((4, 2), dtype = "float32")
# the top-left point has the smallest sum whereas the
# bottom-right has the largest sum
s = pts.sum(axis = 1)
rect[0] = pts[np.argmin(s)]
rect[3] = pts[np.argmax(s)]
# compute the difference between the points -- the top-right
# will have the minumum difference and the bottom-left will
# have the maximum difference
diff = np.diff(pts, axis = 1)
rect[2] = pts[np.argmin(diff)]
rect[1] = pts[np.argmax(diff)]
dst = np.float32([[0,0],[0,100],[200,0],[200,100]])
M = cv2.getPerspectiveTransform(rect, dst)
warp = cv2.warpPerspective(im, M, (200, 100))
img_show_hook("??????", warp)
warp = np.array(warp, dtype=np.uint8)
radius = 10
selem = disk(radius)
#????????OTSU????
local_otsu = rank.otsu(warp, selem)
l_otsu = np.uint8(warp >= local_otsu)
l_otsu *= 255
kernel = cv2.getStructuringElement(cv2.MORPH_RECT,(4, 4))
l_otsu = cv2.morphologyEx(l_otsu, cv2.MORPH_CLOSE, kernel)
img_show_hook("?????OTSU??", l_otsu)
print("?????")
print(pytesseract.image_to_string(Image.fromarray(l_otsu)))
cv2.waitKey(0)
return
def img_tesseract_detect(c_rect, im):
# ????minAreaRect??????-90~0??????????????????
# ???????????????????????????????????????
pts = c_rect.reshape(4, 2)
rect = np.zeros((4, 2), dtype = "float32")
# the top-left point has the smallest sum whereas the
# bottom-right has the largest sum
s = pts.sum(axis = 1)
rect[0] = pts[np.argmin(s)]
rect[3] = pts[np.argmax(s)]
# compute the difference between the points -- the top-right
# will have the minumum difference and the bottom-left will
# have the maximum difference
diff = np.diff(pts, axis = 1)
rect[2] = pts[np.argmin(diff)]
rect[1] = pts[np.argmax(diff)]
width = rect[3][0] - rect[0][0]
height = rect[3][1] - rect[0][1]
width = (int)((50.0 / height) * width)
height = 50
dst = np.float32([[0,0],[0,height],[width,0],[width,height]])
M = cv2.getPerspectiveTransform(rect, dst)
warp = cv2.warpPerspective(im, M, (width, height))
img_show_hook("??????", warp)
warp = np.array(warp, dtype=np.uint8)
radius = 13
selem = disk(radius)
#????????OTSU????
local_otsu = rank.otsu(warp, selem)
l_otsu = np.uint8(warp >= local_otsu)
l_otsu *= 255
kernel = cv2.getStructuringElement(cv2.MORPH_RECT,(2,2))
l_otsu = cv2.morphologyEx(l_otsu, cv2.MORPH_CLOSE, kernel)
img_show_hook("?????OTSU??", l_otsu)
print("?????")
print(pytesseract.image_to_string(Image.fromarray(l_otsu), lang="chi-sim"))
cv2.waitKey(0)
return
def unet_candidates():
cands = glob.glob('%s/*.png' % sys.argv[1]) #"/razberry/workspace/dsb_nodule_detection.109fd54/*.png
#df = pd.DataFrame(columns=['seriesuid','coordX','coordY','coordZ','class'])
data = []
imname = ""
origin = []
spacing = []
nrimages = 0
for name in tqdm(cands):
#image = imread(name)
image_t = imread(name)
image_t = image_t.transpose()
#Thresholding
image_t[image_t<THRESHOLD] = 0
image_t[image_t>0] = 1
#erosion
selem = morphology.disk(1)
image_eroded = image_t
image_eroded = morphology.binary_erosion(image_t,selem=selem)
label_im, nb_labels = ndimage.label(image_eroded)
imname3 = os.path.split(name)[1].replace('.png','')
splitted = imname3.split("_")
slice = splitted[1]
imname2 = splitted[0][:-1]
centers = []
for i in xrange(1,nb_labels+1):
blob_i = np.where(label_im==i,1,0)
mass = center_of_mass(blob_i)
centers.append([mass[1],mass[0]])
if imname2 != imname:
if os.path.isfile("../data/1_1_1mm_512_x_512_annotation_masks/spacings/{0}.pickle".format(imname2)):
with open("../data/1_1_1mm_512_x_512_annotation_masks/spacings/{0}.pickle".format(imname2), 'rb') as handle:
dic = pickle.load(handle)
origin = dic["origin"]
spacing = dic["spacing"]
imname = imname2
nrimages +=1
for center in centers:
# coords = voxel_2_world([int(slice),center[1]+(512-324)*0.5,center[0]+(512-324)*0.5],origin,spacing)
coords = [int(slice),center[1],center[0]]
data.append([imname2,coords[2],coords[1],coords[0],'?'])
#if nrimages == 5:
# break
df = pd.DataFrame(data,columns=CANDIDATES_COLUMNS)
save_candidates('%s/candidates.csv' % work_dir, df)
def get_segmented_lungs(image):
#Creation of the markers as shown above:
marker_internal, marker_external, marker_watershed = generate_markers(image)
#Creation of the Sobel-Gradient
sobel_filtered_dx = ndimage.sobel(image, 1)
sobel_filtered_dy = ndimage.sobel(image, 0)
sobel_gradient = np.hypot(sobel_filtered_dx, sobel_filtered_dy)
sobel_gradient *= 255.0 / np.max(sobel_gradient)
#Watershed algorithm
watershed = morphology.watershed(sobel_gradient, marker_watershed)
#Reducing the image created by the Watershed algorithm to its outline
outline = ndimage.morphological_gradient(watershed, size=(3,3))
outline = outline.astype(bool)
#Performing Black-Tophat Morphology for reinclusion
#Creation of the disk-kernel and increasing its size a bit
blackhat_struct = [[0, 0, 1, 1, 1, 0, 0],
[0, 1, 1, 1, 1, 1, 0],
[1, 1, 1, 1, 1, 1, 1],
[1, 1, 1, 1, 1, 1, 1],
[1, 1, 1, 1, 1, 1, 1],
[0, 1, 1, 1, 1, 1, 0],
[0, 0, 1, 1, 1, 0, 0]]
#blackhat_struct = ndimage.iterate_structure(blackhat_struct, 8)
blackhat_struct = ndimage.iterate_structure(blackhat_struct, 14) # <- retains more of the area, 12 works well. Changed to 14, 12 still excluded some parts.
#Perform the Black-Hat
outline += ndimage.black_tophat(outline, structure=blackhat_struct)
#Use the internal marker and the Outline that was just created to generate the lungfilter
lungfilter = np.bitwise_or(marker_internal, outline)
#Close holes in the lungfilter
#fill_holes is not used here, since in some slices the heart would be reincluded by accident
lungfilter = ndimage.morphology.binary_closing(lungfilter, structure=np.ones((5,5)), iterations=3)
#Apply the lungfilter (note the filtered areas being assigned threshold_min HU)
segmented = np.where(lungfilter == 1, image, threshold_min*np.ones(image.shape))
#return segmented, lungfilter, outline, watershed, sobel_gradient, marker_internal, marker_external, marker_watershed
return segmented
def get_segmented_lungs(image):
#Creation of the markers as shown above:
marker_internal, marker_external, marker_watershed = generate_markers(image)
#Creation of the Sobel-Gradient
sobel_filtered_dx = ndimage.sobel(image, 1)
sobel_filtered_dy = ndimage.sobel(image, 0)
sobel_gradient = np.hypot(sobel_filtered_dx, sobel_filtered_dy)
sobel_gradient *= 255.0 / np.max(sobel_gradient)
#Watershed algorithm
watershed = morphology.watershed(sobel_gradient, marker_watershed)
#Reducing the image created by the Watershed algorithm to its outline
outline = ndimage.morphological_gradient(watershed, size=(3,3))
outline = outline.astype(bool)
#Performing Black-Tophat Morphology for reinclusion
#Creation of the disk-kernel and increasing its size a bit
blackhat_struct = [[0, 0, 1, 1, 1, 0, 0],
[0, 1, 1, 1, 1, 1, 0],
[1, 1, 1, 1, 1, 1, 1],
[1, 1, 1, 1, 1, 1, 1],
[1, 1, 1, 1, 1, 1, 1],
[0, 1, 1, 1, 1, 1, 0],
[0, 0, 1, 1, 1, 0, 0]]
#blackhat_struct = ndimage.iterate_structure(blackhat_struct, 8)
blackhat_struct = ndimage.iterate_structure(blackhat_struct, 14) # <- retains more of the area, 12 works well. Changed to 14, 12 still excluded some parts.
#Perform the Black-Hat
outline += ndimage.black_tophat(outline, structure=blackhat_struct)
#Use the internal marker and the Outline that was just created to generate the lungfilter
lungfilter = np.bitwise_or(marker_internal, outline)
#Close holes in the lungfilter
#fill_holes is not used here, since in some slices the heart would be reincluded by accident
lungfilter = ndimage.morphology.binary_closing(lungfilter, structure=np.ones((5,5)), iterations=3)
#Apply the lungfilter (note the filtered areas being assigned threshold_min HU)
segmented = np.where(lungfilter == 1, image, threshold_min*np.ones(image.shape))
#return segmented, lungfilter, outline, watershed, sobel_gradient, marker_internal, marker_external, marker_watershed
return segmented
def get_segmented_lungs(image):
#Creation of the markers as shown above:
marker_internal, marker_external, marker_watershed = generate_markers(image)
#Creation of the Sobel-Gradient
sobel_filtered_dx = ndimage.sobel(image, 1)
sobel_filtered_dy = ndimage.sobel(image, 0)
sobel_gradient = np.hypot(sobel_filtered_dx, sobel_filtered_dy)
sobel_gradient *= 255.0 / np.max(sobel_gradient)
#Watershed algorithm
watershed = morphology.watershed(sobel_gradient, marker_watershed)
#Reducing the image created by the Watershed algorithm to its outline
outline = ndimage.morphological_gradient(watershed, size=(3,3))
outline = outline.astype(bool)
#Performing Black-Tophat Morphology for reinclusion
#Creation of the disk-kernel and increasing its size a bit
blackhat_struct = [[0, 0, 1, 1, 1, 0, 0],
[0, 1, 1, 1, 1, 1, 0],
[1, 1, 1, 1, 1, 1, 1],
[1, 1, 1, 1, 1, 1, 1],
[1, 1, 1, 1, 1, 1, 1],
[0, 1, 1, 1, 1, 1, 0],
[0, 0, 1, 1, 1, 0, 0]]
#blackhat_struct = ndimage.iterate_structure(blackhat_struct, 8)
blackhat_struct = ndimage.iterate_structure(blackhat_struct, 14) # <- retains more of the area, 12 works well. Changed to 14, 12 still excluded some parts.
#Perform the Black-Hat
outline += ndimage.black_tophat(outline, structure=blackhat_struct)
#Use the internal marker and the Outline that was just created to generate the lungfilter
lungfilter = np.bitwise_or(marker_internal, outline)
#Close holes in the lungfilter
#fill_holes is not used here, since in some slices the heart would be reincluded by accident
lungfilter = ndimage.morphology.binary_closing(lungfilter, structure=np.ones((5,5)), iterations=3)
#Apply the lungfilter (note the filtered areas being assigned threshold_min HU)
segmented = np.where(lungfilter == 1, image, threshold_min*np.ones(image.shape))
#return segmented, lungfilter, outline, watershed, sobel_gradient, marker_internal, marker_external, marker_watershed
return segmented
