In [1]:
%pylab inline
In [2]:
import numpy as np
import matplotlib.pyplot as plt
In [3]:
import os
import nibabel as nib
The MNI templates are here: http://www.bic.mni.mcgill.ca/~vfonov/icbm/2009/mni_icbm152_nlin_asym_09a_nifti.zip
I unzipped these into the directory of this notebook, giving me the subdirectory below.
In [4]:
MNI_PATH = 'mni_icbm152_nlin_asym_09a'
The templates contain probability maps of gray matter, white matter and CSF. The images have values between 0 and 1 giving the probability that the voxel is (for example) within gray mattr.
In [5]:
gm_fname = os.path.join(MNI_PATH, 'mni_icbm152_gm_tal_nlin_asym_09a.nii')
wm_fname = os.path.join(MNI_PATH, 'mni_icbm152_wm_tal_nlin_asym_09a.nii')
csf_fname = os.path.join(MNI_PATH, 'mni_icbm152_csf_tal_nlin_asym_09a.nii')
In [6]:
gm_img = nib.load(gm_fname)
gm_data = gm_img.get_data()
wm_img = nib.load(wm_fname)
wm_data = wm_img.get_data()
csf_img = nib.load(csf_fname)
csf_data = csf_img.get_data()
csf_data.shape
Out[6]:
What do the real T1 and T2 images look like?
In [7]:
t1_fname = os.path.join(MNI_PATH, 'mni_icbm152_t1_tal_nlin_asym_09a.nii')
t2_fname = os.path.join(MNI_PATH, 'mni_icbm152_t2_tal_nlin_asym_09a.nii')
In [8]:
t1_slice = nib.load(t1_fname).get_data()[:, :, 94]
t2_slice = nib.load(t2_fname).get_data()[:, :, 94]
plt.imshow(np.hstack((t1_slice, t2_slice)), cmap="gray")
Out[8]:
We make a pretend, simpler T1-like and T2-like slice from the gray matter, white matter, CSF
In [9]:
# Get binary slices for gray matter, white matter, CSF
mid_gm = gm_data[:, :, 94] > 0.5
mid_wm = wm_data[:, :, 94] > 0.5
mid_csf = csf_data[:, :, 94] > 0.5
In [10]:
plt.imshow(np.hstack((mid_gm, mid_wm, mid_csf)), cmap="gray")
Out[10]:
In [11]:
# Combine GM, WM, CSF in different proportions for T1-like and T2-like
t1_like = mid_gm * 0.7 + mid_wm * 1.0 + mid_csf * 0.1
t2_like = mid_gm * 0.5 + mid_wm * 0.2 + mid_csf * 0.8
plt.imshow(np.hstack((t1_like, t2_like)), cmap="gray")
Out[11]:
In [12]:
plt.plot(t1_like.ravel(), t2_like.ravel(), '.')
plt.axis([-0.1, 1.1, -0.1, 1.1])
Out[12]:
In [13]:
hgram_registered, x_edges, y_edges = np.histogram2d(t1_like.ravel(), t2_like.ravel(), 30)
hgram_registered.max()
Out[13]:
In [14]:
plt.imshow(np.clip(hgram_registered, 0, 1000), cmap='gray', interpolation='nearest')
Out[14]:
Remember the histogram image is flipped up-down compared to the plot, because 0, 0 by convention starts at the top left for the image display.
Mutual information is defined like this:
$$ I(X;Y) = \sum_{y \in Y} \sum_{x \in X} p(x,y) \log{ \left(\frac{p(x,y)}{p(x)\,p(y)} \right) } $$In [15]:
def mutual_information(hgram):
""" Mutual information for joint histogram
"""
pxy = hgram / float(np.sum(hgram)) # Convert to probability
px = np.sum(pxy, 1) # marginal for x over y
py = np.sum(pxy, 0) # marginal for y over x
px_py = px[:, None] * py[None, :] # Broadcast to multiply marginals
nzs = pxy > 0 # Only non-zero pxy values contribute to the sum
return np.sum(pxy[nzs] * np.log(pxy[nzs] / px_py[nzs]))
In [16]:
mutual_information(hgram_registered)
Out[16]:
Move T2-like image 30 pixels down
In [17]:
t2_like_moved = np.zeros_like(t2_like)
t2_like_moved[10:, :] = t2_like[:-10, :]
plt.imshow(t2_like_moved, cmap="gray")
Out[17]:
In [18]:
plt.plot(t1_like.ravel(), t2_like_moved.ravel(), '.')
plt.axis([-0.1, 1.1, -0.1, 1.1])
Out[18]:
In [19]:
hgram_off, x_edges, y_edges = np.histogram2d(t1_like.ravel(), t2_like_moved.ravel(), 30)
plt.imshow(np.clip(hgram_off, 0, 1000), cmap='gray', interpolation='nearest')
Out[19]:
In [20]:
mutual_information(hgram_off)
Out[20]: