This talk outlines a novel class of high-order methods for the efficient numerical evaluation of volume potentials (VPs) defined by volume integrals over complex geometries. Inspired by the Density Interpolation Method (DIM) for boundary integral operators, the proposed methodology leverages Green’s third identity and a local polynomial interpolation of the density function to recast a given VP as a linear combination of surface-layer potentials and a volume integral with a regularized (bounded or smoother) integrand. The layer potentials can be accurately and efficiently evaluated inside and outside the integration domain using existing methods (e.g. DIM), while the regularized volume integral can be accurately evaluated by applying elementary quadrature rules to integrate over structured or unstructured domain decompositions without local numerical treatment at and around the kernel singularity. The proposed methodology is flexible, easy to implement, and fully compatible with well-established fast algorithms such as the Fast Multipole Method and H-matrices, enabling VP evaluations to achieve linearithmic computational complexity. To demonstrate the merits of the proposed methodology, we applied it to the Nyström discretization of the Lippmann-Schwinger volume integral equation for frequency-domain Helmholtz scattering problems in piecewise-smooth variable media. This is joint work with Thomas G. Anderson (Rice University), Luiz Faria (ENSTA Paris), and Marc Bonnet (ENSTA Paris).