fmri_cluster_report.m

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

function [cluster_info, cluster_fig] = fmri_cluster_report(varargin)
%
% USAGE:  cluster_info =  fmri_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) | all(st_origin == 0)

      if all(st_dims == [40 48 1 34]) & all(st_voxel_size == [4 4 4])
         st_origin = [20 29 12];
      elseif all(st_dims == [91 109 1 91]) & all(st_voxel_size == [2 2 2])
         st_origin = [46 64 37];
      else
         % according to SPM: if the origin field contains 0, then the origin is
         % assumed to be at the center of the volume.
         %
         st_origin = floor((dims([1 2 4])+1)/2);
         % st_origin = round(st_dims([1 2 4])/2);
      end;
   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);

   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;

   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;

   if (recursion_cnt >= recursion_limit)
      return;					% too many recursion
   end;

   neighbor_voxels = [  x-1 y   z;   x+1 y   z;
			x   y-1 z;   x   y+1 z;
      			x   y   z-1; x   y   z+1 ];
   for i=1:size(neighbor_voxels,1), 
      % exclude points that are out of boundary
      %
      curr_xyz = neighbor_voxels(i,:);
      within_boundary = and(curr_xyz>=[1 1 1],curr_xyz<=mask_dims);
      if isempty(find(~within_boundary)),
	  fillmask(curr_xyz,cnt);
      end;
   end;

   if (first_iteration)			% done all recursive calls, setup
      if ~isempty(combine_clusters) 
          same_clusters = unique(combine_clusters);
      end;
      updated_mask = mask;		%   return values
      c_size = cluster_size;
   end;

   return;					% fillmask


%-------------------------------------------------------------------------
%
function [xyz, xyz_mm] = coords2xyz(fig, coords, origin, voxel_size, dim);

   sa = getappdata(fig, 'sa');

   if isempty(sa)
      [xyz, xyz_mm] = rri_coord2xyz(coords, dim, origin, voxel_size);
   elseif sa
      origin = origin([3 1 2]);
      voxel_size = voxel_size([3 1 2]);
      [xyz, xyz_mm] = rri_coord2xyz(coords, dim, origin, voxel_size);
   else
      origin = origin([1 3 2]);
      voxel_size = voxel_size([1 3 2]);
      [xyz, xyz_mm] = rri_coord2xyz(coords, dim, origin, voxel_size);
   end

   return;

%
   if isempty(coords),
      xyz = [];
      xyz_mm = [];
      return;
   end;

   y_idx = repmat([1:st_dims(2)],1,st_dims(4)); 
   z_idx = repmat([1:st_dims(4)],st_dims(2),1); 
   
   x = mod(coords-1,st_dims(1)) + 1;
   col_idx = floor((coords-x)/st_dims(1));
%   y = st_dims(2) - y_idx(col_idx+1) + 1;
   y = y_idx(col_idx+1);
   z = z_idx(col_idx+1);

   xyz = [x(:) y(:) z(:)];

   xyz_offset = xyz - ones(length(x),1) * st_origin;
%   xyz_offset = [ x(:)-st_origin(1) ...
%		  (st_dims(2)-st_origin(2)+1)-y(:) ...
%	          z(:)-st_origin(3)];

   xyz_mm = xyz_offset * diag(voxel_size);

   return; 					% coords2xyz


%-------------------------------------------------------------------------
function cluster_fig = show_cluster_report(cluster_info,report_type,mainfig),
%
   cluster_fig = [];

   isbsr = getappdata(mainfig,'ViewBootstrapRatio');
   curr_lv_idx = getappdata(mainfig,'CurrLVIdx');
   cluster_blv = getappdata(mainfig,'cluster_blv');
   cluster_bsr = getappdata(mainfig,'cluster_bsr');

   if isbsr
      cluster_bsr{curr_lv_idx} = cluster_info;