{"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":"<div style=\"text-align:center\"> <img src=\"https://sun9-45.userapi.com/impg/Flbnug2OUli1ecXsoIKeUasIGXGj_5hqjX4cRg/z2nfz8-b3a0.jpg?size=2560x1153&quality=96&sign=0536543610f1655d967af88dbc775e98&type=album\" width=\"800\">","metadata":{}},{"cell_type":"code","source":"import 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-miccai-brain-tumor-radiogenomic-classification/train/'","metadata":{"execution":{"iopub.status.busy":"2021-07-24T13:33:49.457637Z","iopub.execute_input":"2021-07-24T13:33:49.457993Z","iopub.status.idle":"2021-07-24T13:33:50.301565Z","shell.execute_reply.started":"2021-07-24T13:33:49.457964Z","shell.execute_reply":"2021-07-24T13:33:50.300444Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_dirs = os.listdir(train_path)","metadata":{"execution":{"iopub.status.busy":"2021-07-24T13:33:50.303057Z","iopub.execute_input":"2021-07-24T13:33:50.303387Z","iopub.status.idle":"2021-07-24T13:33:50.350607Z","shell.execute_reply.started":"2021-07-24T13:33:50.303357Z","shell.execute_reply":"2021-07-24T13:33:50.349509Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plt.figure(figsize=(14,7))\nplt.subplot(121)\nplt.imshow(pydicom.dcmread(f'{train_path + train_dirs[0]}/T2w/Image-111.dcm').pixel_array)\nplt.subplot(122)\nplt.imshow(pydicom.dcmread(f'{train_path + train_dirs[0]}/T1w/Image-111.dcm').pixel_array)","metadata":{"execution":{"iopub.status.busy":"2021-07-24T13:33:50.352135Z","iopub.execute_input":"2021-07-24T13:33:50.352484Z","iopub.status.idle":"2021-07-24T13:33:50.878468Z","shell.execute_reply.started":"2021-07-24T13:33:50.352451Z","shell.execute_reply":"2021-07-24T13:33:50.87742Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Different spatial orientation of images is inconvinient. But we can make it the same. ","metadata":{}},{"cell_type":"markdown","source":"# 1. (not really important part) Affine matrix and simple resampling","metadata":{}},{"cell_type":"markdown","source":"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). ","metadata":{}},{"cell_type":"markdown","source":"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. ","metadata":{}},{"cell_type":"code","source":"reader = sitk.ImageSeriesReader()\nreader.LoadPrivateTagsOn()","metadata":{"execution":{"iopub.status.busy":"2021-07-24T11:29:20.434484Z","iopub.execute_input":"2021-07-24T11:29:20.434792Z","iopub.status.idle":"2021-07-24T11:29:20.446044Z","shell.execute_reply.started":"2021-07-24T11:29:20.434747Z","shell.execute_reply":"2021-07-24T11:29:20.444745Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"filenamesDICOM = reader.GetGDCMSeriesFileNames(f'{train_path}/{train_dirs[2]}/T1w')\nreader.SetFileNames(filenamesDICOM)\nt1_sitk = reader.Execute()","metadata":{"execution":{"iopub.status.busy":"2021-07-24T11:29:20.447228Z","iopub.execute_input":"2021-07-24T11:29:20.447516Z","iopub.status.idle":"2021-07-24T11:29:20.712751Z","shell.execute_reply.started":"2021-07-24T11:29:20.447489Z","shell.execute_reply":"2021-07-24T11:29:20.71148Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"sitk.WriteImage(t1_sitk,'t1.nii')","metadata":{"execution":{"iopub.status.busy":"2021-07-24T11:29:20.714276Z","iopub.execute_input":"2021-07-24T11:29:20.714578Z","iopub.status.idle":"2021-07-24T11:29:20.757181Z","shell.execute_reply.started":"2021-07-24T11:29:20.71455Z","shell.execute_reply":"2021-07-24T11:29:20.756231Z"},"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":"t1_nib = nib.load('t1.nii')\nt1_nib","metadata":{"execution":{"iopub.status.busy":"2021-07-24T11:29:20.758368Z","iopub.execute_input":"2021-07-24T11:29:20.758883Z","iopub.status.idle":"2021-07-24T11:29:20.770522Z","shell.execute_reply.started":"2021-07-24T11:29:20.758848Z","shell.execute_reply":"2021-07-24T11:29:20.768975Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"t1_nib_array = t1_nib.get_fdata() #the voxel array\nt1_nib_array[:3]","metadata":{"execution":{"iopub.status.busy":"2021-07-24T11:29:20.77249Z","iopub.execute_input":"2021-07-24T11:29:20.773293Z","iopub.status.idle":"2021-07-24T11:29:20.811317Z","shell.execute_reply.started":"2021-07-24T11:29:20.773237Z","shell.execute_reply":"2021-07-24T11:29:20.810036Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"t1_nib_array.shape","metadata":{"execution":{"iopub.status.busy":"2021-07-24T11:29:20.813718Z","iopub.execute_input":"2021-07-24T11:29:20.814091Z","iopub.status.idle":"2021-07-24T11:29:20.821816Z","shell.execute_reply.started":"2021-07-24T11:29:20.814059Z","shell.execute_reply":"2021-07-24T11:29:20.82037Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plt.imshow(t1_nib_array[:,:,t1_nib_array.shape[2]//2])","metadata":{"execution":{"iopub.status.busy":"2021-07-24T11:29:31.059655Z","iopub.execute_input":"2021-07-24T11:29:31.060113Z","iopub.status.idle":"2021-07-24T11:29:31.297943Z","shell.execute_reply.started":"2021-07-24T11:29:31.060063Z","shell.execute_reply":"2021-07-24T11:29:31.296712Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"So, this is the voxel array. Its orientation in scanner space is encoded in the affine matrix:","metadata":{}},{"cell_type":"code","source":"np.set_printoptions(precision=4, suppress=True)\nt1_nib.affine","metadata":{"execution":{"iopub.status.busy":"2021-07-24T11:29:41.760267Z","iopub.execute_input":"2021-07-24T11:29:41.760682Z","iopub.status.idle":"2021-07-24T11:29:41.769153Z","shell.execute_reply.started":"2021-07-24T11:29:41.760645Z","shell.execute_reply":"2021-07-24T11:29:41.767705Z"},"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.","metadata":{}},{"cell_type":"markdown","source":"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).","metadata":{"execution":{"iopub.status.busy":"2021-07-21T19:48:25.58817Z","iopub.execute_input":"2021-07-21T19:48:25.588484Z","iopub.status.idle":"2021-07-21T19:48:25.59415Z","shell.execute_reply.started":"2021-07-21T19:48:25.588456Z","shell.execute_reply":"2021-07-21T19:48:25.593035Z"}}},{"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        $$","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":"2021-07-24T11:31:12.583395Z","iopub.execute_input":"2021-07-24T11:31:12.583796Z","iopub.status.idle":"2021-07-24T11:31:12.59248Z","shell.execute_reply.started":"2021-07-24T11:31:12.583749Z","shell.execute_reply":"2021-07-24T11:31:12.591658Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The 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":"2021-07-24T11:31:13.692475Z","iopub.execute_input":"2021-07-24T11:31:13.692843Z","iopub.status.idle":"2021-07-24T11:31:13.702037Z","shell.execute_reply.started":"2021-07-24T11:31:13.692812Z","shell.execute_reply":"2021-07-24T11:31:13.700612Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Apparently, 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":{"execution":{"iopub.status.busy":"2021-07-21T20:19:28.660919Z","iopub.execute_input":"2021-07-21T20:19:28.661207Z","iopub.status.idle":"2021-07-21T20:19:28.667071Z","shell.execute_reply.started":"2021-07-21T20:19:28.661183Z","shell.execute_reply":"2021-07-21T20:19:28.665763Z"}}},{"cell_type":"code","source":"filenamesDICOM = reader.GetGDCMSeriesFileNames(f'{train_path}/{train_dirs[2]}/FLAIR')\nreader.SetFileNames(filenamesDICOM)\nflair_sitk = reader.Execute()\nsitk.WriteImage(flair_sitk,'flair.nii')\n\nflair_nib = nib.load('flair.nii')\nflair_nib_array = flair_nib.get_fdata()","metadata":{"execution":{"iopub.status.busy":"2021-07-24T11:42:25.240793Z","iopub.execute_input":"2021-07-24T11:42:25.241716Z","iopub.status.idle":"2021-07-24T11:42:28.58729Z","shell.execute_reply.started":"2021-07-24T11:42:25.241656Z","shell.execute_reply":"2021-07-24T11:42:28.585673Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plt.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":{"execution":{"iopub.status.busy":"2021-07-24T11:42:28.589656Z","iopub.execute_input":"2021-07-24T11:42:28.590153Z","iopub.status.idle":"2021-07-24T11:42:29.675126Z","shell.execute_reply.started":"2021-07-24T11:42:28.590097Z","shell.execute_reply":"2021-07-24T11:42:29.673572Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from nilearn.image import resample_img","metadata":{"execution":{"iopub.status.busy":"2021-07-24T11:42:29.678051Z","iopub.execute_input":"2021-07-24T11:42:29.678525Z","iopub.status.idle":"2021-07-24T11:42:31.069898Z","shell.execute_reply.started":"2021-07-24T11:42:29.678474Z","shell.execute_reply":"2021-07-24T11:42:31.06854Z"},"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":"2021-07-24T11:42:31.071973Z","iopub.execute_input":"2021-07-24T11:42:31.072435Z","iopub.status.idle":"2021-07-24T11:42:41.856092Z","shell.execute_reply.started":"2021-07-24T11:42:31.072384Z","shell.execute_reply":"2021-07-24T11:42:41.854774Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plt.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":{"execution":{"iopub.status.busy":"2021-07-24T11:42:41.857515Z","iopub.execute_input":"2021-07-24T11:42:41.857827Z","iopub.status.idle":"2021-07-24T11:42:42.334266Z","shell.execute_reply.started":"2021-07-24T11:42:41.857796Z","shell.execute_reply":"2021-07-24T11:42:42.332777Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"It works rather slowly, but there is a faster way. ","metadata":{"execution":{"iopub.status.busy":"2021-07-21T20:34:39.135758Z","iopub.execute_input":"2021-07-21T20:34:39.136209Z","iopub.status.idle":"2021-07-21T20:34:39.14222Z","shell.execute_reply.started":"2021-07-21T20:34:39.13617Z","shell.execute_reply":"2021-07-21T20:34:39.140951Z"}}},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# 2. SimpleITK resampling","metadata":{"execution":{"iopub.status.busy":"2021-07-21T20:36:15.766371Z","iopub.execute_input":"2021-07-21T20:36:15.766827Z","iopub.status.idle":"2021-07-21T20:36:15.770927Z","shell.execute_reply.started":"2021-07-21T20:36:15.766754Z","shell.execute_reply":"2021-07-21T20:36:15.769916Z"}}},{"cell_type":"code","source":"filenamesDICOM = reader.GetGDCMSeriesFileNames(f'{train_path}/{train_dirs[2]}/T1w')\nreader.SetFileNames(filenamesDICOM)\nt1_sitk = reader.Execute()\nt1_sitk","metadata":{"execution":{"iopub.status.busy":"2021-07-24T11:42:42.335631Z","iopub.execute_input":"2021-07-24T11:42:42.33601Z","iopub.status.idle":"2021-07-24T11:42:42.563158Z","shell.execute_reply.started":"2021-07-24T11:42:42.335973Z","shell.execute_reply":"2021-07-24T11:42:42.562042Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"So, 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:","metadata":{"execution":{"iopub.status.busy":"2021-07-21T20:37:58.207601Z","iopub.execute_input":"2021-07-21T20:37:58.208267Z","iopub.status.idle":"2021-07-21T20:37:58.214299Z","shell.execute_reply.started":"2021-07-21T20:37:58.20823Z","shell.execute_reply":"2021-07-21T20:37:58.212918Z"}}},{"cell_type":"code","source":"t1_sitk.GetOrigin() # which is a translation column from the affine matrix, but negative","metadata":{"execution":{"iopub.status.busy":"2021-07-24T11:45:19.600927Z","iopub.execute_input":"2021-07-24T11:45:19.601529Z","iopub.status.idle":"2021-07-24T11:45:19.608351Z","shell.execute_reply.started":"2021-07-24T11:45:19.601489Z","shell.execute_reply":"2021-07-24T11:45:19.607299Z"},"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":"2021-07-24T11:45:19.890434Z","iopub.execute_input":"2021-07-24T11:45:19.890912Z","iopub.status.idle":"2021-07-24T11:45:19.899266Z","shell.execute_reply.started":"2021-07-24T11:45:19.89087Z","shell.execute_reply":"2021-07-24T11:45:19.897894Z"},"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":"2021-07-24T11:45:20.659936Z","iopub.execute_input":"2021-07-24T11:45:20.660296Z","iopub.status.idle":"2021-07-24T11:45:20.668341Z","shell.execute_reply.started":"2021-07-24T11:45:20.660265Z","shell.execute_reply":"2021-07-24T11:45:20.666826Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We 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).","metadata":{}},{"cell_type":"code","source":"x_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":"2021-07-24T12:01:14.999745Z","iopub.execute_input":"2021-07-24T12:01:15.00013Z","iopub.status.idle":"2021-07-24T12:01:15.009361Z","shell.execute_reply.started":"2021-07-24T12:01:15.000095Z","shell.execute_reply":"2021-07-24T12:01:15.007921Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Also, if we divide each column by the corresponding spacing number, we'll get a cosine matrix.","metadata":{}},{"cell_type":"code","source":"t1_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":"2021-07-24T12:44:21.638965Z","iopub.execute_input":"2021-07-24T12:44:21.63948Z","iopub.status.idle":"2021-07-24T12:44:21.72442Z","shell.execute_reply.started":"2021-07-24T12:44:21.639362Z","shell.execute_reply":"2021-07-24T12:44:21.722016Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"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. ","metadata":{"execution":{"iopub.status.busy":"2021-07-24T11:46:51.653135Z","iopub.execute_input":"2021-07-24T11:46:51.653577Z","iopub.status.idle":"2021-07-24T11:46:51.661329Z","shell.execute_reply.started":"2021-07-24T11:46:51.653534Z","shell.execute_reply":"2021-07-24T11:46:51.660238Z"}}},{"cell_type":"markdown","source":"Btw, this information can be acquired from DICOM files. There's a tag for this:","metadata":{}},{"cell_type":"code","source":"cosine_from_dcm = pydicom.dcmread(filenamesDICOM[1]).ImageOrientationPatient\ncosine_from_dcm","metadata":{"execution":{"iopub.status.busy":"2021-07-24T12:09:29.341923Z","iopub.execute_input":"2021-07-24T12:09:29.34238Z","iopub.status.idle":"2021-07-24T12:09:29.363907Z","shell.execute_reply.started":"2021-07-24T12:09:29.342339Z","shell.execute_reply":"2021-07-24T12:09:29.362652Z"},"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:","metadata":{}},{"cell_type":"code","source":"np.cross(cosine_from_dcm[:3], cosine_from_dcm[3:])","metadata":{"execution":{"iopub.status.busy":"2021-07-24T12:11:41.150964Z","iopub.execute_input":"2021-07-24T12:11:41.151444Z","iopub.status.idle":"2021-07-24T12:11:41.159466Z","shell.execute_reply.started":"2021-07-24T12:11:41.151403Z","shell.execute_reply":"2021-07-24T12:11:41.158476Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Back to resampling. It's a bit more complicated but still pretty straightforward:","metadata":{"execution":{"iopub.status.busy":"2021-07-24T12:12:42.06904Z","iopub.execute_input":"2021-07-24T12:12:42.069427Z","iopub.status.idle":"2021-07-24T12:12:42.074792Z","shell.execute_reply.started":"2021-07-24T12:12:42.069391Z","shell.execute_reply":"2021-07-24T12:12:42.073733Z"}}},{"cell_type":"code","source":"def 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    resamped_image = resampler.Execute(image)\n    \n    return resamped_image","metadata":{"execution":{"iopub.status.busy":"2021-07-24T12:14:18.465612Z","iopub.execute_input":"2021-07-24T12:14:18.46627Z","iopub.status.idle":"2021-07-24T12:14:18.473349Z","shell.execute_reply.started":"2021-07-24T12:14:18.466229Z","shell.execute_reply":"2021-07-24T12:14:18.471959Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"flair_resampled = resample(flair_sitk, t1_sitk)","metadata":{"execution":{"iopub.status.busy":"2021-07-24T12:14:18.696583Z","iopub.execute_input":"2021-07-24T12:14:18.697014Z","iopub.status.idle":"2021-07-24T12:14:18.821114Z","shell.execute_reply.started":"2021-07-24T12:14:18.696979Z","shell.execute_reply":"2021-07-24T12:14:18.819726Z"},"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":"2021-07-24T12:14:18.886928Z","iopub.execute_input":"2021-07-24T12:14:18.887275Z","iopub.status.idle":"2021-07-24T12:14:18.898938Z","shell.execute_reply.started":"2021-07-24T12:14:18.887244Z","shell.execute_reply":"2021-07-24T12:14:18.897882Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plt.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":{"execution":{"iopub.status.busy":"2021-07-24T12:14:19.065148Z","iopub.execute_input":"2021-07-24T12:14:19.065514Z","iopub.status.idle":"2021-07-24T12:14:19.486017Z","shell.execute_reply.started":"2021-07-24T12:14:19.06548Z","shell.execute_reply":"2021-07-24T12:14:19.484911Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"It works **much** faster. ","metadata":{"execution":{"iopub.status.busy":"2021-07-21T20:53:06.836892Z","iopub.execute_input":"2021-07-21T20:53:06.837292Z","iopub.status.idle":"2021-07-21T20:53:06.842321Z","shell.execute_reply.started":"2021-07-21T20:53:06.837259Z","shell.execute_reply":"2021-07-21T20:53:06.841373Z"}}},{"cell_type":"code","source":"def normalize(data):\n    return (data - np.min(data)) / (np.max(data) - np.min(data))","metadata":{"execution":{"iopub.status.busy":"2021-07-24T12:14:19.582171Z","iopub.execute_input":"2021-07-24T12:14:19.582812Z","iopub.status.idle":"2021-07-24T12:14:19.589289Z","shell.execute_reply.started":"2021-07-24T12:14:19.582733Z","shell.execute_reply":"2021-07-24T12:14:19.588321Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\nfilenamesDICOM = reader.GetGDCMSeriesFileNames(f'{train_path}/{train_dirs[2]}/T1w')\nreader.SetFileNames(filenamesDICOM)\nt1_sitk = reader.Execute()\n\nfilenamesDICOM = reader.GetGDCMSeriesFileNames(f'{train_path}/{train_dirs[2]}/FLAIR')\nreader.SetFileNames(filenamesDICOM)\nflair_sitk = reader.Execute()\n\nfilenamesDICOM = reader.GetGDCMSeriesFileNames(f'{train_path}/{train_dirs[2]}/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":"2021-07-24T12:35:34.119686Z","iopub.execute_input":"2021-07-24T12:35:34.120257Z","iopub.status.idle":"2021-07-24T12:35:44.245547Z","shell.execute_reply.started":"2021-07-24T12:35:34.120197Z","shell.execute_reply":"2021-07-24T12:35:44.243939Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"im","metadata":{"execution":{"iopub.status.busy":"2021-07-24T12:35:59.112932Z","iopub.execute_input":"2021-07-24T12:35:59.113441Z","iopub.status.idle":"2021-07-24T12:35:59.170699Z","shell.execute_reply.started":"2021-07-24T12:35:59.113401Z","shell.execute_reply":"2021-07-24T12:35:59.169536Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"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":"reader = sitk.ImageSeriesReader()\nreader.LoadPrivateTagsOn()\nfilenamesDICOM = reader.GetGDCMSeriesFileNames(f'{train_path}/{train_dirs[20]}/T1w')\nreader.SetFileNames(filenamesDICOM)\nt1_reference = reader.Execute()","metadata":{"execution":{"iopub.status.busy":"2021-07-24T12:18:25.641739Z","iopub.execute_input":"2021-07-24T12:18:25.642178Z","iopub.status.idle":"2021-07-24T12:18:26.126667Z","shell.execute_reply.started":"2021-07-24T12:18:25.642136Z","shell.execute_reply":"2021-07-24T12:18:26.125892Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"sitk.GetArrayFromImage(t1_reference).shape","metadata":{"execution":{"iopub.status.busy":"2021-07-24T12:19:07.196699Z","iopub.execute_input":"2021-07-24T12:19:07.1971Z","iopub.status.idle":"2021-07-24T12:19:07.207872Z","shell.execute_reply.started":"2021-07-24T12:19:07.197061Z","shell.execute_reply":"2021-07-24T12:19:07.206605Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plt.imshow(sitk.GetArrayFromImage(t1_reference)[15,:,:])","metadata":{"execution":{"iopub.status.busy":"2021-07-24T12:19:07.408167Z","iopub.execute_input":"2021-07-24T12:19:07.408614Z","iopub.status.idle":"2021-07-24T12:19:07.636092Z","shell.execute_reply.started":"2021-07-24T12:19:07.408572Z","shell.execute_reply":"2021-07-24T12:19:07.635027Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig, 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}/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}/T1wCE')\n    reader.SetFileNames(filenamesDICOM)\n    t1 = reader.Execute()\n\n    filenamesDICOM = reader.GetGDCMSeriesFileNames(f'{train_path}/{folder}/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":{"execution":{"iopub.status.busy":"2021-07-24T12:33:19.303729Z","iopub.execute_input":"2021-07-24T12:33:19.30438Z","iopub.status.idle":"2021-07-24T12:34:41.864077Z","shell.execute_reply.started":"2021-07-24T12:33:19.304336Z","shell.execute_reply":"2021-07-24T12:34:41.862537Z"}},"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). ","metadata":{}}]}