<jupyter_start><jupyter_text>I made this kernel for generating masks using deep learning. I was inspired by the discussion in this thread: https://www.kaggle.com/c/osic-pulmonary-fibrosis-progression/discussion/166123Using a pre-trained model, we can get a lung mask pretty easily. Let's check this out:<jupyter_code>#! pip install git+https://github.com/SoufianeDataFan/lungmask
from os import listdir, mkdir
import numpy as np
import pandas as pd
import pydicom
import seaborn as sns
import matplotlib.pyplot as plt
from tqdm import tqdm

import skimage
from skimage.measure import label,regionprops
from skimage.segmentation import clear_border
import matplotlib.pyplot as plt
from matplotlib import colors
from matplotlib import cm
import cv2
import math
#Interaction user
from ipywidgets import interact, interactive, fixed, interact_manual
import ipywidgets as widgets

from IPython.display import HTML
from IPython.display import clear_output
import warnings
warnings.simplefilter("ignore")
warnings.filterwarnings("ignore", category=DeprecationWarning)
warnings.filterwarnings("ignore", category=UserWarning)
warnings.filterwarnings("ignore", category=FutureWarning)

def selectDicom(x=6):
    print("El Dicom Escogido es: {}".format(x))
    return x
def selectPatient(x=0):
    print("Selected Patient: {}".format(x))
    return x
def selectSeries(x=0):
    print("Selected Serie: {}".format(x))
    return x
def selectDirectory(x="./"):
    print("Selected Directory: {}".format(x))
    return x
widgetSelDirectory = interactive(selectDirectory, x=["CD_2_CLEAN","CD_1_CLEAN"])
display(widgetSelDirectory)
basepathAqp = "../../datasetAqp/"+widgetSelDirectory.result
listEstudios = listdir(basepathAqp)
datasetAqp = pd.DataFrame(listEstudios, columns=["Estudio"])
estudios = []
for estudio in listEstudios:
    estudios.append(basepathAqp + "/"+estudio)
datasetAqp['fullPath'] = estudios
datasetAqp = datasetAqp.sort_values(by="Estudio", ascending=True, key=lambda col: col.str[5:],ignore_index = True)
datasetAqp
slices = []
def load_scans(dcm_path,stringInit):
    if stringInit  == "I":
        #ordenar por nro de orden del archivo -> porque no hay ImagePosition metadata por cual podriamos filtrar
        files = listdir(dcm_path)
        file_nums = [np.int(file.split(".")[0][1:]) for file  in files]
        sorted_file_nums = np.sort(file_nums)[::-1]
        
        #cortes ordenados y leidos se debe agregar un I al inicio del corte dicom (ojo no lleva extension .dcm)
        for file_num in tqdm(sorted_file_nums):
            if (file_num == 0):
                slices.append(pydicom.dcmread(dcm_path + "/I" + str(file_num) ))
            elif (int(math.log10(file_num))+1 == 1):
                slices.append(pydicom.dcmread(dcm_path + "/I" + str(file_num) ))
            else:
                slices.append(pydicom.dcmread(dcm_path + "/I" + str(file_num) ))
    elif stringInit == "Z":
        #ordenar por nro de orden del archivo -> porque no hay ImagePosition metadata por cual podriamos filtrar
        files = listdir(dcm_path)
        file_nums = [np.int(file.split(".")[0][1:]) for file  in files]
        sorted_file_nums = np.sort(file_nums)[::-1]
        #cortes ordenados y leidos se debe agregar un I al inicio del corte dicom (ojo no lleva extension .dcm)
        for file_num in tqdm(sorted_file_nums):
            if int(math.log10(file_num))+1 == 1:
                slices.append(pydicom.dcmread(dcm_path + "/Z0" + str(file_num)))
            else:
                slices.append(pydicom.dcmread(dcm_path + "/Z" + str(file_num)))
    return slices<jupyter_output><empty_output><jupyter_text>Widget Control Parameters<jupyter_code>widgetSelectionDicom = interactive(selectDicom,x=(0,20))
