Time-resolved diffuse scattering experiments have gained increasing attention due to their potential to reveal non-equilibrium dynamics of crystal lattice vibrations with full momentum resolution. Although progress has been made in interpreting experimental data on the basis of one-phonon scattering, understanding the role of individual phonons can be sometimes hindered by multi-phonon excitations. In Ref. [{\it arXiv:2103.10108}] we have introduced a rigorous approach for the calculation of the all-phonon inelastic scattering intensity of solids from first-principles. In the present work, we describe our implementation in detail and show that multi-phonon interactions are captured efficiently by exploiting translational and time-reversal symmetries of the crystal. We demonstrate its predictive power by calculating the scattering patterns of monolayer molybdenum disulfide (MoS$_2$), bulk MoS$_2$, and black phosphorus (bP), and we obtain excellent agreement with our measurements of thermal electron diffuse scattering. Remarkably, our results show that multi-phonon excitations dominate in bP across multiple Brillouin zones, while in MoS$_2$ they play a less pronounced role. We expand our analysis for each system and examine the effect of individual atomic and interatomic vibrational motion on the diffuse scattering signals. We further demonstrate the high-throughput capability of our approach by reporting all-phonon scattering maps of 2D MoSe2, WSe2, WS2, graphene, and CdI2, rationalizing in each case the effect of multi-phonon processes. As a side point, we show that the special displacement method reproduces the thermally distorted configuration that generates precisely the all-phonon diffuse pattern.