%pylab inline
Welcome to pylab, a matplotlib-based Python environment [backend: module://IPython.kernel.zmq.pylab.backend_inline]. For more information, type 'help(pylab)'.
import numpy as np
import matplotlib.pyplot as plt
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.
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.
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')
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
(197, 233, 189)
What do the real T1 and T2 images look like?
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')
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")
<matplotlib.image.AxesImage at 0x421c250>
We make a pretend, simpler T1-like and T2-like slice from the gray matter, white matter, CSF
# 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
plt.imshow(np.hstack((mid_gm, mid_wm, mid_csf)), cmap="gray")
<matplotlib.image.AxesImage at 0x4535310>
# 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")
<matplotlib.image.AxesImage at 0x4837550>
plt.plot(t1_like.ravel(), t2_like.ravel(), '.')
plt.axis([-0.1, 1.1, -0.1, 1.1])
[-0.1, 1.1, -0.1, 1.1]
hgram_registered, x_edges, y_edges = np.histogram2d(t1_like.ravel(), t2_like.ravel(), 30)
hgram_registered.max()
27157.0
plt.imshow(np.clip(hgram_registered, 0, 1000), cmap='gray', interpolation='nearest')
<matplotlib.image.AxesImage at 0x4860a10>
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) } $$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]))
mutual_information(hgram_registered)
1.0374382238573592
Move T2-like image 30 pixels down
t2_like_moved = np.zeros_like(t2_like)
t2_like_moved[10:, :] = t2_like[:-10, :]
plt.imshow(t2_like_moved, cmap="gray")
<matplotlib.image.AxesImage at 0x4c97690>
plt.plot(t1_like.ravel(), t2_like_moved.ravel(), '.')
plt.axis([-0.1, 1.1, -0.1, 1.1])
[-0.1, 1.1, -0.1, 1.1]
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')
<matplotlib.image.AxesImage at 0x51805d0>
mutual_information(hgram_off)
0.41711040301909241