fmri_apply_std.m

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

function fmri_apply_std(isfmri, prefix, plspath)

   if isfmri
      datamatfile = [prefix '_fMRIdatamat.mat'];
   else
      datamatfile = [prefix '_BfMRIdatamat.mat'];
   end;

   datamatfile = fullfile(plspath, datamatfile);

   try
      load(datamatfile, 'st_sessionFile', 'st_coords', 'st_dims', 'st_win_size');
   catch
      msg = ['Datamat file "' datamatfile];
      msg = [msg '" does not exist, please Create Datamat first.'];
      uiwait(msgbox(msg,'Error'));
      return;
   end

   try
      load(st_sessionFile);
   catch
      msg = ['Can not open session file "' st_sessionFile];
      msg = [msg '".'];
      uiwait(msgbox(msg,'Error'));
      return;
   end

   num_runs = session_info.num_runs;
   num_voxels = prod(st_dims);

   %  estimate and check memory
   %
   for run_num = 1:num_runs

      run_info = session_info.run(run_num);
      data_path = run_info.data_path;
      flist = run_info.data_files;
      num_files = length(flist);

      mem = round(num_files * num_voxels * 8 / 1000000);	% 8 came from double data type

      if mem > 32
         msg = ['After expanding image files in Run ', num2str(run_num), ', you need more than '];
         msg = [msg, num2str(mem), 'MB of memory to process it. Do you want to continue?'];
         ans = questdlg(msg,'Warning','No','Yes','No');

         if strcmp(ans,'No')
            return;
         end
      end

   end;

   msg1 = 'Enter maximum standard deviation allowed in the brain voxels:';
   def1 = '2';
   msg2 = 'Enter number of scans to be skipped:';
   def2 = '0';
   msg3 = 'Enter runs to be included:';
   def3 = num2str([1:num_runs]);
   msg4 = 'Enter slices to be ignored:';
   def4 = '';

   tit = 'Apply Max STD value';

   ans = inputdlg({msg1, msg2, msg3, msg4}, tit, 1, {def1, def2, def3, def4});

   if isempty(ans)
      msg = 'value of maximum standard deviation must be inputed';
      uiwait(msgbox(msg,'Error'));
      return;
   end

   sd_thresh = str2num(ans{1});
   num_skipped_scans = str2num(ans{2});
   run_idx = str2num(ans{3});
   ignore_slices = str2num(ans{4});

   if isempty(ignore_slices)
      s_idx = [1:st_dims(4)];			% handle all slices
   else
      s_idx = zeros(1,st_dims(4));
      s_idx(ignore_slices) = 1;
      s_idx = find(s_idx == 0);
   end;

   num_runs = length(run_idx);
   coords = zeros(1, num_voxels);
   coords(st_coords) = 1;
   ischange = 0;
   run_cnt = 0;

   progress_hdl = rri_progress_ui('initialize', 'Loading data ...');

   for run_num = run_idx

      run_cnt = run_cnt + 1;

      msg = ['Loading images in Run ', num2str(run_num)];
      rri_progress_ui(progress_hdl, 'Loading data ...', msg);

      run_info = session_info.run(run_num);
      data_path = run_info.data_path;
      flist = run_info.data_files;
      total_num_files = length(flist);

      used_files = zeros(1, total_num_files);

      if isfmri
         for cond_num = 1:length(run_info.evt_onsets)
            onset = run_info.evt_onsets{cond_num};

            for onset_num = 1:length(onset)
               for win_num = 1:st_win_size
                  used_files(onset(onset_num)+win_num) = 1;
               end
            end
         end
      else
         for cond_num = 1:length(run_info.blk_onsets)
            onset = run_info.blk_onsets{cond_num};
            width = run_info.blk_length{cond_num};

            for onset_num = 1:length(onset)
               for win_num = 1:width(onset_num)
                  used_files(onset(onset_num)+win_num) = 1;
               end
            end
         end
      end

      used_files = find(used_files);
      num_files = length(used_files);

      dataset = zeros(num_files,num_voxels);

      for j=used_files

         progress_step = (run_cnt-1)/num_runs + (j-1)/(total_num_files*num_runs);
         rri_progress_ui(progress_hdl, 'Loading data ...', progress_step);

         if length(flist)>1
            img_file = fullfile(data_path, flist{j-num_skipped_scans});
            img = load_nii(img_file);
         else
            img_file = fullfile(data_path, flist{1});
            img = load_nii(img_file, j-num_skipped_scans);
         end

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

      end;

      datamat = dataset(:,st_coords);

      %  get rid of voxels that have standard deviation "sd_thresh" times larger 
      %  than the grand average standard deviation.
      %
      std_d = std(datamat);
      std_d_max = max(std_d);
      std_d_mean = mean(std_d);

      thresh = std_d_mean * sd_thresh;

      if std_d_max > thresh
         idx = find(std_d > thresh);
         idx = st_coords(idx);			% non brain voxels
         coords(idx) = 0;

         msg0 = ['. The st_coords is changed, and st_datamat is modified'];
         ischange = 1;
      else
         msg0 = ['. There is no change based on Run ', num2str(run_num)];
      end

      msg = ['Maximum STD is ', num2str(std_d_max), ', Average STD is ', num2str(std_d_mean), msg0];
      disp(msg);

   end;

   clear dataset datamat std_d;

   if ischange
      coords = find(coords);

      load(datamatfile);
      [r c] = size(st_datamat);
      st_datamat = reshape(st_datamat, [st_win_size, r, c/st_win_size]);

      [tmp idx] = intersect(st_coords, coords);
      st_datamat = st_datamat(:,:,idx);
      st_datamat = reshape(st_datamat, [r, st_win_size*length(coords)]);
      st_coords = coords;

      try
         save(datamatfile, 'behavdata', 'behavname', 'create_ver', ...
		'normalize_volume_mean', 'st_coords', 'st_datamat', ...
		'st_dims', 'st_evt_list', 'st_origin', 'st_sessionFile', ...
		'st_voxel_size', 'st_win_size');
      catch
         msg = ['Can not save file "', datamatfile, '".'];
         uiwait(msgbox(msg,'Error'));
      end
   end

   close(progress_hdl);

   return;