Spectral color constancy using a maximum entropy approach

S. Skaff and J. J. Clark

Vol. 28, No. 11 / November 2011 / J. Opt. Soc. Am. A

2385

Spectral color constancy using a maximum entropy approach Sandra Skaff1,* and James J. Clark2 2

1 Xerox Research Center Webster, Xerox Corporation, 800 Phillips Road, Webster, New York 14580, USA Centre for Intelligent Machines, McGill University, 3480 University Street, Montreal, Quebec H3A 2A7, Canada *Corresponding author: [email protected]

Received January 12, 2011; revised September 5, 2011; accepted September 8, 2011; posted September 28, 2011 (Doc. ID 140734); published October 31, 2011 This paper proposes a solution to the spectral color constancy problem. The method is based on a statistical model for the surface reflectance spectrum and applies a maximum entropy constraint. Unlike prior methods based on linear models, the solution process does not require a set of basis functions to be defined, nor does it require a database of spectra to be specified in advance. Experiments on simulated and real data show that spectral estimation using the maximum entropy approach is feasible and performs similarly to existing spectral methods in spite of the lower level of a priori information required. © 2011 Optical Society of America OCIS codes: 330.1720, 330.1690, 330.4060.

1. INTRODUCTION Color constancy is the ability of the human visual system to perceive a surface as having a given color independently of the color of the light illuminating the surface. The soughtfor constancy of color is fundamentally predicated on the existence of an invariant property of surfaces, that being the surface’s reflectance function. Therefore, following [1], we approach the color constancy problem by estimating this invariant property directly rather than compensating for the effect of the illuminant on an arbitrary set of measurements. Many of the color constancy algorithms found in the literature take the other road and do not model the spectrum of the incoming light. Instead they are based on simple 3D color space models involving various transformations of the raw sensor signals (e.g., RGB, CMY, YIQ, etc.). However, it has been noted that the most flexible and robust color constancy methods are based on spectral models [2]. Basing color constancy on reflectance spectra instead of on 3D projections onto sensor response curves constitutes a big step toward the resolution of problems arising from the discrepancies in color spaces and sensor and display gamuts in different contexts. For example, capturing and consequently rendering color on displays of different devices such as cameras, computer monitors, iPods, and cell phones can pose color matching problems. Problems may also arise in identifying color printed on different media such as banner and sign material, photographic paper, and coated paper. Finally, it is worth noting that in computer graphics, researchers have suggested that the use of spectra would provide more accurate color rendering [3,4]. Spectral models are also advantageous in that they provide a common structure for integration of multiple sources of information. This information, for example, could arise from sensor responses arising from multiple surfaces that constitute a given scene, or could come from sensors of different spectral sensitivity profiles viewing a given surface [5]. 1084-7529/11/112385-15$15.00/0

Most of the spectral color constancy algorithms in the literature follow the seminal paper of Maloney and Wandell [6], who suggested utilizing linear models for the illuminant spectra and the wavelength dependence of the surface reflectance. These model a spectrum with a weighted sum of predefined spectral basis functions. A major difficulty in applying linear model representations of spectra is that the basis functions need to be specified in advance. In a few cases, the bases used to model reflectance spectra are polynomial [7] or Fourier [8] ones. In other and more common cases, the basis functions are obtained by performing principal components analysis (PCA) on a representative database of spectra. However, such a database is not always available. Even if such a database were available, it might not be consistent with, or representative of, the data used in a given application. In this paper we propose a spectral color constancy method that is not based on a linear model and, as such, does not require specification of a set of basis functions. Instead, we use the Bayesian approach of Brainard and Freeman [2,9] with a maximum entropy prior. This paper is organized as follows. First, an overview of color constancy approaches that represent surfaces and illuminants by spectral models is provided in Section 2. Next, the paper contributes a derivation of statistical models for surface and illuminant spectra employed in the proposed maximum entropy color constancy approach in Section 3. Results of estimating surface and illuminant spectra, in both simulation and experiment, are provided in Section 4. The proposed approach is compared to state-of-the-art approaches, both spectral- and 3D-based, in Section 5.

2. OVERVIEW OF SPECTRAL COLOR CONSTANCY Numerous algorithms that compute surface and illuminant spectra have been developed in the literature. None of these algorithms function perfectly in all situations, and most fail outside of a limited domain of applicability. The difficulty © 2011 Optical Society of America

2386

J. Opt. Soc. Am. A / Vol. 28, No. 11 / November 2011

S. Skaff and J. J. Clark

in solving the color constancy problem lies in the bilinear nature of the underlying relation between the illuminant and surface reflectance spectra and the measurements provided by the camera sensor responses. To explain the spectral color constancy problem, we start by describing the spectrum of light arriving at location x on an array of sensors by I x ðλÞ ¼ EðλÞS x ðλÞ;

ð1Þ

where λ denotes the wavelength taken over the visible range, assumed to take on discrete values, and x denotes the location in the image. S x ðλÞ is the surface spectral reflectance, which is the amount of light, relative to the incident light, reflected from a surface at each wavelength. EðλÞ is the illuminant spectral power distribution, which is the amount of power emitted by the illuminant and incident on the surface at each wavelength. I x ðλÞ is referred to as the product spectrum. Assuming that there are p sensors and the relative wavelength sensitivity of the kth sensor is Rk ðλÞ, the response recorded at each location x is given by P xk ¼

M X

S x ðλÞEðλÞRk ðλÞ;

k ¼ 1; 2; …; p;

ð2Þ

λ¼1

where M is the dimension of the spectra. Note that Eq. (2) holds only for flat and uniformly illuminated Lambertian surfaces, without fluorescence. In previous spectral color constancy approaches of the literature, the models used for S and E are linear with n and m basis functions, respectively: S x ðλÞ ¼

n X j¼1

EðλÞ ¼

m X

σ xj S j ðλÞ;

ð3Þ

ϵi Ei ðλÞ:

ð4Þ

i¼1

The goal in these approaches is to compute the surface and illuminant linear model weights σ xj , j ¼ 1; …; n and ϵi , i ¼ 1; …; m respectively. Linear model representations appear in Brill and West’s work [10–14] on von Kries adaptation [15], and Judd [16] had first hypothesized that this adaptation can be used to attain color constancy. The von Kries adaptation assumes that the three sensor types (long, medium, short) in the human retina react independently to discount the effect of changing illumination. In [10], it is proven that the ratio of products of reflectance, illuminant, and photoreceptor sensitivity spectra for two surfaces are illumination invariant when three classes of colors are considered. It is assumed that the illuminant spectrum is a linear combination of three Fourier basis functions. No constraints were imposed on the number of surface spectral basis functions, which are also considered to be Fourier. It is then shown in [11] that the same hypothesis holds, except that the assumptions on the number of basis functions for each of the surface and illuminant spectra are reversed. In [12], West discussed the limitations of color constancy as carried out by a device mimicking the human visual system. In [13], the authors presented spectral classes of illuminant and surface reflectances for which color constancy would

be possible given the von Kries model. These are sufficient but not necessary conditions for color constancy. In [14], West and Brill presented necessary and sufficient conditions for color constancy to be possible while representing illuminant spectra by linear models. In later work, Worthey [17] studied the limitations of color constancy by examining what features of the illuminant spectrum the eye is best able to discount. Linear model representations also appear in Buchsbaum’s [18] and Gershon et al.’s work [19] proposing variants of the Gray World algorithm, which assumed that the average of the surface reflectances in the scene is gray. In these methods, first the illuminant is estimated. Then the estimated illuminant is used to find the weights of the basis functions of the surface spectrum linear model. Buchsbaum chose to use three basis functions to represent both the illuminant and reflectance spectra following a study performed by Wyszecki and Stiles [20], who showed that such spectra can be represented with third order basis sets. Gershon et al. used three basis functions in the illuminant linear model representations, relying on a study performed by Judd et al., who showed that models with three basis functions provide approximations that are very close to the actual daylight illuminant spectra [21]. They also used three basis functions for the surface reflectance spectra relying on a study performed by Cohen, who showed that the first three characteristic vectors account for 99% of the variance in 150 Munsell chips [22]. Buchsbaum’s choice of basis functions was ad hoc, while the basis functions used in Gershon et al.’s algorithm were the principal components of real sets of illuminant and surface reflectance spectra. However, the assumption that the average of the surface reflectances in a scene is gray can be easily violated; therefore the Gray World approaches do not perform well in practice. Linear model representations also appear in Maloney and Wandell’s [6,23,24] and Yuille’s work [25]. The latter work is similar to the former except that Yuille generalizes the approach to account for cases when the illumination varies through space and when objects are three-dimensional as opposed to flat. Maloney and Wandell determined the plane spanning the sensor responses, and the knowledge of this plane led to recovering the illuminant spectra basis function weights. They then suggest that the surface basis function weights can be consequently computed. In both approaches, the assumption of having more sensor responses than the number of surface spectra basis functions is required. This is a serious limitation, especially in estimating surface spectra. It implies that no more than three basis functions can be used in the event that the rod photoreceptor response is used in addition to those of the cones. In addition, further constraints need to be imposed for the uniqueness of the solutions in both approaches. In [26–30], linear model representations are also employed. In [26] D’Zmura and Lennie developed an algorithm that discounted the illuminant effect relying on light adaptation and making use of eye movements. Three basis functions for each of the reflectance and illuminant spectra were recovered. In [27] D’Zmura proposed an algorithm, which requires three surfaces and two views for each surface with the assumption of uniform illumination, to recover three basis functions for the linear model representations of each of the surface and illuminant spectra in a scene. In [28,29] D’Zmura and Iverson proposed a model check algorithm to determine necessary and

S. Skaff and J. J. Clark

