function [ij_max, V] = find_maxima_dep(im, inten_range)
	% [ij_max, f0] = find_maxima_dep(im)
	%
	% find the location and intensity of local maxima in the 
	% 3D image array im which are brighter than a certain threshold
	% 
	% Local maxima are maxima over the 2 nearest neighbour voxels in depth
	% direction
	%
	% Input:
	%	im				Ni x Nj x Nk	image stack
	% Output:
	%	ij_max		3 x N				indices of local maxima
	%	f0				1 x N				values of local maxima
	
	%% get size
	sz			= size(im);
	assert(all(sz >= 3));
	
	imin		= min(inten_range);
	if length(inten_range) < 2
		imax  = inf;
	else
		imax  = max(inten_range);
	end
	
	%% get local maxima
	middle	= im(2:end-1,2:end-1,2:end-1);
	lap		= del2(middle);
	is_max	= middle > im(2:end-1,2:end-1,1:end-2) ...
				& middle > im(2:end-1,2:end-1,3:end  ) ...
				& middle > imin ...
				& middle < imax ...
				& lap    < 0;
	[i,jk]	= find(is_max);
	[j,k]		= ind2sub(sz(2:3)-2, jk); % remember edges have been cut off...
	
	%% add offset to account for cutting off edges
	ij_max	= [i(:) j(:) k(:)]' + 1;
	V			= middle(sub2ind(sz-2, i, j, k))';
end