21_Transform_and_Resampling.ipynb

{
 "metadata": {
  "name": ""
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import SimpleITK as sitk\n",
      "import numpy as np"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def cShow(img1, img2, title=\"\", margin=0.05, expand=1):\n",
      "    simg1 = sitk.Cast(sitk.RescaleIntensity(img1), sitk.sitkUInt8)\n",
      "    simg2 = sitk.Cast(sitk.RescaleIntensity(img2), sitk.sitkUInt8)\n",
      "    cimg = sitk.Compose(simg1, simg2, simg1/2.+simg2/2.)\n",
      "    if expand != 1:\n",
      "        cimg = sitk.Expand(cimg, [expand]*3)\n",
      "    del simg1, simg2\n",
      "\n",
      "    aimg = sitk.GetArrayFromImage(cimg)\n",
      "    \n",
      "    \n",
      "    if aimg.ndim == 3:\n",
      "        ysize,xsize,c = aimg.shape\n",
      "        if not c in (3,4):\n",
      "            aimg = aimg[c//2,:,:]\n",
      "    else:\n",
      "        ysize,xsize = aimg.shape\n",
      "\n",
      "        \n",
      "    dpi=80\n",
      "    \n",
      "    # Make a figure big enough to accomodate an axis of xpixels by ypixels\n",
      "    # as well as the ticklabels, etc...\n",
      "    figsize = (1 + margin) * ysize / dpi, (1 + margin) * xsize / dpi\n",
      "\n",
      "    fig = plt.figure(figsize=figsize, dpi=dpi)\n",
      "    # Make the axis the right size...\n",
      "    ax = fig.add_axes([margin, margin, 1 - 2*margin, 1 - 2*margin])\n",
      "    \n",
      "    extent = (0,xsize*cimg.GetSpacing()[1],0,ysize*cimg.GetSpacing()[0])\n",
      "    t = ax.imshow(aimg,extent=extent)\n",
      "    \n",
      "    if aimg.ndim == 2:\n",
      "        t.set_cmap(\"gray\")\n",
      "    \n",
      "    if(title):\n",
      "        plt.title(title)\n",
      "    "
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Load sample data from ANTS registration.\n"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%cd ~/src/SimpleITK-Notebooks/\n",
      "\n",
      "img1 = sitk.Cast(sitk.ReadImage(\"Data/B1.tiff\"), sitk.sitkFloat32)\n",
      "img2 = sitk.Cast(sitk.ReadImage(\"Data/B2.tiff\"), sitk.sitkFloat32)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "cShow(img1,img2, title=\"Pre-Regstration\",expand=3)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "<code>\n",
      " antsRegistration -d 2 -r [ B1.tiff, B2.tiff ,1]\n",
      "                        -m mattes[  B1.tiff, B2.tiff , 1 , 32, regular, 0.1 ]\n",
      "                         -t affine[ 0.05 ]\n",
      "                         -c [100x0x0,1.e-8,20]\n",
      "                        -s 4x2x1\n",
      "                        -f 6x4x2 -l 1\n",
      "                        -m mattes[  B1.tiff, B2.tiff , 1 , 32 ]\n",
      "                         -t syn[ .3, 3, 0.0 ]\n",
      "                         -c [10x0x0,1.e-8,20]\n",
      "                        -s 2x1x0\n",
      "                        -f 4x2x1 -l 1 -u 1 -z 1\n",
      "                        --write-composite-transform 1\n",
      "\t\t        -o [${nm},${nm}_diff.nii.gz,${nm}_inv.nii.gz]\n",
      "</code>"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "tx = sitk.ReadTransform(\"Data/B1_fixed_B2_moving0GenericAffine.mat\")\n",
      "print tx"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "img2_a = sitk.Resample(img2, img1, tx, sitk.sitkNearestNeighbor, sitk.sitkFloat32)\n",
      "cShow(img1,img2_a, \"Affine Registered\", expand=3)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "tx_comp = sitk.ReadTransform(\"Data/B1_fixed_B2_movingComposite.h5\")\n",
      "img2_a = sitk.Resample(img2, img1, tx_comp, sitk.sitkGaussian, sitk.sitkFloat32)\n",
      "cShow(img1,img2_a, \"Deformed Registered\", expand=3)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}