Retinex Poisson Equation : a Model for Color Perception

This workshop may eventually be submitted and published.

Overview
References
Online Demo
The PDE-Retinex Model
The Algorithm
Implementation
Examples

Nicolas Limare nicolas.limare@cmla.ens-cachan.fr, CMLA, ENS Cachan, CNRS, UniverSud
Jean-Michel Morel morel@cmla.ens-cachan.fr, CMLA, ENS Cachan, CNRS, UniverSud
Ana Belén Petro anabelen.petro@uib.es, TAMI, Universitat de les Illes Balears
Catalina Sbert catalina.sbert@uib.es, TAMI, Universitat de les Illes Balears

Overview

In 1964 Edwin H. Land (ref. 4) formulated the Retinex theory, the first attempt to simulate and explain how the human visual system perceives color. His theory and an extension, the ``reset Retinex'' (ref. 5) were further formalized by Land and McCann in 1971. Several Retinex algorithms have been developed ever since. These color constancy algorithms modify the RGB values at each pixel to give an estimate of the physical color independent of the shading.

The Retinex original algorithm was both complex and not fully specified. Indeed, this algorithm computes at each pixel an average of a very large and unspecified set of paths on the image. For this reason, Retinex has received several interpretations and implementations which, among other aims, attempt to tune down its excessive complexity.

But, as shown in ref. 1, the original Retinex algorithm can be formalized as a (discrete) partial differential equation. More precisely, it can be shown that if the Retinex paths are interpreted as symmetric random walks, then Retinex is equivalent to a Neumann problem for a Poisson equation. This result gives a fast algorithm involving just one parameter, also present in the original theory.

The Retinex Poisson equation (given below) is very similar to Horn’s and Blake’s equations, which were proposed as alternatives to Retinex. It also is one of the “Poisson editing” equations proposed in Perez et al. (ref. 3). The final principle of the algorithm is extremely simple. Given a color image I, its small gradients (those with magnitude lower than a threshold t) in each channel are replaced by zero. The resulting vector field is no more the gradient of a function, but the Poisson equation reconstructs an image whose gradient is closest to this vector field. Thus, a new image is obtained, where small details and shades of the original have been eliminated. The elimination of the shades creates more homogeneous colors. This fact, according to Land and McCann, models the property of our perception to perceive constant colors regardless of their shading.

The formalization proved in ref. 1 yields a fast implementation of the Land-McCann theory using only two DFT's. You can the theory on line on your own color images.

References

Jean-Michel Morel, Ana Belén Petro and Catalina Sbert, A PDE Formalization of the Retinex Theory. IEEE Trans. on Image Processing, in print, 2010.
Jean-Michel Morel, Ana Belén Petro and Catalina Sbert, Fast Implementation of color constancy algorithms. Color Imaging XIV: Displaying, Processing, Hardcopy and Application. Proc. of Electronic Imaging SPIE, vol 7241. January 2009. (preprint)
P. Pérez, M. Gangnet and A. Blake, Poisson Image Editing. ACM Transactions on Image Processing. Proc. of ACM SIGGRAPH 2003, vol 22, Issue 3 Pages: 313 - 318 (July 2003).
Edwin H. Land, "The retinex". Am. Sci. 52(2): 247-64. (1964).
Edwin H. Land and John J. McCann, Lightness and Retinex Theory. J. Opt. Soc. Am. 61, 1-11 (1971).

Online Demo

An online demo allows you to try Retinex with your own images. The demo has only one parameter, the contrast threshold present in the original theory. The images can be slightly resized for an efficient Fourier transform.

As the input image has to be normalized after the Retinex algorithm, the demo also provides the result of this normalization without Retinex. We use the "Simplest Color Balance" algorithm for this normalization.

The PDE-Retinex Model

In ref. 1 it is proven that the output of the Retinex algorithm proposed by Land and McCann is the solution of the discrete partial differential equation with Neumann boundary conditions

$- \Delta_d u(i,j) = F(i,j)$

where

$\Delta_d u(i,j) = u(i+1,j) + u(i-1,j) + u(i,j+1) + u(i,j-1) - 4 u(i,j)$

is the discrete Laplacian, dim = N M is the size of the image,

F(i,j) = f(I(i,j) - I(i+1,j)) + f(I(i,j) - I(i-1,j))
+ f(I(i,j) - I(i,j+1)) + f(I(i,j) - I(i,j-1))

and f(x) is a threshold function, whose value is zero if |x| < t and the identity in other case and I is the image to process. This function f eliminates the small variations of the intensity image I.

The parameter t (the threshold) is by default t = 4 but you can choose the value depending of the variations you want to eliminate. When I is a gray level image, the algorithm applies to I. When I is a color image, the algorithm is applied to each scalar channel separately.

The Algorithm

The output of the algorithm are two images: the first one is the white balance of the original color image (on each channel the darkest pixels are put to zero and the brightest ones are put to 255); the second image is the result of the Retinex PDE applied to the three channels of the color image, completed by the same white balance.

The discrete partial differential equation is easily solved by Fourier transform. To enforce the Neumann boundary condition, the image is first mirrored across its right, and bottom sides to obtain an four times larger image which is symmetric with respect to its vertical and horizontal medial axes.

The discrete Fourier transform of a two-dimensional function u(n,m) defined on a N x M grid is defined for (k,l) in {0,..., M-1}x{0,..., N-1} by

$\widehat{u}(k,l) = \frac{1}{NM}\sum_{n=0}^{N-1}\sum_{m=0}^{M-1} u(n,m)e^{-i\frac{2\pi k n}{N}} e^{-i\frac{2\pi l m}{M}}$

and the discrete inverse Fourier transform for (m,n) in {0,... , M-1}x{0,..., N-1} by

