1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
|
#!/usr/bin/env python3
# emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*-
# vi: set ft=python sts=4 ts=4 sw=4 et:
__doc__ = """
This is an example where:
1. An sequence of fMRI volumes are simulated
2. A design matrix describing all the effects related to the data is computed
3. A GLM is applied to all voxels
4. A contrast image is created
Requires matplotlib
Author : Bertrand Thirion, 2010
"""
print(__doc__)
import os
import os.path as op
import numpy as np
try:
import matplotlib.pyplot as plt
except ImportError:
raise RuntimeError("This script needs the matplotlib library")
from nibabel import Nifti1Image, save
import nipy.modalities.fmri.design_matrix as dm
from nipy.labs.utils.simul_multisubject_fmri_dataset import surrogate_4d_dataset
from nipy.modalities.fmri.experimental_paradigm import EventRelatedParadigm
from nipy.modalities.fmri.glm import GeneralLinearModel
#######################################
# Simulation parameters
#######################################
# volume mask
shape = (20, 20, 20)
affine = np.eye(4)
# Acquisition parameters: number of scans (n_scans) and volume repetition time
# value in seconds
n_scans = 128
tr = 2.4
# input paradigm information
frametimes = np.linspace(0, (n_scans - 1) * tr, n_scans)
# conditions are 0 1 0 1 0 1 ...
conditions = np.arange(20) % 2
# 20 onsets (in sec), first event 10 sec after the start of the first scan
onsets = np.linspace(5, (n_scans - 1) * tr - 10, 20)
# model with canonical HRF (could also be :
# 'canonical with derivative' or 'fir'
hrf_model = 'canonical'
# fake motion parameters to be included in the model
motion = np.cumsum(np.random.randn(n_scans, 6), 0)
add_reg_names = ['tx', 'ty', 'tz', 'rx', 'ry', 'rz']
########################################
# Design matrix
########################################
paradigm = EventRelatedParadigm(conditions, onsets)
X, names = dm.dmtx_light(frametimes, paradigm, drift_model='cosine',
hfcut=128, hrf_model=hrf_model, add_regs=motion,
add_reg_names=add_reg_names)
#######################################
# Get the FMRI data
#######################################
fmri_data = surrogate_4d_dataset(shape=shape, n_scans=n_scans)[0]
# if you want to save it as an image
data_file = 'fmri_data.nii'
save(fmri_data, data_file)
########################################
# Perform a GLM analysis
########################################
# GLM fit
Y = fmri_data.get_fdata().reshape(np.prod(shape), n_scans)
glm = GeneralLinearModel(X)
glm.fit(Y.T)
# specify the contrast [1 -1 0 ..]
contrast = np.zeros(X.shape[1])
contrast[0] = 1
contrast[1] = - 1
# compute the contrast image related to it
zvals = glm.contrast(contrast).z_score()
contrast_image = Nifti1Image(np.reshape(zvals, shape), affine)
# if you want to save the contrast as an image
contrast_path = 'zmap.nii'
save(contrast_image, contrast_path)
print(f'Wrote the some of the results as images in directory {op.abspath(os.getcwd())}')
h, c = np.histogram(zvals, 100)
# Show the histogram
plt.figure()
plt.bar(c[: - 1], h, width=.1)
plt.title(' Histogram of the z-values')
plt.show()
|
