Hyperspectral imaging in color vision research : tutorial

This tutorial offers an introduction to terrestrial and close-range hyperspectral imaging and some of its uses in human color vision research. The main types of hyperspectral cameras are described together with procedures for image acquisition, postprocessing, and calibration for either radiance or reflectance data. Image transformations are defined for colorimetric representations, color rendering, and cone receptor and postreceptor coding. Several example applications are also presented. These include calculating the color properties of scenes, such as gamut volume and metamerism, and analyzing the utility of color in observer tasks, such as identifying surfaces under illuminant changes. The effects of noise and uncertainty are considered in both image acquisition and color vision applications.


INTRODUCTION
Hyperspectral imaging combines spatial and spectral data in such a way that each pixel of an image of a scene or object represents a continuous radiance or reflectance spectrum [1].Terminology in spectral imaging is not standardized [2], but can be indicative.Thus hyperspectral imaging differs from multispectral imaging, where spectral sampling takes place over multiple nonadjacent wavelength bands or is relatively coarse within a band so that some spectral information is lost.In turn, multispectral imaging differs from conventional trichromatic or RGB color imaging, where spectral sampling is reduced to just three components, from the long-, medium-, and shortwavelength regions of the visible spectrum.
Hyperspectral imaging is relevant to research in vision, especially vision in natural scenes, because hyperspectral images capture the full spatial and spectral properties of the light signal, thereby preserving the relevant degrees of freedom.In the main, these signals are spatially and spectrally complex [3]. Figure 1 shows radiance spectra reflected from individual pixels in a hyperspectral image of a flower scene.
Natural scenes are interpreted broadly to include both undeveloped and developed land cover [6,7].Specifically, scenes may be predominantly vegetated, containing woodland, shrubland, herbaceous vegetation, and cultivated land, or predominantly nonvegetated, containing barren land, urban development, as well as farm outbuildings and painted or treated surfaces [8].Importantly, they embody the everyday visual environment, unlike some arrays and tableaux used in laboratory studies.
In addition to enabling visual function to be realistically analyzed and computationally modeled in the natural world, hyperspectral imaging allows more accurate testing of colorimetric and color appearance formulas and more faithful rendering of scene data on color display devices, for both laboratory research and its practical applications.
Unavoidably, the technical literature spans multiple disciplines.The aim, therefore, of this tutorial is to provide an integrated introduction to terrestrial and close-range hyperspectral imaging and some of its uses in human color vision research.
Section 2 describes the main types of hyperspectral cameras and procedures for image acquisition, postprocessing, and calibration.Section 3 compares spectral radiance and reflectance images and the effects of real and simulated illumination changes.Sections 4 and 5 present image transformations for colorimetric representations, color rendering, and cone receptor and postreceptor coding.Sections 6 and 7 contain example applications.These include the calculation of color gamut, color constancy, and metamerism, and an analysis of the utility of color in tasks requiring an observer to distinguish surfaces or identify them under illuminant changes.Notation is guided by [9][10][11][12][13] but amended where necessary for clarity or consistency with previous work.As with other tutorials [14], no attempt is made at a comprehensive review of the field or to report novel research findings.Citations may not reflect priority.
The color images and associated data in this article are from hyperspectral images acquired with an in-house wavelengthscanning hyperspectral camera, which provided images of size 1344 × 1024 pixels over a spectral range 400-720 nm in 10 nm steps [8].Apart from material in Section 2, however, where different camera types are described, none of this tutorial is specific to any one imaging system.

HYPERSPECTRAL IMAGING
A hyperspectral image can be treated as a three-dimensional set of data, sometimes called a datacube or hypercube, with values representing a spectroradiometric quantity, such as spectral radiance or spectral reflectance, indexed by spatial coordinates u, v and wavelength λ.A two-dimensional section of the datacube, known as a slice, can be defined for a particular passband center wavelength λ, or for a particular spatial coordinate u or v, to within the sensor resolution of the camera.Figure 2 shows an array of wavelength slices.Each is a grayscale intensity image at wavelengths of 400 nm, 410 nm, …, 720 nm making up the hyperspectral radiance image of the flower scene in Fig. 1.The different spectral reflecting properties of the surfaces are evident, with, e.g., the flowers on the left of the scene appearing brighter in the grayscale images from around 600 nm onward.
Depending on the hyperspectral camera, images can be acquired a slice at a time or in a single operation.Significant postprocessing of the raw signal is required either internally by embedded hardware or externally by the user in order to reduce the effects of noise and to obtain accurate radiance or reflectance estimates.

A. Hyperspectral Cameras
In the present context, which hyperspectral cameras are the most suitable for scene imaging?Systems based on wavelength scanning, line scanning, and non-scanning, single-shot imaging are all popular [16,17,42].Non-scanning systems offer greater efficiencies than scanning systems [43] but may entail compromises in resolution.
Active hyperspectral imaging, such as terrestrial laser scanning, delivers only reflectance data [44].

Wavelength Scanning
In wavelength-scanning imaging, also known as area or focalplane scanning, or front staring imaging, the incident light passes through a tunable spectral filter placed in front of the camera lens or between the camera lens and camera sensor.
Spectral sampling may be achieved with a fixed number of mechanically interchangeable bandpass filters [45][46][47], or with an electronically tunable liquid crystal or acousto-optical filter Examples of reflected radiance spectra from a flower scene.Data are from a hyperspectral radiance image of size 1344 × 1024 pixels, corresponding to approximately 14 × 11 deg visual angle at the camera, with spectra sampled at 400 nm, 410 nm, …, 720 nm.The plots show radiance spectra at individual pixels (radiance scales adjusted for range).The small light-gray sphere near the top of the scene is covered in Munsell N7 matte paint [4], and the reflected spectrum (top right plot) follows the typically uneven spectrum of light from the sun and sky [5].The long, thin, light-gray rectangular plate at the bottom of the scene is a reference reflectance surface.The hyperspectral image of the scene was acquired on October 10, 2003, from a garden in Sameiro in the Minho region of Portugal.
or an interferometer [17,42,48,49].The image is constructed either directly, one wavelength at a time, or indirectly from the interferogram.
Wavelength-scanning systems have the advantage that with well-corrected optics, they preserve the geometry of the scene or object during spectral sampling, and exposure time can be adjusted with wavelength according to the spectral properties of the scene.They have the disadvantage that acquisition time can be prolonged at low light levels and that spatial resolution is limited by the size of the sensor array.Additionally, scene movement may generate chromatic artifacts in the rendered image.

Line Scanning
With a line-scanning or pushbroom hyperspectral camera, the incident light is focused by a lens onto a slit aperture from which the light is spectrally dispersed by a grating or other device [50] and then focused onto the camera sensor.The image of the scene is thus constructed a row (or equivalently a column) at a time, and at every point along the row, a complete spectrum, indexed by wavelength λ, is recorded [51].Successive rows, indexed by spatial coordinate v, may be acquired with a moving mirror or by rotating the whole camera [16].
Line-scanning systems have the advantage that spatial sampling in the other spatial dimension u is, in principle, unlimited, but the disadvantage that asymmetric sampling may lead to distortion [52].Exposure time cannot be adjusted with wavelength in the same way as with wavelength-scanning systems, and large variations in incident radiance with wavelength can prove a challenge [53].Line-scanning systems are also vulnerable to scene movement.

Single Shot
With a single-shot or snapshot hyperspectral camera, the incident light is divided simultaneously over multiple spectrally selective sensor arrays, or it is optically coded and spectrally dispersed over a single sensor array, or imaged on a single spectrally patterned sensor array, analogous to the Bayer filter mosaic in a conventional color camera [2,54,55].
Single-shot systems have the advantage that complete spatial and spectral data are acquired simultaneously, but the disadvantage that there is a tradeoff between the number of spatial and spectral samples according to how the incident light is divided or coded [2,43,54,56,57].As the number of spectral samples increases, the number of sensor elements allocated to a particular wavelength generally decreases.The alternative of increasing the size of the sensor array can lead to longer integration or readout times.Other relevant issues include spatial and spectral aliasing [2,58], demosaicing artifacts [56,59], and, with spectrally patterned arrays, the consistency of wavelength tuning across corresponding sensor elements.

