Characterizing real-life optical surfaces usually involves finding the best-fit of an appropriate surface model to a set of discrete measurement data. This process can be greatly simplified by choosing orthonormal polynomials for the surface description. In case of rotationally symmetric aspherical surfaces, new sets of orthogonal polynomials were introduced by Forbes to replace the numerical unstable standard description. From these, for the application of surface retrieval using experimental ray tracing, the sag orthogonal Q(con)-polynomials are of particular interest. However, these are by definition orthogonal over continuous data and may not be orthogonal for discrete data. In this case, the simplified solution is not valid. Hence, a Gram-Schmidt orthonormalization of these polynomials over the discrete data set is proposed to solve this problem. The resulting difference will be presented by a performance analysis and comparison to the direct matrix inversion method.