source: misclib/fmri_analysis.m @ 188

Last change on this file since 188 was 188, checked in by tjniskan, 7 years ago
  • Fixed a bug in fmri analysis. Threshold estimation failed when fmri time series had no variation.

M misclib/fmri_analysis.m
M aedes_revision.m

File size: 8.9 KB
Line 
1function [maps_out,fdrth,H] = fmri_analysis(DATA,varargin)
2%
3% This function does a quick fMRI analysis in SPM style. No motion
4% correction or slice timing at this time...
5%
6
7% Defaults
8TR = [];
9contr = [1];
10onset = [];
11durat = [];
12smooth_kernel = [];
13d_cutoff = [];      % Don't temporally filter data by default
14qFDR = 0.05;        % Default threshold for FDR correction
15globalNorm = false; % Don't do global normalization by default
16omitVols = [];      % Don't omit any volumes from the analysis by default
17show_wbar = true;   % Show waitbar by default
18mask = [];          % Don't mask by default
19regr = [];          % Don't use additional regressors by default
20
21if rem(length(varargin),2)~=0
22  error('parameters/values missing!')
23end
24
25% Parse varargin
26for ii=1:2:length(varargin)
27  param = lower(varargin{ii});
28  value = varargin{ii+1};
29  switch param
30    case 'tr'
31      TR = value;
32    case {'onset','onsets'}
33      onset = value;
34      onset = onset(:);
35    case {'durat','durations','duration'}
36      durat = value;
37      durat = durat(:);
38    case 'smooth'
39      smooth_kernel = value;
40    case 'contrast'
41      contr = value;
42    case 'fdrth'
43      qFDR = value;
44    case {'hipass','hicutoff'}
45      d_cutoff = value;
46    case 'omitvolumes'
47      omitVols = value;
48    case 'wbar'
49      if value==1
50        show_wbar = true;
51      else
52        show_wbar = false;
53      end
54    case 'mask'
55      if isnumeric(value) || islogical(value)
56        mask = value;
57      else
58        error('Invalid mask!')
59      end
60    case 'regressor'
61      regr = value;
62  end
63end
64
65% Check that we have something to fit...
66if isempty(onset) && isempty(regr)
67        warning('No onsets or regressors defined. Only mean will be fitted!')
68end
69
70% Check that TR has been given
71if isempty(TR)
72        error('Repetition time (TR) has not been defined!')
73end
74
75% Check onsets and durations
76if ~isempty(onset)
77        if isempty(durat)
78                error('Durations for onsets have not been defined!')
79        end
80  if length(durat)==1
81    durat = repmat(durat,1,length(onset));
82    durat = durat(:);
83  elseif length(durat)~=length(onset)
84    error('Mismatch in the number of elements in onset and duration!')
85  end
86end
87
88% Check data
89if isstruct(DATA) && isfield(DATA,'FTDATA')
90  data = DATA.FTDATA;
91elseif isnumeric(DATA) && ndims(DATA)>2
92  data=DATA;
93else
94  error('Input data has to be a 3D/4D numeric matrix or Aedes data structure!')
95end
96
97% Permute 3D matrix to 4D
98if ndims(data)==3
99  data = permute(data,[1 2 4 3]);
100end
101
102% Initialize mask
103if isempty(mask)
104  mask = true(size(data,1),size(data,2),size(data,3));
105end
106
107% Check regressors
108if ~isempty(regr)
109  if size(regr,1)~=size(data,4)
110    error('The lengths of the regressors do not match data length!')
111        end
112  regr = detrend(regr);
113end
114
115% Load Rat HRF (8 seconds long) NOTE: DON'T USE FOR HUMAN DATA!!!
116bf = [
1170
118   0.000008876945146
119   0.000387621003097
120   0.003012441669394
121   0.011548214895067
122   0.030056522080168
123   0.061233628208021
124   0.105349995559045
125   0.160158500633323
126   0.221528149631163
127   0.284405279593501
128   0.343761753554314
129   0.395323573307322
130   0.436002661868442
131   0.464045998242252
132   0.478965423658889
133   0.481327079129854
134   0.472473786978796
135   0.454237528221769
136   0.428680153860941
137   0.397883140401152
138   0.363793656845973
139   0.328124931280142
140   0.292303457117353
141   0.257453130446147
142   0.224406054242129
143   0.193730677047637
144   0.165769521647024
145   0.140680556701329
146   0.118477987483032
147   0.099069734765398
148   0.082290068881472
149   0.067926763712502
150   0.055742762013807
151   0.045492744488986
152   0.036935220196777
153   0.029840852217657
154   0.023997740489123
155   0.019214335904674
156   0.015320580887648
157   0.012167779438633
158   0.009627606069476
159   0.007590575511228
160   0.005964217696074
161   0.004671136982604
162   0.003647081075751
163   0.002839102788446
164   0.002203865389816
165   0.001706118264261
166   0.001317352436400
167   0.001014633776781
168   0.000779604140824
169   0.000597636251764
170   0.000457125954290
171   0.000348904854607
172   0.000265756797153
173   0.000202022712087
174   0.000153279812313
175   0.000116082720651
176   0.000087755728843
177   0.000066226941806
178   0.000049896490408
179   0.000037532277385
180];
181
182fprintf(1,'\n******************************************\n');
183fprintf(1,'Starting fMRI analysis.\n');
184if isstruct(DATA) && isfield(DATA,'HDR') && isfield(DATA.HDR,'fpath')
185  filename = [DATA.HDR.fpath,DATA.HDR.fname];
186  fprintf(1,'File: %s\n\n',filename)
187else
188  fprintf(1,'\n');
189end
190
191% Omit volumes from data and stimulus function if requested
192if ~isempty(omitVols)
193 
194  fprintf(1,'Skipping requested volumes...\n');
195  if ~isempty(onset)
196    % Calculate new onsets and durations
197    ton = onset;
198    tof = onset+durat+1;
199    tmp=zeros(1,size(data,4));
200    tmp(ton)=1;tmp(tof)=-1;
201    sf=cumsum(tmp);
202    sf(omitVols)=[];
203    tmp=diff([0 sf]);
204    onset = find(tmp==1);
205    durat = find(tmp==-1);
206    if isempty(durat)
207      % Block stays up...
208      durat=length(tmp)-onset(end)-1;
209    elseif length(durat) < length(onset)
210      durat(1:end-1) = durat(1:end-1)-onset(1:end-1)-1;
211      durat(end) = length(tmp)-onset(end)-1;
212    else
213      durat = durat-onset-1;
214    end
215    onset=onset(:);
216    durat=durat(:);
217   
218    fprintf(1,['New onsets: ',num2str(onset(:)'),'\n']);
219    fprintf(1,['New duration: ',num2str(durat(:)'),'\n']);
220  end
221 
222  if ~isempty(regr)
223    regr(omitVols,:)=[];
224  end
225  data(:,:,:,omitVols) = [];
226   
227end
228
229% Create stimulus function (32 bin offset)
230k = size(data,4); % Number of scans
231T     = 16;
232dt    = TR/T;
233if ~isempty(onset)
234  u     = onset.^0;
235  if ~any(durat)
236    u  = u/dt;
237  end
238  ton       = round(onset*TR/dt) + 32;                  % onsets
239  tof       = round(durat*TR/dt) + ton + 1;                     % offset
240  sf        = zeros((k*T + 128),size(u,2));
241  ton       = max(ton,1);
242  tof       = max(tof,1);
243  for j = 1:length(ton)
244    if numel(sf)>ton(j),
245      sf(ton(j),:) = sf(ton(j),:) + u(j,:);
246    end;
247    if numel(sf)>tof(j),
248      sf(tof(j),:) = sf(tof(j),:) - u(j,:);
249    end;
250  end
251  sf        = cumsum(sf);                                       % integrate
252 
253  % Convolve stimulus with the HRF
254  conv_sf = conv(sf,bf);
255 
256  % Resample the convolved stimulus
257  RR = conv_sf([0:(k-1)]*16+1+32);
258  RR=RR(:);
259else
260  RR=[];
261end
262
263
264% Smooth data if requested
265if ~isempty(smooth_kernel) && ~ismember(smooth_kernel,[1 1 1],'rows')
266        if show_wbar
267                wbh = aedes_calc_wait('Smoothing data...');
268                drawnow
269        end
270        if all(smooth_kernel)
271                tic
272                smooth_data = fmri_smooth(data,smooth_kernel);
273                toc
274        else
275                fprintf(1,'fMRI analysis warning: Could not smooth data!\n');
276                smooth_data = data;
277        end
278        if show_wbar
279                close(wbh);
280        end
281else
282        smooth_data = data;
283end
284
285
286if ~isempty(d_cutoff)
287  % Create filtering matrix
288  nCosine = ceil((2*k*TR)/(d_cutoff + 1));
289  S = sqrt(2/k)*cos([1:nCosine]'*pi*([1:k]/k)).';
290  KKT = eye(size(S,1))-2*S*S'+S*S'*S*S';
291else
292  KKT = eye(k);
293  S = 0;
294end
295H = [RR regr ones(k,1)];
296H = H-S*(S'*H); % Filter design matrix
297nParam = size(H,2);
298maps_out = struct('pmap',[],'tmap',[]);
299maps_out.pmap = zeros(size(smooth_data,1),size(smooth_data,2),size(smooth_data,3),nParam);
300maps_out.tmap = zeros(size(smooth_data,1),size(smooth_data,2),size(smooth_data,3));
301
302% Calculate parametric map(s)
303nPlanes = size(smooth_data,3);
304nCols = size(smooth_data,2);
305nRows = size(smooth_data,1);
306if length(contr)<nParam
307        contr(nParam)=0; % Pad with zeros
308end
309c = repmat(contr,nRows,1);
310HTH = pinv(H'*H);
311R = eye(size(H,1))-H*HTH*H';
312
313if show_wbar
314  wbh = aedes_wbar(0,sprintf('Estimating parameters. Processing plane 0/%d',nPlanes));
315  drawnow
316end
317
318% Process data in columns
319for ii=1:nPlanes
320 
321  %fprintf(1,'Processing plane %d/%d\n',ii,nPlanes);
322  for kk=1:nCols
323    if show_wbar
324      aedes_wbar(ii/nPlanes,wbh,sprintf('Estimating parameters. Processing plane %d/%d, column %d/%d',ii,nPlanes,kk,nCols));
325    end
326    col_data = squeeze(smooth_data(:,kk,ii,:)).';
327    col_data = col_data-S*(S'*col_data);
328    th=H\col_data;
329    r = col_data-H*th;
330    rr=diag(r'*r);
331    rr=rr(:);
332    sig2 = rr./trace(R*KKT);
333    sig2 = repmat(sig2,1,nParam);
334    T = diag(c*th)./sqrt(c.*sig2*HTH*H'*KKT*H*HTH*contr');
335    T(find(mask(:,kk,ii)==0))=0;
336                maps_out.pmap(:,kk,ii,:) = th.';
337                maps_out.tmap(:,kk,ii)=T.';
338        end
339end
340if show_wbar
341  close(wbh);
342end
343
344if show_wbar
345  wbh = aedes_calc_wait('Calculating threshold...');
346  drawnow
347end
348
349% Set NaN:s to zeros
350maps_out.tmap(isnan(maps_out.tmap)) = 0;
351
352% Calculate effective degrees of freedom
353dof = (trace(R*KKT).^2)/trace(R*KKT*R*KKT);
354
355% p-values
356pval_map = 1-tdist(maps_out.tmap,dof);
357
358% Perform FDR (False Discovery Rate) correction
359cV = 1;
360
361pValues = pval_map(:);
362tValues = maps_out.tmap(:);
363
364[pValuesSorted,sortInd] = sort(pValues);
365tValuesSorted = tValues(sortInd);
366nP = length(pValues);
367
368pFDR = [1:nP]'/nP*qFDR/cV; % FDR-correction
369thresFDRind = find(pValuesSorted<=pFDR,1,'last');
370if ~isempty(thresFDRind)
371  thresFDR = tValuesSorted(thresFDRind);
372else
373  thresFDR = [];
374end
375
376if show_wbar
377  close(wbh);
378end
379
380if nargout>=2
381  fdrth = thresFDR;
382end
383
384if ~isempty(thresFDR)
385  fprintf(1,['FDR threshold at p<',num2str(qFDR),': %.3f\n'],thresFDR)
386else
387  fprintf(1,['No significant voxels over FDR threshold at p<',num2str(qFDR),'!\n'])
388end
389
390fprintf(1,'******************************************\n');
391
392
Note: See TracBrowser for help on using the repository browser.

Powered by Trac 1.0.9.Copyright © Juha-Pekka Niskanen 2008