{"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":"# 01 Image registration helper code\nCode based on https://www.kaggle.com/boojum/connecting-voxel-spaces/notebook\n\ntwo approaches are outlined in the reference notebook: \n1. transforming scans to scanner space with affine matrix defined in dcm metadata --> resampling to match one channel\n    this approach is slower, and if there is any disruption to scanning or meta data input, it can lead to incorrect results\n2. simpleitk image registration\n    this approach is more general \n\nBacklog: \n1. **(DONE)** view a sample of img_registration --> scaling --> crop pipeline\n2. **(DONE)** image registration function: take reference channel as argument\n3. **(DONE)** image registration function: re-order voxels if needed \n4. check values distribution based on metadata (MR Acquisition and Pixel encoding)\n5. log scaling","metadata":{}},{"cell_type":"code","source":"import os\nimport sys \n\n#from tqdm import tqdm\nfrom pathlib import Path\n#from PIL import Image\n\nimport pydicom\nimport numpy as np\nnp.set_printoptions(precision=4)\nnp.set_printoptions(suppress=True)\n\nimport pandas as pd\n\nimport plotly.express as px\nfrom ipywidgets import widgets\nfrom IPython.display import display, clear_output, Image\n#import matplotlib.pyplot as plt\n#import matplotlib.colors\n\n#import nibabel as nib\nimport SimpleITK as sitk\nsitk.ProcessObject_SetGlobalWarningDisplay(False)\n\ntrain_path = '../input/rsna-miccai-brain-tumor-radiogenomic-classification/train/'\ntrain_dirs = os.listdir(train_path)\n\nDATASET = 'train'\nscan_types = ['FLAIR','T1w','T1wCE','T2w']\ndata_root = Path(\"../input/rsna-miccai-brain-tumor-radiogenomic-classification\")","metadata":{"execution":{"iopub.status.busy":"2021-10-26T12:26:50.484724Z","iopub.execute_input":"2021-10-26T12:26:50.485042Z","iopub.status.idle":"2021-10-26T12:26:52.673236Z","shell.execute_reply.started":"2021-10-26T12:26:50.484962Z","shell.execute_reply":"2021-10-26T12:26:52.672456Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def resample(image, ref_image):\n    resampler = sitk.ResampleImageFilter()\n    resampler.SetReferenceImage(ref_image) \n    resampler.SetInterpolator(sitk.sitkLinear)\n    resampler.SetTransform(sitk.AffineTransform(image.GetDimension()))\n    resampler.SetOutputSpacing(ref_image.GetSpacing())\n    resampler.SetSize(ref_image.GetSize())\n    resampler.SetOutputDirection(ref_image.GetDirection())\n    resampler.SetOutputOrigin(ref_image.GetOrigin())\n    resampler.SetDefaultPixelValue(image.GetPixelIDValue())\n    resamped_image = resampler.Execute(image)\n    return resamped_image\n\n\ndef image_registration(study_id, ref_channel='FLAIR'):\n    \n    reader = sitk.ImageSeriesReader()\n    reader.LoadPrivateTagsOn()\n    \n    filenamesDICOM = reader.GetGDCMSeriesFileNames(str(data_root.joinpath(DATASET,study_id,'FLAIR')))\n    reader.SetFileNames(filenamesDICOM)\n    flair_sitk = reader.Execute()\n    \n    filenamesDICOM = reader.GetGDCMSeriesFileNames(str(data_root.joinpath(DATASET,study_id,'T1w')))\n    reader.SetFileNames(filenamesDICOM)\n    t1_sitk = reader.Execute()\n    \n    filenamesDICOM = reader.GetGDCMSeriesFileNames(str(data_root.joinpath(DATASET,study_id,'T1wCE')))\n    reader.SetFileNames(filenamesDICOM)\n    t1wce_sitk = reader.Execute()\n\n    filenamesDICOM = reader.GetGDCMSeriesFileNames(str(data_root.joinpath(DATASET,study_id,'T2w')))\n    reader.SetFileNames(filenamesDICOM)\n    t2_sitk = reader.Execute()\n\n    ref_sitk_dct = {\n        'FLAIR': flair_sitk,\n        'T1w': t1_sitk,\n        'T1wCE': t1wce_sitk,\n        'T2w': t2_sitk,\n    }\n    ref_sitk = ref_sitk_dct[ref_channel]\n    \n    flair_resampled = resample(flair_sitk, ref_sitk)\n    t1_resampled = resample(t1_sitk, ref_sitk)\n    t1wce_resampled = resample(t1wce_sitk, ref_sitk)\n    t2_resampled = resample(t2_sitk, ref_sitk)\n    \n    flair_array = sitk.GetArrayFromImage(flair_resampled)\n    t1_array = sitk.GetArrayFromImage(t1_resampled) \n    t1wce_array = sitk.GetArrayFromImage(t1wce_resampled) \n    t2_array = sitk.GetArrayFromImage(t2_resampled)\n    \n    stacked = np.stack([flair_array, t1_array, t1wce_array, t2_array,])\n    \n    ref_dir = ref_sitk.GetDirection()\n    ref_dir = np.abs(np.round(np.array(ref_dir)))\n    \n    if np.array_equal(ref_dir, np.array([1, 0, 0, 0, 0, 1, 0, 1, 0,])):\n        return stacked.transpose((0, 2, 1, 3))\n    elif np.array_equal(ref_dir, np.array([1, 0, 0, 0, 1, 0, 0, 0, 1,])):\n        return stacked\n    elif np.array_equal(ref_dir, np.array([0, 0, 1, 1, 0, 0, 0, 1, 0,])):\n        return stacked.transpose((0, 2, 3, 1))\n    \n\ndef min_max_scaling(img):\n    mins = img.min(axis=(1,2,3), keepdims=True)\n    maxs = img.max(axis=(1,2,3), keepdims=True)\n    scaled = (img - mins) / (maxs-mins)\n    return scaled\n\n\n# def crop(img, min_threshold = 0.0001):\n#     if img.sum() == 0:\n#         return voxel\n    \n#     keep = (img.mean(axis=(0, 1, 2)) > min_threshold)\n#     img = img[:, :, :, keep]\n    \n#     keep = (img.mean(axis=(0, 1, 3)) > min_threshold)\n#     img = img[:, :, keep, :]\n    \n#     keep = (img.mean(axis=(0, 2, 3)) > min_threshold)\n#     img = img[:, keep, :, :]\n#     return img\n\n\ndef crop(img, scan_type):\n    if img.sum() == 0:\n        return voxel\n    \n    channel_map = {\n        'FLAIR': 0,\n        'T1w': 1,\n        'T1wCE': 2,\n        'T2w': 3\n    }\n    \n    c = channel_map[scan_type]\n    \n    keep = (img[c].mean(axis=(0, 1)) > 0)\n    img = img[:, :, :, keep]\n    \n    keep = (img[c].mean(axis=(0, 2)) > 0)\n    img = img[:, :, keep, :]\n    \n    keep = (img[c].mean(axis=(1, 2)) > 0)\n    img = img[:, keep, :, :]\n    return img\n\ndef log_scaling(img):\n    mins = np.log(img.min(axis=(1,2,3), keepdims=True) + 100)\n    maxs = np.log(np.percentile(img, 99, axis=(1,2,3), keepdims=True) + 1)\n    print(mins)\n    print(maxs)\n    scaled = (np.log(img + 1) - mins) / (maxs-mins)\n    clipped = np.clip(scaled, 0, 1)\n    return clipped","metadata":{"execution":{"iopub.status.busy":"2021-10-26T12:42:04.328562Z","iopub.execute_input":"2021-10-26T12:42:04.328909Z","iopub.status.idle":"2021-10-26T12:42:04.357721Z","shell.execute_reply.started":"2021-10-26T12:42:04.328858Z","shell.execute_reply":"2021-10-26T12:42:04.356470Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# from matplotlib import animation, rc\n# rc('animation', html='jshtml')\n\n\n# def display_img(img):\n#     rgb = img[:,img.shape[1]//2,:,:].transpose(1,2,0)\n#     im = Image.fromarray((rgb * 255).astype(np.uint8))\n#     return im\n\n\n# def get_slice(img, n):\n#     rgb = img[:,n].transpose((1,2,0))\n#     return (rgb * 255).astype(np.uint8)\n\n\n# def get_channel(img, n_slice, channel):\n#     rgb = img[channel,n_slice]\n#     return (rgb * 255).astype(np.uint8)\n\n\n# def anim_breakdown(ims, title=''):\n#     fig, (ax1, ax2, ax3, ax4) = plt.subplots(1,4 ,figsize=(24,6))\n    \n#     plt.axis('off')\n#     fig.patch.set_facecolor('white')\n\n#     im1 = ax1.imshow(get_channel(ims, n_slice=0, channel=0))\n#     im2 = ax2.imshow(get_channel(ims, n_slice=0, channel=1))\n#     im3 = ax3.imshow(get_channel(ims, n_slice=0, channel=2))\n#     im4 = ax4.imshow(get_channel(ims, n_slice=0, channel=3))\n    \n#     ax1.title.set_text('FLAIR')\n#     ax2.title.set_text('T1w')\n#     ax3.title.set_text('T1wCE')\n#     ax4.title.set_text('T2w')\n   \n\n#     def animate_func(i):\n#         im1.set_array(get_channel(ims, n_slice=i, channel=0))\n#         im2.set_array(get_channel(ims, n_slice=i, channel=1))\n#         im3.set_array(get_channel(ims, n_slice=i, channel=2))\n#         im4.set_array(get_channel(ims, n_slice=i, channel=3))\n#         return fig\n    \n#     plt.close()\n#     return animation.FuncAnimation(fig, animate_func, frames = ims.shape[1], interval = 1000//24)","metadata":{"execution":{"iopub.status.busy":"2021-10-26T08:08:27.358527Z","iopub.execute_input":"2021-10-26T08:08:27.358993Z","iopub.status.idle":"2021-10-26T08:08:27.378683Z","shell.execute_reply.started":"2021-10-26T08:08:27.35896Z","shell.execute_reply":"2021-10-26T08:08:27.377593Z"},"_kg_hide-input":true,"jupyter":{"source_hidden":true},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#ref_scantype = 'T1wCE'\n#reg_im = image_registration(train_dirs[2], ref_scantype)\n#print(f'{ref_scantype} aligned image {reg_im.shape}')\n#anim_breakdown(crop(min_max_scaling(reg_im)))","metadata":{"execution":{"iopub.status.busy":"2021-10-26T08:08:27.381074Z","iopub.execute_input":"2021-10-26T08:08:27.381875Z","iopub.status.idle":"2021-10-26T08:08:42.553433Z","shell.execute_reply.started":"2021-10-26T08:08:27.381837Z","shell.execute_reply":"2021-10-26T08:08:42.552305Z"},"jupyter":{"outputs_hidden":true,"source_hidden":true},"collapsed":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"cropping based on https://www.kaggle.com/ren4yu/normalized-voxels-align-planes-and-crop\nsince linear interpolation is used, threshold is required instead of checking for zeros","metadata":{}},{"cell_type":"markdown","source":"Animation code based on https://www.kaggle.com/ihelon/brain-tumor-eda-with-animations-and-modeling","metadata":{}},{"cell_type":"code","source":"def show_all_alignments(study_path):\n    \n    for scan_type in scan_types:\n        reg_im = image_registration(study_path, scan_type)\n        im = crop(min_max_scaling(reg_im), scan_type)\n\n        fig = px.imshow(im, \n                        animation_frame=1, \n                        facet_col=0, \n                        binary_string=True, \n                        title=f'Reference Channel = {scan_type}, Dimensions = {im.shape}')\n        \n        for i, label in enumerate(scan_types):\n            fig.layout.annotations[i]['text'] = label\n\n        fig.show()\n\nshow_all_alignments(train_dirs[2])","metadata":{"execution":{"iopub.status.busy":"2021-10-26T12:42:08.990600Z","iopub.execute_input":"2021-10-26T12:42:08.991691Z","iopub.status.idle":"2021-10-26T12:43:03.664169Z","shell.execute_reply.started":"2021-10-26T12:42:08.991634Z","shell.execute_reply":"2021-10-26T12:43:03.662952Z"},"trusted":true},"execution_count":null,"outputs":[]}]}