source: plugins/fmri_plugins/basic_fmri_analysis.m

Last change on this file was 152, checked in by tjniskan, 8 years ago
  • aedes_readvnmr.m almost works with epi data...
  • The ROI averages of fMRI time series are now shown in one column instead of two in plugins/fmri_plugins/basic_fmri_analysis.m

M aedes_readvnmr.m
M plugins/fmri_plugins/basic_fmri_analysis.m
M aedes_revision.m
M vnmr_recon/epi_recon.m

File size: 11.4 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;
73%tmap = 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  fh=figure;
115        nRows = length(ROI);
116        nCols = 1;
117        ax_w = 0.9;
118        ax_l = (1-ax_w)/2;
119        ax_gap = 0.01;
120        ax_h = (1-(length(ROI)+5)*ax_gap)/length(ROI);
121        bgax = axes('parent',fh,...
122                'units','normal','position',[0 0 1 1],...
123                'xlim',[0 1],'ylim',[0 1],'visible','off');
124%   if length(ROI)<=3
125%     nRows = length(ROI);
126%     nCols=1;
127%   elseif length(ROI)==4
128%     nRows = 2;
129%     nCols=2;
130%   elseif length(ROI)<10
131%     nRows = 3;
132%     nCols = ceil(length(ROI)/nRows);
133%   else
134%     nRows = 4;
135%     nCols = ceil(length(ROI)/nRows);
136%   end
137  for kk=1:length(ROI)
138   
139    % Get time-series indices from ROI
140    ind = repmat(ROI(kk).voxels{1}(:,:,:,1),[1 1 1 size(DATA{1}.FTDATA,4)]);
141
142    % Mean EPI time series
143    ts_data = reshape(double(DATA{1}.FTDATA(ind)),[],size(DATA{1}.FTDATA,4));
144    ts_data = ts_data.';
145    ts_data = detrend(ts_data)+repmat(mean(ts_data),size(ts_data,1),1);
146   
147    % Normalize mean to bold-%
148    mean_ts_data = mean(ts_data,2);
149   
150    % estimate BOLD-% for the ROI mean
151    [bold_prc,z,th,z_hat]=estimate_bold(new_onset,new_durat,TR,nScans,mean_ts_data);
152   
153    z_norm = (z./th(2)-1)*100;
154    %z_hat_norm = (z_hat./th(2)-1)*100;
155   
156    z_hat_norm=aedes_trendest(double(z_norm),10);
157   
158    % Normalize raw timeseries to bold-%
159    %ts_data_bold = (ts_data./repmat(mean(ts_data(1:20,:)),size(ts_data,1),1)-1)*100;
160   
161    % Plot results
162    %ax=subplot(nRows,nCols,kk,'align','parent',fh,'fontsize',8);
163                ax=axes('parent',fh,'units','normal',...
164                        'position',[ax_l 1-kk*ax_h-kk*ax_gap ax_w ax_h],...
165                        'yaxislocation','right','layer','top','box','on');
166    line(1:length(z),z_norm,'color','k',...
167      'parent',ax);
168    line(1:length(z),z_hat_norm,...
169      'color','r','linewidth',2,'parent',ax);
170                text(ax_l-0.01,1-kk*ax_h-kk*ax_gap+ax_h/2,...
171                        ['ROI: ',ROI(kk).label],'parent',bgax,...
172                        'rotation',90,'fontsize',8,...
173                        'horizontalalign','center',...
174                        'verticalalign','bottom');
175    %title(['Time series for ROI: ',ROI(kk).label],'fontsize',8)
176    ylabel(ax,'BOLD-%','fontsize',8);
177    set(ax,'xlim',[0 length(z)],...
178      'ylim',[min(z_norm)-min(z_norm)*0.05 ...
179      max(z_norm)+max(z_norm)*0.05],'fontsize',8);
180    for ii=1:length(new_onset)
181      xdata = [new_onset(ii) new_onset(ii)+new_durat(ii) ...
182        new_onset(ii)+new_durat(ii) new_onset(ii)];
183      tmp = get(ax,'ylim');
184      ydata = [tmp(1) tmp(1) tmp(2) tmp(2)];
185      patch('parent',ax,...
186        'xdata',xdata,'ydata',ydata,'FaceColor','r',...
187        'FaceAlpha',0.3,'LineStyle','none');
188                end
189                if kk~=length(ROI)
190                        set(ax,'xticklabel',[])
191                end
192                set(ax,'layer','top')
193                if kk==length(ROI)
194                        xlabel('Time (scans)','fontsize',8)
195                end
196   
197    fprintf(1,'BOLD-%% %s: %.3f\n',ROI(kk).label,bold_prc);
198        end
199        if length(ROI)>=3
200                tmp_pos1 = get(0,'screensize');
201                tmp_pos2 = get(fh,'position');
202                set(fh,'position',[tmp_pos2(1) 100 tmp_pos2(3)*1.5 tmp_pos1(4)-100]);
203        end
204  fprintf(1,'---------------------------------------\n');
205else
206  % Estimate BOLD strength for the time-series corresponding to the
207  % maximum value in T-map.
208  nScans = size(DATA{1}.FTDATA,4);
209  [mx,I_mx]=max(maps_out.tmap(:));
210  [II_mx,JJ_mx,KK_mx]=ind2sub(size(maps_out.tmap),I_mx);
211  [mn,I_mn]=min(maps_out.tmap(:));
212  [II_mn,JJ_mn,KK_mn]=ind2sub(size(maps_out.tmap),I_mn);
213 
214  ts_data_max = double(squeeze(DATA{1}.FTDATA(II_mx,JJ_mx,KK_mx,:)));
215  ts_data_min = double(squeeze(DATA{1}.FTDATA(II_mn,JJ_mn,KK_mn,:)));
216 
217  [bold_prc_max,z_max,th_max,z_hat_max]=estimate_bold(new_onset,new_durat,TR,nScans,ts_data_max);
218  [bold_prc_min,z_min,th_min,z_hat_min]=estimate_bold(new_onset,new_durat,TR,nScans,ts_data_min);
219 
220  z_norm_max = (z_max./th_max(2)-1)*100;
221  z_hat_norm_max=aedes_trendest(double(z_norm_max),10);
222  %z_hat_norm_max = (z_hat_max./th_max(2)-1)*100;
223 
224  z_norm_min = (z_min./th_min(2)-1)*100;
225  z_hat_norm_min=aedes_trendest(double(z_norm_min),10);
226  %z_hat_norm_min= (z_hat_min./th_min(2)-1)*100;
227 
228  % Plot results
229  figure;
230  ax1=subplot(2,1,1,'align');
231  line(1:length(z_max),z_norm_max,'color','k',...
232    'parent',ax1);
233  line(1:length(z_max),z_hat_norm_max,...
234    'color','r','linewidth',2,'parent',ax1);
235  title('Time series for Max T-value');
236  ylabel('BOLD-%');
237  set(ax1,'xlim',[0 length(z_max)],...
238    'ylim',[min(z_norm_max)-min(z_norm_max)*0.05 ...
239    max(z_norm_max)+max(z_norm_max)*0.05]);
240
241  ax2=subplot(2,1,2,'align');
242  line(1:length(z_min),z_norm_min,'color','k',...
243    'parent',ax2);
244  line(1:length(z_min),z_hat_norm_min,...
245    'color','r','linewidth',2,'parent',ax2);
246  title('Time series for Min T-value');
247  ylabel('BOLD-%');
248  set(ax2,'xlim',[0 length(z_min)],...
249    'ylim',[min(z_norm_min)-min(z_norm_min)*0.05 ...
250    max(z_norm_min)+max(z_norm_min)*0.05]);
251 
252  for ii=1:length(new_onset)
253    xdata = [new_onset(ii) new_onset(ii)+new_durat(ii) ...
254      new_onset(ii)+new_durat(ii) new_onset(ii)];
255    tmp = get(ax1,'ylim');
256    ydata = [tmp(1) tmp(1) tmp(2) tmp(2)];
257    patch('parent',ax1,...
258      'xdata',xdata,'ydata',ydata,'FaceColor','r',...
259      'FaceAlpha',0.3,'LineStyle','none');
260    tmp = get(ax2,'ylim');
261    ydata = [tmp(1) tmp(1) tmp(2) tmp(2)];
262    patch('parent',ax2,...
263      'xdata',xdata,'ydata',ydata,'FaceColor','r',...
264      'FaceAlpha',0.3,'LineStyle','none');
265  end
266 
267  % Print percents to command window
268  fprintf(1,'---------------------------------------\n');
269  fprintf(1,'Max T BOLD-%%: %.3f\n',bold_prc_max);
270  fprintf(1,'Min T BOLD-%%: %.3f\n',bold_prc_min);
271  fprintf(1,'---------------------------------------\n');
272end
273
274% Show in Aedes as overlay
275aedes(DATA,[],overlay);
276
277if isempty(fdr_th)
278  warndlg(['No significant voxels over FDR threshold at p < ',...
279    num2str(qFDR),'!'],'No significant voxels!','modal')
280end
281
282
283function [bold_prc,z,th,z_hat]=estimate_bold(onset,durat,TR,num_scans,ts_data)
284%% -------------------------------------------------
285
286% Load Rat HRF (8 seconds long) NOTE: DON'T USE FOR HUMAN DATA!!!
287bf = [
288  0
289  0.000008876945146
290  0.000387621003097
291  0.003012441669394
292  0.011548214895067
293  0.030056522080168
294  0.061233628208021
295  0.105349995559045
296  0.160158500633323
297  0.221528149631163
298  0.284405279593501
299  0.343761753554314
300  0.395323573307322
301  0.436002661868442
302  0.464045998242252
303  0.478965423658889
304  0.481327079129854
305  0.472473786978796
306  0.454237528221769
307  0.428680153860941
308  0.397883140401152
309  0.363793656845973
310  0.328124931280142
311  0.292303457117353
312  0.257453130446147
313  0.224406054242129
314  0.193730677047637
315  0.165769521647024
316  0.140680556701329
317  0.118477987483032
318  0.099069734765398
319  0.082290068881472
320  0.067926763712502
321  0.055742762013807
322  0.045492744488986
323  0.036935220196777
324  0.029840852217657
325  0.023997740489123
326  0.019214335904674
327  0.015320580887648
328  0.012167779438633
329  0.009627606069476
330  0.007590575511228
331  0.005964217696074
332  0.004671136982604
333  0.003647081075751
334  0.002839102788446
335  0.002203865389816
336  0.001706118264261
337  0.001317352436400
338  0.001014633776781
339  0.000779604140824
340  0.000597636251764
341  0.000457125954290
342  0.000348904854607
343  0.000265756797153
344  0.000202022712087
345  0.000153279812313
346  0.000116082720651
347  0.000087755728843
348  0.000066226941806
349  0.000049896490408
350  0.000037532277385
351  ];
352
353
354%TR=2.039;
355%onset=[30 75 120]';
356%durat=[15 15 15]';
357onset = onset(:);
358durat = durat(:);
359
360%bf = spm_hrf(1/16);
361
362% Create stimulus function (32 bin offset)
363%k = 165; % Number of scans
364k = num_scans;
365T     = 16;
366dt    = TR/T;
367u     = ones(size(onset));
368if ~any(durat)
369  u  = u/dt;
370end
371ton = round(onset*TR/dt) + 32;                  % onsets
372tof = round(durat*TR/dt) + ton + 1;                     % offset
373sf = zeros((k*T + 128),size(u,2));
374
375for j = 1:length(ton)
376  if numel(sf)>ton(j),
377    sf(ton(j),:) = sf(ton(j),:) + u(j,:);
378  end;
379  if numel(sf)>tof(j),
380    sf(tof(j),:) = sf(tof(j),:) - u(j,:);
381  end;
382end
383sf        = cumsum(sf);                                 % integrate
384
385% Convolve stimulus with the HRF
386conv_sf = conv(sf,bf);
387
388% Resample the convolved stimulus
389R = conv_sf([0:(k-1)]*16+1+32);
390R=R(:);
391
392% Remove linear trend but keep the mean
393mean_ts_data = mean(ts_data);
394ts_data = detrend(ts_data)+mean_ts_data;
395
396% Estimate baseline from the paradigm
397bl_ind = find(sf==0);
398bl_ind = bl_ind(1:min(length(bl_ind),20));
399if isempty(bl_ind)
400  baseline = mean(ts_data);
401else
402  baseline = mean(ts_data(bl_ind));
403end
404
405mean_ts_data_bold = (ts_data./baseline-1)*100;
406
407
408% Fit block model to the data
409%z=mean_ts_data_bold;
410z=ts_data;
411H = [R ones(size(R))];
412th = H\z;
413z_hat = H*th;
414
415% Calculate results
416bold_prc = ((th(1)*max(R))/th(2))*100;
417
Note: See TracBrowser for help on using the repository browser.

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