Как улучшить сегментацию легких без глубокого обучения

Question

Моя цель состоит в том, чтобы получить логическое ИЛИ ROI легких по оси срезов на воксель с помощью трехмерного КТ-сканирования легких, используя быстрые стандартные подходы .

Я пробовал использовать разные подходы, такие как пороговое значение и морфология и водораздел , но результаты не идеальны, как вы можете видеть на изображении внизу.

Я взял логическое ИЛИ всех масок на воксел для сравнения с результирующим результатом моих манипуляций и взял стандартное отклонение вокселей вдоль оси срезов, чтобы извлечь ROI в зависимости от него.

Мой вопрос:

Как сделать так, чтобы Результат выглядел очень близко к маскам после логического ИЛИ ??

Используемый набор данных (изображения и маски): Здесь .

в следующем коде и пытались объединить пороговое значение и морфологию, но проблема в том, что:

# import necessary libraries
import os
import numpy as np
import matplotlib.pyplot as plt 
# plt.switch_backend('Qt5Agg')
import glob
import nibabel as nib
from skimage.io import imsave
from skimage.filters import threshold_otsu
from skimage.transform import resize
import skimage.morphology as morph
from skimage import measure, exposure , segmentation
from scipy import ndimage
#---------------------------------------------------------------
# Images
path = 'Nifti_DataSet/COVID-19-CT-Seg_20cases'
Dataset = glob.glob( os.path.join(path, '*.gz') )
# Masks
masks_path = 'Nifti_DataSet/Lung_and_Infection_Mask'
masks_Dataset = glob.glob( os.path.join(masks_path, '*.gz') )
#---------------------------------------------------------------
# Source: https://www.kaggle.com/threedb/improved-lung-segmentation-using-watershed
#---------------------------------------------------------------
def generate_markers(image,thr):
    #Creation of the internal Marker
    marker_internal = (image < thr*image.max()).reshape(image.shape)
    marker_internal
    marker_internal = segmentation.clear_border(marker_internal)
    marker_internal_labels = measure.label(marker_internal)
    areas = [r.area for r in measure.regionprops(marker_internal_labels)]
    areas.sort()
    if len(areas) > 2:
        for region in measure.regionprops(marker_internal_labels):
            if region.area < areas[-2]:
                for coordinates in region.coords:                
                       marker_internal_labels[coordinates[0], coordinates[1]] = 0
    marker_internal = marker_internal_labels > 0
    #Creation of the external Marker
    external_a = ndimage.binary_dilation(marker_internal, iterations=10)
    external_b = ndimage.binary_dilation(marker_internal, iterations=55)
    marker_external = external_b ^ external_a
    #Creation of the Watershed Marker matrix
    marker_watershed = np.zeros((512, 512), dtype=np.int)
    marker_watershed += marker_internal * 255
    marker_watershed += marker_external * 128

    return marker_internal, marker_external, marker_watershed
#================================================================================        
def get_best_th(image):
    seg = np.zeros_like(image)
    for th in range(10):
        th = th/10
        m_internal,m_external, _ = generate_markers(image,th)
        if( (np.sum(m_external) > np.sum(m_internal)) and (np.sum(m_internal) > np.sum(seg))):
            seg = m_internal
            #print(th)
    return seg
#=====================================
def generate_watershed_segmentation(img):
    img2 = exposure.equalize_hist(img)
    segm = get_best_th(img2)
    segm = np.array(segm)
    segm = np.clip(segm,0,255)
    return segm
#===============================================================
# 
# Set a Counter 
ctr = 1
#---------------------------------------------------------------
for image in Dataset:
    # Load the images as nd-arrays
    images = nib.load(image).get_fdata()
    # Resize the images
    masks = resize(images, (256,256))
    # Load the masks as nd-arrays
    masks  = nib.load(os.path.join(masks_path,os.path.basename(image))).get_fdata() 
    # Resize the masks
    masks = resize(masks, (256,256))
    # make logical or for all the masks along slices axis
    masks = masks.any(axis=2)
    # print(masks.shape)
    # Take the Std Deviation for all the images along slices axis
    std = np.std(images, axis=2)
    stdcp = std.copy()
    # suppress the low values
    std[std < 0.1 * np.max(std)] = 0.0
    # Normalize The Output
    image = 1. * std / np.max(std)
    # Threshold the result
    val = threshold_otsu(image-std)
    image2 = (image-std) < val
    # Threshold again after masking with the previous result
    val2 = threshold_otsu(image2*image*(image-std))
    result = (image2*image*(image-std)) < val2
    #-------------------------------------------
    #           Morphology Part
    #===========================================
    # Structuring Element
    radius = 5
    # SE = morph.disk(radius)
    SE = morph.star(radius)
    # SE = morph.diamond(radius)
    #-------------------------------------------
    #        Morphological Operations
    #-------------------------------------------
    res2 = morph.binary_erosion(result,selem=SE)
    res2 = morph.remove_small_objects(res2)
    # res2 = morph.binary_opening(res2,selem=SE)
    res2 = morph.label(res2)
    #---------Largest Connected Component--------------------
    lab = np.argmax(np.unique(res2,return_counts=True)[1][1:])
    lab +=1
    res2[res2!=lab] = 0
    res2 = res2.astype(np.uint8)
    # ----More Morphological operations---------
    res2 = morph.remove_small_holes(res2)
    #-------------------------------------------
    res3 = generate_watershed_segmentation(stdcp)
    #-------------------------------------------
    #        Visualize the Results
    #-------------------------------------------
    plt.figure(num='image_%d'%ctr)

    plt.subplot('221')
    plt.imshow(masks)
    plt.title('masks after Logical OR')
    plt.axis('off')
    plt.subplot('222')
    plt.imshow(result)
    plt.title('after first stage')
    plt.axis('off')
    plt.subplot('223')
    plt.imshow(res2)
    plt.title('Result')
    plt.axis('off')   
    plt.subplot('224')
    plt.imshow(res3)
    plt.title('watershed')
    plt.axis('off')  
    # figManager = plt.get_current_fig_manager()
    # figManager.window.showMaximized()
    plt.show()
    #-----------------
    # Increase the counter
    ctr+=1