sufficient conditions for unique recovery of surface and illuminant spectra. They analyzed different combinations of the number of surface and illuminant spectral basis functions, the number of surfaces in a scene, the number of views per surface, and the number of sensor responses. They showed that unique recovery is possible only for a limited number of combinations of the different parameters. They then developed an algorithm that checks whether a given bilinear model provides unique recovery given the different parameters [30]. Another approach is Ho et al.’s, which estimates surface and illuminant spectra using linear model representations from the color signal as compared with the sensor responses [31]. This signal is the product of the illuminant and surface spectra and it needs to be measured using a spectroradiometer if the approach were to be applied on real data. The computations were carried out with three basis functions for each of the surface and illuminant spectra, even though the authors mention that more basis functions could be used.

3. MAXIMUM ENTROPY SPECTRAL-BASED COLOR CONSTANCY In the spectral methods described above, determination of the correct number of spectral basis functions versus sensor responses is required. Moreover, in all these approaches, either the surface spectral reflectance or the illuminant spectrum has at most three basis functions in its linear model representation. Such a small number of basis functions may not be sufficient for an accurate representation of a surface or illuminant spectrum. These issues led Brainard and Freeman to introduce their Bayesian color constancy approach [2,9]. This approach used a linear model, but imposed no restrictions on the number of surface and illuminant spectra basis functions compared to the number of sensor responses. It uses a prior model of the weights to resolve the ambiguity introduced by having more basis functions than measurements. The Bayesian approach is extended in [5,32,33], in which it is shown how surface and illuminant spectra facilitate spectral fusion by gathering sensor responses from different types of sensors. The Brainard and Freeman approach still has the problem of requiring the specification of a set of basis functions. This is easily remedied, however, by modifying the form of the prior model being used in the Bayesian formulation. There is no need to limit the prior to a linear model. A nonlinear model can be used. The approach we take in this paper is to use an (log) entropic prior. The basic idea is to take as the solution the spectrum that jointly minimizes the negative of the entropy of the spectrum considered as a probability distribution and the distance between the actual sensor values and the model sensor values. With regards to the problem of modeling surface reflectance and illuminant spectra, the question arises as to whether such spectra are typically found to have high entropies. Jaynes pointed out that a physical quantity frequently observed in practice will tend to a value that can be produced in the largest number of ways [34]. In the case of physical processes generating spectra, many surfaces observed in our everyday-life surroundings have spectra that have high entropy, as compared with monochromatic surfaces that have low entropy spectra [35]. The case for illuminants is less certain. Applying Jaynes’ argument, we would expect that illuminants observed in our everyday-life surroundings have relatively

Vol. 28, No. 11 / November 2011 / J. Opt. Soc. Am. A

2387

high entropy spectra, which therefore can be produced in the largest number of ways. However, many illuminants encountered in man-made environments are designed to use processes that rarely occur in nature and can have spectra that have quite low entropy. For example, LED or laser light, or the light from fluorescent bulbs, typically have very spiky spectra with low entropy. On the other hand, black body illumination such as sunlight or the light from incandescent bulbs have relatively high entropy. So we expect that our maximum entropy method will perform better under black body illuminants than under illuminants such as LED or fluorescent lighting. Stiles and Wyszecki [36] carried out an investigation on the number of metamers corresponding to given tristimulus values. They assumed a uniform sampling of spectral reflectances and showed both analytically and empirically that this method of sampling gives rise to a distribution of tristimulus values that are peaked in the center of the range. Such peaks correspond to colors with low saturation, which are shown by Clark and Skaff [1] to correspond to high spectral entropy. Thus the Stiles and Wyszecki metamer counting work supports our choice of maximum entropy priors. Morovic and Finlayson [37] suggest applying a Gaussian model to the distribution of linear model weights that define elements of a metamer set and choosing the single metamer that has the maximum probability. Such a choice will be similar, if not identical to, the maximum entropy metamer, for the same reason that underlies the result found by the Stiles and Wyszecki study, namely that frequently encountered spectra (the most probable ones) will be of high entropy. It should be noted that maximum entropy methods have been used previously in the field of color science for many different purposes than that described here. For example, the authors in [38] applied the maximum entropy approach to detect edges in multispectral imagery. Moreover, in [39,40], the authors identified the direction of lighting to perform shadow removal by taking it to be that corresponding to the minimum entropy. In [41], the authors showed how the maximum entropy principle is used in color conversions. In addition, the authors in [42] used maximum entropy models to detect skin in images. Finally, Marin-Franch and Foster computed the entropies of the distributions of the color-code values of the image pixels under two illuminants and consequently the mutual information between these entropies. The goal of this work was to determine how much color information is preserved in surfaces when using different color spaces, different colordifference formulae, or different spectral sharpenings [43]. A. Statistical Spectral Modeling Our approach builds on the work of [35], in which a statistical model of spectra was developed and where maximum entropy models were successfully used to estimate spectra given only photoreceptor responses. However, that work did not consider the problem of estimating the surface reflectance spectrum when the scene illuminant is not white. Entropies are generally defined only for probability density functions (pdf’s), which sum to one. It is worth noting, however, that in signal and image analysis applications, entropy measures are often computed for quantities that are not necessarily pdf’s, but that are normalized to sum to one [44]. Typically these quantities are normalized histograms. In the case of spectral color constancy, the illuminant spectrum can be

2388

J. Opt. Soc. Am. A / Vol. 28, No. 11 / November 2011

S. Skaff and J. J. Clark

considered as a histogram by counting the number of photons emitted from the light source in various wavelength bins. The histogram of the surface spectrum can be obtained by counting the number of photons reflected from a surface in various wavelength bins when the surface is illuminated by white light. In our work, however, we choose to model both the surface and illuminant spectra with statistical representations which identify them with physical quantities. The statistical model of a spectrum is obtained by considering single photons selected from a flux of photons with various wavelengths. We can identify the spectrum of the light to be the probability that a given photon in the flux has a specific wavelength. We will refer to this probability as pðλÞ. The development of a statistical model for a surface reflectance spectrum is more complicated, however, since surfaces act to transform the relative probabilities of the wavelengths of photons in the incident and reflected fluxes. Suppose that a photon source with an associated pdf pE ðλÞ and intensity β is shining on a colored surface as in Fig. 1. The spectrum of the light source is therefore EðλÞ ¼ βpE ðλÞ:

ð5Þ

A photon incident on the surface can either be absorbed, transmitted, or reflected. The fate of a photon after interaction with the surface is denoted with the symbol Q, where states Q ¼ A, Q ¼ T, and Q ¼ R indicate that the photon is absorbed, transmitted, or reflected, respectively. We assume that the wavelength of a photon does not change under reflection or transmission. The joint probability of a given photon after interaction with the surface having wavelength λ and in state Q will be denoted by pS ðλ; QÞ where pS ðλ; Q ¼ RÞ þ pS ðλ; Q ¼ AÞ þ pS ðλ; Q ¼ TÞ ¼ pE ðλÞ:

ð6Þ

The conditional probability of a photon having wavelength λ after interacting with the surface given that it is in state Q is pS ðλjQÞ ¼

pS ðλ; QÞ : pQ ðQÞ

ð7Þ

Using Bayes rule [45] this probability can be expressed as pS ðλjQÞ ¼

Q=R

pS ðQjλÞpE ðλÞ : pQ ðQÞ

ð8Þ

SðλÞ ¼ LðλjQ ¼ RÞ;

ð9Þ

LðλjQ ¼ RÞ ¼ pS ðQ ¼ RjλÞ

ð10Þ

where

is the likelihood of a photon having a wavelength λ given that it was reflected. This likelihood depends only on the properties of the surface. Equations (9) and (10) imply SðλÞ ¼ pS ðQ ¼ RjλÞ:

ð11Þ

In general, the probability that a given incident photon from a light source, with arbitrary associated pdf, is reflected is given by pQ ðQ ¼ RÞ ¼

M X

pS ðQ ¼ R; λÞ

ð12aÞ

pS ðQ ¼ RjλÞpE ðλÞ

ð12bÞ

SðλÞpE ðλÞ;

ð12cÞ

λ¼1

¼

M X λ¼1

¼

M X λ¼1

where M is the number of wavelength bins obtained by discretizing the visible wavelength range, which is the range over which photons can be detected. Equation (12c) is obtained from Eq. (12b) by substituting for the probability pS ðQ ¼ RjλÞ using Eq. (11). We now define the albedo, α, of the surface to be α ¼ pQ ðQ ¼ RÞ:

ð13Þ

Note that the albedo in this definition depends both on the surface reflectance as well as the spectrum of the incident light. If the incident light has all of its photons concentrated in a wavelength range for which the surface is not very reflective, then the albedo will be low. Conversely, if the incident light has all of its photons concentrated in wavelengths where the surface is reflective, then the albedo will be high. The pdf associated with the reflected light is therefore pR ðλÞ ¼ pS ðλjQ ¼ RÞ

ð14aÞ

¼

pS ðQ ¼ RjλÞpE ðλÞ pQ ðQ ¼ RÞ

ð14bÞ

¼

SðλÞ p ðλÞ: α E

ð14cÞ

w w w

w

Next, we define the surface reflectance function SðλÞ as the likelihood function:

Q=A

Q=T

Fig. 1. (Color online) Gray circle represents a photon coming from the illuminant and incident on the surface. The light blue circles represent this photon after it has been either reflected from the surface (Q ¼ T), absorbed by the surface (Q ¼ A), or transmitted by the surface (Q ¼ T). The photon has wavelength λ.

