Multiphoton microscopy (MPM) is a relatively recent tool involved in biological imaging. Its resolution is somewhat limited due to its physical principles, resulting experimentally in a deteriorated planar resolution of about few tenths of micrometers and a reduced axial resolution of a few micrometers. In this communication, we present a numerical approach taking the form of a processing pipeline for image restoration going from the raw images of 0.2 μm diameter fluorescent microbeads to the reconstructed ones. The strategy consists in estimating 3D Point-Spread Function (PSF) shapes in automatically selected volumes of interest, considering a 3D Gaussian profile of intensity. An automatic crop procedure selects areas of a few dozens of individual beads. Our algorithms, based on a variational approach for multivariate Gaussian fitting, allows us to identify evolutions of PSF dimensions along the 3 dimensions of the volume. A proximal alternating optimization method called FIGARO is employed to minimize a cost function serving to quantify the optimality of the process. The difficulty overcome thanks to this novel variational approach for multivariate Gaussian fitting lies in its ability to estimate parametrically Gaussian shapes in a 3D space, which leads to accurate models of the PSF. Thanks to our technique, the PSF measurements show a dimension in (X,Y,Z) of (0.21, 0.27, 1.49) μm. It is the first time to our knowledge that such a direct 3D measurement has been made possible. This process paves the way to a significant improvement of the resolution, perfectly suited to MPM.