B. Image Acquisition
The workflow for image acquisition outdoors depends on the choice of hyperspectral camera, viewing geometry, and whether spectral radiance or reflectance data are required.The following steps are typical for a spectroradiometrically uncalibrated scanning system [8,16,25,32,50,60,61]: 1. Adjust the camera position and focus for the scene (or region or object).
2. To aid later spectral calibration (Sections 2.C.5, 2.C.6), introduce into the scene or field of view one or more small, spectrally neutral, reference surfaces, e.g., a barium sulfate plug or a Munsell gray sphere or plate [4], as in Fig. 1, or a multicolor chart (e.g., ColorChecker, X-Rite Pantone, Grand Rapids MI, USA) for additional control measurements [62,63].Ensure that the reference surface is oriented correctly with respect to the camera axis and incident illumination.Insert or identify any structures in the scene that can provide fiducial marks for later spatial calibration and image registration in postprocessing (Section 2.C.4).
3. Adjust the aperture of the camera, and, if necessary, estimate the exposure timing in a trial acquisition sequence, avoiding point specularities.
4. Acquire a hyperspectral image of the scene.5.If the image is to be processed as a spectral radiance image rather than a reflectance image, record the spectral radiance reflected from the reference surface or surfaces in the scene or field of view with a calibrated spectroradiometer or telespectroradiometer 6. Acquire a hyperspectral image of a uniform, neutral, matte, flat, reflecting plate placed in front of the camera to cover the field of view.This flat-field image is used to estimate fixed-pattern noise (Section 2.C.1), although its effects are more systematic than random.As a spatial measurement, it is distinct from the spectral measurements with the reference surface in steps 2 and 5. Flat-field images can also be obtained in the laboratory from a full spectrum source with an isotropic diffuser placed in front of the camera or with an integrating sphere [32,64].  1, with wavelength indicated in nm at the top left of each slice.The intensity range in each slice image has been stretched for illustration.As evident from the spectral radiance plots in Fig. 1, the surfaces in the scene reflect little energy at very short wavelengths.7. Acquire a hyperspectral image with the lens covered but with the same exposure timing as for the scene image and at the same temperature.This dark-field image is used to estimate the combined effects of dark-current noise, bias, and readout noise.
8. If necessary, repeat the scene image acquisition to reduce uncorrelated image noise by subsequent averaging.Photon noise can be a major source of error.Both flat-field and dark-field image acquisitions should also be repeated several times [32].The plate used for flat-fielding may be moved a little between acquisitions to ensure the uniformity of the estimate.
Image sequences should be examined for defects, most often movement artifacts and sensor element saturation resulting from transient point specularities.Defective sequences may need to be discarded.If the imaging conditions are constant, then common flat-field images and common dark-field images may be used.Details of the scene and imaging conditions should be recorded as part of the acquisition documentation.

C. Image Postprocessing
After acquisition, images obtained with an uncalibrated system need to be corrected for noise errors and other systematic errors and then calibrated spatially to give angular subtense and spectrally to give spectral radiance or reflectance [8].Correcting for achromatic geometric errors and spatial blurring are considered in Refs.[65,66] and for spectral blurring in Refs.[67,68].Notation here differs somewhat from that in Refs.[10,13].

Noise Errors
Errors due to noise can be divided into non-imaging and imaging errors.Non-imaging noise errors include camera dark-current noise, which increases with exposure duration and temperature; bias, which provides a constant offset; and readout noise.The main imaging noise error is fixed-pattern noise, which includes non-uniformity in the sensor array response, i.e., pixel-to-pixel variation in quantum efficiency, gain, and dust [10,31,69], and off-axis vignetting, i.e., darkening at the corners and edges of the image [27,64].These and other sources of noise are analyzed in Refs.[28,31,32].
Let I S u, v; λ be the scene image intensity indexed by spatial coordinates u, v, and wavelength λ; let I F u, v; λ be the corresponding flat-field image intensity; and let I D u, v; λ be the corresponding dark-field image intensity (independent estimates of dark-current noise, bias, and readout noise are usually unneeded [32]).The approximately corrected scene image intensity I CS u, v; λ is then given [32,70] by This expression assumes that after subtraction of dark-field image noise I D u, v; λ, the response of each sensor element is linear.Both CCD and sCMOS sensors can deliver an approximately linear response to light over several decades [32,69], and any departures from linearity can be accommodated by adding correction terms to Eq. (1).Fluctuations in the recorded values of I D u, v; λ may lead to the corrected scene image intensity I CS u, v; λ in Eq. ( 1) taking on numerically small, negative values at some points.These values may be set either to zero or to some mean estimate derived from other areas of the image.
A correction for stray light or light scatter within the camera may also be applied [29,71,72], but in practice, it is necessary only if there are locally very bright reflecting or emissive surfaces within the scene.
Notice that values of I CS u, v; λ do not represent spectroradiometric quantities until the image is calibrated for spectral radiance or reflectance (Sections 2.C.5, 2.C.6).

Wavelength Errors
Errors in the nominal passband center wavelength setting of the system can be quantified with a calibrated precision monochromator.Errors in the passband center wavelength across the field of view can be quantified by repeating a flat-field measurement in which the illumination is from spectral line sources [64].Both kinds of wavelength error ought to be small and capable of being accommodated by corrections to Eq. (1).

Spectral Sampling and Out-of-Band Transmission
The wavelength sampling interval of the imaging system is usually larger than the passband width, so sampling can be approximated by a delta function [28].Alternatively, a routine correction can be applied for bandpass error [13,73].
Spectral leakage, i.e., transmission at wavelengths outside the selected passband of the system [10,25,64], may need to be considered.It may be extremely small, several orders of magnitude less than maximum transmission [74], yet, in principle, its effects can be cumulative at long wavelengths where the radiant energy of the sun is greater, vegetated surfaces reflect more strongly (the "red edge"), and the quantum efficiency of the camera sensor is higher.A correction for spectral leakage can be made as follows.
Let tλ 0 , λ be the spectral transmission of the system as a function of wavelength λ at each passband center wavelength λ 0 .The in-band portion of this function can be measured with a calibrated precision monochromator, but the residual out-ofband transmission may need other methods.Let I u, v; λ and Ĩ u, v; λ be the respective image spectral intensities at the point u, v before and after spectral sampling.Then, The tilde symbol signifies that Ĩu, v; λ is an approximation to I u, v; λ.The discretized form of Eq. ( 2) can be solved by matrix methods to yield the required estimate of I u, v; λ, with due attention paid to the accuracy of the numerical data for t and the stability of the solution.
In practice, with an infrared-blocking filter in place, the correction is usually very small and may be omitted under good imaging conditions.Nevertheless, one or more calibrated surfaces that reflect predominantly at long wavelengths, e.g., Munsell 10YR [4], may be included in the scene to monitor the extent of long-wavelength leakage.

Image Registration
Longitudinal chromatic aberration of the optical system and consequent chromatic image blur are reduced with achromatic or apochromatic lenses [75].Even so, lateral chromatic aberration, i.e., variation in optical image size with wavelength, may still be visible at the extremes of the spectrum, resulting in misalignments of several pixels across the full image width.Figure 3, top left, shows a color image of a stone cottage rendered from a hyperspectral radiance image with an enlarged copy of the square area marked in white, bottom left.The color fringing on the drainpipe is obvious, but also detectable on the branches of trees and the walls of buildings.
The principle of image registration is to take one wavelength slice as a reference or base image, usually from the middle of the spectrum, and apply operations such as scaling and translation to the other wavelength slices to optimize the alignment, quantified, e.g., by an areal measure such as normalized cross correlation [76].Figure 3, top right, shows the effect of registration.An enlarged copy of the square area marked in white is shown bottom right.The color fringing is eliminated.
As an alternative to areal registration, image features such as edges or corners may be identified and local edge fitting applied to define control points for a single scaling and translation operation [76].Non-Euclidean transformations may also be applied [77].
The method of registration needs to be chosen with reference to the composition of the scene, since regions of the image may reverse contrast from one part of the spectrum to another.The effect is evident in Fig. 2. At 520 nm, the flower on the lower right of the scene appears darker than the leaf beneath it in the grayscale image, but at 550 nm, it appears brighter, with the result that the correlation between the slices changes sign.In these parts of the spectrum, regions, edges, and corners may lose their identity or shift spatially.
The problem of aligning wavelength slices differing significantly in structure can be avoided by applying registration iteratively, each slice being compared with the next [78], though with the risk that errors may accumulate across the spectrum.
Whichever method is chosen, the edge of the registered image may need to be trimmed to eliminate pixel artifacts with zero or replicated padded values.

Spectral Radiance Calibration
A spectral radiance image can be obtained from the corrected scene image by exploiting flat-field data (Section 2.B), but the calculation requires additional measurements and contains uncertainties.A more robust way to proceed is to use the spectral radiance data recorded independently from the reference surface or surfaces in the scene or field of view (Section 2.B).
Suppose that l 0 λ is the spectral radiance recorded from a reference surface at a point u 0 , v 0 as a function of wavelength λ.Then, for the selected viewing geometry, the spectral radiance Lu, v; λ at the point u, v is estimated by This estimate, signified by the hat symbol, is called a one-point calibration [26].It can be improved by using data from several reference surfaces [70].Suppose that l 0 λ, l 1 λ, …, l n λ are the spectral radiances recorded from these surfaces at the n 1 points u 0 , v 0 , u 1 , v 1 , …, u n , v n .Perform an mthorder regression of l 0 λ, l 1 λ, …, l n λ on I CS u 0 , v 0 ; λ, I CS u 1 , v 1 ; λ, …, I CS u n , v n ; λ, where 1 ≤ m ≤ n.Thus, if m 1 and the estimated regression coefficients are k0 λ and k1 λ, then Lu, v; λ k0 λ k1 λI CS u, v; λ: Notice the dependence of the coefficients k0 λ and k1 λ on λ; i.e., the regression is performed at each wavelength.The goodness of the fit can be tested by cross-validation against known reflecting surfaces in the scene.The offset term k0 λ ought to be unnecessary if the estimate of the dark-field image intensity I D u, v; λ in Eq. ( 1) is correct.

