fmri_get2_datamat.m

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

%FMRI_GET_DATAMAT will save the run data into 'prefix_run%d' file.
%
%	The run data includes:
%
%		coords
%		dims
%		voxel_size
%		origin
%		sessionFile
%		normalize_volume_mean
%
%	which can also be obtained from FMRI_GEN_DATAMAT,
%	and the followings:
%
%		tmp_datamat
%		baseline_signals
%
%	both are cell arrays for each run, which is indexed by num_cond
%	and num_onsets.
%
%	Usage: fmri_get_datamat(sessionFile, run_num, brain_mask, ...
%		sd_thresh, ignore_slices, normalize_volume_mean, ...
%		num_skipped_scans, win_size)
%
%
%   Input: 
%       sessionFile:  the name of file that contains the session information.
%       brain_mask:  a matrix, with same size as the IMG files, defines the
%                   brain region.  Any element larger than 0 is treated as
%                   part of the brain region.
%       coord_thresh: threshold to generate compute the brain region 
%                    [default: 0.15];
%       sd_thresh: maximum standard deviation (s.d.) for the brain voxels.  Any
%       (optional) voxels with s.d. greater than this value will be removed to
%                 avoid the inclusion of venous and sinus. [default: 2];
%       ignored_slices: a vector specifies the number of slices which should
%       (optional) be skipped to generate the datamat. [default: []]
%       num_skipped_scans:  (optional) [default: num_skipped_scans = 0]
%       win_size:  (optional) [default: win_size = 8]
%                 number of scans for the temporal windows. 
%
%
%   Output Files:
%      {prefix_datamat}_run?.mat:  
%                     (created in the directory specified by dir_path in 
%                      the sessionFile)
%               tmp_datamat: spatial/temporal datamat
%		baseline_signals: average of baseline scans
%               coords: coordinates for the datamat
%               dims: dimension of the IMG file.
%
%   EXAMPLE:
%     fmri_gen_datamat('PLSsession',3);      % generate datamat for run#3
%
%
%	See also FMRI_GEN_DATAMAT, FMRI_GEN_ST_DATAMAT