Bayes rule is used to obtain Eq. (14b) from Eq. (14a). Equations (11) and (13) are used to substitute for pS ðQjλÞ and pQ ðQÞ in Eq. (14b) to obtain Eq. (14c). So far we have shown how each of the surface and illuminant spectra SðλÞ and EðλÞ can be associated with pS ðQ ¼ RjλÞ and pE ðλÞ, respectively. We denote by RðλÞ the reflected light spectrum. The intensity of the reflected light is the product of the incident light intensity and the albedo, α. The spectrum of the reflected light is therefore

S. Skaff and J. J. Clark

Vol. 28, No. 11 / November 2011 / J. Opt. Soc. Am. A

RðλÞ ¼ αβpR ðλÞ

ð15aÞ


SðλÞ pE ðλÞ ¼ αβ α

ð15bÞ

¼ SðλÞ½βpE ðλÞ


ð15cÞ

¼ SðλÞEðλÞ:

ð15dÞ

Substituting pR ðλÞ from Eq. (14c) into Eq. (15a) yields Eq. (15b). Equation (5) is used to deduce Eq. (15d) from Eq. (15c). Note that we now have statistical models for SðλÞ and EðλÞ. In the case of the former the model is a likelihood function and in the latter it is a pdf (scaled by the intensity of the incident illumination). Using these statistical models we can sensibly compute the entropy of SðλÞ and EðλÞ. The definition and computation of entropies of likelihood functions is discussed in [46,47]. B. Maximum Entropy Approach The proposed color constancy approach aims at recovering surface and illuminant spectra by representing them using maximum entropy models, given only camera sensor responses. We had shown in simulation and experiment that such estimation is feasible in [48]. For simplicity, we start by illustrating the approach with the case of one surface patch in a scene illuminated by a single light source in Subsection 3.B.1. Note that color constancy problems are typically addressed for cases where there are at least two surface patches in the scene. We then move on to generalizing our approach to the case of n surface patches in Subsection 3.B.2. 1. One Surface Patch As a reference, a scene comprised of n surface patches and one illuminant is shown in Fig. 2. Our goal is to estimate the surface and illuminant spectra given the constraints imposed by the sensor responses [Eq. (2)] by modeling each of the surface and illuminant spectrum by a maximum entropy approach. Therefore, we seek the joint entropy H [49] of pS and pE . Assuming independence of the surface with the incident light, this joint entropy is given by the sum of the individual entropies:

HðpS ; pE Þ ¼ HðpS Þ þ HðpE Þ ¼−

M X

pS ðλÞ log pS ðλÞ −

λ¼1

M X

pE ðλÞ log pE ðλÞ: ð16Þ

λ¼1

As derived in Subsection 3.A, pS ðλÞ denotes the normalized likelihood function representing the surface reflectance, while pE ðλÞ ¼ EðλÞ=β [following Eq. (5)] is the normalized illuminant spectrum. Note that pS ðλÞ and pE ðλÞ are independent since SðλÞ and EðλÞ are independent. The latter two spectra are independent from each other as they are physical properties of the surface and the illuminant, respectively. The sensor response corresponding to one surface [x ¼ 1 in Eq. (2)] is P ¼ ½ P 1 P 2 P 3
T . SðλÞ and EðλÞ can be estimated by maximizing H as given in Eq. (16) subject to the following constraints: SðλÞ ≥ 0;

λ ¼ 1; …; M;

ð17aÞ

EðλÞ ≥ 0;

λ ¼ 1; …; M;

ð17bÞ

SðλÞ ≤ 1;

λ ¼ 1; …; M;

ð18aÞ

EðλÞ ≤ 1;

λ ¼ 1; …; M;

ð18bÞ

P − RT ðS:  EÞ ¼ 0;

ð19Þ

where R in Eq. (19) is the sensor spectral sensitivity matrix of the camera, and “:” denotes element-wise multiplication. Equations (17a) and (17b) enforce positivity on the surface and illuminant spectra respectively, while Eqs. (18a) and (18b) enforce the constraints that the surface and illuminant spectra are less than one respectively. Note it is assumed that the albedo of the surface α and the light source intensity β are known. This is a typical assumption in spectral color constancy algorithms that aim at recovering only the shape of the spectra. 2. Generalization: n Surface Patches After illustrating the maximum entropy approach with a one surface patch scene, we move on to a scene of n surface patches, where n ≥ 2. The goal is to estimate the surface patch spectra and the illuminant spectrum by extending the maximum entropy approach described in Subsection 3.B.1. The entropy for the all the surface and illuminant spectra can be expressed as Hn ¼ −

M X λ¼1

−

M X λ¼1

Fig. 2. (Color online) Scene with n surface patches of reflectance spectra S 1 ðλÞ, S 2 ðλÞ, S 3 ðλÞ, … S n ðλÞ illuminated by one light source of spectrum EðλÞ.

2389

pS1 ðλÞ log pS1 ðλÞ − pSn ðλÞ log pSn ðλÞ −

M X λ¼1 M X

pS2 ðλÞ log pS2 ðλÞ… pE ðλÞ log pE ðλÞ:

ð20Þ

λ¼1

Note that the reflectance spectrum is a physical property of a surface and is independent of the reflectance spectra of other surfaces in the scene. Consequently, the normalized surface reflectance spectra are independent. The normalized reflectance spectra for the surfaces are denoted by pS1 ðλÞ, pS2 ðλÞ,

2390

J. Opt. Soc. Am. A / Vol. 28, No. 11 / November 2011

and pSn ðλÞ for the first, second, and nth surface, respectively. The constraints imposed by the camera sensor responses are given by [Eq. (2)]: P1 − RT ðS1 :  EÞ ¼ 0;

ð21aÞ

P2 − RT ðS2 :  EÞ ¼ 0;

ð21bÞ

Pn − RT ðSn :  EÞ ¼ 0;

ð21cÞ

where P1 ; P2 ; …; Pn are the sensor responses of the n surfaces. The patches are assumed to be uniform and therefore have the same camera sensor responses at all pixels. These responses can be expressed as P1 ¼ ½ P 11

P 12

P 13
T ;

ð22aÞ

¼ ½ P 21

P 22

P 23
T ;

ð22bÞ

Pn ¼ ½ P n1

P n2

P n3
T :

ð22cÞ

P2

In addition, a set of positivity constraints is imposed on the surface and illuminant spectra as such: S 1 ðλÞ ≥ 0;

λ ¼ 1; …; M;

ð23aÞ

S 2 ðλÞ ≥ 0;

λ ¼ 1; …; M;

ð23bÞ

S n ðλÞ ≥ 0;

λ ¼ 1; …; M;

ð23cÞ

EðλÞ ≥ 0;

λ ¼ 1; …; M:

ð23dÞ

It is also ensured that the spectra estimated are less than one: S 1 ðλÞ ≤ 1;

λ ¼ 1; …; M;

ð24aÞ

S 2 ðλÞ ≤ 1;

λ ¼ 1; …; M;

ð24bÞ

S n ðλÞ ≤ 1;

λ ¼ 1; …; M;

ð24cÞ

EðλÞ ≤ 1;

λ ¼ 1; …; M:

ð24dÞ

The goal is to estimate the surface and illuminant spectra by maximizing H n as given in Eq. (20) subject to the constraints given in Eqs. (21), (23), and (24). Equivalently to maximizing the entropy H n , the negative of H n is minimized subject to the constraints imposed by the sensor responses. The location of the mode of the negative of this entropy function is estimated by a standard nonlinear constrained optimization function, fmincon, from the MATLAB Optimization Toolbox, resulting in a set of estimated surface and light source spectra. The algorithm used in solving the optimization problem was the sequential quadratic programming algorithm. The parameters set are 105 on the number of iterations and 108 on the number of function evaluations. The optimization algorithm converged at similar cost function values with the computed sensor responses matching the measured ones. Since the performance of our approach is demonstrated on four surface patch scenes in all cases considered, we use it to compare execution time as follows. For simulation without adding

S. Skaff and J. J. Clark

noise to the measurements (Subsection 4.B), the execution time is as short as 1 min. In the case of simulation where noise is added to the measurements (Subsection 5), the execution time varies from 1 min to 30 minutes. Finally, in the case of experiment (Subsection 4.C), the runtime can be in the order of 2 hours. The run times are quoted for a 3:2 GHz Pentium 4 processor with 1 GB RAM. This long running time is due to the large optimization dimensionality, which increases with the increase in the number of surfaces in a scene. The spectral models are 31D each. Therefore, in the case of four surface patches and one illuminant in a scene, for example, the solution sought is a 155D vector. For this reason, we run the optimization algorithm with 10 different randomly chosen starting points for each of a sample set of scenes only in each case.

4. RESULTS A. Spectral Estimation Evaluation The performance of the algorithm in estimating surface and illuminant spectra in both simulation and experiment is measured using the root mean square error (RMSE). This error conveys an unbiased difference between spectra without taking the human visual system into account. Such an error measure is suitable for our work, the goal of which is to estimate spectra for use in different applications such as computer graphics, printing and rendering. The actual spectrum of a surface S A ðλÞ and the spectrum obtained from the proposed algorithm S M ðλÞ, referred to as a model spectrum, are each normalized to a maximum of one. This is a common practice in solving color constancy problems as intensity information is not sought. Denoting the normalized actual and model spectra of a surface by S 0A ðλÞ and S 0M ðλÞ respectively, the RMSE is computed as such: vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u M u1 X RMSE ¼ t ðS 0 ðλÞ − S 0M ðλÞÞ2 : M λ¼1 A

ð25Þ

The RMSE between the normalized actual spectrum E0A ðλÞ of an illuminant and the normalized model spectrum obtained from the proposed algorithm E0M ðλÞ is computed as such: vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u M u1 X RMSE ¼ t ðE 0 ðλÞ − E0M ðλÞÞ2 : M λ¼1 A

ð26Þ

