function [ X, Y ]  = peak_find_cm(input, threshold)
% [ X, Y ] = peak_find_cm( input, threshold )
%
% Find the peaks in the two-dimensional input array.  This is done
% by the following scheme:
%
% A threshold (supplied) is applied to the input array to form a
% binary image.  Connected components are identified in this
% thresholded image, and from each connected component we compute
% the "center of mass" by weighting each pixel with its value in
% the original input array.  The result is not necessarily the
% actual peak!  You may want to apply a "hill climbing" algorithm
% to the returned coordinates to get the actual peaks.
%	   
% If you are interested in local maxima (which are not necessarily
% absolute maxima) then apply a high-pass filter to the input
% before using this routine.
%
% 2002-06-24 Tobin Fricke <tobin@splorg.org> Weizmann Institute of Science
% $Id: peak_find_cm.m,v 1.1 2002/06/24 13:18:36 tobin Exp tobin $

% Perform the thresholding

input_t = (input > threshold);

% Identify connected components

% use Hoshen-Kopelman Algorithm for Cluster Identification  (MATLAB "bwlabel")
% see http://splorg.org/people/tobin/projects/hoshenkopelman/
% TODO: Implement parallelized version (http://www.cs.utk.edu/~berry/parhk/)
 
[input_labelled, n_clusters] = bwlabel(input_t,4);
[I,J,V_clusters] = find(input_labelled);  

% Now I is a list of X coordinates of the pixels that are in the connected components
% and J is a list of Y coords.  V_clusters tells us what cluster that pixel belongs to.

% Extract the actual height of each pixel into V_heights

% See [ http://www.mathworks.com/company/digest/sept01/matrix.shtml ]

indicies = sub2ind(size(input),I,J);
V_heights = input(indicies);

% Now weight the (X,Y) coordinates by the height of the corresponding pixel

X_weights = V_heights .* I;
Y_weights = V_heights .* J;

% Now we need to add up the weighted coordinates for each cluster, and divide
% by the total weight of that cluster.

X = [];
Y = [];
weight = [];

for cluster = 1:n_clusters,
  weight(cluster) = sum(V_heights(V_clusters == cluster));
  X(cluster) = sum(X_weights(V_clusters == cluster))/weight(cluster);
  Y(cluster) = sum(Y_weights(V_clusters == cluster))/weight(cluster);
end

% If this works it will be amazing.