display(widgetSelectionDicom)
#muestra de un estudio dado
example = datasetAqp.fullPath.values[widgetSelectionDicom.result]
scans = load_scans(example,"I")
print("El estudio contiene {} cortes".format(len(scans)))
estudio_sel = example[28:]
estudio_sel
df = pd.DataFrame(columns=["PatientID","StudyDescription","SeriesNumber","SeriesDescription","InstanceNumber","SeriesInstanceUID"])
#df = pd.DataFrame(columns=["PatientID"])
for scan in tqdm(scans):
    #serie
    dicSerie = {"PatientID":scan.PatientID,
                "StudyDescription":scan.StudyDescription,
                "SeriesNumber":scan.SeriesNumber,
                "SeriesDescription":scan.SeriesDescription,
                "InstanceNumber":scan.InstanceNumber,
                "SeriesInstanceUID":scan.SeriesInstanceUID
               }
    df = df.append(dicSerie, ignore_index=True)
df
# Oonly one pacient with different type of series
listPatients = df.PatientID.unique()
listPatients
widgetSelPatient = interactive(selectPatient, x=listPatients)
display(widgetSelPatient)
#Series List
listSeries = df[df["PatientID"]==widgetSelPatient.result]["SeriesNumber"].unique()
listSeries
#Description of every serie
listSeriesDescription = df[df["PatientID"]==listPatients[0]]["SeriesDescription"].unique()
listSeriesDescription
widgetSelSerie = interactive(selectSeries, x=listSeries)
display(widgetSelSerie)
#Select Slice AXIAL S/C PULM which is parenquimal slice
dicomAislada = df[(df["PatientID"]==widgetSelPatient.result) & (df["SeriesNumber"]==widgetSelSerie.result)]
dicomAislada
serieID_select = dicomAislada["SeriesInstanceUID"][dicomAislada.index.min()]
serieID_select
print("Min Número de Corte: {}".format(dicomAislada.index.min()))
print("Max Número de Corte: {}".format(dicomAislada.index.max()))
scansFiltrado = scans[dicomAislada.index.min():dicomAislada.index.max()+1]
orderScans = sorted(scansFiltrado, key=lambda corte: corte.InstanceNumber)<jupyter_output><empty_output><jupyter_text>Lectura de pixels y transformacion a HU<jupyter_code>#definiendo el aire con valores cero (pixeles completamente oscuros vacio)
def set_outside_scanner_to_air(raw_pixelarrays):
    #filtrando pintando todo negro el aire
    raw_pixelarrays[raw_pixelarrays <= -1000] = 0
    return raw_pixelarrays
#tranformar a unidades de hounsfield
def transform_to_hu(slices):
    images = np.stack([file.pixel_array for file in slices])
    #convertido a escala de grises de 16 bits
    images = images.astype(np.int16)
    #aplicar pintado de aire de todos los cortes apilados <= -1000
    images = set_outside_scanner_to_air(images)
    
    #convertirmos a HU
    for n in range(len(slices)):
        intercept = slices[n].RescaleIntercept
        slope = slices[n].RescaleSlope
        
        #FORMULA = mb + i donde m slope o pendiente e i es el intercepto
        
        #MULTIPLICAR EL SLOPE
        #solo en caso de q el slope sea 0.5, 2, 0.9 
        if slope != 1:
            #para ser multiplicado por slope debe pasarse a flotante
            images[n] = slope * images[n].astype(np.float64)
            #despues ya convertimos a la escala de grise 16 bits
            images[n] = images[n].astype(np.int16)
        
        #ADICIONAR EL INTERCEPTO
        #directamente el intercepto en caso el slope sea 1
        #trata como escala de 16 bits el intercepto
        images[n] += np.int16(intercept)
        
    return np.array(images, dtype=np.int16)
hu_scans = transform_to_hu(orderScans)
print("Total Nro de Cortes: {}".format(len(hu_scans)))
fix, ax =  plt.subplots(1,4,figsize=(20,5))
ax[0].set_title("Original CT-Scan")
ax[0].imshow(orderScans[30].pixel_array, cmap="bone")
ax[1].set_title("PixelArray Distribution")
sns.distplot(orderScans[30].pixel_array.flatten(),ax=ax[1])

