Nanoparticles, exhibiting functionally relevant structural heterogeneity, are at the forefront of cutting-edge research. Now, high-throughput single-particle imaging (SPI) with x-ray free-electron lasers (XFELs) creates unprecedented opportunities for recovering the shape distributions of millions of particles that exhibit functionally relevant structural heterogeneity. To realize this potential, three challenges have to be overcome: (1) simultaneous parametrization of structural variability in real and reciprocal spaces; (2) efficiently inferring the latent parameters of each SPI measurement; (3) scaling up comparisons between $10^5$ structural models and $10^6$ XFEL-SPI measurements. Here, we describe how we overcame these three challenges to resolve the non-equilibrium shape distributions within millions of gold nanoparticles imaged at the European XFEL. These shape distributions allowed us to quantify the degree of asymmetry in these particles, discover a relatively stable `shape envelope' amongst nanoparticles, discern finite-size effects related to shape-controlling surfactants, and extrapolate nanoparticles' shapes to their idealized thermodynamic limit. Ultimately, these demonstrations show that XFEL SPI can help transform nanoparticle shape characterization from anecdotally interesting to statistically meaningful.