{"cells":[{"metadata":{},"cell_type":"markdown","source":"# Preprocessing: Spatially adjacent RGB images and cropping\nThis is a follow on from my previous notebook https://www.kaggle.com/anjum48/reconstructing-3d-volumes-from-metadata\n\nThis is also a part of our 5th place solution https://www.kaggle.com/c/rsna-intracranial-hemorrhage-detection/discussion/117232#latest-672657"},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"import os\nimport pandas as pd\nimport numpy as np\nimport pydicom\nfrom PIL import Image, ImageFile\nimport matplotlib.pylab as plt\nfrom tqdm import tqdm_notebook\n%matplotlib inline\n\ndata_path = \"../input/rsna-intracranial-hemorrhage-detection/rsna-intracranial-hemorrhage-detection\"\nmetadata_path = \"../input/rsna-ich-metadata\"\nos.listdir(metadata_path)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"os.listdir(data_path)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Load the metadata (see my previous notebook to see how this dataset was generated)\n\nDataset: https://www.kaggle.com/anjum48/rsna-ich-metadata"},{"metadata":{"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","trusted":true},"cell_type":"code","source":"train_metadata = pd.read_parquet(f'{metadata_path}/train_metadata.parquet.gzip')\ntest_metadata = pd.read_parquet(f'{metadata_path}/test_metadata.parquet.gzip')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Create triplets of images\nBased on `StudyInstanceUID` and sorting on `ImagePositionPatient` it is possible to reconstruct 3D volumes for each study. However since each study contained a variable number of axial slices (between 20-60) this makes it difficult to create a architecture that implements 3D convolutions. Instead, triplets of images were created from the 3D volumes to represent the RGB channels of an image, i.e. the green channel being the target image and the red & blue channels being the adjacent images. If an image was at the edge of the volume, then the green channel was repeated. This is essentially a 3D volume but only using 3 axial slices. At this stage no windowing was applied and the image is retained in Hounsfield units."},{"metadata":{"trusted":true},"cell_type":"code","source":"def build_triplets(metadata):\n    metadata.sort_values(by=\"ImagePositionPatient_2\", inplace=True, ascending=False)\n    studies = metadata.groupby(\"StudyInstanceUID\")\n    triplets = []\n\n    for study_name, study_df in tqdm_notebook(studies):\n        padded_names = np.pad(study_df.index, (1, 1), 'edge')\n\n        for i, img in enumerate(padded_names[1:-1]):\n            t = [padded_names[i], img, padded_names[i + 2]]\n            triplets.append(t)\n\n    return pd.DataFrame(triplets, columns=[\"red\", \"green\", \"blue\"])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_triplets = build_triplets(train_metadata)\ntest_triplets = build_triplets(test_metadata)\ntrain_triplets.to_csv(\"train_triplets.csv\")\ntest_triplets.to_csv(\"stage_1_test_triplets.csv\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_triplets.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Creating the RGB image\nLet's construct an RGB image using these triplets"},{"metadata":{"trusted":true},"cell_type":"code","source":"def prepare_dicom(dcm):\n    \"\"\"\n    Converts a DICOM object to a 16-bit Numpy array (in Hounsfield units)\n    :param dcm: DICOM Object\n    :return: Numpy array in int16\n    \"\"\"\n\n    try:\n        # https://www.kaggle.com/jhoward/cleaning-the-data-for-rapid-prototyping-fastai\n        if dcm.BitsStored == 12 and dcm.PixelRepresentation == 0 and dcm.RescaleIntercept > -100:\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\n        pixels = dcm.pixel_array.astype(np.float32) * dcm.RescaleSlope + dcm.RescaleIntercept\n    except ValueError as e:\n        print(\"ValueError with\", dcm.SOPInstanceUID, e)\n        return np.zeros((512, 512))\n\n    # Pad the image if it isn't square\n    if pixels.shape[0] != pixels.shape[1]:\n        (a, b) = pixels.shape\n        if a > b:\n            padding = ((0, 0), ((a - b) // 2, (a - b) // 2))\n        else:\n            padding = (((b - a) // 2, (b - a) // 2), (0, 0))\n        pixels = np.pad(pixels, padding, mode='constant', constant_values=0)\n\n    return pixels.astype(np.int16)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"channels = train_triplets.iloc[3079]\nchannels","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"rgb = []\n\nfor image_id in channels:\n    dcm = pydicom.dcmread(os.path.join(data_path, \"stage_2_train\", image_id + \".dcm\"))\n    rgb.append(prepare_dicom(dcm))\n    \nimg = np.stack(rgb, -1)\nimg = np.clip(img, 0, 255).astype(np.uint8) ","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"On the final line above, we actually clip the image between 0 & 255 Hounfield units. This is ok, since most of the important features in the image are between this range.\n\nSince the image is still technically in Hounsfield units, you can apply a window to it later (e.g. brain, subdural etc.), however since we have already clipped this image, this may impact the bone window.\n\nLet's check out our image:"},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.figure(figsize=(8, 8))\nplt.imshow(img);","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Cropping the image\nAs you can see there is a) a lot of black space which we don't want to waste our valuable GPU operations on and b) the CT scanner headrest. Let's get rid of both of these using `scipy.ndimage.label`"},{"metadata":{"trusted":true},"cell_type":"code","source":"from scipy import ndimage","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"labeled_blobs, number_of_blobs = ndimage.label(img)\nblob_sizes = np.bincount(labeled_blobs.flatten())\nnumber_of_blobs","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"blob_sizes","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"This labels adjacent \"blobs\" in the image, i.e. groups of connected pixels.\n\nIn this example we have 39 blobs, each with different sizes. Lets see the first blob:"},{"metadata":{"trusted":true},"cell_type":"code","source":"blob_0 = labeled_blobs == 0  # label 0 has 593185 connected pixels\nblob_0 = np.max(blob_0, axis=-1)\n\nplt.figure(figsize=(8, 8))\nplt.imshow(blob_0);","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Ok, the largest blob appears to be the background. Let's try the next largest blob"},{"metadata":{"trusted":true},"cell_type":"code","source":"blob_1 = labeled_blobs == 1  # label 1 has 173243 connected pixels\nblob_1 = np.max(blob_1, axis=-1)\n\nplt.figure(figsize=(8, 8))\nplt.imshow(blob_1);","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Great! We've found that the **2nd largest blob** is the head and there is no headrest. We can use this mask to define our cropping extent (`x_min, x_max, y_min, y_max`) and create a cropping function"},{"metadata":{"trusted":true},"cell_type":"code","source":"class CropHead(object):\n    def __init__(self, offset=10):\n        \"\"\"\n        Crops the head by labelling the objects in an image and keeping the second largest object (the largest object\n        is the background). This method removes most of the headrest\n\n        Originally made as a image transform for use with PyTorch, but too slow to run on the fly :(\n        :param offset: Pixel offset to apply to the crop so that it isn't too tight\n        \"\"\"\n        self.offset = offset\n\n    def crop_extents(self, img):\n        try:\n            if type(img) != np.array:\n                img_array = np.array(img)\n            else:\n                img_array = img\n\n            labeled_blobs, number_of_blobs = ndimage.label(img_array)\n            blob_sizes = np.bincount(labeled_blobs.flatten())\n            head_blob = labeled_blobs == np.argmax(blob_sizes[1:]) + 1  # The number of the head blob\n            head_blob = np.max(head_blob, axis=-1)\n\n            mask = head_blob == 0\n            rows = np.flatnonzero((~mask).sum(axis=1))\n            cols = np.flatnonzero((~mask).sum(axis=0))\n\n            x_min = max([rows.min() - self.offset, 0])\n            x_max = min([rows.max() + self.offset + 1, img_array.shape[0]])\n            y_min = max([cols.min() - self.offset, 0])\n            y_max = min([cols.max() + self.offset + 1, img_array.shape[1]])\n\n            return x_min, x_max, y_min, y_max\n        except ValueError:\n            return 0, 0, -1, -1\n\n    def __call__(self, img):\n        \"\"\"\n        Crops a CT image to so that as much black area is removed as possible\n        :param img: PIL image\n        :return: Cropped image\n        \"\"\"\n\n        x_min, x_max, y_min, y_max = self.crop_extents(img)\n\n        try:\n            if type(img) != np.array:\n                img_array = np.array(img)\n            else:\n                img_array = img\n\n            return Image.fromarray(np.uint8(img_array[x_min:x_max, y_min:y_max]))\n        except ValueError:\n            return img\n\n    def __repr__(self):\n        return self.__class__.__name__ + '(offset={})'.format(self.offset)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"crop_head = CropHead()\nimg_cropped = crop_head(img)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# The final product\nLet's check our cropped image"},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.figure(figsize=(8, 8))\nplt.imshow(img_cropped);","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Looks good! Most of the black space is removed and most of the headrest is gone. The image can then be resized into a square (e.g. 224x224) for use in a CNN"},{"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":1}