B. Simulations The performance of the maximum entropy approach in estimating surface and illuminant spectra is tested on matte Munsell patches. These are representative of a range of 40 hues encountered in printing. Their corresponding spectra were obtained from [50]. The illuminant spectra used are those of daylight and skylight, measured by the University of Joensuu Color Group [50], and a set of tungsten light spectra. The set of tungsten light spectra is composed of 10 blackbody spectra with temperatures ranging from 2600 K to 3500 K, in steps of 100 K, obtained from the IES lighting handbook [51]. The surface and illuminant spectra are multiplied to obtain the spectra of the light falling on the sensor. To obtain the simulated sensor responses for the surface patches used, this

1. Results We plot the model and actual spectra for one surface patch and the illuminant in three cases for each of two types of scenes, one illuminated by skylight and one illuminated by tungsten light. The three cases considered are those where there are one, two, and four surface patches in a scene. The plots are shown in Figs. 3 and 4 respectively. Note that each of the surface and illuminant model spectra appears to comprise a product of three exponentials, one in each wavelength range (long, medium, and short). This means that our models are conforming with the product of exponentials form that is representative of a maximum entropy solution. As for the illuminants, since the skylight spectrum has a peak in the low-wavelength range, it is better estimated by our model than the tungsten light spectrum, which has no peaks or valleys. This can be observed through both the relative shapes of the model and actual spectra as well as the corresponding RMSEs. Another observation made from the plots is that the spectral models of both the surfaces and illuminants are better estimates of the actual spectra when there are four surface patches in the scene. This can be attributed to the fact that there is a larger number of sensor responses in the case of four surface patches than in the case of one or two patches. The additional sensor responses contribute more information for improved scene illuminant spectral estimation and, consequently, for improved surface spectral estimation. Next, the average RMSEs between the estimated and actual spectra over all 100 scenes for each of the surface patch repeated in five scenes and the illuminant are plotted in Fig. 5. The figure shows that the spectral estimates generally improve upon combining the sensor responses from multiple surfaces. For the case of the surface patches, the average RMSE drops from 0.2325 in the case of one surface patch in the scene to 0.1537 in the case of five surface patches in the scene. For the case of the illuminants, the average RMSE drops from 0.2300 in the case of one surface patch in the scene

(a)

2391

1 0.9

Normalized Spectral Reflectance

spectrum is multiplied by the sensor spectral sensitivity curves of a Panasonic WV-CP410 camera [Eq. (2)]. These sensitivity curves were obtained from the manufacturer. The wavelength range considered for these spectra is 400 to 700 nm, and is discretized into 10 nm bins. This yields a dimension of M ¼ 31 in Eq. (2) for each spectrum. Noise is added to the responses in Subsection 5.A, where the approach is compared to previous spectral color constancy methods. We analyze the performance of the approach in the case of n surface patches in the scene, where n ¼ 1, 2, 3, 4, 5. These scenes are Mondrians consisting of several overlapping, matte (Lambertian) patches. The light illuminating a Mondrian scene is assumed to be locally constant. This means that the spectral characteristics of the light vary slowly in space. The artificial scenes are constructed as follows. Twenty surface patch and illuminant pairs are chosen at random. Five 5-surface patch scenes are constructed from each of these pairs by selecting the four remaining patches randomly. Therefore, a total of 100 scenes are constructed. The selection of scenes was done in this fashion so as to allow evaluation of the repeatability of spectral estimation in the scenes. The model surface and illuminant spectra are computed for each scene. To evaluate the performance of our approach, we use the RMSE measure as described in Subsection 4.A.

Vol. 28, No. 11 / November 2011 / J. Opt. Soc. Am. A

0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 400

(b)

actual model, 1 surface model, 2 surfaces model, 4 surfaces 450

500

550 600 Wavelength nm

650

700

550 600 Wavelength nm

650

700

1 0.9 0.8

Normalized Radiant Power

S. Skaff and J. J. Clark

0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 400

actual model, 1 surface model, 2 surfaces model, 4 surfaces 450

500

Fig. 3. Model and actual spectra obtained in simulation for a Munsell patch scene with one, two, and four surface patches in the scene: (a) Munsell patch 2.5P 2:5=2 (RMSE ¼ 0:2022, one surface; RMSE ¼ 0:1516, two surfaces; RMSE ¼ 0:1304, four surfaces) illuminated with (b) skylight 20 (RMSE ¼ 0:1501, one surface; RMSE ¼ 0:1352, two surfaces; RMSE ¼ 0:1150, four surfaces).

to 0.1798 in the case of five surface patches. However, the lowest error obtained in the case of the illuminants is 0.1729 in the case of three surface patches. In addition, we can see that the average RMSE over both the surface and illuminant spectral estimates decreases steadily from 0.2312 in the one surface patch case to 0.1668 in the five surface patch case. To gain a better understanding of the comparative performance of the approach in surface and illuminant spectral estimation in the cases of different numbers of surfaces in a scene, the paired T-test [52] is run on the average RMS errors for the 100 scenes [depicted in Figs. 5(a) and 5(b)] for a 95% confidence level. The pairs for which the average RMS errors are employed are the (1) one surface patch case and five surface patch case and (2) two surface patch case and five surface patch case. As a result of this test, we can conclude that the improvements in the surface spectral estimates in the five surface patch case compared with the one (t ¼ 5:5235, p < 0:05) and two (t ¼ 3:8834, p < 0:05) surface patch cases are statistically significant. Next, for the illuminant spectra we can conclude that the improvement in the illuminant spectral estimates in (1) is statistically significant with (t ¼ 2:5675, p < 0:05). However, we cannot conclude that the improvement in the illuminant spectral estimates in (2) is statistically significant (t ¼ 0:3181, p < 0:05). This is attributed to the fact that the average RMS error of the illuminant spectral estimates over all 100 scenes decreases very marginally, which is from 0.1838 in the two surface patch case to 0.1798 in the five surface patch case. We should note that a similar calculation for the median RMS error shows that it increases

2392

J. Opt. Soc. Am. A / Vol. 28, No. 11 / November 2011 (a)

S. Skaff and J. J. Clark (a)

1

0.3

0.25

0.8 Average of the RMSE’s

Normalized Spectral Reflectance

0.9

0.7 0.6 0.5 0.4 0.3 actual model, 1 surface model, 2 surfaces model, 4 surfaces

0.2 0.1 0 400

(b)

450

500

550 600 Wavelength nm

650

0.15

0.1

0.05

0

700

(b)

1

0.2

1

2 3 4 Number of Surfaces in the Scene

5

1

2 3 4 Number of Surfaces in the Scene

5

1

2 3 4 Number of Surfaces in the Scene

5

0.3

0.9 0.25 Average of the RMSE’s

0.7 0.6 0.5 0.4 0.3 actual model, 1 surface model, 2 surfaces model, 4 surfaces

0.2 0.1 0 400

450

500

550 600 Wavelength nm

650

slightly, from 0.1666 in the two surface patch case to 0.1708 in the five surface patch case. We would like to note also that the average RMS errors for the illuminant spectral estimates were found to be less in 60% of the groups for the five surface patch case as compared with the two surface patch case. Next, we run the T-test under the same conditions as above on the averages of the average RMS errors for both the surface and illuminant spectral estimates as depicted in Fig. 5(c) for pairs (1) and (2). We can conclude that the improvements in spectral estimates are statistically significant for (1) with (t ¼ 3:9741, p < 0:05). However, we cannot conclude that the improvements in spectral estimates are statistically significant for (2) with (t ¼ 1:9707, p < 0:05). We would like to note though that the average RMS error decreases from 0.1900 in the two surface patch case to 0.1667 in the five surface patch case. In addition, the average RMS errors were found to be less in 70% of the groups for the five surface patch case as compared with the two surface patch case. The median statistic confirms these results, as the median of the RMS errors over all scenes decreases from 0.1813 in the two surface patch case to 0.1728 in the five surface patch case. Also, the median RMS errors were found to be less in 70% of the groups for the five surface patch case as compared with the two surface patch case. Finally, we would like to note that the optimization problem is difficult and is of high dimensionality as described in Subsection 3.B.2. Therefore convergence at the global minimum cannot be guaranteed. Finally, these results show that fusing sensor responses from multiple surfaces is useful in improving spectral estimation. The use of spectral models facilitates fusion of informa-

0.15

0.1

0

700

Fig. 4. Model and actual spectra obtained in simulation for a Munsell patch scene with one, two, and four surface patches in the scene: (a) Munsell patch 2.5Y 6=2 (RMSE ¼ 0:1813, one surface; RMSE ¼ 0:1416, two surfaces; RMSE ¼ 0:0976, four surfaces) illuminated with (b) tungsten light at temperature 2900 K (RMSE ¼ 0:2434, one surface; RMSE ¼ 0:1861, two surfaces; RMSE ¼ 0:1903, four surfaces).

0.2

0.05

(c)

0.3

0.25 Average of the RMSE’s

Normalized Radiant Power

0.8

0.2

0.15

0.1

0.05

0

Fig. 5. Average RMSEs between the actual and model spectra obtained in simulation over 100 scenes for the (a) Munsell surface patches repeated in a set of five scenes, (b) daylight, skylight, and tungsten illuminants, and (c) the combination of both the surface patches and illuminants.

tion as the approach is easily extended to incorporate more surface patches in a scene. The entropy terms of the spectra of the additional patches are added to the overall entropy term [Eq. (20)] and the sensor responses from the additional surface patches are added to the set of constraints [Eq. (21)]. C. Experiments Real scenes were constructed with two categories of surface patches: Munsell and construction paper. The latter are used to show that the approach presented in this work can estimate spectra of surfaces other than Munsell ones. 1. Munsell Patches Images of the 41 pages in the Munsell Book of Color [53] were captured with a Panasonic WV-CP410 camera. Samples of 50 × 50 pixels in size were obtained from each of the Munsell patches using a segmentation algorithm. The algorithm

