A decorrelation stretch is a linear, pixel-wise operation in
which the specific parameters depend on the values of actual and desired
(target) image statistics. The vector `a` containing
the value of a given pixel in each band of the input image `A` is
transformed into the corresponding pixel `b` in output
image `B` as follows:

`b = T * (a - m) + m_target`.

`a` and `b` are `nBands`-by-1
vectors, `T` is an `nBands`-by-`nBands` matrix,
and `m` and `m_target` are `nBands`-by-1
vectors such that

`m` contains the mean of each band
in the image, or in a subset of image pixels that you specify

`m_target` contains the desired output
mean in each band. The default choice is `m_target = m`.

The linear transformation matrix `T` depends
on the following:

The band-to-band sample covariance of the image, or
of a subset of the image that you specify (the same subset as used
for `m`), represented by matrix `Cov`

A desired output standard deviation in each band.
This is conveniently represented by a diagonal matrix, `SIGMA_target`.
The default choice is `SIGMA_target = SIGMA`, where `SIGMA` is
the diagonal matrix containing the sample standard deviation of each
band. `SIGMA` should be computed from the same pixels
that were used for `m` and `Cov`,
which means simply that:

`SIGMA(k,k) = sqrt(Cov(k,k), k = 1,..., nBands`.

`Cov`, `SIGMA`, and `SIGMA_target` are `nBands`-by-`nBands`,
as are the matrices `Corr`, `LAMBDA`,
and `V`, defined below.

The first step in computing `T` is to perform
an eigen-decomposition of either the covariance matrix `Cov` or
the correlation matrix

`Corr = inv(SIGMA) * Cov * inv(SIGMA)`.

In the correlation-based method, `Corr` is
decomposed: `Corr = V LAMBDA V'`.

In the covariance-based method, `Cov` is
decomposed: `Cov = V LAMBDA V'`.

`LAMBDA` is a diagonal matrix of eigenvalues
and `V` is the orthogonal matrix that transforms
either `Corr` or `Cov` to `LAMBDA`.

The next step is to compute a stretch factor for each band,
which is the inverse square root of the corresponding eigenvalue.
It is convenient to define a diagonal matrix `S` containing
the stretch factors, such that:

`S(k,k) = 1 / sqrt(LAMBDA(k,k))`.

Finally, matrix `T` is computed from either

`T = SIGMA_target V S V' inv(SIGMA) ` (correlation-based
method)

or

`T = SIGMA_target V S V'` (covariance-based
method).

The two methods yield identical results if the band variances
are uniform.

Substituting `T` into the expression for `b`:

`b = m_target + SIGMA_target V S V' inv(SIGMA) * (a
- m)`

or

`b = m_target + SIGMA_target V S V' * (a - m)`

and reading from right to left, you can see that the decorrelation
stretch:

Removes a mean from each band

Normalizes each band by its standard deviation (correlation-based
method only)

Rotates the bands into the eigenspace of `Corr` or `Cov`

Applies a stretch `S` in the eigenspace,
leaving the image decorrelated and normalized in the eigenspace

Rotates back to the original band-space, where the
bands remain decorrelated and normalized

Rescales each band according to `SIGMA_target`

Restores a mean in each band.