{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Creating a Printable Model from a 3D Medical Image\n",
    "\n",
    "## A Tutorial on dicom2stl.py\n",
    "\n",
    "[https://github.com/dave3d/dicom2stl](https://github.com/dave3d/dicom2stl)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "![heads](https://github.com/dave3d/dicom2stl/blob/main/examples/Data/head_diagram.jpg?raw=true)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "import SimpleITK as sitk\n",
    "%matplotlib notebook"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Digital Imaging and Communications in Medicine (DICOM)\n",
    "\n",
    "DICOM is the standard for the communication and management of **medical imaging information** and related data.\n",
    "\n",
    "DICOM is most commonly used for storing and transmitting medical images enabling the **integration of medical imaging devices** such as scanners, servers, workstations, printers, network hardware, and **picture archiving and communication systems (PACS)** from multiple manufacturers\n",
    "\n",
    "[https://en.wikipedia.org/wiki/DICOM](https://en.wikipedia.org/wiki/DICOM)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Imaging Modalities\n",
    "\n",
    " * CT (computed tomography)\n",
    " * MRI (magnetic resonance imaging)\n",
    " * ultrasound\n",
    " * X-ray\n",
    " * fluoroscopy\n",
    " * angiography\n",
    " * mammography\n",
    " * breast tomosynthesis\n",
    " * PET (positron emission tomography)\n",
    " * SPECT (single photon emission computed tomography)\n",
    " * Endoscopy\n",
    " * microscopy and whole slide imaging\n",
    " * OCT (optical coherence tomography)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "ct_image = sitk.ReadImage('Data/ct_example.nii.gz')\n",
    "mri_image = sitk.ReadImage('Data/mri_t1_example.nii.gz')\n",
    "\n",
    "import gui\n",
    "gui.MultiImageDisplay(image_list=[ct_image, mri_image], title_list=['CT Head', 'MRI T1 Head'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# CT Houndsfield Units\n",
    "\n",
    "Hounsfield units (HU) are a dimensionless unit universally used in computed tomography (CT) scanning to express CT numbers in a standardized and convenient form. Hounsfield units are obtained from a linear transformation of the measured attenuation coefficients 1\n",
    "\n",
    " * Water is 0 HU\n",
    " * Air is -1000 HU\n",
    " * Very dense bone is 2000 HU\n",
    " * Metal is 3000 HU\n",
    " \n",
    " [Houndsfield Wikipedia page](https://en.wikipedia.org/wiki/Hounsfield_scale)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Image Segmentation\n",
    "\n",
    "The process of partitioning an image into multiple segments.\n",
    "\n",
    "Typically used to locate objects and boundaries in images.\n",
    "\n",
    "We use thresholding (selecting a range of image intesities), but SimpleITK has a variety of algorithms\n",
    "\n",
    "[SimpleITK Notebooks](https://github.com/InsightSoftwareConsortium/SimpleITK-Notebooks/tree/master/Python)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "from myshow import myshow, myshow3d\n",
    "\n",
    "ct_bone = ct_image>200\n",
    "\n",
    "# To visualize the labels image in RGB with needs a image with 0-255 range\n",
    "ct255_image = sitk.Cast(sitk.IntensityWindowing(ct_bone,0,500.0,0.,255.), \n",
    "                        sitk.sitkUInt8)\n",
    "\n",
    "ct255_bone = sitk.Cast(ct_bone, sitk.sitkUInt8)\n",
    "\n",
    "myshow(sitk.LabelOverlay(ct255_image, ct255_bone), \"Basic Thresholding\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Iso-surface extraction\n",
    "\n",
    "Extract a polygonal surface from a 3D image.  The most well known algorithm is Marching Cubes (Lorenson & Cline, SIGGRAPH 1987).  The 2D version is Marching Squares, shown below\n",
    "\n",
    "![Marching Squares](https://github.com/dave3d/dicom2stl/blob/main/examples/Data/marching_squares.png?raw=true)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Marching Cubes\n",
    "\n",
    "And here is the lookup table for Marching Cubes\n",
    "\n",
    "![Marching Cubes](https://github.com/dave3d/dicom2stl/blob/main/examples/Data/marching_cubes.png?raw=true)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# dicom2stl.py processing pipeline\n",
    "\n",
    "SimpleITK image processing pipeline\n",
    "\n",
    " * **Shrink** the volume to 256^3\n",
    " * Apply **anisotripic smoothing**\n",
    " * **Threshold**\n",
    "    - Preset tissue types: skin, bone, fat, soft tissue\n",
    "    - User specified iso-value\n",
    " * **Median filter**\n",
    " * **Pad** the volume with black\n",
    " \n",
    "VTK mesh pipeline\n",
    "\n",
    " * Run **Marching Cubes** to extract surface\n",
    " * Apply **CleanMesh** filter to merge vertices\n",
    " * Apply **SmoothMesh** filter\n",
    " * Run **polygon reduction**\n",
    " * Write STL\n",
    " \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "import itkwidgets\n",
    "head = sitk.ReadImage(\"Data/ct_head.nii.gz\")\n",
    "itkwidgets.view(head)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "import sys, os\n",
    "\n",
    "# download dicom2stl if it's not here already\n",
    "if not os.path.isdir('dicom2stl'):\n",
    "    !{'git clone https://github.com/dave3d/dicom2stl.git'}\n",
    "    \n",
    "!{sys.executable} dicom2stl/dicom2stl.py -h"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "!{sys.executable} dicom2stl/dicom2stl.py -i 400 -o bone.stl Data/ct_head.nii.gz"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "from dicom2stl.utils import vtkutils\n",
    "mesh = vtkutils.readMesh('bone.stl')\n",
    "itkwidgets.view(head, geometries=[mesh])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "celltoolbar": "Slideshow",
  "kernelspec": {
   "display_name": "Python 3",
   "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.7.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}