pet_plot_rf_task.m

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

       setappdata(gcf,'PlottedDatamat',selected_data);
       setappdata(gcf,'STDatamat',st_datamat);
       setappdata(gcf,'org_coords',org_coords);
       setappdata(gcf,'num_subj',num_subj);
       setappdata(gcf,'subj_name',subj_name);

       set(findobj(gcf,'Tag','MessageLine'),'String','');
   else
       org_coords = getappdata(gcf,'org_coords');
   end;
   
   %  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);
       set(gca,'visible','off');
       return;
   end;

   coord_idx = find(st_coords == cur_coord);
   if isempty(coord_idx)
       msg = 'ERROR: The selected point is outside the brain region.';
       set(findobj(gcf,'Tag','MessageLine'),'String',msg);
       set(gca,'visible','off');
       return;
   end;
   coord_idx = find(org_coords == cur_coord);

   %  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));
         all_voxel(st_coords) = (abs(BSRatio) >= BSThreshold);
         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(org_coords, neighbor_coord);

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

   st_datamat = getappdata(gcf,'STDatamat'); 

%   [nr,nc] = size(st_datamat);
%   cols = nc / st_win_size;
%   st_data = reshape(st_datamat,[nr,st_win_size,cols]);
%   st_data = squeeze(st_data(:,[1:num_pts_plotted],coord_idx));

   st_data = st_datamat;
   
%   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;

   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);
%      if (max_num_stim < condition{i}.num_stim)
%          max_num_stim = condition{i}.num_stim;
%      end;
%   end;

   % generate the plots
   %
   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;
%       plot_dims = [num_rows num_cols];
%   end;

%   pet_plot_cond_stim_ui('STARTUP', st_data, axes_margin, cond_idx);

   setappdata(gcf,'ST_data',st_data);
   setappdata(gcf,'PlotCondIdx',cond_idx);

   if isempty(getappdata(gcf,'CombinePlots'))
      setappdata(gcf,'CombinePlots',0);
   end;

   if isempty(getappdata(gcf,'ShowAverage'))
      setappdata(gcf,'ShowAverage',0);
   end;

   %  the following code is to get an intensity array for the voxel and plot it

   %  find out the coord_idx
   %
%   coord = getappdata(gcf,'Coord');
%   org_coords = getappdata(gcf,'org_coords');
%   coord_idx = find(org_coords == coord);

   num_subj = getappdata(gcf,'num_subj');
   subj_name = getappdata(gcf,'subj_name');
   subjects = [subj_name, {'mean'}];

   for k = cond_idx			% cond
      for n = 1:num_subj		% subj
         j = n+(k-1)*num_subj;		% row number in datamat
         intensity(k,n) = mean(st_data(j, ncoord_idx), 2);
         intensity_hdl(k,n) = 0;	% initialization
      end
      intensity_avg_hdl(k) = 0;
   end

   intensity_avg = mean(intensity,2);

   color_code =[ 'bo';'rx';'g+';'m*';'bs';'rd';'g^';'m<';'bp';'r>'; ...
                 'gh';'mv';'ro';'gx';'m+';'b*';'rs';'gd';'m^';'b<'];

   % need more color
   %
   if num_subj+1 > size(color_code,1)

      tmp = [];

      for i=1:ceil((num_subj+1)/size(color_code,1))
         tmp = [tmp; color_code];
      end

      color_code = tmp;

   end

   cla; grid off; hold on;
   for k = cond_idx

      for n = 1:num_subj
         intensity_hdl(k,n) = plot(k,intensity(k,n), ...
			color_code(n,:));
      end

      intensity_avg_hdl(k) = bar(k,intensity_avg(k));
      set(intensity_avg_hdl(k),'facecolor','none')
      % set(intensity_avg_hdl(k), 'linewidth', 2);

   end


   % normalize with intensity(:), not st_data(:)
   %
   min_y = min(intensity(:)); max_y = max(intensity(:));
   margin_y = abs((max_y - min_y) / 100);

   axis_scale = getappdata(gcf, 'axis_scale');

   if isempty(axis_scale)
      axis_scale = [0 length(cond_idx)+1 min_y-margin_y max_y+margin_y];
   else
      if (min_y-margin_y)<axis_scale(3)
         axis_scale(3) = min_y-margin_y;
      end;

      if (max_y+margin_y)>axis_scale(4)
         axis_scale(4) = max_y+margin_y;
      end;
   end

   setappdata(gcf, 'axis_scale', axis_scale);
   axis(axis_scale);
   set(gca, 'xtick', cond_idx);

   xlabel('Conditions');
   ylabel('Intensities');

   xyz = getappdata(gcf,'XYZ');
   title(['Intensity for subjects and mean at voxel:  [',num2str(xyz),']']);

   hold off;

   l_hdl = [];

   if ~isempty(subjects),

      intensity_legend = [intensity_hdl(1,:), intensity_avg_hdl(1)];

      % remove the old legend to avoid the bug in the MATLAB5
      old_legend = getappdata(gcf,'LegendHdl');
      if ~isempty(old_legend),
        try
          delete(old_legend{1});
        catch
        end;
      end;

      % create a new legend, and save the handles
      [l_hdl, o_hdl] = legend(intensity_legend, subjects, 0);
      legend_txt(o_hdl);
      set(l_hdl,'color',[0.9 1 0.9]);
      setappdata(gcf,'LegendHdl',[{l_hdl} {o_hdl}]);

   else

      setappdata(gcf,'LegendHdl',[]);

   end;

   setappdata(gcf, 'ncoord_idx', ncoord_idx);
   setappdata(gcf, 'ncoord', ncoord);
   set(findobj(gcf,'Tag','DataMenu'),'visible','on');

   return;					% plot_response_fn