Spectral Reflectance Calibration
The procedure for obtaining a spectral reflectance image from the corrected scene image I CS u, v; λ is analogous to that for obtaining a spectral radiance image.But only the spectral reflectance r 0 λ at a point u 0 , v 0 on the reference surface is required as a function of wavelength λ, not its spectral radiance l 0 λ.For the selected viewing geometry, the spectral reflectance Ru, v; λ at the point u, v is estimated by This scaling of I CS u, v; λ by r 0 λ∕I CS u 0 , v 0 ; λ, unlike scaling for spectral radiance in Eq. ( 3), may amplify image errors at short wavelengths, when the incident illumination is weak.But as with spectral radiance estimates, this one-point calibration can be improved by using data from several reference surfaces and performing an mth-order regression of spectral reflectances r 0 λ, r 1 λ, …, r n λ on If m 1 and the estimated regression coefficients are k0 λ and k1 λ, then Ru, v; λ k0 λ k1 λI CS u, v; λ: As with the spectral radiance calibration, the offset term k0 λ ought to be unnecessary if the estimate of the dark-field image intensity I D u, v; λ in Eq. ( 1) is correct.

Spatial Subsampling
Even after registration, data from adjacent pixels in the hyperspectral image of a scene are likely to be correlated because of the limited resolution of the camera optical system.The pointspread function (PSF), which in general varies with wavelength, provides a useful guide.It can be estimated by methods using point imaging [79], edge imaging [80], and random-noise imaging [81].If the PSF is much larger than the pixel spacing, the spectral radiance image Lu, v; λ or spectral reflectance image Ru, v; λ may need to be spatially subsampled before being used in analysis and modeling.

RADIANCE AND REFLECTANCE ESTIMATES
Given the procedural similarities in calibrating for spectral radiance and spectral reflectance, which of the two representations is the more appropriate?It depends on how the data are used, the differences between real and simulated illumination changes, and the variations in natural reflected spectra.

A. Effective Reflectances and Global Illuminants
Recall from Section 2.C.6 that a spectral reflectance image Ru, v; λ, indexed by spatial coordinates u, v and wavelength λ, is obtained by scaling the corrected scene image I CS u, v; λ by the spectral reflectance of one or more reference surfaces in the scene.It is strictly an effective reflectance in the sense that it represents the scene for the selected viewing geometry as if it consists of planar Lambertian reflecting surfaces illuminated by a spatially uniform or global illuminant Eλ, which need not coincide with the true illumination on the scene [8].More precisely, suppose, again, that l 0 λ is the spectral radiance recorded at a point u 0 , v 0 on a reference surface in the scene or field of view (Section 2.C.5) and that r 0 λ is the known spectral reflectance of that surface (Section 2.C.6).Then, the global illuminant Eλ is defined by Eλ l 0 λ∕r 0 λ: The reflected spectral radiance Lu, v; λ may then be represented as the product of Eλ and the effective spectral reflectance Ru, v; λ; thus, It is immaterial whether the reflected spectral radiance is estimated directly by Eq. ( 3) or indirectly through Eq. ( 4), since the outcome is the same.Effective spectra are sometimes referred to as apparent [82] or equivalent [83,84].Any nonlinear contributions to Eq. ( 4) from processes such as mutual reflection and fluorescence are assumed to be negligible.Mutual reflection, where light reflected from one surface illuminates a second [85], also referred to as mutual illumination or interreflection [86], depends on the proximity and orientation of the surfaces in the scene [87,88].Fluorescence, which is commonly associated with vegetation [89], accounts for perhaps 4%-15% of the nominal reflectance in the 670-690 nm region of the spectrum [90].Both contributions, which are difficult to predict for the general case, are probably outweighed by other factors, as illustrated in Section 3.B.
These specific nonlinearities aside, the reflected radiance Lu, v; λ is the result of a fundamentally more complicated interaction than that summarized by Eq. ( 4), since it involves the bidirectional reflectance distribution function (BRDF) [9,91].With an abuse of notation, let Eθ, φ; u, v; λ be the incident spectral irradiance at the point u, v in the direction θ, φ at wavelength λ, where the polar and azimuthal angles θ and φ are defined with respect to a fixed directional coordinate system.Likewise, let Rθ 0 , φ 0 ; θ, φ; u, v; λ be the bidirectional spectral reflectance distribution function at u, v in the direction θ, φ at wavelength λ, for the viewing direction θ 0 , φ 0 , where the points of incidence and reflection are assumed to be the same [9].The reflected spectral radiance Lu, v; λ in the direction θ 0 , φ 0 is then given by where dω is the solid angle element in the direction θ, φ, and the explicit dependence on θ 0 , φ 0 is suppressed.
The role of direction of incidence θ, φ and position u, v in Eθ, φ; u, v; λ and Rθ, φ; u, v; λ is nontrivial.In the product of effective quantities in Eq. ( 4), the variation in the incident spectral irradiance is absorbed into the effective spectral reflectance, whose magnitude may exceed unity at some points because the reference surface is oriented with its normal to the camera and at an angle to the sun.These unreal values of Ru, v; λ lose their significance when multiplied by Eλ in Eq. ( 4) to form the reflected spectral radiance Lu, v; λ.
It is emphasized that in practice the effective spectral reflectance Ru, v; λ recorded from a physical matte object whose spectral reflectance is nominally constant across its surface can still vary markedly from pixel to pixel.In addition to photon and acquisition noise identified in Sections 2.B and 2.C, there are multiple other sources of noise associated with physical surfaces under natural illumination.These sources are enumerated in Section 3.C.

B. Real and Simulated Illumination Changes
Whether a reflected spectral radiance can be represented as a product of a global illuminant and an effective spectral reflectance, as in Eq. ( 4), depends on the separability of the integrand in Eq. ( 5), i.e., with a further abuse of notation, whether Eθ, φ; u, v; λ can be approximated as the product of spatial and spectral factors: In that event, the expression in Eq. ( 5) for the reflected spectral radiance Lu, v; λ becomes where the explicit dependence of Lu, v; λ and Rθ, φ; u, v; λ on θ 0 , φ 0 is again suppressed.This equation reduces to Eq. ( 4).The approximation can be usefully tested with two kinds of natural illumination change, one mainly spectral, the other mainly geometric.The color of the illumination is usually identified with its correlated color temperature in kelvin [13,92], where 4000 K typifies the yellow-orange of the setting sun and 25,000 K the blue of the north or polar sky.For the scene in Fig. 4, the spectrum of the direct illumination on the scene can be calculated from the hyperspectral image radiance at the small neutral matte Munsell N7 rectangular plate at the bottom right of the scene, arrowed in (a).The correlated color temperature of the incident light changes from 4180 K at 18:15 in (a) to 3350 K at 18:40 in (b).For later reference, notional average daylight D65 has a correlated color temperature of approximately 6500 K [22], which on a reciprocal color temperature scale [92], falls roughly midway between 4000 K and 25,000 K. Can the reflected spectral radiances in Figs.4(a) and 4(b) be adequately represented as products, as in Eq. ( 4), with a common effective spectral reflectance?Suppose that E 1 λ and E 2 λ are global illuminant spectra with correlated color temperatures 4180 K and 3350 K, respectively, and that L 1 u, v; λ and L 2 u, v; λ are the corresponding spectral radiances reflected from the scene.Suppose further that R 1 u, v; λ is an effective spectral reflectance defined by Figure 4(c) shows the color image rendered from L2 u, v; λ according to Eq. ( 6).This simulated image and the real one in Fig. 4(b) are closely similar, though there are small differences in the shadows.The normalized root mean square deviation between the spectral radiances is 1.9%, where normalization is with respect to the range.

Geometric Illumination Changes
More frequently, it is geometric effects that dominate natural illumination changes, as illustrated in Figs.
Figure 4(f ) shows the color image rendered from the approximation L2 u, v; λ for L 2 u, v; λ according to Eq. ( 6).The simulated image in (f ) and the real one in (e) are now clearly different.The normalized root mean square deviation between the spectral radiances is 13.7%, about seven times larger than in (b) and (c).
Although the assumption of a common effective spectral reflectance R 1 u, v; λ for the two images is here unsafe, some local statistical properties may nonetheless be preserved [93].
Methods for segmenting the effects of multiple illuminants are described in Refs.[94,95].

C. Variations in Natural Spectra
Reflected spectra recorded outdoors may show much greater variation than in the laboratory for several reasons: 1. Inorganic and organic materials have significant microstructure, which combined with dirt, weathering, biological degradation, and moisture can give non-Lambertian behavior [3], including specularities [96,97].
3. Natural illumination is inconstant.In a clear sky with no visible cloud, the primary source of variation is solar elevation, but unobservable cirrus clouds and aerosols can be highly inhomogeneous, leading to additional variation in the solar beam.Average minute-to-minute temporal fluctuations in irradiance can be on the order of 0.1% around midday and much more around sunrise and sunset (data provided by A. R. D. Smedley and A. R. Webb, University of Manchester).The angular and spectral distributions of light from the sun and sky are also very different [5], and the balance between the direct beam and diffuse radiation and the radiance distribution across the sky change with solar elevation.Nearer to the scene being imaged, object occlusions, mutual reflections between surfaces, and transilluminance all contribute to further variation in the incident light [103].
4. Reflected spectra from adjacent surfaces in scenes suffer local mixing at the camera because of the resolution of the optical system and sensor array.Mixing adds to the variation in the estimated surface spectral reflectances [104,105].
On the other hand, in the laboratory, reflection spectra are typically recorded from relatively homogenous materials with a non-imaging device such as a reflection spectrophotometer, which averages over extended areas under controlled illumination.Where materials and imaging conditions coincide, however, closely similar spectral estimates ought to be obtainable from a correctly calibrated hyperspectral imaging system and a non-imaging spectroradiometer [8, Fig. 1].

