{"cells":[{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"# This Python 3 environment comes with many helpful analytics libraries installed\n# It is defined by the kaggle/python Docker image: https://github.com/kaggle/docker-python\n# For example, here's several helpful packages to load\n\nimport numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\n\n# Input data files are available in the read-only \"../input/\" directory\n# For example, running this (by clicking run or pressing Shift+Enter) will list all files under the input directory\n\nimport os\nimport glob\n\nimport os\nimport numpy as np\nimport pandas as pd\nimport pydicom as dcm\nimport matplotlib\nimport plotly.express as px\nimport matplotlib.pyplot as plt\nimport seaborn as sns\nimport glob\n# import gdcm\nfrom matplotlib import animation, rc\nfrom plotly.subplots import make_subplots\nimport plotly.graph_objs as go\n\nimport pydicom\nimport scipy.ndimage\n# import gdcm\nimport imageio\n\n\nimport os\nimport copy\nfrom datetime import timedelta, datetime\nimport imageio\nimport matplotlib.pyplot as plt\nfrom matplotlib import cm\nimport multiprocessing\nimport numpy as np\nimport os\nfrom pathlib import Path\nimport pydicom\nimport pytest\nimport scipy.ndimage as ndimage\nfrom scipy.ndimage.interpolation import zoom\nfrom skimage import measure, morphology, segmentation\nfrom skimage.transform import resize\nfrom time import time, sleep\nfrom tqdm import trange, tqdm\nimport tensorflow as tf\nfrom tensorflow.keras.models import Sequential\nfrom tensorflow.keras.preprocessing.image import ImageDataGenerator\nfrom tensorflow.keras.layers import *\nfrom tensorflow.data import Dataset\nimport torch\nfrom torch.utils.tensorboard import SummaryWriter\nfrom torchvision import transforms\nimport warnings\nimport seaborn as sns\nimport glob as glob\nimport imageio\nfrom IPython.display import Image\n\n#for masking\nfrom skimage.measure import label,regionprops\nfrom sklearn.cluster import KMeans\nfrom skimage.segmentation import clear_border\n\nimport onnx\n# +++++++++++++++?\n\nroot = \"/kaggle/input/rsna-str-pulmonary-embolism-detection/\"\n\nfor item in os.listdir(root):\n    path = os.path.join(root, item)\n    if os.path.isfile(path):\n        print(path)\n\n# You can write up to 5GB to the current directory (/kaggle/working/) that gets preserved as output when you create a version using \"Save & Run All\" \n# You can also write temporary files to /kaggle/temp/, but they won't be saved outside of the current session","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Training Data Visualization"},{"metadata":{"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","trusted":true},"cell_type":"code","source":"train_data = pd.read_csv(\"/kaggle/input/rsna-str-pulmonary-embolism-detection/train.csv\")\ntrain_data.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_data.shape","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_data.columns","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_data.info()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_data.describe()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"pd.isnull(train_data).any()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def rhead(x, nrow = 6, ncol = 4):\n    pd.set_option('display.expand_frame_repr', False)\n    seq = np.arange(0, len(x.columns), ncol)\n    for i in seq:\n        print(x.loc[range(0, nrow), x.columns[range(i, min(i+ncol, len(x.columns)))]])\n    pd.set_option('display.expand_frame_repr', True)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"rhead(train_data)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"for i in range(train_data.shape[1]-3):\n\n    train_data.hist(column=train_data.columns[i+3])","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### bar chart show train_data only at non-zero bars"},{"metadata":{"trusted":true},"cell_type":"code","source":"import matplotlib.pyplot as plt\nimport plotly.express as px\n\ntrain_data_drop = train_data.drop(['StudyInstanceUID', 'SeriesInstanceUID', 'SOPInstanceUID'], axis=1)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"\ntrain_data_select = train_data_drop.sum(axis=0).sort_values().reset_index()\ntrain_data_select.columns = ['columns', 'nonzero_records']\n\nfig = px.bar(\n    train_data_select, \n    x='nonzero_records', \n    y='columns', \n    orientation='h', \n    title='Columns and non zero samples', \n    height=800, \n    width=600\n)\nfig.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"\ntrain_data_select = train_data_drop.astype(bool).sum(axis=1).reset_index()\ntrain_data_select.columns = ['rows', 'count']\n\ntrain_data_select = train_data_select.groupby(['count'])['rows'].count().reset_index()\nfig = px.pie(\n    train_data_select, \n    values=round((100 * train_data_select['rows'] / len(train_data)), 2), \n    names=\"count\", \n    title='Every sample (Percent)', \n    width=500, \n    height=500\n)\nfig.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"\nf = plt.figure(figsize=(20, 20))\nplt.matshow(train_data_drop.corr(), fignum=f.number)\nplt.xticks(range(train_data_drop.shape[1]), train_data_drop.columns, fontsize=12, rotation=90)\nplt.yticks(range(train_data_drop.shape[1]), train_data_drop.columns, fontsize=12)\ncb = plt.colorbar()\ncb.ax.tick_params(labelsize=13)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_data_dir = \"../input/rsna-str-pulmonary-embolism-detection/train\"\nprint('Total number of dirictories in training set {}'.format(len(os.listdir(train_data_dir))))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Show a training Image"},{"metadata":{"trusted":true},"cell_type":"code","source":"# !pip install dicom","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"import vtk\nfrom vtk.util import numpy_support\nimport cv2\n\nreader = vtk.vtkDICOMImageReader()\ndef get_img(path):\n    reader.SetFileName(path)\n    reader.Update()\n    _extent = reader.GetDataExtent()\n    ConstPixelDims = [_extent[1]-_extent[0]+1, _extent[3]-_extent[2]+1, _extent[5]-_extent[4]+1]\n\n    ConstPixelSpacing = reader.GetPixelSpacing()\n    imageData = reader.GetOutput()\n    pointData = imageData.GetPointData()\n    arrayData = pointData.GetArray(0)\n    ArrayDicom = numpy_support.vtk_to_numpy(arrayData)\n    ArrayDicom = ArrayDicom.reshape(ConstPixelDims, order='F')\n    ArrayDicom = cv2.resize(ArrayDicom,(512,512))\n    return ArrayDicom","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"\nimport matplotlib.pyplot as plt\n\n\ndef show_dicom_images(dcom):\n    f, ax = plt.subplots(1,1, figsize=(12,10))\n    ax.imshow(dcom, cmap=plt.cm.bone)\n    ax.axis('off')\n    ax.set_title('Original DICOM Image')\n    plt.show()\n    \n#test read a dcom file and view it\nimg_path = \"../input/rsna-str-pulmonary-embolism-detection/train/005a0dbcb4b7/4ceaee66edc8/01a737504be7.dcm\"\nimg_get = get_img(img_path)\nshow_dicom_images(img_get)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### know a sample image"},{"metadata":{"trusted":true},"cell_type":"code","source":"import pydicom as dcm\nfig, ax = plt.subplots(figsize=(12, 12))\nax.imshow(dcm.dcmread(\"../input/rsna-str-pulmonary-embolism-detection/train/00c07cd8129d/8877e4d12ce9/00feb47a8d76.dcm\").pixel_array)\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### shape of a train image"},{"metadata":{"trusted":true},"cell_type":"code","source":"train_one_img = dcm.dcmread(\"../input/rsna-str-pulmonary-embolism-detection/train/00c07cd8129d/8877e4d12ce9/00feb47a8d76.dcm\").pixel_array\nprint(train_one_img.shape)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### info of an image"},{"metadata":{"trusted":true},"cell_type":"code","source":"train_one_img_info = dcm.dcmread(\"../input/rsna-str-pulmonary-embolism-detection/train/00c07cd8129d/8877e4d12ce9/00feb47a8d76.dcm\")\nprint(train_one_img_info)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### dcm images loading and preprocessing"},{"metadata":{"trusted":true},"cell_type":"code","source":"data_path = '../input/rsna-str-pulmonary-embolism-detection/train/0003b3d648eb/'\n\noutput_path = '../input/output/'\n\ntrain_image_files = sorted(glob.glob(os.path.join(data_path, '*','*.dcm')))\ntrain_image_files_list = os.listdir(data_path)\n\ntrain_image_files_list.sort()\n\nprint('Some sample ID''s :', len(train_image_files))\nprint(\"\\n\".join(train_image_files[:5]))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def load_scan(path):\n    \"\"\"\n    Loads scans from a folder and into a list.\n    \n    Parameters: path (Folder path)\n    \n    Returns: slices (List of slices)\n    \"\"\"\n    \n    slices = [pydicom.read_file(path + '/' + s) for s in os.listdir(path)]\n    slices.sort(key = lambda x: int(x.InstanceNumber))\n    \n    try:\n        slice_thickness = np.abs(slices[0].ImagePositionPatient[2] - slices[1].ImagePositionPatient[2])\n    except:\n        slice_thickness = np.abs(slices[0].SliceLocation - slices[1].SliceLocation)\n        \n    for s in slices:\n        s.SliceThickness = slice_thickness\n        \n    return slices","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"HUs can be calculated from the pixel data with a DICOM Image using the following formula:\n\nHU = m ∗ P + b\n\nwhere,\n\nm = RescaleSlope attribute of the DICOM image,\n\nb = RescaleIntercept attribute of the DICOM image,\n\nP = Pixel Array"},{"metadata":{"trusted":true},"cell_type":"code","source":"def get_pixels_hu(scans):\n    \"\"\"\n    Converts raw images to Hounsfield Units (HU).\n    \n    Parameters: scans (Raw images)\n    \n    Returns: image (NumPy array)\n    \"\"\"\n    \n    image = np.stack([s.pixel_array for s in scans])\n    image = image.astype(np.int16)\n\n    # Since the scanning equipment is cylindrical in nature and image output is square,\n    # we set the out-of-scan pixels to 0\n    image[image == -2000] = 0\n    \n    \n    # HU = m*P + b\n    intercept = scans[0].RescaleIntercept\n    slope = scans[0].RescaleSlope\n    \n    if slope != 1:\n        image = slope * image.astype(np.float64)\n        image = image.astype(np.int16)\n        \n    image += np.int16(intercept)\n    \n    return np.array(image, dtype=np.int16)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_img_scans = load_scan(data_path + train_image_files_list[0])\ntrain_images = get_pixels_hu(train_img_scans)\n\n#We'll be taking a random slice to perform segmentation:\n\nfor imgs in range(len(train_images[0:5])):\n    f, (ax1, ax2, ax3) = plt.subplots(1, 3, sharey=True, figsize=(15,15))\n    ax1.imshow(train_images[imgs], cmap=plt.cm.bone)\n    ax1.set_title(\"Original Slice\")\n    \n    ax2.imshow(train_images[imgs], cmap=plt.cm.bone)\n    ax2.set_title(\"Original Slice\")\n    \n    ax3.imshow(train_images[imgs], cmap=plt.cm.bone)\n    ax3.set_title(\"Original Slice\")\n    plt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"from IPython.display import Image\n\ndef set_lungwin(img, hu=[-1200., 600.]):\n    lungwin = np.array(hu)\n    newimg = (img-lungwin[0]) / (lungwin[1]-lungwin[0])\n    newimg[newimg < 0] = 0\n    newimg[newimg > 1] = 1\n    newimg = (newimg * 255).astype('uint8')\n    return newimg\n\n\nscans = load_scan(data_path + train_image_files_list[0])\nscan_array = set_lungwin(get_pixels_hu(scans))\n\nimageio.mimsave(\"/tmp/gif.gif\", scan_array, duration=0.0001)\nImage(filename=\"/tmp/gif.gif\", format='png')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### marker controlled watershed transformation"},{"metadata":{"trusted":true},"cell_type":"code","source":"def generate_markers(image):\n    \"\"\"\n    Generates markers for a given image.\n    \n    Parameters: image\n    \n    Returns: Internal Marker, External Marker, Watershed Marker\n    \"\"\"\n    \n    #Creation of the internal Marker\n    marker_internal = image < -400\n    marker_internal = segmentation.clear_border(marker_internal)\n    marker_internal_labels = measure.label(marker_internal)\n    \n    areas = [r.area for r in measure.regionprops(marker_internal_labels)]\n    areas.sort()\n    \n    if len(areas) > 2:\n        for region in measure.regionprops(marker_internal_labels):\n            if region.area < areas[-2]:\n                for coordinates in region.coords:                \n                       marker_internal_labels[coordinates[0], coordinates[1]] = 0\n    \n    marker_internal = marker_internal_labels > 0\n    \n    # Creation of the External Marker\n    external_a = ndimage.binary_dilation(marker_internal, iterations=10)\n    external_b = ndimage.binary_dilation(marker_internal, iterations=55)\n    marker_external = external_b ^ external_a\n    \n    # Creation of the Watershed Marker\n    marker_watershed = np.zeros((512, 512), dtype=np.int)\n    marker_watershed += marker_internal * 255\n    marker_watershed += marker_external * 128\n    \n    return marker_internal, marker_external, marker_watershed","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_img_scans = load_scan(data_path + train_image_files_list[0])\ntrain_images = get_pixels_hu(train_img_scans)\nprint(len(train_img_scans))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"\n\ntest_patient_internal, test_patient_external, test_patient_watershed = generate_markers(train_images[2])\n\nf, (ax1, ax2, ax3) = plt.subplots(1, 3, sharey=True, figsize=(15,15))\n\nax1.imshow(test_patient_internal, cmap='gray')\nax1.set_title(\"Internal Marker\")\nax1.axis('off')\n\nax2.imshow(test_patient_external, cmap='gray')\nax2.set_title(\"External Marker\")\nax2.axis('off')\n\nax3.imshow(test_patient_watershed, cmap='gray')\nax3.set_title(\"Watershed Marker\")\nax3.axis('off')\n\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Mask DICOM images"},{"metadata":{"trusted":true},"cell_type":"code","source":"sample_image = pydicom.dcmread(train_image_files[7])\nimg = sample_image.pixel_array\n\nplt.imshow(img, cmap='gray')\nplt.title('Original Image')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Now we will create a Binary Mask Image with rescale intercept and slope and adjusting values below -400 HU"},{"metadata":{"trusted":true},"cell_type":"code","source":"img = (img + sample_image.RescaleIntercept) / sample_image.RescaleSlope\nimg = img < -400 #HU unit range for lungs CT SCAN\n\nplt.imshow(img, cmap='gray')\nplt.title('Binary Mask Image')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Cleaning Border"},{"metadata":{"trusted":true},"cell_type":"code","source":"img = clear_border(img)\nplt.imshow(img, cmap='gray')\nplt.title('Cleaned Border Image')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Labelling a small region of scan"},{"metadata":{"trusted":true},"cell_type":"code","source":"img = label(img)\nplt.imshow(img, cmap='gray')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"areas = [r.area for r in regionprops(img)]\nareas.sort()\nif len(areas) > 2:\n    for region in regionprops(img):\n        if region.area < areas[-2]:\n            for coordinates in region.coords:                \n                img[coordinates[0], coordinates[1]] = 0\nimg = img > 0\nplt.imshow(img, cmap='gray')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Other Masks"},{"metadata":{"trusted":true},"cell_type":"code","source":"\ndef make_pemask(img, display=False):\n    row_size= img.shape[0]\n    col_size = img.shape[1]\n    \n    mean = np.mean(img)\n    std = np.std(img)\n    img = img-mean\n    img = img/std\n    \n    # Find the average pixel value near the lungs\n        # to renormalize washed out images\n    middle = img[int(col_size/5):int(col_size/5*4),int(row_size/5):int(row_size/5*4)] \n    mean = np.mean(middle)  \n    max = np.max(img)\n    min = np.min(img)\n    \n    # To improve threshold finding, I'm moving the \n    # underflow and overflow on the pixel spectrum\n    img[img==max]=mean\n    img[img==min]=mean\n    \n    # Using Kmeans to separate foreground (soft tissue / bone) and background (lung/air)\n    \n    kmeans = KMeans(n_clusters=2).fit(np.reshape(middle,[np.prod(middle.shape),1]))\n    centers = sorted(kmeans.cluster_centers_.flatten())\n    threshold = np.mean(centers)\n    thresh_img = np.where(img<threshold,1.0,0.0)  # threshold the image\n\n    # First erode away the finer elements, then dilate to include some of the pixels surrounding the lung.  \n    # We don't want to accidentally clip the lung.\n\n    eroded = morphology.erosion(thresh_img,np.ones([3,3]))\n    dilation = morphology.dilation(eroded,np.ones([8,8]))\n\n    labels = measure.label(dilation) # Different labels are displayed in different colors\n    label_vals = np.unique(labels)\n    regions = measure.regionprops(labels)\n    good_labels = []\n    for prop in regions:\n        B = prop.bbox\n        if B[2]-B[0]<row_size/10*9 and B[3]-B[1]<col_size/10*9 and B[0]>row_size/5 and B[2]<col_size/5*4:\n            good_labels.append(prop.label)\n    mask = np.ndarray([row_size,col_size],dtype=np.int8)\n    mask[:] = 0\n\n\n    #  After just the lungs are left, we do another large dilation\n    #  in order to fill in and out the lung mask \n    \n    for N in good_labels:\n        mask = mask + np.where(labels==N,1,0)\n    mask = morphology.dilation(mask,np.ones([10,10])) # one last dilation\n\n    if (display):\n        fig, ax = plt.subplots(3, 2, figsize=[12, 12])\n        ax[0, 0].set_title(\"Original\")\n        ax[0, 0].imshow(img, cmap='gray')\n        ax[0, 0].axis('off')\n        ax[0, 1].set_title(\"Threshold\")\n        ax[0, 1].imshow(thresh_img, cmap='gray')\n        ax[0, 1].axis('off')\n        ax[1, 0].set_title(\"After Erosion and Dilation\")\n        ax[1, 0].imshow(dilation, cmap='gray')\n        ax[1, 0].axis('off')\n        ax[1, 1].set_title(\"Color Labels\")\n        ax[1, 1].imshow(labels)\n        ax[1, 1].axis('off')\n        ax[2, 0].set_title(\"Final Mask\")\n        ax[2, 0].imshow(mask, cmap='gray')\n        ax[2, 0].axis('off')\n        ax[2, 1].set_title(\"Apply Mask on Original\")\n        ax[2, 1].imshow(mask*img, cmap='gray')\n        ax[2, 1].axis('off')\n        \n        plt.show()\n    return mask*img","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Select a sample\npath = \"../input/rsna-str-pulmonary-embolism-detection/train/000f7f114264/9f7378c3b2ab/0003aa3e734b.dcm\"\ndataset = pydicom.dcmread(path)\nimg = dataset.pixel_array\n\n# Masked image\nmask_img = make_pemask(img, display=True)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"import re\ntrain_img_dir = \"../input/rsna-str-pulmonary-embolism-detection/train/000f7f114264/9f7378c3b2ab\"\ndatasets = []\n\n# First Order the files in the dataset\nfiles = []\nfor dcm in list(os.listdir(train_img_dir)):\n    files.append(dcm) \nfiles.sort(key=lambda f: int(re.sub('\\D', '', f)))\n\n# Read in the Dataset\nfor dcm in files:\n    path = train_img_dir + \"/\" + dcm\n    datasets.append(pydicom.dcmread(path))\n    \nimgs = []\nfor data in datasets:\n    img = data.pixel_array\n    imgs.append(img)\n    \n    \n# Show masks\nfig=plt.figure(figsize=(16, 6))\ncolumns = 10\nrows = 3\n\nfor i in range(1, columns*rows +1):\n    img = make_pemask(datasets[i-1].pixel_array)\n    fig.add_subplot(rows, columns, i)\n    plt.imshow(img, cmap=\"gray\")\n    plt.title(i, fontsize = 9)\n    plt.axis('off');","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# To be continued ... "}],"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat":4,"nbformat_minor":4}