source: plugins/fmri_plugins/basic_fmri_analysis.m @ 82

Last change on this file since 82 was 82, checked in by tjniskan, 10 years ago
  • Some minor modifications to aedes_readfid.m
  • Added BOLD strength estimation to

plugins/fmri_plugins/basic_fmri_analysis.m

  • Changed the image in about dialog
  • Added support for negative threshold in misclib/fmri_blob_overlay.m

M misclib/fmri_blob_overlay.m
M aedes_helpabout.m
M plugins/fmri_plugins/basic_fmri_analysis.m
M aedes_readfid.m
M aedes_revision.m

File size: 9.1 KB
Line 
1function basic_fmri_analysis(DATA,ROI,AddInfo)
2%
3% This Aedes plugin does a quick basic fMRI analysis.
4%
5
6if AddInfo.isDataMixed
7  % Make data a 4D-matrix if it is loaded in aedes as a stack of 2D or 3D images
8  DATA2{1}=DATA{1};
9  DATA2{1}.FTDATA = zeros([size(DATA{1}.FTDATA,1),...
10    size(DATA{1}.FTDATA,2),size(DATA{1}.FTDATA,3),length(DATA)],...
11    class(DATA{1}.FTDATA));
12  for ii=1:length(DATA)
13    DATA2{1}.FTDATA(:,:,:,ii)=DATA{ii}.FTDATA;
14  end
15  DATA=DATA2;
16end
17
18% Prompt for design and contrast
19prompt={'Stimulus onsets (scans):','Stimulus durations (scans):',...
20  'Contrast vector','Smooth FWHM','TR','High-pass filter cutoff (sec.)',...
21  'P-value for FDR correction','Volumes to omit from analysis'};
22name='Enter parameters for fMRI analysis';
23numlines=1;
24defaultanswer={'30 75 120','15 15 15','1 0',...
25  '2 2 1','2.039','256','0.05',''};
26
27resp=inputdlg(prompt,name,numlines,defaultanswer);
28if isempty(resp)
29  % Canceled
30  return
31end
32
33% Call fmri_analysis in /misclib
34onset = str2num(resp{1});
35durat = str2num(resp{2});
36contr = str2num(resp{3});
37if length(contr)~=2
38  errordlg('Invalid contrast vector!','Error','modal');
39  return
40end
41smooth_sz = str2num(resp{4});
42TR = str2num(resp{5});
43hipass = str2num(resp{6});
44qFDR = str2num(resp{7});
45omitVols = str2num(strrep(resp{8},'end',num2str(size(DATA{1}.FTDATA,4))));
46
47% if isfield(DATA{1},'PROCPAR') && isfield(DATA{1}.PROCPAR,'readres')
48%   % Remove reference image from VNMR EPI data
49%   DATA{1}.FTDATA = DATA{1}.FTDATA(:,:,:,2:end);
50% end
51
52if length(durat)==1
53  durat = ones(1,length(onset))*durat;
54end
55
56if ~isempty(ROI) && any(strcmpi({ROI(:).label},'mask'))
57  roi_ind = find(strcmpi({ROI(:).label},'mask'));
58  [maps_out,fdr_th] = fmri_analysis(DATA{1},'TR',TR,...
59    'onsets',onset,'durations',durat,...
60    'contrast',contr,'smooth',smooth_sz,'FDRth',qFDR,...
61    'hipass',hipass,'omitvolumes',omitVols,...
62    'mask',ROI(roi_ind).voxels{1}(:,:,:,1));
63  ROI(roi_ind)=[];
64else
65  [maps_out,fdr_th] = fmri_analysis(DATA{1},'TR',TR,...
66    'onsets',onset,'durations',durat,...
67    'contrast',contr,'smooth',smooth_sz,'FDRth',qFDR,...
68    'hipass',hipass,'omitvolumes',omitVols);
69end
70
71% Resize tmap to EPI time series size
72tmap = maps_out.tmap;
73tmap = repmat(tmap,[1,1,1,size(DATA{1}.FTDATA,4)]);
74
75overlay.FTDATA = tmap;
76overlay.ImageOverlayCmapStr = 'hot';
77overlay.ImageOverlayTholdDirPos = 1;
78overlay.ImageOverlayAlpha = 1;
79overlay.fMRIonsets = onset;
80overlay.fMRIdurats = durat;
81if ~isempty(fdr_th)
82  overlay.ImageOverlayThold = fdr_th;
83  overlay.ImOverlayMin = 0;
84else
85  overlay.ImageOverlayAlpha = 0;
86end
87
88% If some volumes were omitted, calculate new onsets and durations
89if ~isempty(omitVols)
90  % create stick function
91  ton = onset;
92  tof = onset+durat+1;
93  tmp=zeros(1,size(DATA{1}.FTDATA,4));
94  tmp(ton)=1;tmp(tof)=-1;
95  sf=cumsum(tmp);
96  sf(omitVols)=[];
97  tmp=diff([0 sf]);
98  new_onset = find(tmp==1);
99  new_durat = find(tmp==-1)-new_onset-1;
100  overlay.fMRIonsets = new_onset;
101  overlay.fMRIdurats = new_durat;
102  DATA{1}.FTDATA(:,:,:,omitVols)=[];
103  overlay.FTDATA(:,:,:,omitVols)=[];
104else
105  new_onset = onset;
106  new_durat = durat;
107end
108
109
110% Estimate BOLD strength for ROIs
111if ~isempty(ROI)
112  nScans = size(DATA{1}.FTDATA,4);
113  fprintf(1,'---------------------------------------\n');
114  for kk=1:length(ROI)
115   
116    % Get time-series indices from ROI
117    ind = repmat(ROI(kk).voxels{1}(:,:,:,1),[1 1 1 size(DATA{1}.FTDATA,4)]);
118
119    % Mean EPI time series
120    ts_data = reshape(double(DATA{1}.FTDATA(ind)),[],size(DATA{1}.FTDATA,4));
121    ts_data = ts_data.';
122    ts_data = detrend(ts_data)+repmat(mean(ts_data),size(ts_data,1),1);
123   
124    % Normalize mean to bold-%
125    mean_ts_data = mean(ts_data,2);
126   
127    % estimate BOLD-% for the ROI mean
128    [bold_prc,z,th,z_hat]=estimate_bold(new_onset,new_durat,TR,nScans,mean_ts_data);
129   
130    z_norm = (z./th(2)-1)*100;
131    z_hat_norm = (z_hat./th(2)-1)*100;
132   
133    % Normalize raw timeseries to bold-%
134    %ts_data_bold = (ts_data./repmat(mean(ts_data(1:20,:)),size(ts_data,1),1)-1)*100;
135   
136    % Plot results
137    figure;
138    ax=axes;
139    line(1:length(z),z_norm,'color','k',...
140      'parent',ax);
141    line(1:length(z),z_hat_norm,...
142      'color','r','linewidth',2,'parent',ax);
143    title(['Time series for ROI: ',ROI(kk).label])
144    ylabel('BOLD-%');
145    set(ax,'xlim',[0 length(z)],...
146      'ylim',[min(z_norm)-min(z_norm)*0.05 ...
147      max(z_norm)+max(z_norm)*0.05]);
148   
149    fprintf(1,'BOLD-%% %s: %.3f\n',ROI(kk).label,bold_prc);
150  end
151  fprintf(1,'---------------------------------------\n');
152else
153  % Estimate BOLD strength for the time-series corresponding to the
154  % maximum value in T-map.
155  nScans = size(DATA{1}.FTDATA,4);
156  [mx,I_mx]=max(maps_out.tmap(:));
157  [II_mx,JJ_mx,KK_mx]=ind2sub(size(maps_out.tmap),I_mx);
158  [mn,I_mn]=min(maps_out.tmap(:));
159  [II_mn,JJ_mn,KK_mn]=ind2sub(size(maps_out.tmap),I_mn);
160 
161  ts_data_max = squeeze(DATA{1}.FTDATA(II_mx,JJ_mx,KK_mx,:));
162  ts_data_min = squeeze(DATA{1}.FTDATA(II_mn,JJ_mn,KK_mn,:));
163 
164  [bold_prc_max,z_max,th_max,z_hat_max]=estimate_bold(new_onset,new_durat,TR,nScans,ts_data_max);
165  [bold_prc_min,z_min,th_min,z_hat_min]=estimate_bold(new_onset,new_durat,TR,nScans,ts_data_min);
166 
167  z_norm_max = (z_max./th_max(2)-1)*100;
168  z_hat_norm_max = (z_hat_max./th_max(2)-1)*100;
169 
170  z_norm_min = (z_min./th_min(2)-1)*100;
171  z_hat_norm_min= (z_hat_min./th_min(2)-1)*100;
172 
173  % Plot results
174  figure;
175  ax1=subplot(2,1,1);
176  line(1:length(z_max),z_norm_max,'color','k',...
177    'parent',ax1);
178  line(1:length(z_max),z_hat_norm_max,...
179    'color','r','linewidth',2,'parent',ax1);
180  title('Time series for Max T-value');
181  ylabel('BOLD-%');
182  set(ax1,'xlim',[0 length(z_max)],...
183    'ylim',[min(z_norm_max)-min(z_norm_max)*0.05 ...
184    max(z_norm_max)+max(z_norm_max)*0.05]);
185 
186  ax2=subplot(2,1,2);
187  line(1:length(z_min),z_norm_min,'color','k',...
188    'parent',ax2);
189  line(1:length(z_min),z_hat_norm_min,...
190    'color','r','linewidth',2,'parent',ax2);
191  title('Time series for Min T-value');
192  ylabel('BOLD-%');
193  set(ax2,'xlim',[0 length(z_min)],...
194    'ylim',[min(z_norm_min)-min(z_norm_min)*0.05 ...
195    max(z_norm_min)+max(z_norm_min)*0.05]);
196 
197  % Print percents to command window
198  fprintf(1,'---------------------------------------\n');
199  fprintf(1,'Max T BOLD-%%: %.3f\n',bold_prc_max);
200  fprintf(1,'Min T BOLD-%%: %.3f\n',bold_prc_min);
201  fprintf(1,'---------------------------------------\n');
202end
203
204% Show in Aedes as overlay
205aedes(DATA,[],overlay);
206
207if isempty(fdr_th)
208  warndlg(['No significant voxels over FDR threshold at p < ',...
209    num2str(qFDR),'!'],'No significant voxels!','modal')
210end
211
212
213function [bold_prc,z,th,z_hat]=estimate_bold(onset,durat,TR,num_scans,ts_data)
214%% -------------------------------------------------
215
216% Load Rat HRF (8 seconds long) NOTE: DON'T USE FOR HUMAN DATA!!!
217bf = [
218  0
219  0.000008876945146
220  0.000387621003097
221  0.003012441669394
222  0.011548214895067
223  0.030056522080168
224  0.061233628208021
225  0.105349995559045
226  0.160158500633323
227  0.221528149631163
228  0.284405279593501
229  0.343761753554314
230  0.395323573307322
231  0.436002661868442
232  0.464045998242252
233  0.478965423658889
234  0.481327079129854
235  0.472473786978796
236  0.454237528221769
237  0.428680153860941
238  0.397883140401152
239  0.363793656845973
240  0.328124931280142
241  0.292303457117353
242  0.257453130446147
243  0.224406054242129
244  0.193730677047637
245  0.165769521647024
246  0.140680556701329
247  0.118477987483032
248  0.099069734765398
249  0.082290068881472
250  0.067926763712502
251  0.055742762013807
252  0.045492744488986
253  0.036935220196777
254  0.029840852217657
255  0.023997740489123
256  0.019214335904674
257  0.015320580887648
258  0.012167779438633
259  0.009627606069476
260  0.007590575511228
261  0.005964217696074
262  0.004671136982604
263  0.003647081075751
264  0.002839102788446
265  0.002203865389816
266  0.001706118264261
267  0.001317352436400
268  0.001014633776781
269  0.000779604140824
270  0.000597636251764
271  0.000457125954290
272  0.000348904854607
273  0.000265756797153
274  0.000202022712087
275  0.000153279812313
276  0.000116082720651
277  0.000087755728843
278  0.000066226941806
279  0.000049896490408
280  0.000037532277385
281  ];
282
283
284%TR=2.039;
285%onset=[30 75 120]';
286%durat=[15 15 15]';
287onset = onset(:);
288durat = durat(:);
289
290%bf = spm_hrf(1/16);
291
292% Create stimulus function (32 bin offset)
293%k = 165; % Number of scans
294k = num_scans;
295T     = 16;
296dt    = TR/T;
297u     = ones(size(onset));
298if ~any(durat)
299  u  = u/dt;
300end
301ton = round(onset*TR/dt) + 32;                  % onsets
302tof = round(durat*TR/dt) + ton + 1;                     % offset
303sf = zeros((k*T + 128),size(u,2));
304
305for j = 1:length(ton)
306  if numel(sf)>ton(j),
307    sf(ton(j),:) = sf(ton(j),:) + u(j,:);
308  end;
309  if numel(sf)>tof(j),
310    sf(tof(j),:) = sf(tof(j),:) - u(j,:);
311  end;
312end
313sf        = cumsum(sf);                                 % integrate
314
315% Convolve stimulus with the HRF
316conv_sf = conv(sf,bf);
317
318% Resample the convolved stimulus
319R = conv_sf([0:(k-1)]*16+1+32);
320R=R(:);
321
322% Remove linear trend but keep the mean
323mean_ts_data = mean(ts_data);
324ts_data = detrend(ts_data)+mean_ts_data;
325
326% Estimate baseline from the paradigm
327bl_ind = find(sf==0);
328bl_ind = bl_ind(1:min(length(bl_ind),20));
329if isempty(bl_ind)
330  baseline = mean(ts_data);
331else
332  baseline = mean(ts_data(bl_ind));
333end
334
335mean_ts_data_bold = (ts_data./baseline-1)*100;
336
337
338% Fit block model to the data
339%z=mean_ts_data_bold;
340z=ts_data;
341H = [R ones(size(R))];
342th = H\z;
343z_hat = H*th;
344
345% Calculate results
346bold_prc = ((th(1)*max(R))/th(2))*100;
347
Note: See TracBrowser for help on using the repository browser.

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