function particles = find_circles(u_axis, v_axis, im, opts)
	% find circles on image using hough transform
	%
	% Input
	%	u_axis	1 x M	u coordinate axis
	%	v_axis	1 x N	v coordinate axis
	%	im		M x N	image to search
	%	opts	struct	options for search
	%
	% Options structure
	%	radius		radius range to search for, pixels
	%	threshold	threshold for identifying edges [0-1]
	%	polarity	polarity of the circular objects, 'bright'/'dark'
	%	binning   	downsample image by this factor to
	%				accelerate search for circles
	%
	% Output
	%	uv_ctr	2 x P	coordinates of circle centres
	%	radius	1 x P	radii of detected circles
	
	assert(ismatrix(im));
	assert(isvector(u_axis));
	assert(isvector(v_axis));
	
	%% find objects in image using OpenCV hough circle transform
	r_search		= opts.radius;	% range of radii to search
	r_min			= min(r_search);	% min distance between circles
	% rescale image to binary grayscale uint8
	im				= single(im);
	imax			= max(im(:));
	imin			= min(im(:));
	im_uint8		= uint8((im - imin)/(imax-imin) * 256);
	% estimate threshold for number of accumulator votes
	vote_threshold	= round(min(r_search)*2*pi * 0.2);
	edge_threshold	= opts.threshold;
	[idx_ctr,radius]= find_circles_ocv(im_uint8', r_search, r_min, edge_threshold, vote_threshold);
	
	%% convert indices provided by imfindcircles to image coordinate system
	du			= u_axis(2)-u_axis(1);
	dv			= v_axis(2)-v_axis(1);
	u0			= u_axis(1);
	v0			= v_axis(1);
	uv_ctr		= diag([du dv]) * idx_ctr + [u0; v0];

	%% convert to structure array
	uv_ary		= num2cell(uv_ctr, 1);
	r_ary		= num2cell(radius, 1);
	particles	= struct('centroid', uv_ary, 'radius', r_ary);
end
					