How does DCT decorrelate images?
[EDIT] In 1991, Nasir Ahmed wrote: &quot;How I Came Up with the Discrete Cosine Transform&quot;. Interesting to read, on how he was inspired by  Chebyshev polynomials, and on how he didn't get funding, for a tool at the heart of JPEG and MP3.
Natural images are not very stationary, but locally, their covariance is often modeled by a first- or second-order process. Being able to &quot;concentrate&quot; those matrices is one key to efficient compression (not the only one, however).
In different words: one aims at projecting, or representing, the pixels in another system of vectors (a basis, a frame), such that the information content of the picture is better represented, or sparser. Classically, one preserves energy (hence orthogonal or close to orthogonal transformations. Two main options:

compute or learn the best adaptive representation from each image. PRO: more efficient in compaction; CON: computationally expensive and storage lost because should be sent for decoding
use a known and fixed transform (or a set of known transforms). PRO: fast; CON: less optimal.

The DCT is option 2, with the performance of option 1: a fixed transform with close to adaptive performance, for most standard images. It worked, people still use JPEG and MP3, at the heart of which one finds the DCT.
The goal is thus to find the eigenvectors of such (covariance matrices, with less computational cost than matrix diagonalisation. The original DCT (DCT-II) from Discrete Cosine Transform (1974) aimed at finding basis vectors that could approximately diagonalize the first-order Markov Toeplitz  matrix:
$$ \left[\begin{array}{c}
1&amp;\rho&amp; \cdots &amp;\rho^{N-1}\\
\rho&amp; 1&amp;\cdots &amp;\rho^{N-2}\\
\vdots&amp; \vdots&amp;\ddots &amp;\vdots\\
\rho^{N-1}&amp; \rho^{N-2}&amp;\ddots &amp;1\end{array}\right] $$

The original motivation for defining the DCT was that its basis set
provided a good approximation to the eigenvectors of the class of
Toeplitz matrices that constitutes the autocovariance matrix of a
first-order stationary Markov process, with the result that it had a
better performance than the discrete Fourier transform (DFT) and some
other transforms [l], 3, 4 with respect to such kinds of
processes. In fact, as shown in 4, the DCT is asymptotically
equivalent to the Karhunen-Loeve transform (KLT) of a first-order
stationary Markov process as $\rho$ tends to 1, where $\rho$ is the
correlation coefficient.

Looking at the basis functions for $\rho=0.9$

the match is pretty good, and it turns out that the processing performance is close as well to that of the Karhunen-Loève transform.
In their 1995 paper Diagonalizing Properties of the Discrete Cosine Transforms, the authors provide an extension of the above property.

In this paper, we obtain the eight types of DCT’s as the complete
orthonormal set of eigenvectors generated by a general form of
matrices in the same way as the discrete Fourier transform (DFT) can
be obtained as the eigenvectors of an arbitrary circulant matrix.
These matrices can be decomposed as the sum of a symmetric Toeplitz
matrix plus a Hankel or close to Hankel matrix scaled by some constant
factors.

In 1991, Nasir Ahmed wrote: &quot;How I Came Up with the Discrete Cosine Transform&quot;. Interesting to read, some excerpts:

What intrigued me was that the KLT was indeed the optimal transform on
the basis of the meansquare- error criterion and the first-order
Markov process model, and yet there was no efficient algorithm
available to compute it. As such, the focus of my research was to
determine whether it would be possible to come up with a good
approximation to the KLT that could be computed efficiently.


to study a “cosine transform” using Chebyshev polynomials


The motivation for looking into such “cosine functions”
was that they closely resembled KLT basis
functions for a range of values of the correlation coefficient
$\rho$ (in the covariance matrix).