pet_cluster_report.m

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

   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;
      setappdata(mainfig,'cluster_bsr',cluster_bsr);
   else
      cluster_blv{curr_lv_idx} = cluster_info;
      setappdata(mainfig,'cluster_blv',cluster_blv);
   end

   cluster_fig = create_report_figure(report_type,mainfig);
   setup_cluster_data_idx(cluster_info);

   setup_header;
   setup_footer;

   show_report_page(1);

   return; 					% show_cluster_report


%-------------------------------------------------------------------------
function show_report_page(page_num),
%
   report_type = getappdata(gcf,'ReportType');
   num_rows = getappdata(gcf,'NumRows');
   top_margin = getappdata(gcf,'TopMargin');
   line_height = getappdata(gcf,'LineHeight');
   cluster_info = getappdata(gcf,'ClusterInfo');
   cluster_data_idx = getappdata(gcf,'ClusterDataIdx');
   standard_text_obj = getappdata(gcf,'StandardTextObj');
   h_list = getappdata(gcf,'HdlList');

   axes_h = findobj(gcf,'Tag','ClusterReportAxes');
   axes_pos = get(axes_h,'Position');
   first_row_pos = 1 - top_margin;
   last_row_pos = first_row_pos - (num_rows-1)*line_height;

   if isempty(cluster_data_idx)
      return;
   end;

   show_page_keys(page_num);

   %  display the cluster information for the current page
   %
   start_row = num_rows * (page_num - 1) + 1;
   if (cluster_data_idx(1,start_row) == 0),	% the first row is empty
	start_row = start_row - 1;		
	num_rows = num_rows - 1;
   end;

   v_pos = [first_row_pos:-line_height:last_row_pos];

   if (report_type == 0)		% BLV value report
      h_pos = [0.06 0.18 0.28 0.44 0.71 0.88];
      h_pvalue_pos = 0;
   else					% Bootstrap ratio report
      h_pos = [0.06 0.18 0.28 0.44 0.67 0.88];
      h_pvalue_pos = 0.75;
   end;

   delete(h_list);
   h_list = [];
   last_row = min(size(cluster_data_idx,2),start_row+num_rows-1);
   for idx=start_row:last_row,

      v_idx = idx-start_row+1;
      cluster_id  = cluster_data_idx(1,idx);
      if (cluster_id ~= 0),
         cluster_lag = cluster_data_idx(2,idx);
         cluster_idx = cluster_data_idx(3,idx);
         c_data = cluster_info.data{cluster_lag};

         h = copyobj(standard_text_obj,axes_h);			% cluster # 
         cluster_num = sprintf('%3d',cluster_id);
         set(h, 'String',cluster_num, ...
		'Position',[h_pos(1) v_pos(v_idx) 0], ...
		'Visible','on'); 
	 h_list = [h_list h];

         h = copyobj(standard_text_obj,axes_h);			% voxel mean
         voxel_mean_str = sprintf('%6.2f',cluster_info.voxel_means_avg(cluster_idx));
         set(h, 'String',voxel_mean_str, ...
		'Position',[h_pos(2) v_pos(v_idx) 0], ...
		'Visible','on'); 
	 h_list = [h_list h];

	 h = copyobj(standard_text_obj,axes_h);			% peak xyz
         peak_xyz_str = sprintf('[%3d %3d %3d]',c_data.peak_xyz(cluster_idx,:));
	 set(h, 'String',peak_xyz_str, ...
		'Position',[h_pos(3) v_pos(v_idx) 0], ...
		'Visible','on'); 
	 h_list = [h_list h];

	 h = copyobj(standard_text_obj,axes_h);			% peak loc
         peak_loc_str = sprintf('[%6.1f %6.1f %6.1f]', ...
					c_data.peak_loc(cluster_idx,:));
	 set(h, 'String',peak_loc_str, ...
		'Position',[h_pos(4) v_pos(v_idx) 0], ...
		'Visible','on'); 
	 h_list = [h_list h];

	 h = copyobj(standard_text_obj,axes_h);			% peak value
         peak_value_str = sprintf('%8.4f',c_data.peak_values(cluster_idx));
	 set(h, 'String',peak_value_str, ...
		'Position',[h_pos(5) v_pos(v_idx) 0], ...
		'Visible','on'); 
	 h_list = [h_list h];

	 if (h_pvalue_pos ~= 0)
	    h = copyobj(standard_text_obj,axes_h);		% P value
	    p_value = ratio2p(abs(c_data.peak_values(cluster_idx)), ...
						0,1 );
            p_value_str = sprintf('(%6.4f)',p_value);
	    set(h, 'String',p_value_str, ...
		'Position',[h_pvalue_pos v_pos(v_idx) 0], ...
		'Visible','on'); 
	    h_list = [h_list h];
         end;

	 h = copyobj(standard_text_obj,axes_h);			% cluster size
         size_str = sprintf('%4d',c_data.size(cluster_idx));
	 set(h, 'String',size_str, ...
		'Position',[h_pos(6) v_pos(v_idx) 0], ...
		'Visible','on'); 
	 h_list = [h_list h];
      end;
	
   end;

   setappdata(gcf,'CurrPage',page_num);
   setappdata(gcf,'HdlList',h_list);

   return; 					% show_report_page


