Let Delta(alpha,Y) be the bounded from above self-adjoint realization in L-2(R-3) of the Laplacian with n point scatterers placed at Y={y(1),...,y(n)} subset of R-3, the parameters (alpha(1),...alpha(n)) equivalent to alpha is an element of R-n being related to the scattering properties of the obstacles. Let u(f epsilon)(alpha,Y) and u(f epsilon)(empty set) denote the solutions of the wave equations corresponding to Delta(alpha,Y) and to the free Laplacian Delta respectively, with a source term given a pulse f(epsilon) supported in epsilon-neighborhoods of the points in X-N={x(1),...,x(N)}, X-N boolean AND Y= empty set. We show that, for any fixed gimel > sup sigma(Delta(alpha,Y)) >= 0, there exists N-circle >= 1such that the locations of the points in Y can be determined by the knowledge of a finite-dimensional scattering data operator F-lambda(N) : R-N -> R-N, N >= N-circle. Such an operator is defined in terms of the limit as epsilon SE arrow 0 of the Laplace transform of u(f epsilon)(alpha,f) (t, x(k)) - u(f epsilon)(empty set) (t, x(k)), k = 1,...,N. We exploit the factorized form of the resolvent difference (-Delta(alpha,Y)+lambda)(-1) - (-Delta+lambda)(-1) and a variation on the finite-dimensional factorization in the MUSIC algorithm. Multiple scattering effects are not neglected; our model can be interpreted as the time-domain version of a frequency-domain scattering from an array of Foldy's point-like obstacles. (c) 2022 Elsevier Inc. All rights reserved.