COLORIMETRIC REPRESENTATION
For colorimetric calculations, psychophysical measurements, image display, or other related purposes, hyperspectral radiance images may be mapped into a variety of color spaces.How, then, to decide on which space? Factors to consider include the effect of different scene illuminants in the space, the uniformity of the space, i.e., whether equal distances correspond to equal color differences, and adjustments for the characteristics of any display device.Color spaces modeled on physiological data are discussed in Section 5.
In colorimetric applications, it is conventional to generate an intermediate representation expressed in CIE XYZ tristimulus values [21,22].Thus, let xλ, ȳλ, zλ be the CIE XYZ [13] color-matching functions for the 2°standard observer, and let Lu, v; λ be a spectral radiance image, indexed by spatial coordinates u, v and wavelength λ.The tristimulus values X , Y , Z at each point or pixel u, v are given by where the constant k is chosen so that Y 100 for a perfectly white surface under full illumination.Because of this normalization, the CIE XYZ representation of the image requires only relative not absolute spectral radiance values [21].
There is generally no easy relationship between the RGB values of a scene recorded by a color camera and the corresponding CIE XYZ values.The camera spectral sensitivities would need to be linear transformations of the CIE colormatching functions, i.e., meet the Luther condition [21,106].
To simplify notation, the explicit dependence of X , Y , Z on position (u, v) is suppressed.

A. Compensating for Different Illuminations
Given a trichromatic representation of a scene, a color image, it is possible to compensate approximately for a pure spectral change in illumination by transforming the representation.The transformation, traditionally called a chromatic adaptation transform, transforms the colors of the scene under one illuminant, referred to as the test, to the corresponding colors under another illuminant, the reference, so that, ideally, they have the same color appearance [21,22,107].
Historically, chromatic adaptation transforms are based on von Kries scaling [22,92].The idea is to compensate for a shift in the color of the illumination, e.g., towards long wavelengths, by reducing the gain of the appropriate cone receptor class, here the long-wavelength-sensitive cones.This adjustment is independent of any adjustments to the medium-and short-wavelength-sensitive cones.In matrix notation, it is a diagonal matrix transformation [108].
In practice, with color images represented by their CIE XYZ tristimulus values, additional transformations are needed to convert the tristimulus values to cone responses or some sharpened version [109] for von Kries scaling and then back again.Such a transformation is incorporated into the standard chromatic adaptation transform CMCCAT2000, which takes as parameters the tristimulus values of the test and reference illuminants X t , Y t , Z t and X r , Y r , Z r , along with some other optional parameters [110].Routines for performing chromatic adaptation transforms are available from several sources, e.g., Ref. [21].
The reason that chromatic adaptation can compensate only approximately for illuminant changes is that for most scenes, the number of degrees of freedom for surface reflectance and illumination spectra is greater than for their trichromatic representations.Chromatic adaptation transforms can, though, be optimized for particular sets of reflecting surfaces.These are often manufactured sets, such as textiles and Munsell papers [110].But unlike hyperspectral reflectance data, they do not capture the distributions of spectral reflectances in the natural world.

B. Representing Hyperspectral Images in Uniform Color Spaces
Although CIE XYZ tristimulus values allow the equality of image colors to be decided under the same conditions, equal differences in tristimulus values need not correspond to equal differences in discriminability or color appearance.What is required is a uniform color space, i.e., one in which equally spaced points correspond to equal perceptual color differences [92,111].Standard examples include the approximately uniform color space CIELAB [13] and the more uniform space CIECAM02 [12], used here as a color space, not an appearance model.The coordinates of CIELAB space are L , a , b , where L represents lightness, ranging from 0 to 100, and a and b represent nominally red-green and yellow-blue chromatic components, respectively.The corresponding coordinates of CIECAM02 space are J, a C , b C [21,22].
It is straightforward to map a hyperspectral radiance image into CIELAB space.Let X , Y , Z be the tristimulus values at each point u, v and X t , Y t , Z t be the tristimulus values of the scene illuminant at the reference surface, the illuminant being assumed global.The transformation of X , Y , Z to CIELAB values L , a , b with respect to X t , Y t , Z t is defined [13] as follows: where f is a compressive function of its argument r, i.e., 841∕108r 4∕29 otherwise: The color difference ΔE across two points u 1 , v 1 and u 2 , v 2 is defined by the Euclidean metric; i.e., if ΔL , Δa , Δb are the respective individual differences in the components . The chromatic projections a , b may be converted into the more intuitive polar form, with chroma C ab a 2 b 2 1∕2 and hue angle h ab tan −1 b ∕a , where 0°corresponds to the red end of the red-green axis and 90°to the yellow end of the yellow-blue axis.
As an illustration, the CIELAB hue angles of the selected samples from the flower scene in Fig. 1 are 140°and 114°f or the green leaves, left and top; 49°for the orange flower, left; 305°(i.e., −55°) for the purple flower; and 71°for the yellow flower, right.Further discussion can be found in Refs.[21,22].
Somewhat inconveniently, CIELAB space is not suited to comparing differences ΔE between surfaces under illuminants other than average daylight [13,22].The reason is that in Eq. ( 8) the illuminant appears in a quotient of tristimulus values, whereas to accurately accommodate chromatic adaptation, the quotient should be of cone responses (Section 4.A).The formula delivers what is sometimes called the wrong von Kries transformation [108].
A better procedure [13,22] for comparing color differences under changes in illuminant is to first apply an accurate chromatic adaptation transform with respect to some standard reference illuminant (not to be confused with the actual illuminant at the reference surface in the scene) and then transform the result to CIELAB coordinates.The chromatic adaptation transform CMCCAT2000 is used in this analysis, but recall from Section 4.A that chromatic adaptation can compensate only approximately for illuminant changes.
As before, let X , Y , Z be the tristimulus values at each point u, v and X t , Y t , Z t be the tristimulus values of the scene illuminant at the reference surface.Then: 1. Let X r , Y r , Z r be the tristimulus values of the reference illuminant, e.g., an equi-energy illuminant or standard daylight with correlated color temperature of 6500 K.
2. Apply a chromatic adaptation transform such as CMCCAT2000 with parameters X t , Y t , Z t and X r , Y r , Z r to transform X , Y , Z to corresponding values X , Ỹ , Z under the reference illuminant. 3.
The procedure is the same for any other scene illuminant.The separate issue of the approximate uniformity of CIELAB space can be accommodated by replacing the Euclidean metric with a more accurate color-difference formula, e.g., CIEDE2000 [13,112,113].
For CIECAM02 space, color differences ΔE are defined by the Euclidean metric, i.e., ΔE ΔJ 2 Δa C 2 Δb C 2 1∕2 , or some other orthogonal combination [114].But the comparison of color differences across scene illuminants is simpler than with CIELAB space because an accurate chromatic adaptation transform is built into CIECAM02 space.Its uniformity is also greater than that of CIELAB space, although it too can be improved [115].

C. Color Rendering Hyperspectral Images
To present a spectral radiance image on a display device requires a color space for rendering and some knowledge of the display device [116].Color rendering for image illustration can be treated as a form of spectral visualization [117].
With the default RGB color space sRGB [118], rendering proceeds as follows.First, radiance values are converted to CIE XYZ tristimulus values X , Y , Z according to Eq. ( 7), but rescaled to range from 0 to 1 rather than the default from 0 to 100.Finally, to satisfy range constraints, sRGB values less than 0 are set to 0 and sRGB values greater than 1 set to 1.
A nonlinear correction can be applied to compensate approximately for the input-output function of the display device, sometimes referred to as gamma encoding [116], despite its ambiguity [118].A typical approximate correction has the form More accurate and comprehensive corrections, including an offset for black, are described in Refs.[21,116,119].
A key distinction between these transformations and those in Section 4.B is that color rendering transformations depend on the characteristics of the display device.
Figure 5, left, shows an sRGB rendering of a hyperspectral radiance image of a yellow flower (adapted from [120]).The image appears dark because the exposure timing used during image acquisition was limited by the specular highlight from the shiny leaf at the top right of the scene.To improve appearance, the sRGB levels of the image can be clipped to the level of a less specular but still bright region, such as the arrowed area on the light-gray sphere in the left image [121,122], and then scaled.Thus if c is the maximum of the sRGB values in the arrowed area, the clipped and scaled sRGB values at each point are given by R c minfR, cg∕c, G c minfG, cg∕c, B c minfB, cg∕c: The result is shown in Fig. 5, right.The detail in the darker areas is now much clearer.Clipping affects about 1% of the pixels in this scene.Clipping can, however, lead to artifacts, and more nuanced approaches exploit gamut-mapping algorithms to preserve the relationships between colors [122] and tone-mapping algorithms to preserve the detail in both light and dark regions [116].

RECEPTOR AND POSTRECEPTOR CODING
Estimating cone receptor and postreceptor responses to a spectral radiance image is more complicated than mapping points into a color space.The receptor mosaic is spatially non-uniform; preceptor absorption in the ocular media varies between individuals and with age; and chromatic adaptation and postreceptor coding take place at multiple levels.Calculations may, though, be simplified with the use of standardized observer data and viewing conditions [13].
As with colorimetric representations, there is generally no easy relationship between the RGB image of a scene recorded by a color camera and the corresponding receptor or postreceptor responses, standardized or not.
To simplify notation, the explicit dependence of cone signals on position (u, v) is suppressed.

