{"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":"In this notebook, we will convert the dicom files into nifti format. Also, we will use the SRI24 dataset to resample the images to a consistent orientation. This is, within the same modality - say T2w, the sequene for different patient appear in different orientation (axial, sagittal etc). We will use the SRI24 axial image as the referece image to resample our data.","metadata":{}},{"cell_type":"markdown","source":"The following notebooks/tutorial were of great help.\n\n1. https://www.kaggle.com/boojum/connecting-voxel-spaces/\n2. https://simpleitk.readthedocs.io/en/master/link_examples.html","metadata":{}},{"cell_type":"markdown","source":"# Imports and paths","metadata":{}},{"cell_type":"code","source":"import warnings\nwarnings.filterwarnings('ignore')","metadata":{"execution":{"iopub.status.busy":"2021-08-29T04:14:19.589621Z","iopub.execute_input":"2021-08-29T04:14:19.589976Z","iopub.status.idle":"2021-08-29T04:14:19.594625Z","shell.execute_reply.started":"2021-08-29T04:14:19.589943Z","shell.execute_reply":"2021-08-29T04:14:19.593664Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import os\nfrom tqdm import tqdm\n\nimport nibabel as nib\nimport SimpleITK as sitk\n\nfrom fastai.medical.imaging import *\nfrom fastai.vision.all import *","metadata":{"execution":{"iopub.status.busy":"2021-08-29T04:14:19.778991Z","iopub.execute_input":"2021-08-29T04:14:19.779377Z","iopub.status.idle":"2021-08-29T04:14:19.784883Z","shell.execute_reply.started":"2021-08-29T04:14:19.779339Z","shell.execute_reply":"2021-08-29T04:14:19.783796Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"path_train =  Path('../input/rsna-miccai-brain-tumor-radiogenomic-classification/train')","metadata":{"execution":{"iopub.status.busy":"2021-08-29T04:14:22.790856Z","iopub.execute_input":"2021-08-29T04:14:22.791234Z","iopub.status.idle":"2021-08-29T04:14:22.795847Z","shell.execute_reply.started":"2021-08-29T04:14:22.791199Z","shell.execute_reply":"2021-08-29T04:14:22.794794Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"path_train_t2w,path_train_t1wce,path_train_t1w,path_train_flair = [],[],[],[]\nfor each in path_train.ls():\n    path_train_t2w.append(each.ls()[0])\n    path_train_t1wce.append(each.ls()[1])\n    path_train_t1w.append(each.ls()[2])\n    path_train_flair.append(each.ls()[3])","metadata":{"execution":{"iopub.status.busy":"2021-08-29T04:14:22.964929Z","iopub.execute_input":"2021-08-29T04:14:22.965307Z","iopub.status.idle":"2021-08-29T04:14:25.081670Z","shell.execute_reply.started":"2021-08-29T04:14:22.965276Z","shell.execute_reply":"2021-08-29T04:14:25.080737Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Functions","metadata":{}},{"cell_type":"code","source":"def get_array(fn):\n    \"opens .nii file and return the array\"\n    img = sitk.ReadImage(str(fn))\n    imgd = sitk.GetArrayFromImage(img)\n    return imgd\n\ndef plot_slice(imgd, sli):\n    \"given an image of shape slices x height x width, plots a slice\"\n    plt.imshow(imgd[sli], cmap='gray')\n    plt.axis('off')","metadata":{"execution":{"iopub.status.busy":"2021-08-29T04:16:25.085370Z","iopub.execute_input":"2021-08-29T04:16:25.085870Z","iopub.status.idle":"2021-08-29T04:16:25.091687Z","shell.execute_reply.started":"2021-08-29T04:16:25.085822Z","shell.execute_reply":"2021-08-29T04:16:25.090752Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def dicom2nifti(image_dir, save=True):\n    \"given a dicom directory, loads them into single file and can save it as .nii file\"\n    reader = sitk.ImageSeriesReader()\n    reader.LoadPrivateTagsOn()\n    filenamesDICOM = reader.GetGDCMSeriesFileNames(str(image_dir))\n    reader.SetFileNames(filenamesDICOM)\n    img = reader.Execute()\n    img = sitk.Cast(img, sitk.sitkFloat32)\n    \n    if save:\n        sitk.WriteImage(img, f'/kaggle/working/T2w/{image_dir.parent.name}.nii')\n    else:\n        return img","metadata":{"execution":{"iopub.status.busy":"2021-08-29T04:14:32.835991Z","iopub.execute_input":"2021-08-29T04:14:32.836572Z","iopub.status.idle":"2021-08-29T04:14:32.843399Z","shell.execute_reply.started":"2021-08-29T04:14:32.836522Z","shell.execute_reply":"2021-08-29T04:14:32.842595Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def resample_nifti(image_dir, ref_image, fn, save=True):\n    \"resample using a reference image\"\n\n    image = dicom2nifti(image_dir, save=False)\n    \n    initial_transform = sitk.CenteredTransformInitializer(ref_image, \n                                                          image, \n                                                          sitk.Euler3DTransform(), \n                                                          sitk.CenteredTransformInitializerFilter.GEOMETRY)\n\n    resampler = sitk.ResampleImageFilter()\n    resampler.SetReferenceImage(ref_image)\n    resampler.SetInterpolator(sitk.sitkLinear)\n    resampler.SetTransform(initial_transform)\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    \n    if save:\n        sitk.WriteImage(resamped_image, fn)\n\n    return resamped_image","metadata":{"execution":{"iopub.status.busy":"2021-08-29T04:17:32.503388Z","iopub.execute_input":"2021-08-29T04:17:32.503821Z","iopub.status.idle":"2021-08-29T04:17:32.511108Z","shell.execute_reply.started":"2021-08-29T04:17:32.503784Z","shell.execute_reply":"2021-08-29T04:17:32.510336Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#code from simpleitk examples\ndef threshold_based_crop_and_bg_median(image):\n    '''\n    Use Otsu's threshold estimator to separate background and foreground. In medical imaging the background is\n    usually air. Then crop the image using the foreground's axis aligned bounding box and compute the background \n    median intensity.\n    Args:\n        image (SimpleITK image): An image where the anatomy and background intensities form a bi-modal distribution\n                                 (the assumption underlying Otsu's method.)\n    Return:\n        Cropped image based on foreground's axis aligned bounding box.\n        Background median intensity value.\n    '''\n    # Set pixels that are in [min_intensity,otsu_threshold] to inside_value, values above otsu_threshold are\n    # set to outside_value. The anatomy has higher intensity values than the background, so it is outside.\n    inside_value = 0\n    outside_value = 255\n    bin_image = sitk.OtsuThreshold(image, inside_value, outside_value)\n\n    # Get the median background intensity\n    label_intensity_stats_filter = sitk.LabelIntensityStatisticsImageFilter()\n    label_intensity_stats_filter.SetBackgroundValue(outside_value)\n    label_intensity_stats_filter.Execute(bin_image,image)\n    bg_median = label_intensity_stats_filter.GetMedian(inside_value)\n    \n    # Get the bounding box of the anatomy\n    label_shape_filter = sitk.LabelShapeStatisticsImageFilter()    \n    label_shape_filter.Execute(bin_image)\n    bounding_box = label_shape_filter.GetBoundingBox(outside_value)\n    # The bounding box's first \"dim\" entries are the starting index and last \"dim\" entries the size\n    return bg_median, sitk.RegionOfInterest(image, bounding_box[int(len(bounding_box)/2):], bounding_box[0:int(len(bounding_box)/2)])","metadata":{"execution":{"iopub.status.busy":"2021-08-29T04:17:49.358668Z","iopub.execute_input":"2021-08-29T04:17:49.359185Z","iopub.status.idle":"2021-08-29T04:17:49.367377Z","shell.execute_reply.started":"2021-08-29T04:17:49.359141Z","shell.execute_reply":"2021-08-29T04:17:49.366469Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Let's see what we are talking about ","metadata":{}},{"cell_type":"markdown","source":"Let's open one and see how it looks like","metadata":{}},{"cell_type":"code","source":"samp = path_train_t2w[0]","metadata":{"execution":{"iopub.status.busy":"2021-08-29T04:19:21.767523Z","iopub.execute_input":"2021-08-29T04:19:21.767886Z","iopub.status.idle":"2021-08-29T04:19:21.771624Z","shell.execute_reply.started":"2021-08-29T04:19:21.767851Z","shell.execute_reply":"2021-08-29T04:19:21.770861Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"samp_img = dicom2nifti(samp, False)\nsamp_imgd = sitk.GetArrayFromImage(samp_img)","metadata":{"execution":{"iopub.status.busy":"2021-08-29T04:20:35.350453Z","iopub.execute_input":"2021-08-29T04:20:35.350823Z","iopub.status.idle":"2021-08-29T04:20:37.398244Z","shell.execute_reply.started":"2021-08-29T04:20:35.350785Z","shell.execute_reply":"2021-08-29T04:20:37.397551Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"samp_imgd.shape","metadata":{"execution":{"iopub.status.busy":"2021-08-29T04:20:51.395273Z","iopub.execute_input":"2021-08-29T04:20:51.395776Z","iopub.status.idle":"2021-08-29T04:20:51.401134Z","shell.execute_reply.started":"2021-08-29T04:20:51.395741Z","shell.execute_reply":"2021-08-29T04:20:51.400229Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plot_slice(samp_imgd, 100)","metadata":{"execution":{"iopub.status.busy":"2021-08-29T04:21:05.255356Z","iopub.execute_input":"2021-08-29T04:21:05.255728Z","iopub.status.idle":"2021-08-29T04:21:05.391032Z","shell.execute_reply.started":"2021-08-29T04:21:05.255697Z","shell.execute_reply":"2021-08-29T04:21:05.389927Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Now lets resample using a SIR24 as referece image","metadata":{}},{"cell_type":"code","source":"ref_image = sitk.ReadImage(\"../input/sri24-dataset/sri24/late.nii\", sitk.sitkFloat32)","metadata":{"execution":{"iopub.status.busy":"2021-08-29T04:23:55.383374Z","iopub.execute_input":"2021-08-29T04:23:55.383871Z","iopub.status.idle":"2021-08-29T04:23:55.420733Z","shell.execute_reply.started":"2021-08-29T04:23:55.383838Z","shell.execute_reply":"2021-08-29T04:23:55.419558Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"samp_resamp = resample_nifti(samp, ref_image, \"\", False)\nsamp_resampd = sitk.GetArrayFromImage(samp_resamp)","metadata":{"execution":{"iopub.status.busy":"2021-08-29T04:25:03.809907Z","iopub.execute_input":"2021-08-29T04:25:03.810284Z","iopub.status.idle":"2021-08-29T04:25:05.930304Z","shell.execute_reply.started":"2021-08-29T04:25:03.810249Z","shell.execute_reply":"2021-08-29T04:25:05.929241Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We have reshaped the samp images to the same dimensions as the reference image","metadata":{}},{"cell_type":"code","source":"sitk.GetArrayFromImage(ref_image).shape, samp_resampd.shape","metadata":{"execution":{"iopub.status.busy":"2021-08-29T04:25:11.488245Z","iopub.execute_input":"2021-08-29T04:25:11.488631Z","iopub.status.idle":"2021-08-29T04:25:11.511503Z","shell.execute_reply.started":"2021-08-29T04:25:11.488593Z","shell.execute_reply":"2021-08-29T04:25:11.510588Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"And we have resampled to the axial view! ","metadata":{}},{"cell_type":"code","source":"plot_slice(samp_resampd, 98)","metadata":{"execution":{"iopub.status.busy":"2021-08-29T04:26:37.182981Z","iopub.execute_input":"2021-08-29T04:26:37.183378Z","iopub.status.idle":"2021-08-29T04:26:37.276195Z","shell.execute_reply.started":"2021-08-29T04:26:37.183342Z","shell.execute_reply":"2021-08-29T04:26:37.275227Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We can also crop them to the region of interest","metadata":{}},{"cell_type":"code","source":"_,samp_resamp_cropped = threshold_based_crop_and_bg_median(samp_resamp)\nsamp_resamp_croppedd = sitk.GetArrayFromImage(samp_resamp_cropped)","metadata":{"execution":{"iopub.status.busy":"2021-08-29T04:28:20.369780Z","iopub.execute_input":"2021-08-29T04:28:20.370130Z","iopub.status.idle":"2021-08-29T04:28:21.339284Z","shell.execute_reply.started":"2021-08-29T04:28:20.370101Z","shell.execute_reply":"2021-08-29T04:28:21.338351Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"samp_resamp_croppedd.shape","metadata":{"execution":{"iopub.status.busy":"2021-08-29T04:28:47.450608Z","iopub.execute_input":"2021-08-29T04:28:47.451136Z","iopub.status.idle":"2021-08-29T04:28:47.457325Z","shell.execute_reply.started":"2021-08-29T04:28:47.451103Z","shell.execute_reply":"2021-08-29T04:28:47.456442Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plot_slice(samp_resamp_croppedd, 80)","metadata":{"execution":{"iopub.status.busy":"2021-08-29T04:28:55.819818Z","iopub.execute_input":"2021-08-29T04:28:55.820197Z","iopub.status.idle":"2021-08-29T04:28:55.908052Z","shell.execute_reply.started":"2021-08-29T04:28:55.820166Z","shell.execute_reply":"2021-08-29T04:28:55.907385Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"As you can see the first few slices have no region of interest. By cropping to the ROI, we remove the not useful slices. ","metadata":{}},{"cell_type":"code","source":"plot_slice(samp_resampd, 0)","metadata":{"execution":{"iopub.status.busy":"2021-08-29T04:31:01.351444Z","iopub.execute_input":"2021-08-29T04:31:01.351904Z","iopub.status.idle":"2021-08-29T04:31:01.436097Z","shell.execute_reply.started":"2021-08-29T04:31:01.351868Z","shell.execute_reply":"2021-08-29T04:31:01.434965Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plot_slice(samp_resamp_croppedd, 0)","metadata":{"execution":{"iopub.status.busy":"2021-08-29T04:31:55.581907Z","iopub.execute_input":"2021-08-29T04:31:55.582305Z","iopub.status.idle":"2021-08-29T04:31:55.665059Z","shell.execute_reply.started":"2021-08-29T04:31:55.582265Z","shell.execute_reply":"2021-08-29T04:31:55.663994Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Resample all T2w images","metadata":{}},{"cell_type":"code","source":"!mkdir t2w_preproc_v2","metadata":{"execution":{"iopub.status.busy":"2021-08-28T03:36:49.072761Z","iopub.execute_input":"2021-08-28T03:36:49.073317Z","iopub.status.idle":"2021-08-28T03:36:52.124824Z","shell.execute_reply.started":"2021-08-28T03:36:49.073271Z","shell.execute_reply":"2021-08-28T03:36:52.123711Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"ref_image = sitk.ReadImage('../input/sri24-dataset/sri24/spgr.nii', sitk.sitkFloat32)","metadata":{"execution":{"iopub.status.busy":"2021-08-28T03:36:52.126426Z","iopub.execute_input":"2021-08-28T03:36:52.126896Z","iopub.status.idle":"2021-08-28T03:36:52.319179Z","shell.execute_reply.started":"2021-08-28T03:36:52.126859Z","shell.execute_reply":"2021-08-28T03:36:52.318131Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"for fn in tqdm(path_train_t2w, total=len(path_train_t2w)):\n    pat_id = str(fn).split('/')[-2]\n    final_fn = f\"/kaggle/working/t2w_preproc_v2/{pat_id}.nii.gz\"\n    resample_nifti(fn, ref_image, final_fn, True)\n","metadata":{"execution":{"iopub.status.busy":"2021-08-28T03:37:52.467486Z","iopub.execute_input":"2021-08-28T03:37:52.467915Z","iopub.status.idle":"2021-08-28T03:38:10.310368Z","shell.execute_reply.started":"2021-08-28T03:37:52.467876Z","shell.execute_reply":"2021-08-28T03:38:10.308752Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from zipfile import ZipFile\nimport os\nimport shutil\n\n\n# iterate over all the files in directory\nfor folderName, subfolders, filenames in os.walk(f'/kaggle/working/t2w_preproc_v2'):\n    # create a ZipFile object\n    #print(folderName)\n    with ZipFile(folderName.split('/')[-1] + '.zip', 'w') as zipObj:\n        for filename in filenames:\n            # create complete filepath of file in directory\n            filePath = os.path.join(folderName, filename)\n            # add file to zip\n            zipObj.write(filePath, os.path.basename(filePath))\n            # delete the file to open space\n            os.remove(filePath)\nshutil.rmtree(f\"/kaggle/working/t2w_preproc_v2\")\n","metadata":{"execution":{"iopub.status.busy":"2021-08-28T03:36:35.426229Z","iopub.status.idle":"2021-08-28T03:36:35.426812Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}