function [solutions, values] = h3d_evaluate_solutions(p, ht1, ht2, solutions)
	% function [solutions, values] = h3d_evaluate_solutions(p, ht1, ht2, solutions)
	% 
	% Evaluates the goodness of a list of solutions.
	%
	% - solutions is a list of tuples {rotation, translation}
	
	% Find good directions for the correlation
	fprintf('h3d_evaluate_solutions: Finding peaks on HT...\n');
	peaks = h3d_find_hs_peak(ht1);
	peaks = peaks(2:4,:);
	
	% Peaks that we use
	peaks_used = zeros(3,0);
	
	for i=1:size(peaks,2)
		peak = peaks(:,i);
		
		% Do not use this direction if it is too close to one that we already used.
		
		if closer_than_deg(peak, peaks_used, 5)
			fprintf('Discarding peak at [%1.2f %1.2f %1.2f] (too close)\n', peak(1),peak(2),peak(3));
			continue
		end
		
		% Do not use this direction if the other spectrum is very poor
			
%		[f1, i1, j1] = h3d_coords_s2_to_cell(peak, ht1.cube_ncells);
%		[f2, i2, j2] = h3d_coords_s2_to_cell(R_est * peak, ht2.cube_ncells);

%		if(ht2.cube_hs(f2,i2,j2) <= 0.01)
%			fprintf('Discarding peak at [%1.2f %1.2f %1.2f] (value of other hs is %lf)\n', peak(1),peak(2),peak(3), ht2.cube_hs(f2,i2,j2) );
%			continue;
%		end

		fprintf('Using peak at [%1.2f %1.2f %1.2f] \n', peak(1),peak(2),peak(3));
		peaks_used(1:3, size(peaks_used,2) + 1 ) = peak;
		
		if(size(peaks_used,2) >= p.guess_tran_max_directions)
			break;
		end
	end
	
	fprintf('Directions to use:\n')
	peaks_used

	show_corr = 0;

	values = [];
	for sol = 1:numel(solutions)
		progress(sol, numel(solutions));
		
		R_est = solutions{sol}{1};
		t_est = solutions{sol}{2};
		
		if show_corr
			figure;
		end
		
		value = 0;

		npeaks = size(peaks_used,2);
		for ipeak=1:npeaks
		
			peak = peaks_used(1:3,ipeak);
		
			[f1, i1, j1] = h3d_coords_s2_to_cell(peak, ht1.cube_ncells);
			[f2, i2, j2] = h3d_coords_s2_to_cell(R_est * peak, ht2.cube_ncells);
		
			s1 = squeeze(ht1.cube_ht(f1,i1,j1,:));
			s2 = squeeze(ht2.cube_ht(f2,i2,j2,:));
			
			lag_est  = ( (R_est * peak)' * t_est )          / ht1.rho_cell_size;
			lag_true = ( (R_est * peak)' * p.debug_true_T ) / ht1.rho_cell_size;

			pad = zeros(floor(lag_est),1);
			if lag_true > 0
				f1 = [pad; s1];
				f2 = [s2; pad];
			else
				f1 = [s1; pad];
				f2 = [pad; s2];				
			end
			
			value = value + f1' * f2;
			
			if show_corr
				subplot(npeaks,1,ipeak); hold on;
				plot(f1,'r.')
				plot(f2,'b.')
			end
		end
		
		values(sol) = value;
	end

	if false
	fprintf('Summary:\n');
	for i=1:numel(solutions)

		R = solutions{i}{1};
		T = solutions{i}{2};
		val = values(i);

		if numel(p.debug_true_T) == 3
			e_r = h3d_rotation_error(R, p.debug_true_R);
			e_t = (T -  p.debug_true_T);
			s = sprintf('error R: %03.1f deg; T: %3.1f cm', e_r, norm(e_t) * 100);
		else
			s= '';
		end
		fprintf('%d Val = %1.3f   R = %s  T = %s   %s\n', i, val, display_rot_aa(R), display_v3(T),s);	
	end
	end
	
function res = closer_than_deg(v, others, max_deg)
	% closer_than_deg(v, others, max_deg)
	%   v: 3x1 vector
	%   others: 3xn matrix
	%   returns 1 if v or -v is closer than max_deg to at least 
	%   one column of others
	
	n = size(others,2);
	
	for i=1:n
		if angle_between_vectors(v, others(:,i)) < deg2rad(max_deg)
			res = 1; return;
		end
	end
	
	res = 0;