{"cells":[{"metadata":{},"cell_type":"markdown","source":"# DICOMs to JPGs using FastAI method\n\n### JPGs are generated for both scan images and scans segmented using Watershed-Segmentation with markers\n\nThe process followed for generating the JPGs is from this amazing kernel by [Jeremy Howard](https://www.kaggle.com/jhoward) - https://www.kaggle.com/jhoward/don-t-see-like-a-radiologist-fastai. I've also explained breifly why this process is required.\n\nFor generating the segemented lung images Marker controlled watershed Segmentation is used from the kernel https://www.kaggle.com/aadhavvignesh/lung-segmentation-by-marker-controlled-watershed by [Aadhav Vignesh](https://www.kaggle.com/aadhavvignesh)\n\n## Datasets Used\n\n**DICOM Metadata** - generated using FastAI\nhttps://www.kaggle.com/gautham11/osic-pulmonary-fibrosis-progression-dicom-metadata\n\n**FastAI2 whls** - https://www.kaggle.com/vijayabhaskar96/fastai2-wheels\n\n## Dataset Generated\n\n**JPG images dataset** - https://www.kaggle.com/gautham11/pulmonaryfibrosisprogressionsegmentationjpgs","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"### Key Insight\n\nWhile normalizing the floating point image array to `np.uint8` for storing as jpg, care must be taken to not scale every image by it's `mean` and `std.dev`, because by this the mapping from DICOM values to normalized values will change from image to image. ","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"import os\nimport cv2\n\nimport pydicom\nimport pandas as pd\nimport numpy as np \nimport torch\nimport tensorflow as tf \nimport matplotlib.pyplot as plt \nfrom pathlib import Path\nimport scipy.ndimage as ndimage\nfrom skimage import measure, morphology, segmentation\nfrom scipy.ndimage.interpolation import zoom\nfrom PIL import Image \n\nfrom tqdm.notebook import tqdm\n%matplotlib inline","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"!pip install ../input/fastai2-wheels/fastcore-0.1.18-py3-none-any.whl -q\n!pip install ../input/fastai2-wheels/fastai2-0.0.17-py3-none-any.whl -q","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"from fastai2.basics import *\nfrom fastai2.medical.imaging import *\nimport cv2\nimport torch","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"data_dir = Path('/kaggle/input/osic-pulmonary-fibrosis-progression/')\npatient_paths = (data_dir/'train').ls()\nsample_patient = patient_paths[0]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"## Generated by\n# dicom_meta = pd.DataFrame()\n# for patient in tqdm(patient_paths):\n#     dicom_meta = pd.concat([dicom_meta, pd.DataFrame.from_dicoms(patient.ls())], ignore_index=True)\n\ndicom_meta = pd.read_pickle('../input/osic-pulmonary-fibrosis-progression-dicom-metadata/train_dicom_df')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"dicom_meta.head().T","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"## Code from https://www.kaggle.com/aadhavvignesh/lung-segmentation-by-marker-controlled-watershed \n\ndef generate_markers(image):\n    \"\"\"\n    Generates markers for a given image.\n    \n    Parameters: image\n    \n    Returns: Internal Marker, External Marker, Watershed Marker\n    \"\"\"\n    \n    #Creation of the internal Marker\n    marker_internal = image < -400\n    marker_internal = segmentation.clear_border(marker_internal)\n    marker_internal_labels = measure.label(marker_internal)\n    \n    areas = [r.area for r in measure.regionprops(marker_internal_labels)]\n    areas.sort()\n    \n    if len(areas) > 2:\n        for region in measure.regionprops(marker_internal_labels):\n            if region.area < areas[-2]:\n                for coordinates in region.coords:                \n                       marker_internal_labels[coordinates[0], coordinates[1]] = 0\n    \n    marker_internal = marker_internal_labels > 0\n    \n    # Creation of the External Marker\n    external_a = ndimage.binary_dilation(marker_internal, iterations=10)\n    external_b = ndimage.binary_dilation(marker_internal, iterations=55)\n    marker_external = external_b ^ external_a\n    \n    # Creation of the Watershed Marker\n    marker_watershed = np.zeros(image.shape, dtype=np.int)\n    marker_watershed += marker_internal * 255\n    marker_watershed += marker_external * 128\n    \n    return marker_internal, marker_external, marker_watershed\n\n\ndef seperate_lungs(image, iterations = 1):\n    \"\"\"\n    Segments lungs using various techniques.\n    \n    Parameters: image (Scan image), iterations (more iterations, more accurate mask)\n    \n    Returns: \n        - Segmented Lung\n        - Lung Filter\n        - Outline Lung\n        - Watershed Lung\n        - Sobel Gradient\n    \"\"\"\n    \n    # Store the start time\n    # start = time.time()\n    \n    marker_internal, marker_external, marker_watershed = generate_markers(image)\n    \n    \n    '''\n    Creation of Sobel Gradient\n    '''\n    \n    # Sobel-Gradient\n    sobel_filtered_dx = ndimage.sobel(image, 1)\n    sobel_filtered_dy = ndimage.sobel(image, 0)\n    sobel_gradient = np.hypot(sobel_filtered_dx, sobel_filtered_dy)\n    sobel_gradient *= 255.0 / np.max(sobel_gradient)\n    \n    \n    '''\n    Using the watershed algorithm\n    \n    \n    We pass the image convoluted by sobel operator and the watershed marker\n    to morphology.watershed and get a matrix matrix labeled using the \n    watershed segmentation algorithm.\n    '''\n    watershed = morphology.watershed(sobel_gradient, marker_watershed)\n    \n    '''\n    Reducing the image to outlines after Watershed algorithm\n    '''\n    outline = ndimage.morphological_gradient(watershed, size=(3,3))\n    outline = outline.astype(bool)\n    \n    \n    '''\n    Black Top-hat Morphology:\n    \n    The black top hat of an image is defined as its morphological closing\n    minus the original image. This operation returns the dark spots of the\n    image that are smaller than the structuring element. Note that dark \n    spots in the original image are bright spots after the black top hat.\n    '''\n    \n    # Structuring element used for the filter\n    blackhat_struct = [[0, 0, 1, 1, 1, 0, 0],\n                       [0, 1, 1, 1, 1, 1, 0],\n                       [1, 1, 1, 1, 1, 1, 1],\n                       [1, 1, 1, 1, 1, 1, 1],\n                       [1, 1, 1, 1, 1, 1, 1],\n                       [0, 1, 1, 1, 1, 1, 0],\n                       [0, 0, 1, 1, 1, 0, 0]]\n    \n    blackhat_struct = ndimage.iterate_structure(blackhat_struct, iterations)\n    \n    # Perform Black Top-hat filter\n    outline += ndimage.black_tophat(outline, structure=blackhat_struct)\n    \n    '''\n    Generate lung filter using internal marker and outline.\n    '''\n    lungfilter = np.bitwise_or(marker_internal, outline)\n    lungfilter = ndimage.morphology.binary_closing(lungfilter, structure=np.ones((5,5)), iterations=3)\n    \n    '''\n    Segment lung using lungfilter and the image.\n    '''\n    segmented = np.where(lungfilter == 1, image, -2000*np.ones(image.shape))\n    \n    return segmented, lungfilter, outline, watershed, sobel_gradient\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#DICOM Read utils\ndef fix_pxrepr(dcm):\n    if dcm.PixelRepresentation != 0 or dcm.RescaleIntercept<-100:\n        return dcm\n    x = dcm.pixel_array + 1000\n    px_mode = 4096\n    x[x>=px_mode] = x[x>=px_mode] - px_mode\n    dcm.PixelData = x.tobytes()\n    dcm.RescaleIntercept = -1000\n    return dcm\n\ndef read_dcm(path):\n    dcm = fix_pxrepr(Path(path).dcmread())\n    if dcm.Rows != 512 or dcm.Columns != 512: \n        dcm.zoom_to((512,512))\n    return dcm\n\n\n# Generate JPGs\ndef get_maskedlung_scan_img(path):\n    dcm = read_dcm(path)\n    _, lungmask, _, _, _ = seperate_lungs(dcm.scaled_px.numpy())\n    masked_lung = torch.tensor(np.where(lungmask, dcm.scaled_px, -2048))\n    return dcm, masked_lung\n\ndef save_maskedlung_lungimg(path, masked_path, img_path, masked_bins=None, bins=None):\n    dcm, masked_lung = get_maskedlung_scan_img(path)\n    save_jpg(masked_lung, masked_path, masked_bins)\n    save_jpg(dcm, img_path, bins)\n    \n    \ndef save_jpg(img, path, bins):\n    windows = [dicom_windows.lungs,dicom_windows.subdural]\n    img.save_jpg(path, windows, bins=bins)\n    \ndef read_jpg(path):\n    return np.array(Image.open(str(path)))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Generating Bins for Histogram scaling\n\nTo convert an image from Housfield values to a [0, 1] range, we split the Housfield pixel values of images into bins, such that each bin has an equal number of pixels. In order to maintain the same mapping from Hounsfield values to normalized value for all images, these bin boundaries are generated from Hounsfield values of multiple sampled images. This process is explained in detail in  https://www.kaggle.com/jhoward/don-t-see-like-a-radiologist-fastai","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"dicom_meta['image_format'] = dicom_meta.PixelRepresentation.astype(str) + '_' + dicom_meta.BitsStored.astype(str)\ndicom_meta['image_format'].value_counts().to_dict()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"format_df = dicom_meta[['image_format', 'fname']]\nselected_files = []\nfor img_format, n in format_df['image_format'].value_counts().to_dict().items():\n    selected_files += list(format_df.loc[format_df['image_format'] == str(img_format)].sample(min(1000, n // 3))['fname'])\n    \nselected_files[0:5]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def get_scan_freqhist_bins():\n    sample_images = []\n    for path in tqdm(selected_files):\n        try:\n            dcm = read_dcm(path)\n            sample_images.append(dcm.scaled_px)\n        except Exception as e:\n            print(e)\n    samples_arr = torch.stack(sample_images)\n    del sample_images\n    bins = samples_arr.freqhist_bins()\n    del samples_arr\n    return bins","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# bins = get_scan_freqhist_bins()\n\n# generated from the above code\nbins = torch.tensor([-4096., -3024., -2048., -2000., -1109., -1025., -1024., -1023., -1020.,\n        -1017., -1014., -1011., -1008., -1005., -1003., -1001.,  -998.,  -996.,\n         -993.,  -991.,  -988.,  -985.,  -981.,  -978.,  -973.,  -969.,  -964.,\n         -958.,  -952.,  -946.,  -940.,  -933.,  -926.,  -919.,  -912.,  -904.,\n         -895.,  -886.,  -875.,  -862.,  -847.,  -829.,  -806.,  -777.,  -738.,\n         -684.,  -607.,  -501.,  -386.,  -301.,  -247.,  -210.,  -183.,  -163.,\n         -148.,  -136.,  -126.,  -117.,  -109.,  -101.,   -93.,   -85.,   -77.,\n          -69.,   -60.,   -51.,   -42.,   -33.,   -24.,   -15.,    -6.,     0.,\n            5.,    13.,    20.,    27.,    34.,    41.,    48.,    55.,    62.,\n           70.,    79.,    90.,   102.,   118.,   137.,   162.,   193.,   237.,\n          301.,   417.,   750.,  1278.])\n\nplt.plot(bins, torch.linspace(0,1,len(bins)));\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def get_maskedscan_freqhist_bins():\n    sample_masked_images = []\n    for path in tqdm(selected_files):\n        try:\n            dcm = read_dcm(path)\n            dcm_np = dcm.scaled_px.numpy()\n            _, lungmask, _, _, _ = seperate_lungs(dcm_np)\n            masked_lung = torch.tensor(np.where(lungmask, dcm_np, -2048))\n            sample_masked_images.append(masked_lung)\n        except Exception as e:\n            print(e)\n    sample_mask_arr = torch.stack(sample_masked_images)\n    del sample_masked_images\n    masked_bins = sample_mask_arr.freqhist_bins()\n    del sample_mask_arr\n    return masked_bins","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# masked_bins = get_maskedscan_freqhist_bins()\n\nmasked_bins = torch.tensor([-2048., -1000.,  -942.,  -908.,  -879.,  -847.,  -808.,  -751.,  -656.,\n         -467.,   -81.,   251.])\n\nplt.plot(masked_bins, torch.linspace(0,1,len(masked_bins)));\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Testing on sample Patient","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"dcom = read_dcm(data_dir/'train/ID00228637202259965313869/100.dcm')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sample_patient","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"fig, axes = plt.subplots(2,2, figsize=[10, 10])\naxes[0][0].imshow(dcom.scaled_px)\naxes[0][0].set_title('Housfield values')\naxes[0][1].imshow(dcom.scaled_px.hist_scaled(bins))\naxes[0][1].set_title('Hist Scaled')\naxes[1][0].hist(dcom.scaled_px.flatten(), bins=bins)\naxes[1][1].hist(dcom.scaled_px.hist_scaled(bins).flatten())","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"_, lungmask, _, _, _ = seperate_lungs(dcom.scaled_px.numpy())\nmasked_lung = torch.tensor(np.where(lungmask, dcom.scaled_px, -2048))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"fig, axes = plt.subplots(2,2, figsize=[10, 10])\naxes[0][0].imshow(masked_lung)\naxes[0][0].set_title('Housfield values')\naxes[0][1].imshow(masked_lung.hist_scaled(masked_bins))\naxes[0][1].set_title('Hist Scaled')\naxes[1][0].hist(masked_lung.flatten(), bins=masked_bins)\naxes[1][1].hist(masked_lung.hist_scaled(bins).flatten())","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Normalized Array to JPGs\n\nIn order to fully utilize the 3 channels of JPG image, the first 2 channels are filled with windowed CT scans - https://radiopaedia.org/articles/windowing-ct, and the last channel is filled with histogram scaled values with the generated `bins` ","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"nchan_img = dcom.to_nchan([dicom_windows.lungs,dicom_windows.subdural], bins)\nnchan_masked_img = masked_lung.to_nchan([dicom_windows.lungs,dicom_windows.subdural], masked_bins)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"fig, axes = plt.subplots(1, 3, figsize=(15, 5))\naxes[0].imshow(nchan_img[0,:,:])\naxes[0].set_title('dicom_windows.lungs')\naxes[1].imshow(nchan_img[1,:,:])\naxes[1].set_title('dicom_windows.subdural')\naxes[2].imshow(nchan_img[2,:,:])\naxes[2].set_title('hist_scaled')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"fig, axes = plt.subplots(1, 3, figsize=(15, 5))\naxes[0].imshow(nchan_masked_img[0,:,:])\naxes[0].set_title('dicom_windows.lungs')\naxes[1].imshow(nchan_masked_img[1,:,:])\naxes[1].set_title('dicom_windows.subdural')\naxes[2].imshow(nchan_masked_img[2,:,:])\naxes[2].set_title('hist_scaled')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Generate Dataset","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"!mkdir 'osic-pulmonary-fibrosis-progression-segmentation-jpgs'\n!mkdir 'osic-pulmonary-fibrosis-progression-segmentation-jpgs/train'\n!mkdir 'osic-pulmonary-fibrosis-progression-segmentation-jpgs/train/scans'\n!mkdir 'osic-pulmonary-fibrosis-progression-segmentation-jpgs/train/masked_scans'","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"output_dir = Path('osic-pulmonary-fibrosis-progression-segmentation-jpgs/train')\noutput_dir_scans = output_dir/'scans'\noutput_dir_mscans = output_dir/'masked_scans'","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"patient_dirs = (data_dir/'train/').ls()\nlen(patient_dirs)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"import traceback\nimport shutil\n\ndef create_a_zip_archive(i, patient_dir_batch):\n    for patient in tqdm(patient_dir_batch):\n        patient_id = str(patient).split('/')[-1]\n        (output_dir_scans/patient_id).mkdir(parents=True, exist_ok=True)\n        (output_dir_mscans/patient_id).mkdir(parents=True, exist_ok=True)\n        for img_path in patient.glob('*'):\n            try:\n                img_id = str(img_path).split('/')[-1][:-4]\n                if not (output_dir_mscans/patient_id/f'{img_id}.jpg').exists():\n                    save_maskedlung_lungimg(\n                        img_path, \n                        output_dir_mscans/patient_id/f'{img_id}.jpg',\n                        output_dir_scans/patient_id/f'{img_id}.jpg',\n                        masked_bins,\n                        bins\n                    )\n            except Exception as e:\n                traceback.print_exc()\n                print(e)\n                print(img_path)\n    shutil.make_archive(f'lung_segmentation_scans_{i}', 'zip', str(output_dir))\n    shutil.rmtree(str(output_dir_scans))\n    shutil.rmtree(str(output_dir_mscans))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"n = 10\nstart_n = 0\nfor i, n_start in enumerate(range(start_n, len(patient_dirs), n)):\n    create_a_zip_archive(i + int(start_n / n), patient_dirs[n_start: n_start + n])\n    # Link to generated dataset added\n    break","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}