A new data processing method based on orthonormal rectangular gradient polynomials is introduced in this work. This methodology is capable of effectively reconstructing surfaces or wavefronts with data obtained from deflectometry systems, especially during fabrication and metrology of high resolution and freeform surfaces. First, we derived a complete and computationally efficient vector polynomial set, called G polynomials. These polynomials are obtained from gradients of Chebyshev polynomials of the first kind-a basis set with many qualities that are useful for modal fitting. In our approach both the scalar and vector polynomials, that are defined and manipulated easily, have a straightforward relationship due to which the polynomial coefficients of both sets are the same. This makes conversion between the two sets highly convenient. Another powerful attribute of this technique is the ability to quickly generate a very large number of polynomial terms, with high numerical efficiency. Since tens of thousands of polynomials can be generated, mid-To-high spatial frequencies of surfaces can be reconstructed from high-resolution metrology data. We will establish the strengths of our approach with examples involving simulations as well as real metrology data from the Daniel K. Inouye Solar Telescope (DKIST) primary mirror.