radius_range	= [20 40];%30 65];
edge_threshold	= 50;
edge_threshold_mat=0.01;
n_circles		= 250;
im_size			= [2560 1600];

addpath([pwd '/../']);

%% create an image with circles in it
fprintf('making synthetic image...\n');
u_axis			= 0 : im_size(1)-1;
v_axis			= 0 : im_size(2)-1;
uv_ref			= diag(im_size) * rand(2, n_circles);
r_ref			= rand(1, n_circles)*(radius_range*[-1;1]) + radius_range(1);
[um,vm]			= ndgrid(u_axis, v_axis);
im				= ones(im_size);
sigma_edge		= 0.5;

for n = 1 : n_circles
	dist		= sqrt((um - uv_ref(1,n)).^2 + (vm - uv_ref(2,n)).^2);
	im_shad		= (erf((dist - r_ref(n)) / sigma_edge) + 1)/2;
	im			= im .* im_shad;
	%b_inside	= dist < r_ref(n);
	%im(b_inside)= 1;
end

%G_blur			= fspecial('gaussian', [5 5]);
%im				= convn(im, G_blur, 'same');
%im				= 1 - im;

%% find centers with MATLAB
tic;
opts			= struct(	'circle_radius', radius_range, ...
							'circle_threshold', edge_threshold_mat, ...
							'circle_polarity', 'dark', ...
							'circle_binning', 1);
[ctr_mat, radius_mat]= find_circles_matlab(u_axis, v_axis, im, opts); 
toc;

%% find centers with openCV
tic;
opts			= struct(	'circle_radius', radius_range, ...
							'circle_threshold', edge_threshold);
particles		= find_circles(u_axis, v_axis, im, opts); 
ctr				= [particles.centroid];
radius			= [particles.radius];
toc;

%% evaluate 
[i_mat,d_mat]	= knnsearch(ctr_mat', uv_ref');
[i_ocv,d_ocv]	= knnsearch(ctr', uv_ref');
% residuals
uv_res_mat		= ctr_mat(:,i_mat) - uv_ref;
uv_res_ocv		= ctr(:,i_ocv) - uv_ref;

d_thresh		= 1;
n_mat			= nnz(d_mat < d_thresh);
n_ocv			= nnz(d_ocv < d_thresh);
fprintf('MATLAB\n');
fprintf('   found %d circles within %0.1f px\n', n_mat, d_thresh);
fprintf('   rms error: %0.3f %0.3f px\n', rms(uv_res_mat, 2));
fprintf('   med error: %0.3f %0.3f px\n', median(abs(uv_res_mat), 2));

fprintf('OCV:\n');
fprintf('   found %d circles within %0.1f px\n', n_ocv, d_thresh);
fprintf('   rms error: %0.3f %0.3f px\n', rms(uv_res_ocv, 2));
fprintf('   med error: %0.3f %0.3f px\n', median(abs(uv_res_ocv), 2));

%% plot
figure(1);
clf;
imagesc(u_axis, v_axis, im');
hold on;
quiver(uv_ref(1,:), uv_ref(2,:), uv_res_ocv(1,:), uv_res_ocv(2,:), 'r');

axis equal tight ij;
colormap gray;
viscircles(ctr_mat', radius_mat', 'color', 'b');
viscircles(ctr', radius', 'color', 'r');