$u(n,m)=\sum_{k=0}^{N-1}\sum_{l=0}^{M-1} \widehat{u}(k,l)e^{i\frac{2\pi k n}{N}} e^{i\frac{2\pi l m}{M}}$ .

The discrete Fourier transform has the following property

$u(n-n_0,m-m_0)=\sum_{k=0}^{N-1}\sum_{l=0}^{M-1} \widehat{g}(k,l)e^{i\frac{2\pi k n}{N}} e^{i\frac{2\pi l m}{M}}\;$ ,

where

$\widehat{g}(k,l)=\widehat{u}(k,l) e^{-i\frac{2\pi k n_0}{N}} e^{-i\frac{2\pi l m_0}{M}} \;$ .

Applying the discrete Fourier transform to the discrete Poisson equation and using this last property yields

$\widehat{u}(k,l)\left(4 -2\cos\frac{2\pi k}{N}-2\cos\frac{2\pi l}{M}\right) = \widehat{F}(k,l) \;$ ,

which entails

$\widehat{u}(k,l)= \frac{\widehat{F}(k,l)}{4 -2\cos\frac{2\pi k}{N}-2\cos\frac{2\pi l}{M}}, \; \hbox{ for } (k,l)\neq (0,0)$ .

Using the inverse Fourier transform we obtain the value of u at each point of the grid, defined up to a constant since the constant Fourier coefficient is arbitrary. The values of u are finally normalized to the interval [0,255], which amounts to apply a white balance. All of the above computations are performed on the extended symmetric image F defined on the 2N x 2M grid. F being symmetric, its Fourier coefficients are real. This property is transferred by the equation to the Fourier coefficients of u, and u is therefore also symmetric and verifies the Neumann boundary condition. All of these operations are performed for each channel of the color image, u being in turn the red, green and blue channel.

The algorithm (applied to each channel) therefore is

Compute F(i,j);
Compute the Fourier transform of F by DFT;
Deduce the Fourier transform of u using the formula above;
Compute the final solution u by the inverse DFT and apply a white balance.

Implementation

The retinex_pde implementation is available : source code, documentation (online).

It should compile on any system since it's only ANSI C. This implementation is used in the online demo.

This code requires libpng for PNG file input/output and libfftw3 to process the Fourier transforms.

Linux. You can install libpng and libfftw3 with your package manager.
Mac OSX. You can get libpng and libfftw3 from the Fink project.
Windows. Precompiled DLLs are available online for libpng and libfftw3.

Compilation and usage instructions are provided in the README.txt file.

Implementation notes:

the white balance is applied to the two output images with 3% saturation, using the "Simplest Color Balance" algorithm;
with the libfftw3 Fourier transform, the image mirroring is emulated by the cosine transform functions.

Examples

Here there are some examples, note that Retinex is not a model for image enhancement or image improvement; it is only an algorithm to mimic our color perception. Thus, Retinex will enhance an effect that our perception does anyway.

The role of the t threshold is to eliminate the small intensity variations due to shading. Thus t cannot be too large (less than 10 typically) to avoid removing significant details. But is must be large enough to eliminate light shading effects. This is not always possible, since shadows can be very contrasted.

Adelson's checker

A first classic example, shows Adelson's checker shadow illusion. In the left image a green cylinder standing on a black and white checker-board casts a diagonal shadow across the board.

The image has been so constructed that the white squares in the shadow, one of which is labeled "B," have actually the very same gray value as the black squares outside the shadow, one of which is labeled "A."

If Retinex is faithful to human perception, it should make B much brighter than A, and it does.

original image (A and B have a gray value of 120)	white balance image (A and B have gray value of 88)

Retinex result with t=3 (A has a gray value of 85 and B of 94)	Retinex result with t=8

Circles on a Gradient

Simultaneous contrast, namely the fact that the appearance of a color depends on the colors surrounding it. The original image (already white balanced) shows a background with a smooth variation and two circles with the same gray value (170). One of them is placed in the darker part of the image, and the other one in the brighter part.

The usual perception is that the circle in the darker part looks conspicuously brighter than the other. If we use a threshold t=3 larger than the background variation, the result is an image with early constant background (45). The left circle gets a 18 gray value and the right circle a 245 gray value, which predicts well with our perception. In this image, the normalization have to done with care.

original image, white balanced	Retinex result

Color Illusion

This example shows the successful simulation of a color illusion. On the left image, two X's with exactly the same color are put on different backgrounds, one yellow and the other purple.

The striking illusion is that the X on the left seems to have a color similar to the purple background color on the right, while the X on the right has an apparent color similar to the yellow background color on the left. Applying the Retinex algorithm the X's really get these illusory colors.

Noisy Image

To understand the effect of the threshold t the next image shows a noisy original and the result of Retinex with increasing threshold values t=1, t=3 and t=5.

The background clutter and the shades are progressively filtered out when t increases, but the main edges are kept. At t=5, however, edges start loosing contrast and low contrasted details could disappear.

original image	White balance	Retinex result with t=1

Retinex result with t=3	Retinex result with t=5

Dark Image

An example with a dark image. There is no significant difference between the white balanced image and the Retinex result. A careful examination shows nonetheless the color on homogeneous regions becoming more constant.

original Image	white Balance image	Retinex result

Shadows removal

These two images demonstrate how Retinex can be used for removing shadows. It is not always effective, but here we can observe two good examples. The shadow removal works only if the boundary of the shadow is blurry, and therefore has a small gradient.

original image	White balance	Retinex result

original image	White balance	Retinex result

Lena

Application to Lena.

The smooth shading variations on the shoulder or the face and in the background fade out when the threshold t increases.

original Image	white balance image	Retinex result with t=8	Retinex result with t=12

credits

Edward H. Adelson
the authors CC-BY
standard test image