A. Receptor Responses
For a 2°standard observer [13], the calculation of normalized cone excitations is exactly analogous to the calculation of CIE XYZ tristimulus values (Section 4), since both are transforms of the color-matching functions.
Let Lu, v; λ be a spectral radiance image, indexed by spatial coordinates u, v and wavelength λ.If S L λ, S M λ, S S λ are the long-, medium-, and short-wavelength-sensitive cone spectral sensitivities, measured at the cornea, i.e., incorporating preceptor absorption, then at each point u, v, the corresponding cone excitations q L , q M , q S are given by The spectral sensitivities S L λ, S M λ, S S λ are usually normalized to a maximum of unity on a linear energy (or quantal) scale [123].Cone excitations may instead be calculated by linear transformation of the corresponding XYZ tristimulus values [19,21].Individual differences in observer preceptor effects may be accommodated by replacing cone spectral sensitivities measured at the cornea by pigment spectral sensitivities combined with individual lens absorption and macular pigment data, as described in Refs.[124,125].
Because of the variation in receptor density and receptor type across the retina, assumptions may need to be made about how the gaze of the observer samples the spectral radiance image.For vision outside the central fovea, the intrusion of rod signals [126,127] and melanopsin-generated signals from intrinsically photosensitive retinal ganglion cells [128,129] may also need to be considered.

B. Postreceptor Transformations
The cone excitations estimated at each point u, v may be transformed in several ways, according to the requirements of the simulation or analysis and whether cone-opponency or color-opponency is relevant [130].
Given the large overlap in long-and medium-wavelengthsensitive cone spectral sensitivities, a frequent computational objective is to reduce redundancy.Other objectives may be to optimize chromatic adaptation or to maximize the information that can be extracted from scenes.The postreceptor transformations are often assumed to be linear, albeit sometimes preceded by a nonlinear stage [131].
An early study by Buchsbaum and Gottschalk [132] provided an influential model [133].The postreceptor transformation is designed to linearly decorrelate cone signals q L , q M , q S in response to random combinations of monochromatic spectra.Cone spectral sensitivities are from the Vos-Walraven fundamentals [134], which differ a little from those due to Stockman and Sharpe [124,125].The result is an achromatic response r A , a red-green chromatic response r RG , and a yellowblue chromatic response r YB ; thus, q L q M q S 3 7

5:
The terms "red-green" and "yellow-blue" are used generically in this context.The number of significant figures in the matrix entries follows author usage.In Ref. [109], the postreceptor transformation maximally constrains responses to narrow bands of wavelengths.These sharpened spectral sensitivities improve the action of von Kries scaling in compensating for the effects of illuminant changes with Munsell reflectance spectra.The cone fundamentals are derived from [134][135][136].The postreceptor responses r # L , r # M , r # S are given by 2 q L q M q S 3 7 5: In Ref. [137], the postreceptor transformation maximizes the Shannon information available from natural scenes under simulated daylight changes.Responses are assumed to be subjected to von Kries scaling and values are therefore normalized to a mean of 1.0.The cone fundamentals are from [124,125].The postreceptor responses r # L , r # M , r # S are given by 2 q L q M q S 3 7

5:
Unsurprisingly, these postreceptor transformations are only loosely correlated with neural function [138,139].Still, with variable coefficients, they enable performance to be optimized with respect to specific functional objectives.By contrast, simpler, fixed, postreceptor transformations may be used to explore possible psychophysical implications of the data, e.g., by mapping cone excitations into a particular color space such as the MacLeod-Boynton chromaticity diagram [140] or the color space due to Derrington, Krauskopf, and Lennie [141], as illustrated in Refs.[18,138].The same caveat about neural correlates applies [19].

PROPERTIES OF SCENE COLORS
Some of the most basic color properties of a scene have to do with the existence and stability of colors, e.g., how many distinguishable colors it contains, the constancy of its colors under changes in illumination, and the presence of aliased colors, i.e., metamerism.
Since color is used here as a proxy for spectral reflectance or radiance, there is no reference to the spatial characteristics of the scene, at least explicitly, though procedures intended to compensate for illumination changes may well do so.Critically, these color properties set a least upper bound on what observers can infer about scenes from surface color alone.
Elsewhere, interest may be in color phenomena that do relate to spatial characteristics.For example, an extension of CIELAB known as S-CIELAB [142,143] incorporates preprocessing by the approximate spatial contrast sensitivity functions of achromatic and chromatic visual pathways.Some of the relationships between the chromatic and spatial properties of natural scenes are analyzed in Refs.[144][145][146][147].
The material in this section is based partly on [148].

A. Number of Distinguishable Colors
Because natural spectra come from a restricted population of reflecting surfaces and illuminants, the number of distinguishable colors in any one scene is likely to be very much smaller than the maximum of 1.7-2.3 million associated with an idealized population, the theoretical object color solid [149][150][151][152]. How, then, to estimate the number of distinguishable colors, N say, in a given scene?With an approximately uniform color space such as CIELAB or CIECAM02 and a nominal observer threshold color difference, ΔE thr say, an estimate of N can be obtained from a hyperspectral radiance image, or from a reflectance image in combination with a global illuminant, as follows: 1. Map each point of the scene into the color space.2. Segment the color space into unit cells, most simply cubes of side equal to the threshold interval ΔE thr .
3. Count the number of occupied cells.
Crucially, this estimate of N takes no account of how often each color appears or where it comes from in the scene.All that is required is that a cell in the color space is occupied [35,153].The number of occupied cells defines the number of points that on average are separated by at least ΔE thr .
The estimate of N does, of course, depend on the metric used to define distances in the color space, the size of the discrimination threshold ΔE thr , and the method used to define the unit cell.In traditional color discrimination measurements, representative values for ΔE thr are typically 1.0 for CIELAB [154,155] and 0.5-0.6 with the more accurate color-difference formula CIEDE2000 [156,157].In practical applications, though, representative values may be set several times larger.Unit cells may be cubes, spheres, or dodecahedra [158], and the distribution of cells need not form a regular array [159].More realistically, thresholds may be defined probabilistically [153].
Incidentally, nomenclature in this area varies.Along with colors being said to be distinguishable [152], they may also be called distinct [160], discernible [149,159], and discriminable [148,161].Their description here as distinguishable is for compatibility with more general applications in information theory [162].
Based on 50 close-up and distant hyperspectral radiance images of rural and urban scenes, the average number of distinguishable colors per scene is found to be approximately 2.7 × 10 5 [159], for spherical unit cells and the color-difference formula CIEDE2000 with discrimination threshold ΔE thr 0.6.This average number is almost an order of magnitude smaller than with the theoretical object color solid.Similar calculations can be performed with other kinds of hyperspectral images, e.g., of artworks [35].
Where illumination can be controlled, the number of distinguishable colors can be evaluated under different illuminants, e.g., incandescent, fluorescent, and metal halide lamps, and LED sources [34], which can be optimized for observers with color deficiency [163].

B. Gamut Volume
The range or gamut of the colors in a scene can be quantified by the number of colors distinguishable by an observer.This measure is dimensionless but dependent on the size of the threshold color difference ΔE thr , as explained in Section 6.A.More commonly, however, the color gamut is expressed as a volume in a color space [122].Other issues need then to be addressed such as whether points are treated as members of a discrete or continuous set and whether the boundary of the set should form a convex hull.Several estimation methods are available [122,164].
Values of the gamut volume also depend on the dimensions of the chosen color space, typically CIELAB or CIECAM02 space [165].In CIECAM02 space, gamut volume is measured in cubic CIECAM02 J, a C , b C units [122].
Given an independent estimate V of the gamut volume for some scene, the number of distinguishable colors N may be estimated not by the cell-counting method in Section 6.A, but by dividing V by the volume of the unit cell, in the same color space; i.e., if the unit cell is a cube of side ΔE thr , then A way of interpreting N and V to accommodate the relative frequencies with which colors occur is developed in Section 7.A.

C. Color Constancy
Color constancy refers to the effect whereby the perceived or apparent color of a surface remains constant despite changes in the intensity and spectral composition of the illumination [120,166,167].This definition is not the same as that sometimes used in computational color science, where it refers to the estimation of a representative illuminant on the scene [84,168], i.e., a global illuminant.The difference between real illumination changes and their approximation by global illuminants is discussed in Section 3.B.

