{"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 class='alert alert-info' style='text-align:center'><h1>Export Images By Plane</h1>\n- yet another MR processing notebook -</div>\n\n#### This notebook exports a series from each study in the training set into the same plane.\n- Since all series are not in the same plane, it makes sense to export them into the same plane.\n- Slices that are not in the specified plane already are reformated with SimpleITK.\n- Images are exported as JPG, ignoring empty slices (all black).\n- A zip file is created that contains a single series from each test study, in the plane you choose.","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2021-08-12T17:23:13.595475Z","iopub.execute_input":"2021-08-12T17:23:13.595923Z","iopub.status.idle":"2021-08-12T17:23:14.662964Z","shell.execute_reply.started":"2021-08-12T17:23:13.595826Z","shell.execute_reply":"2021-08-12T17:23:14.66168Z"}}},{"cell_type":"code","source":"import os\nimport sys \nimport numpy as np\nimport nibabel as nib\nimport SimpleITK as sitk\nimport matplotlib.pyplot as plt\nimport pydicom\nimport cv2\nimport shutil\nfrom IPython.display import FileLink\nfrom PIL import Image","metadata":{"execution":{"iopub.status.busy":"2021-09-24T13:06:53.303518Z","iopub.execute_input":"2021-09-24T13:06:53.304007Z","iopub.status.idle":"2021-09-24T13:06:54.35791Z","shell.execute_reply.started":"2021-09-24T13:06:53.303906Z","shell.execute_reply":"2021-09-24T13:06:54.356509Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Specify a series and plane to export\nPLANE = 'axial'\nSERIES = 'T1w'","metadata":{"execution":{"iopub.status.busy":"2021-09-24T13:06:54.359564Z","iopub.execute_input":"2021-09-24T13:06:54.359943Z","iopub.status.idle":"2021-09-24T13:06:54.364989Z","shell.execute_reply.started":"2021-09-24T13:06:54.359909Z","shell.execute_reply":"2021-09-24T13:06:54.363447Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Get a list of study directories\ntrain_path = '../input/rsna-miccai-brain-tumor-radiogenomic-classification/train/'\ntrain_dirs = os.listdir(train_path)\n\nif not os.path.exists(f'/kaggle/working/train'):\n    os.mkdir('/kaggle/working/train')","metadata":{"execution":{"iopub.status.busy":"2021-09-24T13:06:54.3673Z","iopub.execute_input":"2021-09-24T13:06:54.367704Z","iopub.status.idle":"2021-09-24T13:06:54.420766Z","shell.execute_reply.started":"2021-09-24T13:06:54.36766Z","shell.execute_reply":"2021-09-24T13:06:54.41961Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# This function gives a basic plane from the ImageOrientationPatient tag. It doesn't account for obliqueness. But we don't need to care about it.\n# Will return 'unknown' if the image isn't exactly orthogonal.\ndef get_image_plane(loc):\n    row_x = round(loc[0])\n    row_y = round(loc[1])\n    row_z = round(loc[2])\n    col_x = round(loc[3])\n    col_y = round(loc[4])\n    col_z = round(loc[5])\n    if (row_x, row_y, col_x, col_y) == (1,0,0,0):\n        return \"Coronal\"\n    if (row_x, row_y, col_x, col_y) == (0,1,0,0):\n        return \"Sagittal\"\n    if (row_x, row_y, col_x, col_y) == (1,0,0,1):\n        return \"Axial\"\n    return \"Unknown\"","metadata":{"execution":{"iopub.status.busy":"2021-09-24T13:06:54.422381Z","iopub.execute_input":"2021-09-24T13:06:54.4227Z","iopub.status.idle":"2021-09-24T13:06:54.431689Z","shell.execute_reply.started":"2021-09-24T13:06:54.422661Z","shell.execute_reply":"2021-09-24T13:06:54.430172Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Resample function. Original idea from -> https://www.kaggle.com/boojum/connecting-voxel-spaces\n\ndef 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","metadata":{"execution":{"iopub.status.busy":"2021-09-24T13:06:54.433579Z","iopub.execute_input":"2021-09-24T13:06:54.433965Z","iopub.status.idle":"2021-09-24T13:06:54.443154Z","shell.execute_reply.started":"2021-09-24T13:06:54.433935Z","shell.execute_reply":"2021-09-24T13:06:54.44231Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Setup a simpleITK reader and load a volume as a reference plane. Choose a series that's in the plane you want to export the target series to.\n# We'll use the following series as reference planes since we already know their planes.\n\n# train_dirs[0]/T2w = sagittal plane\n# train_dirs[0]/T1w = axial plane\n# train_dirs[5]/T1wCE = coronal plane\n\nif PLANE == 'axial':\n    train_dir = train_dirs[0]\n    train_series = \"T1w\"\nif PLANE == 'sagittal':\n    train_dir = train_dirs[0]\n    train_series = \"T2w\"   \nif PLANE == 'coronal':\n    train_dir = train_dirs[5]\n    train_series = \"T1wCE\"\n\nreader = sitk.ImageSeriesReader()\nreader.LoadPrivateTagsOn()\n\nref_files = reader.GetGDCMSeriesFileNames(f'{train_path}/{train_dir}/{train_series}')\nreader.SetFileNames(ref_files)\nref_set = reader.Execute()\n     \nref_plane = get_image_plane(ref_set.GetDirection())\nref_array = sitk.GetArrayFromImage(ref_set)\nprint(\"Reference Image Plane: \" + ref_plane)","metadata":{"execution":{"iopub.status.busy":"2021-09-24T13:06:54.444893Z","iopub.execute_input":"2021-09-24T13:06:54.445196Z","iopub.status.idle":"2021-09-24T13:06:56.743269Z","shell.execute_reply.started":"2021-09-24T13:06:54.445166Z","shell.execute_reply":"2021-09-24T13:06:56.741473Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Get a series into an array\ndef get_target_array(study, series):\n\n    target_files = reader.GetGDCMSeriesFileNames(f'{train_path}/{study}/{series}')\n    reader.SetFileNames(target_files)\n    target_set = reader.Execute()\n\n    target_plane = get_image_plane(target_set.GetDirection())\n\n    # Only resample images that aren't already in the reference plane, otherwise just export as-is.\n    if ref_plane != target_plane:\n        target_new = resample(target_set, ref_set)\n    else:\n        target_new = target_set\n\n    return sitk.GetArrayFromImage(target_new)","metadata":{"execution":{"iopub.status.busy":"2021-09-24T13:06:56.744617Z","iopub.execute_input":"2021-09-24T13:06:56.745157Z","iopub.status.idle":"2021-09-24T13:06:56.752292Z","shell.execute_reply.started":"2021-09-24T13:06:56.745116Z","shell.execute_reply":"2021-09-24T13:06:56.751091Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Crunch pixels down to 8 bit\ndef normalize_pixels(pixels):\n    pixels = pixels - np.min(pixels)\n    pixels = pixels / np.max(pixels)\n    pixels = (pixels * 255).astype(np.uint8)\n    return pixels","metadata":{"execution":{"iopub.status.busy":"2021-09-24T13:06:56.754692Z","iopub.execute_input":"2021-09-24T13:06:56.755003Z","iopub.status.idle":"2021-09-24T13:06:56.767941Z","shell.execute_reply.started":"2021-09-24T13:06:56.754972Z","shell.execute_reply":"2021-09-24T13:06:56.766768Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def get_images(study, series):\n    count = 0\n\n    # Get the images in the target series\n    target_array = get_target_array(study, series) \n\n    # Make dirs to save images to\n    if not os.path.exists(f'/kaggle/working/train/{study}'):\n        os.mkdir(f'/kaggle/working/train/{study}')\n    if not os.path.exists(f'/kaggle/working/train/{study}/{series}'):\n        os.mkdir(f'/kaggle/working/train/{study}/{series}')\n\n    # iterate through each image in the array\n    for i in range(0,target_array.shape[0]):\n\n        if len(target_array) > 0:\n            target = target_array[i,:,:]\n            pixels = normalize_pixels(target)\n\n            # Just check the center of the image to ignore the grey border we get with matrices that aren't the same shape (256 vs 194 etc)\n            w = pixels.shape[0] - 100\n            h = pixels.shape[1] - 100\n            pix_mean = np.mean(pixels[100:w,100:h])\n\n            # Export only images that have 'some' pixels. Not sure of the exact threshold here. 30 seems reasonable.\n            # Maybe this should be a count of pixels that have any value other than the background value (~0) instead of mean since it will be different for each series type\n            if pix_mean > 20:\n                filename = f'/kaggle/working/train/{study}/{series}/{i}.jpg'\n                cv2.imwrite(filename, pixels)\n                count += 1","metadata":{"execution":{"iopub.status.busy":"2021-09-24T13:06:56.769729Z","iopub.execute_input":"2021-09-24T13:06:56.770113Z","iopub.status.idle":"2021-09-24T13:06:56.783267Z","shell.execute_reply.started":"2021-09-24T13:06:56.770082Z","shell.execute_reply":"2021-09-24T13:06:56.781951Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Iterate through each study directory\nfor study in train_dirs:\n    get_images(study, SERIES)","metadata":{"execution":{"iopub.status.busy":"2021-09-24T13:06:56.787281Z","iopub.execute_input":"2021-09-24T13:06:56.787702Z","iopub.status.idle":"2021-09-24T13:06:57.404804Z","shell.execute_reply.started":"2021-09-24T13:06:56.787658Z","shell.execute_reply":"2021-09-24T13:06:57.403269Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Zip it up\nshutil.make_archive(f'{SERIES}_{PLANE}', 'zip', '/kaggle/working/train')","metadata":{"execution":{"iopub.status.busy":"2021-09-24T13:06:57.406765Z","iopub.execute_input":"2021-09-24T13:06:57.407182Z","iopub.status.idle":"2021-09-24T13:06:57.445064Z","shell.execute_reply.started":"2021-09-24T13:06:57.407131Z","shell.execute_reply":"2021-09-24T13:06:57.444051Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Create a link to download the zip\nos.chdir(r'/kaggle/working')\n%cd /kaggle/working\nFileLink(f'{SERIES}_{PLANE}.zip')","metadata":{"execution":{"iopub.status.busy":"2021-09-24T13:06:57.446417Z","iopub.execute_input":"2021-09-24T13:06:57.446753Z","iopub.status.idle":"2021-09-24T13:06:57.455951Z","shell.execute_reply.started":"2021-09-24T13:06:57.446721Z","shell.execute_reply":"2021-09-24T13:06:57.454458Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#### Some of my other MR notebooks\n- Tumor Object Detection -> https://www.kaggle.com/davidbroberts/brain-tumor-object-detection\n- Determining MR Slice Orientation -> https://www.kaggle.com/davidbroberts/determining-mr-slice-orientation\n- Determining DICOM Image Order -> https://www.kaggle.com/davidbroberts/determining-dicom-image-order/\n- Determining MR image planes -> https://www.kaggle.com/davidbroberts/determining-mr-image-planes\n- Reference Lines on MR images -> https://www.kaggle.com/davidbroberts/mr-reference-lines\n- Manual VOI LUT on MR images -> https://www.kaggle.com/davidbroberts/manual-voi-lut-on-mr-images\n- Standardizing MR images -> https://www.kaggle.com/davidbroberts/standardizing-mr-images","metadata":{}}]}