%--------------------------------------------------------------------------
function  save_cluster_report()

   report_type = getappdata(gcf,'ReportType');
   cluster_info = getappdata(gcf,'ClusterInfo');

   if (report_type == 0)		% BLV value report
      cluster_info.type = 'BLV Report';
   else
      cluster_info.type = 'Bootstrap Ratio Report';
   end;

   result_name = getappdata(gcbf, 'ResultName');

   [filename, pathname] = ...
           uiputfile('*_PETcluster.mat','Save to .mat file');

   if ischar(filename) & (length(filename)<9 | isempty(findstr(lower(filename),'_petcluster')))
      [tmp filename] = fileparts(filename);
      filename = [filename, '_PETcluster.mat'];
   end
   
   if isequal(filename,0)
         status = 0;
         return;
   end;
   report_file = fullfile(pathname,filename);

   try
      save (report_file, 'cluster_info', 'report_type' );
      msg = ['File ', filename, ' has been saved'];
      msgbox(msg, 'Info');
   catch
      msg = ['File ', filename, ' can not be open to write'];
      msgbox(msg, 'Error');
%      msg = sprintf('Cannot save cluster report to %s',report_file),
 %     set(findobj(gcf,'Tag','MessageLine'),'String',['ERROR: ', msg]);
      status = 0;
      return;
   end;

   return; 					% save_cluster_report


%--------------------------------------------------------------------------
function  save_cluster_txt()

   [filename, pathname] = ...
           uiputfile('*_PETcluster.txt','Save to .txt file');

   if ischar(filename) & (length(filename)<9 | isempty(findstr(lower(filename),'_petcluster.txt')))
      [tmp filename] = fileparts(filename);
      filename = [filename, '_PETcluster.txt'];
   end
   
   if isequal(filename,0)
         status = 0;
         return;
   end;
   report_file = fullfile(pathname,filename);

   fid = fopen(report_file, 'wt');

   if fid == -1
      msg = ['File ', filename, ' can not be open to write'];
      msgbox(msg, 'Error');
      return;
   end

   cluster_info = getappdata(gcf,'ClusterInfo');
   cluster_data_idx = getappdata(gcf, 'ClusterDataIdx');
   report_type = getappdata(gcf,'ReportType');

   str = sprintf('Clu#\tV_Mean\tX\tY\tZ\tX(mm)\tY(mm)\tZ(mm)\t');

   if report_type
      str = [str sprintf('BSR\tAppro.P\tClu_Size(voxels)\n')];
      str = [str '--------------------------------------------------------'];
      str = [str '----------------------------------------'];
   else
      str = [str sprintf('BLV\tClu_Size(voxels)\n')];
      str = [str '--------------------------------------------------------'];
      str = [str '--------------------------------'];
   end

   fprintf(fid, '\n%s\n\n', str);

   for idx = 1 : size(cluster_data_idx, 2)
      cluster_id = cluster_data_idx(1, idx);

      if cluster_id ~= 0
         cluster_lag = cluster_data_idx(2, idx);
         cluster_idx = cluster_data_idx(3, idx);
         c_data = cluster_info.data{cluster_lag};

         cluster_num = sprintf('%d\t', cluster_id);
         voxel_mean_str = ...
		sprintf('%.2f\t',cluster_info.voxel_means_avg(cluster_idx));
         peak_xyz_str = ...
		sprintf('%d\t%d\t%d\t',c_data.peak_xyz(cluster_idx,:));
         peak_loc_str = ...
		sprintf('%.1f\t%.1f\t%.1f\t',c_data.peak_loc(cluster_idx,:));
         peak_value_str = sprintf('%.4f\t',c_data.peak_values(cluster_idx));

         if report_type
            p_value = ...
		ratio2p(abs(c_data.peak_values(cluster_idx)), 0, ...
			1);
            p_value_str = sprintf('%.4f\t',p_value);
         else
            p_value_str = '';
         end

         size_str = sprintf('%d\n',c_data.size(cluster_idx));

         str = [cluster_num,voxel_mean_str,peak_xyz_str,peak_loc_str];
         str = [str peak_value_str,p_value_str,size_str];
         fprintf(fid, '%s', str);
      end
   end

   %  display the footer
   %
   source_file = cluster_info.source;
   source_line = sprintf('Source: %s',source_file);

   lv_str = sprintf('LV Idx: %d',cluster_info.lv_idx);
   thresh_str = sprintf('Threshold: %2.4f',cluster_info.threshold);
   min_size_str = sprintf('Min Size: %d(voxels)',cluster_info.min_size);
   min_dist_str = sprintf('Min Distance: %3.1fmm ',cluster_info.min_dist);
   parameters_line = sprintf('%s, %s, %s, %s', ...
				lv_str,thresh_str,min_size_str,min_dist_str);

   if report_type
      str = ['--------------------------------------------------------'];
      str = [str '----------------------------------------'];
      str = [str sprintf('\n%s\n%s',source_line,parameters_line)];
   else
      str = ['--------------------------------------------------------'];
      str = [str '--------------------------------'];
      str = [str sprintf('\n%s\n%s',source_line,parameters_line)];
   end

   fprintf(fid, '\n%s\n\n', str);

   fclose(fid);

   msg = ['File ', filename, ' has been saved'];
   msgbox(msg, 'Info');

   return; 					% save_cluster_txt


