Skip to content

Commit a8ca43b

Browse files
pzeidmanpzeidman
committed
Add developer example
1 parent a954b7e commit a8ca43b

File tree

2 files changed

+200
-0
lines changed

2 files changed

+200
-0
lines changed

‎example_3T/scripts/Run_developer_demo_analysis.m

+66
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,66 @@
1+
% Example script for demonstrating the example pRF function for developers
2+
3+
% Settings
4+
5+
% Directory of the downloaded example dataset
6+
data_root_dir = 'D:\samdata\';
7+
data_dir = fullfile(data_root_dir,'example','pRF');
8+
9+
% Directory of GLM
10+
glm_dir = fullfile(pwd,'../GLM');
11+
12+
TR = 1; % Repetition time
13+
TE = 0.055; % Echo time
14+
15+
% Which sessions to include
16+
sess = 1:10;
17+
num_sess = length(sess);
18+
19+
%% Prepare inputs
20+
21+
% Build a structure containing which stimulus pixels were illuminated at
22+
% each time step.
23+
load(fullfile(data_dir,'aps_Bars.mat'));
24+
U = prepare_inputs_polar_samsrf(ApFrm,TR);
25+
26+
% The timeseries from each session are stored in a VOI_xx.mat file. Build a
27+
% cell array of the VOI files for each session.
28+
xY = cell(1,num_sess);
29+
for i = 1:num_sess
30+
filename = sprintf('VOI_Mask_%d.mat',sess(i));
31+
xY{i} = fullfile(glm_dir,filename);
32+
end
33+
%% Specify pRF model (all 6669 voxels)
34+
35+
% Load SPM for timing information / image dimensions
36+
SPM = load(fullfile(glm_dir,'SPM.mat'));
37+
SPM = SPM.SPM;
38+
39+
% Update SPM path as we don't know where this example will be saved
40+
SPM.swd = glm_dir;
41+
42+
% Set pRF specification options
43+
options = struct('TE', TE,...
44+
'voxel_wise', true,...
45+
'name', 'developer_example',...
46+
'model', 'spm_prf_fcn_template',...
47+
'B0',3);
48+
49+
% Specify pRF model (.mat file will be stored in the GLM directory)
50+
PRF = spm_prf_analyse('specify',SPM,xY,U,options);
51+
52+
%% Estimate one voxel (global peak) as an example: idx 3439, [-4,-70,-3]
53+
54+
% Model to estimate
55+
prf_file = fullfile(glm_dir,'PRF_developer_example.mat');
56+
57+
% Estimation options
58+
options = struct('voxels',3439,'init','none');
59+
60+
% Estimate
61+
PRF_est = spm_prf_analyse('estimate',prf_file,options);
62+
63+
%% Review (single voxel)
64+
prf_file = fullfile(glm_dir,'PRF_developer_example.mat');
65+
66+
spm_prf_review(prf_file, 3439);

‎toolbox/response_functions/spm_prf_fcn_template.m