ax[2].set_title("Ct_Scan in HU")
ax[2].imshow(hu_scans[30],cmap='bone')
ax[3].set_title("HU distribution values")
sns.distplot(hu_scans[30].flatten(), ax=ax[3])

for m in [0,2]:
    ax[m].grid(False)
# import matplotlib.animation as animation

# fig = plt.figure()

# ims = []
# for image in orderScans:
#     newimage = image.pixel_array
#     im = plt.imshow(newimage, animated=True, cmap="Greys")
#     plt.axis("off")
#     ims.append([im])

# ani = animation.ArtistAnimation(fig, ims, interval=100, blit=False,
#                                 repeat_delay=1000)
# HTML(ani.to_jshtml())<jupyter_output><empty_output><jupyter_text>**Original Lung CT-Scan**<jupyter_code># #Aplicar UH al dicom completo
# def get_pixels_hu(slices):
#     image = np.stack([s.pixel_array for s in slices])
#     # Convert to int16 (from sometimes int16), 
#     # should be possible as values should always be low enough (<32k)
#     image = image.astype(np.int16)

#     # Set outside-of-scan pixels to 0
#     # The intercept is usually -1024, so air is approximately 0
#     image[image == -2000] = 0
    
#     # Convert to Hounsfield units (HU)
#     for slice_number in range(len(slices)):
        
#         intercept = slices[slice_number].RescaleIntercept
#         slope = slices[slice_number].RescaleSlope
        
#         if slope != 1:
#             image[slice_number] = slope * image[slice_number].astype(np.float64)
#             image[slice_number] = image[slice_number].astype(np.int16)
            
#         image[slice_number] += np.int16(intercept)
    
#     return np.array(image, dtype=np.int16)
# #mascara creada a partir del intervalo de UH ya definido y con la
# #imagen ya aplicada la transformacion de UH
# def set_lungwin(img, hu=[-1200., 600.]):
#     lungwin = np.array(hu)
#     newimg = (img-lungwin[0]) / (lungwin[1]-lungwin[0])
#     newimg[newimg < 0] = 0
#     newimg[newimg > 1] = 1
#     newimg = (newimg * 255).astype('uint8')
#     return newimg
# scan_array = set_lungwin(get_pixels_hu(scans))<jupyter_output><empty_output><jupyter_text>**Mask using U-net(R231) Lung CT-Scan**<jupyter_code>from lungmask import mask
import SimpleITK as sitk
import sys
import os
# #agregando una dimension adicional
# arrayImage = np.expand_dims(im, axis=0)
# #pasando a simple itk
# print(arrayImage.shape)
# input_image = sitk.GetImageFromArray(arrayImage)
# #volviendo a array numpy
# arrayImage = sitk.GetArrayFromImage(input_image).astype(np.int32)
# #arrayImage.shape

# z = 0
# sliceImage = sitk.GetArrayFromImage(input_image)[z,:,:,:]
# plt.imshow(sliceImage)
mascarasPred = []
def get_mascaras(slices_hu):
    
    for corte in tqdm(slices_hu):   
        #agregando una dimension adicional
        arrayImage = np.expand_dims(corte.astype(np.float32)  , axis=0)
        #pasando a simple itk
        input_image = sitk.GetImageFromArray(arrayImage)      
        
        mascarasPred.append(mask.apply(input_image)[0])
    return mascarasPred
len(orderScans)
mascaras = get_mascaras(hu_scans[::-1])
fig = plt.figure(figsize=(12, 12))
plt.imshow(mascaras[15],norm=colors.Normalize(0,1))
# import matplotlib.animation as animation

# fig = plt.figure()

# ims = []
# for image in mascaras:
#     im = plt.imshow(image, animated=True, cmap="Greys")
#     plt.axis("off")
#     ims.append([im])

