Keywords

These keywords were added by machine and not by the authors. This process is experimental and the keywords may be updated as the learning algorithm improves.

1 Introduction

The most popular paper about the original Retinex formulation is [19]. Retinex stands for ‘Retina plus Cortex’, which refers to the fact that the mechanisms underlying human color vision depend both on the retinal photoreceptors catches and on the cortex interpretation of this signals. The original Retinex is a computational model with the aim of finding a perceptual correlate of reflectance, called ‘lightness’ by Land, to be tested with psychophysical measurements.

Through a series of groundbreaking experiments, mostly performed with the famous ‘Mondrian tableaux’, Land and McCann proved that human perception of a surface’s color is much more influenced by the spatial distribution of the surrounding surfaces than by the spectral distribution of the light used to illuminate the Mondrian tableau. As underlined by McCann in many papers and conference speeches, spatial locality of color perception is the central concept in the whole Retinex theory. Thus, at least in its original form, the aim of Retinex is not to discard illumination and recover the intrinsic reflectance of surfaces, as several authors claim in their paper even nowadays, but to quantify how the points of the spatial surround cooperate to modify color perception.

In spite of their innovative and important experimental achievements, neither Land nor McCann ‘carved their model into stone’ through a rigorous mathematical formulation. In this paper, we are going to discuss the two major classes of Retinex that can be found in the literature: ratio-reset Retinex and Horn’s Retinex [11]. We will underline how profound is the difference between these two interpretations thanks to variational principles and partial differential equations.

Before entering in the mathematical details of variational formulations, it is worth introducing, in the following section, the basic formalization of the original Retinex formula developed in [25]. This will help us fix the ideas and the notation about many concepts and notations that will be discussed in the following sections.

2 Land and McCann’s Original Retinex Model

As previously commented, the original Retinex model of Land and McCann [19] is based on the assumption that the HVS operates with three retinal-cortical systems, each processing independently the low, middle and high wavelengths of the visible electromagnetic spectrum. Every independent process forms a separate image determining a quantity that they called lightness and denoted with L. Land and McCann found a computational way to reproduce lightness for their Mondrian tableaux by introducing spatial comparisons among intensities, calculated over paths. The comparison is performed through a multiplicative chain of ratios, subjected to these non-linear operations: Threshold mechanism: if the ratio does not differ from 1 more than a fixed threshold value, then it is set to be unitary; Reset mechanism: if the cumulated product of ratios overcomes the value 1 in a certain point of the path, then it is forced to 1, so that the computation restarts from it. In this way, this point becomes a local white reference, so hat the reset mechanism is responsible for the white-patch behavior of Retinex.

Let us now present the mathematical formalization of Land and McCann ratio-threshold-reset Retinex computation provided in [25]. Given a discrete digital image function with normalized range, \(I:\varOmega \subset \mathbb {Z} ^2 \rightarrow [0, 1]\), consider a collection of N oriented paths \({\varvec{\gamma }}=\{\gamma _1,\ldots ,\gamma _N\}\) composed by ordered chains of pixels starting in \(y_k\) and ending in x, \(k=1,\ldots ,N\). Let \(n_k\) be the number of pixels traveled by the path \(\gamma _k\) and let \(t_k=1,\ldots ,n_k\) be its parameter, i.e. \( \gamma _k:\{1,\ldots ,n_k\}\rightarrow \varOmega \subset \mathbb {R} ^2 \), \(\gamma _k(1)=y_k\) and \(\gamma _k(n_k)=x\). Write, for simplicity, two subsequent pixels of the path as \(\gamma _k(t_k)=y_{t_k}\) and \(\gamma _k(t_k + 1)=y_{t_k + 1}\), for \(t_k=1,\ldots ,n_k - 1\). Consider, in every fixed chromatic channel \(c\in \{R,G,B\}\), their intensities \(I(y_{t_k})\), \(I(y_{t_k + 1})\) and then compute the ratio \(R_{t_k}=\frac{I(y_{t_k +1})}{I(y_{t_k})}\) with the initial condition \(R_0=1\).

