{"cells":[{"metadata":{},"cell_type":"markdown","source":"# 1. Pre-processing and Feature Engineering\nIn this new version, we stripped down all unnecessary pre-processing, removing padding, to tensor, zero center/standarization.\nWe must use `MaskWatershed`, so that's the choice here.\nFinally, we implemented the following features:\n- Lung volume\n- Chest circumference\n- Lung height\n- HU histogram features: mean, stddev, skew, kurtosis, and % of HU in each bin","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"# 2. Imports","execution_count":null},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"from datetime import timedelta\nimport numpy as np\nimport pandas as pd\nfrom pathlib import Path\nimport pydicom\nfrom scipy.ndimage.interpolation import zoom\nfrom scipy.stats import skew, kurtosis\nimport scipy.ndimage as ndimage\nfrom skimage.measure import label, regionprops\nfrom skimage.segmentation import clear_border\nfrom skimage import measure, morphology, segmentation\nfrom time import time, sleep\nimport torch\nfrom torch.utils.data import Dataset\nfrom torchvision import transforms\nfrom tqdm import trange\nimport warnings\nimport matplotlib.pyplot as plt\nfrom matplotlib import cm\nimport pytest\nimport multiprocessing as mp\nfrom time import time\nfrom datetime import timedelta\nfrom functools import partial\nfrom tqdm import tqdm","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# 3. Dataset interface","execution_count":null},{"metadata":{"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","trusted":true},"cell_type":"code","source":"class CTScansDataset(Dataset):\n    def __init__(self, root_dir, transform=None):\n        self.root_dir = Path(root_dir)\n        self.patients = [p for p in self.root_dir.glob('*') if p.is_dir()]\n        self.transform = transform\n\n    def __len__(self):\n        return len(self.patients)\n\n    def __getitem__(self, idx):\n        if torch.is_tensor(idx):\n            idx = idx.tolist()\n\n        image, metadata, msbs = self.load_scan(self.patients[idx])\n        sample = {'image': image, 'metadata': metadata, 'msbs': msbs}\n        if self.transform:\n            sample = self.transform(sample)\n\n        return sample\n\n    def save(self, path):\n        t0 = time()\n        Path(path).mkdir(exist_ok=True, parents=True)\n        print('Saving pre-processed dataset to disk')\n        sleep(1)\n        cum = 0\n        var = 0\n\n        bar = trange(len(self))\n        for i in bar:\n            sample = self[i]\n            image, data = sample['image'], sample['metadata']\n            cum += torch.mean(image).item()\n            var += torch.var(image).item()\n\n            bar.set_description(f'Saving CT scan {data.PatientID}')\n            fname = Path(path) / f'{data.PatientID}.pt'\n            torch.save(image, fname)\n\n        sleep(1)\n        bar.close()\n        print(f'Done! Time {timedelta(seconds=time() - t0)}\\n'\n              f'Mean value: {cum / len(self)}\\n'\n              f'Std value: {(var / len(self)) ** (1/2)}')\n\n    def get_patient(self, patient_id):\n        patient_ids = [str(p.stem) for p in self.patients]\n        return self.__getitem__(patient_ids.index(patient_id))\n\n    @staticmethod\n    def load_scan(path):\n        raw_slices = [pydicom.read_file(p) for p in path.glob('*.dcm')]\n        final_slices = [s for s in raw_slices\n                        if hasattr(s, 'ImagePositionPatient')]\n\n        if len(final_slices) > 0:\n            final_slices.sort(key=lambda x: float(x.ImagePositionPatient[2]))\n            msbs = 0\n            for i in range(len(final_slices) - 1):\n                a = final_slices[i].ImagePositionPatient[2]\n                b = final_slices[i + 1].ImagePositionPatient[2]\n                msbs += abs(a - b)\n\n            msbs /= len(final_slices)\n            image = np.stack([s.pixel_array.astype(float) for s in final_slices])\n            metadata = final_slices[0]\n\n        else:\n            # Guess based on eyballing other mean spaces between slices:\n            msbs = None\n            if raw_slices[0].PatientID == 'ID00132637202222178761324':\n                msbs = 0.7\n            elif raw_slices[0].PatientID == 'ID00128637202219474716089':\n                msbs = 5.\n\n            warnings.warn(f'Patient {raw_slices[0].PatientID} CT scan does '\n                          f'not have \"ImagePositionPatient\". Assuming '\n                          f'filenames in the right scan order. Also, assuming'\n                          f'mean space between slices of {msbs}')\n\n            image = np.stack([s.pixel_array.astype(float) for s in raw_slices])\n            metadata = raw_slices[0]\n\n        return image, metadata, msbs","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# 4. Pre-processing functions","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"class CropBoundingBox:\n    @staticmethod\n    def bounding_box(img3d: np.array):\n        mid_img = img3d[int(img3d.shape[0] / 2)]\n        same_first_row = (mid_img[0, :] == mid_img[0, 0]).all()\n        same_first_col = (mid_img[:, 0] == mid_img[0, 0]).all()\n        if same_first_col and same_first_row:\n            return True\n        else:\n            return False\n\n    def __call__(self, sample):\n        image, data = sample['image'], sample['metadata']\n        if not self.bounding_box(image):\n            return sample\n\n        mid_img = image[int(image.shape[0] / 2)]\n        r_min, r_max = None, None\n        c_min, c_max = None, None\n        for row in range(mid_img.shape[0]):\n            if not (mid_img[row, :] == mid_img[0, 0]).all() and r_min is None:\n                r_min = row\n            if (mid_img[row, :] == mid_img[0, 0]).all() and r_max is None \\\n                    and r_min is not None:\n                r_max = row\n                break\n\n        for col in range(mid_img.shape[1]):\n            if not (mid_img[:, col] == mid_img[0, 0]).all() and c_min is None:\n                c_min = col\n            if (mid_img[:, col] == mid_img[0, 0]).all() and c_max is None \\\n                    and c_min is not None:\n                c_max = col\n                break\n\n        image = image[:, r_min:r_max, c_min:c_max]\n        return {'image': image, 'metadata': data, 'msbs': sample['msbs']}","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"class ConvertToHU:\n    def __call__(self, sample):\n        image, data = sample['image'], sample['metadata']\n\n        img_type = data.ImageType\n        is_hu = img_type[0] == 'ORIGINAL' and not (img_type[2] == 'LOCALIZER')\n        # if not is_hu:\n        #     warnings.warn(f'Patient {data.PatientID} CT Scan not cannot be'\n        #                   f'converted to Hounsfield Units (HU).')\n\n        intercept = data.RescaleIntercept\n        slope = data.RescaleSlope\n        image = (image * slope + intercept).astype(np.int16)\n        return {'image': image, 'metadata': data, 'msbs': sample['msbs']}","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"class Resample:\n    def __init__(self, new_spacing=(1, 1, 1)):\n        \"\"\"new_spacing means now much every pixel represent in mm, in each\n        dimension. E.g. 2, 2, 2 means every pixel represent 2mm in  every\n        dimension.\n        \"\"\"\n        assert isinstance(new_spacing, tuple)\n        self.new_spacing = new_spacing\n\n    def __call__(self, sample):\n        image, data, msbs = sample['image'], sample['metadata'], sample['msbs']\n\n        spacing = np.array([msbs] + list(data.PixelSpacing), dtype=np.float32)\n        resize_factor = spacing / self.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\n        image = zoom(image, real_resize_factor, mode='nearest')\n        return {'image': image, 'metadata': data}","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"class Clip:\n    def __init__(self, bounds=(-1000, 500)):\n        self.min = min(bounds)\n        self.max = max(bounds)\n\n    def __call__(self, sample):\n        image, data = sample['image'], sample['metadata']\n        image[image < self.min] = self.min\n        image[image > self.max] = self.max\n        return {'image': image, 'metadata': data}","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"class MaskWatershed:\n    def __init__(self, min_hu, iterations, show_tqdm):\n        self.min_hu = min_hu\n        self.iterations = iterations\n        self.show_tqdm = show_tqdm\n\n    def __call__(self, sample):\n        image, data = sample['image'], sample['metadata']\n\n        stack = []\n        if self.show_tqdm:\n            bar = trange(image.shape[0])\n            bar.set_description(f'Masking CT scan {data.PatientID}')\n        else:\n            bar = range(image.shape[0])\n        for slice_idx in bar:\n            sliced = image[slice_idx]\n            stack.append(self.seperate_lungs(sliced, self.min_hu,\n                                             self.iterations))\n\n        return {\n            'image': np.stack(stack),\n            'metadata': sample['metadata']\n        }\n\n    @staticmethod\n    def seperate_lungs(image, min_hu, iterations):\n        h, w = image.shape[0], image.shape[1]\n\n        marker_internal, marker_external, marker_watershed = MaskWatershed.generate_markers(image)\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        watershed = morphology.watershed(sobel_gradient, marker_watershed)\n\n        outline = ndimage.morphological_gradient(watershed, size=(3,3))\n        outline = outline.astype(bool)\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        lungfilter = np.bitwise_or(marker_internal, outline)\n        lungfilter = ndimage.morphology.binary_closing(lungfilter, structure=np.ones((5,5)), iterations=3)\n\n        segmented = np.where(lungfilter == 1, image, min_hu * np.ones((h, w)))\n\n        return segmented\n\n    @staticmethod\n    def generate_markers(image, threshold=-400):\n        h, w = image.shape[0], image.shape[1]\n\n        marker_internal = image < threshold\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((h, w), dtype=np.int)\n        marker_watershed += marker_internal * 255\n        marker_watershed += marker_external * 128\n\n        return marker_internal, marker_external, marker_watershed","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# 5. Feature computations","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"class LungVolume:\n    def __init__(self, unit, spacing, min_HU):\n        self.unit = unit\n        self.spacing = spacing\n        self.min_HU = min_HU\n\n    def __call__(self, sample):\n        image = sample['image']\n\n        volume = (image > self.min_HU).sum()\n        volumes_in_perc = volume / np.prod(image.shape)\n        volumes_in_litr = volumes_in_perc * np.prod(image.shape) * \\\n                          np.prod(self.spacing) / 10**6\n\n        if self.unit == '%':\n            sample['volume'] = volumes_in_perc\n        elif self.unit == 'l':\n            sample['volume'] = volumes_in_litr\n\n        return sample","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"class ChestCircumference:\n    def __init__(self, spacing, min_HU):\n        self.spacing = spacing\n        self.min_HU = min_HU\n\n    def ellipsis_perimeter(self, img):\n        # trim slice, removing zero rows and cols\n        img = img[~(img == self.min_HU).all(1)].T\n        img = img[~(img == self.min_HU).all(1)]\n        # ellipsis perimeter approximation 2\n        # https://www.mathsisfun.com/geometry/ellipse-perimeter.html\n        a, b = img.shape\n        a /= 2\n        b /= 2\n        p = np.pi * (3 * (a + b) - np.sqrt((3 * a + b) * (a + 3 * b)))\n        return p\n\n    def __call__(self, sample):\n        image = sample['image']\n\n        perimeters = []\n        for slice_idx in range(image.shape[0]):\n            slc = image[slice_idx]\n            if slc.max() != self.min_HU:  # There's something in the slice\n                perimeters.append(self.ellipsis_perimeter(slc))\n            else:\n                perimeters.append(0)\n\n        # Transform in mm\n        cc_in_mm = max(perimeters) * self.spacing[1] * self.spacing[2]\n        sample['chest'] = int(cc_in_mm)\n        return sample","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"class LungHeight:\n    def __init__(self, spacing):\n        self.spacing = spacing\n\n    def __call__(self, sample):\n        image = sample['image']\n        sample['height'] = image.shape[0] * self.spacing[0]\n        return sample","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"class HUHistogram:\n    def __init__(self, bounds, bins):\n        self.bounds = bounds\n        self.bins = bins\n\n    def __call__(self, sample):\n        image = sample['image']\n        pixels = image[image != min(self.bounds)]\n\n        if pixels.shape[0] > 0:\n            sample['mean'] = pixels.mean()\n            sample['std'] = pixels.std()\n            sample['skew'] = skew(pixels)\n            sample['kurtosis'] = kurtosis(pixels)\n\n            hist, _ = np.histogram(pixels, bins=self.bins, range=self.bounds)\n            sample['hist'] = hist / pixels.size  # Normalized\n            # Density = True in np.histogram doesn't work as expected\n        else:\n            sample['mean'] = 0\n            sample['std'] = 0\n            sample['skew'] = 0\n            sample['kurtosis'] = 0\n            sample['hist'] = np.zeros(self.bins)\n\n        return sample","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# 6. Visualizing features","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"def gennerate_row(dataset, i):\n    sample = dataset[i]\n    df = pd.DataFrame(columns=np.arange(8).tolist())\n    df.loc[0] = [sample['metadata'].PatientID, sample[\"volume\"],\n                 sample[\"chest\"], sample[\"height\"], sample[\"mean\"],\n                 sample[\"std\"], sample[\"skew\"], sample[\"kurtosis\"]]\n\n    dfhist = pd.DataFrame()\n    dfhist[0] = sample['hist']\n    dfhist = dfhist.T\n    df = pd.concat([df, dfhist], axis=1, ignore_index=True)\n\n    columns = ['PatientID', 'Volume', 'Chest', 'Height', 'Mean',\n               'Std', 'Skew', 'Kurtosis']\n    columns += [f'bin{b + 1}' for b in range(len(dfhist.columns))]\n    df.columns = columns\n\n    return df","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"t0 = time()\ndataset = CTScansDataset(\n    root_dir='../input/osic-pulmonary-fibrosis-progression/test',\n    transform=transforms.Compose([\n        CropBoundingBox(),\n        ConvertToHU(),\n        Resample(new_spacing=(1, 1, 1)),\n        Clip(bounds=(-1000, 400)),\n        MaskWatershed(min_hu=-1000, iterations=1, show_tqdm=False),\n        # Features\n        LungVolume(unit='l',\n                   spacing=(1, 1, 1),\n                   min_HU=-1000),\n        ChestCircumference(spacing=(1, 1, 1), min_HU=-1000),\n        LungHeight(spacing=(1, 1, 1)),\n        HUHistogram(bounds=(-1000, 400), bins=14)\n    ]))\n\nf = partial(gennerate_row, dataset)\nwith mp.Pool(processes=mp.cpu_count()) as pool:\n    results = list(tqdm(pool.imap_unordered(f, range(len(dataset))),\n                        total=len(dataset)))\n\nfeatures = pd.concat(results, ignore_index=True)\nfeatures.to_csv('features.csv', index=False)\n\nprint(f'\\nDone! Time: {timedelta(seconds=time() - t0)}')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"features","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"pd.read_csv('../input/osic-pulmonary-fibrosis-progression/sample_submission.csv').to_csv('submission.csv', index=False)","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}