# ani = animation.ArtistAnimation(fig, ims, interval=100, blit=False,
#                                 repeat_delay=1000)
# HTML(ani.to_jshtml())<jupyter_output><empty_output><jupyter_text>Saving mask<jupyter_code>print("Reading Dicom directory:", basepathAqp +"/"+ estudio_sel)
print("Con el ID de serie: ",serieID_select)
reader = sitk.ImageSeriesReader()

dicom_names = reader.GetGDCMSeriesFileNames(basepathAqp + "/" + estudio_sel, serieID_select)
#dicom_names = reader.GetGDCMSeriesFileNames(basepathAqp + "/" + estudio_sel)
reader.SetFileNames(dicom_names)

imageDicom = reader.Execute()
size = imageDicom.GetSize()
print("Image size:", imageDicom.GetSize())
print("Image direction:", imageDicom.GetDirection())
print("Image dimension:", imageDicom.GetDimension())
print("Image origin:", imageDicom.GetOrigin())
len(mascaras)
maskAsImage = np.zeros((len(mascaras),512,512))
for idxMask in range(len(mascaras)):
    maskAsImage[idxMask,:,:] = mascaras[idxMask]
newMask = sitk.GetImageFromArray(maskAsImage)
newMask.CopyInformation(imageDicom)
print("Image size:", newMask.GetSize())
print("Image direction:", newMask.GetDirection())
print("Image dimension:", newMask.GetDimension())
print("Image origin:", newMask.GetOrigin())
writer = sitk.ImageFileWriter()
writer.SetFileName("maskAsImage_"+estudio_sel+".nii.gz")
writer.Execute(newMask);
#dice e IOU sin clase

def dice_coef_metric(probabilities: np.ndarray,
                     truth: np.ndarray,
                     treshold: float = 0.5,
                     eps: float = 1e-9) -> np.ndarray:
    """
    Calculate Dice score for data batch.
    Params:
        probobilities: model outputs after activation function.
        truth: truth values.
        threshold: threshold for probabilities.
        eps: additive to refine the estimate.
        Returns: dice score aka f1.
    """
    scores = []
    num = probabilities.shape[0]
    
    predictions = (probabilities >= treshold)
    
    assert(predictions.shape == truth.shape)
    #recorremos el nro de cortes
    for i in range(num):
        prediction = predictions[i]
        truth_ = truth[i]
        intersection = 2.0 * (truth_ * prediction).sum()
        union = truth_.sum() + prediction.sum()
        if truth_.sum() == 0 and prediction.sum() == 0:
            scores.append(1.0)
        else:
            scores.append((intersection + eps) / union)
    return np.mean(scores)


def jaccard_coef_metric(probabilities: np.ndarray,
               truth: np.ndarray,
               treshold: float = 0.5,
               eps: float = 1e-9) -> np.ndarray:
    """
    Calculate Jaccard index for data batch.
    Params:
        probobilities: model outputs after activation function.
        truth: truth values.
        threshold: threshold for probabilities.
        eps: additive to refine the estimate.
        Returns: jaccard score aka iou."
    """
    scores = []
    num = probabilities.shape[0]
    predictions = (probabilities >= treshold)
    assert(predictions.shape == truth.shape)
    
    #recorremos nro de cortes
    for i in range(num):
        prediction = predictions[i]
        truth_ = truth[i]
        intersection = (prediction * truth_).sum()
        union = (prediction.sum() + truth_.sum()) - intersection + eps
        if truth_.sum() == 0 and prediction.sum() == 0:
            scores.append(1.0)
        else:
            scores.append((intersection + eps) / union)
    return np.mean(scores)
