{"cells":[{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"!conda install -c conda-forge gdcm -y","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","trusted":true},"cell_type":"code","source":"import numpy as np \nimport pandas as pd \nimport matplotlib.pyplot as plt\nimport seaborn as sns\n\nimport pydicom\nimport scipy.ndimage\nimport gdcm\n\nfrom skimage import measure \nfrom mpl_toolkits.mplot3d.art3d import Poly3DCollection\nfrom skimage.morphology import disk, opening, closing\nfrom tqdm import tqdm\n\nfrom IPython.display import HTML\nfrom PIL import Image\n\nimport warnings\nwarnings.filterwarnings(\"ignore\", category=DeprecationWarning)\nwarnings.filterwarnings(\"ignore\", category=UserWarning)\nwarnings.filterwarnings(\"ignore\", category=FutureWarning)\n\nfrom os import listdir,mkdir","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"listdir(\"../input/\")\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#basepath = \"../input/osic-pulmonary-fibrosis-progression/\"\n# or if you are taking part in RSNA pulmonary embolism detection:\nbasepath = \"../input/rsna-str-pulmonary-embolism-detection/\"\nlistdir(basepath)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train = pd.read_csv(basepath + \"train.csv\")\ntest = pd.read_csv(basepath + \"test.csv\")\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train.shape","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"if basepath == \"../input/osic-pulmonary-fibrosis-progression/\":\n    train[\"dcm_path\"] = basepath + \"train/\" + train.Patient + \"/\"\nelse:\n    train[\"dcm_path\"] = basepath + \"train/\" + train.StudyInstanceUID + \"/\" + train.SeriesInstanceUID  ","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"> #SORTING DONE TO LOAD DATA "},{"metadata":{"trusted":true},"cell_type":"code","source":"def load_scans(dcm_path):\n    if basepath == \"../input/osic-pulmonary-fibrosis-progression/\":\n        # in this competition we have missing values in ImagePosition, this is why we are sorting by filename number\n        files = listdir(dcm_path)\n        file_nums = [np.int(file.split(\".\")[0]) for file in files]\n        sorted_file_nums = np.sort(file_nums)[::-1]\n        slices = [pydicom.dcmread(dcm_path + \"/\" + str(file_num) + \".dcm\" ) for file_num in sorted_file_nums]\n    else:\n        # otherwise we sort by ImagePositionPatient (z-coordinate) or by SliceLocation\n        slices = [pydicom.dcmread(dcm_path + \"/\" + file) for file in listdir(dcm_path)]\n        slices.sort(key = lambda x: float(x.ImagePositionPatient[2]))\n    return slices","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"example = train.dcm_path.values[0]\nscans = load_scans(example)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"CHECK FIRST IMAGE"},{"metadata":{"trusted":true},"cell_type":"code","source":"scans[0]","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We need to transform to Hounsfield units as the spectral composition of the x-rays depends on the measurement settings like acquisition parameters and tube voltage. By normalizing to values of water and air (water has HU 0 and air -1000) the images of different measurements are becoming comparable"},{"metadata":{"trusted":true},"cell_type":"code","source":"fig, ax = plt.subplots(1,2,figsize=(20,5))\nfor n in range(10):\n    image = scans[n].pixel_array.flatten()\n    rescaled_image = image * scans[n].RescaleSlope + scans[n].RescaleIntercept\n    sns.distplot(image.flatten(), ax=ax[0]);\n    sns.distplot(rescaled_image.flatten(), ax=ax[1])\nax[0].set_title(\"Raw pixel array distributions for 10 examples\")\nax[1].set_title(\"HU unit distributions for 10 examples\");","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"For some examples we can see that there are raw values at -2000. They correspond to images with a circular boundary within the image. The \"outside\" of this circle value is often set to -2000"},{"metadata":{"trusted":true},"cell_type":"code","source":"def set_outside_scanner_to_air(raw_pixelarrays):\n    # in OSIC we find outside-scanner-regions with raw-values of -2000. \n    # Let's threshold between air (0) and this default (-2000) using -1000\n    raw_pixelarrays[raw_pixelarrays <= -1000] = 0\n    return raw_pixelarrays","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def transform_to_hu(slices):\n    images = np.stack([file.pixel_array for file in slices])\n    images = images.astype(np.int16)\n\n    images = set_outside_scanner_to_air(images)\n    \n    # convert to HU\n    for n in range(len(slices)):\n        \n        intercept = slices[n].RescaleIntercept\n        slope = slices[n].RescaleSlope\n        \n        if slope != 1:\n            images[n] = slope * images[n].astype(np.float64)\n            images[n] = images[n].astype(np.int16)\n            \n        images[n] += np.int16(intercept)\n    \n    return np.array(images, dtype=np.int16)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Now all raw values per slice are scaled to H-units."},{"metadata":{"trusted":true},"cell_type":"code","source":"hu_scans = transform_to_hu(scans)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"fig, ax = plt.subplots(1,4,figsize=(20,3))\nax[0].set_title(\"Original CT-scan\")\nax[0].imshow(scans[0].pixel_array, cmap=\"bone\")\nax[1].set_title(\"Pixelarray distribution\");\nsns.distplot(scans[0].pixel_array.flatten(), ax=ax[1]);\n\nax[2].set_title(\"CT-scan in HU\")\nax[2].imshow(hu_scans[0], cmap=\"bone\")\nax[3].set_title(\"HU values distribution\");\nsns.distplot(hu_scans[0].flatten(), ax=ax[3]);\n\nfor m in [0,2]:\n    ax[m].grid(False)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"N = 1000","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":""},{"metadata":{"trusted":true},"cell_type":"code","source":"def get_window_value(feature):\n    if type(feature) == pydicom.multival.MultiValue:\n        return np.int(feature[0])\n    else:\n        return np.int(feature)\n\npixelspacing_r = []\npixelspacing_c = []\nslice_thicknesses = []\npatient_id = []\npatient_pth = []\nrow_values = []\ncolumn_values = []\nwindow_widths = []\nwindow_levels = []\n\nif basepath == \"../input/osic-pulmonary-fibrosis-progression/\":\n    patients = train.Patient.unique()[0:N]\nelse:\n    patients = train.SeriesInstanceUID.unique()[0:N]\n\nfor patient in patients:\n    patient_id.append(patient)\n    if basepath == \"../input/osic-pulmonary-fibrosis-progression/\":\n        path = train[train.Patient == patient].dcm_path.values[0]\n    else:\n        path = train[train.SeriesInstanceUID == patient].dcm_path.values[0]\n    example_dcm = listdir(path)[0]\n    patient_pth.append(path)\n    dataset = pydicom.dcmread(path + \"/\" + example_dcm)\n    \n    window_widths.append(get_window_value(dataset.WindowWidth))\n    window_levels.append(get_window_value(dataset.WindowCenter))\n    \n    spacing = dataset.PixelSpacing\n    slice_thicknesses.append(dataset.SliceThickness)\n    \n    row_values.append(dataset.Rows)\n    column_values.append(dataset.Columns)\n    pixelspacing_r.append(spacing[0])\n    pixelspacing_c.append(spacing[1])\n    \nscan_properties = pd.DataFrame(data=patient_id, columns=[\"patient\"])\nscan_properties.loc[:, \"rows\"] = row_values\nscan_properties.loc[:, \"columns\"] = column_values\nscan_properties.loc[:, \"area\"] = scan_properties[\"rows\"] * scan_properties[\"columns\"]\nscan_properties.loc[:, \"pixelspacing_r\"] = pixelspacing_r\nscan_properties.loc[:, \"pixelspacing_c\"] = pixelspacing_c\nscan_properties.loc[:, \"pixelspacing_area\"] = scan_properties.pixelspacing_r * scan_properties.pixelspacing_c\nscan_properties.loc[:, \"slice_thickness\"] = slice_thicknesses\nscan_properties.loc[:, \"patient_pth\"] = patient_pth\nscan_properties.loc[:, \"window_width\"] = window_widths\nscan_properties.loc[:, \"window_level\"] = window_levels\nscan_properties.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"pixelspacing attribute-e. It tells us how much physical distance is covered by one pixel."},{"metadata":{},"cell_type":"markdown","source":"between patients the pixelspacing can differ due to personal or institutional preferences of doctors and the clinic and it also depends on the scanner type. Consequently if you compare two images in the size of the lungs it does not automatically mean that the bigger one is really larger in the physical size of the organ!"},{"metadata":{"trusted":true},"cell_type":"code","source":"fig, ax = plt.subplots(1,2,figsize=(20,5))\nsns.distplot(pixelspacing_r, ax=ax[0], color=\"Limegreen\", kde=False)\nax[0].set_title(\"Pixel spacing distribution \\n in row direction \")\nax[0].set_ylabel(\"Counts in train\")\nax[0].set_xlabel(\"mm\")\nsns.distplot(pixelspacing_c, ax=ax[1], color=\"Mediumseagreen\", kde=False)\nax[1].set_title(\"Pixel spacing distribution \\n in column direction\");\nax[1].set_ylabel(\"Counts in train\");\nax[1].set_xlabel(\"mm\");","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"PIXEL SPACING IS VARYING A LOT ! NEED TO IMPROVE"},{"metadata":{"trusted":true},"cell_type":"code","source":"counts = scan_properties.groupby([\"rows\", \"columns\"]).size()\ncounts = counts.unstack()\ncounts.fillna(0, inplace=True)\n\n\nfig, ax = plt.subplots(1,2,figsize=(20,5))\nsns.distplot(slice_thicknesses, color=\"orangered\", kde=False, ax=ax[0])\nax[0].set_title(\"Slice thicknesses of all patients\");\nax[0].set_xlabel(\"Slice thickness in mm\")\nax[0].set_ylabel(\"Counts in train\");\n\nfor n in counts.index.values:\n    for m in counts.columns.values:\n        ax[1].scatter(n, m, s=counts.loc[n,m], c=\"midnightblue\")\nax[1].set_xlabel(\"rows\")\nax[1].set_ylabel(\"columns\")\nax[1].set_title(\"Pixel area of ct-scan per patient\");","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The slice thickness tells us how much distance is covered in Z-direction by one slic"},{"metadata":{},"cell_type":"markdown","source":"pixel_array of raw values covers a specific area given by row and column values"},{"metadata":{},"cell_type":"markdown","source":"THIN SLICES IN CT SCAN MEANS MORE DETAIL"},{"metadata":{"trusted":true},"cell_type":"code","source":"scan_properties[\"r_distance\"] = scan_properties.pixelspacing_r * scan_properties.rows\nscan_properties[\"c_distance\"] = scan_properties.pixelspacing_c * scan_properties[\"columns\"]\nscan_properties[\"area_cm2\"] = 0.1* scan_properties[\"r_distance\"] * 0.1*scan_properties[\"c_distance\"]\nscan_properties[\"slice_volume_cm3\"] = 0.1*scan_properties.slice_thickness * scan_properties.area_cm2","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We have some images with extreme large sliche areas and volumes"},{"metadata":{"trusted":true},"cell_type":"code","source":"fig, ax = plt.subplots(1,2,figsize=(20,5))\nsns.distplot(scan_properties.area_cm2, ax=ax[0], color=\"purple\")\nsns.distplot(scan_properties.slice_volume_cm3, ax=ax[1], color=\"magenta\")\nax[0].set_title(\"CT-slice area in $cm^{2}$\")\nax[1].set_title(\"CT-slice volume in $cm^{3}$\")\nax[0].set_xlabel(\"$cm^{2}$\")\nax[1].set_xlabel(\"$cm^{3}$\");","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"max_path = scan_properties[\n    scan_properties.area_cm2 == scan_properties.area_cm2.max()].patient_pth.values[0]\nmin_path = scan_properties[\n    scan_properties.area_cm2 == scan_properties.area_cm2.min()].patient_pth.values[0]\n\nmin_scans = load_scans(min_path)\nmin_hu_scans = transform_to_hu(min_scans)\n\nmax_scans = load_scans(max_path)\nmax_hu_scans = transform_to_hu(max_scans)\n\nbackground_water_hu_scans = max_hu_scans.copy()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def set_manual_window(hu_image, custom_center, custom_width):\n    w_image = hu_image.copy()\n    min_value = custom_center - (custom_width/2)\n    max_value = custom_center + (custom_width/2)\n    w_image[w_image < min_value] = min_value\n    w_image[w_image > max_value] = max_value\n    return w_image","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Smallest and larges CT-slice area"},{"metadata":{"trusted":true},"cell_type":"code","source":"fig, ax = plt.subplots(1,2,figsize=(20,10))\nax[0].imshow(set_manual_window(min_hu_scans[np.int(len(min_hu_scans)/2)], -700, 255), cmap=\"YlGnBu\")\nax[1].imshow(set_manual_window(max_hu_scans[np.int(len(max_hu_scans)/2)], -700, 255), cmap=\"YlGnBu\");\nax[0].set_title(\"CT-scan with small slice area\")\nax[1].set_title(\"CT-scan with large slice area\");\nfor n in range(2):\n    ax[n].axis(\"off\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"fig, ax = plt.subplots(1,2,figsize=(20,5))\nsns.distplot(max_hu_scans[np.int(len(max_hu_scans)/2)].flatten(), kde=False, ax=ax[1])\nax[1].set_title(\"Large area image\")\nsns.distplot(min_hu_scans[np.int(len(min_hu_scans)/2)].flatten(), kde=False, ax=ax[0])\nax[0].set_title(\"Small area image\")\nax[0].set_xlabel(\"HU values\")\nax[1].set_xlabel(\"HU values\");","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"we can see that the large one has a lot of useless region covered. We could crop it."},{"metadata":{"trusted":true},"cell_type":"code","source":"max_path = scan_properties[\n    scan_properties.slice_volume_cm3 == scan_properties.slice_volume_cm3.max()].patient_pth.values[0]\nmin_path = scan_properties[\n    scan_properties.slice_volume_cm3 == scan_properties.slice_volume_cm3.min()].patient_pth.values[0]\n\nmin_scans = load_scans(min_path)\nmin_hu_scans = transform_to_hu(min_scans)\n\nmax_scans = load_scans(max_path)\nmax_hu_scans = transform_to_hu(max_scans)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"fig, ax = plt.subplots(1,2,figsize=(20,10))\nax[0].imshow(set_manual_window(min_hu_scans[np.int(len(min_hu_scans)/2)], -700, 255), cmap=\"YlGnBu\")\nax[1].imshow(set_manual_window(max_hu_scans[np.int(len(max_hu_scans)/2)], -700, 255), cmap=\"YlGnBu\");\nax[0].set_title(\"CT-scan with small slice volume\")\nax[1].set_title(\"CT-scan with large slice volume\");\nfor n in range(2):\n    ax[n].axis(\"off\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"fig, ax = plt.subplots(1,2,figsize=(20,5))\nsns.distplot(max_hu_scans[np.int(len(max_hu_scans)/2)].flatten(), kde=False, ax=ax[1])\nax[1].set_title(\"Large slice volume\")\nsns.distplot(min_hu_scans[np.int(len(min_hu_scans)/2)].flatten(), kde=False, ax=ax[0])\nax[0].set_title(\"Small slice volume\")\nax[0].set_xlabel(\"HU values\")\nax[1].set_xlabel(\"HU values\");","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"3D-reconstruction of CT-scans"},{"metadata":{"trusted":true},"cell_type":"code","source":"def plot_3d(image, threshold=700, color=\"navy\"):\n    \n    # Position the scan upright, \n    # so the head of the patient would be at the top facing the camera\n    p = image.transpose(2,1,0)\n    \n    verts, faces,_,_ = measure.marching_cubes_lewiner(p, threshold)\n\n    fig = plt.figure(figsize=(10, 10))\n    ax = fig.add_subplot(111, projection='3d')\n\n    # Fancy indexing: `verts[faces]` to generate a collection of triangles\n    mesh = Poly3DCollection(verts[faces], alpha=0.2)\n    mesh.set_facecolor(color)\n    ax.add_collection3d(mesh)\n\n    ax.set_xlim(0, p.shape[0])\n    ax.set_ylim(0, p.shape[1])\n    ax.set_zlim(0, p.shape[2])\n\n    plt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plot_3d(max_hu_scans)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"old_distribution = max_hu_scans.flatten()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"example = train.dcm_path.values[0]\nscans = load_scans(example)\nhu_scans = transform_to_hu(scans)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plot_3d(hu_scans)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":""},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.figure(figsize=(20,5))\nsns.distplot(old_distribution, label=\"weak 3d plot\", kde=False)\nsns.distplot(hu_scans.flatten(), label=\"strong 3d plot\", kde=False)\nplt.title(\"HU value distribution\")\nplt.legend();","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"print(len(max_hu_scans), len(hu_scans))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Resampling the voxel size"},{"metadata":{},"cell_type":"markdown","source":"We can resize, crop, blur and shift intensities"},{"metadata":{"trusted":true},"cell_type":"code","source":"def resample(image, scan, new_spacing=[1,1,1]):\n    # Determine current pixel spacing\n    spacing = np.array([scan[0].SliceThickness] + list(scan[0].PixelSpacing), dtype=np.float32)\n    \n    resize_factor = spacing / new_spacing\n    new_shape = np.round(image.shape * resize_factor)\n    \n    # recompute the resize factor and spacing such that we match the rounded new shape above\n    rounded_resize_factor = new_shape / image.shape\n    rounded_new_spacing = spacing / rounded_resize_factor\n    \n    # zoom with resize factor\n    image = scipy.ndimage.interpolation.zoom(image, rounded_resize_factor, mode='nearest')\n    \n    return image, rounded_new_spacing","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"img_resampled, spacing = resample(max_hu_scans, scans, [1,1,1])\nprint(\"Shape before resampling\\t\", max_hu_scans.shape)\nprint(\"Shape after resampling\\t\", img_resampled.shape)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plot_3d(img_resampled)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"image_sizes = scan_properties.groupby([\"rows\", \"columns\"]).size().sort_values(ascending=False)\nimage_sizes","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.figure(figsize=(8,8))\nfor n in counts.index.values:\n    for m in counts.columns.values:\n        plt.scatter(n, m, s=counts.loc[n,m], c=\"dodgerblue\", alpha=0.7)\nplt.xlabel(\"rows\")\nplt.ylabel(\"columns\")\nplt.title(\"Pixel area of ct-scan per patient\");\nplt.plot(np.arange(0,1400), '-.', c=\"purple\", label=\"squared\")\nplt.plot(888 * np.ones(1400), '-.', c=\"crimson\", label=\"888 rows\");\nplt.legend();","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"class ImageObserver:\n    \n    def __init__(self, scan_properties, batch_size):\n        self.scan_properties = scan_properties\n        self.batch_size = batch_size\n    \n    def select_group(self, group=(512,512)):\n        self.group = group\n        self.name = \"rows {}, columns {}\".format(group[0], group[1])\n        self.batch_shape = (self.batch_size, group[0], group[1])\n        self.selection = self.scan_properties[\n            (self.scan_properties[\"rows\"]==group[0]) & (self.scan_properties[\"columns\"]==group[1])\n        ].copy()\n        self.patient_pths = self.selection.patient_pth.unique()\n    \n    \n    def get_loader(self):\n        \n        idx=0\n        images = np.zeros(self.batch_shape)\n        \n        for path in self.patient_pths:\n            \n            scans = load_scans(path)\n            hu_scans = transform_to_hu(scans)\n            images[idx,:,:] = hu_scans[0]\n            \n            idx += 1\n            if idx == self.batch_shape[0]:\n                yield images\n                images = np.zeros(self.batch_shape)\n                idx = 0\n        if idx > 0:\n            yield images","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"my_choice = image_sizes.index.values[0]\nprint(my_choice)\nto_display = 4","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"observer = ImageObserver(scan_properties, to_display)\nobserver.select_group(my_choice)\nobserver_iterator = observer.get_loader()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"images = next(observer_iterator)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"fig, ax = plt.subplots(1,to_display,figsize=(20,5))\n\n\nfor m in range(to_display):\n    image = images[m]\n    ax[m].imshow(set_manual_window(image, -500, 1000), cmap=\"YlGnBu\")\n    ax[m].set_title(observer.name)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"scan_properties.shape","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"scan_properties.head(1)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def resize_scan(scan, new_shape):\n    # read slice as 32 bit signed integers\n    img = Image.fromarray(scan, mode=\"I\")\n    # do the resizing\n    img = img.resize(new_shape, resample=Image.LANCZOS)\n    # convert back to 16 bit integers\n    resized_scan = np.array(img, dtype=np.int16)\n    return resized_scan","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def crop_scan(scan):\n    img = Image.fromarray(scan, mode=\"I\")\n    \n    left = (scan.shape[0]-512)/2\n    right = (scan.shape[0]+512)/2\n    top = (scan.shape[1]-512)/2\n    bottom = (scan.shape[1]+512)/2\n\n    img = img.crop((left, top, right, bottom))\n    # convert back to 16 bit integers\n    cropped_scan = np.array(img, dtype=np.int16)\n    return cropped_scan","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def crop_and_resize(scan, new_shape):\n    img = Image.fromarray(scan, mode=\"I\")\n    \n    left = (scan.shape[0]-512)/2\n    right = (scan.shape[0]+512)/2\n    top = (scan.shape[1]-512)/2\n    bottom = (scan.shape[1]+512)/2\n    \n    img = img.crop((left, top, right, bottom))\n    img = img.resize(new_shape, resample=Image.LANCZOS)\n    \n    cropped_resized_scan = np.array(img, dtype=np.int16)\n    return cropped_resized_scan","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def preprocess_to_hu_scans(scan_properties, my_shape, output_dir):\n    \n    for i, patient in enumerate(tqdm(scan_properties.patient.values)):\n        pth = scan_properties.loc[scan_properties.patient==patient].patient_pth.values[0]\n        scans = load_scans(pth)\n        hu_scans = transform_to_hu(scans) \n        prepared_scans = np.zeros((hu_scans.shape[0], my_shape[0], my_shape[1]), dtype=np.int16)\n        \n        # if squared:\n        if hu_scans.shape[1] == hu_scans.shape[2]:\n            \n            # if size is as desired\n            if hu_scans.shape[1] == my_shape[0]:\n                continue\n            # else resize:\n            else:\n               # as we have not converted to jpeg to keep all information, we need to do a workaround\n                hu_scans = hu_scans.astype(np.int32)\n                for s in range(hu_scans.shape[0]): \n                    prepared_scans[s] = resize_scan(hu_scans[s,:,:], my_shape)\n\n        # if non-squared - do a center crop to 512, 512 and then resize to desired shape\n        else:\n            hu_scans = hu_scans.astype(np.int32)\n            for s in range(hu_scans.shape[0]):\n                # if desired shape is 512x512:\n                if my_shape[0]==512:\n                    prepared_scans[s] = crop_scan(hu_scans[s,:,:])\n                else:\n                    prepared_scans[s] = crop_and_resize(hu_scans[s,:,:], my_shape)\n                \n        # save the prepared scans of patient:\n        np.save(output_dir + \"/\" + patient + '_hu_scans', prepared_scans)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"generate_512_512 = False\n\nif generate_512_512:\n    output_dir = \"scans_512x512\"\n    mkdir(output_dir)\n    my_shape = (512, 512)\n    preprocess_to_hu_scans(scan_properties, my_shape, output_dir)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"generate_224_224 = False\n\nif generate_224_224:\n    output_dir = \"scans_224x224\"\n    mkdir(output_dir)\n    my_shape = (224, 224)\n    preprocess_to_hu_scans(scan_properties, my_shape, output_dir)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"generate_128_128 = False\n\nif generate_128_128:\n    output_dir = \"scans_128x128\"\n    mkdir(output_dir)\n    my_shape = (128, 128)\n    preprocess_to_hu_scans(scan_properties, my_shape, output_dir)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"generate_64_64 = False\n\nif generate_64_64:\n    output_dir = \"scans_64x64\"\n    mkdir(output_dir)\n    my_shape = (64, 64)\n    preprocess_to_hu_scans(scan_properties, my_shape, output_dir)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"import matplotlib.pyplot as plt\nimport numpy as np\nimport random\nimport os\nimport glob\nimport pandas as pd \nfrom tqdm import tqdm\nimport sys\nimport glob\nimport cv2\n\nimport pydicom\nfrom sklearn.utils import shuffle\n\nimport albumentations as A\n\nimport torch\nimport torch.nn as nn\nfrom torch import optim\nfrom torch.utils.data import DataLoader\nfrom torch.utils.data.dataset import Dataset\nfrom torch.optim.lr_scheduler import  ReduceLROnPlateau","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sys.path.append('../input/efficientnet-pytorch/EfficientNet-PyTorch-master')\nsys.path.append('../input/pretrainedmodels/pretrainedmodels-0.7.4/')\nsys.path.append('../input/segmentation-models-pytorch/')\nimport segmentation_models_pytorch as smp","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"root_path = '../input/osiclungmask100'\nimg_path = sorted(glob.glob(root_path+'/*/img.png'))\nmask_path = sorted(glob.glob(root_path+'/*/post_label.png'))\n\nimgpaths,maskpaths = shuffle(img_path,mask_path, random_state=0)\n\ntrain_images_path = imgpaths[:int(len(imgpaths)*0.8)]\ntrain_masks_path = maskpaths[:int(len(imgpaths)*0.8)]\nval_images_path = imgpaths[int(len(imgpaths)*0.8):]\nval_masks_path = maskpaths[int(len(maskpaths)*0.8):]\n\ntransform = A.Compose([\n    A.Rotate(p=0.2,limit=30),\n    A.HorizontalFlip(p=0.2),\n    A.OneOf([\n        A.GridDistortion(p=0.1,distort_limit=0.2),\n        A.ElasticTransform(sigma=10, alpha=1,  p=0.1)\n    ]),\n])\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"batch = 8\nlr = 0.0003\nwd = 5e-4\nepochs = 80\noutput_path = './'\ndevice =  torch.device('cuda:0')\nexperiment_name = 'lung_Unet_densenet121'","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"class Data_Generate(Dataset):\n    def __init__(self,img_paths,seg_paths=None,transform=None):\n        self.img_paths = img_paths\n        self.seg_paths = seg_paths\n        self.transform = transform\n        \n    def __getitem__(self,index):\n        if self.seg_paths is not None:\n            img_path = self.img_paths[index]\n            mask_path = self.seg_paths[index]\n            \n            mask = cv2.imread(mask_path,0)/255\n            img = cv2.imread(img_path,0)/255\n\n            if self.transform != None:\n                aug = transform(image=img,mask=mask)\n                img = aug['image']\n                mask = aug['mask']\n                \n            img = img[None,:,:]\n            img = img.astype(np.float32)\n            mask = mask[None,:,:]\n            mask = mask.astype(np.float32)\n            \n            return img,mask\n        \n        else:\n            img = cv2.imread(self.img_paths[index],0)/255\n            img = img[None,:,:]\n            img = img.astype(np.float32)\n            return img\n        \n    def __len__(self):\n        return len(self.img_paths)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_db = Data_Generate(train_images_path,train_masks_path,transform=transform)\ntrain_loader = DataLoader(train_db, batch_size=batch, shuffle=True, num_workers=4)\nval_db = Data_Generate(val_images_path,val_masks_path,transform=None)\nval_loader = DataLoader(val_db, batch_size=batch, shuffle=False, num_workers=4)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"f,ax = plt.subplots(4,4,figsize=(16,16))\nfor i in range(16):\n    img = train_db[i][0]\n    ax[i//4,i%4].imshow(img[0])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"class EarlyStopping:\n    \"\"\"Early stops the training if validation loss doesn't improve after a given patience.\"\"\"\n    def __init__(self, patience=7, verbose=False, delta=0, path='checkpoint.pt'):\n        \"\"\"\n        Args:\n            patience (int): How long to wait after last time validation loss improved.\n                            Default: 7\n            verbose (bool): If True, prints a message for each validation loss improvement. \n                            Default: False\n            delta (float): Minimum change in the monitored quantity to qualify as an improvement.\n                            Default: 0\n            path (str): Path for the checkpoint to be saved to.\n                            Default: 'checkpoint.pt'\n        \"\"\"\n        self.patience = patience\n        self.verbose = verbose\n        self.counter = 0\n        self.best_score = None\n        self.early_stop = False\n        self.val_loss_min = np.Inf\n        self.delta = delta\n        self.path = path\n\n    def __call__(self, val_loss, model):\n\n        score = -val_loss\n\n        if self.best_score is None:\n            self.best_score = score\n            self.save_checkpoint(val_loss, model)\n        elif score < self.best_score + self.delta:\n            self.counter += 1\n            print(f'EarlyStopping counter: {self.counter} out of {self.patience}')\n            if self.counter >= self.patience:\n                self.early_stop = True\n        else:\n            self.best_score = score\n            self.save_checkpoint(val_loss, model)\n            self.counter = 0\n\n    def save_checkpoint(self, val_loss, model):\n        '''Saves model when validation loss decrease.'''\n        if self.verbose:\n            print(f'Validation loss decreased ({self.val_loss_min:.6f} --> {val_loss:.6f}).  Saving model ...')\n        torch.save(model.state_dict(), self.path)\n        self.val_loss_min = val_loss","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"model = smp.Unet('densenet121', classes=1, in_channels=1,activation='sigmoid',encoder_weights='imagenet').to(device)\n    \noptimizer = optim.Adam(filter(lambda p: p.requires_grad, model.parameters()), lr=lr, weight_decay=wd)\nscheduler = ReduceLROnPlateau(optimizer, mode='min', factor=0.5, patience=3, min_lr=1e-8, verbose=True)\n\ncriterion = smp.utils.losses.DiceLoss(eps=1.)\niou = smp.utils.metrics.IoU()\nearly_stopping = EarlyStopping(patience=6, verbose=True,path=os.path.join(output_path, f'best_{experiment_name}.pth'))\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"num_train_loader = len(train_loader)\nnum_val_loader = len(val_loader)\nfor epoch in range(epochs):\n    train_losses,train_score,val_losses,val_score = 0,0,0,0\n    model.train()\n\n    for idx, sample in enumerate(train_loader):\n        image, label = sample\n        image, label = image.to(device), label.to(device)\n        out = model(image)\n        loss = criterion(out, label)\n        optimizer.zero_grad()\n        loss.backward()\n        optimizer.step()\n        train_losses += loss/num_train_loader\n        train_score += iou(out,label)/num_train_loader\n    \n    model.eval()\n    for idx, sample in enumerate(val_loader):\n        image, label = sample\n        image, label = image.to(device), label.to(device)\n        with torch.no_grad():\n            out = model(image)\n        loss = criterion(out, label)\n        val_losses += loss/num_val_loader\n        val_score += iou(out,label)/num_val_loader\n    print('epoch {}/{}\\t LR:{}\\t train_loss:{}\\t train_score:{}\\t val_loss:{}\\t val_score:{}' \\\n          .format(epoch+1, epochs, optimizer.param_groups[0]['lr'], train_losses, train_score, val_losses, val_score))\n    scheduler.step(val_losses)\n    \n    early_stopping(val_losses, model)\n    if early_stopping.early_stop:\n        print(\"Early stopping\")\n        break","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"class Test_Generate(Dataset):\n    def __init__(self,img_paths):\n        self.img_paths = img_paths\n        \n    def __getitem__(self,index):\n        dicom = pydicom.dcmread(self.img_paths[index])\n        slice_img = dicom.pixel_array\n        slice_img = (slice_img-slice_img.min())/(slice_img.max()-slice_img.min())\n        slice_img = (slice_img*255).astype(np.uint8)\n        if slice_img.shape[0] != 512:\n            slice_img = cv2.resize(slice_img,(512,512))\n            \n        slice_img = slice_img[None,:,:]\n        slice_img = (slice_img/255).astype(np.float32)\n        return slice_img\n        \n    def __len__(self):\n        return len(self.img_paths)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"dicom_root_path = '../input/osic-pulmonary-fibrosis-progression/train/*/*'\ndicom_paths = glob.glob(dicom_root_path)\ndicom_paths = random.sample(dicom_paths,16)\n\ntest_db = Test_Generate(dicom_paths)\ntest_loader = DataLoader(test_db, batch_size=batch, shuffle=False, num_workers=0)\n\nmodel.load_state_dict(torch.load('./best_lung_Unet_densenet121.pth'))\nmodel.eval()\n\nouts = []\nfor idx, sample in enumerate(test_loader):\n    image = sample\n    image = image.to(device)\n    with torch.no_grad():\n        out = model(image)\n    out = out.cpu().data.numpy()\n    out = np.where(out>0.5,1,0)\n    out = np.squeeze(out)\n    outs.append(out)\n    \nouts = np.concatenate(outs)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"f,ax = plt.subplots(4,4,figsize=(16,16))\naxes = ax.flatten()\nfor idx in range(len(outs)//2):\n    axes[idx*2].imshow(test_db[idx][0])\n    axes[idx*2+1].imshow(outs[idx])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"","execution_count":null,"outputs":[]}],"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}