{
 "metadata": {
  "name": "",
  "signature": "sha256:5b44004e97625e733a6588149bed4fc0c9eda5b70ed3c9420b3c41b63e6842a0"
 },
 "name": "smoothing",
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# - compatibility with Python 3\n",
      "from __future__ import print_function  # print('me') instead of print 'me'\n",
      "from __future__ import division  # 1/2 == 0.5, not 0"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# - show figures inside the notebook\n",
      "%matplotlib inline"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# - import common modules\n",
      "import numpy as np  # the Python array package\n",
      "import matplotlib.pyplot as plt  # the Python plotting package"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# - set gray colormap and nearest neighbor interpolation by default\n",
      "plt.rcParams['image.cmap'] = 'gray'\n",
      "plt.rcParams['image.interpolation'] = 'nearest'"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Smoothing and convolution"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Start by smoothing in 1D.  We first load our new fixed functional image, and get the first volume:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# - get first volume from fixed functional\n",
      "import nibabel as nib\n",
      "img = nib.load('fds114_sub009_t2r1.nii')\n",
      "data = img.get_data()\n",
      "vol0 = data[..., 0]"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Now get a slice over the third dimension:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# - get slice over third dimension\n",
      "mid_slice = vol0[:, :, 14]\n",
      "plt.imshow(mid_slice)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Get one line of this slice, from left to right, from the middle of the brain:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# - get middle line from middle slice\n",
      "mid_line = mid_slice[:, 25]\n",
      "plt.plot(mid_line)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "We are going to start by smoothing this line (1D smooth)."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Copy the `fwhm2sigma` function from the [smoothing as convolution](http://perrin.dynevor.org/smoothing_as_convolution.html) notebook:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# fwhm2sigma"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# - check fwhm2sigma\n",
      "assert np.allclose(fwhm2sigma(4), 1.698643600576)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "We are going to smooth by 8 pixels FWHM, so we need the corresponding `sigma`:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# sigma for FWHM 8"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "We use a probability density function of a normal Gaussian distribution for our kernel:\n",
      "\n",
      "$$\n",
      "f(x, \\mu, \\sigma) = p e^q\n",
      "$$\n",
      "\n",
      "where:\n",
      "\n",
      "$$\n",
      "p = \\frac{1}{\\sqrt{2\\pi\\sigma^2}} \\\\\n",
      "q = -\\frac{(x-\\mu)^2}{2\\sigma^2}\n",
      "$$\n",
      "\n",
      "See : https://en.wikipedia.org/wiki/Normal_distribution\n",
      "\n",
      "To make our smoothing kernel, we can sample the PDF, for a given mu and sigma.\n",
      "\n",
      "We can use our old friend `scipy.stats` to sample the PDF."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# - import the normal distribution function from scipy.stats\n",
      "from scipy.stats import norm\n",
      "norm_pdf = norm.pdf"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Now we need to make our kernel.\n",
      "\n",
      "* Work out the +/- limit for the kernel x values with `limit = round(sigma * 4)`\n",
      "* Make a vector of integers `x` to sample the PDF from `-limit` to `limit` - *inclusive*;\n",
      "* Sample the PDF at these x values, for $\\mu = 0$ and given $\\sigma$, to get the `kernel` vector;\n",
      "* Plot the x values against the y values;\n",
      "* Work out how many indices from the start of the `kernel` vector, to the value corresponding to x=0; (this is the `kernel_offset`."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Work out the +/- limit for the kernel x values;\n",
      "# Make a vector `x` to sample the PDF;\n",
      "# Get `kernel` vector by sampling the PDF at these x values (mu=0);\n",
      "# Work out `kernel_offset`."
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Now we can convolve our image line with the `kernel`.  We take the result starting at `kernel_offset`:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Convolve image line with kernel\n",
      "# Slice out the result we want using kernel_offset\n",
      "# Plot smoothed with unsmoothed line"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Now we are going to smooth a whole slice.\n",
      "\n",
      "* Make a new array the same shape as `mid_slice` called `smoothed_slice`;\n",
      "* Loop over the first dimension of `mid_slice`, and fill each line of `smoothed_slice` with the smoothed version;\n",
      "* Loop over the second dimension of `smoothed_slice` and replace each line with the smoothed version;\n",
      "* Show."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Make a new array `smoothed_slice` for smoothed output\n",
      "# Loop over first dimension applying kernel as above\n",
      "# Loop over second dimension of smoothed image, applying kernel\n",
      "# Show with `imshow`"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Now time to do the same thing with `scipy.ndimage.gaussian_filter`."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# - import gaussian_filter\n",
      "from scipy.ndimage import gaussian_filter"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Have a look at the help for `gaussian_filter`.\n",
      "\n",
      "Use `gaussian_filter` to make another smoothed version of `mid_slice`."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Use gaussian filter to smooth `mid_slice`"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Subtract the `gaussian_filter` version from your handcrafted version.\n",
      "\n",
      "Can you explain the difference?"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Subtract two versions of the smoothed slice, show"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}