In many image processing applications, the discrete values of an image can be embedded in a continuous function. This type of representation can be useful for interpolation, geometrical transformations or special features extraction. Given a rectangular M multiplied by N discrete image (or sub-image), it is shown how to compute a continuous polynomial function that guarantees an exact fit at the considered pixel locations. The polynomials coefficients can be expressed as a linear one-to-one separable transform of the pixels. The transform matrices can be computed using a fast recursive algorithm which enables efficient inversion of a Vandermonde matrix. It is also shown that the least square polynomial approximation with M' multiplied by N' coefficients, in the separable formulation, involves the inversion of two M' multiplied by M' and N' multiplied by N' Hankel matrices.