fmri_plot_rf.m

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

           fname = sprintf('%s_fMRIdatamat.mat',session_info.datamat_prefix);
        end

        st_filename{cnt}.name = fname;
        st_filename{cnt}.fullname = fullfile(session_info.pls_data_path,fname);
        st_filename{cnt}.group = i;
        st_filename{cnt}.profile = sessionFileList{i}{j};
     end;
  end;

  setappdata(gcf,'STFiles',st_filename);

  return;                                       % get_st_datamat_filename


%--------------------------------------------------------------------------
%
function  [st_datamat, org_coords,st_win_size,st_evt_list] =  ...
                                         load_plotted_datamat,

   old_pointer = get(gcf,'Pointer');
   set(gcf,'Pointer','watch');

   selected_files = get_selected_filename;

   if (length(selected_files) == 1), 		
       load(selected_files{1});

       common_coords = getappdata(gcf,'STCoords');
       [tmp, coord_idx] = intersect(st_coords, common_coords);

       nr = length(st_evt_list);
       nc = length(st_coords);
       win_size = st_win_size;
       num_voxels = length(st_coords);
       num_cols = win_size * num_voxels;
       tmp_datamat = reshape(st_datamat,[nr,win_size,nc]); 
       tmp_datamat = tmp_datamat(:,:,coord_idx);
       num_voxels = length(common_coords);
       num_cols = win_size * num_voxels;
       st_datamat = reshape(tmp_datamat,[nr,num_cols]);
       st_coords = common_coords;
   else					%  merge files together
       [st_datamat,st_coords,st_win_size,st_evt_list] =  ...
                                   merge_st_datamat(selected_files);
   end;

   org_coords = st_coords;

   cond_selection = getappdata(gcf,'cond_selection');

   [mask, st_evt_list, evt_length] = ...
	fmri_mask_evt_list(st_evt_list, cond_selection);

   st_datamat = st_datamat(mask,:);

   set(gcf,'Pointer',old_pointer);

   return;                                       % load_plotted_datamat


%--------------------------------------------------------------------------
%
function  [selected_files] = get_selected_filename(select_all_flg),

   st_filename = getappdata(gcf,'STFiles');

   if exist('select_all_flg','var') & (select_all_flg == 1)
      data_option = 3;
   else
      h = findobj(gcf,'Tag','STDatamatPopup');
      selected_idx = get(h,'Value');
      data_option = get(h,'Userdata');
   end;

   switch (data_option)
      case {1}, 				% individual file
	 selected_files{1} = st_filename{selected_idx}.fullname;
      case {2},					% group data
         cnt = 0;
         for i=1:length(st_filename),
            if (st_filename{i}.group == selected_idx)
               cnt = cnt+1;
               selected_files{cnt} = st_filename{i}.fullname;
            end;
         end;
      case {3},					% all data
         cnt = 0;
         for i=1:length(st_filename),
            selected_files{i} = st_filename{i}.fullname;
         end;
   end; 

   return;					%  get selected filenames


%--------------------------------------------------------------------------
%
function  [grp_st_datamat,coords,win_size,evt_list] =  ...
                                   merge_st_datamat(selected_files),

   num_files = length(selected_files);

   %  compute the common coords first and the total number of events
   %  (don't need to compute the common coords, because it has been
   %   computed during run PLS)
   %
   load(selected_files{1},'st_dims','st_win_size'); 
   total_num_evt = 0;
%   m = zeros(st_dims);
   for i=1:num_files,
      load(selected_files{i},'st_coords','st_evt_list');
      total_num_evt = total_num_evt + length(st_evt_list); 
%      m(st_coords) = m(st_coords) + 1; 
   end;
%   coords = find(m == num_files); 
   coords = getappdata(gcf,'STCoords');

   %  ready to merge the st_datamat together now ...
   %
   win_size = st_win_size;
   num_voxels = length(coords);
   num_cols = win_size * num_voxels;
   grp_st_datamat = zeros(total_num_evt,num_cols);
   evt_list = [];
   first_row = 1;
   for i=1:num_files,

      load(selected_files{i});

%      coord_idx = find(m(st_coords) == num_files);
      [tmp, coord_idx] = intersect(st_coords, coords);

      nr = length(st_evt_list);
      nc = length(st_coords);
      last_row = nr + first_row - 1;

      num_voxels = length(st_coords);
      num_cols = win_size * num_voxels;

      tmp_datamat = reshape(st_datamat,[nr,win_size,nc]); 
      tmp_datamat = tmp_datamat(:,:,coord_idx);
      num_voxels = length(coords);
      num_cols = win_size * num_voxels;
      tmp_datamat = reshape(tmp_datamat,[nr,num_cols]);

      grp_st_datamat(first_row:last_row,:) = tmp_datamat;
      evt_list = [evt_list st_evt_list];

      first_row = last_row + 1;

      clear st_datamat tmp_datamat;
   end; 

   return;                               	% merge_st_datamat


