Resample multiple MRI sequences of the kaggle MICCAI challenge to axial view

resample_mri.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "6076bb78",
   "metadata": {},
   "source": [
    "# Resample a 3D Image"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "id": "2ad4d83f",
   "metadata": {},
   "outputs": [],
   "source": [
    "import SimpleITK as sitk\n",
    "import numpy as np\n",
    "from matplotlib import pyplot as plt\n",
    "from ipywidgets import interact, fixed"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "id": "0e7a7b31",
   "metadata": {},
   "outputs": [],
   "source": [
    " def read_dicom_series(fn:str):\n",
    "    \"load a DICOM series as `sitk.Image`\"\n",
    "    SeriesReader = sitk.ImageSeriesReader()\n",
    "    dicom_names = SeriesReader.GetGDCMSeriesFileNames(str(fn))\n",
    "    SeriesReader.SetFileNames(dicom_names)\n",
    "    im = SeriesReader.Execute()\n",
    "    return sitk.Cast(im, sitk.sitkInt16)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "id": "842b0013",
   "metadata": {},
   "outputs": [],
   "source": [
    "def resize_image(im, new_size=(128, 128, 64)):\n",
    "    \"Resample a sitk.Image to `new_size` maintaining oriantation and direction\"\n",
    "    # get old image meta data\n",
    "    old_size = im.GetSize()\n",
    "    old_spacing = im.GetSpacing()\n",
    "    old_orig = im.GetOrigin()\n",
    "    old_dir = im.GetDirection()\n",
    "    old_pixel_id = im.GetPixelIDValue()\n",
    "    \n",
    "    # calculate new spacing\n",
    "    new_spacing = [spc * (old_sz/new_sz) for spc, old_sz, new_sz in zip(old_spacing, old_size, new_size)]\n",
    "    \n",
    "    # create reference plane to resample image\n",
    "    ref_image = sitk.Image(new_size, old_pixel_id)\n",
    "    ref_image.SetSpacing(new_spacing)\n",
    "    ref_image.SetOrigin(old_orig)\n",
    "    ref_image.SetDirection(old_dir)\n",
    "    \n",
    "    # resample image to `new_size`\n",
    "    return sitk.Resample(im, ref_image)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 64,
   "id": "4a56e042",
   "metadata": {},
   "outputs": [],
   "source": [
    "def resample_to_axial(image, spacing=None):\n",
    "    \"resample an image to direction (1,0,0,0,1,0,0,0,1)\"\n",
    "    \n",
    "    # define the Euler Transformation for Image Rotation\n",
    "    euler3d = sitk.Euler3DTransform() # define transform for rotation of the image\n",
    "    image_center = (np.array(image.GetSize())/2.0)  # set rotation center to image center\n",
    "    image_center_as_sitk_point = image.TransformContinuousIndexToPhysicalPoint(image_center)\n",
    "    euler3d.SetCenter(image_center_as_sitk_point)\n",
    "    \n",
    "    # get index of volume edges\n",
    "    w,h,d = image.GetSize()\n",
    "    extreme_points = [(0, 0, 0), (w, 0, 0), (0, h, 0), (0, 0, d), \n",
    "                      (w, h, 0), (w, 0, d), (0, h, d), (w, h, d)]\n",
    "    # transform edges to physical points in the global coordinate system\n",
    "    extreme_points = [image.TransformIndexToPhysicalPoint(pnt) for pnt in extreme_points]\n",
    "    inv_euler3d = euler3d.GetInverse()\n",
    "    extreme_points_transformed = [inv_euler3d.TransformPoint(pnt) for pnt in extreme_points]\n",
    "    \n",
    "    # get new min and max coordinates of image edges\n",
    "    min_x = min(extreme_points_transformed)[0]\n",
    "    min_y = min(extreme_points_transformed, key=lambda p: p[1])[1]\n",
    "    min_z = min(extreme_points_transformed, key=lambda p: p[2])[2]\n",
    "        \n",
    "    max_x = max(extreme_points_transformed)[0]\n",
    "    max_y = max(extreme_points_transformed, key=lambda p: p[1])[1]\n",
    "    max_z = max(extreme_points_transformed, key=lambda p: p[2])[2]\n",
    "        \n",
    "    # define new direction\n",
    "    # take spacing and size from original image\n",
    "    # calculate new origin from extree points\n",
    "    output_spacing = spacing if spacing else image.GetSpacing()\n",
    "    output_direction = tuple(np.identity(3).flatten())\n",
    "    output_origin = [min_x, min_y, min_z]\n",
    "    output_size = [int((max_x-min_x)/output_spacing[0]), \n",
    "                   int((max_y-min_y)/output_spacing[1]), \n",
    "                   int((max_z-min_z)/output_spacing[2])]\n",
    "        \n",
    "    resampled_image = sitk.Resample(image, output_size, euler3d, sitk.sitkLinear, output_origin, output_spacing, output_direction)\n",
    "    return resampled_image "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 66,
   "id": "36277436",
   "metadata": {},
   "outputs": [],
   "source": [
    "def show_image(image, slice_id):\n",
    "    img_slice = sitk.GetArrayViewFromImage(image)[slice_id, :, :]\n",
    "    plt.imshow(img_slice, cmap='bone')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 67,
   "id": "c4337437-ab0e-4d0b-990a-91112f05004f",
   "metadata": {},
   "outputs": [],
   "source": [
    "def read_and_resample(fn, plot=False):\n",
    "    im = read_dicom_series(fn)\n",
    "    im_resampled = resample_to_axial(im)\n",
    "    if plot:\n",
    "        interact(show_image, image=fixed(im_resampled), slice_id = (0, im_resampled.GetDepth()-1));\n",
    "    else: return im_resampled"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 68,
   "id": "1be20060-2a44-4b09-83b1-9c34b8b0a2b0",
   "metadata": {},
   "outputs": [],
   "source": [
    "i = 114"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 69,
   "id": "33f79c71-aa22-4f0d-972c-a76e0cd7deaa",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "bcc4dcf7094b4acaa0ec0a7dcd8341a9",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(IntSlider(value=13, description='slice_id', max=27), Output()), _dom_classes=('widget-in…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "4c8b6e8b8d2c4df884a1b39da74c5986",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(IntSlider(value=13, description='slice_id', max=27), Output()), _dom_classes=('widget-in…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "5c94bdc2956a45089c60cb997ee7cd43",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(IntSlider(value=95, description='slice_id', max=191), Output()), _dom_classes=('widget-i…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "ae5134ba104b4a1da18729fb8791a132",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(IntSlider(value=13, description='slice_id', max=27), Output()), _dom_classes=('widget-in…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "flair = read_and_resample(f'images/{i:05}/FLAIR', plot=True)\n",
    "t1w = read_and_resample(f'images/{i:05}/T1w', plot=True)\n",
    "t2w = read_and_resample(f'images/{i:05}/T1wCE', plot=True)\n",
    "t1wce = read_and_resample(f'images/{i:05}/T2w', plot=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b51aed08-a5a4-4f69-81cc-4215c2c9d877",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}