S. Skaff and J. J. Clark

Vol. 28, No. 11 / November 2011 / J. Opt. Soc. Am. A

assumed prior knowledge of the location of the center of the possible patches and segmented the samples around the centers. Since we assume Mondrian scenes where the illumination is locally constant, the RGB (red, green, blue) responses of all the pixels in a segmented patch were averaged to obtain one RGB response per patch. The algorithm assumed prior knowledge of the RGB of a white patch. Therefore if the response of a sample was greater than that of a white patch, it was discarded based on the assumption that it did not correspond to a patch. We did not have the γ specification of the camera. Therefore, we assume linear images, which means that the pixel values are assumed to be gamma-corrected with a factor of 1 (γ ¼ 1) in the experiments. The algorithm was run on the same 100 scenes used in simulation. The surface and illuminant spectra were taken to be over the same wavelength range and with the same discretizations as in simulation. The only exception was the illumination for which a tungsten light bulb of temperature 2800 K was used. Since we assume flat scenes, where there is no inter-reflection, the average responses corresponding to the patches in the scene are inputted into our algorithm. Experimental Results The model and actual spectra for one surface patch and the illuminant from a scene are plotted in Fig. 6. The RMSEs are computed as described in Subsection 4.B. We can see from the figure that the approach provides a reasonable estimate for the surface patch spectrum shown when there are four surface patches in the scene. The RMSE (a)

1

Normalized Spectral Reflectance

0.9 0.8 0.7 0.6 0.5 0.4 0.3 actual model, 1 surface model, 2 surfaces model, 4 surfaces

0.2 0.1 0 400

(b)

450

500

550 600 Wavelength nm

650

700

1 0.9

Normalized Radiant Power

0.8 0.7 0.6 0.5 0.4 0.3 actual model, 1 surface model, 2 surfaces model, 4 surfaces

0.2 0.1 0 400

450

500

550 600 Wavelength nm

650

700

Fig. 6. Model and actual spectra obtained in experiment for a Munsell patch scene with one, two, and four surface patches in the scene: (a) Munsell patch 2.5Y 6=2 (RMSE ¼ 0:2058, one surface; RMSE ¼ 0:1960, two surfaces; RMSE ¼ 0:1099, four surfaces) illuminated with (b) tungsten light at 2900 K (RMSE ¼ 0:5896, one surface; RMSE ¼ 0:5192, two surfaces; RMSE ¼ 0:4361, four surfaces).

2393

decreases from 0.2058 in the one surface patch case to 0.1099 in the four surface patch case. However, the approach does not perform as well in estimating the illuminant spectrum. As mentioned in Subsection 4.B.1, each of the surface and illuminant spectra is in the form of a product of exponentials. Since a tungsten light spectrum has the form of a straight line, a product of exponentials model would not serve as its best approximation. The average RMSEs between the estimated and actual spectra over all 100 scenes for each of the surface patch repeated in five scenes and the illuminant are plotted in Fig. 7. The figure shows that the spectral estimates generally improve upon combining the sensor responses from multiple surfaces. For the case of the surface patches, the average RMSE drops from 0.1795 in the case of one surface patch in the scene to 0.1554 in the case of five surface patches in the scene. For the case of the illuminants, the average RMSE drops from 0.5182 in the case of one surface patch in the scene to 0.4460 in the case of five surface patches in the scene. In addition, we can see that the average RMSE over both the surface and illuminant spectral estimates decreases steadily from 0.3488 in the one surface patch case to 0.3007 in the five surface patch case. To gain a better understanding of the comparative performance of the approach in surface and illuminant spectral estimation in the cases of different numbers of surfaces in a scene, the paired T test [52] is run on the average RMS errors for the 100 scenes [depicted in Figs. 7(a) and 7(b)] for a 95% confidence level. The pairs for which the average RMS errors are employed are the (1) one surface patch case and five surface patch case and (2) two surface patch case and five surface patch case. As a result of this test, we cannot conclude that the improvements in the surface spectral estimates in the five surface patch case compared with the one (t ¼ 1:8793, p < 0:05) and two (t ¼ 1:1165, p < 0:05) surface patch cases are statistically significant. However, the average RMS error of the surface spectral estimates over all 100 scenes decreases from 0.1795 and 0.1655 in the one and two surface patch cases, respectively, to 0.1554 in the five surface patch case. Moreover, the average RMS errors for the surface spectral estimates were found to be less in 65% of the groups for the five surface patch case compared with the one surface patch case. This average error was also found to be less in 55% of the groups for the five surface patch case compared with the two surface patch case. Next, for the illuminant spectra we conclude that the improvements in the illuminant spectral estimates in both (1) and (2) are statistically significant with (t ¼ 14:3559, p < 0:05) and (t ¼ 10:0063, p < 0:05) respectively. Next, we run the T test under the same conditions as above on the averages of the average RMS errors for both the surface and illuminant spectral estimates as depicted in Fig. 7(c) for pairs (1) and (2). We can conclude that the improvements in spectral estimates are statistically significant for cases (1) with (t ¼ 5:4344) p < 0:05) and (2) with (t ¼ 4:2083, p < 0:05). Finally, these results show that fusing sensor responses from multiple surfaces is useful in improving spectral estimation. In addition, the use of spectral models facilitates this fusion as the approach is easily extended from the one surface patch case to the five surface patch case. The entropy terms of the spectra of the additional patches are added to the overall

2394

J. Opt. Soc. Am. A / Vol. 28, No. 11 / November 2011 (a)

S. Skaff and J. J. Clark (a)

0.2

Normalized Spectral Reflectance

Average of the RMSE’s

0.16 0.14 0.12 0.1 0.08 0.06 0.04

(b)

0.7 0.6 0.5 0.4 0.3 0.2

1

2 3 4 Number of Surfaces in the Scene

0 400

5

0.6

(b)

actual spectrum model spectrum 450

500

550 600 Wavelength nm

Normalized Spectral Reflectance

0.4

0.3

0.2

0.1

(c)

1

2 3 4 Number of Surfaces in the Scene

5

0.8 0.7 0.6 0.5 0.4 0.3 0.2

0 400

actual spectrum model spectrum 450

0.4

500

550 600 Wavelength nm

650

700

Fig. 8. Model and actual spectra of construction paper samples, (a) light blue (RMSE ¼ 0:1238) and (b) green (RMSE ¼ 0:1155).

0.35 Average of the RMSE’s

700

1

0.1

0

650

0.9

0.5 Average of the RMSE’s

0.8

0.1

0.02 0

1 0.9

0.18

0.3

the entropy of 3.4340 for a flat spectrum. These computed entropies correspond to logarithms of base e.

0.25 0.2 0.15 0.1 0.05 0

1

2 3 4 Number of Surfaces in the Scene

5

Fig. 7. Average RMSEs between the actual and model spectra obtained in experiment over 100 scenes for the (a) Munsell surface patches repeated in a set of five scenes, (b) tungsten illuminants, and (c) the combination of both the surface patches and illuminants.

entropy term [Eq. (20)], and the sensor responses from the additional surface patches are added to the set of constraints [Eq. (21)]. 2. Construction Paper We segment four samples of 30 × 30 pixels in size from an image of construction paper patches. We average the camera responses of all the pixels in a patch to obtain one RGB response per patch. We feed these averages into our algorithm to obtain the corresponding spectral estimates and thereby constructing a four surface patch scene. Surface spectral estimates for two patches are shown in Fig. 8. We can see that the approach provides good spectral estimates of construction paper, which is expected as construction paper spectra have entropies similar to those of the Munsell patches. The average entropy of the spectra of nine construction paper samples is 3.2818, while that of all the Munsell patches is 3.3632. These averages are very close to

D. Comparison between Simulation and Experimental Results We compare the results obtained in simulation and experiment using both the average and median of the RMSEs between the model and actual spectra as shown in Table 1 for the case of four surface patches in a scene. On the first row, we show the average followed by the median of the RMSEs for the surface patch model spectra in all scenes. We then show the same statistics for the illuminant model spectra, for the simulations case. The second and third rows show equivalent results for the set of simulated scenes using tungsten light, and for the set of real scenes respectively. Both the average and the median statistics allow us to draw the same conclusions. As shown by the table, the performance of the proposed approach is similar in simulation and experiment in estimating surface patch spectra. However, our approach does not do as well at estimating the illuminant Table 1. Average and Median of the RMSEs in Surface Patch and Illuminant Spectral Estimates in All Munsell Scenes for Simulation and Experimental Data Scene Type Simulation (all lights) Simulation (tungsten) Experiment (tungsten)

Surfaces Average

Surfaces Median

Illuminants Average

Illuminants Median

0.1704

0.1539

0.1747

0.1650

0.1800

0.1048

0.2197

0.2422

0.1652

0.1379

0.4577

0.4590

S. Skaff and J. J. Clark

Vol. 28, No. 11 / November 2011 / J. Opt. Soc. Am. A

spectra of tungsten light for real scenes as compared with simulated scenes. Moreover, our approach does better at estimating daylight and skylight spectra than tungsten light spectra. This is evident in the lower average and median RMSE for the illuminants in the case of simulation in the all-lights case as compared with the tungsten illuminant case. The improvement in performance is due to the product of exponentialsshaped spectrum given by the maximum entropy approach for each of the surface and illuminant. Such a model spectrum is more aligned in shape with daylight and skylight spectra, which have more peaks and valleys as compared with tungsten light spectra.

