{"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":"![](https://raw.githubusercontent.com/fepegar/torchio/main/docs/source/favicon_io/for_readme_2000x462.png)\n\n[TorchIO](http://torchio.org/) is a Python library for loading, preprocessing, augmentation and sampling for multidimensional medical images in deep learning.\n\nWe can leverage the new [`RSNACervicalSpineFracture`](https://torchio.readthedocs.io/datasets.html#rsnacervicalspinefracture) class to avoid dealing with the complex [DICOM](https://www.dicomstandard.org/) and [NIfTI](https://nifti.nimh.nih.gov/) formats for medical images.","metadata":{}},{"cell_type":"code","source":"%pip install --quiet torchio","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2022-08-02T23:55:47.69815Z","iopub.execute_input":"2022-08-02T23:55:47.698936Z","iopub.status.idle":"2022-08-02T23:56:03.417322Z","shell.execute_reply.started":"2022-08-02T23:55:47.698838Z","shell.execute_reply":"2022-08-02T23:56:03.415772Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import copy\n\nimport torchio as tio\nimport torch\nfrom tqdm import tqdm\nimport matplotlib as mpl\nimport matplotlib.pyplot as plt\n\ntorch.manual_seed(20220802)\nmpl.rcParams['figure.figsize'] = 12, 8","metadata":{"execution":{"iopub.status.busy":"2022-08-02T23:56:03.42022Z","iopub.execute_input":"2022-08-02T23:56:03.420671Z","iopub.status.idle":"2022-08-02T23:56:06.363251Z","shell.execute_reply.started":"2022-08-02T23:56:03.42063Z","shell.execute_reply":"2022-08-02T23:56:06.361573Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Dataset\n\nLet's create an instance of `RSNACervicalSpineFracture`, which inherits directly from [`SubjectsDataset`](https://torchio.readthedocs.io/data/dataset.html) and indirectly from [`torch.utils.data.Dataset`](https://pytorch.org/docs/stable/data.html#torch.utils.data.Dataset).","metadata":{}},{"cell_type":"code","source":"root_dir = '/kaggle/input/rsna-2022-cervical-spine-fracture-detection'\ndataset = tio.datasets.RSNACervicalSpineFracture(root_dir, add_segmentations=True, add_bounding_boxes=True)\nprint('Number of subjects:', len(dataset))","metadata":{"execution":{"iopub.status.busy":"2022-08-02T23:56:06.36596Z","iopub.execute_input":"2022-08-02T23:56:06.367793Z","iopub.status.idle":"2022-08-02T23:56:09.480991Z","shell.execute_reply.started":"2022-08-02T23:56:06.367719Z","shell.execute_reply":"2022-08-02T23:56:09.479779Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Subject\n\nLet's inspect the first [`Subject`](https://torchio.readthedocs.io/data/subject.html) in the dataset.","metadata":{}},{"cell_type":"code","source":"first_subject = dataset[0]\nfirst_subject","metadata":{"execution":{"iopub.status.busy":"2022-08-02T23:56:09.48376Z","iopub.execute_input":"2022-08-02T23:56:09.484258Z","iopub.status.idle":"2022-08-02T23:56:15.327811Z","shell.execute_reply.started":"2022-08-02T23:56:09.484224Z","shell.execute_reply":"2022-08-02T23:56:15.326595Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"for key, value in first_subject.items():\n    print(f'{key}: {value}')","metadata":{"execution":{"iopub.status.busy":"2022-08-02T23:56:15.329411Z","iopub.execute_input":"2022-08-02T23:56:15.329901Z","iopub.status.idle":"2022-08-02T23:56:15.353785Z","shell.execute_reply.started":"2022-08-02T23:56:15.329855Z","shell.execute_reply":"2022-08-02T23:56:15.352097Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We can see that the patient's C1 and C2 vertebrae are fractured, according to the labels. We can also see that the CT scan is an instance of [`ScalarImage`](https://torchio.readthedocs.io/data/image.html#torchio.ScalarImage) with $512 \\times 512 \\times 243$ voxels, in-plane pixel spacing of 0.44 mm, slice thickness of 0.80 mm, [LPS+ orientation](https://nipy.org/nibabel/image_orientation.html) data type \"short tensor\" (or signed 16-bit integers) and takes 121.5 MiB of memory.\n\nLet's take a look at the 3D image in the subject.","metadata":{}},{"cell_type":"code","source":"first_subject.plot()","metadata":{"execution":{"iopub.status.busy":"2022-08-02T23:56:15.355292Z","iopub.execute_input":"2022-08-02T23:56:15.35572Z","iopub.status.idle":"2022-08-02T23:56:16.526794Z","shell.execute_reply.started":"2022-08-02T23:56:15.355681Z","shell.execute_reply":"2022-08-02T23:56:16.525715Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The [CT scan](https://www.nhs.uk/conditions/ct-scan/) looks black where there is no reconstruction data, and the contrast in the region of interest is not great. Let's look at the intensity distribution.","metadata":{}},{"cell_type":"code","source":"first_subject.ct.hist()","metadata":{"execution":{"iopub.status.busy":"2022-08-02T23:56:16.527966Z","iopub.execute_input":"2022-08-02T23:56:16.528563Z","iopub.status.idle":"2022-08-02T23:56:23.6089Z","shell.execute_reply.started":"2022-08-02T23:56:16.52852Z","shell.execute_reply":"2022-08-02T23:56:23.607581Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Intensity preprocessing\n\nWe can look at the [Hounsfield scale](https://en.wikipedia.org/wiki/Hounsfield_scale) to keep only intensity values within a sensible range using the [`Clamp`](https://torchio.readthedocs.io/transforms/preprocessing.html#clamp) transform.","metadata":{}},{"cell_type":"code","source":"HOUNSFIELD_AIR, HOUNSFIELD_BONE = -1000, 1900\nclamp = tio.Clamp(out_min=HOUNSFIELD_AIR, out_max=HOUNSFIELD_BONE)\nfirst_subject_clamped = clamp(first_subject)\nfirst_subject_clamped.ct.hist()\nfirst_subject_clamped.plot()","metadata":{"execution":{"iopub.status.busy":"2022-08-02T23:56:23.610946Z","iopub.execute_input":"2022-08-02T23:56:23.611708Z","iopub.status.idle":"2022-08-02T23:56:28.020366Z","shell.execute_reply.started":"2022-08-02T23:56:23.611655Z","shell.execute_reply":"2022-08-02T23:56:28.018968Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"That looks much better! Let's now normalize the values to [0, 1], a common practice in neural network training. We will also saturate values in the first and last half percentiles, to remove potential outliers. Clamping and rescaling can be concatenated using [`Compose`](https://torchio.readthedocs.io/transforms/augmentation.html#compose) to create an intensity preprocessing transform.","metadata":{}},{"cell_type":"code","source":"rescale = tio.RescaleIntensity(percentiles=(0.5, 99.5))\npreprocess_intensity = tio.Compose([\n    clamp,\n    rescale,\n])","metadata":{"execution":{"iopub.status.busy":"2022-08-02T23:56:28.022017Z","iopub.execute_input":"2022-08-02T23:56:28.022422Z","iopub.status.idle":"2022-08-02T23:56:28.029793Z","shell.execute_reply.started":"2022-08-02T23:56:28.022388Z","shell.execute_reply":"2022-08-02T23:56:28.028209Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"first_subject_preprocessed = preprocess_intensity(first_subject)\nfirst_subject_preprocessed.ct.hist(show=False), plt.ylim(0, 1e6)\nfirst_subject_preprocessed.plot()","metadata":{"execution":{"iopub.status.busy":"2022-08-02T23:56:28.034209Z","iopub.execute_input":"2022-08-02T23:56:28.035279Z","iopub.status.idle":"2022-08-02T23:56:34.310006Z","shell.execute_reply.started":"2022-08-02T23:56:28.035231Z","shell.execute_reply":"2022-08-02T23:56:34.308737Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Let's look for a subject with an associated segmentation. We will use the [`dry_iter` method of `SubjectsDataset`](https://torchio.readthedocs.io/data/dataset.html#torchio.data.SubjectsDataset.dry_iter) as we don't want to load any data here.","metadata":{}},{"cell_type":"code","source":"for subject in dataset.dry_iter():\n    if 'seg' in subject:\n        subject_with_seg = subject\n        break\nsubject_with_seg.plot(reorient=False)","metadata":{"execution":{"iopub.status.busy":"2022-08-02T23:56:34.311984Z","iopub.execute_input":"2022-08-02T23:56:34.312958Z","iopub.status.idle":"2022-08-02T23:56:42.261408Z","shell.execute_reply.started":"2022-08-02T23:56:34.312907Z","shell.execute_reply":"2022-08-02T23:56:42.26Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Spatial preprocessing\n\nThe CT scan and the corresponding segmentation are not aligned! Let's look at their orientations.","metadata":{}},{"cell_type":"code","source":"print('CT orientation:', subject_with_seg.ct.orientation)\nprint('Segmentation orientation:', subject_with_seg.seg.orientation)","metadata":{"execution":{"iopub.status.busy":"2022-08-02T23:56:42.262866Z","iopub.execute_input":"2022-08-02T23:56:42.263261Z","iopub.status.idle":"2022-08-02T23:56:42.273749Z","shell.execute_reply.started":"2022-08-02T23:56:42.263228Z","shell.execute_reply":"2022-08-02T23:56:42.272143Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"It looks like voxels along the second spatial dimension grow towards the back (Posterior) in the CT and towards the front (Anterior) in the segmentation. That's why the images looked flipped with respect to the [coronal](https://en.wikipedia.org/wiki/Coronal_plane) plane.\n\nWe can normalize the orientation (to RAS+) using [`ToCanonical`](https://torchio.readthedocs.io/transforms/preprocessing.html#tocanonical). Also, as the resolution is quite high, we will downsample to a sensible value (1 mm isotropic) for faster computations using [`Resample`](https://torchio.readthedocs.io/transforms/preprocessing.html#resample). These two transforms are implemented on top of [NiBabel](https://nipy.org/nibabel/) and [SimpleITK](https://simpleitk.org/), respectively.\n\nAs before, we use `Compose` to concatenate preprocessing transforms.","metadata":{}},{"cell_type":"code","source":"normalize_orientation = tio.ToCanonical()\ndownsample = tio.Resample(1)\npreprocess_spatial = tio.Compose([\n    normalize_orientation,\n    downsample,\n])\npreprocess = tio.Compose([\n    preprocess_intensity,\n    preprocess_spatial,\n])","metadata":{"execution":{"iopub.status.busy":"2022-08-02T23:56:42.276026Z","iopub.execute_input":"2022-08-02T23:56:42.277474Z","iopub.status.idle":"2022-08-02T23:56:42.285982Z","shell.execute_reply.started":"2022-08-02T23:56:42.27742Z","shell.execute_reply":"2022-08-02T23:56:42.284419Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"subject_with_seg_preprocessed = preprocess(subject_with_seg)\nsubject_with_seg_preprocessed.plot(reorient=False)","metadata":{"execution":{"iopub.status.busy":"2022-08-02T23:56:42.287856Z","iopub.execute_input":"2022-08-02T23:56:42.288298Z","iopub.status.idle":"2022-08-02T23:56:50.243202Z","shell.execute_reply.started":"2022-08-02T23:56:42.288261Z","shell.execute_reply":"2022-08-02T23:56:50.24215Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"That's better! We actually didn't need to pass `reorient=False` to the `plot` method, as all that does is apply `ToCanonical` internally, which we have already done.","metadata":{}},{"cell_type":"markdown","source":"## Data augmentation\n\nThere are also many [augmentation transforms](https://torchio.readthedocs.io/transforms/augmentation.html) available in TorchIO. Let's compose some appropriate ones. If we were using MRI, we should definitely leverage some of the [MRI k-space artifact simulation transforms](https://torchio.readthedocs.io/transforms/augmentation.html#randommotion). Here, we will use simpler ones.","metadata":{}},{"cell_type":"code","source":"augment = tio.Compose([\n    tio.RandomAnisotropy(p=0.25),\n    tio.RandomAffine(),\n    tio.RandomFlip(),\n    tio.RandomNoise(p=0.25),\n    tio.RandomGamma(p=0.5),\n])","metadata":{"execution":{"iopub.status.busy":"2022-08-02T23:56:50.244611Z","iopub.execute_input":"2022-08-02T23:56:50.244934Z","iopub.status.idle":"2022-08-02T23:56:50.251685Z","shell.execute_reply.started":"2022-08-02T23:56:50.244906Z","shell.execute_reply":"2022-08-02T23:56:50.250387Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"In this case, we always apply an affine (anisotropic scaling + rotation) transform. The others are applied with a probability of 0.25 or 0.5.\n\nInternally, [`RandomFlip`](https://torchio.readthedocs.io/transforms/augmentation.html#randomflip) uses 0.5). By default, flipping is applied only around the sagittal plane, i.e., left-right. But if you feel adventurous, feel free to explore other axes as well!\n\nLet's set our training and validation transforms. We want augmentation during training only; our validation transform is simply our deterministic preprocessing.","metadata":{}},{"cell_type":"code","source":"train_transform = tio.Compose([\n    preprocess,\n    augment,\n])\nval_transform = preprocess","metadata":{"execution":{"iopub.status.busy":"2022-08-02T23:56:50.253317Z","iopub.execute_input":"2022-08-02T23:56:50.253752Z","iopub.status.idle":"2022-08-02T23:56:50.262549Z","shell.execute_reply.started":"2022-08-02T23:56:50.253672Z","shell.execute_reply":"2022-08-02T23:56:50.261205Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Let's look at some preprocessed and randomly augmented variations of our subject.","metadata":{}},{"cell_type":"code","source":"for _ in range(5):\n    subject_with_seg_augmented = train_transform(subject_with_seg)\n    subject_with_seg_augmented.plot()","metadata":{"execution":{"iopub.status.busy":"2022-08-02T23:56:50.264319Z","iopub.execute_input":"2022-08-02T23:56:50.265705Z","iopub.status.idle":"2022-08-02T23:57:32.326816Z","shell.execute_reply.started":"2022-08-02T23:56:50.26565Z","shell.execute_reply":"2022-08-02T23:57:32.325828Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Intensity transforms are only applied to instances of `ScalarImage`, whereas spatial transforms are applied also to instances of [`LabelMap`](https://torchio.readthedocs.io/data/image.html#torchio.LabelMap). Of course, the same random parameters for spatial transforms are applied to all images within the same subject!","metadata":{}},{"cell_type":"markdown","source":"## Training and validation splits","metadata":{}},{"cell_type":"markdown","source":"Let's use 80% of our data for training and 20% for validation.","metadata":{}},{"cell_type":"code","source":"train_to_val_ratio = 0.8\nnum_subjects = len(dataset)\nnum_train_subjects = int(train_to_val_ratio * num_subjects)\nnum_val_subjects = num_subjects - num_train_subjects\nprint('Number of subjects for training:  ', num_train_subjects)\nprint('Number of subjects for validation: ', num_val_subjects)","metadata":{"execution":{"iopub.status.busy":"2022-08-02T23:57:32.328418Z","iopub.execute_input":"2022-08-02T23:57:32.329127Z","iopub.status.idle":"2022-08-02T23:57:32.335439Z","shell.execute_reply.started":"2022-08-02T23:57:32.32909Z","shell.execute_reply":"2022-08-02T23:57:32.334513Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We will assume we won't need the segmentations or bounding boxes during training. We create a new instance of `RSNACervicalSpineFracture` and random split it between training and validation sets using PyTorch.","metadata":{}},{"cell_type":"code","source":"no_segs_dataset = tio.datasets.RSNACervicalSpineFracture(root_dir)\ntrain_set, val_set = torch.utils.data.random_split(no_segs_dataset, [num_train_subjects, num_val_subjects])","metadata":{"execution":{"iopub.status.busy":"2022-08-02T23:57:32.336777Z","iopub.execute_input":"2022-08-02T23:57:32.337356Z","iopub.status.idle":"2022-08-02T23:57:34.051384Z","shell.execute_reply.started":"2022-08-02T23:57:32.337324Z","shell.execute_reply":"2022-08-02T23:57:34.050377Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We use the training and validation transforms compute above for our training and validation sets.","metadata":{}},{"cell_type":"code","source":"# Both subsets share the \"dataset\" attribute, so the last statement would\n# overwrite the previous one unless we deep-copy one of the subsets\nval_set = copy.deepcopy(val_set)\ntrain_set.dataset.set_transform(train_transform)\nval_set.dataset.set_transform(val_transform)","metadata":{"execution":{"iopub.status.busy":"2022-08-02T23:57:34.052435Z","iopub.execute_input":"2022-08-02T23:57:34.052845Z","iopub.status.idle":"2022-08-02T23:57:34.472123Z","shell.execute_reply.started":"2022-08-02T23:57:34.052809Z","shell.execute_reply":"2022-08-02T23:57:34.470587Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Let's look at some of the images in our training set. They are randomly transformed by our `train_transform`.","metadata":{}},{"cell_type":"code","source":"for i in range(5):\n    train_set[i].plot()","metadata":{"execution":{"iopub.status.busy":"2022-08-02T23:57:34.473608Z","iopub.execute_input":"2022-08-02T23:57:34.473976Z","iopub.status.idle":"2022-08-02T23:58:38.914326Z","shell.execute_reply.started":"2022-08-02T23:57:34.473945Z","shell.execute_reply":"2022-08-02T23:58:38.912971Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Let's now look at some of the images in the validation set, which are not augmented (only preprocessed).","metadata":{}},{"cell_type":"code","source":"for i in range(5):\n    val_set[i].plot()","metadata":{"execution":{"iopub.status.busy":"2022-08-02T23:58:38.915943Z","iopub.execute_input":"2022-08-02T23:58:38.916537Z","iopub.status.idle":"2022-08-02T23:59:54.222995Z","shell.execute_reply.started":"2022-08-02T23:58:38.916462Z","shell.execute_reply":"2022-08-02T23:59:54.221778Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Conclusion\n\nWe have used TorchIO to easily load, inspect, visualize, preprocess and augment complex DICOM and NIfTI 3D CT scans. TorchIO empowers researchers to iterate quickly and improve their models without worrying about (re)implementing preprocessing and augmentation for medical images.\n\n## Next\n\n- As these images are large and have different fields of view, it might be a good idea to use [patch-based training and inference](https://torchio.readthedocs.io/patches/index.html).\n- TorchIO can be used with other popular frameworks such as [MONAI and PyTorch Lightning](https://colab.research.google.com/github/fepegar/torchio-notebooks/blob/main/notebooks/Brain_parcellation_with_TorchIO_and_HighRes3DNet.ipynb).\n\n## Feedback\n\n- If you find this notebook useful, please upvote and comment!\n- Feel free to [request features on GitHub](https://github.com/fepegar/torchio/issues/new?assignees=&labels=enhancement&template=feature_request.md&title=) that might be useful for this challenge.","metadata":{}}]}