I recently found myself squinting at two images at work and wondering if one is actually “better” than the other. The two images were showing the same object, but they had slightly but noticeably different value ranges. That made it hard for me to understand which of the differences were due to display brightness and contrast settings and which were part of the actual structure of the image. Then I remembered a simple technique that I’ve been using on and off since my PhD days.

The Problem

So let’s say we have two images \(\boldsymbol{I}\) and \(\boldsymbol{J}\) we want to normalize to a common brightness range. Despite showing the same structures, maybe we have generated the two images using slightly different processing techniques. The structures of interest might present a little differently, but the value ranges might also be shifted noticeably between the images. Using the minimum and maximum values to normalize the value ranges is often a bad idea because those properties can be very sensitive to noise. On the other hand, using the mean or mode to scale one image to the other is often good enough, but gives us only one degree of freedom. Note that we assume that the images are spatially aligned and the only thing we want to correct is a brightness transform.

A Simple Brightness Transform

Now let’s find a simple brightness transform that maps image \(\boldsymbol{J}\) to image \(\boldsymbol{I}\). Let’s allow a linear brightness transform with coefficients because such a transform is typically so constrained that it won’t alter the structures of interest. It will only help us shift one image to the value range of the other image to make it comparable. Every pixel \(J_k\) at index \(k\) of image \(\boldsymbol{J}\) is transformed as¹:

\[J_k \mapsto a J_k + b, \tag{1} \label{transform}\]

where \(a, b \in \mathbb{R}\) are the coefficients of the transformation. So, what are the best coefficients \(a, b\) such that the brightness variations between the images are minimized? Well, we can, e.g., minimize the pixelwise deviations of the images in a least squares sense²:

\[\min_{a,b} \frac{1}{N} \sum_{k=1}^{N} \left( a J_k +b - I_k \right)^2 \tag{2} \label{ols}\]

The neat thing is that there’s an analytical solution for the coefficients, which can be found by simply calculating the point at which the partial derivatives \(\partial/\partial a, \partial/\partial b\) of the expression both vanish. After some rewriting, this leads us to:

\[\begin{eqnarray} a &=& \frac{\text{cov}(\boldsymbol{I},\boldsymbol{J})}{\text{var}(\boldsymbol{J})} \label{a-ols} \tag{3a}\\ b &=& \text{mean}(\boldsymbol{I}) - a\cdot \text{mean}(\boldsymbol{J}) \label{b-ols} \tag{3b}, \end{eqnarray}\]

where \(\text{mean}\) calculates the mean value of an image, \(\text{cov}\) is the covariance, and \(\text{var}\) is the variance:

\[\begin{eqnarray} \text{mean}(\boldsymbol{I}) &=& \frac{1}{N} \sum_{k=1}^{N} I_k \tag{4a} \label{mean}\\ \text{cov}(\boldsymbol{I},\boldsymbol{J}) &=& \frac{1}{N-1} \sum_{k=1}^{N} (I_k-\text{mean}(\boldsymbol{I}))\cdot(J_k- \text{mean}(\boldsymbol{J})) \tag{4b} \label{cov}\\ \text{var}(\boldsymbol{J}) &=& \text{cov}(\boldsymbol{J},\boldsymbol{J}) \tag{4c} \label{var} \end{eqnarray}\]

That means using eq. \(\eqref{transform}\) with coefficients from \(\eqref{a-ols}\) and \(\eqref{b-ols}\) to adjust the brightness values in image \(\boldsymbol{J}\) will adjust the image to fit the brightness range of image \(\boldsymbol{I}\).

Weighted Least Squares

There’s one more improvement to the technique we can make quite easily, which is to use weighted least squares as the minimization objective rather than the ordinary least squares we used above. This allows us to give each pixel \(k\) a corresponding weight \(w_k\) that modifies its contribution to the objective function. For example, if our images contain large areas of near-zero background, we could force areas with higher signal to contribute more strongly by setting the weights as \(w_k = I_k\); we could also set the weights for pixels below a threshold to \(0\). There are many possible ways to improve the optimization with weighting.

For weighted least squares, our minimization objective becomes

\[\min_{a,b} \frac{1}{N} \sum_{k=1}^{N} w_k \left( a J_k +b - I_k \right)^2 \tag{5} \label{wls},\]

and we can find very similar analytical solutions for the coefficients

\[\begin{eqnarray} a &=& \frac{\text{cov}_{\boldsymbol{w}}(\boldsymbol{I},\boldsymbol{J})}{\text{var}_{\boldsymbol{w}}(\boldsymbol{J})} \label{a-wls} \tag{6a}\\ b &=& \text{mean}_{\boldsymbol{w}}(\boldsymbol{I}) - a\cdot \text{mean}_{\boldsymbol{w}}(\boldsymbol{J}) \label{b-wls} \tag{6b}. \end{eqnarray}\]

Now we have to use the weighted mean, covariance, and variance defined as

\[\begin{eqnarray} \text{mean}_{\boldsymbol{w}}(\boldsymbol{I}) &=& \frac{1}{\sum_{k=1}^N w_k} \sum_{k=1}^{N} w_k I_k \tag{7a} \label{mean-w}\\ \text{cov}_{\boldsymbol{w}}(\boldsymbol{I},\boldsymbol{J}) &=& \frac{N}{N-1}\frac{\sum_{k=1}^{N} w_k (I_k-\text{mean}_{\boldsymbol{w}}(\boldsymbol{I}))\cdot(J_k- \text{mean}_{\boldsymbol{w}}(\boldsymbol{J}))}{\sum_{k=1}^N w_k} \tag{7b} \label{cov-w}\\ \text{var}_{\boldsymbol{w}}(\boldsymbol{J}) &=& \text{cov}_{\boldsymbol{w}}(\boldsymbol{J},\boldsymbol{J}) \tag{7c} \label{var-w} \end{eqnarray}\]

Prior Art and Further Reading

Obviously, I wasn’t the first to come up with this least-squares-based linear brightness normalization technique. In remote sensing it’s known as radiometric normalization and described, e.g., by Zhang et al. ³. The relevant results are given in eqn. (1) and (2) of that paper. The authors also propose an iterative reweighting that goes beyond the scope of this article. For a much more sophisticated framework specific to this field, see Canty et al. ⁴

Outlook

At the beginning, I said that I wanted to normalize two images to a common brightness range, but what I’ve actually presented is a technique to transform one image into the brightness range of another image. How do we know which image to pick? This might seem pedantic for two images, but the problem becomes more obvious if we have a set of more than two images that we want to transform into the same brightness range. We can certainly choose one image as a template and transform all others into its brightness range, and often that is just fine. But what if the one image is an outlier and corrupted by noise or artifacts? Least squares isn’t famous for dealing well with outliers anyway, but in this case it’s more obvious that the choice of template image can introduce an unfavorable bias.

Those considerations quickly lead us into the territory of joint optimizations and latent images, where the math and algorithms get really interesting really fast, despite our simple model. I might tackle this in a follow-up article.

Endnotes