fmri_gen_st_datamat.m

绝对经典,老外制作的功能强大的matlab实现PLS_TOOBOX

function fmri_gen_st_datamat(sessionFile,num_skipped_scans,win_size, ...
		run_idx,coords_info,merge_across_runs_flg,behavdata,behavname)
%
% USAGE: 
%    fmri_gen_st_datamat(sessionFile,num_skipped_scans,win_size,run_idx,
%					coord_file, merge_across_runs_flg)
%  or
%    fmri_gen_st_datamat(sessionFile,num_skipped_scans,win_size,run_idx,
%					coord_struct, merge_across_runs_flg)
%
%   Generate the spatial/temporal datamat from the PLS/datamat_run? file
%   using PLS/combined_coords for the coords.
%   
%   Input: 
%      sessionFile:  the name of file that contains the session information.
%      num_skipped_scans:  (optional) [default: num_skipped_scans = 0]
%      win_size:  (optional) [default: win_size = 8]
%                 number of scans for the temporal windows. 
%      run_idx:   (optional) [default: run_idx contains all the runs]
%		  The indices of runs that will be used to generate st_datamat.
%      coord_file: a file contains the common coords and the indices of the
%                  each individual coords for the runs to make the common coords
%      coord_struct:  a structure contains the same information as the above
%                     coord_file
%         note: if neither coord_file or coord_struct is specified, the coords
%               information is read from "{datamat_prefix}_common_coords.mat"
%               in the PLS data directory.
%
%      merge_across_runs_flg:  1 - merge the same condition across all runs
%      (optional)                  to become one row of the ST datamat
%                              0 - merge the same condition only within a run
%                  	       [default: merge_across_runs_flg = 1] 
%
%   Input File:
%      ${pls_data_path}/${datamat_prefix}_run?.mat: 
%                 the datamat file for the runs.  The files are
%                 created by the "pls_gen_datamat" script.
%
%   Output File:
%      ${pls_data_path}/${datamat_prefix}_fMRIdatamat.mat:  
%           st_datamat: the combined spatial/temporal datamat 
%           st_evt_list:   the type of event associated with each row 
%           
%
%   EXAMPLE:
%       fmri_gen_st_datamat('PLSsession',0,8,[1:3],coords_info,1);
%

   progress_hdl = ShowProgress('initialize');

   load(sessionFile);		% load the session information

   pls_data_path = session_info.pls_data_path;
   datamat_prefix = session_info.datamat_prefix;

   % handle session file without condition baseline specification
   %
   if isfield(session_info,'condition_baseline')  
      cond_baseline = session_info.condition_baseline;
   else
      cond_baseline = [];
   end;


   num_runs = session_info.num_runs;
   if ~exist('run_idx','var') | isempty(run_idx),
      run_idx = [1:num_runs];
   end;

   if ~exist('num_skipped_scans','var') | isempty(num_skipped_scans),
      num_skipped_scans = 0;
   end;

   if ~exist('win_size','var') | isempty(win_size),
      win_size = 8;
   end;

   if ~exist('merge_across_runs_flg','var') | isempty(merge_across_runs_flg),
      merge_across_runs_flg = 1;
   end;

   if ~exist('coords_info','var') | isempty(coords_info),
       coords_file = sprintf('%s_common_coords.mat',datamat_prefix);
       coords_file = fullfile(pls_data_path,coords_file);
       load(coords_file);
   elseif ischar(coords_info(1))		% coords file name is specified
       coords_file = coords_info;
       load(coords_file);
   end;

   if (isfield(coords_info,'coords_idx') == 1),
       all_runs_have_the_same_coords = 0;
   else
       all_runs_have_the_same_coords = 1;
       num_cols = length(coords_info.common) * win_size;
   end;

   st_datamat = [];
   st_evt_list = [];
   st_evt_cnt = [];
   num_runs = length(run_idx);
   for i=1:num_runs,
      run_start_row = size(st_datamat,1)+1;

      msg = sprintf('Loading datamat of run#%d',run_idx(i));
      ShowProgress(progress_hdl,msg)

      run_info = session_info.run(run_idx(i));

      % load the datamat, coords and dims for each run
      %
      datamat_file = sprintf('%s_run%d.mat',datamat_prefix,run_idx(i));

      load(fullfile(pls_data_path,datamat_file));

      if (all_runs_have_the_same_coords == 0),
          c = coords_info.coords_idx{i};	% use only the common voxels
          datamat = datamat(:,c);
          num_cols = length(c) * win_size;
      else
          num_cols = length(coords_info.common) * win_size;
      end;

      % get rid of the skipped scans
      %
      datamat = datamat(num_skipped_scans+1:end,:);
      num_scans = size(datamat,1);

      num_conds = length(run_info.evt_onsets);

      for j=1:num_conds,		% for each condition

          ShowProgress(progress_hdl,(i-1)+j/num_conds);
          msg = sprintf('Working on condition %d of run#%d... ',j,run_idx(i));
          ShowProgress(progress_hdl,msg)

          %  extract the datamat scans that matches the current condition
          %
          row_idx = 0;
          num_onsets = length(run_info.evt_onsets{j});
          evt_datamat = zeros(num_onsets,num_cols);

          for k=1:num_onsets;
             start_scan = floor(run_info.evt_onsets{j}(k)-num_skipped_scans)+1;
             end_scan = start_scan+win_size-1;

	     if isempty(cond_baseline)
                baseline_start = start_scan;
                baseline_scans = start_scan;
                baseline_length = 1;
	     else
	        baseline_info = cond_baseline{j};
                baseline_start = start_scan + baseline_info(1);
                baseline_scans = [baseline_start:baseline_start+baseline_info(2)-1];
	     end;

             if (start_scan <= 0) | (num_scans < end_scan)	% out of bound
                 disp(sprintf('Scans %03d for the condition %d of run%d are not included due to out of bound',floor(run_info.evt_onsets{j}(k)),j,run_idx(i)));
             elseif (baseline_start <= 0)
                 disp(sprintf('Scans %03d for the condition %d of run%d are not included due to out of bound baseline',floor(run_info.evt_onsets{j}(k)),j,run_idx(i)));
             else
		 % convert to the percent signal change
		 %
		 baseline_signals = mean(datamat(baseline_scans,:),1);
		 baseline_signals = baseline_signals(ones(1,win_size),:);
                 tmp_datamat = datamat(start_scan:end_scan,:);
		 tmp_datamat = (tmp_datamat - baseline_signals ) * 100 ./ baseline_signals;
                 row_idx = row_idx+1;  
                 evt_datamat(row_idx,:) = reshape(tmp_datamat,[1 num_cols]);
             end;
          end;

          if (merge_across_runs_flg == 0 & row_idx == 0),
             errmsg = sprintf('No onsets for condition #%d in the run #%d',j,i);
             errordlg(['ERROR: ' errmsg],'Generating ST Datamat Error');
             waitfor(gcf);
             return;
          end;

          if (row_idx ~= 0),	  % some stimuli for the current condition

             if (merge_across_runs_flg == 0)	% merge data within run only
                st_evt_list = [st_evt_list j];
                st_datamat = [st_datamat; mean(evt_datamat(1:row_idx,:),1)];
             else
                if isempty(st_evt_list) 
                   idx = [];
                else
                   idx = find(st_evt_list == j);
                end;

                if isempty(idx),			% new event
                   st_evt_list = [st_evt_list j];
	           st_evt_cnt = [st_evt_cnt row_idx];
                   st_datamat = [st_datamat; sum(evt_datamat(1:row_idx,:),1)];
                else
	           st_evt_cnt(idx) = st_evt_cnt(idx) + row_idx;
                   st_datamat(idx,:) = st_datamat(idx,:) + ...
                                           sum(evt_datamat(1:row_idx,:),1);
                end;

             end; 
          end;   % for (row_idx~=0)

      end;   % for each in condition list

      %  save some memory usage
      %
      clear evt_datamat 

   end;	 % for each run

   ShowProgress(progress_hdl,(num_runs+1));
   if (merge_across_runs_flg == 1)	% merge across all runs 
      for j=1:length(st_evt_list)
         st_datamat(j,:) = st_datamat(j,:) / st_evt_cnt(j);
      end;
   end;

   [new_evt_list, order_idx] = reorder_evt_list(st_evt_list,num_conds);
   if isempty(new_evt_list),
     errmsg = sprintf('ERROR: Some conditions have no trials at all.\n');
     errordlg(errmsg,'Creating ST Datamat Error');
     waitfor(gcf);
     return;
   end;

   st_datamat = st_datamat(order_idx,:);
   st_evt_list = new_evt_list;


   %  everything is done, ready to save the information
   %
   st_dims = dims;
   st_voxel_size = voxel_size;
   st_origin = origin;
   st_coords = coords_info.common;
   st_win_size = win_size;

   st_sessionFile = sessionFile;

   if (st_sessionFile ~= filesep),			
      curr = pwd;
      if isempty(curr)
         curr = filesep;
      end

      st_sessionFile = fullfile(curr,st_sessionFile);
   end;

   fname = sprintf('%s_fMRIdatamat.mat',datamat_prefix);
   st_datafile = fullfile(pls_data_path,fname);

   ShowProgress(progress_hdl,(num_runs+2));
   ShowProgress(progress_hdl,['Saving ST Datamat into the file: ' fname]);

   if(exist(st_datafile,'file')==2)  % datamat file with same filename exist
      dlg_title = 'Confirm File Overwrite';
      msg = ['File ',fname,' exist. Are you sure you want to overwrite it?'];
      response = questdlg(msg,dlg_title,'Yes','No','Yes');

      if(strcmp(response,'No'))
         cd(pls_data_path);
         putfile_filter = [datamat_prefix,'*_fMRIdatamat.mat'];
         [filename, pathname] = uiputfile(putfile_filter,'Save As');
         if isequal(filename,0)
            close(progress_hdl);
            msg1 = ['WARNING: No file is saved.'];
