1. pna2015
  2. day12
In [ ]:
# - compatibility with Python 3
from __future__ import print_function  # print('me') instead of print 'me'
from __future__ import division  # 1/2 == 0.5, not 0
In [ ]:
# - show figures inside the notebook
%matplotlib inline
In [ ]:
# - import common modules
import numpy as np  # the Python array package
import matplotlib.pyplot as plt  # the Python plotting package
In [ ]:
# - set gray colormap and nearest neighbor interpolation by default
plt.rcParams['image.cmap'] = 'gray'
plt.rcParams['image.interpolation'] = 'nearest'

Smoothing and convolution

Start by smoothing in 1D. We first load our new fixed functional image, and get the first volume:

In [ ]:
# - get first volume from fixed functional
import nibabel as nib
img = nib.load('fds114_sub009_t2r1.nii')
data = img.get_data()
vol0 = data[..., 0]

Now get a slice over the third dimension:

In [ ]:
# - get slice over third dimension
mid_slice = vol0[:, :, 14]
plt.imshow(mid_slice)

Get one line of this slice, from left to right, from the middle of the brain:

In [ ]:
# - get middle line from middle slice
mid_line = mid_slice[:, 25]
plt.plot(mid_line)

We are going to start by smoothing this line (1D smooth).

Copy the fwhm2sigma function from the smoothing as convolution notebook:

In [ ]:
# fwhm2sigma
In [ ]:
# - check fwhm2sigma
assert np.allclose(fwhm2sigma(4), 1.698643600576)

We are going to smooth by 8 pixels FWHM, so we need the corresponding sigma:

In [ ]:
# sigma for FWHM 8

We use a probability density function of a normal Gaussian distribution for our kernel:

$$ f(x, \mu, \sigma) = p e^q $$

where:

$$ p = \frac{1}{\sqrt{2\pi\sigma^2}} \\ q = -\frac{(x-\mu)^2}{2\sigma^2} $$

See : https://en.wikipedia.org/wiki/Normal_distribution

To make our smoothing kernel, we can sample the PDF, for a given mu and sigma.

We can use our old friend scipy.stats to sample the PDF.

In [ ]:
# - import the normal distribution function from scipy.stats
from scipy.stats import norm
norm_pdf = norm.pdf

Now we need to make our kernel.

  • Work out the +/- limit for the kernel x values with limit = round(sigma * 4)
  • Make a vector of integers x to sample the PDF from -limit to limit - inclusive;
  • Sample the PDF at these x values, for $\mu = 0$ and given $\sigma$, to get the kernel vector;
  • Plot the x values against the y values;
  • Work out how many indices from the start of the kernel vector, to the value corresponding to x=0; (this is the kernel_offset.
In [ ]:
# Work out the +/- limit for the kernel x values;
# Make a vector `x` to sample the PDF;
# Get `kernel` vector by sampling the PDF at these x values (mu=0);
# Work out `kernel_offset`.

Now we can convolve our image line with the kernel. We take the result starting at kernel_offset:

In [ ]:
# Convolve image line with kernel
# Slice out the result we want using kernel_offset
# Plot smoothed with unsmoothed line

Now we are going to smooth a whole slice.

  • Make a new array the same shape as mid_slice called smoothed_slice;
  • Loop over the first dimension of mid_slice, and fill each line of smoothed_slice with the smoothed version;
  • Loop over the second dimension of smoothed_slice and replace each line with the smoothed version;
  • Show.
In [ ]:
# Make a new array `smoothed_slice` for smoothed output
# Loop over first dimension applying kernel as above
# Loop over second dimension of smoothed image, applying kernel
# Show with `imshow`

Now time to do the same thing with scipy.ndimage.gaussian_filter.

In [ ]:
# - import gaussian_filter
from scipy.ndimage import gaussian_filter

Have a look at the help for gaussian_filter.

Use gaussian_filter to make another smoothed version of mid_slice.

In [ ]:
# Use gaussian filter to smooth `mid_slice`

Subtract the gaussian_filter version from your handcrafted version.

Can you explain the difference?

In [ ]:
# Subtract two versions of the smoothed slice, show