We propose a diagrammatic Monte Carlo approach for quantum impurity models, which can be regarded as a generalization of the strong-coupling expansion for fermionic impurity models. The algorithm is based on a self-consistently computed three-point vertex and a stochastically sampled four-point vertex, and it allows one to obtain numerically exact results in a wide parameter regime. The performance of the algorithm is demonstrated with applications to a spin-boson model representing an emitter in a waveguide. As a function of the coupling strength, the spin exhibits a delocalization-localization crossover at low temperatures, signaling a qualitative change in the real-time relaxation. In certain parameter regimes, the response functions of the emitter coupled to the electromagnetic continuum can be described by an effective Rabi model with appropriately defined parameters. We also discuss the spatial distribution of the photon density around the emitter.