  [Abridged] We present a novel technique, dubbed FiEstAS, to estimate the
underlying density field from a discrete set of sample points in an arbitrary
multidimensional space. FiEstAS assigns a volume to each point by means of a
binary tree. Density is then computed by integrating over an adaptive kernel.
As a first test, we construct several Monte Carlo realizations of a Hernquist
profile and recover the particle density in both real and phase space. At a
given point, Poisson noise causes the unsmoothed estimates to fluctuate by a
factor ~2 regardless of the number of particles. This spread can be reduced to
about 1 dex (~26 per cent) by our smoothing procedure. [...] We conclude that
our algorithm accurately measure the phase-space density up to the limit where
discreteness effects render the simulation itself unreliable. Computationally,
FiEstAS is orders of magnitude faster than the method based on Delaunay
tessellation that Arad et al. employed, making it practicable to recover
smoothed density estimates for sets of 10^9 points in 6 dimensions.
