function [picks_X, picks_Y, picks_angle, picks_height, picks_antiheight, ...
	 picks_segments, n_segments ] = ppcc(image, ref, result, ...
	 anti_result,  best_angles, cutoff, threshold, angle_inc, verbose)
% Particle Picking by Cross Correlation
%
%   function [X, Y, A, H, H', S, n_S] = ppcc(image, ref, result, anti_result,
%   best_angles, cutoff, threshold,
%   angle_inc, verbose)
%
% This works by magic.
%
% The following parameters are used:
%
%   image      the image to be searched
%   ref        the reference image for which we are searching
%   cutoff     the frequency cutoff for the high-pass filter
%   threshold  threshold applied to filtered cross-correlation image
%
% suggested parameter values: 
%  cutoff = 0.03, threshold = 0.60, angle_inc = 4
%
% Note: the process will be much faster if the dimensions of the
% input images (image,ref) are composed of small prime numbers.
% For example, a 500x500 image is much better than a 501x501 pixel
% image, because 500 = 5*5*5*2*2 while 501 = 3*167.  You can use
% the matlab command "factor(n)" to show the prime factor
% decomposition of any integer.  
%
% 2002-06-16 Tobin Fricke <tobin@splorg.org> Weizmann Institute of Science
% $Id: ppcc.m,v 1.10 2002/07/02 12:02:13 tobin Exp tobin $

% if verbose, disp('Performing cross correlations.'); end
%
% [result, anti_result, best_angles] = cc_rotations(image, ref, angle_inc); 
if verbose, disp('Scaling result.'); end

[result, anti_result] = rescale2(result, anti_result);

if verbose, disp('Filtering result.'); end
  
% do we take abs or real?
% we should take real, as the real parts can be negative, and abs()
% does weird things in that case

result_filtered = rescale(abs(apply_highpass_filter(result, cutoff)));  

if verbose, disp('Finding peaks.'); end

[ picks_X, picks_Y ] = peak_find_max(result_filtered, threshold);

n_picks = length(picks_X);

if verbose, disp(sprintf('%d picks have been picked', n_picks)); end
  
% Extract information from the picks: the best correlation value
% acheived (picks_height), the orientation angle of the reference
% that produced that best correlation (picks_angle), and the
% correlation value resulting from the angle 180 degrees opposite
% the best angle (picks_antiheight).

picks_angle = [];
picks_height = [];
picks_antiheight = [];

% TODO: we should interpolate, but instead we'll just use a
% nearby value

indicies = sub2ind(size(result), round(picks_X), round(picks_Y));

for n = 1:n_picks,
  picks_angle(n) = best_angles(indicies(n));
  picks_height(n) = result(indicies(n));
  picks_antiheight(n) = anti_result(indicies(n));
end

if verbose, disp('Grouping by segment.'); end

% Group picks by segment -- if two picks are within a certain
% distance, they are considered to be members of the same segment;
% this relationship is transitive, so we compute equivalence classes.

pairwise_dist = pairwise_distances(picks_X, picks_Y);
pairwise_dist_t = pairwise_dist < 70;
[picks_segments, n_segments] = eclazz(pairwise_dist_t);















