source: misclib/fmri_analysis.m @ 94

Last change on this file since 94 was 94, checked in by tjniskan, 10 years ago
  • Fixed a bug in reading Varian FID-files (regression bug)
  • Added a possibility of Using 3D overlays on 4D data in Aedes
  • misclib/fmri_analysis.m now accepts general regressors

M misclib/fmri_analysis.m
M aedes.m
M aedes_readfid.m
M aedes_revision.m

File size: 8.4 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 = 2.039;
9contr = [1 0];
10onset = [30 75 120].';
11durat = [15 15 15].';
12smooth_kernel = [2 2 1];
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'}
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
65if ~isempty(onset)
66  if length(durat)==1
67    durat = repmat(durat,1,length(onset));
68    durat = durat(:);
69  elseif length(durat)~=length(onset)
70    error('Mismatch in the number of elements in onset and duration!')
71  end
72end
73
74
75% Check data
76if isstruct(DATA) && isfield(DATA,'FTDATA')
77  data = DATA.FTDATA;
78elseif isnumeric(DATA) && ndims(DATA)>2
79  data=DATA;
80else
81  error('Input data has to be 3D or 4D numeric matrix!')
82end
83
84if ndims(data)==3
85  data = permute(data,[1 2 4 3]);
86end
87
88if isempty(mask)
89  mask = true(size(data,1),size(data,2),size(data,3));
90end
91
92% Check regressors
93if ~isempty(regr)
94  if size(regr,1)~=size(data,4) && size(regr,2)~=size(data,4)
95    error('The lengths of the regressors do not match data length!')
96  end
97  if size(regr,1)~=size(data,4)
98    regr = regr.';
99  end
100  regr = detrend(regr);
101end
102
103% Load Rat HRF (8 seconds long) NOTE: DON'T USE FOR HUMAN DATA!!!
104bf = [
1050
106   0.000008876945146
107   0.000387621003097
108   0.003012441669394
109   0.011548214895067
110   0.030056522080168
111   0.061233628208021
112   0.105349995559045
113   0.160158500633323
114   0.221528149631163
115   0.284405279593501
116   0.343761753554314
117   0.395323573307322
118   0.436002661868442
119   0.464045998242252
120   0.478965423658889
121   0.481327079129854
122   0.472473786978796
123   0.454237528221769
124   0.428680153860941
125   0.397883140401152
126   0.363793656845973
127   0.328124931280142
128   0.292303457117353
129   0.257453130446147
130   0.224406054242129
131   0.193730677047637
132   0.165769521647024
133   0.140680556701329
134   0.118477987483032
135   0.099069734765398
136   0.082290068881472
137   0.067926763712502
138   0.055742762013807
139   0.045492744488986
140   0.036935220196777
141   0.029840852217657
142   0.023997740489123
143   0.019214335904674
144   0.015320580887648
145   0.012167779438633
146   0.009627606069476
147   0.007590575511228
148   0.005964217696074
149   0.004671136982604
150   0.003647081075751
151   0.002839102788446
152   0.002203865389816
153   0.001706118264261
154   0.001317352436400
155   0.001014633776781
156   0.000779604140824
157   0.000597636251764
158   0.000457125954290
159   0.000348904854607
160   0.000265756797153
161   0.000202022712087
162   0.000153279812313
163   0.000116082720651
164   0.000087755728843
165   0.000066226941806
166   0.000049896490408
167   0.000037532277385
168];
169
170fprintf(1,'\n******************************************\n');
171fprintf(1,'Starting fMRI analysis.\n');
172if isstruct(DATA) && isfield(DATA,'HDR') && isfield(DATA.HDR,'fpath')
173  filename = [DATA.HDR.fpath,DATA.HDR.fname];
174  fprintf(1,'File: %s\n\n',filename)
175else
176  fprintf(1,'\n');
177end
178
179% Omit volumes from data and stimulus function if requested
180if ~isempty(omitVols)
181 
182  fprintf(1,'Skipping requested volumes...\n');
183  if ~isempty(onset)
184    % Calculate new onsets and durations
185    ton = onset;
186    tof = onset+durat+1;
187    tmp=zeros(1,size(data,4));
188    tmp(ton)=1;tmp(tof)=-1;
189    sf=cumsum(tmp);
190    sf(omitVols)=[];
191    tmp=diff([0 sf]);
192    onset = find(tmp==1);
193    durat = find(tmp==-1);
194    if isempty(durat)
195      % Block stays up...
196      durat=length(tmp)-onset(end)-1;
197    elseif length(durat) < length(onset)
198      durat(1:end-1) = durat(1:end-1)-onset(1:end-1)-1;
199      durat(end) = length(tmp)-onset(end)-1;
200    else
201      durat = durat-onset-1;
202    end
203    onset=onset(:);
204    durat=durat(:);
205   
206    fprintf(1,['New onsets: ',num2str(onset(:)'),'\n']);
207    fprintf(1,['New duration: ',num2str(durat(:)'),'\n']);
208  end
209 
210  if ~isempty(regr)
211    regr(omitVols,:)=[];
212  end
213  data(:,:,:,omitVols) = [];
214   
215end
216
217% Create stimulus function (32 bin offset)
218k = size(data,4); % Number of scans
219T     = 16;
220dt    = TR/T;
221if ~isempty(onset)
222  u     = onset.^0;
223  if ~any(durat)
224    u  = u/dt;
225  end
226  ton       = round(onset*TR/dt) + 32;                  % onsets
227  tof       = round(durat*TR/dt) + ton + 1;                     % offset
228  sf        = zeros((k*T + 128),size(u,2));
229  ton       = max(ton,1);
230  tof       = max(tof,1);
231  for j = 1:length(ton)
232    if numel(sf)>ton(j),
233      sf(ton(j),:) = sf(ton(j),:) + u(j,:);
234    end;
235    if numel(sf)>tof(j),
236      sf(tof(j),:) = sf(tof(j),:) - u(j,:);
237    end;
238  end
239  sf        = cumsum(sf);                                       % integrate
240 
241  % Convolve stimulus with the HRF
242  conv_sf = conv(sf,bf);
243 
244  % Resample the convolved stimulus
245  R = conv_sf([0:(k-1)]*16+1+32);
246  R=R(:);
247else
248  R=[];
249end
250
251
252% Smooth data
253if show_wbar
254  wbh = aedes_calc_wait('Smoothing data...');
255  drawnow
256end
257if all(smooth_kernel)
258  smooth_data = fmri_smooth(data,smooth_kernel);
259else
260  fprintf(1,'fMRI analysis warning: Could not smooth data!\n');
261  smooth_data = data;
262end
263
264if ~isempty(d_cutoff)
265  % Create filtering matrix
266  nCosine = ceil((2*k*TR)/(d_cutoff + 1));
267  S = sqrt(2/k)*cos([1:nCosine]'*pi*([1:k]/k)).';
268  KKT = eye(size(S,1))-2*S*S'+S*S'*S*S';
269else
270  KKT = eye(k);
271  S = 0;
272end
273H = [R regr ones(k,1)];
274H = H-S*(S'*H); % Filter design matrix
275
276% Calculate parametric map(s)
277nPlanes = size(smooth_data,3);
278nCols = size(smooth_data,2);
279nRows = size(smooth_data,1);
280pmap1 = zeros(size(smooth_data,1),size(smooth_data,2),size(smooth_data,3));
281pmap2 = zeros(size(pmap1));
282tmap = zeros(size(pmap1));
283c = repmat(contr,nRows,1);
284HTH = pinv(H'*H);
285R = eye(size(H,1))-H*HTH*H';
286if show_wbar
287  close(wbh);
288end
289if show_wbar
290  wbh = aedes_wbar(0,sprintf('Estimating parameters. Processing plane 0/%d',nPlanes));
291  drawnow
292end
293for ii=1:nPlanes
294 
295  %fprintf(1,'Processing plane %d/%d\n',ii,nPlanes);
296  for kk=1:nCols
297    if show_wbar
298      aedes_wbar(ii/nPlanes,wbh,sprintf('Estimating parameters. Processing plane %d/%d, column %d/%d',ii,nPlanes,kk,nCols));
299    end
300    col_data = squeeze(smooth_data(:,kk,ii,:)).';
301    col_data = col_data-S*(S'*col_data);
302    th=H\col_data;
303    pmap1(:,kk,ii)=th(1,:).';
304    pmap2(:,kk,ii)=th(2,:).';
305   
306    r = col_data-H*th;
307    rr=diag(r'*r);
308    rr=rr(:);
309    sig2 = rr./trace(R*KKT);
310    sig2 = repmat(sig2,1,2);
311    T = diag(c*th)./sqrt(c.*sig2*HTH*H'*KKT*H*HTH*contr');
312    T(find(mask(:,kk,ii)==0))=0;
313    tmap(:,kk,ii)=T.';
314  end
315end
316if show_wbar
317  close(wbh);
318end
319
320if show_wbar
321  wbh = aedes_calc_wait('Calculating threshold...');
322  drawnow
323end
324
325% Calculate effective degrees of freedom
326dof = (trace(R*KKT).^2)/trace(R*KKT*R*KKT);
327
328% p-values
329pval_map = 1-tdist(tmap,dof);
330
331% Perform FDR (False Discovery Rate) correction
332cV = 1;
333
334pValues = pval_map(:);
335tValues = tmap(:);
336
337[pValuesSorted,sortInd] = sort(pValues);
338tValuesSorted = tValues(sortInd);
339nP = length(pValues);
340
341pFDR = [1:nP]'/nP*qFDR/cV; % FDR-correction
342thresFDRind = find(pValuesSorted<=pFDR,1,'last');
343%if isempty(thresFDRind)
344%       fprintf(1,'Warning: No significant values with the FDR-corrected threshold %f\n',...
345%               qFDR);
346%       return
347%end
348if ~isempty(thresFDRind)
349  thresFDR = tValuesSorted(thresFDRind);
350else
351  thresFDR = [];
352end
353
354maps_out.pmap1 = pmap1;
355maps_out.pmap2 = pmap2;
356%maps_out.pmap3 = pmap3;
357maps_out.tmap = tmap;
358
359if show_wbar
360  close(wbh);
361end
362
363if nargout>=2
364  fdrth = thresFDR;
365end
366
367if ~isempty(thresFDR)
368  fprintf(1,['FDR threshold at p<',num2str(qFDR),': %.3f\n'],thresFDR)
369else
370  fprintf(1,['No significant voxels over FDR threshold at p<',num2str(qFDR),'!\n'])
371end
372
373fprintf(1,'******************************************\n');
374
375
Note: See TracBrowser for help on using the repository browser.

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