pet_cluster_report.m

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

function [cluster_info, cluster_fig] = pet_cluster_report(varargin)
%
% USAGE:  cluster_info =  pet_cluster_report(min_size, min_dist)
%
%   Input:  
%	min_size - minimum number of voxels within a cluster (default = 5)
%	min_dist - minimum distant (in mm) between the peaks of any two
%                  clusters (default = 4mm)
%

   cluster_info = [];
   cluster_fig = [];

   if (nargin == 0) | ~ischar(varargin{1})
       if (nargin >= 1)
          min_size = varargin{1};
       else
          min_size = 5;
       end;

       if (nargin >= 2)
          min_dist = varargin{2};
       else
          min_dist = 10;
       end;

       mainfig = gcbf;

       cluster_info = get_cluster_info(min_size,min_dist);

       if ~isempty(cluster_info)
          report_type = getappdata(gcbf,'ViewBootstrapRatio');
          cluster_fig = show_cluster_report(cluster_info,report_type,mainfig);
       end;
       return;
   end;

   action = varargin{1};
   switch (action),
     case {'initialize'}
	 cluster_info = varargin{2};
	 header_type = varargin{3};
         show_cluster_report(cluster_info,header_type,[]);
     case {'PrevPage'}
         curr_page = getappdata(gcbf,'CurrPage');
         show_report_page(curr_page-1);
     case {'NextPage'}
         curr_page = getappdata(gcbf,'CurrPage');
         show_report_page(curr_page+1);
      case {'LoadClusterReport'}
	 cluster_fig = load_cluster_report(varargin{2});
      case {'SaveClusterReport'}
	 save_cluster_report;
      case {'SaveClusterTXT'}
	 save_cluster_txt;
      case {'SaveAllLocation'}
         save_all_location;
      case {'SavePeakLocation'}
         save_peak_location;
      case {'delete_fig'}
         delete_fig;
   end;  

   return;					% find_active_cluster
   