%--------------------------------------------------------------------------
function  save_peak_location()

   result_file = get(findobj(getappdata(gcf,'mainfig'),'tag','ResultFile'),'UserData');
   cluster_info = getappdata(gcf,'ClusterInfo');

   xyz = cluster_info.data{1}.peak_xyz;
   xyz = double(xyz2xyzmm(xyz, result_file));

   [filename, pathname] = ...
           uiputfile('*_PETcluster_peak.txt','Export Peak Location');

   if ischar(filename) & (length(filename)<9 | isempty(findstr(lower(filename),'_petcluster_peak.txt')))
      [tmp filename] = fileparts(filename);
      filename = [filename, '_PETcluster_peak.txt'];
   end
   
   if isequal(filename,0)
         status = 0;
         return;
   end;
   report_file = fullfile(pathname,filename);

   try
      save (report_file, '-ascii', 'xyz' );
      msg = ['File ', filename, ' has been saved'];
      msgbox(msg, 'Info');
   catch
      msg = ['File ', filename, ' can not be open to write'];
      msgbox(msg, 'Error');
%      msg = sprintf('Cannot export peak location to %s',report_file),
 %     set(findobj(gcf,'Tag','MessageLine'),'String',['ERROR: ', msg]);
      status = 0;
      return;
   end;

   return; 					% save_peak_location


%--------------------------------------------------------------------------
function  save_all_location()

   result_file = get(findobj(getappdata(gcf,'mainfig'),'tag','ResultFile'),'UserData');
   dims = getappdata(getappdata(gcf,'mainfig'),'Dims');
   cluster_info = getappdata(gcf,'ClusterInfo');

   xyz = cluster_info.data{1}.idx;
   xyz = rri_coord2xyz(xyz, dims);

   cnum = cluster_info.data{1}.mask;
   for i=1:length(cluster_info.data{1}.id)
      cnum(find(cluster_info.data{1}.mask==cluster_info.data{1}.id(i)))=i; 
   end;

   coord = getappdata(getappdata(gcf,'mainfig'),'BSRatioCoords');

   if isempty(coord)
      coord = getappdata(getappdata(gcf,'mainfig'),'BLVCoords');
   end

   if getappdata(gcf,'ReportType')
      data = getappdata(getappdata(gcf,'mainfig'), 'BSRatio');
   else
      data = getappdata(getappdata(gcf,'mainfig'), 'BLVData');
   end

   data = data(:, cluster_info.lv_idx);

   coord_idx = ones(size(cluster_info.data{1}.idx));

   for i=1:length(coord_idx)
      coord_idx(i) = find(coord==cluster_info.data{1}.idx(i));