1. pna-notebooks
  2. diagnostics.ipynb

Diagnosing the FSL data

In [1]:
%pylab inline

Populating the interactive namespace from numpy and matplotlib

We make sure we have some useful libraries loaded:

In [2]:
import numpy as np # array manipulation
import matplotlib.pyplot as plt # plotting library

Set some defaults for plotting:

In [3]:
plt.rcParams['image.cmap'] = 'gray'
plt.rcParams['image.interpolation'] = 'nearest'

We need the library to load and analyze images:

In [4]:
from nipy import load_image
from nipy.core.api import Image, drop_io_dim, rollimg
from nipy.algorithms.utils.pca import pca_image
from nipy.algorithms.diagnostics import screens

Get the image to look at it:

In [5]:
img = load_image('fmri.nii.gz')
img.shape

/Users/mb312/usr/local/lib/python2.7/site-packages/nipy/io/files.py:53: FutureWarning: Please use ``dataobj`` instead of ``_data``; We will remove this wrapper for ``_data`` soon
  ni_img = nib.Nifti1Image(img._data, img.get_affine(), img.get_header())
Out[5]:
(64, 64, 21, 180)
In [6]:
print(img.metadata['header'])

<class 'nibabel.nifti1.Nifti1Header'> object, endian='>'
sizeof_hdr      : 348
data_type       : 
db_name         : 
extents         : 0
session_error   : 0
regular         : r
dim_info        : 0
dim             : [  4  64  64  21 180   1   1   1]
intent_p1       : 0.0
intent_p2       : 0.0
intent_p3       : 0.0
intent_code     : none
datatype        : int16
bitpix          : 16
slice_start     : 0
pixdim          : [ 1.  4.  4.  6.  1.  1.  1.  1.]
vox_offset      : 352.0
scl_slope       : 0.0
scl_inter       : 0.0
slice_end       : 0
slice_code      : unknown
xyzt_units      : 10
cal_max         : 25500.0
cal_min         : 3.0
slice_duration  : 0.0
toffset         : 0.0
glmax           : 0
glmin           : 0
descrip         : FSL3.2beta
aux_file        : 
qform_code      : unknown
sform_code      : unknown
quatern_b       : 0.0
quatern_c       : 0.0
quatern_d       : 0.0
qoffset_x       : 0.0
qoffset_y       : 0.0
qoffset_z       : 0.0
srow_x          : [ 0.  0.  0.  0.]
srow_y          : [ 0.  0.  0.  0.]
srow_z          : [ 0.  0.  0.  0.]
intent_name     : 
magic           : n+1
In [7]:
t0_img = rollimg(img, 't') # Put time axis first
time_mean_array = t0_img.get_data().mean(axis=0)
In [8]:
plt.imshow(time_mean_array[:, :, 10])
Out[8]:
<matplotlib.image.AxesImage at 0x10060f950>

Try PCA:

In [9]:
img_pca_results = pca_image(t0_img, axis='t', ncomp=10)

/Users/mb312/usr/local/lib/python2.7/site-packages/nipy/algorithms/utils/pca.py:143: RuntimeWarning: divide by zero encountered in divide
  return np.where(rmse<=0, 0, 1. / rmse)
In [10]:
img_pca_results.keys()
Out[10]:
['basis_vectors over t',
 'axis',
 'basis_projections',
 'basis_vectors',
 'pcnt_var']
In [11]:
img_pca_results['basis_vectors'].shape
Out[11]:
(180, 179)
In [12]:
fig, axes = plt.subplots(10, 1)
for i, ax in enumerate(axes):
    ax.plot(img_pca_results['basis_vectors'][:, i])
In [13]:
pca0 = img_pca_results['basis_projections']
plt.imshow(pca0.get_data()[0, :, :, 5])
Out[13]:
<matplotlib.image.AxesImage at 0x107f54510>

Do some more diagostic checks:

In [14]:
dx = screens.screen(img, slice_axis=2)
In [15]:
fig, axes = plt.subplots(4, 1, figsize=(12, 6))
screens.plot_tsdiffs(dx['ts_res'], axes)
Out[15]:
array([<matplotlib.axes.AxesSubplot object at 0x107f5b090>,
       <matplotlib.axes.AxesSubplot object at 0x10abfea10>,
       <matplotlib.axes.AxesSubplot object at 0x10ac24650>,
       <matplotlib.axes.AxesSubplot object at 0x10ac82b50>], dtype=object)