function [im_pos, f0, im_samp] = find_particles_3d(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 particles 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 six nearest neighbours
	% 
	% 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
	%	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
	%	samp			3x3x3 x Np		sample of particle image
	
	im_size					= size(im);
	
	%% identify local maxima
	[ij_part, f0]			= find_maxima_3d(im, inten_range);
	n_part					= size(ij_part, 2);
	
	%% sub-pixel peak location
	fi				= [im(sub2ind(im_size, ij_part(1,:)-1, ij_part(2,:)  , ij_part(3,:)  ));
						im(sub2ind(im_size, ij_part(1,:)+1, ij_part(2,:)  , ij_part(3,:)  ))];
	fj				= [im(sub2ind(im_size, ij_part(1,:)  , ij_part(2,:)-1, ij_part(3,:)  ));
						im(sub2ind(im_size, ij_part(1,:)  , ij_part(2,:)+1, ij_part(3,:)  ))];
	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_fi			= log(max(fi, 1));
	ln_fj			= log(max(fj, 1));
	ln_fk			= log(max(fk, 1));
	ln_f0			= log(max(f0, 1));
	
	% subpixel shift
	di				= 1/2 * (ln_fi(1,:) - ln_fi(2,:)) ./ (ln_fi(1,:) - 2*ln_f0 + ln_fi(2,:));
	dj				= 1/2 * (ln_fj(1,:) - ln_fj(2,:)) ./ (ln_fj(1,:) - 2*ln_f0 + ln_fj(2,:));
	dk				= 1/2 * (ln_fk(1,:) - ln_fk(2,:)) ./ (ln_fk(1,:) - 2*ln_f0 + ln_fk(2,:));
	
	%% sampling of image
	[irel,jrel,krel]	= ndgrid(-1:1, -1:1, -1:1);
	iabs					= bsxfun(@plus, irel, shiftdim(ij_part(1,:),-2));
	jabs					= bsxfun(@plus, jrel, shiftdim(ij_part(2,:),-2));
	kabs					= bsxfun(@plus, krel, shiftdim(ij_part(3,:),-2));
	im_samp				= im(sub2ind(im_size, iabs, jabs, kabs));
	im_samp				= reshape(im_samp, [3 3 3 n_part]);
	
	%% 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,:) + di - 1) * im_di + i_axis(1); 
						(ij_part(2,:) + dj - 1) * im_dj + j_axis(1);
						(ij_part(3,:) + dk - 1) * im_dk + k_axis(1)];	
end