fmri_gen_datamat.m

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

function fmri_gen_datamat(sessionFile,run_num,varargin)

%
% USAGE: fmri_gen_datamat(sessionFile,brain_mask,sd_thresh,ignored_slices)
%    or  fmri_gen_datamat(sessionFile,coord_thresh,sd_thresh,ignore_slices)
%        
%   Generate the datamat based on the information in the sessionFile.  
%
%   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: 7];
%      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: []]
%
%
%   Output Files:
%      {prefix_datamat}_run?.mat:  
%                     (created in the directory specified by dir_path in 
%                      the sessionFile)
%               datamat: spatial/temporal datamat
%               coords: coordinates for the datamat
%               dims: dimension of the IMG file.
%
%   EXAMPLE:
%     fmri_gen_datamat('PLSsession',3);      % generate datamat for run#3
%

   progress_hdl = ShowProgress('initialize');

   brain_mask = [];
   coord_thresh = 0.15;
   sd_thresh = 2;
   ignore_slices = [];
   normalize_volume_mean = 0;
   
   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;
   end
   
   
   load(sessionFile);		% load the session information
   
   pls_data_path = session_info.pls_data_path;
   datamat_prefix = session_info.datamat_prefix;

   %  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));

   run_info = session_info.run(run_num);
   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;

   if 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;   

   %    create a matrix to store all image data
   %
   progress_step = 1 / (num_files + 3);

   num_voxels = prod(dims);
   dataset = zeros(num_files,num_voxels);
   for j=1:num_files,
      img_file = [data_path filesep flist{j}];
      img = rri_load_img(img_file);
      img = img(:,:,:,s_idx);
      dataset(j,:) = img(:)';
      ShowProgress(progress_hdl,j*progress_step);
   end;

   %    determine the coords of the brain region
   %
   if (use_brain_mask == 0)		% compute coords from the dataset
      coords = rri_make_coords(dataset,1/coord_thresh);
      datamat = dataset(:,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);
      datamat = datamat(:,idx);
      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;
      datamat = dataset(:,coords);

   end;

   if (normalize_volume_mean == 1),
      num_voxels = length(coords);
      mean_datamat = mean(datamat,2);
      datamat = datamat ./ mean_datamat(:,ones(1,num_voxels));
   end;

   ShowProgress(progress_hdl,(num_files+1)*progress_step);

   %  save the datamat, which is a matrix contains only the brain voxels.
   %
   fname = sprintf('%s_run%d.mat',datamat_prefix,run_num);
   datamat_file = fullfile(pls_data_path,fname);

   ShowProgress(progress_hdl,['Saving datamat into the file: ' fname ]);
   try 
      save(datamat_file,'datamat','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;


%-------------------------------------------------------------------------
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