%            uiwait(msgbox(msg1,'Uncomplete','modal'));
            return;
         else
            st_datafile = fullfile(pathname, filename);
            session_info.datamat_prefix = strrep(filename,'_fMRIdatamat.mat','');
            save(sessionFile,'session_info');
         end
      end
   end

   try
      save(st_datafile,'st_datamat','st_coords','st_dims','st_voxel_size', ...
                'st_origin','st_evt_list', 'st_win_size','st_sessionFile', ...
		'normalize_volume_mean','behavdata','behavname');
   catch
      errmsg = sprintf('ERROR: Cannot save ST datamat to the file: \n   %s.',st_datafile);
      errordlg(errmsg,'Creating ST Datamat Error');
      waitfor(gcf);
      return;
   end;
   
   return;


%-------------------------------------------------------------------------
%
function detrended_datamat = detrend(datamat);

  num_scans = size(datamat,1);
  linear_term = [1:num_scans]';
  linear_term = linear_term - mean(linear_term);
  quad_term = linear_term .* linear_term;
  
  detrend_terms = [ones(num_scans,1) linear_term quad_term]; 

  beta = detrend_terms \ datamat;
  pred = detrend_terms * beta;

  detrended_datamat = datamat - pred;

  %% add the mean back to the detrended datamat for computing the
  %% percent signal change in later stage.   [by Barry Giesbrecht]
  
  voxel_mean = mean(datamat,1);
  detrended_datamat = detrended_datamat + voxel_mean(ones(1,num_scans),:);

  return;					% detrend


