Finding faces, with processed data
%pylab inline
import numpy as np
import matplotlib.pyplot as plt
Populating the interactive namespace from numpy and matplotlib
You'll need the processed data from the website. Use your web browser to download http://nipy.bic.berkeley.edu/practical_neuroimaging/srabold.nii.gz into this current directory, or do:
curl -O http://nipy.bic.berkeley.edu/practical_neuroimaging/srabold.nii.gz
in the terminal / shell (for OSX and Linux)
import nibabel as nib
img = nib.load('srabold.nii.gz')
print(img.header)
<class 'nibabel.nifti1.Nifti1Header'> object, endian='<' sizeof_hdr : 348 data_type : db_name : extents : 0 session_error : 0 regular : dim_info : 0 dim : [ 4 40 64 64 121 1 1 1] intent_p1 : 0.0 intent_p2 : 0.0 intent_p3 : 0.0 intent_code : none datatype : float64 bitpix : 64 slice_start : 0 pixdim : [-1. 3.5 3.75 3.75 1. 1. 1. 1. ] vox_offset : 352.0 scl_slope : 1.0 scl_inter : 0.0 slice_end : 0 slice_code : unknown xyzt_units : 10 cal_max : 0.0 cal_min : 0.0 slice_duration : 0.0 toffset : 0.0 glmax : 0 glmin : 0 descrip : aux_file : qform_code : scanner sform_code : scanner quatern_b : 0.0 quatern_c : 1.0 quatern_d : 0.0 qoffset_x : 68.25 qoffset_y : -118.125 qoffset_z : -118.125 srow_x : [ -3.5 0. 0. 68.25] srow_y : [ 0. 3.75 0. -118.125] srow_z : [ 0. 0. 3.75 -118.125] intent_name : magic : n+1
Defaults for plotting
plt.rcParams['image.cmap'] = 'gray'
plt.rcParams['image.interpolation'] = 'nearest'
Loading the image may take a little while:
data = img.get_data()
data.shape
(40, 64, 64, 121)
plt.imshow(data.mean(axis=-1)[:, :, 30])
<matplotlib.image.AxesImage at 0x100614190>
Why is this image blurry?
So what was the subject doing here?
https://openfmri.org/dataset/ds000105
event_types = ['house', 'scrambled', 'cat', 'shoe', 'bottle', 'scissors', 'chair', 'face']
You should have the task lists for this subject in your ~/data/ds105
directory. Let's check:
import os
HOME = os.path.expanduser('~')
DATA_PATH = os.path.join(HOME, 'data', 'ds105')
SUB1_TASK_PATH = os.path.join(DATA_PATH, 'sub001', 'model', 'model001', 'onsets', 'task001_run001')
os.path.isdir(SUB1_TASK_PATH)
True
If you got "False" above - you can download the zip file of the task files you need from:
http://nipy.bic.berkeley.edu/practical_neuroimaging/ds105_sub001_tasks.zip
Download the file, then cd ~/data
and unzip ~/Downloads/ds105_sub001_tasks.zip
(if you downloaded the file to your Downloads
folder.
Recheck:
os.path.isdir(SUB1_TASK_PATH)
True
Look in this path:
os.listdir(SUB1_TASK_PATH)
['cond001.txt', 'cond002.txt', 'cond003.txt', 'cond004.txt', 'cond005.txt', 'cond006.txt', 'cond007.txt', 'cond008.txt']
What are in these files?
COND1_FNAME = os.path.join(SUB1_TASK_PATH, 'cond001.txt')
print(open(COND1_FNAME, 'rt').read())
156.000 0.5 1 158.000 0.5 1 160.000 0.5 1 162.000 0.5 1 164.000 0.5 1 166.000 0.5 1 168.000 0.5 1 170.000 0.5 1 172.000 0.5 1 174.000 0.5 1 176.000 0.5 1 178.000 0.5 1
This is a common format for specifying events - any guesses what the columns are?
I want to read this file to give me something I can use to specify my design. Something like this:
def evdefs_from_fname(fname):
""" Get onset, duration, amplitude for event from file
Parameters
----------
fname : str
Filename to read data from
Returns
-------
evdefs : list of lists
list, one element per event, where the elements are lists containing
(onset, duration, amplitude)
"""
evdefs = []
with open(fname, 'rt') as fobj:
for line in fobj:
entries = [float(e) for e in line.split()]
evdefs.append(entries)
return evdefs
# Of course we need to test
evdefs = evdefs_from_fname(COND1_FNAME)
assert len(evdefs) == 12
assert evdefs[0] == [156.000, 0.5, 1]
assert evdefs[-1] == [178.000, 0.5, 1]
We can actually do something neater than that:
dtype_def = [('onset', 'f'), ('duration', 'f'), ('amplitude', 'f')]
evs = np.loadtxt(COND1_FNAME, dtype_def)
evs
array([(156.0, 0.5, 1.0), (158.0, 0.5, 1.0), (160.0, 0.5, 1.0), (162.0, 0.5, 1.0), (164.0, 0.5, 1.0), (166.0, 0.5, 1.0), (168.0, 0.5, 1.0), (170.0, 0.5, 1.0), (172.0, 0.5, 1.0), (174.0, 0.5, 1.0), (176.0, 0.5, 1.0), (178.0, 0.5, 1.0)], dtype=[('onset', '<f4'), ('duration', '<f4'), ('amplitude', '<f4')])
This is nice because:
evs[0]['onset']
156.0
evs['onset']
array([ 156., 158., 160., 162., 164., 166., 168., 170., 172., 174., 176., 178.], dtype=float32)
We can load our events for all types:
event_defs = {}
for eno in range(8):
fname = os.path.join(os.path.join(SUB1_TASK_PATH, 'cond%03d.txt' % (eno + 1)))
print(fname)
event_name = event_types[eno]
event_defs[event_name] = np.loadtxt(fname, dtype_def)
event_defs
/Users/mb312/data/ds105/sub001/model/model001/onsets/task001_run001/cond001.txt /Users/mb312/data/ds105/sub001/model/model001/onsets/task001_run001/cond002.txt /Users/mb312/data/ds105/sub001/model/model001/onsets/task001_run001/cond003.txt /Users/mb312/data/ds105/sub001/model/model001/onsets/task001_run001/cond004.txt /Users/mb312/data/ds105/sub001/model/model001/onsets/task001_run001/cond005.txt /Users/mb312/data/ds105/sub001/model/model001/onsets/task001_run001/cond006.txt /Users/mb312/data/ds105/sub001/model/model001/onsets/task001_run001/cond007.txt /Users/mb312/data/ds105/sub001/model/model001/onsets/task001_run001/cond008.txt
{'bottle': array([(228.0, 0.5, 1.0), (230.0, 0.5, 1.0), (232.0, 0.5, 1.0), (234.0, 0.5, 1.0), (236.0, 0.5, 1.0), (238.0, 0.5, 1.0), (240.0, 0.5, 1.0), (242.0, 0.5, 1.0), (244.0, 0.5, 1.0), (246.0, 0.5, 1.0), (248.0, 0.5, 1.0), (250.0, 0.5, 1.0)], dtype=[('onset', '<f4'), ('duration', '<f4'), ('amplitude', '<f4')]), 'cat': array([(84.0, 0.5, 1.0), (86.0, 0.5, 1.0), (88.0, 0.5, 1.0), (90.0, 0.5, 1.0), (92.0, 0.5, 1.0), (94.0, 0.5, 1.0), (96.0, 0.5, 1.0), (98.0, 0.5, 1.0), (100.0, 0.5, 1.0), (102.0, 0.5, 1.0), (104.0, 0.5, 1.0), (106.0, 0.5, 1.0)], dtype=[('onset', '<f4'), ('duration', '<f4'), ('amplitude', '<f4')]), 'chair': array([(264.0, 0.5, 1.0), (266.0, 0.5, 1.0), (268.0, 0.5, 1.0), (270.0, 0.5, 1.0), (272.0, 0.5, 1.0), (274.0, 0.5, 1.0), (276.0, 0.5, 1.0), (278.0, 0.5, 1.0), (280.0, 0.5, 1.0), (282.0, 0.5, 1.0), (284.0, 0.5, 1.0), (286.0, 0.5, 1.0)], dtype=[('onset', '<f4'), ('duration', '<f4'), ('amplitude', '<f4')]), 'face': array([(48.0, 0.5, 1.0), (50.0, 0.5, 1.0), (52.0, 0.5, 1.0), (54.0, 0.5, 1.0), (56.0, 0.5, 1.0), (58.0, 0.5, 1.0), (60.0, 0.5, 1.0), (62.0, 0.5, 1.0), (64.0, 0.5, 1.0), (66.0, 0.5, 1.0), (68.0, 0.5, 1.0), (70.0, 0.5, 1.0)], dtype=[('onset', '<f4'), ('duration', '<f4'), ('amplitude', '<f4')]), 'house': array([(156.0, 0.5, 1.0), (158.0, 0.5, 1.0), (160.0, 0.5, 1.0), (162.0, 0.5, 1.0), (164.0, 0.5, 1.0), (166.0, 0.5, 1.0), (168.0, 0.5, 1.0), (170.0, 0.5, 1.0), (172.0, 0.5, 1.0), (174.0, 0.5, 1.0), (176.0, 0.5, 1.0), (178.0, 0.5, 1.0)], dtype=[('onset', '<f4'), ('duration', '<f4'), ('amplitude', '<f4')]), 'scissors': array([(12.0, 0.5, 1.0), (14.0, 0.5, 1.0), (16.0, 0.5, 1.0), (18.0, 0.5, 1.0), (20.0, 0.5, 1.0), (22.0, 0.5, 1.0), (24.0, 0.5, 1.0), (26.0, 0.5, 1.0), (28.0, 0.5, 1.0), (30.0, 0.5, 1.0), (32.0, 0.5, 1.0), (34.0, 0.5, 1.0)], dtype=[('onset', '<f4'), ('duration', '<f4'), ('amplitude', '<f4')]), 'scrambled': array([(192.0, 0.5, 1.0), (194.0, 0.5, 1.0), (196.0, 0.5, 1.0), (198.0, 0.5, 1.0), (200.0, 0.5, 1.0), (202.0, 0.5, 1.0), (204.0, 0.5, 1.0), (206.0, 0.5, 1.0), (208.0, 0.5, 1.0), (210.0, 0.5, 1.0), (212.0, 0.5, 1.0), (214.0, 0.5, 1.0)], dtype=[('onset', '<f4'), ('duration', '<f4'), ('amplitude', '<f4')]), 'shoe': array([(120.0, 0.5, 1.0), (122.0, 0.5, 1.0), (124.0, 0.5, 1.0), (126.0, 0.5, 1.0), (128.0, 0.5, 1.0), (130.0, 0.5, 1.0), (132.0, 0.5, 1.0), (134.0, 0.5, 1.0), (136.0, 0.5, 1.0), (138.0, 0.5, 1.0), (140.0, 0.5, 1.0), (142.0, 0.5, 1.0)], dtype=[('onset', '<f4'), ('duration', '<f4'), ('amplitude', '<f4')])}
Are these onsets in TRs or in seconds?
OK - how to make our design?
We could do this by hand (right?) - but let's let nipy
do this for us:
import nipy.modalities.fmri.design as fmrid
Try fmrid.block_design?
at the prompt.
So - we need start and stop for our events:
nipy_dtype = [('start', 'f'), ('end', 'f'), ('type', 'S10')]
n_events = sum([len(evs) for e in event_defs])
nipy_block_specs = np.zeros((n_events,), dtype=nipy_dtype)
event_no = 0
for ev_type in event_defs:
ev_def = event_defs[ev_type]
for e in ev_def:
nipy_block_specs[event_no]['start'] = e['onset']
nipy_block_specs[event_no]['end'] = (e['onset'] + e['duration'])
nipy_block_specs[event_no]['type'] = ev_type
event_no += 1
nipy_block_defs = np.sort(nipy_block_specs)
nipy_block_specs
array([(156.0, 156.5, 'house'), (158.0, 158.5, 'house'), (160.0, 160.5, 'house'), (162.0, 162.5, 'house'), (164.0, 164.5, 'house'), (166.0, 166.5, 'house'), (168.0, 168.5, 'house'), (170.0, 170.5, 'house'), (172.0, 172.5, 'house'), (174.0, 174.5, 'house'), (176.0, 176.5, 'house'), (178.0, 178.5, 'house'), (84.0, 84.5, 'cat'), (86.0, 86.5, 'cat'), (88.0, 88.5, 'cat'), (90.0, 90.5, 'cat'), (92.0, 92.5, 'cat'), (94.0, 94.5, 'cat'), (96.0, 96.5, 'cat'), (98.0, 98.5, 'cat'), (100.0, 100.5, 'cat'), (102.0, 102.5, 'cat'), (104.0, 104.5, 'cat'), (106.0, 106.5, 'cat'), (192.0, 192.5, 'scrambled'), (194.0, 194.5, 'scrambled'), (196.0, 196.5, 'scrambled'), (198.0, 198.5, 'scrambled'), (200.0, 200.5, 'scrambled'), (202.0, 202.5, 'scrambled'), (204.0, 204.5, 'scrambled'), (206.0, 206.5, 'scrambled'), (208.0, 208.5, 'scrambled'), (210.0, 210.5, 'scrambled'), (212.0, 212.5, 'scrambled'), (214.0, 214.5, 'scrambled'), (120.0, 120.5, 'shoe'), (122.0, 122.5, 'shoe'), (124.0, 124.5, 'shoe'), (126.0, 126.5, 'shoe'), (128.0, 128.5, 'shoe'), (130.0, 130.5, 'shoe'), (132.0, 132.5, 'shoe'), (134.0, 134.5, 'shoe'), (136.0, 136.5, 'shoe'), (138.0, 138.5, 'shoe'), (140.0, 140.5, 'shoe'), (142.0, 142.5, 'shoe'), (228.0, 228.5, 'bottle'), (230.0, 230.5, 'bottle'), (232.0, 232.5, 'bottle'), (234.0, 234.5, 'bottle'), (236.0, 236.5, 'bottle'), (238.0, 238.5, 'bottle'), (240.0, 240.5, 'bottle'), (242.0, 242.5, 'bottle'), (244.0, 244.5, 'bottle'), (246.0, 246.5, 'bottle'), (248.0, 248.5, 'bottle'), (250.0, 250.5, 'bottle'), (12.0, 12.5, 'scissors'), (14.0, 14.5, 'scissors'), (16.0, 16.5, 'scissors'), (18.0, 18.5, 'scissors'), (20.0, 20.5, 'scissors'), (22.0, 22.5, 'scissors'), (24.0, 24.5, 'scissors'), (26.0, 26.5, 'scissors'), (28.0, 28.5, 'scissors'), (30.0, 30.5, 'scissors'), (32.0, 32.5, 'scissors'), (34.0, 34.5, 'scissors'), (264.0, 264.5, 'chair'), (266.0, 266.5, 'chair'), (268.0, 268.5, 'chair'), (270.0, 270.5, 'chair'), (272.0, 272.5, 'chair'), (274.0, 274.5, 'chair'), (276.0, 276.5, 'chair'), (278.0, 278.5, 'chair'), (280.0, 280.5, 'chair'), (282.0, 282.5, 'chair'), (284.0, 284.5, 'chair'), (286.0, 286.5, 'chair'), (48.0, 48.5, 'face'), (50.0, 50.5, 'face'), (52.0, 52.5, 'face'), (54.0, 54.5, 'face'), (56.0, 56.5, 'face'), (58.0, 58.5, 'face'), (60.0, 60.5, 'face'), (62.0, 62.5, 'face'), (64.0, 64.5, 'face'), (66.0, 66.5, 'face'), (68.0, 68.5, 'face'), (70.0, 70.5, 'face')], dtype=[('start', '<f4'), ('end', '<f4'), ('type', 'S10')])
Nice - so let's make the design:
n_trs = img.shape[-1]
TR = 2.5
t = np.arange(n_trs) * TR
X_exper, cons_exper = fmrid.block_design(nipy_block_defs, t, level_contrasts=True)
plt.figure(figsize=(10, 20))
plt.imshow(X_exper)
<matplotlib.image.AxesImage at 0x143c14690>
np.set_printoptions(suppress=True)
cons_exper
{'constant_0': array([ 1., 1., 1., 1., 1., 1., 1., 1.]), 'type_bottle_0': array([ 1., -0., -0., -0., -0., -0., -0., -0.]), 'type_cat_0': array([ 0., 1., -0., 0., -0., -0., -0., -0.]), 'type_chair_0': array([-0., -0., 1., -0., -0., 0., -0., 0.]), 'type_face_0': array([-0., 0., -0., 1., -0., -0., -0., -0.]), 'type_house_0': array([-0., -0., 0., -0., 1., 0., -0., -0.]), 'type_scissors_0': array([-0., -0., -0., -0., 0., 1., -0., 0.]), 'type_scrambled_0': array([ 0., -0., -0., 0., -0., -0., 1., 0.]), 'type_shoe_0': array([-0., -0., 0., 0., -0., -0., -0., 1.])}
Next we make the drift terms. Here we use a spline set:
drift = fmrid.natural_spline(t, [t.mean()])
# Make confounds a little easier to display
drift = drift / drift.max(axis=0)
plt.figure(figsize=(10, 20))
plt.imshow(drift)
<matplotlib.image.AxesImage at 0x143b14410>
We combine that with our event design matrix:
X, cons = fmrid.stack_designs((X_exper, cons_exper),
(drift, {}))
plt.figure(figsize=(20, 10))
plt.imshow(X)
<matplotlib.image.AxesImage at 0x143b83410>
Fit the design, the by-hand way:
# Reshape the data for estimation
slice_size = np.prod(data.shape[:-1])
Y = np.reshape(data, (slice_size, n_trs)).T
Y.shape
(121, 163840)
B = np.linalg.pinv(X).dot(Y)
B.shape
(13, 163840)
Let's make a contrast map. But, what are our contrasts?
cons
{'constant_0': array([ 1., 1., 1., 1., 1., 1., 1., 1., 0., 0., 0., 0., 0.]), 'type_bottle_0': array([ 1., -0., -0., -0., -0., -0., -0., -0., 0., 0., 0., 0., 0.]), 'type_cat_0': array([ 0., 1., -0., 0., -0., -0., -0., -0., 0., 0., 0., 0., 0.]), 'type_chair_0': array([-0., -0., 1., -0., -0., 0., -0., 0., 0., 0., 0., 0., 0.]), 'type_face_0': array([-0., 0., -0., 1., -0., -0., -0., -0., 0., 0., 0., 0., 0.]), 'type_house_0': array([-0., -0., 0., -0., 1., 0., -0., -0., 0., 0., 0., 0., 0.]), 'type_scissors_0': array([-0., -0., -0., -0., 0., 1., -0., 0., 0., 0., 0., 0., 0.]), 'type_scrambled_0': array([ 0., -0., -0., 0., -0., -0., 1., 0., 0., 0., 0., 0., 0.]), 'type_shoe_0': array([-0., -0., 0., 0., -0., -0., -0., 1., 0., 0., 0., 0., 0.])}
contrast = cons['type_face_0'] - cons['type_house_0']
contrast
array([ 0., 0., -0., 1., -1., -0., 0., 0., 0., 0., 0., 0., 0.])
con_data = contrast.dot(B)
con_data.shape
(163840,)
con_img = np.reshape(con_data, img.shape[:-1])
con_img.shape
(40, 64, 64)
plt.imshow(con_img[..., 30])
<matplotlib.image.AxesImage at 0x143bb2c10>
Can you find any face activity you believe?