%-------------------------------------------------------------------------
function [cluster_info] = get_cluster_info(min_size,min_dist),
%
%  This function uses data from the "pls_result_ui" figure

   %  make sure the threshold value are used for the current plotting
   %
   view_state = getappdata(gcbf,'ViewBootstrapRatio');
   if (view_state == 0),
      thresh_field = str2num(get(findobj(gcbf,'Tag','Threshold'),'String'));
      plotted_thresh = getappdata(gcbf,'BLVThreshold');
   else
      thresh_field = str2num(get(findobj(gcbf,'Tag','BSThreshold'),'String'));
      plotted_thresh = getappdata(gcbf,'BSThreshold');
   end;

   thresh_diff = abs(plotted_thresh - thresh_field);

   if (thresh_diff > 0.0001)
      msg = 'Threshold field have been changed.  Press PLOT button first';
      set(findobj(gcbf,'Tag','MessageLine'),'String',['ERROR: ' msg]);
      return;
   end;

   
   %  extract the brain LV data or Bootstrap Ratio
   %

   lv_idx = getappdata(gcbf,'CurrLVIdx');
   st_dims = getappdata(gcbf,'STDims');
   win_size = getappdata(gcbf,'WinSize');
   voxel_size = getappdata(gcbf,'STVoxelSize');
   st_origin = getappdata(gcbf,'STOrigin');
   PLSResultFile = get(findobj(gcbf,'Tag','ResultFile'),'UserData');

   if isempty(st_origin)
      st_origin = round(st_dims([1 2 4])/2);
   end;

   if (view_state == 0),
      data = getappdata(gcbf,'BLVData');
      std_ratio = [];
      st_coords = getappdata(gcbf,'BLVCoords');
   else
      data = getappdata(gcbf,'BSRatio');
      curr_bs_ratio = data(:,lv_idx);
      curr_bs_ratio(isnan(curr_bs_ratio)) = 0;

      % avoid the outliers
      idx = find(abs(curr_bs_ratio) < std(curr_bs_ratio) * 5); 
      std_ratio = std(curr_bs_ratio(idx));
      st_coords = getappdata(gcbf,'BSRatioCoords');
   end;

   %  extract voxels that have values above the threshold
   %
   pos_active_list = cell(1,win_size);
   neg_active_list = cell(1,win_size);
   for i=1:win_size,
      voxels = data(i:win_size:end,lv_idx);

      pos_active_voxel_idx = find(voxels >= plotted_thresh);
      pos_active_voxels = voxels(pos_active_voxel_idx);

      neg_active_voxel_idx = find(voxels <= -plotted_thresh);
      neg_active_voxels = voxels(neg_active_voxel_idx);

      [dummy_voxels, sorted_order] = sort(pos_active_voxels);
      sorted_order = sorted_order(end:-1:1);
      pos_active_list{i}.voxels = pos_active_voxels(sorted_order)';
      pos_active_list{i}.coords = st_coords(pos_active_voxel_idx(sorted_order));

      [dummy_voxels, sorted_order] = sort(neg_active_voxels);
      neg_active_list{i}.voxels = neg_active_voxels(sorted_order)';
      neg_active_list{i}.coords = st_coords(neg_active_voxel_idx(sorted_order));
   end;


   %  now compute the cluster information of the active voxels
   %
   cluster_info.source = PLSResultFile;
   cluster_info.min_size = min_size;
   cluster_info.min_dist = min_dist;
   cluster_info.threshold = plotted_thresh;
   cluster_info.std_ratio = std_ratio;
   cluster_info.lv_idx = lv_idx;
   cluster_info.data = cell(1,win_size);


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


   for i=1:win_size
      pos_cluster = compute_cluster_stat(gcbf, pos_active_list{i}, st_origin, ...
					voxel_size,st_dims,min_size,min_dist);
      neg_cluster = compute_cluster_stat(gcbf, neg_active_list{i}, st_origin, ...
					voxel_size,st_dims,min_size,min_dist);

      % put the positive and negative cluster info together
      %
      curr_cluster.id = [ pos_cluster.id -neg_cluster.id ];
      curr_cluster.size = [ pos_cluster.size neg_cluster.size ];
      curr_cluster.peak_values = ...
		   	[ pos_cluster.peak_values neg_cluster.peak_values ];
      curr_cluster.peak_coords = ...
			[ pos_cluster.peak_coords neg_cluster.peak_coords ];
      curr_cluster.peak_xyz = [ pos_cluster.peak_xyz; neg_cluster.peak_xyz ];
      curr_cluster.peak_loc = [ pos_cluster.peak_loc; neg_cluster.peak_loc ];
      curr_cluster.idx = [ pos_cluster.idx neg_cluster.idx ];
      curr_cluster.mask = [ pos_cluster.mask -neg_cluster.mask ];

      sort_by_cluster=[];
      for j=1:length(curr_cluster.id)
         sort_by_cluster = [sort_by_cluster find(curr_cluster.mask==curr_cluster.id(j))];
      end
      curr_cluster.idx = curr_cluster.idx(sort_by_cluster);
      curr_cluster.mask = curr_cluster.mask(sort_by_cluster);

      cluster_info.data{i} = curr_cluster;
   end;

   % calculate voxel mean for the subjects in a task/scan
   %
   datamat_file_lst = getappdata(gcbf,'DatamatFileList');
   datamat_file = load(datamat_file_lst{1},'coords','datamat','session_info');



    origin_pattern = getappdata(gcbf,'origin_pattern');

    if ~isempty(origin_pattern)
       dims = getappdata(gcbf,'STDims');
       new_coord = zeros(dims);
       new_coord(datamat_file.coords) = 1;
       new_coord = find(new_coord(origin_pattern));

       datamat_file.datamat = ...
		rri_xy_orient_data(datamat_file.datamat', ...
		datamat_file.coords, new_coord, dims, origin_pattern);
       datamat_file.datamat = datamat_file.datamat';

       datamat_file.coords = new_coord;
    end



   cond_selection = [];

   warning off;
   try
      load(PLSResultFile,'cond_selection');
   catch
   end
   warning on;

   num_subj = datamat_file.session_info.num_subjects;

   if isempty(cond_selection)
      num_cond = datamat_file.session_info.num_conditions;
      cond_selection = ones(1,num_cond);
   else
      num_cond = sum(cond_selection);
   end

   selected_subjects = ones(num_subj,1);
   bmask = selected_subjects * cond_selection;
   bmask = find(bmask(:));
   datamat_file.datamat = datamat_file.datamat(bmask,:);

   % coord is the coord of all active voxels
   % peak_coord is the coord of peak active voxels
   %
   coord = curr_cluster.idx;
   peak_coord = curr_cluster.peak_coords;

   % calculate index of coord & peak_coord
   %
   for i=1:length(coord)
      coord_idx(i) = find(datamat_file.coords == coord(i));
   end

   % for i=1:length(peak_coord)
   %    peak_coord_idx(i) = find(datamat_file.coords == peak_coord(i));
   % end

   % find location of peak_coord_idx inside coord_idx
   %
   % [tmp,peak_idx] = intersect(coord,peak_coord);

   for i=1:length(peak_coord)
      peak_idx(i) = find(peak_coord(i) == coord);
   end


   % calc intensity
   %
   for v = 1:length(coord_idx)			% traverse active voxels
      for k = 1:num_cond			% cond
         for n = 1:num_subj			% subj
            j = n+(k-1)*num_subj;		% row number in datamat
            intensity(v,k,n) = ...
		datamat_file.datamat(j, coord_idx(v));
         end
      end
   end

