Time-harmonic scattering in variable media can often be modeled by linear elliptic partial differential equations. Such equations pose several numerical challenges. For example, they can be intrinsically ill-conditioned, necessitate the imposition of radiation conditions, and produce pollution errors when discretized with standard finite difference or finite element methods. To avoid these issues, it is often convenient to first reformulate the differential equation as an integral equation. The tradeoffs are that an integral operator with a singular kernel must be discretized and that the resulting linear system that must be inverted is dense. Sometimes, the latter issue can be handled using a fast matrix-vector product algorithm (e.g., the fast Fourier transform or the fast multipole method) paired with an iterative solver (e.g., GMRES). However, this approach can be prohibitively slow when there is a large amount of backscattering and when multiple incident fields are of interest. In these cases, it is better to use direct solvers. In this talk, I describe some recent projects on developing direct solvers in this regime with applications to acoustic scattering and metasurface design.