Changeset 82 for plugins


Ignore:
Timestamp:
Jun 16, 2009, 10:54:45 AM (10 years ago)
Author:
tjniskan
Message:
  • 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:
1 edited

Legend:

Unmodified
Added
Removed
  • plugins/fmri_plugins/basic_fmri_analysis.m

    r78 r82  
    55
    66if AddInfo.isDataMixed
    7   return
     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;
    816end
    917
     
    5361    'hipass',hipass,'omitvolumes',omitVols,...
    5462    'mask',ROI(roi_ind).voxels{1}(:,:,:,1));
     63  ROI(roi_ind)=[];
    5564else
    5665  [maps_out,fdr_th] = fmri_analysis(DATA{1},'TR',TR,...
     
    93102  DATA{1}.FTDATA(:,:,:,omitVols)=[];
    94103  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');
    95202end
    96203
     
    104211
    105212
     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 TracChangeset for help on using the changeset viewer.

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