Diagnosing the FSL data
%pylab inline
Populating the interactive namespace from numpy and matplotlib
We make sure we have some useful libraries loaded:
import numpy as np # array manipulation
import matplotlib.pyplot as plt # plotting library
Set some defaults for plotting:
plt.rcParams['image.cmap'] = 'gray'
plt.rcParams['image.interpolation'] = 'nearest'
We need the library to load and analyze images:
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:
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())
(64, 64, 21, 180)
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
t0_img = rollimg(img, 't') # Put time axis first
time_mean_array = t0_img.get_data().mean(axis=0)
plt.imshow(time_mean_array[:, :, 10])
<matplotlib.image.AxesImage at 0x10060f950>
Try PCA:
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)
img_pca_results.keys()
['basis_vectors over t', 'axis', 'basis_projections', 'basis_vectors', 'pcnt_var']
img_pca_results['basis_vectors'].shape
(180, 179)
fig, axes = plt.subplots(10, 1)
for i, ax in enumerate(axes):
ax.plot(img_pca_results['basis_vectors'][:, i])
pca0 = img_pca_results['basis_projections']
plt.imshow(pca0.get_data()[0, :, :, 5])
<matplotlib.image.AxesImage at 0x107f54510>
Do some more diagostic checks:
dx = screens.screen(img, slice_axis=2)
fig, axes = plt.subplots(4, 1, figsize=(12, 6))
screens.plot_tsdiffs(dx['ts_res'], axes)
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)