%	called by FMRI_CREATE_DATAMAT_UI
%
%   Modified on 12-May-2003 by Jimmy Shen
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function fmri_get_datamat(sessionFile,run_num,varargin)

   progress_hdl = ShowProgress('initialize');

   brain_mask = [];
   coord_thresh = 0.15;
   sd_thresh = 2;
   ignore_slices = [];
   normalize_volume_mean = 0;
   num_skipped_scans = 0;
   win_size = 8;

   mask_dims = [];

   if (nargin == 2),
      use_brain_mask = 0;
   elseif (nargin >= 2),
      if  (ischar(varargin{1})) ,	% use predefined brain region
          use_brain_mask = 1;
          brain_mask_file = varargin{1};
          brain_mask = rri_load_img(brain_mask_file);	% use brain mask file
          mask_dims = size(brain_mask);
      else
          use_brain_mask = 0;		% compute brain region from the data
          coord_thresh = varargin{1};
      end;
   
      if (nargin >= 4)
         sd_thresh = varargin{2};
      end;
   
      if (nargin >= 5)
         ignore_slices = varargin{3};
      end;

      if (nargin >= 6)
         normalize_volume_mean = varargin{4};
      end;

      if (nargin >= 7)
         num_skipped_scans = varargin{5};
      end;

      if (nargin >= 8)
         win_size = varargin{6};
      end;
   end

   load(sessionFile);		% load the session information

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

   %  save run data of each run, into datamat_file
   %
   fname = sprintf('%s_run%d.mat',datamat_prefix,run_num);
   datamat_file = fullfile(pls_data_path,fname);

   %  generate the datamat for each run
   %
   if (sessionFile ~= filesep)		
      curr = pwd;
      if isempty(curr)
         curr = filesep;
      end

      sessionFile = fullfile(curr,sessionFile);
   end;

   ShowProgress(progress_hdl,sprintf('Creating datamat for run#%d ... ',run_num));

   %  get stuff from session_info first
   %
   cond_baseline = session_info.condition_baseline;
   num_conds = session_info.num_conditions;
   run_info = session_info.run(run_num);

   % get maximum number of onsets within this run
   %
   max_onsets = find_max_onsets(run_info);

   %  get image files path/file name
   %
   data_path = run_info.data_path;
   flist = run_info.data_files;
   num_files = length(flist);

   img_file = fullfile(data_path,flist{1});
   [dims,voxel_size,origin] = rri_imginfo(img_file);
   dims = [dims(1) dims(2) 1 dims(3)];

   if (use_brain_mask==1) & ~isequal(dims,mask_dims),
       errmsg ='ERROR: Dimensions of the data do not match that of the brain mask!';
       errordlg(errmsg,'Brain Mask Error');
       waitfor(gcf);
       return;
   end;

   %  skipped scans handling
   %
   if ~exist('num_skipped_scans','var') | isempty(num_skipped_scans)
      num_skipped_scans = 0;
   end;

   %  ignored slices handling
   %
   if ~exist('ignore_slices','var') | isempty(ignore_slices)
      s_idx = [1:dims(4)];		% handle all slices
   else
      idx = zeros(1,dims(4));
      idx(ignore_slices) = 1;
      s_idx = find(s_idx == 0);
   end;

   num_voxels = prod(dims);

   %  initial tmp_datamat and baseline_signals
   %
   tmp_datamat = cell(num_conds, max_onsets);		% will be saved
   baseline_signals = cell(num_conds, max_onsets);	% will be saved
   datamat = [];					% temp datamat

   % create step for progress bar
   %
   progress_step = 1 / (num_conds * max_onsets + 3);

   for j = 1:num_conds
      num_onsets = length(run_info.evt_onsets{j});
      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_end = baseline_start;
         else
            baseline_info = cond_baseline{j};
            baseline_start = start_scan + baseline_info(1);
            baseline_end = baseline_start+baseline_info(2)-1;
         end;

         if (start_scan <= 0) | (num_files < 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

            %  read in image files
            %
            dataset = zeros(length([start_scan:end_scan]), num_voxels);

            for scan = start_scan:end_scan
               img_file = [data_path filesep flist{scan}];
               img = rri_load_img(img_file);
               img = img(:,:,:,s_idx);
               dataset(scan - start_scan + 1,:) = img(:)';
            end;

            tmp_datamat{j,k} = dataset;
            datamat = [datamat; tmp_datamat{j,k}];

            %  read in baseline_signals image files
            %
            dataset = zeros(length([baseline_start:baseline_end]), num_voxels);

            for scan = baseline_start:baseline_end
               img_file = [data_path filesep flist{scan}];
               img = rri_load_img(img_file);
               img = img(:,:,:,s_idx);
               dataset(scan - baseline_start + 1,:) = img(:)';
            end;

            baseline_signals{j,k} = mean(dataset,1);
            datamat = [datamat; baseline_signals{j,k}];

            ShowProgress(progress_hdl, ((j-1)*num_onsets+k)*progress_step);

         end;
      end;	% for num_onsets
   end;		% for num_conds

   %    determine the coords of the brain region
   %
   if (use_brain_mask == 0)		% compute coords from (temp) datamat
      coords = rri_make_coords(datamat, 1/coord_thresh);
      datamat = datamat(:,coords);

      %  get rid of voxels that have standard deviation "sd_thresh" times larger 
      %  than the grand average standard deviation.
      %
      std_d = std(datamat);
      thresh = mean(std_d) * sd_thresh;
      idx = find(std_d < thresh);
      coords = coords(idx);
   else
      coords = find( brain_mask(:) > 0)';
      if ~isempty(ignore_slices)	% skip the slices if any
         m = zeros(dims);		
         m(coords) = 1;
         m(:,:,:,s_idx);
         coords = find(m == 1)';
      end;
   end;

   %    apply the coords of the brain region
   %
   for j = 1:num_conds
      num_onsets = length(run_info.evt_onsets{j});
      for k = 1:num_onsets
         tmp_datamat{j,k} = tmp_datamat{j,k}(:,coords);
         baseline_signals{j,k} = baseline_signals{j,k}(:,coords);

         if (normalize_volume_mean == 1),
            num_voxels = length(coords);

            mean_tmp_datamat = mean(tmp_datamat{j,k},2);
            tmp_datamat{j,k} = tmp_datamat{j,k} ./ ...
		mean_tmp_datamat(:,ones(1,num_voxels));

            mean_baseline_signals = mean(baseline_signals{j,k},2);
            baseline_signals{j,k} = baseline_signals{j,k} ./ ...
		mean_baseline_signals(:,ones(1,num_voxels));
         end;
      end;
   end;

   ShowProgress(progress_hdl, (num_conds * max_onsets + 1)*progress_step);

   %  save the datamat, which is a matrix contains only the brain voxels.
   %
   ShowProgress(progress_hdl,['Saving datamat into the file: ' fname ]);

   try 
      save(datamat_file,'tmp_datamat','baseline_signals','coords', ...
	'dims','voxel_size','origin','sessionFile','normalize_volume_mean');
   catch
      errmsg = sprintf('ERROR: Cannot save datamat to file: \n   %s.', ...
                       datamat_file);
      errordlg(errmsg,'Save Datamat Error');
      waitfor(gcf);
      return;
   end

   return;						% fmri_get_datamat


%-------------------------------------------------------------------------

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


%----------------------------------------------------------------------------
%
%  get the maximum number of onsets within a run
%
%----------------------------------------------------------------------------

function max_onsets = find_max_onsets(run_info)

   num_cond = length(run_info.evt_onsets);
   max_onsets = 0;

   for i = 1:num_cond
      if length(run_info.evt_onsets{i}) > max_onsets
         max_onsets = length(run_info.evt_onsets{i});
      end;
   end;

   return;						% find_max_onsets