This paper presents a general high-order kernel regularization technique
applicable to all four integral operators of Calder\'on calculus associated
with linear elliptic PDEs in two and three spatial dimensions. Like previous
density interpolation methods, the proposed technique relies on interpolating
the density function around the kernel singularity in terms of solutions of the
underlying homogeneous PDE, so as to recast singular and nearly singular
integrals in terms of bounded (or more regular) integrands. We present here a
simple interpolation strategy which, unlike previous approaches, does not
entail explicit computation of high-order derivatives of the density function
along the surface. Furthermore, the proposed approach is kernel- and
dimension-independent in the sense that the sought density interpolant is
constructed as a linear combination of point-source fields, given by the same
Green's function used in the integral equation formulation, thus making the
procedure applicable, in principle, to any PDE with known Green's function. For
the sake of definiteness, we focus here on Nystr\"om methods for the (scalar)
Laplace and Helmholtz equations and the (vector) elastostatic and time-harmonic
elastodynamic equations. The method's accuracy, flexibility, efficiency, and
compatibility with fast solvers are demonstrated by means of a variety of
large-scale three-dimensional numerical examples