fmri_get_datamat.m

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

%FMRI_GET_DATAMAT will create the Event Related or Blocked fMRI datamat
%
%	Usage: fmri_get_datamat(sessionFile, run_idx, brain_mask, ...
%		sd_thresh, ignore_slices, normalize_volume_mean, ...
%		num_skipped_scans, win_size, merge_across_runs_flg, ...
%		behavdata, behavname)
%
%   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];
%       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.
%       run_idx:   (optional) [default: run_idx contains all the runs]
%		  The indices of runs that will be used to generate st_datamat.
%       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]
%
%	See also FMRI_GEN_DATAMAT, FMRI_GEN_ST_DATAMAT

%
%	the following paramater is removed:
%
%       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];
%

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

function fmri_get_datamat(sessionFile, run_idx, varargin)

   progress_hdl = rri_progress_ui('initialize');
   brain_mask_file = '';
   brain_mask = [];
   coord_thresh = [];
   sd_thresh = 2;
   ignore_slices = [];
   normalize_volume_mean = 0;
   num_skipped_scans = 0;
   win_size = 8;
   for_batch = 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 = load_nii(brain_mask_file, 1);
          mask_dims = [brain_mask.hdr.dime.dim(2:3) 1 brain_mask.hdr.dime.dim(4)];
          brain_mask = reshape(int8(brain_mask.img), [brain_mask.hdr.dime.dim(2:3) 1 brain_mask.hdr.dime.dim(4)]);
      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;

      if (nargin >= 9)
         merge_across_runs_flg = varargin{7};
      end;

      if (nargin >= 10)
         behavdata = varargin{8};
      end;

      if (nargin >= 11)
         behavname = varargin{9};
      end;

      if (nargin >= 12)
         session_win_hdl = varargin{10};
      end;

      if (nargin >= 13)
         normalize_signal_mean = varargin{11};
      end;

      if (nargin >= 14)
         ConsiderAllVoxels = varargin{12};
      end;

      if (nargin >= 15)
         SingleSubject = varargin{13};
      end;

      if (nargin >= 16)
         orient = varargin{14};
      end;

      if (nargin >= 17)
         for_batch = varargin{15};
      end;

   end

   create_datamat_info.brain_mask_file = brain_mask_file;
   create_datamat_info.brain_coord_thresh = coord_thresh;
   create_datamat_info.consider_all_voxels_as_brain = ConsiderAllVoxels;
   create_datamat_info.num_skipped_scans = num_skipped_scans;
   create_datamat_info.run_idx = run_idx;
   create_datamat_info.ignore_slices = ignore_slices;
   create_datamat_info.temporal_window_size = win_size;
   create_datamat_info.normalize_volume_mean = normalize_volume_mean;
   create_datamat_info.normalize_with_baseline = normalize_signal_mean;
   create_datamat_info.merge_across_runs = merge_across_runs_flg;
   create_datamat_info.single_subject_analysis = SingleSubject;

   load(sessionFile);		% load the session information

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

   %  save temp data into datamat_file
   %