%   set(gcf,'Pointer',old_pointer);

   % voxel_means are the mean for subjects in a task/scan
   % voxel_means_avg are mean across all scan, for display
   %
   voxel_means = mean(intensity,3);
   voxel_means = voxel_means([peak_idx],:);

   voxel_means_avg = mean(voxel_means,2);
%   voxel_means_avg = voxel_means_avg';
%   voxel_means_avg = voxel_means_avg([peak_idx])';

   cluster_info.voxel_means = voxel_means;
   cluster_info.voxel_means_avg = voxel_means_avg;

   return;					% get_cluster_info


%-------------------------------------------------------------------------
function [cluster_info] = compute_cluster_stat(fig, active_list,st_origin, ...
					voxel_size,st_dims,min_size,min_dist),
%   
   num_active_voxels = length(active_list.voxels);

   mask = zeros(st_dims([1 2 4]));
   mask(active_list.coords) = -1;

   y_idx = repmat([1:st_dims(2)],1,st_dims(4)); 
   z_idx = repmat([1:st_dims(4)],st_dims(2),1); 

   coords = active_list.coords;
   cluster_id = zeros(1,num_active_voxels);
   cluster_size = zeros(1,num_active_voxels);
   for i=1:num_active_voxels,
       x = mod(coords(i)-1,st_dims(1)) + 1;
       col_idx = floor((coords(i)-x)/st_dims(1));
       y = y_idx(col_idx+1);
       z = z_idx(col_idx+1);

       [updated_mask,csize,same_clusters] = fillmask([x y z],i,mask);

       if ~isempty(updated_mask)
          cluster_id(i) = i;
          cluster_size(i) =  csize;
	  mask = updated_mask;

          if ~isempty(same_clusters)
	     same_clusters = [same_clusters i];
             min_cnt = min(cluster_id(same_clusters));
	     cluster_id(same_clusters) = min_cnt;
          end;
       end;
   end;

   %  merge clusters when needed
   %
   active_voxels = find(mask > 0);
   active_clusters = [];
   for i=1:num_active_voxels,
      if (cluster_id(i) ~= 0),

         curr_cluster_id = i;
	 cluster_set = find(cluster_id == curr_cluster_id);
         total_cluster_size = sum(cluster_size(cluster_set));

	 cluster_id(cluster_set) = 0;
	 cluster_size(cluster_set) = 0;
         if (total_cluster_size > min_size)
            active_clusters = [active_clusters curr_cluster_id];
	    cluster_id(curr_cluster_id) = curr_cluster_id;
   	    cluster_size(curr_cluster_id) = total_cluster_size;
         else
            curr_cluster_id = 0;
         end;

	 for j=1:length(cluster_set),
	    idx = find(mask(active_voxels) == cluster_set(j));
            mask(active_voxels(idx)) = curr_cluster_id;
         end;
      end;
   end;


   cluster_id = cluster_id(active_clusters);
   cluster_size = cluster_size(active_clusters);

   peak_values = active_list.voxels(cluster_id); peak_values(:)';
   peak_coords = active_list.coords(cluster_id); peak_coords(:)';
   [peak_xyz, peak_loc] = coords2xyz(fig, peak_coords,st_origin,voxel_size,st_dims);
   mask_idx = find(mask > 0);
   mask_value = mask(mask_idx);

   %  combine clusters that are close together (i.e. their peaks are less
   %  than the distance specified by min_dist);
   %
   min_dist_sq = min_dist^2;
   num_clusters = length(cluster_id);
   is_not_merged = ones(1,num_clusters);
   for i=1:num_clusters-1,
     if (is_not_merged(i)),
        for j=i+1:num_clusters,
         if (is_not_merged(j)),
          dist = (peak_xyz(i,:) - peak_xyz(j,:)) .* voxel_size;
          if ((dist * dist') <= min_dist_sq),	   % cluster needed to be merged
	     cluster_size(i) = cluster_size(i)+cluster_size(j);
	     mask_value(find(mask_value==cluster_id(j))) = cluster_id(i);
	     is_not_merged(j) = 0;
          end;
         end
        end;
     end;
   end;


   %  setup the cluster record 
   %
   final_cluster_idx = find(is_not_merged == 1);
   cluster_info.id = cluster_id(final_cluster_idx);
   cluster_info.size = cluster_size(final_cluster_idx);
   cluster_info.peak_values = peak_values(final_cluster_idx);
   cluster_info.peak_coords = peak_coords(final_cluster_idx);
   cluster_info.peak_xyz = peak_xyz(final_cluster_idx,:);
   cluster_info.peak_loc = peak_loc(final_cluster_idx,:);
   cluster_info.idx = mask_idx(:)';
   cluster_info.mask = mask_value(:)';

   return;					% compute_cluster_stat


%-------------------------------------------------------------------------
%
function [updated_mask,c_size,same_clusters] = fillmask(xyz,cnt,initial_mask),
%   
   persistent mask
   persistent mask_dims
   persistent cluster_size
   persistent start_xzy
   persistent recursion_cnt
   persistent recursion_limit
   persistent combine_clusters

   updated_mask = []; 
   same_clusters = [];
   c_size = 0;

   if exist('initial_mask','var')  
      first_iteration = 1;
      mask = initial_mask;
      mask_dims = size(mask);
      cluster_size = 0;
      start_xyz = xyz;
      combine_clusters = [];
      recursion_cnt = 1;
      recursion_limit = get(0,'RecursionLimit');
   else
      recursion_cnt = recursion_cnt+1;
      first_iteration = 0;
   end;

   x=xyz(1); y=xyz(2); z=xyz(3);
   curr_voxel_value = mask(x,y,z);
   if (curr_voxel_value == 0 | curr_voxel_value == cnt)
      return;
   end;

   %  handle the missing voxels for the previous clusters due to 
   %  run out of stack space (i.e. recursive limitation)
   if (curr_voxel_value > 0)
      combine_clusters = [combine_clusters curr_voxel_value];
      return;
   end;						  
   
   mask(x,y,z) = cnt;
   cluster_size = cluster_size+1;