Example of running PCA on a 4D image¶
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]:
# - import numpy.linalg with a shorter name
import numpy.linalg as npl
In [5]:
# - set gray colormap and nearest neighbor interpolation by default
plt.rcParams['image.cmap'] = 'gray'
plt.rcParams['image.interpolation'] = 'nearest'
Load the image 'ds114_sub009_t2r1.nii' with nibabel
In [6]:
# Get the data array from the image
import nibabel as nib
img = nib.load('ds107_sub012_t1r2.nii')
data = img.get_data()
Reshape to volumes by voxels
In [7]:
vol_shape = data.shape[:-1]
n_vols = data.shape[-1]
n_voxels = np.prod(vol_shape)
X = data.reshape(n_voxels, n_vols).T
Subtract the mean across columns (row means):
In [8]:
X = X - X.mean(axis=1).reshape(n_vols, 1)
Calculate the unscaled covariance (covariance not divided by N-1):
In [9]:
unscaled_cov = X.dot(X.T)
Get the principal components (array V) using SVD:
In [10]:
U, S, V = npl.svd(unscaled_cov)
Plot the first principal component:
In [11]:
plt.plot(V[0])
Out[11]:
Get the projection coefficients (lengths) for each principal component:
In [12]:
c_values = V.dot(X)
c_values.shape
Out[12]:
Reshape the (number of volumes) by (voxels) array to x by y by z by number of volumes:
In [13]:
c_images = c_values.T.reshape(data.shape)
Show the projection coefficients for the first principal component:
In [14]:
plt.imshow(c_images[:, :, 15, 0])
Out[14]:
Show second principal component:
In [15]:
plt.plot(V[1])
Out[15]:
Show projection coefficients for second principal component:
In [16]:
plt.imshow(c_images[:, :, 15, 1])
Out[16]:
