function [im_pos, f0] = find_particles_dep(i_axis, j_axis, k_axis, im, inten_range)
	% im_pos, bright, fit = find_particles_3d(i_axis, j_axis, k_axis, im, inten_range)
	%
	% find locations of local depth-maxima in a 3D image stack to sub-pixel accuracy
	% using 3-point Gaussian fit
	%
	% particles are defined by:
	%   pixel being above threshold
	%   pixel being a local maximum over its two nearest neighbours in depth
	% 
	% Input:
	%	i_axis		1 x Ni			i axis coordinate of image
	%	j_axis		1 x Nj			j axis coordinate of image
	%	k_axis		1 x Nk			k axis coordinate of image (subpx accuracy)
	%	im				Ni x Nj x Nk	image stack
	%	inten_range	scalar/vec		brightness threshold for particle
	%										if scalar, this is a lower bound on bright
	%										if 2-vector, these are lower + upper bounds
	%
	% Output:
	%	im_pos		3 x Np			coordinates of particles in image space
	%	bright		1 x Np			brightness of particles
	
	im_size					= size(im);
	
	%% identify local maxima in depth direction
	[ij_part, f0]			= find_maxima_dep(im, inten_range);
	n_part					= size(ij_part, 2);
	
	%% sub-pixel peak location
	fk				= [im(sub2ind(im_size, ij_part(1,:)  , ij_part(2,:)  , ij_part(3,:)-1));
						im(sub2ind(im_size, ij_part(1,:)  , ij_part(2,:)  , ij_part(3,:)+1))];
	
	% logarithm of brightness
	ln_fk			= log(max(fk, eps));
	ln_f0			= log(f0);
	
	% subpixel shift
	dk				= 1/2 * (ln_fk(1,:) - ln_fk(2,:)) ./ (ln_fk(1,:) - 2*ln_f0 + ln_fk(2,:));
	
	%% translate to image coords
	im_di			= (i_axis(end) - i_axis(1)) / (im_size(1) - 1);
	im_dj			= (j_axis(end) - j_axis(1)) / (im_size(2) - 1);
	im_dk			= (k_axis(end) - k_axis(1)) / (im_size(3) - 1);
	im_pos		= [(ij_part(1,:) + 0  - 1) * im_di + i_axis(1); 
						(ij_part(2,:) + 0  - 1) * im_dj + j_axis(1);
						(ij_part(3,:) + dk - 1) * im_dk + k_axis(1)];	
end