%--------------------------------------------------------------------------
%
function  plot_response_fn()

   neighbor_numbers = 1;

   emp_st_data = 0;
   set(findobj(gcf,'tag','Stable'),'Userdata',0,'Label','&Stable off');

   st_files = getappdata(gcf,'STFiles');
   conditions = getappdata(gcf,'Conditions');

   h = findobj(gcf,'Tag','PlotButton');
   if strcmp(lower(get(h,'Enable')),'off'),  return; end;

   h = findobj(gcf,'Tag','STDatamatPopup');
   popup_string = get(h,'String');
   selected_idx = get(h,'Value');
   selected_data = popup_string{selected_idx};

   %  load the datamat if not loaded yet
   %
   last_datamat = getappdata(gcf,'PlottedDatamat');
   if strcmp(last_datamat,selected_data) == 0 

       set(findobj(gcf,'Tag','MessageLine'),'String','Loading data ... ');

       [st_datamat, org_coords, st_win_size, st_evt_list] = ...
						load_plotted_datamat;

       setappdata(gcf,'PlottedDatamat',selected_data);
       setappdata(gcf,'STDatamat',st_datamat);
       setappdata(gcf,'org_coords',org_coords);
       setappdata(gcf,'STWinSize',st_win_size);
       setappdata(gcf,'STEvtList',st_evt_list);

       set(findobj(gcf,'Tag','MessageLine'),'String','');
   else
       st_win_size = getappdata(gcf,'STWinSize');
       st_evt_list = getappdata(gcf,'STEvtList');
       org_coords = getappdata(gcf,'org_coords');
   end;

   st_coords = getappdata(gcf,'STCoords');	% common_coord from result file

   num_pts_plotted = getappdata(gcf,'NumPtsPlotted');
   if isempty(num_pts_plotted)
      num_pts_plotted = st_win_size;
      setappdata(gcf,'NumPtsPlotted',num_pts_plotted);
   elseif (num_pts_plotted > st_win_size),
      num_pts_plotted = st_win_size;
   end;

   xyz = getappdata(gcf,'XYZ');
%   lag = getappdata(gcf,'lag');
   
   %  extract the currect ploting data
   %
   cur_coord = getappdata(gcf,'Coord');
   if isempty(cur_coord)
       msg = 'ERROR: No point has been seleted to plot.';
       set(findobj(gcf,'Tag','MessageLine'),'String',msg);
%       return;
       emp_st_data = 1;
   end;

   coord_idx = find(st_coords == cur_coord);	% idx of cur_coord in common

   if isempty(coord_idx)
       msg = 'ERROR: The selected point is outside the brain region.';
       set(findobj(gcf,'Tag','MessageLine'),'String',msg);
%       return;
       emp_st_data = 1;
   end;


   %  Get Neighborhood Size
   h = findobj(gcf,'Tag','neighborhoodEdit');
   neighbor_size = round(str2num(get(h,'String')));

   if isempty(neighbor_size) | ~isnumeric(neighbor_size)
      neighbor_size = 0;
   end

   %  Do neighborhood mean only if there is any neighborhood
   %
   if neighbor_size > 0

      CallBackFig = getappdata(gcf,'CallBackFig');
      CurrLVIdx = getappdata(CallBackFig,'CurrLVIdx');
      dims = getappdata(CallBackFig,'STDims');

      %  Get neighborhood XYZs
      %
      xyz = rri_coord2xyz(cur_coord, dims);

      x1 = xyz(1) - neighbor_size;
      if x1 < 1, x1 = 1; end;

      x2 = xyz(1) + neighbor_size;
      if x2 > dims(1), x2 = dims(1); end;

      y1 = xyz(2) - neighbor_size;
      if y1 < 1, y1 = 1; end;

      y2 = xyz(2) + neighbor_size;
      if y2 > dims(2), y2 = dims(2); end;

      z1 = xyz(3) - neighbor_size;
      if z1 < 1, z1 = 1; end;

      z2 = xyz(3) + neighbor_size;
      if z2 > dims(4), z2 = dims(4); end;

      %  Get neighborhood coords relative to whole volume
      %
      neighbor_coord = [];

      for k = z1:z2
         for j=y1:y2
            for i=x1:x2
               neighbor_coord = [neighbor_coord, rri_xyz2coord([i j k], dims)];
            end
         end
      end

      %  If "Cluster Mask" is checked, cluster masked voxels will be used
      %  as a criteria to select surrounding voxels
      %
      isbsr = getappdata(CallBackFig,'ViewBootstrapRatio');

      if isbsr
         cluster_info = getappdata(CallBackFig, 'cluster_bsr');
      else
         cluster_info = getappdata(CallBackFig, 'cluster_blv');
      end

      %  Get cluster voxels coords relative to whole volume
      %
      if length(cluster_info) < CurrLVIdx
         cluster_info = [];
      else
         cluster_info = cluster_info{CurrLVIdx};
         cluster_info = cluster_info.data{1}.idx;
      end

      %  If "Bootstrap" is computed, voxels that meet the bootstrap ratio
      %  threshold will be used as a criteria to select surrounding voxels
      %
      BSThreshold = getappdata(CallBackFig,'BSThreshold');

      if ~isempty(BSThreshold)
         BSRatio = getappdata(CallBackFig,'BSRatio');
         BSRatio = BSRatio(:,CurrLVIdx);

         all_voxel = zeros(1,prod(dims));
         bsr_gt_thresh = (abs(BSRatio) >= BSThreshold);
         bsr_gt_thresh = reshape(bsr_gt_thresh, [st_win_size, round(length(bsr_gt_thresh)/st_win_size)]);
