function [ X, Y ]  = peak_find_max(input, threshold)
% [ X, Y ] = peak_find_max( 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 for each connected component we simply
% find the absolute maximum value within that component.
%
% 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_max.m,v 1.1 2002/06/25 14:49:52 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);

for cluster = 1:n_clusters,
  % OBS! This will only work is V_height is guaranteed to be greater than zero
  [z, l] = max(V_heights .* (V_clusters == cluster));
  
  X(cluster) = I(l);
  Y(cluster) = J(l);
  
end