+134
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,134 @@
1+
function varargout = spm_prf_fcn_template(P,M,U,varargin)
2+
% Template pRF response function. This example models neuronal
3+
% activity as a scaled version of the inputs:
4+
%
5+
% z(t) = alpha * u(t)
6+
%
7+
% Where alpha is a parameter estimated from the data.
8+
%
9+
%
10+
% Inputs:
11+
%
12+
% P - parameter structure
13+
% M - model structure
14+
% U - experimental timing
15+
% action - (optional) the action to performm. If not provided, the
16+
% predicted BOLD and neuronal timeseries are returned.
17+
%
18+
% -------------------------------------------------------------------------
19+
% FORMAT [y,Z] = spm_prf_fcn_template(P,M,U)
20+
% Return the BOLD and neuronal predicted timeseries
21+
%
22+
% P parameters
23+
% M,U model, inputs
24+
%
25+
% y fMRI time series
26+
% Z neuronal response
27+
% -------------------------------------------------------------------------
28+
% FORMAT P = spm_prf_fcn_template(P,M,U,'get_parameters')
29+
% Return the given parameters corrected for display
30+
%
31+
% P parameters
32+
% M,U model, inputs
33+
% -------------------------------------------------------------------------
34+
% FORMAT tf = spm_prf_fcn_template(P,M,U,'is_above_threshold',Cp,v)
35+
% Return whether the model with parameters P and covariance Cp passes an
36+
% arbitrary threshold for display
37+
%
38+
% P parameters
39+
% M,U model, inputs
40+
% Cp parameter covariance matrix
41+
% v voxel index
42+
% -------------------------------------------------------------------------
43+
% FORMAT x = spm_prf_fcn_template(P,M,U,'get_response',xy)
44+
% Return the instantaneous response of the PRF at coordinates xy
45+
%
46+
% P parameters
47+
% M,U model, inputs
48+
% xy [2xn] vector of coordinates to evaluate the PRF
49+
% -------------------------------------------------------------------------
50+
% FORMAT [pE,pC] = spm_prf_fcn_template(P,M,U,'get_priors')
51+
% Return the priors for the model. Importantly, this defines which
52+
% parameters are in the model.
53+
%
54+
% pE structure or vector of prior expectations
55+
% pC prior covariance maitrx
56+
%
57+
% ---------------------------------------------------------------------
58+
% Copyright (C) 2016 Peter Zeidman
59+
% This program is free software: you can redistribute it and/or modify
60+
% it under the terms of the GNU General Public License as published by
61+
% the Free Software Foundation, either version 3 of the License, or
62+
% (at your option) any later version.
63+
%
64+
% This program is distributed in the hope that it will be useful,
65+
% but WITHOUT ANY WARRANTY; without even the implied warranty of
66+
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
67+
% GNU General Public License for more details.
68+
%
69+
% You should have received a copy of the GNU General Public License
70+
% along with this program. If not, see <http://www.gnu.org/licenses/>.
71+
% ---------------------------------------------------------------------
72+
73+
if nargin < 4
74+
% Integrate the model over time and return neuronal timeseries z and
75+
% BOLD timeseries y. As an example, here we have neuronal model
76+
% z(t) = alpha, where alpha is an estimated parameter.
77+
78+
% Number of volumes = inputs
79+
n = length(U);
80+
81+
% Neural timeseries. A vector with one entry per microtime bin. The
82+
% U.nbins field is injected automatially by spm_prf_analyse
83+
z = zeros(1,U(1).nbins);
84+
85+
for t = 1:n
86+
% Microtime index for this volume
87+
ind = U(t).ind;
88+
89+
% pRF response
90+
z(ind) = P.alpha;
91+
end
92+
93+
% Integrate the BOLD model
94+
Z.u=z';
95+
Z.dt=M.dt;
96+
y=spm_int(P,M,Z);
97+
98+
varargout{1} = y;
99+
varargout{2} = Z;
100+
else
101+
% This section of the code provides information on the model, primarily
102+
% for plotting purposes in spm_prf_review()
103+
104+
action = varargin{1};
105+
106+
switch action
107+
case 'get_parameters'
108+
% Get the parameters with any corrections needed for
109+
% display
110+
varargout{1} = P;
111+
case 'is_above_threshold'
112+
% Return binary vector identifying whether each voxel is
113+
% above some threshold for display
114+
varargout{1} = 1;
115+
case 'get_response'
116+
% Return the prediction of the model at coordinates xy
117+
xy = varargin{2};
118+
varargout{1} = ones(1,length(xy));
119+
case 'get_priors'
120+
% Return a structure containing the priors
121+
pE.alpha = 1;
122+
pC.alpha = 1;
123+
varargout{1} = pE;
124+
varargout{2} = pC;
125+
case 'glm_initialize'
126+
% (Optional) Return parameters initialized using some
127+
% rapid initial search on timeseries y
128+
y = varargin{2};
129+
130+
varargout = P;
131+
otherwise
132+
error('Unknown action');
133+
end
134+
end

0 commit comments

Comments
 (0)