{"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":"# Normalized Voxels: Align Planes, Adjust Contrast, and Crop\n\nAs shown in several notebooks, MRI plane type (Axial, Coronal, and Sagittal) is not consistent among patients or MRI scan types (FLAIR, T1w, T1wCE, T2w).\nWhile augmentations might alleviate this inconsistency, it is better to train models using MRI voxels that are consistent in terms of plane type.\nThis notebook shows we can obtain normalized voxels by appropriately rotating MRI voxels.\nI found that simply rotating MRI voxels is not enough because the order of planes is also inconsistent in some cases.\nFor example, even with the same Sagittal type, some of scans were in left-to-right order, while others were the other way around.\nAs Instance Number does not help, Image Position (Patient) is used in this notebook to reorder stacked images.\nAfter normalized voxels with respect to planes, contrast is adjusted and then voxels are cropped.\nFinally, the voxel is resized to arbitrary fixed size.\n\n## Normalized Voxel Datasets\n\nThe normalized voxels created the above procedure were stored as a dataset. Please refer to the second half of this notebook.","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19"}},{"cell_type":"code","source":"from pathlib import Path\nimport numpy as np\nimport cv2\nimport pydicom\nimport matplotlib.pyplot as plt\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-08-28T15:44:23.681524Z","iopub.execute_input":"2021-08-28T15:44:23.681850Z","iopub.status.idle":"2021-08-28T15:44:23.687075Z","shell.execute_reply.started":"2021-08-28T15:44:23.681811Z","shell.execute_reply":"2021-08-28T15:44:23.685976Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# https://www.kaggle.com/arnabs007/part-1-rsna-miccai-btrc-understanding-the-data\n# https://www.kaggle.com/davidbroberts/determining-mr-image-planes\ndef get_image_plane(data):\n    x1, y1, _, x2, y2, _ = [round(j) for j in data.ImageOrientationPatient]\n    cords = [x1, y1, x2, y2]\n\n    if cords == [1, 0, 0, 0]:\n        return 'Coronal'\n    elif cords == [1, 0, 0, 1]:\n        return 'Axial'\n    elif cords == [0, 1, 0, 0]:\n        return 'Sagittal'\n    else:\n        return 'Unknown'","metadata":{"execution":{"iopub.status.busy":"2021-08-28T15:44:24.051729Z","iopub.execute_input":"2021-08-28T15:44:24.052104Z","iopub.status.idle":"2021-08-28T15:44:24.058076Z","shell.execute_reply.started":"2021-08-28T15:44:24.052074Z","shell.execute_reply":"2021-08-28T15:44:24.057062Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# diff from dataset\ndef crop_voxel(voxel):\n    if voxel.sum() == 0:\n        return voxel\n    keep = (voxel.mean(axis=(0, 1)) > 1)\n    voxel = voxel[:, :, keep]\n    keep = (voxel.mean(axis=(0, 2)) > 1)\n    voxel = voxel[:, keep]\n    keep = (voxel.mean(axis=(1, 2)) > 1)\n    voxel = voxel[keep]\n    return voxel","metadata":{"execution":{"iopub.status.busy":"2021-08-28T15:44:24.544512Z","iopub.execute_input":"2021-08-28T15:44:24.545008Z","iopub.status.idle":"2021-08-28T15:44:24.551067Z","shell.execute_reply.started":"2021-08-28T15:44:24.544978Z","shell.execute_reply":"2021-08-28T15:44:24.550116Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def get_voxel(study_id, scan_type):\n    imgs = []\n    dcm_dir = data_root.joinpath(DATASET, study_id, scan_type)\n    dcm_paths = sorted(dcm_dir.glob(\"*.dcm\"), key=lambda x: int(x.stem.split(\"-\")[-1]))\n    positions = []\n    \n    for dcm_path in dcm_paths:\n        img = pydicom.dcmread(str(dcm_path))\n        imgs.append(img.pixel_array)\n        positions.append(img.ImagePositionPatient)\n        \n    plane = get_image_plane(img)\n    voxel = np.stack(imgs)\n    \n    # diff from dataset\n    voxel = crop_voxel(voxel)\n    \n    # reorder planes if needed and rotate voxel\n    if plane == \"Coronal\":\n        if positions[0][1] < positions[-1][1]:\n            voxel = voxel[::-1]\n            print(f\"{study_id} {scan_type} {plane} reordered\")\n        voxel = voxel.transpose((1, 0, 2))\n    elif plane == \"Sagittal\":\n        if positions[0][0] < positions[-1][0]:\n            voxel = voxel[::-1]\n            print(f\"{study_id} {scan_type} {plane} reordered\")\n        voxel = voxel.transpose((1, 2, 0))\n        voxel = np.rot90(voxel, 2, axes=(1, 2))\n    elif plane == \"Axial\":\n        if positions[0][2] > positions[-1][2]:\n            voxel = voxel[::-1]\n            print(f\"{study_id} {scan_type} {plane} reordered\")\n        voxel = np.rot90(voxel, 2)\n    else:\n        raise ValueError(f\"Unknown plane {plane}\")\n    return voxel, plane","metadata":{"execution":{"iopub.status.busy":"2021-08-28T15:44:24.906197Z","iopub.execute_input":"2021-08-28T15:44:24.906704Z","iopub.status.idle":"2021-08-28T15:44:24.917563Z","shell.execute_reply.started":"2021-08-28T15:44:24.906671Z","shell.execute_reply":"2021-08-28T15:44:24.916620Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def normalize_contrast(voxel):\n    if voxel.sum() == 0:\n        return voxel\n    voxel = voxel - np.min(voxel)\n    voxel = voxel / np.max(voxel)\n    voxel = (voxel * 255).astype(np.uint8)\n    return voxel","metadata":{"execution":{"iopub.status.busy":"2021-08-28T15:44:25.265824Z","iopub.execute_input":"2021-08-28T15:44:25.266388Z","iopub.status.idle":"2021-08-28T15:44:25.271941Z","shell.execute_reply.started":"2021-08-28T15:44:25.266356Z","shell.execute_reply":"2021-08-28T15:44:25.270854Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Sample planes along the longest axis and resize the sampled planes.\nBy sampling along the longest axis, the degradation due to sampling is minimized.\nThe best way is to resize twice (e.g. (x, y) axis then (y, z) axis) but it is computationally expensive.","metadata":{}},{"cell_type":"code","source":"def resize_voxel(voxel, sz=64):\n    output = np.zeros((sz, sz, sz), dtype=np.uint8)\n\n    if np.argmax(voxel.shape) == 0:\n        for i, s in enumerate(np.linspace(0, voxel.shape[0] - 1, sz)):\n            output[i] = cv2.resize(voxel[int(s)], (sz, sz))\n    elif np.argmax(voxel.shape) == 1:\n        for i, s in enumerate(np.linspace(0, voxel.shape[1] - 1, sz)):\n            output[:, i] = cv2.resize(voxel[:, int(s)], (sz, sz))\n    elif np.argmax(voxel.shape) == 2:\n        for i, s in enumerate(np.linspace(0, voxel.shape[2] - 1, sz)):\n            output[:, :, i] = cv2.resize(voxel[:, :, int(s)], (sz, sz))\n\n    return output","metadata":{"execution":{"iopub.status.busy":"2021-08-28T15:44:26.030817Z","iopub.execute_input":"2021-08-28T15:44:26.031307Z","iopub.status.idle":"2021-08-28T15:44:26.040934Z","shell.execute_reply.started":"2021-08-28T15:44:26.031277Z","shell.execute_reply":"2021-08-28T15:44:26.039687Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"case_id = '00386'\nscan_t = 'T1wCE'\nsz = 128\nvoxel, _ = get_voxel(case_id, scan_t)\nvoxel = normalize_contrast(voxel)\n# voxel = crop_voxel(voxel)\nvoxel1 = resize_voxel(voxel, sz)\n\nfname = f'../input/rsna-miccai-voxel-{sz}-dataset/voxel/{DATASET}/{case_id}/{scan_t}.npy'\nvoxel2 = np.load(fname)\nprint(f'Normalized voxel\\'s shape: {voxel1.shape}')\nprint(f'Saved voxel\\'s shape: {voxel2.shape}')\nassert np.array_equiv(voxel1, voxel2), 'Normalized voxel differ from array saved on dataset.'","metadata":{"execution":{"iopub.status.busy":"2021-08-28T15:44:26.661081Z","iopub.execute_input":"2021-08-28T15:44:26.661445Z","iopub.status.idle":"2021-08-28T15:44:27.789173Z","shell.execute_reply.started":"2021-08-28T15:44:26.661414Z","shell.execute_reply":"2021-08-28T15:44:27.788224Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}