www.pudn.com > colorseg.zip > calculatePXTheta.m, change:2003-03-11,size:3398b

function M = calculatePXTheta(im,prevClusters)
% CALCULATEPXTHETA - Calculate the probability of each pixel being
% its color conditioned on all of the clusters that were found at
% the previous (coarser) iteration.
%
% Input:
%      IM - MxNx3 image, where each pixel is a point in a
%      perceptually uniform color space.
%
%      PREVCLUSTERS - MxNx3 matrix that contains binary masks
%      representing a pixels membership in a certain cluster as of
%      the previous scale.
%
% Output:
%      M - MxNxC matrix, where the (i,j,k) entry is the probability
%      of pixel (i,j) being its color given that it is a member of
%      cluster k.
%
% Notes:
%
%      The internal variable BINS is rather important because it
%      controls the courseness of the pdf estimation based on 3D
%      histograms in color space. The colorspace cube is divided
%      into BINS^3 sub-cubes for the estimation.  If BINS is large,
%      then the probability of a given pixel given the clusters
%      tends to be 1 for one cluster and 0 for all the rest (thus
%      not really adding any information).  
  
  BINS = 50;  % the number of bins to use for the histrograms along
  imSize = size(im);
  imSize = imSize(1:2);
  numPixels = prod(imSize);
  
  im = reshape(im,prod(imSize),3);
  
  CSPACE_MIN = -0.1; % minimum value allowed in the color space
  CSPACE_MAX = 1.0; % maximum value allowed in the color space 
  C = length(prevClusters(1,1,:)); % the number clusters
  
  % Create 3D color histograms for each of the clusters
  edges = repmat([-inf linspace(CSPACE_MIN,CSPACE_MAX,BINS-1)]',1,3);
  
  %prepare a matrix to hold all of the histogram data
  hists = zeros(BINS,BINS,BINS,C);
  % reshape the image so that we have three 1D vectors
  for i=1:C
    cluster_linear = reshape(prevClusters(:,:,i), ...
			     prod(size(prevClusters(:,:,i))),1);
    ind = find(cluster_linear);
    % im_cluster will be a Px3 matrix, where each row is the color
    % point (in R3) of a point within the image that is in the i^th
    % cluster, and P is the total number of pixels in that cluster
    im_cluster = im(ind,:);
    cluster_size = length(im_cluster(:,1));
    
    [n,x,nbins] = histmulti5(im_cluster,edges);
    % the n that was returned is the BINSxBINSxBINS vector that
    % contains counts for each bin
    n = n / cluster_size;
    hists(:,:,:,i) = n;
    % normalize the histograms to turn them into probabilities
  end
   
   
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   % now, for each pixel in the image, figure out which
   % histogram bin it falls into and use this information to assign
   % a probability to this pixel coming from each of the clusters
   
   M = zeros(numPixels,C);
   edges = edges(:,1);

   % there's probably a better way to do this than the dreaded
   % for-loop, but i can't figure it out
   for i=1:numPixels
     for j=1:C
       M(i,j) = hists(max(find(im(i,1)>edges)),(max(find(im(i,2)>edges))),(max(find(im(i,3)>edges))),j);
     end
   end
   
   M = reshape(M,imSize(1),imSize(2),C);
   
   % there is a normalization problem, so fix it here
   sums = sum(M,3);
   
   %%% hack alert!
   [ii,jj] = find((sums  eps));
   
   if (length(ii) > 0)
     sums(sub2ind(size(sums),ii,jj)) = 1;
     disp(sprintf('calculatePXTheta: found %d zero probabilities', ...
		  length(ii)));
   end
     
   for i=1:C
     M(:,:,i) = M(:,:,i) ./ sums;
   end