Estimating Diracs in continuous two or higher dimensions is a fundamental problem in imaging. Previous approaches extended one dimensional methods, like the ones based on finite rate of innovation (FRI) sampling, in a separable manner, e.g., along the horizontal and vertical dimensions separately in 2D. The separate estimation leads to a sample complexity of O(K^D) for K Diracs in D dimensions, despite that the total degrees of freedom only increase linearly with respect to D. We propose a new method that enforces the continuous-domain sparsity constraints simultaneously along all dimensions, leading to a reconstruction algorithm with linear sample complexity O(K), or a gain of O(K^{D-1}) over previous FRI-based methods. The multi-dimensional Dirac locations are subsequently determined by the intersections of hypersurfaces (e.g., curves in 2D), which can be computed algebraically from the common roots of polynomials. We first demonstrate the performance of the new multi-dimensional algorithm on simulated data: multi-dimensional Dirac location retrieval under noisy measurements. Then we show results on real data: radio astronomy point source reconstruction (from LOFAR telescope measurements) and the direction of arrival estimation of acoustic signals (using Pyramic microphone arrays).