{"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_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# RSNA'23 | I. Series Normalisation\n\nThis notebook resamples the series and segmentations in the training dataset.\n\nOther notebooks and datasets in this analysis:\n1. [RSNA'23 | I. Series Normalisation](https://www.kaggle.com/code/bvinning/rsna-23-i-series-normalisation) -> [RSNA'23 | Normalised Series](https://www.kaggle.com/datasets/bvinning/rsna-2023-normalised-series)\n2. [RSNA'23 | II. Predicting Segmentations with VNet](https://www.kaggle.com/code/bvinning/rsna-23-ii-predicting-segmentations-with-vnet) -> [RSNA'23 | Predicted Segmentations](https://www.kaggle.com/datasets/bvinning/rsna-2023-predicted-segmentations)\n3. [RSNA'23 | III. Applying Organ Segmentations](https://www.kaggle.com/code/bvinning/rsna-23-iii-applying-organ-segmentations) -> [RSNA'23 | Segmented Series](https://www.kaggle.com/datasets/bvinning/rsna-2023-segmented-series)\n4. [RSNA'23 | IV. Predicting Injuries with a 3D CNN](https://www.kaggle.com/code/bvinning/rsna-23-iv-predicting-injuries-with-a-3d-cnn) -> [RSNA'23 | Injury Classifier Model](https://www.kaggle.com/datasets/bvinning/rsna-2023-injury-classifier-model)\n5. [RSNA'23 | V. Final Submission](https://www.kaggle.com/code/bvinning/rsna-23-v-final-submission)","metadata":{}},{"cell_type":"code","source":"# Load libraries.\nimport os\nimport glob\n\nimport joblib\nimport pandas as pd\nimport numpy as np\nimport nibabel as nib\nimport matplotlib.pyplot as plt\nimport pydicom\n\n\nRANDOM_SEED = 777\nRANDOM_STATE = np.random.seed(RANDOM_SEED)","metadata":{"execution":{"iopub.status.busy":"2023-09-20T10:30:42.973662Z","iopub.execute_input":"2023-09-20T10:30:42.974105Z","iopub.status.idle":"2023-09-20T10:30:43.286031Z","shell.execute_reply.started":"2023-09-20T10:30:42.974062Z","shell.execute_reply":"2023-09-20T10:30:43.284876Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Load the train and test meta data.\nmeta_dtypes = {\n    'series_id': int,\n    'patient_id': int,\n    'aortic_hu': float,\n    'incomplete_organ': bool,\n}\n\n\ndef load_meta(train):\n    if train:\n        path = '/kaggle/input/rsna-2023-abdominal-trauma-detection/train_series_meta.csv'\n    else:\n        path = '/kaggle/input/rsna-2023-abdominal-trauma-detection/test_series_meta.csv'\n    df = pd.read_csv(path, dtype=meta_dtypes)\n    assert not df.isnull().any().any()\n    return df\n\n\ntrain_meta = load_meta(train=True)","metadata":{"execution":{"iopub.status.busy":"2023-09-20T10:30:43.288221Z","iopub.execute_input":"2023-09-20T10:30:43.288578Z","iopub.status.idle":"2023-09-20T10:30:43.323029Z","shell.execute_reply.started":"2023-09-20T10:30:43.288546Z","shell.execute_reply":"2023-09-20T10:30:43.322023Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import warnings\nimport shutil\n\n\nclass FileSystem:\n    def __init__(self):\n        self.working_dir = '/kaggle/working/'\n        self.cache_dir = os.path.join(self.working_dir, 'cache')\n        self.input_dir = '/kaggle/input/rsna-2023-abdominal-trauma-detection'\n        self.setup()\n    \n    @property\n    def local(self):\n        return os.path.join(self.working_dir, 'local')\n\n    @property\n    def local_series(self):\n        return os.path.join(self.local, 'series')\n\n    @property\n    def local_segmentations(self):\n        return os.path.join(self.local, 'segmentations')\n\n    @property\n    def local_segmented_series(self):\n        return os.path.join(self.local, 'segmented_series')\n\n    def setup(self):\n        pass\n\n    def teardown(self):\n        pass","metadata":{"execution":{"iopub.status.busy":"2023-09-20T10:30:43.32497Z","iopub.execute_input":"2023-09-20T10:30:43.325425Z","iopub.status.idle":"2023-09-20T10:30:43.335149Z","shell.execute_reply.started":"2023-09-20T10:30:43.325383Z","shell.execute_reply":"2023-09-20T10:30:43.333908Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import functools\n\n\nclass SeriesInfo:\n    def __init__(self, fs):\n        self.fs = fs\n\n    @functools.cached_property\n    def train_segmentations(self):\n        return self.get_series_ids(root_dir=f'{FILESYSTEM.input_dir}/segmentations')\n\n    @functools.cached_property\n    def train_series(self):\n        paths = glob.glob('*/*', root_dir=os.path.join(self.fs.input_dir, 'train_images'))\n        return {int(os.path.split(path)[-1]) for path in paths}\n    \n    @functools.cached_property\n    def test_series(self):\n        paths = glob.glob('*/*', root_dir=os.path.join(self.fs.input_dir, 'test_images'))\n        return {int(os.path.split(path)[-1]) for path in paths}\n\n    @functools.cached_property\n    def series(self):\n        return self.train_series + self.test_series\n\n    @property\n    def local_series(self):\n        return self.get_series_ids(root_dir=self.fs.local_series)\n\n    @property\n    def local_segmentations(self):\n        return self.get_series_ids(root_dir=self.fs.local_segmentations)\n\n    @property\n    def local_segmented_series(self):\n        return self.get_series_ids(root_dir=self.fs.local_segmented_series)\n    \n    @staticmethod\n    def get_series_ids(root_dir):\n        return {int(os.path.splitext(path)[0]) for path in glob.glob('*', root_dir=root_dir)}","metadata":{"execution":{"iopub.status.busy":"2023-09-20T10:30:43.338501Z","iopub.execute_input":"2023-09-20T10:30:43.339126Z","iopub.status.idle":"2023-09-20T10:30:43.352165Z","shell.execute_reply.started":"2023-09-20T10:30:43.339084Z","shell.execute_reply":"2023-09-20T10:30:43.351118Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import yaml\nimport dataclasses\nimport functools\n\n\n@dataclasses.dataclass(frozen=True)\nclass SeriesNormalisationConfig:\n    _series_dtype_out: str = 'uint16'\n    _segmentation_dtype_out: str = 'uint8'\n    shape: tuple[int] = dataclasses.field(default=(128, 128, 128)) # the shape after resampling\n    roi: tuple[int] = dataclasses.field(default=(350, 300, 400)) # region of interest in mm\n\n    @classmethod\n    def load(cls, path):\n        with open(path, 'r') as f:\n            data = yaml.safe_load(f)\n        return cls(**data)\n\n    @functools.cached_property\n    def series_dtype_in(self):\n        return np.dtype(self._series_dtype_in).type\n    \n    @functools.cached_property\n    def series_dtype_out(self):\n        return np.dtype(self._series_dtype_out).type\n    \n    @functools.cached_property\n    def segmentation_dtype_out(self):\n        return np.dtype(self._segmentation_dtype_out).type\n\n    @functools.cached_property\n    def series_dtype_out_max(self):\n        return np.iinfo(self.series_dtype_out).max\n\n    def dump(self):\n        data = dataclasses.asdict(self)\n        with open('config.yaml', 'w') as f:\n            yaml.safe_dump(data, f)","metadata":{"execution":{"iopub.status.busy":"2023-09-20T10:30:43.355568Z","iopub.execute_input":"2023-09-20T10:30:43.356035Z","iopub.status.idle":"2023-09-20T10:30:43.372231Z","shell.execute_reply.started":"2023-09-20T10:30:43.356Z","shell.execute_reply":"2023-09-20T10:30:43.370918Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"FILESYSTEM = FileSystem()\nSERIES_INFO = SeriesInfo(FILESYSTEM)\nCONFIG = SeriesNormalisationConfig()","metadata":{"execution":{"iopub.status.busy":"2023-09-20T10:30:43.373945Z","iopub.execute_input":"2023-09-20T10:30:43.374337Z","iopub.status.idle":"2023-09-20T10:30:43.388916Z","shell.execute_reply.started":"2023-09-20T10:30:43.374306Z","shell.execute_reply":"2023-09-20T10:30:43.387517Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"tags_dtypes = {    \n    'BitsAllocated': int,\n    'BitsStored': int,\n    'Columns': int,\n    'ImageOrientationPatient': str,\n    'ImagePositionPatient': str,\n    'PatientID': int,\n    'PixelRepresentation': int,\n    'PixelSpacing': str,\n    'Rows': int,\n    'SeriesInstanceUID': str,\n    'SliceThickness': float,\n    'path': str,\n}\n\ndef load_tags(train):\n    if train:\n        path = '/kaggle/input/rsna-2023-abdominal-trauma-detection/train_dicom_tags.parquet'\n    else:\n        path = '/kaggle/input/rsna-2023-abdominal-trauma-detection/test_dicom_tags.parquet'\n    \n    df = pd.read_parquet(path, columns=list(tags_dtypes))\n    for k, v in tags_dtypes.items():\n        df[k] = df[k].astype(v)\n    assert not df.isnull().any().any()\n    df['SeriesID'] = df.SeriesInstanceUID.str.split('.').str[-1].astype(int)\n    # Cast strings as arrays.\n    for c in ['ImageOrientationPatient', 'ImagePositionPatient', 'PixelSpacing']:\n        expanded_column = df.pop(c).str.strip('[]').str.split(',', expand=True).astype(np.float32)\n        expanded_column_names = [f'{c}{i}' for i in range(len(expanded_column.columns))]\n        df[expanded_column_names] = expanded_column\n    # Make `path` absolute path.\n    df['path'] = df.path.apply(lambda x: os.path.join(FILESYSTEM.input_dir, x))\n    return df\n\n\ntrain_tags = load_tags(train=True)","metadata":{"execution":{"iopub.status.busy":"2023-09-20T10:30:43.390866Z","iopub.execute_input":"2023-09-20T10:30:43.39166Z","iopub.status.idle":"2023-09-20T10:31:23.551193Z","shell.execute_reply.started":"2023-09-20T10:30:43.391623Z","shell.execute_reply":"2023-09-20T10:31:23.550083Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def sort_tags(df):\n    a = df[\n        ['ImageOrientationPatient0', 'ImageOrientationPatient1', 'ImageOrientationPatient2']\n    ].values\n    b = df[\n        ['ImageOrientationPatient3', 'ImageOrientationPatient4', 'ImageOrientationPatient5']\n    ].values\n    ia = np.cross(a, b) # The imaging axis.\n    ipp = ['ImagePositionPatient0', 'ImagePositionPatient1', 'ImagePositionPatient2']\n    df['DistanceProjection'] = np.sum(df[ipp].values * ia, axis=1)\n    df['SeriesRank'] = df.groupby('SeriesID').DistanceProjection.rank('first').astype('int')\n    return df.sort_values(['SeriesID', 'SeriesRank']).set_index(['SeriesID', 'SeriesRank'])","metadata":{"execution":{"iopub.status.busy":"2023-09-20T10:31:23.552668Z","iopub.execute_input":"2023-09-20T10:31:23.553004Z","iopub.status.idle":"2023-09-20T10:31:23.56086Z","shell.execute_reply.started":"2023-09-20T10:31:23.552976Z","shell.execute_reply":"2023-09-20T10:31:23.559734Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from typing import NamedTuple, List, Tuple, Any\n\n\nclass SeriesMetadata(NamedTuple):\n    origin: Tuple[float, float, float]\n    spacing: Tuple[float, float, float]\n    direction: Tuple[float, float, float, float, float, float]\n\n    @property\n    def sitk_direction(self):\n        # SITK expects directions as 3x3 matrix in row-major order.\n        a_2 = (1 - self.direction[0] ** 2 - self.direction[3] ** 2)\n        a_2 = max(a_2, 0.0)\n        a_5 = (1 - self.direction[1] ** 2 - self.direction[4] ** 2)\n        a_5 = max(a_5, 0.0)\n        a_8 = (1 - self.direction[2] ** 2 - self.direction[5] ** 2)\n        a_8 = max(a_8, 0.0)\n        return [\n            self.direction[0],\n            self.direction[3],\n            a_2,\n            self.direction[1],\n            self.direction[4],\n            a_5,\n            self.direction[2],\n            self.direction[5],\n            a_8,\n        ]","metadata":{"execution":{"iopub.status.busy":"2023-09-20T10:31:23.56237Z","iopub.execute_input":"2023-09-20T10:31:23.562964Z","iopub.status.idle":"2023-09-20T10:31:23.577446Z","shell.execute_reply.started":"2023-09-20T10:31:23.562932Z","shell.execute_reply":"2023-09-20T10:31:23.576427Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import functools\n\n\nclass SeriesNotExistsError(Exception):\n    pass\n\n\ndef load_dcm(path, bit_shift):\n    try:\n        dcm = pydicom.dcmread(path)\n        img = dcm.pixel_array\n        if bit_shift is not None:\n            dtype = img.dtype \n            img = (img << bit_shift).astype(dtype) >> bit_shift\n        img = pydicom.pixel_data_handlers.apply_rescale(img, dcm)\n        img = pydicom.pixel_data_handlers.apply_windowing(img, dcm)\n        return img.T\n    except Exception:\n        return None\n\nload_dcm_delayed = joblib.delayed(load_dcm) \n\n\ndef load_series(tags, series_id, n_jobs=-1):\n\n    try:\n        df = tags.loc[series_id]\n    except KeyError:\n        raise SeriesNotExistsError(f'{series_id}')\n\n    metadata = df.iloc[0]\n\n    if metadata.PixelRepresentation == 1:\n        bit_shift = metadata.BitsAllocated - metadata.BitsStored\n    else:\n        bit_shift = None\n\n    if n_jobs is not None:\n        load_dcm_ = functools.partial(load_dcm_delayed, bit_shift=bit_shift)\n        imgs = joblib.Parallel(n_jobs=n_jobs)(load_dcm_(path) for path in df.path.values)\n    else:\n        load_dcm_ = functools.partial(load_dcm, bit_shift=bit_shift)\n        imgs = [load_dcm_(path) for path in df.path.values]\n\n    fill_value = None\n    failed_ixs = set()\n    for i, img in enumerate(imgs):\n        if fill_value is not None and failed_ixs:\n            break\n        elif img is None:\n            failed_ixs.add(i)\n        elif fill_value is None:\n            fill_value = np.zeros_like(img)\n\n    if failed_ixs and fill_value is None:\n        raise SeriesNotExistsError(f'All files found loading errors in series_id: {series_id}.')\n    elif failed_ixs:\n        for i in failed_ixs:\n            imgs[i] = fill_value\n    \n    img = np.c_[imgs]\n    img = img.transpose((1, 2, 0))\n\n    img -= img.min()\n    img *= CONFIG.series_dtype_out_max / img.max()\n    img = img.astype(CONFIG.series_dtype_out)\n\n    if len(df) > 1:\n        dz = abs(df.iloc[1].ImagePositionPatient2 - metadata.ImagePositionPatient2)\n    else:\n        dz = metadata.SliceThickness\n\n    spacing = [float(metadata.PixelSpacing0), float(metadata.PixelSpacing1), float(dz)]\n    origin = [float(metadata[f'ImagePositionPatient{i}']) for i in range(3)]\n    direction = [float(metadata[f'ImageOrientationPatient{i}']) for i in range(6)]\n    metadata = SeriesMetadata(\n        spacing=spacing,\n        origin=origin,\n        direction=direction,\n    )\n    return img, metadata","metadata":{"execution":{"iopub.status.busy":"2023-09-20T10:31:23.580837Z","iopub.execute_input":"2023-09-20T10:31:23.581227Z","iopub.status.idle":"2023-09-20T10:31:23.600536Z","shell.execute_reply.started":"2023-09-20T10:31:23.581197Z","shell.execute_reply":"2023-09-20T10:31:23.599196Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import SimpleITK as sitk\n\nTARGET_SPACING = np.round(np.divide(CONFIG.roi, CONFIG.shape), 1).tolist()\n\n\ndef resample_image_sitk(image, interpolator):\n    factor = np.divide(image.GetSpacing(), TARGET_SPACING)\n    factor_size = np.asarray(image.GetSize() * factor, dtype=float)\n    new_size = np.max([factor_size, CONFIG.shape], axis=0).astype(int).tolist()\n    return sitk.Resample(\n        image1=image,\n        size=new_size,\n        transform=sitk.Transform(),\n        interpolator=interpolator,\n        outputOrigin=image.GetOrigin(),\n        outputSpacing=TARGET_SPACING,\n        outputDirection=image.GetDirection(),\n        defaultPixelValue=0,\n        outputPixelType=image.GetPixelID(),\n    )\n\n\ndef transform_img(img, metadata, interpolator):\n    img_sitk = sitk.GetImageFromArray(img.T)\n    img_sitk.SetOrigin(metadata.origin)\n    img_sitk.SetDirection(metadata.sitk_direction)\n    img_sitk.SetSpacing(metadata.spacing)\n\n    img_sitk_rs = resample_image_sitk(img_sitk, interpolator)\n    img_rs = sitk.GetArrayFromImage(img_sitk_rs).T\n\n    start = np.subtract(img_rs.shape, CONFIG.shape) // 2\n    end = start + CONFIG.shape\n    return img_rs[start[0]:end[0], start[1]:end[1], start[2]:end[2]]","metadata":{"execution":{"iopub.status.busy":"2023-09-20T10:31:23.602048Z","iopub.execute_input":"2023-09-20T10:31:23.602421Z","iopub.status.idle":"2023-09-20T10:31:24.047105Z","shell.execute_reply.started":"2023-09-20T10:31:23.602367Z","shell.execute_reply":"2023-09-20T10:31:24.045948Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def load_segmentations(series_id):\n    seg = nib.load(f'/kaggle/input/rsna-2023-abdominal-trauma-detection/segmentations/{series_id}.nii')\n    seg = nib.as_closest_canonical(seg) # Ensures output is RAS, DICOM is LPS.\n    return seg.get_fdata().astype(CONFIG.segmentation_dtype_out)[::-1, ::-1, :]","metadata":{"execution":{"iopub.status.busy":"2023-09-20T10:31:24.050056Z","iopub.execute_input":"2023-09-20T10:31:24.050416Z","iopub.status.idle":"2023-09-20T10:31:24.056017Z","shell.execute_reply.started":"2023-09-20T10:31:24.050385Z","shell.execute_reply":"2023-09-20T10:31:24.05525Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from multiprocessing import Pool, cpu_count\n\nimport gc\nimport tqdm\nimport joblib\nimport functools\nfrom tqdm.contrib.concurrent import process_map\n\n\ndef transform_series(series_id, tags, n_jobs=-1):\n    img, metadata = load_series(tags, series_id, n_jobs)\n    img_t = transform_img(img, metadata, interpolator=sitk.sitkLinear)\n    img_path = os.path.join(FILESYSTEM.local_series, f'{series_id}.npz')\n    with open(img_path, 'wb') as f:\n        np.save(f, img_t)\n    # Normalise raw segmentations if exist.\n    if series_id in SERIES_INFO.train_segmentations:\n        img_segs = load_segmentations(series_id)\n        img_segs_t = transform_img(img_segs, metadata, interpolator=sitk.sitkNearestNeighbor)\n        img_segs_path = os.path.join(FILESYSTEM.local_segmentations, f'{series_id}.npz')\n        with open(img_segs_path, 'wb') as f:\n            np.save(f, img_segs_t)\n\n    \ndef transform_data(meta, tags, max_n=None, n_jobs=None, distribute_series_wise=True, chunksize=1):\n    os.makedirs(FILESYSTEM.local_series, exist_ok=True)\n    os.makedirs(FILESYSTEM.local_segmentations, exist_ok=True)\n    \n    tags = sort_tags(tags)\n\n    series_ids = meta.series_id.values\n    if max_n is not None:\n        series_ids = series_ids[:max_n]\n    if n_jobs is None or (n_jobs is not None and not distribute_series_wise):\n        for series_id in tqdm.tqdm(series_ids):\n            transform_series(series_id, tags, n_jobs=n_jobs)\n    else:\n        if n_jobs == -1:\n            n_jobs = cpu_count()\n        f = functools.partial(transform_series, tags=tags, n_jobs=None)\n        r = process_map(f, series_ids, tqdm_class=tqdm.tqdm, max_workers=n_jobs, chunksize=chunksize)","metadata":{"execution":{"iopub.status.busy":"2023-09-20T10:31:24.057061Z","iopub.execute_input":"2023-09-20T10:31:24.057425Z","iopub.status.idle":"2023-09-20T10:31:24.181551Z","shell.execute_reply.started":"2023-09-20T10:31:24.057397Z","shell.execute_reply":"2023-09-20T10:31:24.180382Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"transform_data(train_meta, train_tags, n_jobs=cpu_count() + 4, distribute_series_wise=True, chunksize=5)","metadata":{"execution":{"iopub.status.busy":"2023-09-20T10:31:24.183604Z","iopub.execute_input":"2023-09-20T10:31:24.184046Z","iopub.status.idle":"2023-09-20T10:36:20.611174Z","shell.execute_reply.started":"2023-09-20T10:31:24.184004Z","shell.execute_reply":"2023-09-20T10:36:20.609827Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"CONFIG.dump()\n!mv {FILESYSTEM.local_series} ./\n!mv {FILESYSTEM.local_segmentations} ./\n!rm -r {FILESYSTEM.local}","metadata":{"execution":{"iopub.status.busy":"2023-09-20T10:36:20.613553Z","iopub.execute_input":"2023-09-20T10:36:20.614297Z","iopub.status.idle":"2023-09-20T10:36:23.991911Z","shell.execute_reply.started":"2023-09-20T10:36:20.614244Z","shell.execute_reply":"2023-09-20T10:36:23.990405Z"},"trusted":true},"execution_count":null,"outputs":[]}]}