function particles = find_template(u_axis, v_axis, im, kernel, threshold, mask)
	% particles = find_template(u, v, im, kernel)
	% 
	% find objects matching a template in an image using convolution
	%
	% Input
	%	u			1 x M		u coordinate of image
	%	v			1 x N		v coordinate of image
	%	im			M x N		image to search
	%	kernel		K x K		kernel for convolution
	%	threshold	1 x 1		threshold for identifying objects
	%				or 1 x 2	if scalar, threshold is lower bound for 
	%							intensity of filtered image. 
	%							if 2 element vector, threshold(1) is 
	%							lower bound, threshold(2) is upper bound	
	%	mask		M x N		(optional) image mask, logical array
	%							where mask==1, do not detect center
	%
	% Output
	%	particles	1 x P		structure array of 2d particles with
	%	centroid	2 x 1		(u,v) coordinates of matches
	%	intensity	1 x 2		amplitude of match
	
	assert(isrow(u_axis));
	assert(isrow(v_axis));
	du					= mean(diff(u_axis));
	dv					= mean(diff(v_axis));

	% if only lower bound
	if numel(threshold) == 1
		threshold		= [threshold inf];
	end
	if nargin < 6
		mask			= false(size(im));
	end

	%% convolve image with kernel, find maxima
	im_size				= size(im);
	im					= single(im);
	edgemask			= true(im_size);
	edgemask([1 end],:) = false;
	edgemask(:,[1 end]) = false;
	im_filt				= convn(im, kernel, 'same');
	% find localised maxima, within the radius of the kernel
	im_local_max		= maxfiltnd(single(im_filt), size(kernel));
	[i_ctr,j_ctr]		= find(im_filt == im_local_max & edgemask);
	A_00				= im_filt(sub2ind(im_size, i_ctr  , j_ctr  ));
	b_masked			=    mask(sub2ind(im_size, i_ctr  , j_ctr  ));
	% get brightness of neighbouring pixels to fit sub-pixel position
	A_m0				= im_filt(sub2ind(im_size, i_ctr-1, j_ctr  ));
	A_p0				= im_filt(sub2ind(im_size, i_ctr+1, j_ctr  ));
	A_0m				= im_filt(sub2ind(im_size, i_ctr  , j_ctr-1));
	A_0p				= im_filt(sub2ind(im_size, i_ctr  , j_ctr+1));
	% check if this is a strict maximum above threshold and not masked
	b_ok				= A_00 >= threshold(1) & A_00 <= threshold(2) ...
						& A_00 > A_m0 & A_00 > A_p0 ...
						& A_00 > A_0m & A_00 > A_0p ...
						& ~b_masked;
	% select only pixels which are strict maxima
	i_ctr				= i_ctr(b_ok,1).';
	j_ctr				= j_ctr(b_ok,1).';
	A_00				= A_00(b_ok,1).';
	A_m0				= A_m0(b_ok,1).';
	A_p0				= A_p0(b_ok,1).';
	A_0m				= A_0m(b_ok,1).';
	A_0p				= A_0p(b_ok,1).';

	%% find sub-pixel position using Gaussian peak finder
	fi					= [A_m0; A_00; A_p0];
	fj					= [A_0m; A_00; A_0p];
	% logarithm of brightness
	ln_fi				= log(max(fi, eps));
	ln_fj				= log(max(fj, eps));
	
	% subpixel shift
	part_du				= 1/2 * du * (ln_fi(1,:) - ln_fi(3,:)) ./ (ln_fi(1,:) - 2*ln_fi(2,:) + ln_fi(3,:));
	part_dv				= 1/2 * dv * (ln_fj(1,:) - ln_fj(3,:)) ./ (ln_fj(1,:) - 2*ln_fj(2,:) + ln_fj(3,:));
	% add sub-pixel shift
	uv_ctr0				= [u_axis(1,i_ctr); v_axis(1,j_ctr)];
	uv_ctr				= uv_ctr0 + [part_du; part_dv];

	%% convert to structure array
	uv_ary				= num2cell(uv_ctr, 1);
	A_00				= reshape(A_00, [1 numel(A_00)]);
	A_ary				= num2cell(A_00 , 1);
	particles			= struct('centroid', uv_ary, 'intensity', A_ary);
end