1. pna2015
  2. day12
In [1]:
# - 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 [2]:
# - show figures inside the notebook
%matplotlib inline
In [3]:
# - import common modules
import numpy as np  # the Python array package
import matplotlib.pyplot as plt  # the Python plotting package
In [4]:
# - 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 [5]:
# - 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 [6]:
# - get slice over third dimension
mid_slice = vol0[:, :, 14]
plt.imshow(mid_slice)
Out[6]:
<matplotlib.image.AxesImage at 0x109727610>

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

In [7]:
# - get middle line from middle slice
mid_line = mid_slice[:, 25]
plt.plot(mid_line)
Out[7]:
[<matplotlib.lines.Line2D at 0x109870290>]

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

Copy the fwhm2sigma function from the smoothing as convolution notebook:

In [8]:
# fwhm2sigma
def fwhm2sigma(fwhm):
    return fwhm / np.sqrt(8 * np.log(2))
In [9]:
# - 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 [10]:
# sigma for FWHM 8
sigma = fwhm2sigma(8)
sigma
Out[10]:
3.3972872011520763

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 [11]:
# - 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 [12]:
# 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`.

limit = round(sigma * 4)
# Make an x range between -sigma * 3 and +sigma * 3
x_for_kernel = np.arange(-limit, limit+1)
# Calculate kernel
kernel = norm_pdf(x_for_kernel, 0, sigma)
# Plot kernel
plt.plot(x_for_kernel, kernel)
# Calculate `kernel_offset`
kernel_offset = int(limit)
kernel_offset
Out[12]:
14

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

In [13]:
# Convolve image line with kernel
# Slice out the result we want using kernel_offset
# Plot smoothed with unsmoothed line
smoothed = np.convolve(mid_line, kernel)
smoothed = smoothed[kernel_offset:len(mid_line)+kernel_offset]
plt.plot(mid_line)
plt.plot(smoothed)
Out[13]:
[<matplotlib.lines.Line2D at 0x10a0823d0>]

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 [14]:
# 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`
smoothed_slice = np.zeros(mid_slice.shape)
n_x, n_y = mid_slice.shape
for x in range(n_x):
    line = mid_slice[x, :]
    smoothed = np.convolve(line, kernel)
    smoothed_slice[x, :] = smoothed[kernel_offset:n_y+kernel_offset]
for y in range(n_y):
    line = smoothed_slice[:, y]
    smoothed = np.convolve(line, kernel)
    smoothed_slice[:, y] = smoothed[kernel_offset:n_x+kernel_offset]
plt.imshow(smoothed_slice)
Out[14]:
<matplotlib.image.AxesImage at 0x10a0dcc50>

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

In [15]:
# - 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 [16]:
# Use gaussian filter to smooth `mid_slice`
scipy_slice = gaussian_filter(mid_slice, sigma)
plt.imshow(scipy_slice)
Out[16]:
<matplotlib.image.AxesImage at 0x10a320d50>

Subtract the gaussian_filter version from your handcrafted version.

Can you explain the difference?

In [17]:
# Subtract two versions of the smoothed slice, show
diff_slice = smoothed_slice - scipy_slice
plt.imshow(diff_slice)
Out[17]:
<matplotlib.image.AxesImage at 0x10a4706d0>