Last year, I had the opportunity to visit Tate Britain which at the time hosted The EY Exhibition: Van Gough and Britain. It was a creative tour de force capturing Vincent Van Gough’s ability to transform his agony and pain into ecstatic beauty. Embellishing one of the otherwise unostentatious walls of the exhibit was his masterpiece, Starry Night Over the Rhône, which shows a couple walking near the quay by the river Rhône. The unique depiction of both the celestial and terrestrial gas lamps, and their reflection is stunning. What stands out is his expressive *impasto* (see Wikipedia), a stark contrast to prior Impressionist and Renaissance styles. A picture I took during my visit of this painting is shown below.

This painting motivates this blog post, the intent of which is to demonstrate the predictive power of higher order polynomials for interpolation. I am not aware of papers that focus on image-based interpolation with polynomials, perhaps owing to a large volume of literature using other, arguably more efficient methods. That said, the objective of this post is not to assert that polynomials should be used for this task. It is simply an ode to two things I like: orthogonal polynomials and paintings by Van Gough. Now lets get to the math!

**Pixels and images**

First some image preliminaries. A standard *greyscale* image is a matrix with entries of integers between 0 and 255 (inclusive), where each entry is referred to as a pixel. A pixel value of 0 corresponds to one that is colored black, whilst a pixel value of 255 corresponds to a white one.

*RGB* images are similar, but they have three channels (compared to the one for grayscale image) for the colors red, green and blue as shown below. Here too, pixels values are restricted between 0 and 255. To approximate an RGB image, one typically repeats the same calculation thrice—one for each channel.

**Polynomial Approximation**

We now focus on polynomial approximation. Weierstrauss’s approximation theorem tells us that polynomials can approximate every continuous function on a bounded interval to an arbitrary accuracy (see Chapter 8 in Trefethen). This even includes extremely non-smooth functions. One can think of an image as a discretized version of such a non-smooth function (although it is strictly speaking incorrect to think of an image as being a continuous function).

Rather than attempt to approximate the entire image with a single multivariate (2D) polynomial, we construct k univariate polynomials that interpolate the pixel values within a specific set of *rows groups*. In the image below, rows are divided into k=24 row groups, and within each row group vectorize the pixel values into a single array of 3100 elements.

This approach can be thought of as a very specific form of piecewise polynomial regression, albeit where the support of each polynomial is manually set a priori. If we restrict ourselves to univariate Legendre polynomials with an order of 200, we get the following results. The image below compares the true pixel values over the last row group; the image below that shows the resulting image formed via approximation over all the row groups.

Clearly, this approximation is far from perfect, but its a good start. We improve over this by increasing the number of quadrature points and reducing the number of row groups. For the results below we set k=8 and the order to be 2000.

Better, but not quite there yet. We now set k=48 such that each row group comprises of 5 rows with 1500 pixel entries in each. We set the quadrature order to be 700. This yields the results below.

Specific results for the last row group are also shown below.

So what next? I expect that inspection of the polynomial coefficients within each row group will reveal some information on the nature of the brushstrokes, but I will leave this for another post. Codes for replicating the above results can be found below.