{"cells":[{"metadata":{"_kg_hide-input":true,"_kg_hide-output":true,"trusted":true},"cell_type":"code","source":"!conda install -c conda-forge gdcm -y","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true,"_kg_hide-input":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":{},"cell_type":"markdown","source":"## Table of contents\n\n\n1. [References](#references)\n    * [Data Science Bowl 2017 - Preprocessing Tutorial by Guido Zuidhof](#bowl_2017)\n    * [Papers](#papers)\n2. [Prepare to start](#prepare)\n3. [Working with dicom files](#dicom)\n    * [Loading CT-scans per patient](#ct_scans)\n    * [Transforming to Hounsfield Units](#hunits)\n    * [The voxel size](#voxel)\n    * [CT-scan slice area and volume - EDA](#scan_eda)\n    * [3D-reconstruction of CT-scans](#reconstruction)\n    * [Tissue segmentation](#segmentation)\n4. [Generating a dataset for preprocessed files](#datagenerator)\n    * [Link to public dataset](https://www.kaggle.com/allunia/osic-pulmonary-fibrosis-progression-huscans)"},{"metadata":{},"cell_type":"markdown","source":"# References <a class=\"anchor\" id=\"references\"></a>"},{"metadata":{},"cell_type":"markdown","source":"## Data Science Bowl 2017 - Preprocessing Tutorial by Guido Zuidhof <a class=\"anchor\" id=\"bowl_2017\"></a>\n\nOnce upon a time there was a data science bowl about lung cancer detection that had a fantastic preprocessing tutorial: https://www.kaggle.com/gzuidhof/full-preprocessing-tutorial ;-)\n\n* My notebook heavily uses the concepts and code snippets of this notebook by Guido Zuidhof. \n* With my work I like to bring back his ideas to this competition and I like to add a few more explanations while writing it. \n* Furthermore I hope that I can generate a dataset in the end that contains the preprocessed data to feed into your models."},{"metadata":{},"cell_type":"markdown","source":"## Papers <a class=\"anchor\" id=\"papers\"></a>\n\nTo understand what and why we are doing these concepts here, I started to read a few papers:\n\n* [Intrinsic dependencies of CT radiomic features on voxel size and number of gray levels](https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5462462/)"},{"metadata":{},"cell_type":"markdown","source":"# Prepare to start <a class=\"anchor\" id=\"prepare\"></a>"},{"metadata":{"trusted":true},"cell_type":"code","source":"listdir(\"../input/\")","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":{},"cell_type":"markdown","source":"Let's load the csv-files:"},{"metadata":{"trusted":true},"cell_type":"code","source":"train = pd.read_csv(basepath + \"train.csv\")\ntest = pd.read_csv(basepath + \"test.csv\")","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":"# Working with dicom files <a class=\"anchor\" id=\"dicom\"></a>"},{"metadata":{},"cell_type":"markdown","source":"## Loading CT-scans per patient <a class=\"anchor\" id=\"ct_scans\"></a>\n\n* To load the full 3D-scan we need to order the single dicom files/slices by the ImagePosition: "},{"metadata":{"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","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":"Let's take a look at the first dicom file of our example patient:"},{"metadata":{"trusted":true},"cell_type":"code","source":"scans[0]","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"If you haven't worked with dicom so far, I can recommend this video. If you like to speed up, start at 7 min:"},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"HTML('<iframe width=\"600\" height=\"400\" src=\"https://www.youtube.com/embed/KZld-5W99cI\" frameborder=\"0\" allow=\"accelerometer; autoplay; encrypted-media; gyroscope; picture-in-picture\" allowfullscreen></iframe>')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Adding further information\n\n1. The CT-scan captures information about the radiodensity of an object or tissue exposed to x-rays. A transversal slice of a scan is reconstructed after taking measurements from several different directions.\n2. 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.\n3. A ct-scanner yields roughly 4000 grey values that can't be captured by our eyes. This is why windowing is performed. This way the image is displayed in a HU range that suites most to the region of interest. "},{"metadata":{},"cell_type":"markdown","source":"## Transforming to Hounsfield Units <a class=\"anchor\" id=\"hunits\"></a>"},{"metadata":{},"cell_type":"markdown","source":"Before starting, let's plot the pixelarray distribution of some dicom files to get an impression of the raw data:"},{"metadata":{"_kg_hide-input":true,"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 (or in other competitions I found also -3000) by default."},{"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":{"trusted":true},"cell_type":"code","source":"hu_scans = transform_to_hu(scans)","execution_count":null,"outputs":[]},{"metadata":{"_kg_hide-input":true,"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":{},"cell_type":"markdown","source":"Ok great. Now all raw values per slice are scaled to H-units."},{"metadata":{},"cell_type":"markdown","source":"## The voxel size <a class=\"anchor\" id=\"voxel\"></a>\n\nThe voxel stands for the 3D-pixel that is given in a CT-scan. As far as I know it is spanned by the 2d-plane of the pixelspacing attribute in x- and y-direction and the slice thickness in z-direction."},{"metadata":{},"cell_type":"markdown","source":"### Pixelspacing\n\n* The pixelspacing attribute you can find in the dicom files is an important one. It tells us how much physical distance is covered by one pixel. You can see that there are only 2 values that describe the x- and y-direction in the plane of a transversal slice. \n* For one patient this pixelspacing is usually the same for all slices.\n* But 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!\n\nLet's explore the distributions of the patients pixelspacing widths and heights of this competition:"},{"metadata":{"trusted":true},"cell_type":"code","source":"N = 1000","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Just to speed up the computation, I have selected N patients to consider. Use N = train.shape[0] to do this for all patients in the dataset:"},{"metadata":{"_kg_hide-input":true,"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":{"_kg_hide-input":true,"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":"We can see that the values really vary a lot from patient to patient! As they are given in mm and ct-scans usually cover 512 row and column values... **oh wait! we need to check this!** ... we can compute the minimum and maximum distance that is covered by the images:"},{"metadata":{},"cell_type":"markdown","source":"### Slice thickness and pixel area\n\nThe slice thickness tells us how much distance is covered in Z-direction by one slice. Let's plot the distribution of it as well. Furthermore the pixel_array of raw values covers a specific area given by row and column values. Let's take a look at it as well: "},{"metadata":{"_kg_hide-input":true,"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":"* Very thin slices allow more details to be shown. On the other hand thick slices contain less noise but are more prone to artifacts. Hmm... I'm very excited to see some examples here as well. \n* Even though it is common to have 512x512 pixel size areas, we can see that this is not always true!! We can find a lot of exceptions and even one or a few very large pixel areas (1300x1300)!!! (OSIC)\n* Proper preprocessing of these scans might be very important... we have to check it."},{"metadata":{},"cell_type":"markdown","source":"## Physical area & slice volume covered by a single ct-scan"},{"metadata":{"trusted":true},"cell_type":"markdown","source":"Now, we know some important quantities to compute the physical distance covered by a ct-scan!"},{"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":{"trusted":true,"_kg_hide-input":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":{},"cell_type":"markdown","source":"Ohoh! We have some images with extreme large sliche areas and volumes! I think it's time to do some EDA regarding these features! "},{"metadata":{},"cell_type":"markdown","source":"## CT-scan slice area and volume - EDA <a class=\"anchor\" id=\"scan_eda\"></a>"},{"metadata":{},"cell_type":"markdown","source":"### Smallest and larges CT-slice area"},{"metadata":{"trusted":true,"_kg_hide-input":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":{"trusted":true,"_kg_hide-input":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":{},"cell_type":"markdown","source":"### Insights\n\n* Taking a look at one slice of a scan with smallest and largest slice area, we can see that the large one has a lot of useless region covered. We could crop it.\n* Strange... in the second image with the large area the outside region of the scanner tube is not set to the value of air but rather to some value in the middle of the range -1000 to 1000."},{"metadata":{"_kg_hide-input":true,"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":"* For the OSIC example we can find: Ihh... it was set to water by default in the large image... why?! That's bad! We need to find some strategy to deal with this problem. It's not good that we sometimes have \"water\"-like outside regions and sometimes \"air\"-like regions."},{"metadata":{},"cell_type":"markdown","source":"### Smallest and largest CT-slice volume"},{"metadata":{"trusted":true,"_kg_hide-input":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,"_kg_hide-input":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":{},"cell_type":"markdown","source":"Hmm I can't see a great difference. Perhaps the one with the large slice volume looks a bit more blurred. But as above there is an outside-scanner region that was set to the walue of water (HU value of 0) instead of air."},{"metadata":{"_kg_hide-input":true,"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 <a class=\"anchor\" id=\"reconstruction\"></a>"},{"metadata":{},"cell_type":"markdown","source":"The plot_3d works well in the Data Science Bowl 2017 but in our case the results are not so well. It depends on the threshold but so far I don't know why our reconstructions often look blurred or show tube regions as well. :-("},{"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":{},"cell_type":"markdown","source":"You can see that this looks quite different to this one:"},{"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":"Compare to above this one looks far better. Let's plot the distributions. Perhaps we can understand what's going wrong by taking a look at them:"},{"metadata":{"_kg_hide-input":true,"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":{},"cell_type":"markdown","source":"Hmm"},{"metadata":{"trusted":true},"cell_type":"code","source":"print(len(max_hu_scans), len(hu_scans))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"I think we need to understand the marching_cubes_lewiner algorithm to understand why the plot sometimes works nice and sometimes not. But I think this is not really important for the competition itself. At the moment I don't like to spend more time on that topic. Perhaps it's more important to keep in mind that the overall distributions can be different."},{"metadata":{},"cell_type":"markdown","source":"## Resampling the voxel size"},{"metadata":{},"cell_type":"markdown","source":"I'm not sure if it's really important to do this resampling, because we can already do a lot with augmentations. We can resize, crop, blur and shift intensities. Also the background-problem with outside scanner regions set to water HU-values can be solved. For this reason I don't like to do the resampling. ;-) Perhaps you like to do it. For this purpose I leave the code from Guidos notebook: "},{"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":{},"cell_type":"markdown","source":"Looks the same as the max_hu_scans above!"},{"metadata":{},"cell_type":"markdown","source":"## Tissue segmentation <a class=\"anchor\" id=\"segmentation\"></a>"},{"metadata":{},"cell_type":"markdown","source":"I built this solution upon Guidos that did not work well in this competition. Here is code I won't use but I kept it to show by example why I found it not helpful:"},{"metadata":{"trusted":true,"_kg_hide-input":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 fill_lungs(binary_image):\n    image = binary_image.copy()\n    # For every slice we determine the largest solid structure\n    for i, axial_slice in enumerate(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            image[i][labeling != l_max] = 1\n    return image\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"This is now the segmentation function I am using:"},{"metadata":{"trusted":true},"cell_type":"code","source":"def segment_lung_mask(image):\n    segmented = np.zeros(image.shape)   \n    \n    for n in range(image.shape[0]):\n        binary_image = np.array(image[n] > -320, dtype=np.int8)+1\n        labels = measure.label(binary_image)\n        \n        bad_labels = np.unique([labels[0,:], labels[-1,:], labels[:,0], labels[:,-1]])\n        for bad_label in bad_labels:\n            binary_image[labels == bad_label] = 2\n    \n        #We have a lot of remaining small signals outside of the lungs that need to be removed. \n        #In our competition closing is superior to fill_lungs \n        selem = disk(2)\n        binary_image = opening(binary_image, selem)\n    \n        binary_image -= 1 #Make the image actual binary\n        binary_image = 1-binary_image # Invert it, lungs are now 1\n        \n        segmented[n] = binary_image.copy() * image[n]\n    \n    return segmented","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Understanding the segmentation step by step:"},{"metadata":{},"cell_type":"markdown","source":"First of all we are separating between \"potentially lung\" with values smaller or equal to 320. Why 320? I just kept the value of Guidos old notebook. Try out your own values if you like. ;-) Before you do so, you make like to take a look at the hu_value distribution again:"},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"plt.figure(figsize=(20,5))\nsns.distplot(hu_scans[20], kde=False)\nplt.title(\"Example HU value distribution\");\nplt.xlabel(\"HU-value\")\nplt.ylabel(\"count\")","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"With -320 we are separating between lungs (-700) /air (-1000) and tissue with values close to water (0)."},{"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":{},"cell_type":"markdown","source":"Lungs have values of 1 as well as air backgrounds. In contrast water-like default backgrounds and many other organic tissues or fluids have values of 2. As we only like to segment the lungs, we need to remove the background. In the case of air and air-like default values we need to do set their values manually to 2. We can do so by labelling connected regions in the binary image and extracting the labels of each region that corresponds to the corners or edges of the 2D image slice:"},{"metadata":{"trusted":true},"cell_type":"code","source":"labels = measure.label(binary_image)\n\nbad_labels = np.unique([labels[0,:], labels[-1,:], labels[:,0], labels[:,-1]])","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Now we can set all labelled regions of the binary image that correspond to these corner-labels to the \"not-lung\"-value 2:"},{"metadata":{"trusted":true},"cell_type":"code","source":"binary_image_2 = binary_image.copy()\nfor bad_label in bad_labels:\n    binary_image_2[labels == bad_label] = 2","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The result of these steps looks as follows:"},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"fig, ax = plt.subplots(1,3,figsize=(20,7))\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')\n\nax[0].set_title(\"Binary image\")\nax[1].set_title(\"Labelled image\");\nax[2].set_title(\"Binary image - background removed\");","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Insights\n\n1. The first image shows the raw binary. In this case we find air as background and we need to set it to the \"not-lung\" value of 2.\n2. For this purpose we label all connected regions in the binary image. There are many different labelled regions but the only ones we are interested in are the 4 corner regions at (0,0), (0,500), (500,0) and (500,500).\n3. Knowing the related labels help us to manually set the background to the value 2 (black).\n4. In the end we can see that the lungs are white (1) but we still find a lot of remaining signals that correspond to body tissue that we still need to remove."},{"metadata":{},"cell_type":"markdown","source":"To remove the signals we can use morphological closing. If you like play with the disk value. It must be large enough to cancel the body signals but small enough to keep enough details inside the lungs:"},{"metadata":{"trusted":true},"cell_type":"code","source":"selem = disk(2)\nclosed_binary_2 = closing(binary_image_2, selem)\n\nclosed_binary_2 -= 1 #Make the image actual binary\nclosed_binary_2 = 1-closed_binary_2 # Invert it, lungs are now 1","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Let's compare with Guidos fill-lungs method:"},{"metadata":{"trusted":true},"cell_type":"code","source":"filled_lungs_binary = fill_lungs(binary_image_2)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"And with his air-pocket removal:"},{"metadata":{"trusted":true},"cell_type":"code","source":"air_pocket_binary = closed_binary_2.copy()\n# Remove other air pockets insided body\nlabels_2 = measure.label(air_pocket_binary, background=0)\nl_max = largest_label_volume(labels_2, bg=0)\nif l_max is not None: # There are air pockets\n    air_pocket_binary[labels_2 != l_max] = 0","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"fig, ax = plt.subplots(1,3,figsize=(20,7))\n\nax[0].imshow(closed_binary_2, cmap=\"binary\", interpolation='nearest')\nax[1].imshow(filled_lungs_binary, cmap=\"binary\", interpolation='nearest')\nax[2].imshow(air_pocket_binary, cmap=\"binary\", interpolation='nearest')\n\n\nax[0].set_title(\"Morphological closing\");\nax[1].set_title(\"Guidos filling lung structures\");\nax[2].set_title(\"Guidos air pocket removal\");","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Insights\n\n* The morphological closing worked better than Guidos method. But we are often missing a lot of information inside the lungs just to remove these body signals. There is definitely room for improvements! ;-)\n* In contrast the fill-lung method has troubles with the remaining singals of the body and yields not what we like to obtain.\n* Furthermore the air-pocket removal does also not work well.\n\nProbably I'm having some kind of bug. Can't imagine that Guidos methods didn't work at all?! For this reason I just stayed with the simplest solution of morphological closing."},{"metadata":{},"cell_type":"markdown","source":"### Final solution"},{"metadata":{},"cell_type":"markdown","source":"And here is how it looks like if we mask the original image with the segmented binary lungs in our 2D example:"},{"metadata":{"trusted":true},"cell_type":"code","source":"segmented = segment_lung_mask(np.array([hu_scans[20]]))\n\nfig, ax = plt.subplots(1,2,figsize=(20,10))\nax[0].imshow(set_manual_window(hu_scans[20], -700, 255), cmap=\"Blues_r\")\nax[1].imshow(set_manual_window(segmented[0], -700, 255), cmap=\"Blues_r\");","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"And we can also check how it looks in the 3D case:"},{"metadata":{"trusted":true},"cell_type":"code","source":"segmented_lungs = segment_lung_mask(hu_scans)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"fig, ax = plt.subplots(6,5, figsize=(20,20))\nfor n in range(6):\n    for m in range(5):\n        ax[n,m].imshow(set_manual_window(segmented_lungs[(n+1)*5+m], -700, 255), cmap=\"Blues_r\")","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Interesting that we still have some signals outside of the lungs:"},{"metadata":{"trusted":true},"cell_type":"code","source":"plot_3d(segmented_lungs, threshold=-600)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"If you like try to improve the segmentation and play with your own ideas and concepts! :-) I like to create two datasets - one with segmentation already done and one without. So you can decide on your own if you like to build a custom augmentation that performs the segmentation on the fly when loading your batches."},{"metadata":{},"cell_type":"markdown","source":"# Generating a dataset for preprocessed files <a class=\"anchor\" id=\"datagenerator\"></a>\n\n\n## Dealing with different image-scan sizes <a class=\"anchor\" id=\"image_sizes\"></a>\n\nTo generate the data, we should take a look again at the different images sizes: For OSIC we have two major size groups and some minor outliers. For example we could manually resize or crop the outliers and find a strategy for the two major groups. Let's take a look at the sizes again:"},{"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":{},"cell_type":"markdown","source":"For OSIC it's interesting that we have two different kind of patterns. In the RSNA competition all training images are of shape 512, 512. So you don't need to explore the image sizes further. But perhaps it's still useful to browse through the images to gain ideas for good augmentations."},{"metadata":{"trusted":true,"_kg_hide-input":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":{},"cell_type":"markdown","source":"Such kind of obvious patterns are always worth to take a look at. Perhaps we can find some kind of rule that allows us to create a good resizing strategy. Here is a small observer that allows you to browse through the files by just running the code cell again:"},{"metadata":{"trusted":true,"_kg_hide-input":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":{},"cell_type":"markdown","source":"Just run the following cell again, to observe the next batch of images:"},{"metadata":{"trusted":true},"cell_type":"code","source":"images = next(observer_iterator)","execution_count":null,"outputs":[]},{"metadata":{"_kg_hide-input":true,"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":{},"cell_type":"markdown","source":"### Insights\n\nBrowse through the different groups and you will find out for OSIC that...\n\n* The large squared image sizes often show resolutions higher than with 512 rows and 512 columns. Here resizing to the two major groups (512, 512) or (768, 768) makes sense.\n* In the non-squared cases a center crops should work best as they only have more background values but still belong to the major groups in the inside-scanner-region."},{"metadata":{},"cell_type":"markdown","source":"## Generating non-segmented, original raw image data <a class=\"anchor\" id=\"non_segmented_original\"></a>\n\nhttps://www.kaggle.com/allunia/osic-pulmonary-fibrosis-progression-huscans"},{"metadata":{},"cell_type":"markdown","source":"The least we can do is to resize and crop all those images that don't have a shape of (512,512) and to transform their values to Hounsfield units. Furthermore we can create smaller images by resizing to a desired shape. As far as I know this is not that easy as it sounds as we won't use common Image dtypes but rather original HU values without converting them. Most image preprocessing packages I found can't easily deal with numpy arrays. For this reason I used a workaround:\n\n* I changed the dytpe of my hu-transformed numpy-array scans to Int32\n* This way I can pass them to PILs Image class using fromarray with mode=\"I\"\n* Then I can use the resize and crop method of this package\n* After that the array turns back to Int16 as I like to reduce memory cost\n\nHmm... still feels not good and I'm also a bit confused. You can see that dicom can also be covered with most common Image preprocessing packages this way:\n\n* https://www.kaggle.com/onealbao/dicom-to-jpeg-conversion-kernel\n\nI still need to read about converting dicom or \"how to loose or not loose information\"! ;-)"},{"metadata":{"trusted":true},"cell_type":"code","source":"scan_properties.shape","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"176 different patients in train!"},{"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":{},"cell_type":"markdown","source":"Have to rework this part for RSNA!"},{"metadata":{"_kg_hide-input":true,"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":{},"cell_type":"markdown","source":"TODO!"},{"metadata":{"trusted":true,"_kg_hide-input":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,"_kg_hide-input":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":{"_kg_hide-input":true,"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,"_kg_hide-input":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":{"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}