Changeset 94


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

Files:
4 edited

Legend:

Unmodified
Added
Removed
  • aedes.m

    r80 r94  
    37073707end
    37083708
     3709% If the data underlying the 3D overlay data is 4D, use only the first
     3710% volume of the overlay...
     3711if size(Dat.ImageOverlay,4)==1
     3712  volInd = 1;
     3713else
     3714  volInd = Dat.CurrentVol;
     3715end
    37093716
    37103717% Extract slice data
     
    37523759  if ~Dat.isOverlayRGB
    37533760    % Convert indexed image to RGB image
    3754     slice1_ind=Dat.ImageOverlay(:,:,Dat.Slices(1),Dat.CurrentVol);
     3761    slice1_ind=Dat.ImageOverlay(:,:,Dat.Slices(1),volInd);
    37553762    slice1_ind=double(slice1_ind);
    37563763       
     
    37853792if updateY
    37863793  if ~Dat.isOverlayRGB
    3787     slice2_ind=reshape(Dat.ImageOverlay(:,Dat.Slices(2),:,Dat.CurrentVol),...
     3794    slice2_ind=reshape(Dat.ImageOverlay(:,Dat.Slices(2),:,volInd),...
    37883795          Dat.ImageDim(Dat.DataInd,1), ...
    37893796          Dat.ImageDim(Dat.DataInd,3));
     
    38203827if updateZ
    38213828  if ~Dat.isOverlayRGB
    3822     slice3_ind=reshape(Dat.ImageOverlay(Dat.Slices(3),:,:,Dat.CurrentVol),...
     3829    slice3_ind=reshape(Dat.ImageOverlay(Dat.Slices(3),:,:,volInd),...
    38233830          Dat.ImageDim(Dat.DataInd,2),Dat.ImageDim(Dat.DataInd,3));
    38243831    slice3_ind=double(slice3_ind);
     
    41454152else
    41464153  if ~Dat.isOverlayRGB && ~ismember(sz1,sz2,'rows')
    4147         clear tmp_data;
    4148         h=errordlg({'Data/overlay size mismatch!',...
    4149           '',['Data size: ',num2str(sz2(1)),'x',num2str(sz2(2)),'x',num2str(sz2(3))],...
    4150           ['Overlay size: ',num2str(sz1(1)),'x',num2str(sz1(2)),'x',num2str(sz1(3))]},...
    4151           'Error loading overlay!','modal');
    4152         return
     4154    % Allow Overlay loading if only 4th dimension differs
     4155    if ~ismember(sz1(1:3),sz2(1:3),'rows')
     4156      clear tmp_data;
     4157      h=errordlg({'Data/overlay size mismatch!',...
     4158        '',['Data size: ',num2str(sz2(1)),'x',num2str(sz2(2)),'x',num2str(sz2(3))],...
     4159        ['Overlay size: ',num2str(sz1(1)),'x',num2str(sz1(2)),'x',num2str(sz1(3))]},...
     4160        'Error loading overlay!','modal');
     4161      return
     4162    end
    41534163  elseif Dat.isOverlayRGB && ~ismember(sz1(1:3),sz2(1:3),'rows')
    4154         clear tmp_data;
    4155         h=errordlg({'Data/overlay size mismatch!',...
    4156           '',['Data size: ',num2str(sz2(1)),'x',num2str(sz2(2)),'x',num2str(sz2(3))],...
    4157           ['Overlay size: ',num2str(sz1(1)),'x',num2str(sz1(2)),'x',num2str(sz1(3))]},...
    4158           'Error loading overlay!','modal');
    4159         return
     4164    clear tmp_data;
     4165    h=errordlg({'Data/overlay size mismatch!',...
     4166      '',['Data size: ',num2str(sz2(1)),'x',num2str(sz2(2)),'x',num2str(sz2(3))],...
     4167      ['Overlay size: ',num2str(sz1(1)),'x',num2str(sz1(2)),'x',num2str(sz1(3))]},...
     4168      'Error loading overlay!','modal');
     4169    return
    41604170  end
    41614171end
  • aedes_readfid.m

    r93 r94  
    14511451    elseif length(size(kspace))==3
    14521452      kspace = reshape(kspace,procpar.np/2,procpar.nv,[]);
    1453       if Dat.Sorting
     1453      if Dat.Sorting && ~isempty(Dat.phasetable)
    14541454        kspace(:,Dat.phasetable(:),:) = kspace;
    14551455      end
  • aedes_revision.m

    r93 r94  
    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

    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