#reading manual segmentation
reader = sitk.ImageFileReader()
reader.SetImageIO("NiftiImageIO")
reader.SetFileName("../../manualSegm/nifti/"+estudio_sel+".nii")
manualSegm = reader.Execute()
manualSegmarray = sitk.GetArrayFromImage(manualSegm)
print(manualSegmarray.shape)
print(maskAsImage.shape)
print("Dice Coeff: {}".format(dice_coef_metric(maskAsImage,manualSegmarray)))
print("Jaccard Coeff: {}".format(jaccard_coef_metric(maskAsImage,manualSegmarray)))<jupyter_output>Jaccard Coeff: 0.8536947183981938<jupyter_text>Lectura<jupyter_code>#Lectura de las segmentaciones manuales
metrics = []
initName_manualSegm = "COD00"
estudios_sel = []
not_studies = [3,4,8,10,12,13,17,18,20,26,28]
for i in range(1,30):
    if (i not in not_studies):
        if(int(math.log10(i))+1 == 1):
            estudios_sel.append(initName_manualSegm+"0"+str(i))
        else:
            estudios_sel.append(initName_manualSegm+str(i))


reader = sitk.ImageFileReader()
reader.SetImageIO("NiftiImageIO")
### lectura de la segm manual
for estudio in tqdm(estudios_sel):
    
    reader.SetFileName("../../manualSegm/nifti/"+estudio+".nii")
    manualSemgs = sitk.GetArrayFromImage(reader.Execute())

    reader.SetFileName("saveSegm/maskAsImage_"+estudio+".nii.gz")
    unetSegms = sitk.GetArrayFromImage(reader.Execute())

    metrics.append((estudio[-2:],dice_coef_metric(unetSegms,manualSemgs),jaccard_coef_metric(unetSegms,manualSemgs)))
import seaborn as sns

def show_values_on_bars(axs, h_v="v", space=0.4):
    def _show_on_single_plot(ax):
        if h_v == "v":
            for p in ax.patches:
                _x = p.get_x() + p.get_width() / 2
                _y = p.get_y()+0.01 + p.get_height()
                value = '{:.2f}'.format(p.get_height())
                ax.text(_x, _y, value, ha="center", fontsize=12) 
        elif h_v == "h":
            for p in ax.patches:
                _x = p.get_x() + p.get_width() + float(space)
                _y = p.get_y() + p.get_height()
                value = p.get_width()
                ax.text(_x, _y, value, ha="left")

    if isinstance(axs, np.ndarray):
        for idx, ax in np.ndenumerate(axs):
            _show_on_single_plot(ax)
    else:
        _show_on_single_plot(axs)
        
valDicom = []
valMetric_dice = []
valMetric_jaccard = []
for metric in metrics:
    valDicom.append(metric[0])
    valMetric_dice.append(metric[1])
    valMetric_jaccard.append(metric[2])
dataMetrics = {"dicom": valDicom, "dice": valMetric_dice, "jaccard": valMetric_jaccard}
#data
dfMetric = pd.DataFrame(dataMetrics)
#dfOrderMetric = dfMetric.sort_values(by=['dice_metrics'])

plt.figure(figsize=(12,5))
ax = sns.barplot(x="dicom", y="dice", data=dfMetric,palette="flare")
ax = sns.barplot(x="dicom", y="jaccard", data=dfMetric,palette="crest")
plt.setp(ax.get_xticklabels(), rotation=0)
ax.tick_params(labelsize=14)
ax.tick_params(length=0, axis='x')
ax.set_ylabel("Dice & Jaccard Coefficient", fontsize=12)
ax.set_xlabel("ID Dicoms", fontsize=12)
ax.set_title("Dice & Jaccard Metrics between Unet Segmentation and Manual Segmentation")
show_values_on_bars(ax, "v", 0.3)
plt.tight_layout()
plt.show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.12526539278131635, 'fraction_numerical': 0.015493200206576003, 'mean_word_length': 3.8024249104436483, 'pattern_counts': {'":': 13, '<': 32, '<?xml version=': 0, '>': 35, 'https://': 2, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '7', 'dir': 'jupyter-structured-clean-dedup', 'id': '561568', 'license': 'Apache-2.0', 'n_tokens_mistral': 6200, 'n_tokens_neox': 5949, 'n_words': 1264, 'path': 'lung-mask-using-u-net-r231_v3.ipynb', 'repo_name': 'christiansucavelando/lungmask', 'size': '17427'}