With this notation in mind, the value of lightness provided by the ratio-threshold-reset Retinex algorithm for a generic pixel \(x\in \varOmega \), in every fixed chromatic channel c (that we avoid specifying for the same of a clearer notation), is given by:

$$\begin{aligned} L_{\varepsilon ,{\varvec{\gamma }}}(x)= \frac{1}{N} \sum _{k=1}^N \prod _{t_k = 1}^{n_k - 1} \delta _k(R_{t_k}) \end{aligned}$$
(1)

where \(\delta _k: \mathbb {R} ^+ \rightarrow \mathbb {R} ^+\), \(k=1,\ldots ,N\), are functions defined in this way: \(\delta _k(R_0)=1\) and, for \(t_k=1,\ldots ,n_k-1\),

$$\begin{aligned} \delta _k(R_{t_k})= \left\{ \begin{array}{lll} R_{t_k} &{} \text { if }\; 0< R_{t_k} \le 1 - \varepsilon \\ 1 &{} \mathrm{if }\; 1-\varepsilon< R_{t_k} < 1 + \varepsilon \\ R_{t_k} &{} \text { if } 1+\varepsilon \le R_{t_k} \le \frac{1+\varepsilon }{\prod _{m_k = 0}^{t_k - 1} \delta _k(R_{m_k})} \\ \frac{1}{\prod _{m_k =0}^{t_k - 1} \delta _k(R_{m_k})} &{} \text { if }\; R_{t_k} > \frac{1 + \varepsilon }{\prod _{m_k = 0}^{t_k - 1} \delta _k(R_{m_k})} \end{array} \right. \end{aligned}$$
(2)

being \(\varepsilon >0\) a fixed threshold. The second option is the mathematical implementation of the threshold mechanism while the fourth implements the reset mechanism (and so the white patch behavior) of the algorithm.

It is useful to write the contribution of the single path \(\gamma _k\) to \(L_{\varepsilon ,{\varvec{\gamma }}}(x)\) as:

$$\begin{aligned} L_{\varepsilon ,\gamma _k}(x)=\prod _{t_k = 1}^{n_k - 1} \delta _k(R_{t_k}), \end{aligned}$$
(3)

so that formula (1) reduces simply to the average of these contributions, i.e. \(L_{\varepsilon ,{\varvec{\gamma }}}(x)= \frac{1}{N} \sum \limits _{k=1}^N L_{\varepsilon ,\gamma _k}(x)\).

2.1 The Limit Behavior \(\varepsilon \rightarrow 0\)

The analytical formula to describe the ratio-threshold-reset Retinex algorithm just introduced allowed making predictions about the model. As explained in [25], this can be done if the threshold mechanism is disregarded, or, equivalently, by considering the case \(\varepsilon \rightarrow 0\).

As \(\varepsilon \rightarrow 0\), the functions \(\delta _k\) become much simpler:

$$\begin{aligned} \delta _k(R_{t_k})= \left\{ \begin{array}{lll} R_{t_k} &{} \text { if }\; 0 < R_{t_k} \prod \limits _{m_k = 0}^{t_k - 1} \delta _k (R_{m_k}) \le 1 \\ \frac{1}{\prod \limits _{m_k = 0}^{t_k - 1} \delta _k(R_{m_k})} &{} \text { if }\; R_{t_k} \prod \limits _{m_k = 0}^{t_k - 1} \delta _k(R_{m_k}) > 1 \end{array} \right. \end{aligned}$$
(4)

hence, when \(\varepsilon \rightarrow 0\), \(\delta _k\) behaves either as the identity or the reset function.

In [25] it was proven that this implies the following formula:

$$\begin{aligned} L_{0,{\varvec{\gamma }}}(x)=\frac{1}{N}\sum _{k=1}^N \frac{I(x)}{I(y_{H_k})}, \end{aligned}$$
(5)

where \(y_{H_k}\) is the pixel with highest intensity traveled by \(\gamma _k\). From now on, we will refer to formula (5) as describing the ‘ratio-reset Retinex algorithm’.

Notice that he presence of paths makes the ratio-reset Retinex a local algorithm, where locality is intrinsically represented by the geometry of paths used. However, when \(n_k\rightarrow |\varOmega |\) or \(N\rightarrow \infty \), the ratio-reset Retinex loses its local properties and reduces, see [25], to the global diagonal von Kries model [16]. On the other hand, if we use small values of \(n_k\) or N, the resulting lightness images are affected by a lot of noise, see again [25].

Finally, it is important to underline that, since intensity values are normalized, \(0<I(y_{H_k})\le 1\) for every \(k=1,\ldots ,N\) and then \(\sum _{k=1}^N \frac{1}{I(y_{H_k})}\ge N\). It follows that \(L(x)\ge I(x)\) for every pixel i and this proves that an image filtered with the ratio-reset Retinex is always brighter or equal to the original one. This shows an important limitation of this algorithm: an over-exposed picture can only be worsened by the application of the ratio-reset Retinex used as a color corrector.

2.2 From Paths to Pixel Sprays: RSR and Related Algorithms

The information obtain thanks to the mathematical formulation of Retinex has important consequences on the structure of \(\mathcal{P}_x(\varOmega )\), the set of paths embedded in the image domain \(\varOmega \) and ending in the point x. After formula (5), on this set it is natural to define the following equivalence relation: given \(\gamma , \eta \in \mathcal{P}_x(\varOmega )\),

$$\begin{aligned} \gamma \sim \eta \quad \Leftrightarrow \quad \max _{y\in \gamma ^*} \{I(y)\} = \max _{y\in \eta ^*} \{I(y)\} \end{aligned}$$
(6)

where \(\gamma ^*\) and \(\eta ^*\) are the co-domains of the paths, i.e. the collections of pixels traveled by \(\gamma \) and \(\eta \), respectively.

Paths belonging to different equivalence classes give different contributions to the lightness computation, while every path in a given equivalence class gives rise to the same value of \(L_{0,\gamma _k}(x)\). It follows immediately that \(\mathcal{P}_x(\varOmega )\) contains redundant paths and that the correct set of paths to consider is given by the quotient set \(\mathcal{P}_x(\varOmega )/{\sim }\), whose elements are the equivalence classes of paths with respect to the equivalence relation defined in (6).

In each equivalence class one can choose a single representative path to compute \(L_{0,\gamma _k}(x)\), in particular, the more efficient one is the two-points path whose co-domain is simply given by \(\{y_{H_k},x\}\). Thus, the ordering operations needed to generate the paths are totally unnecessary for the final lightness computation. Moreover, by a mathematical point of view, paths are topological manifolds of dimension 1 embedded in the image, which is a topological manifold of dimension 2, so paths do not really scan local neighborhoods of a pixel, rather particular directions in these neighborhoods. This directional extraction of information can lead to halos or artifacts in the filtered image, see e.g. [7].

These considerations led the authors of [26] to consider 2-dimensional objects such as areas instead of 1-dimensional paths to analyze image locality for an efficient color correction. Roughly speaking, their idea is to implement spatial locality by selecting a fraction of pixels from these areas with a density sample that changes according to a given function of their distance with respect to the target pixel x. Each function generates a different kind of pixel selection around x, leading to different kind of ‘sprays’, each of which shows different local filtering properties. The new implementation of the ratio-reset Retinex that follows this idea is called RSR for ‘Random Sprays Retinex’.

In RSR the role of a path \(\gamma _k\) traveling \(n_k\) pixels and ending in the target x is played by \(S_k(x)\), a spray with \(n_k\) pixels centered in x. Actually, once the number of points per spray is chosen, there is no need to vary it with k, hence, from now on, we will write n instead of \(n_k\) to denote the number of pixels per spray. The ratio-reset operation along a path is substituted by the search of the pixel with highest intensity in the whole spray. The functional expression of formula (5) to compute the lightness remains exactly the same in both algorithms, so the ratio-reset Retinex and RSR share the same intrinsic properties.

In [1, 2] two techniques have been proposed to reduce noise generation also decreasing the computational time of RSR.

In [27] the spray technique was used to fuse RSR with ACE [28], another perceptually-inspired color correction algorithm that makes use of the gray-world hypothesis [5]. The hybrid algorithm is called RACE and it is able to color correct both under and over exposed images.

A more recent proposal to fuse WP and GW features in a single algorithm is that presented in [15] and called STRESS for Spatio-Temporal Retinex-like Envelope with Stochastic Sampling. As Retinex, STRESS computes, for each pixel, the local white reference, but also the black reference in each chromatic channel. This is done through calculating the maximum and minimum envelope functions, denoted as \(E_{\max } (x)\) and \(E_{\min } (x)\), respectively.

Finally, let us mention that, in [9], the RSR sampling technique has been studied from a probabilistic point of view, resulting in the algorithm QBRIX and, in [10], further comparisons among Retinex models are discussed.

3 A Variational Framework for the Ratio-Reset Retinex

The similarities between the ACE formula [28] and the gradient descent equations for histogram equalization obtained in [29], led to the discovery of a variational interpretation of ACE in the paper [4]. The framework were further extended in [23] and, finally, in [3] a variational framework for (an anti-symmetric version of the) ratio-reset Retinex has been discussed. In order to understand how this is possible, let us come back to the lightness formula (5).

Land and McCann proposed a further Retinex mechanism, the scaling, implemented via a strictly increasing function \(f:(0,1]\rightarrow (0,1]\) such that \(f(r)\ge r\) for all \(r\in (0,1]\) applied to the ratio \(r=\frac{I(x)}{I(y_{H_k})}\), so that the Retinex lightness formula becomes:

$$\begin{aligned} L_{0,{\varvec{\gamma }},f}(x)=\frac{1}{N}\sum _{k=1}^N f\left( \frac{I(x)}{I(y_{H_k})}\right) . \end{aligned}$$
(7)

The reset mechanism of Retinex and the scaling operation can be merged: in fact, we can extend f to \((0,+\infty )\) preserving its continuity by defining

$$\begin{aligned} \hat{f} (r) = {\left\{ \begin{array}{ll} f(r) &{} \text {if }r\in (0,1] \\ 1 &{} \text {if } r\in [1,+\infty ). \end{array}\right. } \end{aligned}$$

It is clear that applying this new scaling function \(\hat{f}\) to the ratios I(x)/I(y), with x fixed and y that varies in \(\varOmega \), jointly implements the scaling and the reset mechanism.

Now we have all the elements to introduce the continuous version of the Retinex algorithm presented in [3] under the name ‘Kernel-Based Retinex’, or KBR for short. Given \(x\in \varOmega \), let \(Y_{w,x}\) be the random variable modeling the selection of a pixel in the neighborhood of x according to the density w(x, y).

The output \(L^\text {KBR}_w(x)\) of the KBR algorithm at the pixel x is defined as the conditional expectation of \(\hat{f}\left( \frac{I(x)}{I(Y_{w,x})}\right) \) with respect to the distribution w of pixels around x, i.e.

$$\begin{aligned} L^\text {KBR}_w(x) = \mathbb {E}_{Y_{w,x}}\left[ \hat{f} \left( \frac{I(x)}{I(Y_{w,x})} \right) \right] . \end{aligned}$$
(8)

This formula is used independently for each color channel and can be written more explicitly as

$$\begin{aligned} \begin{array}{ll} \displaystyle L^\text {KBR}_w(x)=\sum _{\{y\in \varOmega : I(y) \ge I(x)\}} w(x,y) \, f\left( \frac{I(x)}{I(y)}\right) + \sum _{\{y\in \varOmega : I(y) < I(x)\}} w(x,y). \end{array} \end{aligned}$$
(9)

All the basic properties of the ratio-reset Retinex are faithfully implemented in (9): KBR is founded on the propagation of a two-pixel ratio comparison between the fixed target x and the generic pixel y that runs across the image; these comparisons are then subjected to the reset and scaling performed by \(\hat{f}\) and, finally, locally averaged with weight w, in order to produce the value of \(L^\text {KBR}_w(x)\).

To study the action of KBR of pixel intensities, it is useful to rewrite (9) introducing the functions

$$\mathrm{sign}^+(\xi ):=\left\{ \begin{array}{ll} 1 &{} \text {if } \xi > 0,\\ \frac{1}{2} &{} \text {if } \xi = 0,\\ 0 &{} \text {if } \xi < 0, \end{array} \right. \qquad \mathrm{sign}^-(\xi ) = 1 - \mathrm{sign}^+(\xi ),$$

so that Eq. (9) can be re-written as

$$\begin{aligned} \displaystyle L^\text {KBR}_w(x) = \sum _{y \in \varOmega } w(x,y) \, f\left( \frac{I(x)}{I(y)}\right) \, \mathrm{sign}^+(I(y)-I(x)) + \sum _{y\in \varOmega } w(x,y) \, \mathrm{sign}^-(I(y)-I(x)). \end{aligned}$$
(10)

Thanks to Eq. (10) we can verify that KBR always increases brightness as the original Retinex implementation. In fact, since \(f(r) \ge r\) for all \(r\in (0,1]\), then \(f\left( \frac{I(x)}{I(y)}\right) \ge \frac{I(x)}{I(y)} \ge I(x)\), so

$$\begin{aligned} L^\text {KBR}_w(x) \ge \sum _{y \in \varOmega } w(x,y) \, I(x) \, \mathrm{sign}^+(I(y)-I(x)) + \sum _{y\in \varOmega } w(x,y) \, \mathrm{sign}^-(I(y)-I(x)) \end{aligned}$$
(11)

moreover, being \(I(x)\le 1\), we can write

$$\begin{aligned} L^\text {KBR}_w(x)&\ge \sum _{y \in \varOmega } w(x,y) \, I(x) \, \mathrm{sign}^+(I(y)-I(x)) + \sum _{y\in \varOmega } w(x,y) \, I(x) \, \mathrm{sign}^-(I(y)-I(x)) \nonumber \\&= I(x) \sum _{y \in \varOmega } w(x,y) \, \left[ \mathrm{sign}^+(I(y)-I(x)) + \mathrm{sign}^-(I(y)-I(x)) \right] \nonumber \\&= I(x) \sum _{y \in \varOmega } w(x,y) = I(x), \end{aligned}$$
(12)

having used the fact that the kernel is normalized. As in the original formulation, this property implies that over-exposed pictures could not be enhanced with Retinex unless we use a post-processing stage and that further iterations of Retinex keep on increasing the intensity until a white image is reached.

This equation of KBR does not correspond to the minimization of an energy functional. However, let us consider the sum of the function \(f\left( \frac{I(x)}{I(y)}\right) \, \mathrm{sign}^+(I(y)-I(x))\) and of its the anti-symmetrized version on the region \(\{x\in \varOmega : I(y) \le I(x)\}\), i.e.

$$\begin{aligned} {\begin{matrix} L^\text {aKBR}_w(x) &{} = \sum _{y \in \varOmega } w(x,y) \,f\left( \frac{I(x)}{I(y)}\right) \, \mathrm{sign}^+(I(y)-I(x)) \\ &{} - \sum _{y\in \varOmega } w(x,y) \, f\left( \frac{I(y)}{I(x)}\right) \, \mathrm{sign}^-(I(y)-I(x)) \end{matrix}} \end{aligned}$$
(13)

where aKBR stands for anti-symmetrized KBR.

In [3] it was proven that the right-hand side of the previous equation can be interpreted as the minimization of the energy functional given by:

$$\begin{aligned} C^f_w(I)=\sum _{x\in \varOmega } \sum _{y\in \varOmega } w(x,y) f\left( \frac{\min (I(x),I(y))}{\max (I(x),I(y))}\right) \!. \end{aligned}$$
(14)

Minimizing \(C^f_w(I)\) corresponds to maximizing the contrast in a local (due to the presence of the weight w) and non linear way (due to the ratio and to the presence of f). This explained in a quantitative and qualitative way how and why the somewhat involved ratio-reset mechanism of Retinex allows for a unilateral contrast enhancement, always directed towards the highest intensity.

KBR, ACE, RACE and STRESS corrected this unilateral behavior. In [24] the spatially-based variational framework was translated into a wavelet-based setting.

4 Retinex: A ‘melody’ that Everyone Plays Differently

In image processing it is hard to find a model whose name has been interpreted in so many different ways as ‘Retinex’. In this section, we present a synthetic description of the evolution of the Retinex interpretation.

Path-wise Retinex share a local WP nature and mostly differ from each other by the path geometry used to explore spatial locality: Land and McCann used piecewise linear paths in [19]. In [6], [21], and [30] those paths were substituted by double spirals, Brownian paths and traces of a specialized swarm of termites, respectively.

Center/surround Retinex are local GW algorithms originated from [18], where Land noticed that he could reproduce Mach bands originated by a spinning white square on a black background by using a different Retinex formulation. Precisely, for every image point, the intensity of the center x is replaced by the ratio between I(x) and the average value of the surround, sampled with a density that decays as the inverse of the square distance from the center. Writing with \(L^\mathrm{CS}\) this ‘center/surround lightness’, we have: \(L^\mathrm{CS}(x) = I(x)/{<}\{I(y), \, y\in \mathrm{Surround}\}{>}\), where \(< \cdot>\) represents the average operator. Comparing this last formula with (5), it can be seen that there is a fundamental difference between this formulation and the original one: there the ratio is performed over the pixel with highest intensity, while in this formulation it is implemented over the mean value of the surround. In practice, this last formulation can be seen as a gray-world method to remove the illuminant component of the image [5].

In 1997, Johbson, Rahman and Woodell [13] re-elaborated Land’s idea presented in [18]: they worked with logarithmic data, approximating the average of the surround by convolving the image function I with a normalized kernel function F, usually a Gaussian. If we use again, for simplicity, the symbol \(L^\mathrm{CS}\), we can write this model as follows: \(L^\mathrm{CS}(x)= \log (I(x)) - \log ((F*I)(x))\), \(\forall x\in \varOmega \).

Multilevel Retinex algorithms were pioneered by Frankle and McCann in [7] and further refined in [8]. In these works a multilevel version of the original local WP Retinex is presented, the authors abandon paths and consider a computation that takes into account all pixels. The input image is progressively sub-sampled averaging a number of pixel that grows as increasing powers of 2. On each sub-sample level a ratio-reset computation (without threshold) is operated a certain number of times, from the coarser sub-sample level to the finest one. Because of the sub-sampling, as we go far from the target pixel, we do not consider actual pixel values, but average values of macroareas of increasing size. A rigorous mathematical formulation of these multilevel algorithms is still lacking.

Based on this idea, Marini, Rizzi and De Carli [21] constructed a local WP multilevel version of Brownian path Retinex that reduced the amount of noise in the output images. A different multilevel proposal has been pointed out by Johbson, Rahman and Woodell in [12]: they introduced a certain number S of scales where performing the convolutions with normalized Gaussian functions \(F_s\), \(s=1,\ldots ,S\). Each scale is associated to a suitable weight \(w_s\), which gives more importance to finer scales than to coarser ones.

Finally there are WP Retinex versions based on solving a Poisson equation. They rely on a work of Horn [11], in which he remarkably pointed out, for the first time, the need for a spatially isotropic two-dimensional version of Retinex. Horn considered, as Land, only Mondrian tableaux illuminated by a smoothly varying light. However, differently from Land, he explicitly tackled the ill-posed problem of inverting the equation \(I_c(x) = S_c(x) L_c(x)\), \(c\in \{R,G,B\}\), with respect to \(S_c(x)\), the reflectance of the point x, knowing only the image intensity \(I_c(x)\) and not the illumination \(L_c(x)\). If we pass to logarithmic values, i.e. \(\log I_c(x) = \log S_c(x) + \log L_c\) or, equivalently, \(\log S_c(x) = \log I_c(x) - \log L_c\) and we apply a differential operator D to both sides, then \(D(\log L_c(x))\) will be small but finite everywhere, while \(D(\log S_c(x))\) will be different from zero only if x is close to sharp edges.

If we apply a threshold operator \(\delta _T\) defined as follows:

$$\begin{aligned} \delta _T(s) = {\left\{ \begin{array}{ll} s &{} \text {if } |s|>T \\ 0 &{} \text {elsewhere}, \end{array}\right. } \end{aligned}$$

for all \(s\in \mathbb {R} \) and if the threshold \(T>0\) is small enough, then we obtain \(D(\log S_c(x))=\delta _T(D(\log I_c(x)))\). Horn insisted on the choice of the Laplacian for D instead of the gradient, arguing that first order derivatives are one-dimensional, while the second order derivatives involved in the Laplacian are isotropic and thus more suited for the topology of an image. By substituting D with the Laplacian operator \(\varDelta \), the last formula becomes a Poisson equation:

$$\begin{aligned} \varDelta (\log S_c(x))=\delta _T(\varDelta (\log I_c(x))), \end{aligned}$$
(15)

whose solution allows to recover the logarithmic reflectance \(\log S_c(x)\). It is clear that Horn’s method is based on quite restrictive hypotheses: smoothness of illumination (violated by scenes with deep shadows, for instance) and a Mondrian-like world (violated each time edges are not sharp).

5 Mathematical Formalizations of Horn’s Interpretation

Besides the variational framework described in Sect. 3, in the literature there exist alternative variational models of Retinex-like algorithms and also formalizations based on partial differential equations (PDE). The aim of this section is not to give an exhaustive list, rather to discuss the main features of the most famous alternative mathematical formalizations of Retinex-like algorithms present in the literature.

The first authors to embed a Retinex-like algorithm in a variational framework were Kimmel and colleagues in [14]. They did not considered the original Land’s ratio-threshold-reset Retinex, but Horn’s interpretation. In fact, they started from the logarithmic equation \(\log I_c(x) = \log S_c(x) + \log L_c(x)\), \(c\in \{R,G,B\}\) and tried to solve it with respect to \(\log L_c(x)\) by imposing the hypothesis of smoothness on the illuminant part of the logarithmic image. Once obtained an estimation of the illumination, they could infer the reflectance information \(S_c(x)\). This one then undergoes suitable transformations and gives an illuminant-invariant version of the original image.

It is important to underline a fundamental difference between this variational technique and the one presented in the previous sections: here contrast enhancement of the original image \(\log I_c(x)\) is obtained by decreasing the contrast of the illuminant image \(\log L_c(x)\). In fact, \(\log I_c(x)\) is measured by the camera and so it is a fixed data, \(\log L_c(x)\) is estimated by using a smoothness prior, thus the estimated reflectance \(\log S_c(x) = \log I_c(x) - \log L_c(x)\), or \(S_c(x)=I_c(x)/L_c(x)\) is forced to have a stronger contrast than the original image data. Instead, the variational principles previously discussed act directly on the contrast of the original image, without taking into account the separation between reflectance and illuminant and related approximations and priors.

Avoiding the subscript c, the functional proposed in [14], with the notations of this paper, can be expressed as follows:

$$\begin{aligned} \small {E_{\alpha ,\beta }(\log L)=\sum _{x\in \varOmega } [|\nabla \log L(x)|^2+\alpha (\log L(x)-\log I(x))^2+\beta |\nabla (\log L(x)-\log I(x))|^2]} \end{aligned}$$
(16)

with the constraints \(\log L(x) \ge \log I(x)\), because the reflectance S(x) is always between 0 and 1, and the boundary condition \(\langle \nabla \log L, {\varvec{n}} \rangle =0\) on \(\partial \varOmega \), i.e. \(\log L\) orthogonal to the normal \({\varvec{n}}\) to the boundary \(\partial \varOmega \) of \(\varOmega \).

The first term of the functional forces spatial smoothness on the illumination L. The authors chose that particular analytical form because the Euler-Lagrange equation associated to \(\sum _{x\in \varOmega } |\nabla \log L(x)|^2\) is the Laplace PDE \(\varDelta \log L=0\), whose steepest descent solution is equivalent to a Gaussian smoothing. The second penalty term forces a proximity between \(\log L\) and \(\log I\), so that their difference \(\log S\), the logarithmic reflectance, tends to 0, i.e. the reflectance R tends to 1, or white. The authors declare that the principal objective of this term is to regularize the problem, so that it is better conditioned in view of a numerical solution and they set the constant \(\alpha \) to a very small value not to force too much \(\log L\) towards \(\log I\). The third term represents a Bayesian penalty, which forces reflectance gradients to be smooth. The authors declared to have introduced it to force R to be visually pleasing, without abrupt variations.

Morel, Petro and Sbert [22] analyzed Land’s original Retinex model [17] without the reset mechanism. They showed that, if the Retinex paths are interpreted as symmetric random walks, then Retinex is equivalent to the following Neumann problem for a Poisson equation:

$$\begin{aligned} {\left\{ \begin{array}{ll} -\varDelta L(x)=F(x) &{} x\in \varOmega \\ \frac{\partial L(x)}{\partial {\varvec{n}}}=0 &{} x\in \partial \varOmega , \end{array}\right. } \end{aligned}$$

where F is a suitable scalar field, see [22] page 2830.

Let us now consider the algorithm STRESS. We recall that the basic information needed by STRESS is given by the two envelope functions \(E_{\min }\) and \(E_{\max }\) which, in the original formulation, are computed through the same random spray technique of RSR [26]. To avoid the typical noise problems related to this technique, in [31], the authors proposed to compute the envelope functions via the minimization of a functional based on the total variation, instead of using the random spray technique. For this reason the corresponding algorithm is called STRETV and corresponds to the minimization of the following functional for E (in this case E denotes the envelope and not the energy functional):

$$\begin{aligned} \sum _{x\in \varOmega } \left[ |\nabla E(x)| + \frac{\lambda }{2} |E(x)-I(x)|^2\right] \end{aligned}$$
(17)

subjected to \(E(x)\ge I(x)\) to compute \(E_{\max }\) and to \(E(x)\le I(x)\) for \(E_{\min }\).

The minimization of the first (total variation) term, assures the spatial smoothness of the envelope functions, the second term is a fidelity term used not to depart too much from the original image values. The authors declare that the coefficient \(\lambda \) must be \({\ll }1\) for good results. The authors do not specify if they consider a spatial kernel to localize their computation or not.

The last variational formalization that we discuss here is that presented in [20] relative to the termite Retinex. Here an energy functional is taken into account to determine the geometry of the paths used by Retinex. Fixed a pixel \(x\in \varOmega \), the authors search for the path \(\gamma :[0, 1]\rightarrow \varOmega \), \(\gamma (0)=x\), the minimizes the energy functional defined by:

$$\begin{aligned} E(\gamma )=\int _0^1 \left[ \frac{1}{1+(D^2-\Vert x-\gamma (s)\Vert ^2)\Vert \nabla I(\gamma (s))\Vert ^2}+\theta (\gamma (s))\right] ds, \end{aligned}$$
(18)

where D is the diagonal of \(\varOmega \) and 1 is introduced to avoid singularities in the case the denominator is 0. The paths that minimize \(E(\gamma )\) are those which balance the fact to remain as close to x as possible and, simultaneously, to explore image areas with high values of the gradient. Both features maximize the denominator of the first term. If x lies in an area with a high density of edges, \(\gamma \) will not go too far from x, instead, if x lies in a rather homogeneous area, \(\gamma \) will be forced to explore the image points far away from x to find the important edge information. \(\theta (\gamma (s))\), the so-called ‘poison term’, is set to zero at the beginning, and it increases each time a path has been traveled, to prevent from exploring the same image area all the time. Once a set of N path has been selected, the intensity I(x) of the pixel x in each separate chromatic channel is modified with the Retinex formula 5.

6 Conclusions

In the past fifteen years, variational methods have been used to formalize color correction algorithms. This permitted to point out similarities and differences among several models that were difficult to detect just looking at their direct equations. In this paper we have described, in a self-contained way, both the direct and the variational versions of several color enhancement algorithms inspired by the seminal Retinex theory of color vision. A particular emphasis has been put in highlighting the very different variational formulations of the original Retinex of Land and McCann and those referring to Horn’s interpretation, which are often misleadingly mixed in the literature.