DStretch is really PCA

Comments
Last updated 2020-12-16 20:47:50 SGT

I happened to chance across a post on the UPenn Language Log a while back that makes reference to an image enhancement algorithm called DStretch, apparently developed at JPL. Reading through the description of the technique, however, it's pretty clear to me that the method is basically PCA. Let's see how this works on the images supplied with the Language Log article.

import numpy as np
import matplotlib.pyplot as plt
from imageio import imread

The input image is assumed to be of shape [(M, N, C)], which, for the purposes of this exercise, we will reshape into an [MN] by [C] design matrix (i.e. assume we have [MN] samples of [C] random variables). I apply autoscaling to the input data, find the corresponding principal components in the usual manner (via SVD of the design matrix), and return the image with the variances normalised to unity.

def process(data):
    X = data.reshape(-1, 3).astype(float)
    μ = np.mean(X, axis=0)
    σ = np.std(X, axis=0)
    Z = (X - μ) / σ
    U, Λ, V = np.linalg.svd(Z, full_matrices=False)
    return U @ V, μ, σ

I also define some convenience functions for display purposes.

def rescale(x):
    return (x - np.min(x)) / (np.max(x) - np.min(x))

def shape(x, data):
    if len(x.shape) == 1:
        return x.reshape(data.shape[:-1])
    else:
        return x.reshape(data.shape)

Let's apply this to the two images supplied by the Language Log post, comparing the results side by side with the originals:

for image in ["sarada1.jpg", "sarada2.jpg"]:
    data = imread(image)

    I, μ, σ = process(data)

    f, ax = plt.subplots(1,2)
    plt.sca(ax[0])
    plt.imshow(shape(rescale(I), data))
    plt.sca(ax[1])
    plt.imshow(data)
    plt.gcf().set_size_inches(12, 12)
    plt.show()

Understandably, this looks far less visually appealing than the sample result shown in the LL post, but then this is a fairly naive implementation. Partly, this is because we are performing this decomposition, and reconstructing the displayed images, in different colourspaces than the original implementation. In my tests of this method on other images (i.e. my collection of desktop wallpapers), these reconstructed images tend to be highly saturated and vibrant. This doesn't necessarily translate into better perceptual contrast, which I surmise is the purpose of using these custom colourspaces as an additional layer of processing.

For quantitative purposes, I think it's also instructive to examine the results of individual colour channels of the processed images:

for image in ["sarada1.jpg", "sarada2.jpg"]:
    data = imread(image)

    I, μ, σ = process(data)

    f, ax = plt.subplots(1,2)
    plt.sca(ax[0])
    plt.imshow(shape(rescale(I[:,1]), data))
    plt.sca(ax[1])
    plt.imshow(data)
    plt.gcf().set_size_inches(12, 12)
    plt.show()

To my eye this representation does appear to yield basically the same improvement in contrast as the full-colour DStretch reconstruction shown in the LL post. It's possible that this might be because the Viridis colourmap is engineered for perceptual uniformity (thereby standing in for the custom output colourspaces used by DStretch), although I don't claim to know for certain.


comments powered by Disqus