{"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":"code","source":"# This Python 3 environment comes with many helpful analytics libraries installed\n# It is defined by the kaggle/python Docker image: https://github.com/kaggle/docker-python\n# For example, here's several helpful packages to load\n\nimport numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\n\n# Input data files are available in the read-only \"../input/\" directory\n# For example, running this (by clicking run or pressing Shift+Enter) will list all files under the input directory\n\nimport os\nfor dirname, _, filenames in os.walk('/kaggle/input'):\n    for filename in filenames:\n        print(os.path.join(dirname, filename))\n\n# You can write up to 20GB to the current directory (/kaggle/working/) that gets preserved as output when you create a version using \"Save & Run All\" \n# You can also write temporary files to /kaggle/temp/, but they won't be saved outside of the current session","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","_kg_hide-input":true,"_kg_hide-output":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#NiBabel Neuroimaging in Python \n\nNiBabel is free-software (beer and speech) and covered by the MIT License.\n\n![](https://nilearn.github.io/stable/_static/nilearn-logo.png)https://nipy.org/nibabel/coordinate_systems.html\n\n\"Coordinate systems and affines\n\n\"A nibabel (and nipy) image is the association of three things:\n\nThe image data array: a 3D or 4D array of image data\n\nAn affine array that tells you the position of the image array data in a reference space.\n\nImage metadata (data about the data) describing the image, usually in the form of an image header.\n\n\"This document describes how the affine array describes the position of the image data in a reference space. On the way we will define what we mean by reference space, and the reference spaces that Nibabel uses.\"\n\nhttps://nipy.org/nibabel/coordinate_systems.html","metadata":{}},{"cell_type":"markdown","source":"#All script by Michael Beregov https://www.kaggle.com/code/boojum/connecting-voxel-spaces\n\n#Vote  Michael's (boojum) work!  ","metadata":{}},{"cell_type":"code","source":"#Code by Michael Beregov https://www.kaggle.com/code/boojum/connecting-voxel-spaces\n\nimport os\nimport sys \nfrom tqdm import tqdm\nimport numpy as np\nimport pandas as pd\nimport matplotlib.pyplot as plt\nfrom PIL import Image\nimport pydicom\n\nimport numpy as np\nimport nibabel as nib\nimport matplotlib.pyplot as plt\nimport SimpleITK as sitk\n\ntrain_path = '../input/rsna-2022-cervical-spine-fracture-detection/train_images/'","metadata":{"execution":{"iopub.status.busy":"2022-08-02T18:54:22.836499Z","iopub.execute_input":"2022-08-02T18:54:22.837712Z","iopub.status.idle":"2022-08-02T18:54:23.732421Z","shell.execute_reply.started":"2022-08-02T18:54:22.837552Z","shell.execute_reply":"2022-08-02T18:54:23.73114Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#Voxels\n\n\"In 3D computer graphics, a voxel represents a value on a regular grid in three-dimensional space. As with pixels in a 2D bitmap, voxels themselves do not typically have their position (i.e. coordinates) explicitly encoded with their values. Instead, rendering systems infer the position of a voxel based upon its position relative to other voxels (i.e., its position in the data structure that makes up a single volumetric image).\"\n\nhttps://en.wikipedia.org/wiki/Voxel","metadata":{}},{"cell_type":"markdown","source":"#How to Automate Voxel Modelling of 3D Point Cloud with Python\n\nBy Florent Poux, Ph.D.\n\n![](https://miro.medium.com/max/1400/1*aIWdWJJ_oxlagxFyrBkEsQ.gif)\nhttps://towardsdatascience.com/how-to-automate-voxel-modelling-of-3d-point-cloud-with-python-459f4d43a227","metadata":{}},{"cell_type":"code","source":"train_dirs = os.listdir(train_path)","metadata":{"execution":{"iopub.status.busy":"2022-08-02T18:54:36.087427Z","iopub.execute_input":"2022-08-02T18:54:36.087766Z","iopub.status.idle":"2022-08-02T18:54:36.225564Z","shell.execute_reply.started":"2022-08-02T18:54:36.08774Z","shell.execute_reply":"2022-08-02T18:54:36.22406Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#Voxel Size\n\nLast revised by Dr Daniel J Bell on 04 Sep 2018\n\n\"Voxel size is an important component of image quality. Voxel is the 3-D analog of a pixel. Voxel size is related to both the pixel size and slice thickness. Pixel size is dependent on both the field of view and the image matrix. The pixel size is equal to the field of view divided by the matrix size. The matrix size is typically 128x, 256x or 512x. Pixel size is typically between 0.5 and 1.5 mm. The smaller the pixel size, the greater the image spatial resolution.\"\n\n\"Increased voxel size results in an increased signal-to-noise ratio. The trade-off for increased voxel size is decreased spatial resolution. Voxel size can be influenced by receiver coil characteristics. For examples, surface coils indirectly improve resolution by enabling a smaller voxel size for the same signal-to-noise ratio.\"\n\n\"Voxel size can contribute to artifacts in MRI. Many MR artifacts are attributable to errors in the underlying spatial encoding of the radiofrequency signals arising from image voxels. Motion artefacts can occur in the phase-encoding direction because a specific tissue voxel may change location between acquisition cycles, leading to phase encoding errors. This manifests as a streak or ghost in the final image, and can be reduced with image gating and regional presaturation techniques.\"\n\nhttps://radiopaedia.org/articles/voxel-size-1","metadata":{}},{"cell_type":"markdown","source":"#Print the 5 first filenames","metadata":{}},{"cell_type":"code","source":"# Print out the first 5 file names to verify we're in the right folder.\nprint (f'Total of {len(train_dirs)} DICOM images.\\nFirst 5 filenames:' )\ntrain_dirs[:5]","metadata":{"execution":{"iopub.status.busy":"2022-08-02T18:54:41.015248Z","iopub.execute_input":"2022-08-02T18:54:41.015645Z","iopub.status.idle":"2022-08-02T18:54:41.025401Z","shell.execute_reply.started":"2022-08-02T18:54:41.015619Z","shell.execute_reply":"2022-08-02T18:54:41.024115Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#Effect of voxel size in CT simulations\n\nAuthors: Andrew L Goertzen, Freek J Beekman, Simon R Cherry - DOI:10.1109/NSSMIC.2000.949326\n\n\"In computer simulations of X-ray CT systems one can either use continuous geometrical descriptions for phantoms or a voxelized representation. The voxelized approach allows arbitrary phantoms to be defined without being confined to geometrical shapes.\"\n\n\"The disadvantage of the voxelized approach is that inherent errors are introduced due to the phantom voxelization. To study effects of phantom discretization, analytical CT simulations were run for a fan-beam geometry with phantom voxel sizes ranging from 0.0625 to 2 times the reconstructed pixel size and noise levels corresponding to 103 to 107 photons per detector pixel prior to attenuation. Differences in the filtered back-projection (FBP) images caused by different phantom matrix sizes were assessed by calculating the difference between reconstructions based on the finest matrix and coarser matrix simulations.\"\n\n\"In noise free simulations, all phantom matrix sizes produced a measurable difference from the almost continuous case. When even a small amount of noise was added to the projection data, the differences due to the phantom discretization were masked by the noise, and in all cases there was almost no improvement by using a phantom matrix that was more than twice as fine as the reconstruction matrix.\"\n\nhttps://www.researchgate.net/publication/3914486_Effect_of_voxel_size_in_CT_simulations","metadata":{}},{"cell_type":"code","source":"#Code by Michael Beregov https://www.kaggle.com/code/boojum/connecting-voxel-spaces\n\nplt.figure(figsize=(14,7))\nplt.subplot(121)\nplt.imshow(pydicom.dcmread(f'{train_path + train_dirs[0]}/1.dcm').pixel_array)\n\nplt.subplot(122)\nplt.imshow(pydicom.dcmread(f'{train_path + train_dirs[0]}/10.dcm').pixel_array);","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-08-02T18:57:23.785238Z","iopub.execute_input":"2022-08-02T18:57:23.785576Z","iopub.status.idle":"2022-08-02T18:57:24.247342Z","shell.execute_reply.started":"2022-08-02T18:57:23.78555Z","shell.execute_reply":"2022-08-02T18:57:24.244325Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#Affine matrix and simple resampling\n\n\"Each DICOM file stores information about its orientation in the scanner space (which is basically the real world, with the center of the coordinate system in the magnet isocenter).\"\n\n\"Let's convert one of the DICOM-series to a NIfTI file to see it most clearly. There are numerous ways to do it, we'll use the functionality of the SimpleITK library.\"\n\nBy Michael Beregov https://www.kaggle.com/code/boojum/connecting-voxel-spaces","metadata":{}},{"cell_type":"code","source":"#Code by Michael Beregov https://www.kaggle.com/code/boojum/connecting-voxel-spaces\n\nreader = sitk.ImageSeriesReader()\nreader.LoadPrivateTagsOn()","metadata":{"execution":{"iopub.status.busy":"2022-08-02T18:57:41.048922Z","iopub.execute_input":"2022-08-02T18:57:41.049247Z","iopub.status.idle":"2022-08-02T18:57:41.062865Z","shell.execute_reply.started":"2022-08-02T18:57:41.049223Z","shell.execute_reply":"2022-08-02T18:57:41.061629Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#Voxel-Wise Mapping and Deep Learning\n\nAutomatic detection and voxel-wise mapping of lumbar spine Modic changes with deep learning\n\nAuthors: Kenneth T. Gao, Radhika Tibrewala, Madeline Hess, Upasana U. Bharadwaj, Gaurav Inamdar, Thomas M. Link, Cynthia T. Chin, Valentina Pedoia, Sharmila Majumdar\n\n\"Modic changes (MCs) are the most prevalent classification system for describing magnetic resonance imaging (MRI) signal intensity changes in the vertebrae. However, there is a growing need for novel quantitative and standardized methods of characterizing these anomalies, particularly for lesions of transitional or mixed nature, due to the lack of conclusive evidence of their associations with low back pain. This retrospective imaging study aims to develop an interpretable deep learning-based detection tool for voxel-wise mapping of MCs.\"\n\n\"The model successfully identified the presence of changes in 85.7% of samples in the unseen test set with a sensitivity of 0.71 (±0.072), specificity of 0.95 (±0.022), and Cohen's kappa score of 0.63. In the AI-assisted experiment, the agreement between the junior radiologist and the senior neuroradiologist significantly improved from Cohen's kappa score of 0.52 to 0.58 (p < 0.05).\"\n\n\"That deep learning-based approach demonstrates substantial agreement with radiologists and may serve as a tool to improve inter-rater reliability in the assessment of MCs.\"\n\nhttps://onlinelibrary.wiley.com/doi/10.1002/jsp2.1204","metadata":{}},{"cell_type":"code","source":"#Code by Michael Beregov https://www.kaggle.com/code/boojum/connecting-voxel-spaces\n\nfilenamesDICOM = reader.GetGDCMSeriesFileNames(f'{train_path}/{train_dirs[2]}')\nreader.SetFileNames(filenamesDICOM)\nt1_sitk = reader.Execute()","metadata":{"execution":{"iopub.status.busy":"2022-08-02T18:58:25.933964Z","iopub.execute_input":"2022-08-02T18:58:25.934297Z","iopub.status.idle":"2022-08-02T18:58:33.410408Z","shell.execute_reply.started":"2022-08-02T18:58:25.93427Z","shell.execute_reply":"2022-08-02T18:58:33.408662Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#Voxel and Dentistry\n\n\"A voxel is the smallest 3D element of the volume,2 and is typically represented as a cube or a box, with height, width and depth. Just as 2D images are made of several pixels (represented as squares, with height and width) and the smaller the pixel the better the quality of the picture, the same concept applies to a 3D data volume. Each three-dimensional voxel represents specific x-ray absorption.\"\n\n\"The voxel size on CBCT images is isotropic, which means that all the sides are the same dimension with uniform resolution in all directions. In contrast, an MDCT voxel is in general nonisotropic meaning that one side of the voxel is different in dimension. This is considered an advantage of the CBCT because if a certain structure needs to be measured, the measurement will be exact in all the three orthogonal planes. There are different voxel sizes depending on the capabilities of each unit. The small field of view units may use a small voxel size of 0.076 mm, which enables visualization of very small changes to structures. Other voxel sizes available for CBCT units are variable, such as 0.2 mm, 0.3 mm, and 0.4 mm. It is important to note that the larger the voxel size, the less resolution the image will have and less capability to differentiate between small structures. The voxel size is dependent of the imaging objective and the size of the unit detector.\"\n\nhttps://www.dentalcare.com/en-us/ce-courses/ce531/voxel","metadata":{}},{"cell_type":"code","source":"sitk.WriteImage(t1_sitk,'t1.nii')#We don' t have T1","metadata":{"execution":{"iopub.status.busy":"2022-08-02T18:59:12.43236Z","iopub.execute_input":"2022-08-02T18:59:12.432692Z","iopub.status.idle":"2022-08-02T18:59:12.703647Z","shell.execute_reply.started":"2022-08-02T18:59:12.432668Z","shell.execute_reply":"2022-08-02T18:59:12.701959Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#Now we can load it with Nibabel. There's a lot of stuff in the nibabel.nifti1.Nifti1Image object, but the two essential things are voxel array and affine matrix.","metadata":{}},{"cell_type":"code","source":"#Code by Michael Beregov https://www.kaggle.com/code/boojum/connecting-voxel-spaces\n\nt1_nib = nib.load('t1.nii')\nt1_nib","metadata":{"execution":{"iopub.status.busy":"2022-08-02T18:59:47.686864Z","iopub.execute_input":"2022-08-02T18:59:47.687177Z","iopub.status.idle":"2022-08-02T18:59:47.697187Z","shell.execute_reply.started":"2022-08-02T18:59:47.687147Z","shell.execute_reply":"2022-08-02T18:59:47.69513Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"t1_nib_array = t1_nib.get_fdata() #the voxel array. I don't have T1,that's why the array isn't as expected\nt1_nib_array[:3]","metadata":{"_kg_hide-output":true,"execution":{"iopub.status.busy":"2022-08-02T19:00:50.213066Z","iopub.execute_input":"2022-08-02T19:00:50.213406Z","iopub.status.idle":"2022-08-02T19:00:50.353131Z","shell.execute_reply.started":"2022-08-02T19:00:50.213379Z","shell.execute_reply":"2022-08-02T19:00:50.351924Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#NiBabel and Voxels\n\nNiBabel: Access a cacophony of neuro-imaging file formats\n\n\"A nibabel (and nipy) image is the association of three things:\n\n\"The image data array: a 3D or 4D array of image data. An affine array that tells you the position of the image array data in a reference space. image metadata (data about the data) describing the image, usually in the form of an image header.\"\n\nThis document describes how the affine array describes the position of the image data in a reference space. On the way we will define what we mean by reference space, and the reference spaces that Nibabel uses.\n\nhttps://nipy.org/nibabel/image_orientation.html","metadata":{}},{"cell_type":"markdown","source":"#NiBabel Snippet:","metadata":{}},{"cell_type":"code","source":"#https://nipy.org/nibabel/image_orientation.html\n\nimport numpy as np\n\nimport nibabel as nib\n\naffine = np.eye(4) # identity affine\nvoxel_data = np.random.normal(size=(10, 11, 12))\nimg = nib.Nifti1Image(voxel_data, affine)","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#NiBabel Tutorials\n\nNiBabel - Neuroimaging in Python\n\nhttps://nipy.org/nibabel/tutorials.html\n\n#Introduction to Dicoms\n\nhttps://nipy.org/nibabel/dicom/dicom_intro.html","metadata":{}},{"cell_type":"code","source":"t1_nib_array.shape","metadata":{"execution":{"iopub.status.busy":"2022-08-02T19:04:27.616467Z","iopub.execute_input":"2022-08-02T19:04:27.617326Z","iopub.status.idle":"2022-08-02T19:04:27.628306Z","shell.execute_reply.started":"2022-08-02T19:04:27.617292Z","shell.execute_reply":"2022-08-02T19:04:27.62697Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plt.imshow(t1_nib_array[:,:,t1_nib_array.shape[2]//2]);","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-08-02T19:04:32.651139Z","iopub.execute_input":"2022-08-02T19:04:32.652515Z","iopub.status.idle":"2022-08-02T19:04:32.827115Z","shell.execute_reply.started":"2022-08-02T19:04:32.652466Z","shell.execute_reply":"2022-08-02T19:04:32.826393Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#Voxel array\n\nSo, this is the voxel array. Its orientation in scanner space is encoded in the affine matrix:","metadata":{}},{"cell_type":"code","source":"#Code by Michael Beregov https://www.kaggle.com/code/boojum/connecting-voxel-spaces\n\nnp.set_printoptions(precision=4, suppress=True)\nt1_nib.affine","metadata":{"_kg_hide-output":true,"execution":{"iopub.status.busy":"2022-08-02T19:14:15.650631Z","iopub.execute_input":"2022-08-02T19:14:15.651032Z","iopub.status.idle":"2022-08-02T19:14:15.659464Z","shell.execute_reply.started":"2022-08-02T19:14:15.651005Z","shell.execute_reply":"2022-08-02T19:14:15.658203Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"\"The first 3х3 part of the matrix provides information about rotation and scaling. The fourth column tells us about translation.\"\n\n\"So, if we want to know the coordinates of a voxel b=[2,5,10] in the scanner space, we can just calculate its dot product with upper left 3х3 corner of the affine matrix, and then sum the result with the translation vector 𝑡 (which is the forth column of the affine matrix).\"\n\n#Code by Michael Beregov https://www.kaggle.com/code/boojum/connecting-voxel-spaces","metadata":{}},{"cell_type":"markdown","source":"$$ 1) \\,\\,\\,\\,\\, \n    \\begin{bmatrix} \n        { A }_{ 00 } & { A }_{ 01 } & { A }_{ 02 }\\\\ \n        { A }_{ 10 } & { A }_{ 11 } & { A }_{ 12 }\\\\ \n        { A }_{ 20 } & { A }_{ 21 } & { A }_{ 22 }\\end{bmatrix} \\cdot\n   \\begin{bmatrix} \n        { b_0 } \\\\ { b_1 }  \\\\ { b_2 } \\end{bmatrix} =\n   \\begin{bmatrix}      \n        { A }_{ 00 } \\\\ { A }_{ 10 } \\\\ { A }_{ 20 } \\end{bmatrix} \\cdot { b_0 } + \n   \\begin{bmatrix}      \n        { A }_{ 01 } \\\\ { A }_{ 11 } \\\\ { A }_{ 21 } \\end{bmatrix} \\cdot { b_1 } +    \n   \\begin{bmatrix}      \n        { A }_{ 02 } \\\\ { A }_{ 12 } \\\\ { A }_{ 22 } \\end{bmatrix} \\cdot { b_2 } =    \n   \\begin{bmatrix}   \n   { A }_{ 00 } { b_0 } + { A }_{ 01 } { b_1 } + { A }_{ 02 } { b_2 } \\\\\n   { A }_{ 10 } { b_0 } + { A }_{ 11 } { b_1 } + { A }_{ 12 } { b_2 } \\\\\n   { A }_{ 20 } { b_0 } + { A }_{ 21 } { b_1 } + { A }_{ 22 } { b_2 } \\end{bmatrix} = \n   \\begin{bmatrix} \n        { x_0 } \\\\ { x_1 }  \\\\ { x_2 } \\end{bmatrix} $$ <br/>\n\n\n$$ 2) \\,\\,\\,\\,\\, \n    \\begin{bmatrix} \n        { x_0 } \\\\ { x_1 }  \\\\ { x_2 } \\end{bmatrix}  + \n    \\begin{bmatrix} \n        { t_0 } \\\\ { t_1 }  \\\\ { t_2 } \\end{bmatrix} = \n    \\begin{bmatrix} \n        { x_0+t_0 } \\\\ { x_1+t_1 }  \\\\ { x_2+t_2 } \\end{bmatrix} = \n    \\begin{bmatrix} \n        { a } \\\\ { b }  \\\\ { c } \\end{bmatrix}\n        $$       \n        ","metadata":{}},{"cell_type":"code","source":"t1_nib.affine[:3,:3] @ np.array([2,5,10]) + t1_nib.affine[:3,3]","metadata":{"execution":{"iopub.status.busy":"2022-08-02T19:14:38.516551Z","iopub.execute_input":"2022-08-02T19:14:38.516887Z","iopub.status.idle":"2022-08-02T19:14:38.525497Z","shell.execute_reply.started":"2022-08-02T19:14:38.516863Z","shell.execute_reply":"2022-08-02T19:14:38.524485Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#Affine matrix\n\nThe last row in the affine matrix is always the [0,0,0,1] like in the identity matrix, it's just there to make the matrix square so we could use it as a linear operator. Therefore we can also just do this:","metadata":{}},{"cell_type":"code","source":"t1_nib.affine @ np.array([2,5,10,1]) #(add 1 as a fourth coordinate)","metadata":{"execution":{"iopub.status.busy":"2022-08-02T19:14:43.908686Z","iopub.execute_input":"2022-08-02T19:14:43.909076Z","iopub.status.idle":"2022-08-02T19:14:43.917219Z","shell.execute_reply.started":"2022-08-02T19:14:43.909041Z","shell.execute_reply":"2022-08-02T19:14:43.916099Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#We don't have Flair now\n\nApparently, as all series in the study have their affines linked to the same scanner space, we can resample them all into the same voxel space.","metadata":{}},{"cell_type":"code","source":"#Code by Michael Beregov https://www.kaggle.com/code/boojum/connecting-voxel-spaces\n\nfilenamesDICOM = reader.GetGDCMSeriesFileNames(f'{train_path}/{train_dirs[2]}')\nreader.SetFileNames(filenamesDICOM)\nflair_sitk = reader.Execute()\nsitk.WriteImage(flair_sitk,'flair.nii')\n\nflair_nib = nib.load('flair.nii')#We don't have FLAIR!\nflair_nib_array = flair_nib.get_fdata()","metadata":{"execution":{"iopub.status.busy":"2022-08-02T19:14:53.992492Z","iopub.execute_input":"2022-08-02T19:14:53.993515Z","iopub.status.idle":"2022-08-02T19:14:57.259776Z","shell.execute_reply.started":"2022-08-02T19:14:53.993485Z","shell.execute_reply":"2022-08-02T19:14:57.258464Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#Code by Michael Beregov https://www.kaggle.com/code/boojum/connecting-voxel-spaces\n\nplt.figure(figsize=(12,6))\nplt.subplot(121)\nplt.imshow(t1_nib_array[:,:,t1_nib_array.shape[2]//2])\nplt.subplot(122)\nplt.imshow(flair_nib_array[:,:,flair_nib_array.shape[2]//2]);","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-08-02T19:15:05.543533Z","iopub.execute_input":"2022-08-02T19:15:05.543891Z","iopub.status.idle":"2022-08-02T19:15:05.941645Z","shell.execute_reply.started":"2022-08-02T19:15:05.543865Z","shell.execute_reply":"2022-08-02T19:15:05.938953Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from nilearn.image import resample_img","metadata":{"execution":{"iopub.status.busy":"2022-08-02T19:15:35.319963Z","iopub.execute_input":"2022-08-02T19:15:35.320417Z","iopub.status.idle":"2022-08-02T19:15:36.842655Z","shell.execute_reply.started":"2022-08-02T19:15:35.320389Z","shell.execute_reply":"2022-08-02T19:15:36.841282Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time \nflair_resampled = resample_img(flair_nib, target_affine=t1_nib.affine, target_shape=t1_nib.shape)\nflair_resampled_array = flair_resampled.get_fdata()","metadata":{"execution":{"iopub.status.busy":"2022-08-02T19:15:41.32016Z","iopub.execute_input":"2022-08-02T19:15:41.320577Z","iopub.status.idle":"2022-08-02T19:15:41.329161Z","shell.execute_reply.started":"2022-08-02T19:15:41.320551Z","shell.execute_reply":"2022-08-02T19:15:41.32811Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#Code by Michael Beregov https://www.kaggle.com/code/boojum/connecting-voxel-spaces\n\nplt.figure(figsize=(12,6))\nplt.subplot(121)\nplt.imshow(t1_nib_array[:,:,t1_nib_array.shape[2]//2])\nplt.subplot(122)\nplt.imshow(flair_resampled_array[:,:,flair_resampled_array.shape[2]//2]);","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-08-02T19:15:48.224209Z","iopub.execute_input":"2022-08-02T19:15:48.224913Z","iopub.status.idle":"2022-08-02T19:15:48.74943Z","shell.execute_reply.started":"2022-08-02T19:15:48.224876Z","shell.execute_reply":"2022-08-02T19:15:48.748517Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#Above, since we don't Have Flair it didn't change anything.","metadata":{}},{"cell_type":"markdown","source":"#SimpleITK.SimpleITK.Image\n\nSo, the SimpleITK.SimpleITK.Image also contains information about voxel space orientation in the real world. It's not stored by the means of affine matrix though. Here, it goes in a few pieces:\n\nBy Michael Beregov https://www.kaggle.com/code/boojum/connecting-voxel-spaces","metadata":{}},{"cell_type":"code","source":"t1_sitk.GetOrigin() # which is a translation column from the affine matrix, but negative","metadata":{"execution":{"iopub.status.busy":"2022-08-02T19:16:03.878226Z","iopub.execute_input":"2022-08-02T19:16:03.879154Z","iopub.status.idle":"2022-08-02T19:16:03.89841Z","shell.execute_reply.started":"2022-08-02T19:16:03.879112Z","shell.execute_reply":"2022-08-02T19:16:03.897392Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"t1_sitk.GetSpacing() # which is how far away the voxel centers are from one another along each of the axes","metadata":{"execution":{"iopub.status.busy":"2022-08-02T19:16:09.837687Z","iopub.execute_input":"2022-08-02T19:16:09.838032Z","iopub.status.idle":"2022-08-02T19:16:09.847725Z","shell.execute_reply.started":"2022-08-02T19:16:09.838005Z","shell.execute_reply":"2022-08-02T19:16:09.846109Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"t1_sitk.GetDirection() # a flatten cosine matrix which shows rotation of voxel space axes relative to scanner space","metadata":{"execution":{"iopub.status.busy":"2022-08-02T19:16:14.840053Z","iopub.execute_input":"2022-08-02T19:16:14.842033Z","iopub.status.idle":"2022-08-02T19:16:14.852582Z","shell.execute_reply.started":"2022-08-02T19:16:14.841982Z","shell.execute_reply":"2022-08-02T19:16:14.850184Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#Spacing (Scaling)\n\nWe could actually extract all that information from the previously seen affine matrix. For example, knowing that each column affects one of the resulting voxel coordinates, we could get information about scaling (aka spacing, in this case).\n\nBy Michael Beregov https://www.kaggle.com/code/boojum/connecting-voxel-spaces","metadata":{}},{"cell_type":"code","source":"#Code by Michael Beregov https://www.kaggle.com/code/boojum/connecting-voxel-spaces\n\nx_scale = np.linalg.norm(t1_nib.affine[:,0])\ny_scale = np.linalg.norm(t1_nib.affine[:,1])\nz_scale = np.linalg.norm(t1_nib.affine[:,2])\nprint(x_scale, y_scale, z_scale)","metadata":{"execution":{"iopub.status.busy":"2022-08-02T19:16:27.059434Z","iopub.execute_input":"2022-08-02T19:16:27.05981Z","iopub.status.idle":"2022-08-02T19:16:27.067587Z","shell.execute_reply.started":"2022-08-02T19:16:27.059782Z","shell.execute_reply":"2022-08-02T19:16:27.066352Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#Cosine Matrix\n\nAlso, if we divide each column by the corresponding spacing number, we'll get a cosine matrix.","metadata":{}},{"cell_type":"code","source":"#Code by Michael Beregov https://www.kaggle.com/code/boojum/connecting-voxel-spaces\n\nt1_pure_rotation = np.hstack((t1_nib.affine[:,0].reshape(-1,1)/x_scale,\n                   t1_nib.affine[:,1].reshape(-1,1)/y_scale,\n                   t1_nib.affine[:,2].reshape(-1,1)/z_scale,\n                   t1_nib.affine[:,3].reshape(-1,1)))\nt1_pure_rotation[:3,:3]","metadata":{"execution":{"iopub.status.busy":"2022-08-02T19:16:34.475438Z","iopub.execute_input":"2022-08-02T19:16:34.476077Z","iopub.status.idle":"2022-08-02T19:16:34.483636Z","shell.execute_reply.started":"2022-08-02T19:16:34.476048Z","shell.execute_reply":"2022-08-02T19:16:34.482793Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#SimpleITK Cosine Matrix\n\n\"The first two rows differ in sign from the SimpleITK cosine matrix, I'm not sure why. It has something to do with a rotation direction.\"\n\n\"Btw, this information can be acquired from DICOM files. There's a tag for this:\"\n\nBy Michael Beregov https://www.kaggle.com/code/boojum/connecting-voxel-spaces","metadata":{}},{"cell_type":"code","source":"#Code by Michael Beregov https://www.kaggle.com/code/boojum/connecting-voxel-spaces\n\ncosine_from_dcm = pydicom.dcmread(filenamesDICOM[1]).ImageOrientationPatient\ncosine_from_dcm","metadata":{"execution":{"iopub.status.busy":"2022-08-02T19:16:47.943895Z","iopub.execute_input":"2022-08-02T19:16:47.944258Z","iopub.status.idle":"2022-08-02T19:16:47.958694Z","shell.execute_reply.started":"2022-08-02T19:16:47.944232Z","shell.execute_reply":"2022-08-02T19:16:47.957657Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"\"As you can see, it has information only about two rows (6 numbers instead of 9). We can calculate the third row though. It's perpendicular two the first two row vectors, so we can calculate their cross product:\"\n\nBy Michael Beregov https://www.kaggle.com/code/boojum/connecting-voxel-spaces","metadata":{}},{"cell_type":"code","source":"np.cross(cosine_from_dcm[:3], cosine_from_dcm[3:])","metadata":{"execution":{"iopub.status.busy":"2022-08-02T19:16:56.478007Z","iopub.execute_input":"2022-08-02T19:16:56.478346Z","iopub.status.idle":"2022-08-02T19:16:56.484159Z","shell.execute_reply.started":"2022-08-02T19:16:56.47832Z","shell.execute_reply":"2022-08-02T19:16:56.483507Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#Resampling\n\nBack to resampling. It's a bit more complicated but still pretty straightforward:","metadata":{}},{"cell_type":"code","source":"#Code by Michael Beregov https://www.kaggle.com/code/boojum/connecting-voxel-spaces\n\ndef resample(image, ref_image):\n\n    resampler = sitk.ResampleImageFilter()\n    resampler.SetReferenceImage(ref_image)\n    resampler.SetInterpolator(sitk.sitkLinear)\n    \n    resampler.SetTransform(sitk.AffineTransform(image.GetDimension()))\n\n    resampler.SetOutputSpacing(ref_image.GetSpacing())\n\n    resampler.SetSize(ref_image.GetSize())\n\n    resampler.SetOutputDirection(ref_image.GetDirection())\n\n    resampler.SetOutputOrigin(ref_image.GetOrigin())\n\n    resampler.SetDefaultPixelValue(image.GetPixelIDValue())\n\n    resampled_image = resampler.Execute(image)\n    \n    return resampled_image","metadata":{"execution":{"iopub.status.busy":"2022-08-02T19:17:04.217936Z","iopub.execute_input":"2022-08-02T19:17:04.218448Z","iopub.status.idle":"2022-08-02T19:17:04.22654Z","shell.execute_reply.started":"2022-08-02T19:17:04.218423Z","shell.execute_reply":"2022-08-02T19:17:04.22496Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#We don't have Flair","metadata":{}},{"cell_type":"code","source":"flair_resampled = resample(flair_sitk, t1_sitk)#No FLAIR!","metadata":{"execution":{"iopub.status.busy":"2022-08-02T19:17:10.052232Z","iopub.execute_input":"2022-08-02T19:17:10.052566Z","iopub.status.idle":"2022-08-02T19:17:11.372684Z","shell.execute_reply.started":"2022-08-02T19:17:10.052541Z","shell.execute_reply":"2022-08-02T19:17:11.371053Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"t1_sitk_array = sitk.GetArrayFromImage(t1_sitk)\nflair_resampled_array = sitk.GetArrayFromImage(flair_resampled)","metadata":{"execution":{"iopub.status.busy":"2022-08-02T19:17:30.757748Z","iopub.execute_input":"2022-08-02T19:17:30.758112Z","iopub.status.idle":"2022-08-02T19:17:30.936408Z","shell.execute_reply.started":"2022-08-02T19:17:30.758087Z","shell.execute_reply":"2022-08-02T19:17:30.934911Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#Code by Michael Beregov https://www.kaggle.com/code/boojum/connecting-voxel-spaces\n\nplt.figure(figsize=(12,6))\nplt.subplot(121)\nplt.imshow(t1_sitk_array[t1_sitk_array.shape[0]//2,:,:])\nplt.subplot(122)\nplt.imshow(flair_resampled_array[flair_resampled_array.shape[0]//2,:,:]);","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-08-02T19:17:36.322495Z","iopub.execute_input":"2022-08-02T19:17:36.322868Z","iopub.status.idle":"2022-08-02T19:17:36.84187Z","shell.execute_reply.started":"2022-08-02T19:17:36.322835Z","shell.execute_reply":"2022-08-02T19:17:36.836689Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def normalize(data):\n    return (data - np.min(data)) / (np.max(data) - np.min(data))","metadata":{"execution":{"iopub.status.busy":"2022-08-02T19:17:48.544142Z","iopub.execute_input":"2022-08-02T19:17:48.544464Z","iopub.status.idle":"2022-08-02T19:17:48.550223Z","shell.execute_reply.started":"2022-08-02T19:17:48.544439Z","shell.execute_reply":"2022-08-02T19:17:48.549198Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\nfilenamesDICOM = reader.GetGDCMSeriesFileNames(f'{train_path}/{train_dirs[2]}')#We don't have /T1w\nreader.SetFileNames(filenamesDICOM)\nt1_sitk = reader.Execute()\n\nfilenamesDICOM = reader.GetGDCMSeriesFileNames(f'{train_path}/{train_dirs[2]}')#We don't have /FLAIR\nreader.SetFileNames(filenamesDICOM)\nflair_sitk = reader.Execute()\n\nfilenamesDICOM = reader.GetGDCMSeriesFileNames(f'{train_path}/{train_dirs[2]}')#We don't have /T2w\nreader.SetFileNames(filenamesDICOM)\nt2_sitk = reader.Execute()\n\nflair_resampled = resample(flair_sitk, t1_sitk)\nt2_resampled = resample(t2_sitk, t1_sitk)\n\nt1_sitk_array = normalize(sitk.GetArrayFromImage(t1_sitk))\nflair_resampled_array = normalize(sitk.GetArrayFromImage(flair_resampled))\nt2_resampled_array = normalize(sitk.GetArrayFromImage(t2_resampled))\n\nstacked = np.stack([t1_sitk_array, t2_resampled_array, flair_resampled_array,])\n\nto_rgb = stacked[:,t1_sitk_array.shape[0]//2,:,:].transpose(1,2,0)\nim = Image.fromarray((to_rgb * 255).astype(np.uint8))","metadata":{"execution":{"iopub.status.busy":"2022-08-02T19:17:59.386737Z","iopub.execute_input":"2022-08-02T19:17:59.387068Z","iopub.status.idle":"2022-08-02T19:18:10.740984Z","shell.execute_reply.started":"2022-08-02T19:17:59.387043Z","shell.execute_reply":"2022-08-02T19:18:10.739405Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#I didn't even expect to have an image giving that we don't have FLAIR, T1w and T2w","metadata":{}},{"cell_type":"code","source":"im","metadata":{"execution":{"iopub.status.busy":"2022-08-02T19:18:22.242377Z","iopub.execute_input":"2022-08-02T19:18:22.24324Z","iopub.status.idle":"2022-08-02T19:18:22.316353Z","shell.execute_reply.started":"2022-08-02T19:18:22.243189Z","shell.execute_reply":"2022-08-02T19:18:22.315077Z"},"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#Let's resample some more volumes and look at some more pictures.","metadata":{}},{"cell_type":"code","source":"#Code by Michael Beregov https://www.kaggle.com/code/boojum/connecting-voxel-spaces\n\nreader = sitk.ImageSeriesReader()\nreader.LoadPrivateTagsOn()\nfilenamesDICOM = reader.GetGDCMSeriesFileNames(f'{train_path}/{train_dirs[20]}')#We don't have T1w\nreader.SetFileNames(filenamesDICOM)\nt1_reference = reader.Execute()","metadata":{"execution":{"iopub.status.busy":"2022-08-02T19:18:42.237431Z","iopub.execute_input":"2022-08-02T19:18:42.237768Z","iopub.status.idle":"2022-08-02T19:18:48.703857Z","shell.execute_reply.started":"2022-08-02T19:18:42.237742Z","shell.execute_reply":"2022-08-02T19:18:48.702402Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"sitk.GetArrayFromImage(t1_reference).shape","metadata":{"execution":{"iopub.status.busy":"2022-08-02T19:18:51.554649Z","iopub.execute_input":"2022-08-02T19:18:51.555014Z","iopub.status.idle":"2022-08-02T19:18:51.593226Z","shell.execute_reply.started":"2022-08-02T19:18:51.554989Z","shell.execute_reply":"2022-08-02T19:18:51.591789Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plt.imshow(sitk.GetArrayFromImage(t1_reference)[15,:,:]);","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-08-02T19:18:56.973501Z","iopub.execute_input":"2022-08-02T19:18:56.973861Z","iopub.status.idle":"2022-08-02T19:18:57.16424Z","shell.execute_reply.started":"2022-08-02T19:18:56.973835Z","shell.execute_reply":"2022-08-02T19:18:57.163268Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#Below, We don't have FLAIR, T1wCE, T2w, therefore we can't resample\n\nError: Non uniform sampling or missing slices detected,  maximum nonuniformity:0.001","metadata":{}},{"cell_type":"code","source":"#Code by Michael Beregov https://www.kaggle.com/code/boojum/connecting-voxel-spaces\n\nfig, axs = plt.subplots(5,4, figsize=(12, 18), facecolor='w', edgecolor='k', dpi=100)\naxs = axs.ravel()\n\nfor i, folder in enumerate(train_dirs[:20]):\n    filenamesDICOM = reader.GetGDCMSeriesFileNames(f'{train_path}/{folder}')#We don't have FLAIR\n    reader.SetFileNames(filenamesDICOM)\n    flair = reader.Execute()\n    \n    flair_resampled = resample(flair, t1_reference)\n    flair_resampled = normalize(sitk.GetArrayFromImage(flair_resampled))\n        \n    filenamesDICOM = reader.GetGDCMSeriesFileNames(f'{train_path}/{folder}')#We don't have T1wCE\n    reader.SetFileNames(filenamesDICOM)\n    t1 = reader.Execute()\n\n    filenamesDICOM = reader.GetGDCMSeriesFileNames(f'{train_path}/{folder}')#We don't have T2w\n    reader.SetFileNames(filenamesDICOM)\n    t2 = reader.Execute()\n        \n    t1_resampled = resample(t1, t1_reference)\n    t1_resampled = normalize(sitk.GetArrayFromImage(t1_resampled))\n\n    t2_resampled = resample(t2, t1_reference)\n    t2_resampled = normalize(sitk.GetArrayFromImage(t2_resampled))\n        \n    stacked = np.stack([t1_resampled, t2_resampled, flair_resampled])\n         \n    to_rgb = stacked[:,18,:,:].transpose(1,2,0)\n    im = Image.fromarray((to_rgb * 255).astype(np.uint8))\n    axs[i].imshow(im)","metadata":{"_kg_hide-output":true,"execution":{"iopub.status.busy":"2022-08-02T19:19:28.867491Z","iopub.execute_input":"2022-08-02T19:19:28.867839Z","iopub.status.idle":"2022-08-02T19:24:29.110874Z","shell.execute_reply.started":"2022-08-02T19:19:28.867813Z","shell.execute_reply":"2022-08-02T19:24:29.109738Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"\"You can also resample these images into coronal plane or saggital plane. Or resample them into axial plane, but using another patient as reference (who knows, maybe it's a good way to augment the data).\"\n\nBy Michael Beregov https://www.kaggle.com/code/boojum/connecting-voxel-spaces","metadata":{}},{"cell_type":"markdown","source":"#Acknowledgements:\n\nMichael Beregov https://www.kaggle.com/code/boojum/connecting-voxel-spaces\n\nNiBabel:\n\nMost work on NiBabel so far has been by Matthew Brett (MB), Chris Markiewicz\n(CM), Michael Hanke (MH), Marc-Alexandre Côté (MC), Ben Cipollini (BC), Paul\nMcCarthy (PM), Chris Cheng (CC), Yaroslav Halchenko (YOH), Satra Ghosh (SG),\nEric Larson (EL), Demian Wassermann, Stephan Gerhard and Ross Markello (RM).\n\nhttps://github.com/nipy/nibabel/releases\n\nhttps://nipy.org/nibabel/#","metadata":{}}]}