%--------------------------------------------------------------------------
%
function  save_response_fn()

   neighbor_numbers = 1;

   %  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

   st_coords = getappdata(gcf,'STCoords');
   org_coords = getappdata(gcf,'org_coords');
   xyz = getappdata(gcf,'XYZ');

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

   coord_idx = find(st_coords == cur_coord);
   if isempty(coord_idx)
       msg = 'ERROR: The selected point is outside the brain region.';
       set(findobj(gcf,'Tag','MessageLine'),'String',msg);
       return;
   end;
%   coord_idx = find(org_coords == cur_coord);
   ncoord_idx = getappdata(gcf, 'ncoord_idx');
   neighbor_numbers = length(ncoord_idx);

   CallBackFig = getappdata(gcf,'CallBackFig');
   dims = getappdata(CallBackFig,'STDims');
   ncoord = getappdata(gcf, 'ncoord');
   xyz = rri_coord2xyz(ncoord,dims);

   %  get selected file names 
   %
   st_files = getappdata(gcf,'STFiles');
   h = findobj(gcf,'Tag','STDatamatPopup');
   selected_idx = get(h,'Value');
   selected_files = st_files{selected_idx}.fullname;

   %  extract the data of the selected voxel
   %
   st_data = getappdata(gcf,'STDatamat');
   data = mean(st_data(:, ncoord_idx), 2);
   standard_deviation = std(st_data(:, ncoord_idx), 0, 2);

   [filename, pathname] = ...
	uiputfile('*_PET_rf_plot.mat','Save the Response Functions');

   if ischar(filename) & isempty(findstr(lower(filename),'_pet_rf_plot'))
      [tmp filename] = fileparts(filename);
      filename = [filename, '_PET_rf_plot.mat'];
   end

   if isequal(filename,0)
	return;
   end;
   rf_plot_file = fullfile(pathname,filename);

   xyz_str = 'xyz';
 
   sa = getappdata(gcf, 'sa');
   
   if ~isempty(sa) & sa == 1
      yzx = xyz;
      xyz_str = 'yzx';
   elseif ~isempty(sa) & sa == 0
      xzy = xyz;
      xyz_str = 'xzy';
   end

   try
     save (rf_plot_file, 'selected_files', 'data', 'standard_deviation', 'xyz', 'neighbor_size', 'neighbor_numbers' );
   catch
     msg = sprintf('Cannot save the response function data to %s',rf_plot_file);
     set(findobj(gcf,'Tag','MessageLine'),'String',['ERROR: ', msg]);
     status = 0;
     return;
   end;

   return;                                      % save_response_fn


%--------------------------------------------------------------------------
%
function  save_response_fn2()

   st_coords = getappdata(gcf,'STCoords');
   org_coords = getappdata(gcf,'org_coords');
   xyz = getappdata(gcf,'XYZ');

   h = findobj(gcf,'Tag','STDatamatPopup');
   selected_idx = get(h,'Value');

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

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

   %  extract the data of the selected voxel
   %
   st_data = getappdata(gcf,'STDatamat');
   data = double(st_data(:,coord_idx));

   pattern = ['*_PET_grp', num2str(selected_idx), '_behavdata.txt'];
   pattern_suffix = ['_PET_grp', num2str(selected_idx), '_behavdata.txt'];

   [filename, pathname] = uiputfile(pattern,'Save Behav Data');

   if ischar(filename) & isempty(findstr(lower(filename), lower(pattern_suffix)))
      [tmp filename] = fileparts(filename);
      filename = [filename, pattern_suffix];
   end

   if isequal(filename,0)