{"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)\nimport matplotlib.pyplot as plt\nimport seaborn as sns\n\nimport pydicom\nimport scipy.ndimage\n\nfrom skimage import measure \nfrom mpl_toolkits.mplot3d.art3d import Poly3DCollection\n\nfrom IPython.display import HTML\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\nfor dirname, _, filenames in os.walk('/kaggle/input'):\n    for filename in filenames:\n        print(os.path.join(dirname, filename))\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":{"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","trusted":true},"cell_type":"code","source":"from os import listdir\n\nbasepath = \"../input/osic-pulmonary-fibrosis-progression/\"\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\ndef load_scans(dcm_path):\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":"print( train)\nprint( train.Patient.values)\nprint( train.Patient.values[0])\n#../input/osic-pulmonary-fibrosis-progression/train\n#../input/osic-pulmonary-fibrosis-progression/train/ID00007637202177411956430/1.dcm\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"example = basepath + \"train/\" + train.Patient.values[0]\nscans = load_scans(example)\n\nscans[0]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"print(scans[0].PatientID)\n#医療機器製造会社\nprint(scans[0].Manufacturer)\n#スライス厚さ\nprint(scans[0][0x0018,0x0050],\"mm\")\n\nct_image = scans[0].pixel_array\nplt.imshow(ct_image)\n\n\n#../input/osic-pulmonary-fibrosis-progression/train/ID00007637202177411956430/10.dcm","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"import math\nax=len(scans)\nx=int(math.sqrt(ax))\ny=x+1\n\nprint(x,y,ax)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"\n\"\"\"\nax=len(scans)\nx=int(math.sqrt(ax))\ny=x+1\n\n\nfig=plt.figure(figsize=(12, 12))\nfor i in range(ax):\n    ct_image = scans[i].pixel_array\n    fig.add_subplot(x,y,i)\n\n    plt.imshow(scans[i].pixel_array, cmap='gray')#グレースケールで出力\n    \nplt.show()\n\"\"\"\n","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    # convert ouside pixel-values to air:\n    # I'm using <= -1000 to be sure that other defaults are captured as well\n    images[images <= -1000] = 0\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)\n","execution_count":null,"outputs":[]},{"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":"","execution_count":null,"outputs":[]},{"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\nfor patient in train.Patient.values:\n    patient_id.append(patient)\n    example_dcm = listdir(basepath + \"train/\" + patient + \"/\")[0]\n    patient_pth.append(basepath + \"train/\" + patient)\n    dataset = pydicom.dcmread(basepath + \"train/\" + patient + \"/\" + 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\n\n\nscan_properties.head()\n","execution_count":null,"outputs":[]},{"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":{"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":{"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\n\nfig, 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\n\ndef 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\n\n\nfig, ax = plt.subplots(1,2,figsize=(20,10))\nax[0].imshow(set_manual_window(min_hu_scans[np.int(len(min_hu_scans)/2)], -500, 1000), cmap=\"YlGnBu\")\nax[1].imshow(set_manual_window(max_hu_scans[np.int(len(max_hu_scans)/2)], -500, 1000), 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":{"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)\n\n\nfig, ax = plt.subplots(1,2,figsize=(20,10))\nax[0].imshow(set_manual_window(min_hu_scans[np.int(len(min_hu_scans)/2)], -500, 1000), cmap=\"YlGnBu\")\nax[1].imshow(set_manual_window(max_hu_scans[np.int(len(max_hu_scans)/2)], -500, 1000), 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\");\n","execution_count":null,"outputs":[]},{"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.5)\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()\n    \n    \nplot_3d(max_hu_scans)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"old_distribution = max_hu_scans.flatten()\n\nexample = basepath + \"train/\" + train.Patient.values[0]\nscans = load_scans(example)\nhu_scans = transform_to_hu(scans)\n\nplot_3d(hu_scans)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.figure(figsize=(20,5))\nsns.distplot(old_distribution, label=\"weak 3d plot\")\nsns.distplot(hu_scans.flatten(), label=\"strong 3d plot\")\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":{"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_real_shape = image.shape * resize_factor\n    new_shape = np.round(new_real_shape)\n    real_resize_factor = new_shape / image.shape\n    new_spacing = spacing / real_resize_factor\n    \n    image = scipy.ndimage.interpolation.zoom(image, real_resize_factor, mode='nearest')\n    \n    return image, new_spacing","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"pix_resampled, spacing = resample(max_hu_scans, scans, [1,1,1])\nprint(\"Shape before resampling\\t\", hu_scans.shape)\nprint(\"Shape after resampling\\t\", pix_resampled.shape)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Guidos code:\n\ndef largest_label_volume(im, bg=-1):\n    vals, counts = np.unique(im, return_counts=True)\n\n    counts = counts[vals != bg]\n    vals = vals[vals != bg]\n\n    if len(counts) > 0:\n        return vals[np.argmax(counts)]\n    else:\n        return None\n    \ndef segment_lung_mask(image, fill_lung_structures=True):\n    \n    # not actually binary, but 1 and 2. \n    # 0 is treated as background, which we do not want\n    binary_image = np.array(image > -320, dtype=np.int8)+1\n    labels = measure.label(binary_image)\n    \n    # Pick the pixel in the very corner to determine which label is air.\n    #   Improvement: Pick multiple background labels from around the patient\n    #   More resistant to \"trays\" on which the patient lays cutting the air \n    #   around the person in half\n    background_label_1 = labels[0,0]\n    background_label_2 = labels[0,-1]\n    background_label_3 = labels[-1,0]\n    background_label_4 = labels[-1,-1]\n    \n    #Fill the air around the person\n    binary_image[background_label_1 == labels] = 2\n    binary_image[background_label_2 == labels] = 2\n    binary_image[background_label_3 == labels] = 2\n    binary_image[background_label_4 == labels] = 2\n    \n    # Method of filling the lung structures (that is superior to something like \n    # morphological closing)\n    if fill_lung_structures:\n        # For every slice we determine the largest solid structure\n        for i, axial_slice in enumerate(binary_image):\n            axial_slice = axial_slice - 1\n            labeling = measure.label(axial_slice)\n            l_max = largest_label_volume(labeling, bg=0)\n            \n            if l_max is not None: #This slice contains some lung\n                binary_image[i][labeling != l_max] = 1\n\n    \n    binary_image -= 1 #Make the image actual binary\n    binary_image = 1-binary_image # Invert it, lungs are now 1\n    \n    # Remove other air pockets insided body\n    labels = measure.label(binary_image, background=0)\n    l_max = largest_label_volume(labels, bg=0)\n    if l_max is not None: # There are air pockets\n        binary_image[labels != l_max] = 0\n \n    return binary_image","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"binary_image = np.array((hu_scans[20]>-320), dtype=np.int8) + 1\nnp.unique(binary_image)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"labels = measure.label(binary_image)\n\nbackground_label_1 = labels[0,0]\nbackground_label_2 = labels[0,-1]\nbackground_label_3 = labels[-1,0]\nbackground_label_4 = labels[-1,-1]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"binary_image_2 = binary_image.copy()\nbinary_image_2[background_label_1 == labels] = 2\nbinary_image_2[background_label_2 == labels] = 2\nbinary_image_2[background_label_3 == labels] = 2\nbinary_image_2[background_label_4 == labels] = 2","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"fig, ax = plt.subplots(1,4,figsize=(20,5))\nax[0].imshow(binary_image, cmap=\"binary\", interpolation='nearest')\nax[1].imshow(labels, cmap=\"jet\", interpolation='nearest')\nax[2].imshow(binary_image_2, cmap=\"binary\", interpolation='nearest')\nax[0].set_title(\"Binary image\")\nax[1].set_title(\"Labelled image\");","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"segmented_lungs = segment_lung_mask(hu_scans, fill_lung_structures=False)\n#segmented_lungs_fill = segment_lung_mask(hu_scans, fill_lung_structures=True)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.imshow(segmented_lungs[20])","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}