%         all_voxel(st_coords) = bsr_gt_thresh(lag+1,:);
         all_voxel(st_coords) = sum(bsr_gt_thresh, 1);

         bsridx = find(all_voxel);
      else
         bsridx = [];
      end

      %  Only including surrounding voxels that meet the bootstrap ratio
      %  threshold, or are part of cluster masked voxels
      %
      bsr_cluster_coords = unique([cluster_info bsridx cur_coord]);

      %  Intersect of neighborhood coord "neighbor_coord" and 
      %  "bsr_cluster_coords"
      %
      neighbor_coord = intersect(neighbor_coord, bsr_cluster_coords);

   else

      neighbor_coord = cur_coord;

   end;		% if neighbor_size > 0

   %  find out neighborhood indices in st_datamat
   %
   [ncoord ncoord_idx] = intersect(st_coords, neighbor_coord);

   neighbor_numbers = length(ncoord_idx);
   h = findobj(gcf,'Tag','neighbornumberEdit');
   set(h,'String',num2str(neighbor_numbers));

   st_datamat = getappdata(gcf,'STDatamat');	% indexed by common_coord

   [nr,nc] = size(st_datamat);
   cols = nc / st_win_size;

   if emp_st_data
      st_data = [];
   else
      st_data = reshape(st_datamat,[nr,st_win_size,cols]);
      st_data = mean(st_data(:,[1:num_pts_plotted],ncoord_idx),3);
   
      h = findobj(gcf,'Tag','NormalizePlots');
      normalize_flg = get(h,'UserData');
      if (normalize_flg == 1),
         max_st_data = max(st_data,[],2);
         min_st_data = min(st_data,[],2);

         max_mtx = max_st_data(:,ones(1,num_pts_plotted));
         min_mtx = min_st_data(:,ones(1,num_pts_plotted));
         scale_factor = max_mtx - min_mtx;
         st_data = (st_data - min_mtx) ./ scale_factor;
      end;
   end

   bsr = getappdata(gcf,'org_bsr');
   bs_thresh = getappdata(gcf,'bs_thresh');

   if ~isempty(bsr)
      [r c] = size(bsr);
      bsr = reshape(bsr, [st_win_size, r/st_win_size, c]);
      bsr = squeeze(bsr([1:num_pts_plotted],coord_idx,:));
      set(findobj(gcf,'tag','Stable'), 'enable', 'on');
   end 


   selected_condition = getappdata(gcf,'SelectedCondition');
   num_conditions = length(selected_condition);
   condition = cell(1,num_conditions);  

   max_num_stim = 0;
   for i=1:num_conditions,
      condition{i}.st_row_idx  = find(st_evt_list == i);
      condition{i}.num_stim = length(condition{i}.st_row_idx);
      condition{i}.name = conditions{i};
      if (max_num_stim < condition{i}.num_stim)
          max_num_stim = condition{i}.num_stim;
      end;
   end;


   % generate the plots
   %
   plot_cond_idx = find(selected_condition == 1);

   plot_dims = getappdata(gcf,'PlotDims');
   if isempty(plot_dims) 
       if (num_conditions < 5)
         num_rows = num_conditions;
       else
         num_rows = 5;
       end;
       if (max_num_stim < 4),
         num_cols = max_num_stim;
       else
         num_cols = 4;
       end;