Changeset 82 for plugins/fmri_plugins/basic_fmri_analysis.m
 Timestamp:
 Jun 16, 2009, 10:54:45 AM (10 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

plugins/fmri_plugins/basic_fmri_analysis.m
r78 r82 5 5 6 6 if AddInfo.isDataMixed 7 return 7 % Make data a 4Dmatrix 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; 8 16 end 9 17 … … 53 61 'hipass',hipass,'omitvolumes',omitVols,... 54 62 'mask',ROI(roi_ind).voxels{1}(:,:,:,1)); 63 ROI(roi_ind)=[]; 55 64 else 56 65 [maps_out,fdr_th] = fmri_analysis(DATA{1},'TR',TR,... … … 93 102 DATA{1}.FTDATA(:,:,:,omitVols)=[]; 94 103 overlay.FTDATA(:,:,:,omitVols)=[]; 104 else 105 new_onset = onset; 106 new_durat = durat; 107 end 108 109 110 % Estimate BOLD strength for ROIs 111 if ~isempty(ROI) 112 nScans = size(DATA{1}.FTDATA,4); 113 fprintf(1,'\n'); 114 for kk=1:length(ROI) 115 116 % Get timeseries 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'); 152 else 153 % Estimate BOLD strength for the timeseries corresponding to the 154 % maximum value in Tmap. 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 Tvalue'); 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 Tvalue'); 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'); 95 202 end 96 203 … … 104 211 105 212 213 function [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!!! 217 bf = [ 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]'; 287 onset = onset(:); 288 durat = durat(:); 289 290 %bf = spm_hrf(1/16); 291 292 % Create stimulus function (32 bin offset) 293 %k = 165; % Number of scans 294 k = num_scans; 295 T = 16; 296 dt = TR/T; 297 u = ones(size(onset)); 298 if ~any(durat) 299 u = u/dt; 300 end 301 ton = round(onset*TR/dt) + 32; % onsets 302 tof = round(durat*TR/dt) + ton + 1; % offset 303 sf = zeros((k*T + 128),size(u,2)); 304 305 for 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; 312 end 313 sf = cumsum(sf); % integrate 314 315 % Convolve stimulus with the HRF 316 conv_sf = conv(sf,bf); 317 318 % Resample the convolved stimulus 319 R = conv_sf([0:(k1)]*16+1+32); 320 R=R(:); 321 322 % Remove linear trend but keep the mean 323 mean_ts_data = mean(ts_data); 324 ts_data = detrend(ts_data)+mean_ts_data; 325 326 % Estimate baseline from the paradigm 327 bl_ind = find(sf==0); 328 bl_ind = bl_ind(1:min(length(bl_ind),20)); 329 if isempty(bl_ind) 330 baseline = mean(ts_data); 331 else 332 baseline = mean(ts_data(bl_ind)); 333 end 334 335 mean_ts_data_bold = (ts_data./baseline1)*100; 336 337 338 % Fit block model to the data 339 %z=mean_ts_data_bold; 340 z=ts_data; 341 H = [R ones(size(R))]; 342 th = H\z; 343 z_hat = H*th; 344 345 % Calculate results 346 bold_prc = ((th(1)*max(R))/th(2))*100; 347
Note: See TracChangeset
for help on using the changeset viewer.