Colorimetric Approaches
Consider the scene in Fig. , possibly with a colordifference formula such as CIEDE2000.The color difference ΔE is sometimes called the color inconstancy index [169].
How good is the color constancy of the sample point on the flower if the constancy operation is supplied by the chromatic adaptation transform CMCCAT2000?The CIELAB color difference ΔE under the 6500 K daylight is 2.2, which is larger than the representative CIELAB threshold value of 1.0.With the color-difference formula CIEDE2000, a more accurate estimate of ΔE is 1.6, which is also larger than the corresponding representative threshold value of 0.5-0.6 (Section 6.A).Color constancy with CMCCAT2000 is therefore good under a moderate illuminant change, but not perfect.
More usually, the color difference ΔE is converted into a dimensionless color constancy index CI by scaling by the distance, ΔE 1,2 say, between the two illuminants in the same color space and subtracting from unity [170], i.e., CI 1 − ΔE∕ΔE 1,2 .Perfect constancy corresponds to an index of unity, and the complete absence of constancy corresponds to an index of zero.With the chromatic adaptation transform CMCCAT2000 and color-difference formula CIEDE2000, the color constancy index of the sample point on the flower is 0.92.With CIECAM02 space, it is 0.91.
Although the scene and its contents are not typical of those used in psychophysical measurements, these index values are compatible with the upper limit recorded for observers with a variety of experimental methods, materials, and illuminant changes [120,167].Other formulations for color constancy indices are described in Refs.[120,[171][172][173].
To extend measures of color constancy from individual surfaces or points to the whole scene, global statistics may be deployed such as the mean and median of the distribution of index values.

Preserving Identity
Instead of using colorimetric measures of scene color constancy, there is a more utilitarian way of proceeding, namely, by measuring how far color constancy allows points to be identified by their color under a change in illuminant [137].
Consider, again, the sample point on the flower under the 4000 K daylight, arrowed in yellow in Fig. 6, bottom left.As already explained, under the 6500 K daylight, the same point arrowed in yellow, bottom right, has tristimulus values X 2 , Y 2 , Z 2 , and that with chromatic adaptation provided by CMCCAT2000, the estimated tristimulus values X 2 , Ỹ 2 , Z 2 differ from the true values X 2 , Y 2 , Z 2 .But the color difference might be acceptable if X 2 , Ỹ 2 , Z 2 were closer to X 2 , Y 2 , Z 2 than the tristimulus values of any other points in the scene-in other words, if X 2 , Ỹ 2 , Z 2 constitutes a nearest-neighbor match according to some colordifference formula.In that event, color can still be used to identify the sample point uniquely despite the inaccuracy of the constancy operation.Notice that it is not the size of the color difference that counts but the risk of confusion.
Color constancy can then be quantified by the number of points in the scene that preserve their identity under the change in illuminant, i.e., the number of illuminant-invariant points.For a given constancy operation, this number obviously varies with the scene and the illuminant change.
A way to estimate the number of illuminant-invariant points is developed in Section 7.C.

D. Metamerism
Metamerism refers to the effect whereby two or more lights appear the same but have different spectral compositions [13,92].It can arise in several forms, though the most important in practice is associated with surfaces.Metameric surfaces have different spectral reflectances but match in color under some illuminant, i.e., have the same tristimulus values.Since the match is likely to fail under some other illuminant [92], it follows that surface color is not a reliable guide to material identity.As with several other color phenomena, metamerism is a consequence of the number of degrees of freedom for reflectance and illumination spectra being greater than for their trichromatic representations.
The extent to which spectral reflectances are classified as the same or different is an empirical issue.Here, metameric surfaces are taken to be those for which color differences are subthreshold under one illuminant and suprathreshold by a certain multiple of threshold-the criterion degree of metamerism -under another illuminant [8].
The proportion or relative frequency of metameric surfaces in a scene may be calculated by taking a random sample of pairs of points from a hyperspectral reflectance image.Notice that the sample takes no account of what physically or semantically defines a surface, e.g., the flowers in Fig. 1 or the stone cottage in Fig. 3. Therefore, points from the same object in shadow and in direct illumination are treated distinctly.This situation differs from that described in Section 3.B in that pairs of points are considered as isolated samples undergoing a change in illuminant spectrum [8, Appendix A].That said, omitting shadowed regions in these calculations seems to have little effect on relative frequency estimates [8, Section 3.F].
Two kinds of relative frequencies may be defined, both with respect to a nominal discrimination threshold.The relative frequency of metamerism is the proportion of times a pair of points chosen at random are indiscriminable under one illuminant but discriminable under another illuminant.The conditional relative frequency of metamerism is the proportion of times a pair of points chosen at random from pairs of points that are indiscriminable under one illuminant are discriminable under another illuminant.
How should these frequencies be estimated for a scene?Given a hyperspectral reflectance image, choose a pair of global illuminants, say daylights of 25,000 K and 4000 K. Choose an approximately uniform color space, CIELAB or CIECAM02, and a nominal threshold color difference ΔE thr .With CIELAB space, apply a chromatic adaptation transform such as CMCCAT2000 to first transform tristimulus values to a reference illuminant such as a 6500 K daylight, as explained in Section 4.B.Then proceed [8] as follows: 1. Choose a random sample of, say, N pairs of points in the scene (the symbol N is used generically).
2. From this sample, find the number of pairs N 0 whose color differences under the 25,000 K daylight are less than the threshold ΔE thr .
3. From this subsample, find the number of pairs N 1 whose color differences under the 4000 K daylight are greater than or equal to ΔE thr (or some multiple, n say, of ΔE thr , the criterion degree of metamerism).
4. The relative frequency of metamerism is N 1 ∕N . 5. The conditional relative frequency of metamerism is N 1 ∕N 0 .
Each of these frequencies is a probability estimate and subject to sampling error.Even so, based on 50 close-up and distant hyperspectral reflectance images of urban and rural scenes under 25,000 K and 4000 K daylights, CIEDE2000 color difference thresholds ΔE thr 0.5 and 1.0, and criterion degrees of metamerism n 1, …, 4 [8], the relative frequency of metamerism appears to be low, from about 10 −6 to 10 −4 , depending on the criterion.In contrast, the conditional relative frequency is very much higher, from about 10 −2 to 10 −1 , sufficiently high to influence visual judgments.Predictably, daylights that have closer correlated color temperatures yield lower relative frequencies [8].
That there should be so large a difference between relative frequency and conditional relative frequency is to be expected, for in a scene with a moderately large gamut of colors, any two surfaces chosen at random are unlikely to match.Thus, for the 50 natural scenes used in Ref. [8], the median CIEDE2000 color difference ΔE between randomly chosen points under a 6500 K daylight is about 19.0, which is much greater than the threshold color difference ΔE thr of 0.5 or 1.0.On the other hand, given a pair of points that do match, i.e., have a color difference less than ΔE thr , they may well fail to match with a change in illuminant.
Broadly similar estimates of the relative frequency of metamerism are obtained by sampling from pooled collections of individual reflectance spectra taken from different natural and manufactured materials illuminated by natural and artificial lights [161].

UTILITY OF COLOR
Basic color properties of a scene such as the number of distinguishable colors and their constancy under changes in illumination do not depend on how often different reflecting surfaces appear, as long as the colors exist and are stable.Yet for color to be of practical use to an observer in distinguishing surfaces and identifying them under illumination changes, their differing relative frequencies do need to be taken into account.
The reason is fundamental.Relative frequencies determine to a greater or lesser extent the composition of any random sample of points from a scene and therefore their distinguishability.How distinguishability should be quantified can be framed in different ways [148], even without introducing cognitive factors [174].For the present purposes, however, information-theoretic methods [162] offer a fruitful guide to what is possible, in particular, the relationship between the number of distinguishable colors in a scene under constant illumination, the number of surfaces or points that can be distinguished by virtue of their color, and the number of points that can be identified under a change in illumination.
To simplify the exposition, technical notions are introduced only when necessary.
The material in this section is based partly on [137,148,175].

A. Entropy and Effective Gamut Volume
Consider the scene in Fig. 6, top right.Histogram estimates of the relative frequency distributions of its CIECAM02 color components are plotted in Fig. 7.The plot on the left is for the lightness component J, in the middle for the red-green component a C , and on the right for the yellow-blue component b C , with the arrows indicating the maximum value of each component.All three histograms are peaked with the long tails in the positive direction for a C and b C due to the orange of the flowers.
These three histograms can be viewed as the marginal estimates of an underlying probability density function f a of a continuous random variable A, say.The color triplets J, a C , b C are the individual instances a of A. The uncertainty in A may be quantified by an information-theoretic measure known as the Shannon differential entropy hA, defined [162] thus: The differential entropy hA is expressed in bits if the logarithm has base 2, which is the base used in this tutorial.As a measure of uncertainty, the differential entropy is unique, up to a constant, according to certain axioms [176].Other kinds of information measures are described in Ref. [177].Methods of evaluating the integral in Eq. ( 10) are described in Section 7.D.It is emphasized that estimates of the differential entropy thereby obtained refer to the underlying continuous distributions, not the discrete hyperspectral samples that represent them.
Given a differential entropy hA, its logarithmic inverse 2 hA (otherwise known as its exponential or antilogarithm) can be interpreted as the volume that defines the effective support set size of the random variable A or, equivalently, of its probability density function f [162,178,179].In the present context, it defines the effective volume of the color gamut, which, as explained in Section 7.B, may be used to estimate the number of points that can be distinguished by their color.Still, as a consistency check, it is reasonable to ask whether for a suitable probability density function f in Eq. ( 10) the definition of effective gamut volume coincides with the conventional definition of gamut volume.
Suppose that the colors in a scene are equally probable, unlike the colors in the scene in Fig. 6, top right.The marginal distributions, in particular, are then flat, by contrast with those in Fig. 7. Let V be the gamut volume.Since the value f a in Eq. ( 10) is constant over all a, its integral is unity, and so f 1∕V .Equation (10) then simplifies to which, after logarithmic inversion, gives 2 hA V : So, as required, for a suitable probability density function, the definitions coincide: effective gamut volume 2 hA equals gamut volume V .
Since the differential entropy of an arbitrary distribution on a bounded set is always less than or equal to the differential entropy of a uniform distribution, the effective gamut volume is always less than or equal to the gamut volume.
As with gamut volume, the effective gamut volume depends on the dimensions of the color space.For the scene in Fig. 6, top right, its differential entropy is approximately 12.5 bits, corresponding to an effective CIECAM02 gamut volume of approximately 2 12.5 5.8 × 10 3 , much smaller than the conventional CIECAM02 gamut volume of 4.7 × 10 4 , calculated by empirically flattening the probability density function f .

B. Number of Distinguishable Surfaces
Differential entropy provides a way to take into account the relative frequencies of colors in a scene.But how does it enable an estimate to be made of the number of surfaces or points that can be distinguished by their color?
As before, treat the CIECAM02 color triplets J, a C , b C as instances of a three-dimensional continuous random variable, A say, and the corresponding observer responses as instances of another three-dimensional continuous random variable, B say.Notation here differs from [153].The amount of information that B provides about A is given [162] by the mutual information, written I A; B. It can be defined in several equivalent ways.The following uses the differential entropies hA and hB and the differential entropy hA, B of the distribution of A and B taken jointly: Because mutual information is a difference of differential entropies, it does not depend on the dimensions in which it is measured, as long as they are related to each other by an invertible linear transformation [162].Importantly, there is a simple relationship between mutual information and the approximate number of distinguishable points that can be sampled from a scene [162,178], and therefore retain their identity in the presence of uncertainty or internal noise in the observer (or color camera).This number, N say, is given by the inverse logarithm of IA; B, i.e., N 2 I A;B : Estimating N requires evaluation of the right-hand side of Eq. ( 12), but it can be simplified.Suppose that observer internal noise is represented by a three-dimensional continuous random variable W such that B A W. With a hard threshold defined by the interval ΔE thr , the probability density function of W is the uniform distribution on a cube with side ΔE thr (alternatives include the uniform distribution on a ball and a three-dimensional Gaussian noise distribution [153]).Providing that W is not too large, then to a good approximation [153], the mutual information in Eq. ( 12) can be written as From Section 7.A, the differential entropy hA is the logarithm of the effective gamut volume, and, by the same kind of calculation that led to Eq. ( 11), the differential entropy hW is the logarithm of the volume of the unit cell ΔE thr 3 .Substituting these values into Eq.( 14) gives Taking the inverse logarithm of both sides and substituting the result into the right-hand side of Eq. ( 13) yields the following for the number of distinguishable points: This expression is precisely analogous to the expression for the number of distinguishable colors in Eq. ( 9).Each is the quotient of a gamut volume by the volume of a unit cell.Indeed Eq. ( 9) is a special case of Eq. ( 15) in which the differential entropy hA is derived for a flattened version of the probability density function f , as in Section 7.A.
The number of distinguishable points can now be estimated for the scene in Fig. 6, top right.From Section 7.A, the differential entropy of the colors is approximately 12.5 bits, corresponding to an effective CIECAM02 gamut volume of approximately 5.8 × 10 3 .If this volume is divided by the volume of the unit cell ΔE thr 3 , with ΔE thr 0.6, then by Eq. ( 15), the number of distinguishable points is approximately 2.7 × 10 4 .This is much smaller than the estimate of 2.2 × 10 5 obtained by taking its CIECAM02 gamut volume of 4.7 × 10 4 (Section 7.A) and dividing by the same volume of the unit cell ΔE thr 3 .
These estimates are for only one scene.Values for 50 natural scenes are given in Ref. [148,Fig. 4].
In other applications, differential entropy is used to predict the relative frequency and magnitude of metamerism across scenes and the changes in relative color cues or color relations between surfaces under changes in illuminant [148,180].In its discrete form [162], entropy is used to classify observers' colornaming behavior with individual colored surfaces [181].Mutual information is additionally used to estimate the relationship between differences in luminance and differences in color within natural scenes [145] and to estimate the optimum spectral locations of cone photopigments [182].

C. Number of Illuminant-Invariant Points
A utilitarian way of measuring color constancy in a scene is outlined in Section 6.C.2, namely, as the number of points that preserve their identity under a change in illuminant.But how to estimate this number?Consider once more the scene in Fig. 6.One approach is to perform an exhaustive nearestneighbor search for color matches, i.e., by counting all those points whose transformed CIE tristimulus values X 2 , Ỹ 2 , Z 2 are closer to their true values X 2 , Y 2 , Z 2 than for any other points in the scene, with respect to some color-difference formula.But the result depends on the chromatic adaptation transform.Different transforms define different nearest-neighbor matches; moreover, nearest-neighbor matches themselves may not always be optimum [183].
Instead, information-theoretic methods may be used in a way similar to that for estimating the number of distinguishable points under constant illumination described in Section 7.B.Such an approach provides a least upper bound on all possible chromatic adaptation transforms and more generally on any color-constancy operation.
Again treat the CIECAM02 color triplets J, a C , b C as instances of a three-dimensional continuous random variable, A 1 say, under the first illuminant and as a different threedimensional continuous random variable, A 2 say, under the second illuminant.Suppose that f 1 a 1 and f 2 a 2 are the corresponding probability density functions of A 1 and A 2 .The amount of information that A 2 retains about A 1 is measured by the mutual information I A 1 ; A 2 .In Section 7.B, it is used to determine the number of points that preserve their identity in the presence of observer internal noise.Here, it determines the number of points that retain their identity in the presence of what might be thought of as illuminant-change noise.This number, N say, is given by the inverse logarithm of For the scene in Fig. 6, top, under 4000 K and 6500 K daylights, the estimated mutual information IA 1 ; A 2 , calculated as in Sections 7.A and 7.B, is approximately 21.9 bits, corresponding to approximately 4.0 × 10 6 illuminant-invariant points.
As before, these estimates are for one scene.Means and standard deviations for 50 natural scenes under different daylight illuminants are given in Ref. [137, Fig. 3(a)].
Although I A 1 ; A 2 does not depend on linear transformations of the color space, it does depend on the type of color space, since the relationship between one space and another is usually nonlinear.In fact, with the chromatic adaptation transform CMCCAT2000 and color-difference formula CIEDE2000, the CIELAB estimate of IA 1 ; A 2 for the scene in Fig. 6 is approximately 22.4 bits, just 0.5 bits larger than with CIECAM02.The corresponding number of illuminant-invariant points is approximately 5.4 × 10 6 .
These estimates of the number of illuminant-invariant points for the scene in Fig. 6 are large but plausible, since they exclude observer uncertainty.To help interpret the magnitude of the illuminant-change noise, it can be expressed in terms of equivalent Gaussian variables [137].Suppose that the noise is represented by a three-dimensional continuous random variable V, such that A 2 A 1 V. Suppose further that the J, a C , b C components of A 1 and of V are distributed normally and independently with respective common variances σ 2 and σ 2 V .Then, as shown in Ref. [162], the mutual information I A 1 ; A 2 3∕2 log1 σ 2 ∕σ 2 V .It follows that σ∕σ V 2 2IA 1 ;A 2 ∕3 − 1 1∕2 ≈ 2 IA 1 ;A 2 ∕3 .Therefore, given that I A 1 ; A 2 is approximately 21.9 bits, the equivalent Gaussian signal-to-noise amplitude ratio σ∕σ V is approximately 158, or, conversely, the percentage noise amplitude 100 σ V ∕σ is about 0.63%.This value is smaller than measures of observer internal noise such as the Fechner fraction, which is on the order of 2% [92].
Suppose, now, that the estimate of the number of illuminant-invariant points includes observer uncertainty, as in Section 7.B.The mutual information I A 1 ; B between A 1 and the observer response B to A 2 decreases, since there are now two sources of noise that are largely independent of each other, one due to an illuminant change and one due to the observer, i.e., I A 1 ; B < IA 1 ; A 2 : Since the number of illuminant-invariant points N is given by N 2 I A 1 ;B , it too decreases.
Thus, to continue with the example of the scene in Fig. 6, top, suppose that the probability density function for observer internal noise W is the uniform distribution on a cube of side ΔE thr .Then, the estimated mutual information I A 1 ; B, calculated as in Sections 7.A and 7.B, falls to approximately 11.8 bits, corresponding to approximately 3.5 × 10 3 illuminant-invariant points, three orders of magnitude smaller than without observer uncertainty.Similar estimates are obtained by replacing the uniform distribution for observer internal noise by a three-dimensional Gaussian noise distribution with the same covariance.
In addition to spectral changes, non-spectral changes in illumination, e.g., changes in the distribution of shadows, specularities, and mutual reflections, along with other uncertainties (Section 3.C) can be interpreted as sources of noise, leading to further reductions in I A 1 ; A 2 and in I A 1 ; B. The complexities of environmental illumination may also affect observer performance in other ways [184].
Notwithstanding the usefulness of mutual information and its equivalents, can colorimetric measures such as the color constancy index provide a similar guide to identifying surfaces or points under illuminant changes?Unfortunately, global statistics based on the mean or median of point-wise color differences do not suffice.The reason is evident with the scene in Fig. 6.With the chromatic adaptation transform CMCCAT2000, the tristimulus values X 1 , Y 1 , Z 1 of the sample point arrowed in yellow, bottom left, are transformed to the values X 2 , Ỹ 2 , Z 2 , which, as noted in Section 6.C.1, differ from the true values X 2 , Y 2 , Z 2 of the sample point, arrowed in yellow, bottom right.These estimates X 2 , Ỹ 2 , Z 2 are actually closer to the tristimulus values of each of the two points arrowed in white (and to the tristimulus values of many other points not indicated in the image), according to the colordifference formula CIEDE2000.A more detailed example is described in Ref. [137].
These kinds of erroneous matches occur with other points in the scene and with other chromatic adaptation transforms.Patently, identifying a point by its color involves not only the color difference between its estimated and true tristimulus values X 2 , Ỹ 2 , Z 2 and X 2 , Y 2 , Z 2 , or their equivalents, but also the color differences between X 2 , Ỹ 2 , Z 2 and the tristimulus values of all the other points in the scene.
Although mutual information sets a least upper bound on the number of illuminant-invariant points in a scene, it does not automatically follow that it is possible to find a chromatic adaptation transform that can reach that bound, even with illumination changes that are purely spectral.In practice [137,175], performance tends to fall several bits below the least upper bound.Locally weighted polynomial regression may offer a better description [185], though it risks overfitting [107].

D. Estimation
Numerical estimates of differential entropy may be obtained directly from probability density functions, e.g., from empirical histograms.Examples of histogram estimates of the marginal distributions are shown in Fig. 7.But using histograms in place of the unknown probability density function f in Eq. ( 10) can lead to marked biases in the estimates [186].Fortunately, differential entropy can be estimated from nearest-neighbor statistics, a safer approach than constructing empirical histograms.The Kozachenko-Leonenko kth-nearest-neighbor estimator [187,188] is asymptotically unbiased and converges reasonably rapidly with increasing sample size.Its convergence is improved with an offset method, details of which are given in Ref. [175].The essence of the method is to decompose the estimate into two components: one the mutual information between equivalent Gaussian variables with known variance-covariance structure, and the other an offset obtained by applying the nearest-neighbor estimator to normalized versions of the variables.
The Kozachenko-Leonenko estimator can be implemented by library routines with efficient C++ code for exact nearestneighbor search [135,136].

COMMENT AND CONCLUSION
Terrestrial and close-range hyperspectral imaging makes possible the realistic analysis and modeling of color vision in the environment.Cone receptor and postreceptor responses can be estimated, scene color statistics calculated, and the utility of color assessed in requiring an observer to distinguish surfaces or identify them under changes in illumination.In addition, colorimetric formulas can be tested and color rendering transformations optimized.
Caution does, though, need to be exercised in securing images and in analyzing them: 1.To be valid, hyperspectral images should be linearly related to the physical quantities they represent.The workflow for image acquisition, processing, and calibration demands careful, accurate management.
2. The choice of spectral radiance or reflectance images leads to different tradeoffs.Spectral reflectance data allow a wider range of questions to be addressed but entail more assumptions or constraints than spectral radiance data.
3. In modeling receptor responses to images, variations in receptor density and type may need to be considered, along with the intrusion of signals from rods and intrinsically photosensitive retinal ganglion cells.
4. In assessing the utility of color in different observer tasks, average measures of scene properties such as the mean or median of color signals may be less relevant than information measures such as Shannon differential entropy.But estimating differential entropy with little or no bias depends on the goodness of the chosen numerical estimator.
More generally, each spectral radiance or reflectance image constitutes a single sample from the real world in which reflecting properties may vary from point to point and reflected light from instant to instant.Inferences drawn from such data are inevitably probabilistic.
Similar caveats apply to extensions of hyperspectral imaging, including stereo, goniometric, and video hyperspectral imaging.Constraints on acquisition time may force reductions in the spatial or spectral sampling density available to the imaging system.Nevertheless, preserving the relevant degrees of freedom associated with natural scene data remains a prerequisite.
Fig.1.Examples of reflected radiance spectra from a flower scene.Data are from a hyperspectral radiance image of size 1344 × 1024 pixels, corresponding to approximately 14 × 11 deg visual angle at the camera, with spectra sampled at 400 nm, 410 nm, …, 720 nm.The plots show radiance spectra at individual pixels (radiance scales adjusted for range).The small light-gray sphere near the top of the scene is covered in Munsell N7 matte paint[4], and the reflected spectrum (top right plot) follows the typically uneven spectrum of light from the sun and sky[5].The long, thin, light-gray rectangular plate at the bottom of the scene is a reference reflectance surface.The hyperspectral image of the scene was acquired on October 10, 2003, from a garden in Sameiro in the Minho region of Portugal.

Fig. 2 .
Fig.2.Grayscale wavelength slices from the hyperspectral radiance image of the flower scene in Fig.1, with wavelength indicated in nm at the top left of each slice.The intensity range in each slice image has been stretched for illustration.As evident from the spectral radiance plots in Fig.1, the surfaces in the scene reflect little energy at very short wavelengths.

Fig. 3 .
Fig. 3. Effect of registration across wavelength.The color images of the stone cottage are rendered from a hyperspectral radiance image without registration, top left, and with registration, top right.Color fringing due to lateral chromatic aberration and its removal by registration can be seen more easily in the enlarged copies of the square areas marked in white, bottom left and bottom right, respectively.The long, thin, light-gray rectangular plate at the bottom of the full scene is a reference reflectance surface.The hyperspectral image of the cottage was acquired on June 4, 2003, under an overcast sky in Ruivães in the Minho region of Portugal.

1 .
Spectral Illumination Changes An example of where spectral illumination effects dominate is shown in Figs.4(a) and 4(b).The color images are rendered from time-lapse hyperspectral radiance images of a rock face acquired in the late afternoon, about 25 min apart.Although there are small differences in the distribution of shadows in (a) and (b), the obvious change is in the shift of the mean scene color towards red, reflecting the change in prevailing illumination.

Fig. 4 .
Fig. 4. Real and simulated changes in illumination on a rock face.The color images in (a) and (b) are rendered from two time-lapse hyperspectral radiance images acquired at 18:15 and 18:40 [93].The image in (c) is rendered from a simulated version of the spectral radiance image in (b) in which a global illuminant is applied to an effective spectral reflectance image derived from the spectral radiance image in (a).The real image in (b) and the simulated one in (c) are closely similar.The images in (d), (e), and (f ) are analogous, except that (d) and (e) are from earlier in the day, at 13:21 and 15:15, and show marked changes in the distribution of shadows.The real image in (e) and the simulated one in (f ) are clearly different.The small light-gray rectangular plate at the bottom of the scene, arrowed in (a), is a reference reflectance surface.The hyperspectral images of the rock face were acquired on October 6, 2003, in Sete Fontes in the Minho region of Portugal [93].
4(d)  and 4(e).The color images are rendered from time-lapse hyperspectral radiance images of the same rock face as in (a) and (b) but acquired at two earlier times of day.The correlated color temperature of the light incident on the rectangular plate at the bottom right of the scene, arrowed in (a), decreases from 5700 K at 13:21 in (d) to 5510 K at 15:15 in (e).The spectral difference is negligible, but there are marked changes in the distribution of shadows.Can the reflected spectral radiances in Figs.4(d) and 4(e), like those in Figs.4(a) and 4(b), still be adequately represented as products with a common effective spectral reflectance?The analysis follows Section 3.B.1.Suppose that E 1 λ and E 2 λ are global illuminant spectra with correlated color temperatures 5700 K and 5510 K, respectively, and that L 1 u, v; λ and L 2 u, v; λ are the corresponding reflected spectral radiances.
Next, sRGB tristimulus values R, G, B are obtained by linear transformation:

Fig. 5 .
Fig. 5. Color rendering of a hyperspectral radiance image of a yellow flower.An sRGB image is shown left and a clipped and scaled version shown right, with the clip level taken from the arrowed area on the light-gray sphere in the left image.The percentage of clipped pixels is about 1%.The hyperspectral image of the flower was acquired on July 31, 2002, in Gualtar in the Minho region of Portugal.

6 .
Physical limits on color constancy under a global illuminant change.The images top left and top right are color renderings of a hyperspectral reflectance image of a terrace with flowers under a global daylight illuminant with respective correlated color temperatures 4000 K and 6500 K. Enlarged copies of the square areas marked in white are shown bottom left and bottom right.By applying a standard chromatic adaptation transform CMCCAT2000 [110], the point on the flower arrowed in yellow, bottom left, transforms to a point that is a closer color match to each of the two points arrowed in white, bottom right, than to the correct point arrowed in yellow.The small light-gray rectangular plate at the bottom of the full scene is a reference reflectance surface.The hyperspectral image of the terrace was acquired on October 7, 2003, in Sameiro in the Minho region of Portugal.

Fig. 7 .
Fig. 7. Histogram estimates of the relative frequency distributions of the lightness component J, left, the red-green component a C , middle, and the yellow-blue component b C , right, of the CIECAM02 representation of the scene in Fig. 6, top right.The arrows indicate the maximum values of the components, which have very small relative frequencies.
6.It shows a color rendering of a hyperspectral reflectance image of a terrace with flowers under a 4000 K daylight, top left, and under a 6500 K daylight, top right.Enlarged copies of the square areas marked in white are shown bottom left and bottom right.Suppose that under the 4000 K daylight, the sample point on the flower, arrowed in yellow, bottom left, has CIE tristimulus values X 1 , Y 1 , Z 1 and that under the 6500 K daylight, the same point arrowed in yellow, bottom right, has tristimulus values X 2 , Y 2 , Z 2 (the white arrowed points are discussed later).Suppose, next, that as the result of some color-constancy operation applied to a trichromatic representation of the scene under the 4000 K daylight, the sample point takes on tristimulus values X 2 , Ỹ 2 , Z 2 under the 6500 K daylight.In general, these values only approximate the true values X 2 , Y 2 , Z 2 .To quantify the error (Section 4.B), transform X 2 , Ỹ 2 , Z 2 and X 2 , Y 2 , Z 2 to corresponding CIELAB values L