Changeset 94 for misclib


Ignore:
Timestamp:
Oct 12, 2009, 1:20:14 PM (10 years ago)
Author:
tjniskan
Message:
  • Fixed a bug in reading Varian FID-files (regression bug)
  • Added a possibility of Using 3D overlays on 4D data in Aedes
  • misclib/fmri_analysis.m now accepts general regressors

M misclib/fmri_analysis.m
M aedes.m
M aedes_readfid.m
M aedes_revision.m

File:
1 edited

Legend:

Unmodified
Added
Removed
  • misclib/fmri_analysis.m

    r83 r94  
    1 function [maps_out,fdrth] = fmri_analysis(DATA,varargin)
     1function [maps_out,fdrth,H] = fmri_analysis(DATA,varargin)
    22%
    33% This function does a quick fMRI analysis in SPM style. No motion
     
    1717show_wbar = true;   % Show waitbar by default
    1818mask = [];          % Don't mask by default
     19regr = [];          % Don't use additional regressors by default
    1920
    2021if rem(length(varargin),2)~=0
     
    5758        error('Invalid mask!')
    5859      end
    59   end
    60 end
    61 
    62 if length(durat)==1
    63   durat = repmat(durat,1,length(onset));
    64   durat = durat(:);
    65 elseif length(durat)~=length(onset)
    66   error('Mismatch in the number of elements in onset and duration!')
     60    case 'regressor'
     61      regr = value;
     62  end
     63end
     64
     65if ~isempty(onset)
     66  if length(durat)==1
     67    durat = repmat(durat,1,length(onset));
     68    durat = durat(:);
     69  elseif length(durat)~=length(onset)
     70    error('Mismatch in the number of elements in onset and duration!')
     71  end
    6772end
    6873
     
    8388if isempty(mask)
    8489  mask = true(size(data,1),size(data,2),size(data,3));
     90end
     91
     92% Check regressors
     93if ~isempty(regr)
     94  if size(regr,1)~=size(data,4) && size(regr,2)~=size(data,4)
     95    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
     100  regr = detrend(regr);
    85101end
    86102
     
    164180if ~isempty(omitVols)
    165181 
    166   % Calculate new onsets and durations
    167   ton = onset;
    168   tof = onset+durat+1;
    169   tmp=zeros(1,size(data,4));
    170   tmp(ton)=1;tmp(tof)=-1;
    171   sf=cumsum(tmp);
    172   sf(omitVols)=[];
    173   tmp=diff([0 sf]);
    174   onset = find(tmp==1);
    175   durat = find(tmp==-1);
    176   if isempty(durat)
    177     % Block stays up...
    178     durat=length(tmp)-onset(end)-1;
    179   elseif length(durat) < length(onset)
    180     durat(1:end-1) = durat(1:end-1)-onset(1:end-1)-1;
    181     durat(end) = length(tmp)-onset(end)-1;
    182   else
    183     durat = durat-onset-1;
    184   end
    185   onset=onset(:);
    186   durat=durat(:);
     182  fprintf(1,'Skipping requested volumes...\n');
     183  if ~isempty(onset)
     184    % Calculate new onsets and durations
     185    ton = onset;
     186    tof = onset+durat+1;
     187    tmp=zeros(1,size(data,4));
     188    tmp(ton)=1;tmp(tof)=-1;
     189    sf=cumsum(tmp);
     190    sf(omitVols)=[];
     191    tmp=diff([0 sf]);
     192    onset = find(tmp==1);
     193    durat = find(tmp==-1);
     194    if isempty(durat)
     195      % Block stays up...
     196      durat=length(tmp)-onset(end)-1;
     197    elseif length(durat) < length(onset)
     198      durat(1:end-1) = durat(1:end-1)-onset(1:end-1)-1;
     199      durat(end) = length(tmp)-onset(end)-1;
     200    else
     201      durat = durat-onset-1;
     202    end
     203    onset=onset(:);
     204    durat=durat(:);
     205   
     206    fprintf(1,['New onsets: ',num2str(onset(:)'),'\n']);
     207    fprintf(1,['New duration: ',num2str(durat(:)'),'\n']);
     208  end
    187209 
    188   fprintf(1,'Skipping requested volumes...\n');
    189   fprintf(1,['New onsets: ',num2str(onset(:)'),'\n']);
    190   fprintf(1,['New duration: ',num2str(durat(:)'),'\n']);
     210  if ~isempty(regr)
     211    regr(omitVols,:)=[];
     212  end
    191213  data(:,:,:,omitVols) = [];
     214   
    192215end
    193216
     
    196219T     = 16;
    197220dt    = TR/T;
    198 u     = onset.^0;
    199 if ~any(durat)
    200   u  = u/dt;
    201 end
    202 ton       = round(onset*TR/dt) + 32;                    % onsets
    203 tof       = round(durat*TR/dt) + ton + 1;                       % offset
    204 sf        = zeros((k*T + 128),size(u,2));
    205 ton       = max(ton,1);
    206 tof       = max(tof,1);
    207 for j = 1:length(ton)
    208   if numel(sf)>ton(j),
    209     sf(ton(j),:) = sf(ton(j),:) + u(j,:);
    210   end;
    211   if numel(sf)>tof(j),
    212     sf(tof(j),:) = sf(tof(j),:) - u(j,:);
    213   end;
    214 end
    215 sf        = cumsum(sf);                                 % integrate
    216 
    217 % Convolve stimulus with the HRF
    218 conv_sf = conv(sf,bf);
    219 
    220 % Resample the convolved stimulus
    221 R = conv_sf([0:(k-1)]*16+1+32);
    222 R=R(:);
     221if ~isempty(onset)
     222  u     = onset.^0;
     223  if ~any(durat)
     224    u  = u/dt;
     225  end
     226  ton       = round(onset*TR/dt) + 32;                  % onsets
     227  tof       = round(durat*TR/dt) + ton + 1;                     % offset
     228  sf        = zeros((k*T + 128),size(u,2));
     229  ton       = max(ton,1);
     230  tof       = max(tof,1);
     231  for j = 1:length(ton)
     232    if numel(sf)>ton(j),
     233      sf(ton(j),:) = sf(ton(j),:) + u(j,:);
     234    end;
     235    if numel(sf)>tof(j),
     236      sf(tof(j),:) = sf(tof(j),:) - u(j,:);
     237    end;
     238  end
     239  sf        = cumsum(sf);                                       % integrate
     240 
     241  % Convolve stimulus with the HRF
     242  conv_sf = conv(sf,bf);
     243 
     244  % Resample the convolved stimulus
     245  R = conv_sf([0:(k-1)]*16+1+32);
     246  R=R(:);
     247else
     248  R=[];
     249end
     250
    223251
    224252% Smooth data
     
    243271  S = 0;
    244272end
    245 H = [R ones(size(R))];
     273H = [R regr ones(k,1)];
    246274H = H-S*(S'*H); % Filter design matrix
    247275
     
    333361end
    334362
    335 if nargout==2
     363if nargout>=2
    336364  fdrth = thresFDR;
    337365end
Note: See TracChangeset for help on using the changeset viewer.

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