%-------------------------------------------------------------------------
function [new_evt_list,reorder_idx] = reorder_evt_list(evt_list,num_conditions)
%  make sure the row order of the st_datamat is repeat order of
%  conditions, i.e cond#1, cond#2,... cond#1, cond#2, .... etc
%

  num_rep = length(evt_list) / num_conditions;

  [new_evt_list, reorder_idx] = sort(evt_list);
  if ~isequal(reshape(new_evt_list,[num_rep,num_conditions]), ...
             repmat([1:num_conditions],num_rep,1))
     new_evt_list = [];
     reorder_idx = [];
  end;

  return;					% reorder evt_list


%-------------------------------------------------------------------------
function hdl = ShowProgress(progress_hdl,info)

  %  'initialize' - return progress handle if any
  %
  if ischar(progress_hdl) & strcmp(lower(progress_hdl),'initialize'),
     if ~isempty(gcf) & isequal(get(gcf,'Tag'),'ProgressFigure'),
         hdl = gcf;
     else
         hdl = [];
     end;
     return;
  end;


  if ~isempty(progress_hdl)
     if ischar(info)
         rri_progress_status(progress_hdl,'Show_message',info);
     else
         rri_progress_status(progress_hdl,'Update_bar',info);
     end;
     return;
  end;

  if ischar(info),
     disp(info)
  end;

  return;					% ShowProgress