{"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\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":"%matplotlib inline\n\nimport numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\nimport pydicom\nimport os\nimport scipy.ndimage\nimport matplotlib.pyplot as plt\n\nfrom skimage import measure, morphology\nfrom mpl_toolkits.mplot3d.art3d import Poly3DCollection\n\n# Some constants \nINPUT_FOLDER = '../input/osic-pulmonary-fibrosis-progression/train/'\npatients = os.listdir(INPUT_FOLDER)\npatients.sort()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def load_scan(path):\n    slices = [pydicom.dcmread(path + '/' + s) for s in os.listdir(path)]\n    slices.sort(key = lambda x: float(x.ImagePositionPatient[2]))\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":{"trusted":true},"cell_type":"code","source":"def get_pixels_hu(slices):\n    image = np.stack([s.pixel_array for s in slices])\n    print((slices[10].pixel_array).shape)\n    # Convert to int16 (from sometimes int16), \n    # should be possible as values should always be low enough (<32k)\n    image = image.astype(np.int16)\n\n    # Set outside-of-scan pixels to 0\n    # The intercept is usually -1024, so air is approximately 0\n    image[image == -2000] = 0\n    \n    # Convert to Hounsfield units (HU)\n    for slice_number in range(len(slices)):\n        \n        intercept = slices[slice_number].RescaleIntercept\n        slope = slices[slice_number].RescaleSlope\n        \n        if slope != 1:\n            image[slice_number] = slope * image[slice_number].astype(np.float64)\n            image[slice_number] = image[slice_number].astype(np.int16)\n            \n        image[slice_number] += np.int16(intercept)\n    \n    return np.array(image, dtype=np.int16)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"first_patient = load_scan(INPUT_FOLDER + patients[0])\nprint(type(first_patient))\nfirst_patient_pixels = get_pixels_hu(first_patient)\n# plt.hist(first_patient_pixels.flatten(), bins=80, color='c')\n# plt.xlabel(\"Hounsfield Units (HU)\")\n# plt.ylabel(\"Frequency\")\n# plt.show()\n\n# Show some slice in the middle\n# for i in range(len(patients[0])):\n#     print(i,'\\n')\n#     plt.imshow(first_patient_pixels[i], cmap=plt.cm.gray)\n#     plt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.imshow(first_patient_pixels[10], cmap=plt.cm.gray)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"s_patient = load_scan(INPUT_FOLDER + patients[2])\ns_patient_pixels = get_pixels_hu(s_patient)\nplt.hist(s_patient_pixels.flatten(), bins=80, color='c')\nplt.xlabel(\"Hounsfield Units (HU)\")\nplt.ylabel(\"Frequency\")\nplt.show()\n\n# Show some slice in the middle\nfor i in range(len(patients[2])):\n    plt.imshow(s_patient_pixels[i], cmap=plt.cm.gray)\n    plt.show()","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] + scan[0].PixelSpacing, dtype=np.float32)\n    spacing = map(float, ([scan[0].SliceThickness] + list(scan[0].PixelSpacing)))\n    spacing = np.array(list(spacing))\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":"print(300.0/334)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# import shutil\n# for file in os.listdir('./'):\n#     if os.path.isdir(file):\n#         shutil.rmtree(file)\n#     else:\n#         os.remove(file)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"import datetime\ndef get_now_string():\n  x = str(datetime.datetime.now())\n  date=x[:10].split('-')\n  time=x[11:].split(':')\n  time[2]=time[2].split('.')[0]\n  dates=\"\".join(date)+\"_\"+\"\".join(time)\n  return dates\n\nOUTPUT_DIR=get_now_string()+'/'\nos.mkdir(OUTPUT_DIR)\nprint(\"saving models to :\",OUTPUT_DIR)\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"np.savez(OUTPUT_DIR+'data', s_patient_pixels)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"from scipy.ndimage import zoom\n\npix_resampled, spacing = resample(first_patient_pixels, first_patient, [1,1,1])\na1 =  300.0/pix_resampled.shape[0]\na2 =  300.0/pix_resampled.shape[1]\na3 =  300.0/pix_resampled.shape[2]\nnew_size=(300,300,300)\n\nnew_array = zoom(pix_resampled, (a1, a2, a3))\nprint(\"Shape before resampling\\t\", new_array.shape)\nprint(\"Shape before resampling\\t\", first_patient_pixels.shape)\nprint(\"Shape after resampling\\t\", pix_resampled.shape)\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def plot_3d(image, threshold=-300):\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_classic(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.70)\n    face_color = [0.45, 0.45, 0.75]\n    mesh.set_facecolor(face_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(first_patient_pixels, 400)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"print(np.array(pix_resampled[15][250]>-320,dtype=np.int8))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def 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    # not actually binary, but 1 and 2. \n    # 0 is treated as background, which we do not want\n    binary_image = np.array(image > -500, 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 = labels[0,0,0]\n    \n    #Fill the air around the person\n    binary_image[background_label == labels] = 2\n    \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":"segmented_lungs = segment_lung_mask(pix_resampled, False)\nsegmented_lungs_fill = segment_lung_mask(pix_resampled, True)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# plot_3d(segmented_lungs, 0)","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}