The application of Gram polynomial basis functions to hierarchical image processing is presented in this research work. In the early 1960s Ming-kuei Hu used continuous geometrical moments for pattern recognition; since then there has been much research work done on moment invariants in the field of image processing. An accurate set of basis function is required to reconstruct the image of larger order with minimum error. The use of Chebychev moments produces numerical instabilities for moments of large order. The polynomial basis functions used in the method proposed here are discrete orthogonal basis, being unary discrete polynomial basis of order n. These basis functions are numerically better conditioned than discrete cosine transform, which leads us to a new method of image compression. A new multiresolution image analysis technique is presented based on hierarchies of images. The structure of the hierarchy is adapted to the image information and artefacts of each sequential images are reduced by Gram polynomial decimation. A major improvement is achieved by implementing suitable amount of decimation at each level; this decimation is implemented via Gram polynomial bases. Both global and local polynomial approximation are considered and compared with the Fourier basis. The issue of Gibbs error in polynomial decimation is examined. It is shown that the Gram basis is superior when applied to signals with strong gradient, i.e., a gradient which generated a significant Gibbs error with Fourier basis. The modified functions are used to compute spectra whereby the Gibbs error associated with local gradients in the image are reduced. The present work in the field of image registration also presents the first direct linear solution to weighted tensor product polynomial approximation. This method is used to regularize the patch coordinates, the solution is equivalent to a Galerkin type solution to a partial differential equations. The new solution is applied to published standard data sets and to data acquired in a production environment. The speed of the new solution justifies explicit reference: the present solution, implemented in MATLAB, requires approximately 1.3s to register an image of size 800 x 500 pixels. This is approximately a factor 10 to 100 time faster than previously published results for the same data set. The proposed algorithm is applied to non-rigid elastic registration of hyper spectral imaging data for the automatic quality control of decorative foils.