From image processing to computational neuroscience: a neural model based on histogram equalization

1. Introduction

The human visual system works in many ways in order to efficiently encode the visual information coming from natural environments, reducing its inherent redundancy, as proposed in the seminal work of Barlow (1961) (see Olshausen and Field, 2000 for a review). For instance, while natural scenes have luminance distributions which are very lopsided, with a high peak and a very rapid fall-off, photoreceptors encode this information with signals that have a much more even distribution: indeed, photoreceptors perform histogram equalization, as demonstrated by Laughlin (1981). And the receptive fields of visual neurons, both retinal and post-retinal, compensate the 1/f² decay of the power spectrum of natural images, whitening the spectrum of the resulting signal and thus minimizing interpixel redundancies and increasing coding efficiency (see Atick, 1992; Dan et al., 1996 where the existence of whitening at the local geniculate nucleus is demonstrated for natural images).

Apart from efficiency in coding, another very important aspect is that of biological efficiency in terms of wiring. The resolution of retinal mosaics is limited by the number of axons that can pass through the optic nerve, which acts as a bottleneck (Olshausen, 2003). But the visual system is able to achieve a visual acuity beyond the limit imposed by the number of photoreceptors at the retina: in their classical paper on contrast constancy, Georgeson and Sullivan (1975) suggest that there are cortical mechanisms for contrast enhancement that compensate for the blurring produced in early stages of the visual process. Very recently Martinez et al. (2014) have confirmed that contrast enhancement takes place at the lateral geniculate nucleus (LGN) and, remarkably, the authors point out that this contrast enhancement procedure is very much alike the common techniques used in image processing.

Alongside mechanisms for coding and wiring efficiency, there is neural activity in region V1 of the human visual cortex that correlates with the perceptual phenomenon of lightness induction, as proven by Pereverzeva and Murray (2008). The term lightness induction or achromatic induction designates the visual phenomenon by which the perceived reflectance of an object depends on its surround. It can take the form of lightness contrast, when the object's lightness shifts away from that of its surroundings: a dark object on a light background appears even darker, or a light object in a dark surround becomes even lighter. The reverse is called lightness assimilation, in which case the appearance of the object shifts in the direction of the lightness of its surround. As pointed out by Shevell (2003), lightness assimilation occurs in situations of high spatial frequency while lightness contrast is associated with relatively lower spatial frequencies.

Our contribution in this paper is to propose a neural activity model, a partial differential equation (PDE) in the form of a Wilson-Cowan equation (Wilson and Cowan, 1972), which takes care simultaneously of the four aspects mentioned above: it performs histogram equalization, spectrum whitening, contrast enhancement, and it also predicts lightness induction. The proposed model is based on a state of the art method for color and contrast enhancement from the image processing literature, so we start the following section reviewing some key image processing concepts.

2. Image Processing for Contrast Enhancement

2.1. Histogram Equalization