5. COMPARISONS TO PREVIOUS WORK The approach introduced in Section 3 is compared to two types of color constancy approaches in the case of four surface patch scenes. The first type are the spectral approaches in Subsection 5.A and the second type are the approaches that perform color correction and employ 3D color representations in Subsection 5.B. A. Spectral Approaches As mentioned earlier, all approaches employing spectral models use linear combinations of spectral basis functions to represent surface and illuminant spectra. The objective in these approaches is to compute the weights of the basis functions. The limitations of these approaches as discussed in Section 2 motivated Brainard and Freeman to introduce the Bayesian algorithm for color constancy [2]. This algorithm imposes no restrictions on the number of surface and illuminant basis functions, nor on the number of photoreceptor responses. This leads us to choose to compare our proposed maximum entropy approach to the Bayesian approach, which is described briefly below. 1. Bayesian Approach The Bayesian approach regularizes the problem of computing the weights for the basis functions of the surface and illuminant spectra [2]. The posterior distribution function is computed over the set of surface and illuminant spectral weights given the photoreceptor responses. The solution weights vector is the one that maximizes the posterior if the maximum a posteriori (MAP) rule is used to minimize the Bayesian expected loss function. Our choice of eight basis functions comes from a study performed by Parkkinen et al. [54], who showed that as many as eight basis functions are needed to achieve an accurate representation of a Munsell patch spectrum. If we denote the vectors of surface and illuminant spectral model weights by a and b, respectively, and the sensor responses by RGB, we can write the posterior function pða; bjRGBÞ as follows: pða; bjRGBÞ ¼

pðRGBja; bÞpða; bÞ pðRGBÞ

∝ pðRGBja; bÞpða; bÞ ¼ pðRGBja; bÞpðaÞpðbÞ:

ð27Þ

pðRGBja; bÞ is the likelihood function that models the relationship between the weights and the sensor responses. pðRGBÞ denotes the probability of the sensor responses

2395

and is a normalization term that does not affect the shape of the posterior distribution. pða; bÞ denotes the prior information on the model weights. The surface and illuminant spectra are assumed to be independent; therefore, the prior distributions over the corresponding basis function weights are statistically independent (pða; bÞ ¼ pðaÞpðbÞ). The posterior, likelihood, and prior functions in Eq. (27) are computed in the same way as in Brainard and Freeman’s paper [2] as described in the next section. The solution for the surface and illuminant spectra basis function weights is obtained by maximizing the posterior function. 2. Numerical Results The simulated sensor responses are computed in the same way as in Subsection 4.B. In this case, however, noise is added to these responses to simulate the real world settings. Possible sources of noise in the real world could be the flickering of the light source while capturing the images, dust on the camera lens, and electronic noise in the camera. The added noise is modeled by a normal distribution having a standard deviation of 5% of the entire range of simulated responses for each of the long-, medium-, and short-wavelength range channels. The likelihood is computed using the model predictions of the sensor responses. The surface spectra are represented by eight basis functions while the illuminant spectra are represented by five basis functions. These basis functions are obtained by performing PCA on the set of Munsell patch spectra. The prior distributions are assumed to be Gaussian with means and variances computed from given databases of surface and illuminant spectra. The approach is simulated for 100 scenes where the surface patches and illuminants are chosen at random. It is assumed that there are four surface patches in the scene. For each scene, both the average and the median of the RMSEs of all the surface patch and the illuminant model spectra are computed as shown in Table 2. Both the mean and the median statistics shown in the table allow us to draw the same conclusions. They show that the performance of our proposed approach is similar to that of the Bayesian approach in estimating surface spectra even though the former requires no a priori information for modeling spectra. On the other hand, the use of basis functions allows the Bayesian approach to better model the illuminant spectra. B. Color Correction Approaches We compare the proposed approach to the prominent color correction algorithms following the framework proposed by Barnard et al. [55,56]. The illuminant estimate is referred to as the sensor response to a white surface under a canonical Table 2. Average and Median of the RMSEs in Surface Patch and Illuminant Spectral Estimates in All Scenes for the Maximum Entropy Approach and the Bayesian Approacha Scene Component Surface Patches Illuminants a

Maximum Entropy Average

Maximum Entropy Median

Bayesian Average

Bayesian Median

0.2930

0.2498

0.3063

0.2969

0.3376

0.3365

0.2047

0.1678

Simulations were carried out in the presence of noise.

2396

J. Opt. Soc. Am. A / Vol. 28, No. 11 / November 2011

illuminant. The canonical illuminant is the one selected from the database of illuminant spectra for which the camera is most balanced. We would like the reader to note that there are other approaches in the literature that perform color correction on images, such as that of Gehler et al.’s [57], which are not included in this study. The first algorithm is based on the Gray World algorithm as introduced by Buchsbaum [18] and modified by Gershon et al. [19]. The first version of the algorithm is similar to Buchsbaum’s approach. It is called GW (Gray World) and it assumes that the scene average is identical to the camera response to a surface of 50% uniform reflectance under a canonical illuminant. The second version of this algorithm, DB-GW, is similar to Gershon et al.’s approach. The scene average is the camera response to a surface of reflectance that is the average of a database of reflectances under canonical illumination. The database represents the reflectances of surfaces that are likely to be observed in a scene. The illuminant RGB, which is equivalent to the camera response to white, is taken to be double the response to gray, again under canonical illumination. The second algorithm, scale-by-max, is a limiting case of one version of Land and McCann’s Retinex theory [58,59]. This theory is based on explaining human color constancy and has been extended and analyzed by other researchers as well [60–63]. The goal of the Retinex approach is to estimate the 3D descriptors corresponding to the long-, medium-, and short-wavelength ranges of the lightness of a surface by comparing to the surface which has maximum lightness. The method can be extended to computing red, green, and blue descriptors as opposed to the lightness ones by multiplying the spectrum of light arriving at the sensor by the corresponding photoreceptor sensitivity curves. As introduced by Barnard, scale-by-max takes the illuminant RGB estimate to be the maximum in each channel (red, green, blue) over all the camera responses to all the surfaces in a scene. The third algorithm is color by correlation proposed by Finlayson et al. [64,65] as an improvement to the color in perspective algorithm [66], which is discussed in the paragraph after the next one. In the color by correlation algorithm, a correlation matrix is built in which each column corresponds to a different training illuminant and each row corresponds to a different image chromaticity. The chromaticities are computed from the surface reflectances in the training set. The convex hull of the chromaticities corresponding to one illuminant is computed. All the entries of the chromaticity bins falling in the convex hull are set to one, while all other elements are set to zero. The correlation matrix is multiplied by a vector containing information about the image chromaticities. This vector has the same discretization as the rows in the correlation matrix. An entry of this vector is set to one if the chromaticity occurs in the image and to zero otherwise. The scene illuminant selected is the one corresponding to the maximum correlation between the correlation matrix and the vector of chromaticities. This algorithm is labeled as C-by-C-01. There is a second version of the color by correlation approach in which the correlation matrix contains the probabilities computed based on the values in the correlation matrix of the previous version. The algorithm as described in [64,65] works well for synthetic data, as noted by Barnard et al. in [56]. In order for the approach to work on real data, Barnard

S. Skaff and J. J. Clark

applies a Gaussian filter to values in the correlation matrix corresponding to each illuminant. This compensates for noise and smoothens out any holes in the distribution. The illuminant is chosen based on the posterior probability distribution, in three different ways. The MAP estimate corresponds to the illuminant with maximum posterior probability. The minimum mean squared error (MMSE) estimate corresponds to the mean of the chromaticities weighted by the posterior distribution. The maximum local mass (MLM) estimate corresponds to the maximum of the function that is the convolution of the posterior distribution with a Gaussian mask. The three variants of this second version of the color by correlation algorithm are denoted by C-by-C-01-MAP, C-by-C-01-MMSE, and C-by-C-01-MLM, respectively. After discussing GW, scale-by-max, and C-by-C-01, we describe the last two algorithms that we compare our approach to. These two algorithms are based on Forsyth’s gamut mapping approach [67], which generally comprises two steps. The first step is to form two convex sets of RGBs represented by their respective convex hulls. One set comprises all possible surface RGBs when illuminated under a canonical illuminant. The other set comprises the set of surface RGBs under the unknown illuminant, and therefore the convex hull is the set of observed RGBs in this case. The two hulls are a unique diagonal mapping of each other under the diagonal assumption of illumination change. The set of possible maps is computed and the solution diagonal map is sought. This map corresponds to the solution set, which corresponds to the set of the surface RGBs under the canonical illuminant. This is the second step of the algorithm. Finlayson added two ideas to Forsyth’s approach to introduce the color in perspective approach [66]. First, he showed how to perform gamut mapping in 2D (chromaticity) instead of 3D (RGB) space as this simplifies gamut mapping and renders it usable in the real world, where specularities and varying illumination are commonlyarising problems. Second, he suggested that the diagonal maps can be restricted to only those corresponding to expected illuminants. Since the illumination constraint is nonconvex, Barnard [68] introduced the assumption of the convex hull for this constraint (ECRULE). As for obtaining the solution set, which is the second step of the algorithm, there are three methods. These are applied to both 2D and 3D gamut mapping, respectively. The first method is choosing the solution set with the maximum volume [66,67], referred to as 2D gamut mapping—MV and 3D gamut mapping—MV. The second method is taking the average of the constraint set [68], referred to as 2D gamut mapping—AVE and 3D gamut mapping—AVE. The third method is averaging over the constraint set, which is nonconvex, using Monte Carlo integration [69]. 1. Numerical Results Our goal is to compare our algorithm to the approaches discussed in the previous section. This implies that the chosen data sets and priors needed for these approaches are chosen in ways that yield them comparable to the method proposed in this paper. First, the scenes used for testing are the same as the ones used in Subsection 4.C.1. They therefore consist of four surface patches illuminated by a single light source. Second, the database of light source spectra used for training in the color by correlation approach is the same as that used in

