function [part_uv, f0, sxx, syy, sxy, im_filt] = find_particles_2d(u_axis, v_axis, im, threshold, filter_scale, finder)
	% part_uv, part_bright = find_particles(u_axis, v_axis im, threshold)
	%
	% find locations of particles in the image im, which is a 2D array
	% with dimensions Ni x Nj
	%
	% particle locations are returned in the 2 x Np array part_ij
	% where Np is the number of particles that were found. The location is
	% determined to sub-pixel accuracy using Gaussian peak finding
	%
	% part_bright is a 1 x Np vector of the peak brightness of the particle image
	%
	% i_axis and j_axis are Ni- and Nj- vectors respectively
	% which define the coordinate system of the image
	%
	% particle images are identified by:
	%   pixel being above threshold
	%   pixel being a local maximum over its four nearest neighbours

	im_size					= size(im);
	
	if nargin < 6
		finder				= 3;
	end
	
	%% filter
	if filter_scale > 0
		hsize					= ceil(3 * filter_scale)*2 + 1;
		G						= fspecial('gaussian', hsize, filter_scale);
		im_filt				= conv2(im, G, 'same');
	else
		im_filt				= im;
	end
	[part_i, part_j, f0] = find_maxima(im_filt, threshold);
	
	%% exclude points near edges
	ok							= part_i > 2 & part_i < im_size(1)-1 ...
								& part_j > 2 & part_j < im_size(2)-1;
	part_i					= part_i(ok);
	part_j					= part_j(ok);
	f0							= f0(ok);
	n_part					= nnz(ok);
	
	%% stop if no particles
	if n_part < 1
		part_uv				= zeros(2,0);
		f0						= [];
		sxx					= [];
		syy					= [];
		sxy					= [];
		return;
	end
	
	%% get sub-pixel location
	im_ary					= zeros(5,5,n_part);
	for k = 1 : n_part
		im_ary(:,:,k)		= im_filt(part_i(k)+(-2:2), part_j(k)+(-2:2));
	end
	peaks						= peakfit(im_ary, finder);
	
	%% translate to image coords
	part_di					= peaks(1,:);
	part_dj					= peaks(2,:);
	f0							= peaks(3,:);
	sxx						= peaks(4,:);
	syy						= peaks(5,:);
	sxy						= peaks(6,:);
	
	du							= (u_axis(end) - u_axis(1)) / (im_size(1) - 1);
	dv							= (v_axis(end) - v_axis(1)) / (im_size(2) - 1);
	part_uv					= [(part_i' + part_di - 1) * du + u_axis(1); 
									(part_j' + part_dj - 1) * dv + v_axis(1)];	
end