Histogram equalization is a classical, very basic image processing technique dating at least to the early 1970s (see Pratt, 2007 and references therein), aiming at enhancing the contrast and improving the appearance of images by way of re-distributing their levels uniformly accross the available range. In this sense, an image would be optimal if its histogram were flat or “equalized,” meaning that all the range is used and all levels are represented by the same amount of pixels. Therefore, when an image has a flat histogram its cumulative histogram is simply a ramp, and this allows for a very straightforward computation for the histogram equalization procedure: assuming we are working on a graylevel image in the range [0,1], we have to substitute each level g in the original image by the value of its normalized cumulative histogram, H(g). The solution is computed very fast using a look-up table (LUT). An example result can be seen in Figure 1 (notice that, while the range has been expanded and the resulting image has a more even histogram, it's not actually uniform).

FIGURE 1

Figure 1. Left: image and associated histogram. Right: after histogram equalization.

While in Figure 1 histogram equalization improves the visual appearance of the image, Figure 2 shows an example where the image is actually made worse, which Pratt (2007) points out is often the case when the image is overexposed, as it is here. This is aggravated by the fact that the equalization procedure is a one-shot technique, that only produces a final result, without any “in-between,” so if the resulting image shows any type of unpleasant artifact there is nothing to do about it. This issue was addressed by Sapiro and Caselles (1997), who proved that the minimization of the energy functional

$\begin{array}{cc}E(I)=2{\displaystyle \sum _x{\left(I(x)-\frac{1}{2}\right)}^2}-\frac{1}{AB}{\displaystyle \sum _x{\displaystyle \sum _y|}}I(x)-I(y)|& (1)\end{array}$

produces an image I with a flat histogram. The range of I is [0,1], x, y are pixels and A, B are the image dimensions. While the result of histogram equalization is very often unsatisfactory and can't be altered, Sapiro and Caselles (1997) propose to start with an input image I₀ and apply to it step after step of the minimization of Equation 1, letting the user decide when to stop. If the user lets the minimization run to convergence, she'll get the same result as with a LUT, but otherwise a better result can be obtained if the iterative procedure stops before the appearance of severe artifacts. The squared differences in the first term of Equation 1 and the absolute differences in the second one are required to ensure that the minimization yields an image with equalized histogram, see Sapiro and Caselles (1997) for details. The energy in Equation 1 can be interpreted as the difference between two positive and competing terms,

$\begin{array}{cc}E(I)=D(I)-C(I).& (2)\end{array}$

The first term measures the dispersion around the average value of ½, as in the gray world hypothesis for color constancy, stating that our visual system estimates the illuminant as one half the average of the colors of the scene, an observation made by Judd (1940, 1979a) and formalized by Buchsbaum (1980). The second term measures the contrast as the sum of the absolute value of the pixel differences.

FIGURE 2

Figure 2. Left: original image. Right: after histogram equalization.

2.2. Perceptually-Based Contrast Enhancement

The abovementioned measure of contrast is global, not local, i.e., the differences are computed regardless of the spatial locations of the pixels. This is not consistent with how we perceive contrast, which is in a localized manner, at each point having neighbors exert a higher influence than far-away points. Using the concepts introduced by the popular perceptually-based color correction method ACE of Rizzi et al. (2003), the authors of Bertalmío et al. (2007) propose an adapted version of the functional of Equation 1 that complies with some very basic visual perception principles, namely those of locality, color constancy and white patch (the latter stating that the brightest spot in the image is perceived as white, an observation that is often attributed, incorrectly, to the Retinex theory of Land (1977), but which has a long history that dates back at least to the works of Helmholtz, as explained by Judd (1979b,c)):

$\begin{array}{cc}\begin{array}{l}E(I)=\frac{\alpha }{2}{\displaystyle \sum _x{\left(I(x)-\frac{1}{2}\right)}^2}-\gamma {\displaystyle \sum _x{\displaystyle \sum _{y}w}}(x,y)|I(x)\\ -I(y)|+\frac{\beta }{2}{\displaystyle \sum _x(I(}x)-I_0(x){)}^2,\end{array}& (3)\end{array}$

where w is a distance function such that its value decreases as the distance between x and y increases, I₀ is the original image and α, β and γ are positive weights (which can be chosen so as to guarantee the white patch property, see Bertalmío et al., 2007 for details). The gradient descent equation for the functional in Equation 3 is the following, and its numerical implementation is essentially equivalent to the method of Rizzi et al. (2003):

$\begin{array}{cc}\begin{array}{l}{I}_t(x)=-\alpha \left(I(x)-\frac{1}{2}\right)+\gamma {\displaystyle \sum _{y}w}(x,y)sgn(I(x)-I(y))\\ -\beta (I(x)-I_0(x)).\end{array}& (4)\end{array}$

Starting from I = I₀, we iterate Equation 4 until we reach a steady state, which will be the result of this algorithm.

By minimizing the energy in Equation 3 we are locally enhancing contrast (second term) and promoting color constancy by discounting the illuminant (first term), while preventing the image from departing too much from its original values (third term). We could also say that the minimization of Equation 3 approximates local histogram equalization.

The method of Bertalmío et al. (2007) has several good properties:

1. It yields very good color constancy results, being able to remove strong color casts and to deal with non-uniform illumination (a challenging scenario for most color constancy algorithms, see Bertalmío, 2014).

2. It increases the dynamic range of the image (i.e., it tends to “flatten” its histogram).

3. It has a very good local contrast enhancement performance, producing results without halos, spurious colors or any other kind of visual artifact.

4. It can deal with both underexposed and overexposed pictures.

5. It reproduces visual perception phenomena such as simultaneous contrast and the Mach Band effect.

But regarding color constancy, there is also a very interesting and close connection with the classical approach of Retinex. In their kernel-based Retinex (KBR) formulation, Bertalmío et al. (2009) take all the essential elements of the Retinex theory of Land (1977) (channel independence, the ratio reset mechanism, local averages, non-linear correction) and propose an implementation that is intrinsically 2D, and therefore free of the issues associated with the 1D paths used in the original Retinex algorithm. The results obtained with this algorithm comply with all the expected properties of Retinex (such as performing color constancy while being unable to deal with overexposed images) but don't suffer from the usual shortcomings such as sensitivity to noise, appearance of halos, etc. In Bertalmío et al. (2009) it is proven that there isn't any energy that is minimized by the iterative application of the KBR algorithm, and this fact is linked to its limitations regarding overxposed pictures. Using the analysis of contrast performed by Palma-Amestoy et al. (2009), Bertalmío et al. (2009) are able to determine how to modify the basic KBR equation so that it can also handle overexposed images, and the resulting, modified KBR equation turns out to be essentially the gradient descent of the energy given by Equation 3. In other words, the method of Bertalmío et al. (2007) can be seen as an iterative application of Retinex, although in a modified version that allows to produce good results also in the case of overexposed images.

3. A New Neural Model

3.1. Connection with Neuroscience

The activity of a population of neurons in the region V1 of the visual cortex evolves in time according to the Wilson-Cowan equations (see Wilson and Cowan, 1972, 1973; Bressloff et al., 2002). Treating V1 as a planar sheet of nervous tissue, the state a(r, ϕ, t) of a population of cells with cortical space coordenates r ∈ ℝ² and orientation preference ϕ ∈ [0, π) can be modeled with the following PDE (Bressloff et al., 2002):

$\begin{array}{cc}\begin{array}{l}\frac{\partial a(r,\varphi ,t)}{\partial t}=-\alpha a(r,\varphi ,t)\\ +\mu {\displaystyle {\int }_0^{\pi }{\displaystyle {\int }_{{ℝ}^2}\omega }}(r,\varphi \Vert r^{\prime },{\varphi }^{\prime })\sigma (a(r^{\prime },{\varphi }^{\prime },t))d{r}^{\prime }d{\varphi }^{\prime }\\ +h(r,\varphi ,t),\end{array}& (5)\end{array}$

where α, μ are coupling coefficients, h(r, ϕ, t) is the external input (visual stimuli), ω(r, ϕ||r′, ϕ′) is a kernel that decays with the differences |r − r′|, |ϕ − ϕ′| and σ is a sigmoid function. If we ignore the orientation ϕ and assume that the input h is constant in time, it can be shown that Equation 5 is closely related to the gradient descent Equation 3, where neural activity a plays the role of image value I, sigmoid function σ behaves as the derivative of the absolute value function, and the visual input h is the initial image I₀. This connection was already pointed out by Bertalmío et al. (2007), and Bertalmío and Cowan (2009) use it to argue that the Wilson-Cowan equations could therefore be the gradient descent of a certain energy, and also that there would appear to be a physical substrate at the cortex for the Retinex theory.

3.2. Lightness Induction

Looking closely at Equation 4, we can see that the spatial arrangement of the image data plays no role in it. Therefore, we can expect that the local contrast enhancement procedure of Bertalmío et al. (2007) will always produce lightness contrast, not assimilation, since as we mentioned earlier assimilation is linked to high spatial frequencies (Shevell, 2003). Figure 3 confirms this: (Figure 3A) produces lightness assimilation, because all gray bars have the same value but they are perceived darker when surrounded by black and lighter when surrounded by white; on the other hand, the result produced by Bertalmío et al. (2007) in Figure 3B actually emulates lightness contrast rather than assimilation, as the line profiles in Figures 3C,D show.

FIGURE 3

Figure 3. (A) Original image, example of lightness assimilation: the gray bars have all the same value but appear different over black and white backgrounds. (B) Result of applying the model of Bertalmío et al. (2007) to image (A). (C) Profile of a line from image (A). (D) Profile of a line from image (B): notice how the model of Bertalmío et al. (2007) actually emulates lightness contrast rather than assimilation.

Rudd (2010) studies lightness induction using a disk-and-ring (DAR) display for matching experiments, see Figure 4 (left). The intensity of the background B, of the left ring R_M and of the disk on the right D_T is kept constant; the intensity of the right ring R_T is modified, and the observer has to adjust the intensity of the left disk D_M so as to match the appearance of the right disk D_T. Using the model of Bertalmío et al. (2007), the predicted value of D_M as a function of the varying R_T can be computed, and it is shown in Figure 4 (right). We can see that as R_T increases D_M always decreases, so according to this model we should only have lightness contrast in this situation. But the data from the perceptual experiments of Rudd (2010) says otherwise, see Figure 5: as R_T increases D_M also increases (lightness assimilation) until R_T reaches some value, beyond which D_M decreases (lightness contrast). These plots are well approximated by parabolas and, as the ring widths become larger, the resulting parabolas have their curvature decrease, implying that “assimilation is more likely to be observed with narrow surrounds” (Rudd, 2010).

FIGURE 4

Figure 4. Left: diagram of the disk-and-ring display used by Rudd (2010), taken from that paper. Right: value of D_M as a function of R_T, predicted with the model of Bertalmío et al. (2007).

FIGURE 5

Figure 5. Value of D_M as a function of R_T: results of perceptual matches for different ring widths, from Rudd (2010).

3.3. Proposed Model

In order to overcome the intrinsic limitations of Bertalmío et al. (2007) with respect to lightness induction, we should introduce spatial frequency in the energy functional. We propose a new model consisting in the following PDE, a modification of the gradient descent Equation (4):

$\begin{array}{cc}\begin{array}{l}{I}_t(x)=-\alpha (I(x)-\mu (x))+\gamma \left((1+{(\sigma (x))}^c\right){\displaystyle \sum _{y}w(x,y)sgn(I(x)}\\ -I(y))-\beta (I(x)-I_0(x)),\end{array}& (6)\end{array}$

where μ(x) is the mean average of the original image data computed over a neighborhood of x, σ(x) is the standard deviation of the image data computed over a small neighborhood of x, and the exponent c is a positive constant. The differences with respect to Equation 4 are that now the average in the first term is no longer global (the 1/2 value of Equation 4) but local, and that the weight for the second term is no longer a constant, but it changes both spatially and with each iteration, according to the local standard deviation σ: if the neighborhood over which it is computed is sufficiently small, standard deviation can provide a simple estimate of spatial frequency. But also, the standard deviation is commonly used in the vision literature as an estimate of local contrast. We have this contrast σ(x) raised to a power c, and this is also the case with other neural models where a power law is applied to the contrast, as we will briefly discuss later.

Again, this is a Wilson-Cowan type of neural activity model, where I₀ is the visual input. We take I₀ as a non-linear modification of the radiance stimulus, e.g., I₀ could be the result of applying the Naka-Rushton equation, which models photoreceptor responses (see Shapley and Enroth-Cugell, 1984), to the radiance stimuli. As we did with Equation 4, we start with an image I = I₀ and iterate Equation 6 until convergence, obtaining a result which we'll see is able to predict perceptual phenomena as well as improve the efficiency of the representation.

4. The Proposed Model Predicts Induction and Improves Efficiency

4.1. Predicting Lightness Induction

Using this new model, now we can qualitatively predict the results of Rudd (2010). We fix B, D_T, R_M and for each value of R_T we find the steady state of Equation 6 at the center of the right disk, at the middle of the right ring and at the middle of the left ring: these would be our model's predictions of the perceived values for D_T, R_T and R_M. Next, we compute the difference between the first two values and add it to the third, yielding the prediction of the perceived (lightness) value for D_M, from which we can recover the actual luminance value D_M using again Equation 6 (see Appendix for implementation details). From these results we can derive the following conclusions, that corroborate the findings of the perceptual experiments of Rudd (2010):

• As shown in Figure 6 (left), the predicted match luminance plots are no longer linear but quadratic, with an initial lightness assimilation regime for low values of R_T followed by a lightness contrast part.

• The curvature of these parabolas decreases with increasing ring width.

• The previous experiments were for a double-decrement display: B > R_T > D_T. If we now make the ring have a value in between disk and background, B > R_T < D_T, the plot curvature remains negative as in the double-decrement case, see Figure 6 (right). As reported by Rudd (2010), this behavior can not be predicted with some other neural models, like that of Rudd and Arrington (2001), with which the curvature changes sign in this situation.

FIGURE 6

Figure 6. Left: predicted value of D_M as a function of R_T, using the proposed model, for different ring widths; the plots are well approximated by parabolas, whose curvature decreases as the ring width increases. Right: prediction when the disk is a luminance increment with respect to the ring; the sign of the curvature remains negative.

We can also predict lightness assimilation in the previous example of the alternating gray bars of Figure 3, as we now show in Figure 7.

FIGURE 7

Figure 7. (A) Original image, example of lightness assimilation: the gray bars have all the same value but appear different over black and white backgrounds. (B) Result of applying the proposed model to image (A). (C) Profile of a line from image (A). (D) Profile of a line from image (B): notice how the proposed model is capable of emulating lightness assimilation.

It is interesting to note that the shape of the curves in Figure 6 does vary with extent of the neighborhood over which the standard deviation is computed, as Figure 8 shows: when the neighborhood covers disk, ring and some background we have an inverted parabola as before (red curve), but if we decrease the neighborhood size so that it only covers disk and ring but no background then the parabola concavity is reversed (green curve), and if the neighborhood is further reduced so that it only covers the disk then the curve is no longer parabolic but linear.

FIGURE 8

Figure 8. The shape of the curves that predict the value of D_M as a function of R_T, using the proposed model, depend on the extent of the neighborhood over which the standard deviation is computed. In red: the neighborhood covers disk, ring and some background. In green: the neighborhood covers disk and ring but no background. In blue: the neighborhood covers just the disk.

Finally, we may point out that some recent models which also predict lightness induction based on neural attributes of the visual system can be found in Otazu et al. (2008) and Penacchio et al. (2013).

4.2. Efficiency: Redundancy Reduction and Contrast Enhancement

In this section we argue that the proposed model of Equation 6, which as we have seen has the form of a Wilson-Cowan equation, performs local contrast enhancement and is closely related to basic image processing techniques, is a good candidate for a neural model providing the contrast constancy effects described by Georgeson and Sullivan (1975) and Martinez et al. (2014). But furthermore we will now see how this new neural model, applied to signals already encoded by photoreceptors, further improves efficiency by reducing redundancy: flattening the histogram and whitening the power spectrum.

Figure 9A shows a high dynamic range (HDR) image or radiance map, linearly scaled to the range [0,1]. Clearly this kind of mapping is useless, which is a way to explain the need for light adaptation and gain control mechanisms in our photoreceptors. Figure 9B shows the result of applying the Naka-Rushton equation to the previous HDR image. Figure 9C shows the result of applying our proposed model to Figure 9B. As we can see, the original radiance image has a very lopsided histogram, which is made considerably more uniform by applying the Naka-Rushton equation and even more flat if we apply our proposed method to the Naka-Rushton output. Local contrast is clearly enhanced as well, see for instance the window frames, the book cases behind the windows, etc. For the implementation details we refer the reader to the Appendix.

FIGURE 9

Figure 9. (A) High dynamic range (HDR) image, linearly scaled. (B) Result of applying the Naka-Rushton equation to the HDR image. (C) Result of applying proposed model to image (B). (D) Histogram of (A). (E) Histogram of (B). (F) Histogram of (C). Original image courtesy of Max Planck Institute.

Figure 10A shows the result of applying the Naka-Rushton equation to a high dynamic range image. Figure 10B shows the result of applying the model of Bertalmío et al. (2007) to Figure 10A (this is roughly equivalent to the tone mapping approach proposed by Ferradans et al. (2011) in an image processing/computer graphics context). Figure 10C shows the result of applying our proposed model to Figure 10A. Figure 10D compares the power spectrum of the three previous images. We can see that our model improves spectrum whitening over the other two results. In this image the contrast enhancement is more subtle but still noticeable, especially in the interior of the tree-trunk and on the leaves and grass in the foreground.

FIGURE 10

Figure 10. (A) Result of applying the Naka-Rushton equation to a high dynamic range image. (B) Result of applying the model of Bertalmío et al. (2007) to image (A). (C) Result of applying proposed model to image (A). (D) Power spectrum of images (A–C). Original image property of Industrial Light and Magic.

An interesting aspect is given by the constant c in Equation 6 and its relationship to the whitening of the power spectrum. Given a Naka-Rushton output, we compute the rotational average of its power spectrum which, in log-log coordinates, can be fit by a line with a certain slope. We do the same for the output of our model, that has been applied to the Naka-Rushton output using some value for c, and obtain a new linear fit with a new slope. Let us estimate the “increase in whitening” provided by our model as the difference between these two slopes, call it W. The value of W is a function of the constant c used in our model. If we now vary c in the interval [0, 1] we can plot the resulting function W(c), as shown in Figure 11. Disregarding the spikes for low values of c, we can see that there is an optimum value for c, with which we can obtain the maximum power spectrum whitening that our model can provide. In our model c is the power to which we raise the local standard deviation σ(x), and this standard deviation is one of the possible measures that are commonly used to estimate local contrast. Because of this reason we are currently investigating the possible connections of our model with the works of Mante et al. (2008) and Kay et al. (2013), since both of them apply a static power-law non-linearity to the contrast. In particular, Kay et al. (2013) computed the value of the exponent of this power law and found that while it varies accross the visual cortex, it is in the range [0, 0.35] with a value of around 1/3 in the case of V1: this is all consistent with the tests we have performed so far for different high dynamic range images.

FIGURE 11

Figure 11. Increase in spectrum whitening as a function of the constant c in Equation 6, for the image in Figure 10.

5. Conclusion and Future Work

In this paper we have proposed a neural model, in the form of a Wilson-Cowan equation, that is derived from an image processing technique for local histogram equalization. This new model is able to predict lightness induction phenomena, and improves the efficiency of the representation by flattening both the histogram and the power spectrum of the image signal and increasing local contrast.

We are very much interested in finding evidence of neural responses following our proposed model. Our method performs contrast enhancement, so we would like to explore whether there is any relationship with the work of Martinez et al. (2014), who have very recently confirmed that contrast enhancement takes place at the LGN and is much alike the common techniques used in image processing. Our model has a term where a power law is applied to the contrast, and we can optimize the exponent of this power law so as to maximize the whitening of the spectrum; for the limited tests that we have performed so far, our results appear to be in agreement with what is reported by Mante et al. (2008) and Kay et al. (2013), so we also want to investigate possible connections with those works. And as immediate future work, we will extend our formulation to the color case in order to predict color induction as well.

Last but not least, we believe we can use our proposed model to go back to some image processing and computer vision applications, which could benefit from the insights gained in the visual neuroscience domain. In particular, we are currently working in extending this new model for problems such as tone mapping, gamut mapping and computational color constancy.

Conflict of Interest Statement

The author declares that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

This work was supported by European Research Council, Starting Grant ref. 306337, by the ICREA Acadèemia Award and by Spanish grants ref. TIN2011-15954-E and ref. TIN2012-38112.

References

Atick, J. (1992). Could information theory provide an ecological theory of sensory processing? Network 3, 213–251. doi: 10.1088/0954-898X/3/2/009

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Barlow, H. B. (1961). Possible principles underlying the transformation of sensory messages. Sens. Commun. 217–234.

Bertalmío, M. (2014). Image Processing for Cinema. Boca Raton, FL: CRC Press; Taylor and Francis.

Bertalmío, M., Caselles, V., and Provenzi, E. (2009). Issues about retinex theory and contrast enhancement. Int. J. Comput. Vis. 83, 101–119. doi: 10.1007/s11263-009-0221-5

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Bertalmío, M., Caselles, V., Provenzi, E., and Rizzi, A. (2007). Perceptual color correction through variational techniques. IEEE Trans. Image Process. 16, 1058–1072. doi: 10.1109/TIP.2007.891777

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Bertalmío, M., and Cowan, J. (2009). Implementing the Retinex algorithm with Wilson-Cowan equations. J. Physiol. 103, 69–72. doi: 10.1016/j.jphysparis.2009.05.001

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Bressloff, P., Cowan, J., Golubitsky, M., Thomas, P., and Wiener, M. (2002). What geometric visual hallucinations tell us about the visual cortex. Neural Comput. 14, 473–491. doi: 10.1162/089976602317250861

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Buchsbaum, G. (1980). A spatial processor model for object colour perception. J. Frankl. Instit. 310, 337–350. doi: 10.1016/0016-0032(80)90058-7

CrossRef Full Text

Dan, Y., Atick, J. J., and Reid, R. C. (1996). Efficient coding of natural scenes in the lateral geniculate nucleus: experimental test of a computational theory. J. Neurosci. 16, 3351–3362.

Pubmed Abstract | Pubmed Full Text

Ferradans, S., Bertalmío, M., Provenzi, E., and Caselles, V. (2011). An analysis of visual adaptation and contrast perception for tone mapping. IEEE Trans. Patt. Analy. Mach. Intelli. 33, 2002–2012. doi: 10.1109/TPAMI.2011.46

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Georgeson, M., and Sullivan, G. (1975). Contrast constancy: deblurring in human vision by spatial frequency channels. J. Physiol. 252, 627–656.

Pubmed Abstract | Pubmed Full Text

Judd, D. (1940). Hue, saturation, and lightness of surface colors with chromatic illumination. J. Opt. Soc. Am. 30, 2–32. doi: 10.1364/JOSA.30.000002

CrossRef Full Text

Judd, D. (1979a). “The unsolved problem of color perception,” in Contributions to Color Science, Vol. 545, ed D. L. MacAdam (Washington, DC: National Bureau of Standards), 516–522.

Judd, D. (1979b). “Appraisal of Land's work on two-primary color projections,” Contributions to Color Science, Vol. 545, ed D. L. MacAdam (Washington, DC: National Bureau of Standards), 471–486.

Pubmed Abstract | Pubmed Full Text

Judd, D. (1979c). “Color appearance,” Contributions to Color Science, Vol. 545, ed D. L. MacAdam (Washington, DC: National Bureau of Standards), 539–564.

Kay, K., Winawer, J., Mezer, A., and Wandell, B. (2013). Compressive spatial summation in human visual cortex. J. Neurophysiol. 110, 481–494. doi: 10.1152/jn.00105.2013

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Land, E. (1977). The Retinex theory of color vision. Sci. Am. 237, 108–128. doi: 10.1038/scientificamerican1277-108

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Laughlin, S. (1981). A simple coding procedure enhances a neurons information capacity. Z. Naturforsch 36:51.

Pubmed Abstract | Pubmed Full Text

Mante, V., Bonin, V., and Carandini, M. (2008). Functional mechanisms shaping lateral geniculate responses to artificial and natural stimuli. Neuron 58, 625–638. doi: 10.1016/j.neuron.2008.03.011

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Martinez, L., Molano-Mazón, M., Wang, X., Sommer, F., and Hirsch, J. (2014). Statistical wiring of thalamic receptive fields optimizes spatial sampling of the retinal image. Neuron 81, 943–956. doi: 10.1016/j.neuron.2013.12.014

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Olshausen, B. (2003). “Principles of image representation in visual cortex,” in The Visual Neurosciences, Vol. 2, eds L. M. Chalupa and J. S. Werner (Cambridge, MA: MIT Press), 1603–1615.

Olshausen, B., and Field, D. (2000). Vision and the coding of natural images. Am. Sci. 88, 238–245. doi: 10.1511/2000.3.238

CrossRef Full Text

Otazu, X., Vanrell, M., and Alejandro Parraga, C. (2008). Multiresolution wavelet framework models brightness induction effects. Vis. Res. 48, 733–751. doi: 10.1016/j.visres.2007.12.008

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Palma-Amestoy, R., Provenzi, E., Bertalmío, M., and Caselles, V. (2009). A perceptually inspired variational framework for color enhancement. IEEE Trans. Patt. Analy. Mach. Intell. 31, 458–474. doi: 10.1109/TPAMI.2008.86

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Penacchio, O., Otazu, X., and Dempere-Marco, L. (2013). A neurodynamical model of brightness induction in v1. PLoS ONE 8:e64086. doi: 10.1371/journal.pone.0064086

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Pereverzeva, M., and Murray, S. (2008). Neural activity in human v1 correlates with dynamic lightness induction. J. Vis. 8:8. doi: 10.1167/8.15.8

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Pratt, W. (2007). Digital Image Processing, 4th Edn. Hoboken, NJ: John Wiley & Sons. doi: 10.1002/0470097434

CrossRef Full Text

Rizzi, A., Gatta, C., and Marini, D. (2003). A new algorithm for unsupervised global and local color correction. Patt. Recognit. Lett. 24, 1663–1677. doi: 10.1016/S0167-8655(02)00323-9

CrossRef Full Text

Rudd, M. (2010). How attention and contrast gain control interact to regulate lightness contrast and assimilation: a computational neural model. J. Vision 10:40. doi: 10.1167/10.14.40

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Rudd, M., and Arrington, K. (2001). Darkness filling-in: a neural model of darkness induction. Vis. Res. 41, 3649–3662. doi: 10.1016/S0042-6989(01)00216-4

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Sapiro, G., and Caselles, V. (1997). Histogram modification via differential equations. J. Diff. Eq. 135, 238–266. doi: 10.1006/jdeq.1996.3237

CrossRef Full Text

Shapley, R., and Enroth-Cugell, C. (1984). Visual adaptation and retinal gain controls. Prog. Retin. Res. 3, 263–346. doi: 10.1016/0278-4327(84)90011-7

CrossRef Full Text

Shevell, S. K. (ed.). (2003). “Color appearance,” The Science of Color (Oxford, UK: Elsevier Science Ltd.), 149–190.

Wilson, H., and Cowan, J. (1972). Excitatory and inhibitory interactions in localized populations of model neurons. Biophys. J. 12, 1–24. doi: 10.1016/S0006-3495(72)86068-5

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Wilson, H., and Cowan, J. (1973). A mathematical theory of the functional dynamics of cortical and thalamic nervous tissue. Biol. Cybern. 13, 55–80.

Pubmed Abstract | Pubmed Full Text

Appendix

In this section we give the implementation details for the results reported in the paper.

For the examples in Figure 6 we computed explicitly the steady state solution of Equation 6, with α = β = γ = 1, the second term of the equation with a weight γ(1 + 5σ^c), and a compression constant of c = 1/3 (Figure 6, left) or c = 0.75 (Figure 6, right).

For the example in Figure 7 we have adapted the fast numerical implementation of Bertalmío et al. (2007), with a polynomial approximation of degree 7 for the sign function, time step Δt = 0.15 and the stopping condition being fulfilled when the difference between the images of the current and the previous iteration falls below 0.1%. The original image is of size 100 × 200 and the parameter values were: α = β = γ = 1, c = 1/3, stencil size for the computation of the standard deviation σ(x) : 21 × 21, effective radius of the locality kernel w(x, y) : 75, and of the neighborhood over which we compute the mean average μ(x) : 19.

For the examples in Figures 9, 10 we have also used the same fast numerical implementation of Bertalmío et al. (2007), where now the stopping condition is fulfilled when the difference between the images of the current and the previous iteration falls below 0.5%. The parameter values were: α = β = γ = 1, c = 1/3, stencil size for the computation of the standard deviation σ(x) : 3 × 3, effective radius of the locality kernel w(x, y), effective radius of the neighborhood over which we compute the mean average μ(x): 1/3 of the numbers of rows or columns, whichever is larger.