S. Skaff and J. J. Clark

Vol. 28, No. 11 / November 2011 / J. Opt. Soc. Am. A

constructing scenes in our approach. This is a set of daylight and skylight spectra, measured by the University of Joensuu Color Group [50], and a set of tungsten light spectra [51]. Next, the canonical illuminant selected from this set is such that the camera is most balanced, and is therefore tungsten light at 3500 K. Note that the canonical illuminant choice can be arbitrary according to Barnard et al. [56]. In addition, the set of RGB measurements used to obtain the convex hull in the 2D and 3D gamut mapping approaches are the computed measurements for all 1269 patches from the Munsell Book of Color [53] under the illumination of that scene. The reflectance spectra of these patches are used in the Gray World approach as well. All the spectra employed in the previous approaches are 31D and over the same wavelength range as in our algorithm. Our algorithm estimates both illuminant and surface spectra. However the previous approaches perform color correction after carrying out an illumination estimation step. Therefore we focus on comparing the performance of illuminant estimation of the different approaches. Since we are interested in the chromaticity of the illuminant, we choose the second error measure in Barnard’s framework [56] for comparison. To obtain this error measure, the chromaticities ðr; gÞ for each of the actual and model illuminants given their RGBs is first computed: ðr a ; ga Þ ¼ ðRa =ðRa þ Ga þ Ba Þ; ðr m ; gm Þ ¼ ðRm =ðRm þ Gm þ Bm Þ;

Ga =ðRa þ Ga þ Ba ÞÞ; Gm =ðRm þ Gm þ Bm ÞÞ; ð28Þ

where (Ra , Ga , Ba ) is the RGB of the actual illuminant and (Rm , Gm , Bm ) is the RGB of the model illuminant. The error E is the Euclidean distance between the two chromaticities: E¼

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðr m − r a Þ2 þ ðgm − ga Þ2 :

ð29Þ

The proposed algorithm estimates an illuminant spectral model, which we project onto the three sensitivity curves of the camera sensor to obtain the equivalent three-dimensional RGB measurement. From this RGB, the corresponding chromaticity coordinates and then E are computed for the illuminant of each scene. To measure the performance of the algorithm in estimating illuminant spectra for all 100 scenes, the median of the Euclidean distance errors is computed. The median, instead of the average, is computed as Hordley and Finlayson showed in [70] that the median statistic is more suitable for comparing the algorithms in [56]. Next, the median of the Euclidean distances is computed over the illuminant chromaticity estimates of the same 100 scenes as given by each of the five main previous approaches and their corresponding variations, where applicable, described at the beginning of Subsection 5.B. The training data used comprised our surface patch and illuminant sets, which are used in Subsection 4.C. The canonical illuminant was tungsten as used in our experiments. The median errors are shown in Table 3. 2. Discussion First, we would like to note that the 2D gamut mapping approach did not find a solution for 28 out of the 100 scenes. This may be due to the fact that the canonical gamut, formed

2397

Table 3. Median of the Euclidean Distances for the Chromaticities of the Illuminant Estimates of 100 Scenes (72 Scenes for 2D Gamut Mapping) Given by Five Main Previous Approaches with Their Corresponding Variations (When Applicable) and Our Maximum Entropy Approacha Algorithm 2D Gamut Mapping—MV 2D Gamut Mapping—ICA GW Scale-by-Max 2D Gamut Mapping—AVE DB-GW C-by-C-MLM 3D Gamut Mapping—ICA 3D Gamut Mapping—AVE C-by-C-MAP 3D Gamut Mapping—MV C-by-C-MMSE C-by-C-01 Maximum Entropy approach a

Median Euclidean Distance Error 0.0607 0.0554 0.0490 0.0478 0.0455 0.0436 0.0424 0.0423 0.0422 0.0420 0.0409 0.0406 0.0288 0.0310

Median errors of the previous approaches are listed in decreasing order.

by the RGBs of the surfaces viewed under the canonical illuminant, is limited as it does not take into account all possible surfaces. Only the surface patches from our data set are considered. The median Euclidean distance error given in the table for 2D gamut mapping is, therefore, over 72 scenes. From the table, it can be concluded that the Gray World and the scale-by-max approaches perform poorly, as expected, since the assumptions imposed by these approaches are severe and are not always valid for real world images. It is worth noting that recently Funt and Shi showed that the scale-bymax approach performs very well when run on high dynamic images [71]. The color by correlation algorithm with all its variations performs relatively well compared to the other four previous approaches, as the former does not impose as strong a set of assumptions as compared with other approaches. Hence the low median Euclidean distance errors in color by correlation. As obtained by Barnard et al., the MMSE estimate of the illuminant is the best among the variants of the second version of the color by correlation approach. C-by-C-01 yields the best results among all the approaches because the scene illuminant belongs to the set of training illuminants. In other words, the sought scene illuminant is one of the training illuminants and it is the set of training illuminants that is used to build the correlation matrix. Therefore the selection of the scene illuminant, based on the maximum correlation between the correlation matrix and the vector of chromaticities, can be very precise. In addition, C-by-C-01 does not impose as strong a set of assumptions as compared with other approaches. Finally, we conclude that the proposed algorithm yields illuminant spectral estimates that are better than most of the approaches considered, except for C-by-C-01. This can be explained by the fact that our approach has an advantage over the gamut mapping approaches in that it does not require a canonical gamut, which should be formed by taking into account as many surface patches as possible. Moreover, our approach does not impose any assumptions on the scene under study, which is the case in the Gray World and scale-by-max

2398

J. Opt. Soc. Am. A / Vol. 28, No. 11 / November 2011

algorithms. Therefore, our approach demonstrates flexibility of use.

6. CONCLUSIONS In this work, we have proposed a maximum entropy approach to spectral color constancy. The novelty in this approach lies in providing a solution to spectral color constancy such that only weak a priori information is required to model spectra. In previous methods, this a priori information is in the form of a pre-specified spectral database used to compute a set of basis functions that may be difficult to obtain or that may not be representative of the spectra being modeled. The performance of the introduced approach was shown to be successful in estimating primarily surface reflectance spectra in both simulation and experiment. Identifying colors with the surface reflectance function casts the color constancy problem in a different light than was done previously. In this approach, obtaining the surface reflectance spectra solves the color constancy problem directly as these spectra are by definition illumination invariant [1]. In order to compute a reflectance spectrum with a maximum entropy approach, a statistical model is needed. The appropriate model was shown to be a normalized likelihood function. The derivation was based on studying the interaction between the spectrum of incident light and the spectrum of light reflected from a surface. The experiments performed on the scenes considered in this paper show that the performance of the proposed approach in estimating surface spectra is similar to that of the best spectral color constancy approaches. It was also shown that our approach outperforms most of the methods that employ 3D models of color. However, the reader should note that the purpose of these comparisons was to show that spectralbased approaches can perform as well as 2D- and 3D-based approaches in solving for color constancy. In fact, the use of spectral models does not necessarily enhance color constancy [72]. The primary advantage of using spectral models lies in their wide applicability, as their use is not contingent on a particular color space. In addition, spectral models facilitate fusion of sensor responses from multiple surfaces of the same scene as was shown throughout this paper.

REFERENCES 1. J. J. Clark and S. Skaff, “A spectral theory of color perception,” J. Opt. Soc. Am. A 26, 2488–2502 (2009). 2. D. H. Brainard and W. T. Freeman, “Bayesian color constancy,” J. Opt. Soc. Am. A 14, 1393–1411 (1997). 3. R. Hall, “Comparing spectral color computation methods,” IEEE Comp. Grap. Appl. 19, 36–45 (1999). 4. K. Devlin, A. Chalmers, A. Wilkie, and W. Purgathofer, “Tone reproduction and physically based spectral rendering,” in Computer Graphics Forum (Eurographics, 2002), pp. 101–123. 5. S. Skaff, T. Arbel, and J. J. Clark, “A sequential Bayesian approach to color constancy using multiple sensors,” Computer Vision and Image Understanding 113, 993–1004 (2009). 6. L. T. Maloney and B. A. Wandell, “Color constancy: a method for recovering surface spectral reflectance,” J. Opt. Soc. Am. A 3, 29–33 (1986). 7. P. Moon and D. E. Spencer, “Polynomial representation of reflectance curves,” J. Opt. Soc. Am. 35, 597–598 (1945). 8. R. Schettini, “Deriving spectral reflectance functions of computer-simulated object colors,” Comput. Graph. Forum 13, 211–217 (1994). 9. W. T. Freeman and D. H. Brainard, “Bayesian decision theory, the maximum local mass estimate, and color constancy,” in

S. Skaff and J. J. Clark

10. 11. 12. 13. 14. 15. 16. 17. 18. 19.

20. 21. 22. 23. 24. 25. 26. 27. 28. 29. 30. 31.

32.

33. 34. 35.

36. 37. 38.

