Changeset 164


Ignore:
Timestamp:
May 23, 2011, 9:18:21 AM (8 years ago)
Author:
tjniskan
Message:
  • Added option for using external regressors in misclib/fmri_analysis.m

M misclib/fmri_smooth.m
M misclib/fmri_analysis.m
M aedes_revision.m

Files:
3 edited

Legend:

Unmodified
Added
Removed
  • aedes_revision.m

    r163 r164  
    9393% bash-script every time it is called so that this file "aedes_revision.m" is
    9494% always in the list of committed files. DO NOT EDIT THE NEXT LINE!!!
    95 % - SVN Hook -
     95% - Svn Hook -
  • misclib/fmri_analysis.m

    r119 r164  
    66
    77% Defaults
    8 TR = 2.039;
    9 contr = [1 0];
    10 onset = [30 75 120].';
    11 durat = [15 15 15].';
    12 smooth_kernel = [2 2 1];
     8TR = [];
     9contr = [1];
     10onset = [];
     11durat = [];
     12smooth_kernel = [];
    1313d_cutoff = [];      % Don't temporally filter data by default
    1414qFDR = 0.05;        % Default threshold for FDR correction
     
    3333      onset = value;
    3434      onset = onset(:);
    35     case {'durat','durations'}
     35    case {'durat','durations','duration'}
    3636      durat = value;
    3737      durat = durat(:);
     
    6363end
    6464
     65% Check that we have something to fit...
     66if isempty(onset) && isempty(regr)
     67        warning('No onsets or regressors defined. Only mean will be fitted!')
     68end
     69
     70% Check that TR has been given
     71if isempty(TR)
     72        error('Repetition time (TR) has not been defined!')
     73end
     74
     75% Check onsets and durations
    6576if ~isempty(onset)
     77        if isempty(durat)
     78                error('Durations for onsets have not been defined!')
     79        end
    6680  if length(durat)==1
    6781    durat = repmat(durat,1,length(onset));
     
    7286end
    7387
    74 
    7588% Check data
    7689if isstruct(DATA) && isfield(DATA,'FTDATA')
     
    7992  data=DATA;
    8093else
    81   error('Input data has to be 3D or 4D numeric matrix!')
    82 end
    83 
     94  error('Input data has to be a 3D/4D numeric matrix or Aedes data structure!')
     95end
     96
     97% Permute 3D matrix to 4D
    8498if ndims(data)==3
    8599  data = permute(data,[1 2 4 3]);
    86100end
    87101
     102% Initialize mask
    88103if isempty(mask)
    89104  mask = true(size(data,1),size(data,2),size(data,3));
     
    92107% Check regressors
    93108if ~isempty(regr)
    94   if size(regr,1)~=size(data,4) && size(regr,2)~=size(data,4)
     109  if size(regr,1)~=size(data,4)
    95110    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
     111        end
    100112  regr = detrend(regr);
    101113end
     
    243255 
    244256  % Resample the convolved stimulus
    245   R = conv_sf([0:(k-1)]*16+1+32);
    246   R=R(:);
    247 else
    248   R=[];
    249 end
    250 
    251 
    252 % Smooth data
    253 if show_wbar
    254   wbh = aedes_calc_wait('Smoothing data...');
    255   drawnow
    256 end
    257 if all(smooth_kernel)
    258   tic
    259   smooth_data = fmri_smooth(data,smooth_kernel);
    260   toc
    261 else
    262   fprintf(1,'fMRI analysis warning: Could not smooth data!\n');
    263   smooth_data = data;
    264 end
     257  RR = conv_sf([0:(k-1)]*16+1+32);
     258  RR=RR(:);
     259else
     260  RR=[];
     261end
     262
     263
     264% Smooth data if requested
     265if ~isempty(smooth_kernel) && ~ismember(smooth_kernel,[1 1 1],'rows')
     266        if show_wbar
     267                wbh = aedes_calc_wait('Smoothing data...');
     268                drawnow
     269        end
     270        if all(smooth_kernel)
     271                tic
     272                smooth_data = fmri_smooth(data,smooth_kernel);
     273                toc
     274        else
     275                fprintf(1,'fMRI analysis warning: Could not smooth data!\n');
     276                smooth_data = data;
     277        end
     278        if show_wbar
     279                close(wbh);
     280        end
     281else
     282        smooth_data = data;
     283end
     284
    265285
    266286if ~isempty(d_cutoff)
     
    273293  S = 0;
    274294end
    275 H = [R regr ones(k,1)];
     295H = [RR regr ones(k,1)];
    276296H = H-S*(S'*H); % Filter design matrix
     297nParam = size(H,2);
     298maps_out = struct('pmap',[],'tmap',[]);
     299maps_out.pmap = zeros(size(smooth_data,1),size(smooth_data,2),size(smooth_data,3),nParam);
     300maps_out.tmap = zeros(size(smooth_data,1),size(smooth_data,2),size(smooth_data,3));
    277301
    278302% Calculate parametric map(s)
     
    280304nCols = size(smooth_data,2);
    281305nRows = size(smooth_data,1);
    282 pmap1 = zeros(size(smooth_data,1),size(smooth_data,2),size(smooth_data,3));
    283 pmap2 = zeros(size(pmap1));
    284 tmap = zeros(size(pmap1));
     306if length(contr)<nParam
     307        contr(nParam)=0; % Pad with zeros
     308end
    285309c = repmat(contr,nRows,1);
    286310HTH = pinv(H'*H);
    287311R = eye(size(H,1))-H*HTH*H';
    288 if show_wbar
    289   close(wbh);
    290 end
     312
    291313if show_wbar
    292314  wbh = aedes_wbar(0,sprintf('Estimating parameters. Processing plane 0/%d',nPlanes));
    293315  drawnow
    294316end
     317
     318% Process data in columns
    295319for ii=1:nPlanes
    296320 
     
    303327    col_data = col_data-S*(S'*col_data);
    304328    th=H\col_data;
    305     pmap1(:,kk,ii)=th(1,:).';
    306     pmap2(:,kk,ii)=th(2,:).';
    307    
    308329    r = col_data-H*th;
    309330    rr=diag(r'*r);
    310331    rr=rr(:);
    311332    sig2 = rr./trace(R*KKT);
    312     sig2 = repmat(sig2,1,2);
     333    sig2 = repmat(sig2,1,nParam);
    313334    T = diag(c*th)./sqrt(c.*sig2*HTH*H'*KKT*H*HTH*contr');
    314335    T(find(mask(:,kk,ii)==0))=0;
    315     tmap(:,kk,ii)=T.';
    316   end
     336                maps_out.pmap(:,kk,ii,:) = th.';
     337                maps_out.tmap(:,kk,ii)=T.';
     338        end
    317339end
    318340if show_wbar
     
    329351
    330352% p-values
    331 pval_map = 1-tdist(tmap,dof);
     353pval_map = 1-tdist(maps_out.tmap,dof);
    332354
    333355% Perform FDR (False Discovery Rate) correction
     
    335357
    336358pValues = pval_map(:);
    337 tValues = tmap(:);
     359tValues = maps_out.tmap(:);
    338360
    339361[pValuesSorted,sortInd] = sort(pValues);
     
    343365pFDR = [1:nP]'/nP*qFDR/cV; % FDR-correction
    344366thresFDRind = find(pValuesSorted<=pFDR,1,'last');
    345 %if isempty(thresFDRind)
    346 %       fprintf(1,'Warning: No significant values with the FDR-corrected threshold %f\n',...
    347 %               qFDR);
    348 %       return
    349 %end
    350367if ~isempty(thresFDRind)
    351368  thresFDR = tValuesSorted(thresFDRind);
     
    354371end
    355372
    356 maps_out.pmap1 = pmap1;
    357 maps_out.pmap2 = pmap2;
    358 %maps_out.pmap3 = pmap3;
    359 maps_out.tmap = tmap;
    360 
    361373if show_wbar
    362374  close(wbh);
  • misclib/fmri_smooth.m

    r119 r164  
    3434% Calculate standard deviations using FWHM
    3535if fwhmInPixels
    36    stds = fwhm_sz/sqrt(8*log(2));
     36        stds = fwhm_sz/sqrt(8*log(2));
    3737else
    3838  stds = (fwhm_sz/sqrt(8*log(2)))./voxsize;
     
    4444  kernel_sz(3)=floor((size(data,3)-1)/2);
    4545end
    46 kernel_sz
    4746
    4847% Construct the smoothing kernel
     
    6867  tmp_sz = size(data(:,:,:,1));
    6968  %tmp_sz(3)=tmp_sz(3)+fwhm_sz(3);
    70   ind_hi=ceil(tmp_sz/2)+floor(size(s_kernel)/2);
    71   ind_lo=ind_hi-(size(s_kernel)-1);
    72   ind_lo(ind_lo<1)=1;
    73   tmp_kernel = zeros(tmp_sz);
    74   tmp_kernel((ind_lo(1):ind_hi(1))+1,(ind_lo(2):ind_hi(2))+1,...
    75     (ind_lo(3):ind_hi(3)))=s_kernel;
     69  %ind_hi=ceil(tmp_sz/2)+floor(size(s_kernel)/2);
     70  %ind_lo=ind_hi-(size(s_kernel)-1);
     71  %ind_lo(ind_lo<1)=1;
     72  %tmp_kernel = zeros(tmp_sz);
     73  %tmp_kernel((ind_lo(1):ind_hi(1))+1,(ind_lo(2):ind_hi(2))+1,...
     74  %  (ind_lo(3):ind_hi(3)))=s_kernel;
     75        tmp_kernel = s_kernel;
     76        tmp_kernel(tmp_sz(1),tmp_sz(2),tmp_sz(3))=0; % Pad with zeros
    7677end
    7778for ii=1:nVols
     
    9293    %tmp_data(pad_sz(1),pad_sz(2),pad_sz(3))=0;
    9394    %tmp_data = cat(3,tmp_data(:,:,1),tmp_data,tmp_data(:,:,end));
    94     tmp_smooth = fftshift(ifftn(fftn(tmp_data).*fftn(tmp_kernel)));
    95     smooth_data(:,:,:,ii) = tmp_smooth(:,:,1:size(data,3));
     95    tmp_smooth = fftshift(ifftn(fftn(tmp_data).*fftn(tmp_kernel)),3);
     96    smooth_data(:,:,:,ii) = tmp_smooth;%tmp_smooth(:,:,1:size(data,3));
    9697  end
    9798end
Note: See TracChangeset for help on using the changeset viewer.

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