The irregular shapes of small solar system bodies are modelled by lognormal statistics, i.e., assuming that the shapes are realisations of the so-called Gaussian random sphere. The Gaussian sphere is fully described by the mean radius and the covariance function of the logarithmic radius. The stochastic shape is thus given by the covariance function, or the discrete spectrum of its Legendre coefficients. A maximum likelihood estimator is here provided for inverting the covariance function from three-dimensional sample shapes. The inverse method is applied to sophisticated shape data on altogether 14 small solar system bodies: the asteroids 4 Vesta, 243 Ida, 951 Gaspra, 1620 Geographos, 4179 Toutatis, and 4769 Castalia; the Martian satellites Phobos and Deimos; the Jovian satellite Amalthea; the Saturnian satellites Hyperion, Epimetheus, Janus, and Prometheus; and the Neptunian satellite Proteus. Inversion yields sigma = 0.245 for the relative standard deviation of radius, shows that most of the spectral power lies in the second-degree spherical harmonics, and gives Gamma = 32.7 degrees for the correlation angle. Even though the first results are promising, caution is recommended because the number of sample shapes is still small. Omitting one sample shape at a time and repeating the inversion shows that the results are not too sensitive to any one sample shape. As an example application, thermal light curves are simulated for 1000 Gaussian sample spheres in order to study the uncertainties in diameters and masses derived for asteroids. As compared to the Standard Thermal Model that assumes spherical asteroids, the irregular shape is shown to cause a 5 % systematic effect with 10 % scatter in diameter estimation whereas, in mass estimation, the respective numbers are larger at 17 % and 33 %.
