Three-dimensional reconstruction of the corneal surface provides a powerful tool for managing corneal diseases. This study proposes a novel method for reconstructing the corneal surface from elevation point clouds, using modal schemes capable of reproducing corneal shapes using surface polynomial functions. The multivariable polynomial fitting was performed using a non-dominated sorting multivariable genetic algorithm (NS-MVGA). Standard reconstruction methods using least-squares discrete fitting (LSQ) and sequential quadratic programming (SQP) were compared with the evolutionary algorithm-based approach. The study included 270 corneal surfaces of 135 eyes of 102 patients (ages 11–63) sorted in two groups: control (66 eyes of 33 patients) and keratoconus (KC) (69 eyes of 69 patients). Tomographic information (Sirius, Costruzione Strumenti Oftalmici, Italy) was processed using Matlab. The goodness of fit for each method was evaluated using mean squared error (MSE), measured at the same nodes where the elevation data were collected. Polynomial fitting based on NS-MVGA improves MSE values by 86% compared to LSQ-based methods in healthy patients. Moreover, this new method improves aberrated surface reconstruction by an average value of 56% if compared with LSQ-based methods in keratoconus patients. Finally, significant improvements were also found in morpho-geometric parameters, such as asphericity and corneal curvature radii.