Proceedings of the IEEE International Conference on Computer Vision (IEEE, 1995), pp. 210–217. M. H. Brill, “A device performing illuminant-invariant assessment of chromatic relations,” J. Theor. Biol. 71, 473–478 (1978). M. H. Brill, “Further features of the illuminant-invariant trichromatic photosensor,” J. Theor. Biol. 78, 305–308 (1979). G. West, “Color perception and the limits of color constancy,” J. Math. Biol. 8, 47–53 (1979). M. H. Brill and G. West, “Contributions to the theory of invariance of colour under the condition of varying illumination,” J. Math. Biol. 11, 337–350 (1981). G. West and M. H. Brill, “Necessary and sufficient conditions for Von Kries chromatic adaptation to give color constancy,” J. Math. Biol. 15, 249–258 (1982). J. von Kries, “Beitrag zur Physiologie der Gesichtsempfinding,” Arch. Anat. Physiol. 2, 5050–5524 (1878). D. B. Judd, “Hue saturation and lightness of surface colors with chromatic illumination,” J. Opt. Soc. Am. 30, 2–32 (1940). J. A. Worthey, “Limitations of color constancy,” J. Opt. Soc. Am. A 2, 1014–1026 (1985). G. Buchsbaum, “A spatial processor model for object colour perception,” J. Franklin Inst. 310, 1–26 (1980). R. Gershon, A. D. Jepson, and J. K. Tsotsos, “From [R, G, B] to surface reflectance: computing color constant descriptors in images,” in Proceedings of the International Joint Conference on Artificial Intelligence (AAAI/MIT, 1987), pp. 755–758. G. Wyszecki and W. S. Stiles, Color Science: Concepts and Methods, Quantitative Data and Formulas (Wiley, 1967). D. B. Judd, D. L. MacAdam, and G. Wyszecky, “Spectral distribution of typical daylight as a function of correlated color temperature,” J. Opt. Soc. Am. 54, 1031–1040 (1964). J. Cohen, “Dependency of the spectral reflectance curves of the Munsell color chips,” Psychon. Sci. 1, 369–370 (1964). L. T. Maloney, “Computational approaches to color constancy,” Ph.D. thesis (Stanford University, 1984). B. A. Wandell, “The synthesis and analysis of color images,” IEEE Trans. Pattern Anal. Machine Intell. 9, 2–13 (1987). A. Yuille, “A method for computing spectral reflectance,” Biol. Cybern. 56, 195–201 (1987). M. D’Zmura and P. Lennie, “Mechanisms of color constancy,” J. Opt. Soc. Am. A 3, 1662–1672 (1986). M. D. D’Zmura, “Color constancy: surface color from changing illumination,” J. Opt. Soc. Am. A 9, 490–493 (1992). M. D. D’Zmura and G. Iverson, “Color Constancy II. Results for two-stage linear recovery of spectral descriptions for lights and surfaces,” J. Opt. Soc. Am. A 10, 2166–2180 (1993). M. D. D’Zmura and G. Iverson, “Color Constancy I. Basic theory of two-stage linear recovery of spectral descriptions for lights and surfaces,” J. Opt. Soc. Am. A 10, 2148–2165 (1993). M. D. D’Zmura and G. Iverson, “Color Constancy III. General linear recovery of spectral descriptions for lights and surfaces,” J. Opt. Soc. Am. A 11, 2389–2400 (1994). J. Ho, B. V. Funt, and M. S. Drew, “Separating a color signal into illumination and surface reflectance components: theory and applications,” IEEE Trans. Pattern Anal. Machine Intell. 12, 966–977 (1990). S. Skaff, T. Arbel, and J. J. Clark, “Active Bayesian color constancy with non-uniform sensors,” in Proceedings of the IEEE International Conference on Pattern Recognition (IEEE, 2002), pp. 681–684. S. Skaff, “Active Bayesian color constancy with non-uniform sensors,” master’s thesis (McGill University, 2002). E. T. Jaynes, “Prior probabilities,” IEEE Trans. Syst. Sci. Cyber. 4, 227–241 (1968). J. J. Clark and S. Skaff, “Maximum entropy models of surface reflectance spectra,” in Proceedings of the IEEE Instrumentation and Measurement Technology Conference (IEEE, 2005), pp. 1557–1560. W. S. Stiles and G. Wyszecki, “Counting metameric object colors,” J. Opt. Soc. Am. 52, 313–322 (1962). P. Morovic and G. D. Finlayson, “Metamer-set-based approach to estimating surface reflectance from camera RGB,” J. Opt. Soc. Am. A 23, 1814–1822 (2006). N. H. Younan, K. Ponnala, and N. Alapati, “Edge detection in multispectral imagery via maximum entropy,” presented at the

S. Skaff and J. J. Clark

39. 40. 41. 42.

43.

44. 45. 46. 47. 48.

49. 50. 51. 52. 53. 54.

International Symposium on Remote Sensing of Environment, San José, Costa Rica (25–29 June 2007). G. D. Finlayson, M. S. Drew, and C. Lu, “Intrinsic images by entropy minimization,” in Proceedings of the European Conference on Computer Vision (Springer, 2004), pp. 582–595. C. Lu, “Removing shadows from color images,” Ph.D. thesis (School of Computing, Simon Fraser University, 2006). M. R. Gupta and R. M. Gray, “Color conversions using maximum entropy estimation,” in Proceedings of the IEEE International Conference on Image Processing (IEEE, 2001), pp. 118–121. B. Jedynak, H. Zheng, M. Daoudi, and D. Barret, “Maximum entropy models for skin detection,” in Proceedings of the IEEE International Workshop on Energy Minimization Methods in Computer Vision and Pattern Recognition (IEEE, 2003), pp. 180–193. I. Marin-Franch and D. H. Foster, “Estimating the information available from coloured surfaces in natural scenes,” in Proceedings of the European Conference on Color in Graphics, Imaging and Vision (The Society for Imaging Science and Technology, 2006), pp. 44–47. E. L. Hall, “TV and color coding,” in Computer Image Processing and Recognition, W. Rheinboldt, ed. (Academic, 1979), pp. 214–239. R. O. Duda, P. E. Hart, and D. G. Stork, Pattern Classification and Scene Analysis (Wiley, 2001). M. Hardy, “Entropies of likelihood functions,” in Maximum Entropy and Bayesian Methods, C. R. Smith, G. J. Erickson, and P. O. Neudorfer, eds. (Academic, 1991), pp. 127–130. M. Toews and T. Arbel, “Entropy-of-likelihood feature selection for image correspondence,” in Proceedings of the IEEE Conference on Computer Vision (IEEE, 2003), pp. 1041–1047. S. Skaff and J. J. Clark, “Maximum entropy spectral models for color constancy,” in Proceedings of the Color Imaging Conference (The Society for Imaging Science and Technology, 2007), pp. 100–105. T. M. Cover and J. A. Thomas, Elements of Information Theory (Wiley, 1991). University of Joensuu Color Group, “Spectral database,” http:// spectral.joensuu.fi. J. E. Kaufman, IES Lighting Handbook (Illuminating Engineering Society of North America, 1981). D. J. Sheskin, Handbook of Parametric and Nonparametric Statistical Procedures (Chapman and Hall, 2000). Munsell Color Corporation, Munsell Book of Color-Matte Finish Collection (Munsell Color, 1976). J. P. S. Parkkinen, J. Hallikainen, and T. Jaaskelainen, “Characteristic spectra of Munsell colors,” J. Opt. Soc. Am. A 6, 318–322 (1989).

Vol. 28, No. 11 / November 2011 / J. Opt. Soc. Am. A

2399

55. K. Barnard, “Practical colour constancy,” Ph.D. thesis (School of Computing, Simon Fraser University, 1999). 56. K. Barnard, V. Cardei, and B. Funt, “A comparison of computational color constancy algorithms. I.: methodology and experiments with synthesized data,” IEEE Trans. Image Process. 11, 972–984 (2002). 57. P. V. Gehler, C. Rother, A. Blake, T. Sharp, and T. Minka, “Bayesian color constancy revisited,” in Proceedings of the IEEE International Conference on Computer Vision and Patter Recognition (IEEE, 2008), pp. 1–8. 58. E. H. Land and J. J. McCann, “Lightness and Retinex theory,” J. Opt. Soc. Am. 61, 1–11 (1971). 59. E. H. Land, “The Retinex theory of color vision,” Sci. Am. 237, 108–128 (1977). 60. B. K. P. Horn, “Determining lightness from an image,” Comput. Vision Graph. 3, 277–299 (1974). 61. A. Blake, “Boundary conditions for lightness computation in Mondrian World,” Lect. Notes Comput. Sci. 32, 314–327 (1985). 62. D. H. Brainard and B. A. Wandell, “Analysis of the Retinex theory of color vision,” J. Opt. Soc. Am. A 3, 1651–1661 (1986). 63. A. Hurlbert, “Formal connections between lightness algorithms,” J. Opt. Soc. Am. A 3, 1684–1693 (1986). 64. G. D. Finlayson, P. M. Hubel, and S. Hordley, “Color by correlation,” in Proceedings of the Fifth Color Imaging Conference (The Society for Imaging Science and Technology, 1997), pp. 6–11. 65. G. D. Finlayson, S. D. Hordley, and P. M. Hubel, “Color by correlation: a simple unifying framework for color constancy,” IEEE Trans. Pattern Anal. Machine Intell. 23, 1209–1221 (2001). 66. G. D. Finlayson, “Color in perspective,” IEEE Trans. Pattern Anal. Machine Intell. 18, 1034–1038 (1996). 67. D. A. Forsyth, “A novel algorithm for color constancy,” Int. J. Comput. Vis. 5, 5–36 (1990). 68. K. Barnard, “Computational colour constancy: taking theory into practice,” Master’s thesis (School of Computing, Simon Fraser University, 1995). 69. G. Finlayson and S. Hordley, “Selection for gamut mapping colour constancy,” Image Vis. Comp. 17, 597–604 (1999). 70. S. D. Hordley and G. D. Finlayson, “Reevaluation of color constancy algorithm performance,” J. Opt. Soc. Am. A 23, 1008–1020 (2006). 71. B. Funt and L. Shi, “The Rehabilitation of MaxRGB,” in Proceedings of the Eighteenth Color Imaging Conference (The Society for Imaging Science and Technology, 2010), pp. 256–259. 72. M. Mosny and B. Funt, “Multispectral color constancy,” in Proceedings of the Color Imaging Conference (The Society for Imaging Science and Technology, 2006), pp. 309–313.