%   fname = sprintf('%s_run.mat',datamat_prefix);
%   datamat_file = fullfile(pls_data_path,fname);

      curr = pwd;
      if isempty(curr)
         curr = filesep;
      end

   if (sessionFile ~= filesep)		

      sessionFile = fullfile(curr,sessionFile);
   end;

   %  get stuff from session_info first
   %
   cond_baseline = session_info.condition_baseline0;
   num_conds = session_info.num_conditions0;

   % get maximum number of onsets
   %
   max_onsets = 0;

   for i = run_idx
      run_info = session_info.run(i);
      tmp = find_max_onsets(run_info);

      if tmp > max_onsets
         max_onsets = tmp;
      end
   end

   %  get image files path/file name <for read only one image file>
   %
   data_path = run_info.data_path;
   flist = run_info.data_files;

   img_file = fullfile(data_path,flist{1});

   if ~isempty(orient) & ~isempty(orient.orient_pattern)
       dims = orient.dims;
       voxel_size = orient.voxel_size;
       origin = orient.origin;
   else
      [dims,voxel_size,origin] = rri_imginfo(img_file);
   end

   dims = [dims(1) dims(2) 1 dims(3)];

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

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

   %  default is merge_across_runs
   %
   if ~exist('merge_across_runs_flg','var') | isempty(merge_across_runs_flg),
      merge_across_runs_flg = 1;
   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;

   if (use_brain_mask == 1)			% coords from brain_mask
      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;

   num_voxels = prod(dims);
   num_cols = num_voxels * win_size;

   if (use_brain_mask == 0)
      coords = zeros(1, num_voxels);		% initial st_coords
   end

   datamat = [];				% st datamat before coords
   st_datamat = [];				% st datamat
   st_evt_list = [];
   st_evt_cnt = [];

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

   for i = run_idx

      run_info = session_info.run(i);

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

      if num_files==1
         img_file = fullfile(data_path,flist{1});
         num_files = get_nii_frame(img_file);
      end

      for j = 1:num_conds

         num_onsets = length(run_info.evt_onsets{j});

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

         if SingleSubject
            evt_datamat = [];
         else
            evt_datamat = zeros(1,num_cols);
         end

         rri_progress_ui(progress_hdl, '', ...
            sprintf('Creating datamat for run#%d condition#%d ... ',i,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
                  if length(flist)>1
                     img_file = fullfile(data_path, flist{scan});
                     img = load_nii(img_file);
                  else
                     img_file = fullfile(data_path, flist{1});
                     img = load_nii(img_file, scan);
                  end

                  v7 = version;
                  if str2num(v7(1))<7
                     img = reshape(double(img.img), [img.hdr.dime.dim(2:3) 1 img.hdr.dime.dim(4)]);
                  else
                     img = reshape(single(img.img), [img.hdr.dime.dim(2:3) 1 img.hdr.dime.dim(4)]);
                  end

                  if ~isempty(orient) & ~isempty(orient.orient_pattern)
                     img = img(orient.orient_pattern);
                     img = reshape(img,dims);
                  end

                  img = img(:,:,:,s_idx);
                  dataset(scan - start_scan + 1,:) = img(:)';
               end;

               tmp_datamat = dataset;

               %  find brain voxel coords for each onset, and
               %  accumulated to find common coords for all onsets
               %
               if (use_brain_mask == 0)
                  coords = coords + find_onset_coords(dataset,coord_thresh,ConsiderAllVoxels);
               end

               if (normalize_volume_mean == 1)
                  if (use_brain_mask == 0)
                     mean_dataset = dataset(:,find(coords == 0));
                  else
                     mean_dataset = dataset(:,coords);
                  end

                  mean_dataset = mean(mean_dataset, 2);
                  mean_dataset = mean_dataset(:,ones(1,num_voxels));
                  tmp_datamat = dataset ./ mean_dataset;
               end

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

               for scan = baseline_start:baseline_end
                  if length(flist)>1
                     img_file = fullfile(data_path, flist{scan});
                     img = load_nii(img_file);
                  else
                     img_file = fullfile(data_path, flist{1});
                     img = load_nii(img_file, scan);
                  end

                  img = reshape(double(img.img), [img.hdr.dime.dim(2:3) 1 img.hdr.dime.dim(4)]);

                  if ~isempty(orient) & ~isempty(orient.orient_pattern)
                     img = img(orient.orient_pattern);
                     img = reshape(img,dims);
                  end

                  img = img(:,:,:,s_idx);
                  dataset(scan - baseline_start + 1,:) = img(:)';
               end;

               %  this is not a duplicate of normalize_volume_mean,
               %  it is doing for baseline dataset
               %
               if (normalize_volume_mean == 1)
                  if (use_brain_mask == 0)
                     mean_dataset = dataset(:,find(coords == 0));
                  else
                     mean_dataset = dataset(:,coords);
                  end

                  mean_dataset = mean(mean_dataset, 2);
                  mean_dataset = mean_dataset(:,ones(1,num_voxels));
                  dataset = dataset ./ mean_dataset;
               end

               baseline_signals = mean(dataset,1);

               %  repmat for win_size of rows
               %
               baseline_signals = baseline_signals(ones(1,win_size),:);

               %  if 0, make it less significant, so it can be removed later
               %
		zero_baseline_idx = find(baseline_signals==0);
               baseline_signals(find(baseline_signals==0)) = 99999;

               if normalize_signal_mean
                  tmp_datamat = ...
                     (tmp_datamat - baseline_signals ) * 100 ./ baseline_signals;
               end

		tmp_datamat(zero_baseline_idx) = 0; %99999

               row_idx = row_idx + 1;

               if SingleSubject
                  evt_datamat = [evt_datamat; single(reshape(tmp_datamat,[1 num_cols]))];
               else
                  evt_datamat = single(double(evt_datamat) + reshape(tmp_datamat,[1 num_cols]));
               end

               stp = (i-1)*num_conds*max_onsets + (j-1)*max_onsets + k;
               rri_progress_ui(progress_hdl, '', stp * progress_step);

            end;
         end;				% for num_onsets

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

         %  Generate the spatial/temporal datamat from evt_datamat
         %  depending on 'across run' or 'within run'
         %
         if (row_idx ~= 0),		% some stimuli for the current condition

            if (merge_across_runs_flg == 0)	% merge data within run only
               if SingleSubject
%                  st_evt_list = [st_evt_list j*ones(1, size(evt_datamat,1))];
                  st_evt_list = [st_evt_list ((i-1)*num_conds+j)*ones(1, size(evt_datamat,1))];
                  datamat = [datamat; evt_datamat];
               else
%                  st_evt_list = [st_evt_list j];
                  st_evt_list = [st_evt_list (i-1)*num_conds+j];
                  datamat = [datamat; single(double(evt_datamat) / num_onsets)];
               end
            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];

                  if SingleSubject
                     datamat = [datamat; evt_datamat(:)'];