Designing Color Filters that Make Cameras MoreColorimetric
Graham D. Finlayson Yuteng Zhu
cs.CV










Designing Color Filters that Make Cameras More Colorimetric
    Graham D. Finlayson (g.finlayson@uea.ac.uk)
and Yuteng Zhu* (yuteng.zhu@uea.ac.uk) 

(co-first authors) 

Graham D. Finlayson and Yuteng Zhu are with the School of Computing Sciences, University of East Anglia, NR4 7TJ, Norwich, UK (g.finlayson@uea.ac.uk, corresponding author: yuteng.zhu@uea.ac.uk).
    December 30, 2023
================================================================================================================================================================================================================================================================================================================




When we place a colored filter in front of a camera the effective camera response functions are equal to the given camera spectral sensitivities multiplied by the filter spectral transmittance. In this paper, we solve for the filter which returns the modified sensitivities as close to being a linear transformation from the color matching functions of human visual system as possible.  When this linearity condition - sometimes called the Luther condition - is approximately met,  the  `camera+filter' system can be used for accurate color measurement. Then, we reformulate our filter design optimisation for making the sensor responses as close to  the CIEXYZ tristimulus values as possible given the knowledge of real measured surfaces and illuminants spectra data. 
This data-driven method in turn is extended to incorporate constraints on the filter (smoothness and bounded transmission). Also, because how  the optimisation is initialised is shown to impact on the performance of the solved-for filters, a multi-initialisation optimisation is developed. 

Experiments demonstrate that, by taking pictures through our optimised color filters we can make cameras significantly more colorimetric.



Color filter, camera sensitivity functions, color measurement 



 EDICS Category: ELI-COL, ELI-SDP 

 




§ INTRODUCTION


Digital cameras are designed in analogy to the trichromatic human visual system which has three cone sensors. If a camera is to capture colors like a human observer, arguably the camera sensors should equal the  cone fundamentals <cit.>. Practically, however, engineering cameras having spectral sensitivities similar to the cone fundamentals is only required  if we wished to construct a biologically plausible model of how we see<cit.>.  For most practical applications - e.g. photography and video - it is more important that we can transform the recorded device RGBs to drive a display so that the image captured by a camera either looks the same to a human observer or records triplets of numbers - e.g. CIE XYZ coordinates - that are referenced to the human visual system<cit.>. We say that a digital camera is colorimetric if it meets the so-called Luther condition <cit.>, i.e. its spectral sensitivity functions are linearly related to the CIE XYZ color matching functions (CMFs).

The Luther condition places a very strong constraint on the shape of camera spectral sensitivities. A strong constraint is required  because the Luther condition effectively assumes that any and all spectral stimuli are possible. However, many studies have shown that the actual spectra (measured in the real world) are far from being arbitrary. Indeed, reflectance spectra tend to be quite smooth<cit.> and as a consequence can be fit with low dimensional linear basis<cit.>. Indeed, spectral basis with dimensions from six to eight, for different applications, are often proposed as adequate models of spectral reflectance. Illuminants by contrast are much less describable by small parameter models. Indeed, for artificial lights such as fluorescent and LED lights, the light spectrum can be very spiky and the number and position of the spikes can vary considerably. And yet, illumninants are also far from being arbitrary. They are designed to have colors near the Planckian locus<cit.>, a requirement to score highly on color rendering indices <cit.>.

Possibly, a more practically useful variant of the Luther condition would be one that is data-driven. That is, where camera RGBs can be mapped to XYZs for the spectral data that are likely to be encountered in the real world. Equally, in principle, we might consider whether a non-linear mapping could or indeed, should be used.



It is a classical result <cit.> that if reflectances were exactly modelled by a 3-dimensional linear model then for a given light spectrum there would be a specific 3 × 3 transform matrix  taking camera RGBs to XYZs. While reflectance spectra are not adequately described by a 3-dimensional model, RGBs can be approximately mapped to corresponding XYZs using a 3 × 3 matrix.  Indeed, this regression approach is adopted in almost all cameras with good results (we are mostly happy with the colors a camera records). But, as we shall see later the `fit' is not sufficient from a color measurement point of view. 

Of course, rather than using a linear matrix to  map  RGBs to XYZs, we could  use a non-linear transform instead. Possible non-linear methods include Polynomial and Root-polynomial  regressions  and Look-up-tables <cit.>.
However, the linear transform method of using a 3 × 3 matrix - even though it is not optimal in terms of fitting error -  has two  advantages compared to most non-linear methods.  First,  the transform scales linearly with exposure. If the scene is made twice as bright (e.g. by doubling the quantity of incoming light), the same matrix correctly maps the camera measurements to XYZs (because the magnitude of camera RGBs and  XYZs both double). Typically, non-linear methods do not have this exposure-invariant property (one exception is <cit.>).

The second advantage is that a linear transform is, well, linear. The human eye measures color stimuli linearly: at the cone quantal catch level, the response to the sum of two spectral stimuli is equal to the sum of the responses to the two stimuli viewed individually <cit.>. This can be an important physical consideration. As an example, when we view a surface that has highlights, the recorded color is a convex sum of the so-called body color (the color name we would assign to the object) and the color of the highlight<cit.>, sometimes called the interface color. If, we are viewing a red shiny plastic surface the body color is red and if the viewing light is white then the interface color is also white (i.e. the same as the color of the light). As we move from pure body to pure highlight color, the measured XYZs lie on a 2D plane in the color space. Equally the camera, which at the sensor level has a linear response, will also make measurements that lie on a 2D plane. But, a non-linear correction will distort the plane and the result will be an image that is not physically accurate or even physically plausible. This problem is discussed in detail in <cit.>.



Practically, the closer the spectral sensitivities of a camera are to being linearly related to the XYZ color matching functions the better it will perform as a tool for color measurement, i.e. the more colorimetric it will become. Interestingly, when we linearly regress RGBs to XYZ tristimuli, we can interpret this as linearly transforming the sensors themselves. That is, a new camera whose sensitivities are modified by the linear regression transform approximately measures the desired XYZs. It follows that one strategy to improving the color measurement capability of a camera would be to change the camera sensitivities (to ones that are more linearly related to the XYZ color matching functions). However,  there are constraints on sensitivities achievable in physically realisable cameras which means we can never get 100% colorimetric accuracy. 

In this paper, we make an easy modification to the camera spectral sensitivities. We propose simply to place a specially designed filter in front of a camera with the goal of making the filtered RGB measurements it records are as linearly related to the target XYZs as possible. 


How we design the best filter to make a camera colorimetric is the central concern of this paper. 
We develop optimisation methods to solve for the optimal filter that either makes the camera best satisfy the Luther condition or - in a data-driven approach - can best predict measured XYZs for a range of real reflectances and lights. In the rest of this paper we respectively discuss Luther-condition and Data-driven filters.

The methods we develop are not closed-form but adopt the alternating least-squares paradigm<cit.>. For the Luther-condition filter optimisation problem we aim to find a filter so that  the filtered camera response functions are a linear transform from the XYZ CMFs <cit.>. In the first step, we find the filter that best maps the spectral sensitivities of the camera to the XYZ CMFs directly. Then we find the best linear combination of the filtered camera sensitivities that approximate the XYZ sensitivities. Holding this mapping fixed we can solve for a new  best filter. Then we hold the new filter fixed to solve for the best linear transform. We iterate in this way until the procedure converges <cit.>. Each individual step in the optimisation can be solved, in closed-form, using simple linear least-squares.
In the Data-driven approach we analogously find the filter based on actual measured RGBs and XYZs following the alternating least-squares technique <cit.>. For the Luther- and Data-driven techniques, the constraint that the recovered filters must be positive is considered.



 
Clearly, the filter shown in Fig. 1e is not desireable. In the short wavelengths there is a sharp change in transmittance and as a whole the filter is not smooth. For most of the wavelength range the filter transmits little light (<20%).
Thus, we extend our optimisation framework to incorporate minimum and maximum bounds on transmittance and also that the filters are smooth <cit.>.


Experiments demonstrate that  we can find Luther-condition and Data-driven  filters that can dramatically increase the colorimetric accuracy across a large set of commercial cameras.

The rest of the paper is structured as follows. In Section II we review the color matching and color image formation, both ideas underpin our filter design method. The mathematical optimisations for the optimal color filter for a given camera are presented in Section III. The experimental results are reported in Section IV. The paper concludes in Section V.




§ BACKGROUND




 §.§ Color Matching Functions

Color matching functions  provide a quantitative link between the physical light stimuli and the colors perceived by the human vision system. 
Figure <ref> shows a typical  setup for the color matching experiment. The observer views a bipartite field where one side is lit by a test light while the other side is lit by the light mixtures of three primaries (i.e. monochromatic red, green, blue lights). The intensities of three primary lights are adjusted  by the observer  to make a visual match, i.e. the two stimuli on each side  the bipartite field match if they look visually indistinguishable  to the observer. Sometimes, no match is possible. In this case one of the primary lights should be added to the test light field. Mathematically, we can model this as if we were subtracting some of the primary light. See <cit.> for more discussion.

By successively setting the test light to be a unit power monochromatic light at sample points across the visible spectrum, we can measure the Color Matching Functions<cit.>. That is we find the  RGB mixture that affords a match on a per wavelength basis. Because color matching is linear (the sum of two test spectral lights is matched by the sum of their individual matches) then given the color matching functions we can compute the match for any arbitrary test light <cit.>. We simply integrate any test light spectrum with the Color Matching Functions. It can also be shown that the color matching functions are necessarily linearly related to the cone sensitivities <cit.>. 

Assuming monochromatic primary lights at 650nm, 530nm and 460nm, the resulting CIE RGB CMFs are shown to the right of Fig. <ref>.
The XYZ CMFs are a linear combination of the RGB CMFs, see in Fig. 1b. That XYZ CMFs are used (as opposed to RGB CMFs), in part, because they have, by design, no negative sensitivities. Standardised in 1931, the lack of negatives made pencil and paper calculations with matching curves easy. 

The X, Y and Z scalar values we compute when we integrate a test light spectrum with the XYZ CMFs are called XYZ tristimuli.
In this paper we are interested in using a camera to measure XYZ tristimuli.

The color matching experiment is of direct practical importance. Indeed, suppose a display device can output three colors equal to the R, G and B stimuli used in the color matching experiment. It follows that a camera that had sensitivities equal to RGB Color Matching Functions would see the correct RGBs to drive the display that would result in a perceptual match to the test light.  Equally, it would suffice that camera's sensitivities were linearly related to the CMFs since we could linearly transform the camera measurements to the correct RGB to drive the display.  Of course, we note that - unlike the RGB color matching functions - we can only drive the putative display with positive numbers. Consequently, there are  real world colors that cannot be reproduced on displays. 



 §.§ Continuous Color Image Formation


The color of a pixel recorded by a digital camera  depends on three physical factors: the spectral power distribution of the illuminant, the spectral reflectance of the object, and the sensor's spectral sensitivities. The color formation at a pixel, under the Lambertian surface model, can be modeled as
    
    ρ_k = ∫_ωE(λ) S(λ)Q_k(λ) λ ,    k ∈{R, G, B}

where the subscript k indicates the color channel, λ denotes wavelength variable defined on the visible spectrum  ω (approximately, 400nm to 700nm). The functions  E(λ) and S(λ) respectively denote the spectral power distribution of the illuminant and the spectral reflectance. The function  Q_k(λ) is the spectral sensitivity of the kth camera color channel.

Let us define the color signal C(λ):


    C(λ) = E(λ)S(λ).

Substituting Equation (<ref>) into Equation (<ref>):
    
    ρ_k = ∫_ω C(λ)Q_k(λ) λ ,    k 		∈{R, G, B}.

    


 §.§ Discrete Color Image Formation


In the optimisations that will be presented in Section III, it is useful to recast the continuous integrated responses using discrete representations where spectra are represented as sampled vectors of measurements. Typically and justifiably from a color accuracy point of view<cit.>, a spectrum can be represented by 31  measurements made between 400nm and 700nm by 10nm sampling interval. Note, the methods we set forth - since they are designed for the discrete domain - are agnostic about the sampling interval. If the data is given at a 5nm sampling distance then each spectrum would be represented as a 61-component vector. Henceforth,  we will talk about spectra being 31-vectors (and this corresponds to the format of most available measured spectral data).

Equation (<ref>) is now, equivalently, written as:

    ρ_k= C·Q_k,   k ∈{R, G, B}.

Respectively C and Q_k denote the sampled version of the color signal spectrum and the kth spectral sensitivity function. We assume the sampling distance is incorporated in the spectral sensitivity vectors. Here `·' denotes the dot product of two vectors.

One advantage of transcribing image formation into vector notation is that we can usefully deploy the tools of linear algebra. Let the  matrices C and Q denote 31× n and 31× 3 matrices whose columns respectively contain n color signal spectra and the 3 device spectral sensitivities. The n× 3 set of RGB responses can be written as a single concise expression:


    P = C^TQ,

where ^T denotes the matrix transpose.

Denoting, X as a 31× 3 matrix whose colmuns contain the discrete XYZ color matiching functions, the  XYZ tristimulus responses can be written as: 

    X = C^T X.




 §.§ The Luther Condition for Camera Spectral Sensitivities


Let us denote individual camera and XYZ responses as ρ and x. We can use a camera to exactly measure colors if and only if there exists a function such that g(ρ)=x for all spectra. It follows that if a pair of spectra  integrates to  the same RGB then this pair must also integrate to the same XYZ:


    C^T_1Q= C^T_2Q  ⇒ C^T_1 X= C^T_2 X

 
This implies C^T_1-C^T_2 is simultaneously in the null space of Q and X. Since any spectrum in the null-space of Q is a physically plausible spectral difference this implies that the null-spaces of Q and X are the same and this in turn implies the Luther condition 


    X=QM


 where M is a 3× 3 matrix. See <cit.> for the original proof of the Luther condition (which we precis above). 
 
 

 §.§ The Data-driven Luther Condition 

 
The Luther condition for Spectral Sensitivities  presupposes that any and all physical spectra are plausible. However, we know that reflectance spectra are smooth<cit.>. Lights while more arbitrary are designed to integrate to fall near or close to the Planckian locus<cit.> and to score well on measures such as the Color Rendering Index<cit.>. Pragmatically, a camera is colorimetric if its responses are a linear transformation from the XYZ tristimulus values

    C^T X = C^TQM

where M is a 3× 3 matrix.




 §.§ Filter Design


There are many papers in the literature where the spectral sensitivities that a camera should have had are designed <cit.>. For example, we can solve for the camera spectral sensitivities with respect to which the RGBs mapped with respect to many illuminant spectra can be mapped to corresponding XYZs for a single target illuminant. The procedure where we  measure colors under a changing light - e.g. in the real world -  but then reference (that is, map) these colors back to a fixed illuminant viewing condition is a standard methodology in color science. Curiously, the best sensors that solve this problem are not linearly related to the XYZ CMFs <cit.>.


Perhaps, the closest work to our study is <cit.>. Here two images are captured. The first with the native spectral sensitivities and the second through a colored filter. The emphasis of that work was to increase the spectral dimensionality of the capture device.  Since we make two captures we have 6 measurements per pixel. Effectively we have a 6-dimensional sensor set. We match target XYZs by applying a 6× 3 correction matrix. In <cit.>, the best filter was  chosen from a  set of commercially available filters. 

There are many other works which propose recovering spectral information by capturing multiple exposures of the same scene through multiple filters, e.g. <cit.>. A disadvantage of all these methods is that the capture process is longer and more laborious. Filters need to be changed between exposures. The multi-exposure process then need to be registered. Image registration remains a far from solved problem.  Moreover, scene elements between exposures may move (making registration impossible). 

The method we propose here is much simpler. We simply place a specially designed filter in front of a camera and then take conventional single exposure images.




§ DESIGNING A FILTER TO MAKE A CAMERA MORE COLORIMETRIC
 




 §.§ Luther-condition Filters


The Luther condition states that a camera system is colorimetric if its sensitivities are a linear transform from the XYZ color matching functions. We propose a modified  Luther condition where a camera is said to be colorimetric if there exists a physically realisable filter which, when placed in front of the camera, generates effective sensitivities which are a linear transform from the XYZ CMFs. 

Physically, the role of a filter, which absorbs light on a per wavelength basis, is multiplicative. If f(λ) is a transmittance filter and Q(λ) denotes the camera sensitivities then f(λ)Q(λ) is a physically accurate model of the effect of placing a filter in front of the camera sensor.

In Equation (<ref>) we write an optimisation statement for the Filter-based Luther condition:


    min_f,M∥diag(f)QM -  X∥^2_F s.t. f>0

Here Q and X are 31 × 3 matrices encoding respectively the spectral sensitivities of a digital camera and the CIE standard XYZ color matching curves. The 31-vector f is the sampled equivalent of the continuous filter function  f(λ). The diag() operation converts a vector into a diagonal matrix (the vector components appear along the diagonal). The meaning of diag(f)Q  is the same as f(λ)Q(λ), i.e. the diagonal matrix allows us to express component-wise multiplication. M denotes a 3× 3 matrix. We minimise the square of Frobenius norm ∥ ∥^2_F (we minimise the sum of squares error). Notice the constraint that the filter value is larger than 0  (physically, we cannot have a filter that has negative transmittance).
 
We do not have to constrain the maximum transmittance because
we can only solve for f and M up to an unknown scaling factor. Indeed, suppose the filter f is returned where the max transmittance is larger than 1. The fitting error in Equation
(<ref>) is  unchanged if we divide f by its maximum value (resulting in a max transmittance of 100%) so long as we multiply the corresponding correction matrix M by the same value.


We minimise Equation (<ref>)  using an  alternating least-squares (ALS) procedure given in the following algorithm.



where ∏ denotes component-wise matrix multiplication. 
Steps 4 and 5 - where we find the filter and then the linear transform - are solved using simple, closed-form least-squares estimation. For completeness we provide details of how these calculations are made in the Appendix.


At each iteration, the filter and linear transform -  f^i and M^i - are calculated relative to the previous i-1 filters and matrices. It follows in step <ref> that the final solution is the multiplication of all the per-iteration solutions: f^Luther = ∏_s=1^i f^s (component-wise vector multiplication) and M_j = ∏_s=1^i M^s (component-wise matrix multiplication).


Notice nowhere in the above procedure do we constrain the filter transmittance to be larger than 0 (even although this constraint is in the optimisation statement). Empirically, we found that the optimised filter is always positive for all the cameras we tested (see experimental section). Moreover, Theorem 3.1, presented below, proves that when there exists a filter which makes the camera sensors perfectly colorimetric that the filter has to be everywhere positive.

The theorem is presented for continuous spectral sensitivity functions. As such, we write the XYZ CMFs and camera sensitivities as vector functions: X(λ) and Q(λ). In this representation, we, effectively, have taken transposes of the matrices Q and X. So, here, we write M^T for the 3 × 3 matrix. Of course, matrix M in the proof and in the algorithm presented above signifies the same linear transform.





Assuming there exists an exact solution that f(λ)>0 for M^TQ(λ)f(λ)= X(λ), the variable λ is defined over the domain where  X(λ)>0 are continuous  and full rank (no one spectral  sensitivity  can be written as a linear combination of the other two) and Q(λ) are also continuous then  f(λ)>0.






 Proof: First we remark on the continuity of the camera and XYZ functions. Both are the result of physical processes which are continuous in nature. To our knowledge it is not possible to make a physical sensor system that captures light which has discontinuous sensitivities. And, in terms of physiological systems, biological sensor response functions are always continuous. 
 
 Next, if f(λ)<0 across all wavelengths and   M^TQ(λ)f(λ)= X(λ), then -M^TQ(λ)*(-f(λ))= X(λ). In this case -f(λ) must be all positive and so an all-positive filter can be found. The interesting case to consider is when the filter has both negative and positive values.


Clearly the 3× 3 matrix M^T must be full rank otherwise the mapped camera sensitivities would be rank deficient and therefore could not model the CMFs. Equally,  multiplying by a filter does not change the rank of the sensor set. Because, by assumption Q(λ) are continuous it follows that f(λ) must also be a continuous function since otherwise M^TQ(λ)f(λ) would be discontinuous (multiplying a continuous and discontinuous functions together, save for the case where one of the functions is everywhere 0, results in a discontinuous function).


As f(λ) is continuous if the function has both negative and positive values there must be at least one wavelength  λ_v where f(λ_v)=0 and so M^TQ(λ_v)f(λ_v)=0. But, this cannot be the case since the XYZ color matching functions are not all zero at any given wavelength within the defined domain.

 QED



 §.§ Data-driven Filters
 

Simple Case: in the simple Data-driven approach, we look for a  color filter and the 3 × 3 color correction matrix that, in a least-squares sense, best maps camera measurements for a training color signal data set to the corresponding ground-truth XYZ tristimulus values. Denoting a collection of n color signals in the n× 31 color signal matrix C,  the Data-driven optimisation is written as:

    min_f, M∥ C^Tdiag(f)QM - C^T X∥_F^2  s.t. f>0


Solving Equation (<ref>) depends on the structure of the color signal matrix. If we choose C= I_31 (the 31× 31 identity matrix) then we can solve this optimisation using Algorithm 1 (in this case, Equations (<ref>) and (<ref>) are the same). This assumption is related to the Maximum Ignorance assumption <cit.> where all possible spectra are considered equally likely. 



General Case:
Let us develop a more general optimisation statement. One,  where we have cnt  color signal matrices - denoted C_j (j∈{1,2,⋯, cnt}) and the corresponding cnt color correction matrices M_j. Each color signal matrix typically corresponds to a training set of surface reflectances illuminated by a single spectrum of light E_j, thus the color signal matrix is C_j=diag(E_j)S, where S is a 31× n matrix of reflectances, one reflectance spectrum per column. But, the different light assumption is not a necessary assumption. As an example, we might mix color signals for the Maximum Ignorance assumption with measured data (i.e. two color signal matrices) where both measurement sets contain multiple lights. 

The general Data-driven filter optimisation problem is written as:


    min_f, M_j∑_j=1^cnt∥ C^T_jdiag(f)QM_j - C^T_k X∥_F^2   s.t. f>0


and is solved using Algorithm 2.


Finally,  k could denote some other privileged standard reference viewing condition (where the reference viewing illuminant is not in the set of training lights). For example, in color measurement we are often interested in the XYZ tristimuli for a daylight illuminant D65 which has a prescribed but not easily physically replicable spectrum.



We are going to solve Equation (<ref>) for the filter f using an  alternating least-squares procedure. Notice that the input to the optimisation is an initial filter guess denoted by f^seed. 
Let us consider 3 candidate minimisations corresponding to 3 common scenarios for color measurement each of which can be solved using Algorithm <ref>.

1) Multiple Lights: Here we assume that j indexes over cnt illuminants and k=j (per illuminant we make the target XYZs using the same color signals).  We find a single filter which given per-illuminant  optimal least-square 3× 3 correction  matrices will best fit camera data to the corresponding multi-illuminant XYZs.

2) Multiple measurement lights, single target light: Again j indexes over cnt illuminants. But, the target is a single illuminant, for example CIE D65<cit.>.

3) Single Light. This case is, in effect, the simple restriction of the general case, cnt=1. We have one measurement light and one target light. Like the Luther-condition optimisation, we solve for a single correction matrix.
 
In Algorithm 2,  it is straightforward to solve for the jth color correction matrix at iteration step i, M^i_j, in step <ref> using the Moore-Penrose inverse. Step <ref>, where we find the filter, can also be solved directly using simple least-squares, although the basic equations need to be rearranged. Details of the least-squares computation are given in the Appendix. Here, to ensure that the filter is all positive we can also solve for the filter by solving the optimisation subject to the positivity constraints, we solve a quadratic programming problem <cit.> (unlike the Luther-condition case there is no a prior physical reason why the best filter should be all positive). 

Quadratic programming allows linear least-squares problem subject to linear constraints to be solved rapidly and, crucially, a global optimum is found.



 §.§ Adding Filter Constraints
 

By default the filter found using Algorithm 2 can be arbitrarily non-smooth and might also be very non-transmissive. Non-smoothness limits manufacturability (at least with dye based filters) and a filter that absorbs most of the incoming light would, perforce, have limited practical utility. Both these problems can be solved by placing constraints on the filter optimisation.

Let us now constrain the target filter f according to:

    f=Bc,  s.t.  f_min≤f≤ f_max


here B denotes a 31× m basis matrix with each column representing a basis and the vector c denotes an m-component coefficient vector. The scalars f_min and f_max denote lower and upper thresholds on the desired transmittance  of the filter; specifically, f_max is set to 1 by default as fully transmissive and f_min is a positive value between 0 and 1. Equation (<ref>) forces the optimised filter  to be in the span of the column vectors of  B and to meet the min and max transmittance constraints.
By judicious choice of the basis we can effectively bound the smoothness of the filter. For example, we could choose the first m terms of the discrete cosine transform basis expansion<cit.>. 


With respect to the new filter representation, we rewrite the new overall  filter design optimisation in Equation (<ref>) as

    min_c, M_j∑_j=1^cnt∥ C^T_j diag(Bc)QM_j - C^T_k X∥_F^2,    f_min≤  Bc≤ f_max


The current minimisation can be solved using the same alternating least-squares paradigm of Algorithm 2. But, in step <ref>, we need to substitute  the constraint f^i>0 with

    f_min≤ diag(∏_s=0^i-1f^s)f^i = Bc^i ≤ f_max

That is, the filter we find at the ith iteration when multiplied by all the filters from the previous iterations is constrained to be in the basis B. 

Looking at Equation (<ref>) we see that 

    f^i = [diag(∏_s=0^i-1f^s)]^-1Bc^i


That is, effectively the basis for the ith filter changes at each iteration. Again we can solve step <ref> subject to the constraints of Equation (<ref>) using Quadratic Programming. 

Finally, we note that we could, of course, rewrite Algorithm 2 so that at each iteration we solve for a filter defined by a coefficient weight directly (we could solve for an m-term  coefficient vector rather than a 31-component filter). The two formalisms are equivalent. Here, we chose to solve for the per iteration filter for notational convenience: we can use the general Data-driven algorithm and simply change how we calculate the per iteration filter optimisation.




 §.§ Initialising the  Data-driven Optimisation
 

Alternating least-squares is guaranteed to converge but it will not necessarily return the  global optimal result. But, it is deterministic. So, given the state of the correction matrices and filter at the ith iteration we will ultimately arrive at the same solution. Equally, if we change the initialisation condition, f^seed, in Algorithm 2, we may end up solving for a different filter.  Empirically, we observed that the filter returned by the algorithm depends strongly on the initial filter that seeds the optimisation.

Let us consider 3 different ways to seed the Data-driven optimisation:

1) Default:  f^seed= 1. This uniform unit vector denotes a fully transmissive filter over the spectrum. 

2) Luther filter: f^seed=f^Luther. That is we seed the Data-driven optimisation with the optimal Luther-condition filter found using Algorithm <ref>. 

3) By sampling: Here we find a set of sample filters F (which meet our smoothness and transmittance boundedness constraints) and for each filter sample, f∈ F^seed, we will run Algorithm 2.

Algorithm 3  generates #filters (number of initial filters)—subject to bounded smoothness constraints—by uniformly and randomly sampling the filter coefficient space. Before sampling, the algorithm first finds the min and max values of the coefficients (which are calculated in each of the m dimensions individually).  Explicitly, for the ith component in vector c, we denote its minimum and maximum values as c_i^min and c_i^max.

In Algorithm 3, for the minimum value of the ith coefficient, we write: 
c^min_i = minc_i , s.t. 	f_min≤ Bc≤ f_max, i=1,2,...,m. That is, over all possible vectors c, which satisfy the transmittance constraints, we take note of the minimum value of the ith component. That is we find the minimum value that c_i can be over the set of all possible solutions. The maximum of the ith coefficient term is written similarly (see step 3, Algorithm 3).
The minimum and maximum values of c_i can be solved using linear programming<cit.>.


All m min and max components, taken together, make  the two  vectors  c^min and c^max. These vectors together define the extremal values in each dimension of an m-dimensional hypercube. A vector that lies outside the hypercube is guaranteed not to satisfy the boundedness and transmittance constraints we have placed on our filters. This hypercube usefully delimits our search space (of the sample set of solutions).






To generate a set of filters for initialising the optimisation (solved in Algorithm 2) we will sample uniformly and randomly this hypercube.  We use the notation  c_i ∼U(c^min_i,c^max_i  )  to denote sampling a number in the interval [c^min_i,c^max_i] uniformly and randomly.  A filter constructed from the corresponding vector c=[c_1 c_2 ⋯ c_m]^T (f=Bc), will be added into the initial filter set F only if it lies within the transmittance bounds and it is sufficiently far from those filters already in the set, see step 7. In algorithm 3 sufficiently far means at least θ degrees from the other set members. The function cardinality() returns the number of members in a set.




§ RESULTS




 §.§ Experiments for Luther-condition Optimised Filters
 

Let us return to the example shown in Fig. 1.  The optimal Luther-condition filter solved using Algorithm <ref> is shown in Fig. 1e. We multiply the camera sensors by this filter and then find the linear least-squares transform mapping the new effective sensitivities to the XYZ matching functions. The fitted filtered camera sensitivities (to the XYZ color matching functions) are shown  Fig. 1f. In contrast, Fig. 1d, shows the native camera spectral sensitivities fitted to the XYZs. Visually, the addition of our derived filter makes the camera much more colorimetric. 




The reader will notice that the filter, Fig. 1e, absorbs more than 80% of the light  except at the shortest wavelengths where it is maximally transmissive. This need not be a problem for color measurement as we can increase exposure time, for example. Though, it does mean that the camera with and without the filter  would, for the same recording conditions, capture significantly less light. And, this could result in an increase in the  noise in the final image output by a camera reproduction pipeline. Indeed, if we deploy this filter - and keep the capture conditions the same - we would need to `scale up' the recorded values and this operation also scales up the noise. Effectively, we capture an image at a higher ISO number.



 §.§ Vora-Values


To quantitatively measure the spectral match between the filtered and linear transformed camera sensors and the XYZ color matching functions, we calculate the Vora-Value<cit.>.
The Vora-Value measures the  closeness between the spaces spanned by a set of filter sensitivities set Q  and that by the color matching functions X. It is defined as 

    ν( X,Q) =
    	Trace( QQ^+ XX^+)/3

where Trace() returns the sum of the diagonal of a matrix and ^+ is the Moore-Penrose inverse (see Appendix). The Vora-Value is a number between 0 and 1 where 1 indicates the two sensitivity sets span the same space, i.e. the Luther condition is fully satisfied. While there is not a guide on what different Vora-Values mean, empirically we have found when the Vora-Value is respectively larger than 0.95 and 0.99 then we have acceptable and very good color measurement performance. An explanation of why the Vora-Value is useful for quantifying the color measurement potential of a set of color filters together with its derivation  can be found in <cit.>. 



The  Vora-Value performance for a set of 28  camera spectral sensitivities <cit.>  - with and without  their optimised filters - are shown in Fig. <ref>.  The Vora-Value for the unfiltered, NATive sensitivities are shown in blue and for the LUTHer-condition optimised sensitivities in red. On average, the native Vora-Value  is 0.918 but when the optimised filter is added it increases to 0.961. This digital cameras data set <cit.> comprises of diverse camera types, including professional DSLRs, point-and-shoot, industrial and mobile cameras. 




Note, we make a distinction between using a camera for color measurement and for making attractive looking images. Here, we are interested in using a camera to measure XYZ tristimuli (or measures like CIE Lab values that are derived from tristimuli <cit.>). From a measurement perspective, we need higher tolerances and a higher Vora-Value. Clearly, from the point of view of making attractive images cameras that have Vora-Values less than 0.95 can work very well. Indeed, many of the 28 cameras with Vora-Values less than 0.95 can still take images that look appealing from a photographic perspective. But, commercial cameras are not suitable vehicles for accurate color measurement. Quantitative color measurement results are presented in the next section. 



 §.§ Color Measurement Experiments


Now let us evaluate the derived Luther-condition and Data-driven filters in terms of a perceptually relevant color measurement/image reproduction metric. The CIELAB Δ E_ab^* <cit.> - the Euclidean distance calculated between two tristimuli - is  a single number that roughly correlates with  the perceived perceptual distance  between two colors. One  Δ E_ab^* corresponds approximately to a `Just Noticeable Difference' to a human observer <cit.>. When Δ E_ab^* is less than 1 we say that the difference is imperceptible to the visual system. 





  §.§.§ Single Light Case


Let us use the  Canon 5D Mark II camera spectral sensitivities as a putative measurement device and quantify how well it can measure colors - with and without a filter.
In this experiment  the measurement light is either
 a CIE D65 (bluish) or a  CIE A (yellowish) illuminant.   For reflectances we use the SFU set of 1995 spectra <cit.> (itself a composite of many widely used sets).
 The 1995 XYZs for these reflectance and lights are the ground-truth with respect to which we measure color measurement error. 
 
 Using the Canon camera sensitivities, the reflectance spectra and either the CIE D65 and A lights we numerically calculate two sets of 1995 RGBs. 
 Now, we linearly regress the  RGBs for each light to their corresponding ground-truths (we map the native RGBs for CIE D65 and A to respectively the XYZs under the same lights). We call these color corrected RGBs the NATive camera predictions (and we adopt this notation in the results shown in Table <ref>). Rows 1 and 4 of Table <ref> report the mean, median, 90, 95 and 99 percentile and the maximum CIELAB Δ E  error for CIE D65 and CIE A lights.
 
 Now let us  place a filter in front of the camera. 
 Again, we calculate two sets of RGBs (one for each light) for the camera spectral sensitivities multiplied  by the filter found using the Luther-condition optimisation (Algorithm 1). The recorded filtered RGBs are mapped  best  to corresponding XYZs using linear regression. The LUTHer Δ E_ab^*  error statistics are shown in rows 2 and 5. It is clear that placing a Luther-condition Filter substantially increases the ability of the camera to measure colors accurately. Across all metrics the Δ E_ab^*  errors reported are about a third of those found when a filter is not used.
 
 We repeat the experiment for a Data-driven color filter (found using Algorithm 2 where the seed filter for the optimisation is the Luther-condition optimised filter).  Again, the two sets of filtered RGBs are linearly mapped to corresponding XYZs to minimise a least-squares error. Results for the corrected DATA-driven filtered RGBs for the two lights are reported in rows 3 and 6 of Table I. Clearly, the camera plus filter can measure colors more accurately compared to the case where a filter is not used. Across all metrics the Δ E_ab^* errors reported are about a quarter of those found when a filter is not used.
 
 Significantly,  the best Data-driven filter also delivers improved performance compared with the results reported for the  Luther-condition optimised filter. The errors are further reduced by about a quarter. Incorporating knowledge of typical lights and surfaces into the optimisation leads to improved color measurement performance.
 
 
 

  §.§.§ Multiple Lights Case


We now repeat the experiment for a  set of 102 measured lights<cit.>. The results of this second experiment are  shown in the last 3 rows of Table <ref>. Here, the reported error statistics are averages. For each illuminant - as described in the single light case above -  we calculate the mean, median, 90 percentile, 95 percentile, 99 percentile and maximum Δ E_ab^*. That is, we calculate 6 error measures for 102 lights. We then take the mean of each error statistic over all the lights. The aggregate illuminant set performance is reported in rows 7, 8 and 9 of Table I for respectively unfiltered RGBs and RGBs measured with respect to Luther-condition and Data-driven optimised filters.

In terms of the reported errors of the raw RGBs compared to the filtered RGBs for the Luther- and Data-driven filters we see the same data trend for the multiple lights case as we saw previously for single lights. A Luther-condition derived filter reduces the measurement error by 2/3 and for the Data-driven filter the measurement error is reduced by 3/4, on average.



  §.§.§ Multiple Cameras


Now, we calculate the mean Δ E_ab^* error (for the 102 lights and 1995 reflectances) for each of 28 cameras <cit.>. For each camera we calculate the optimal Luther-condition and Data-driven optimal filters (where as before the Luther-condition filter seeds the Data-driven optimisation). Per camera, Figure <ref>
summarises the per camera mean and 95 percentile Δ E_ab^* performance.

Grey bars  in Figs. 4a and 4b show respectively the mean and 95 percentile error performance of native (un-filtered) color corrected RGBs for the 28 cameras. Respectively, the dashed green and solid red lines record the performance of the best Luther-condition and  Data-driven filters. 

It is evident that the optimised filters support improved color measurement for all 28 cameras and on average the performance increment is significant. For many cameras the Data-driven optimised filter delivers significantly better color measurement performance compared with using the  Luther-condition optimised filter. 
  







 §.§ Smooth and Bounded Transmittance Filters
 




Both the Luther-condition and Data-driven filters 
absorb much of the available light (low transmittance values) and are far from being smooth, e.g. see the derived Luther-condition filter in Fig. 1e. Here across much of the visible spectrum the filter transmits little (below 20%) of the available light. 
When a strongly light-absorptive filter is placed in front of a camera then we need to either increase the   exposure time (or widen the  aperture) or apply a  higher ISO (which can increase the conspicuity of noise) to obtain an image with the same range of brightnesses (as when a filter is not used). Plus the filter in Fig. 1e is not smooth it may be difficult to manufacture.

In this subsection we will constrain our optimisation so that the calculated filters are smooth and transmit, per wavelength, a minimum amount of incident light (here, 20%). We enforce smoothness indirectly by assuming that our filters lie within the span of either a 6- or 8-dimensional Cosine basis.


Let us visualise the 20%  bounded transmittance constraint using the  Canon 5D Mark II camera sensitivities and our Data-driven optimisation. First, we initialise the optimisation (Algorithm 2) using  the uniform vector 1 (100% transmittance across the spectrum) as the seed  filter. The optimisation returns the filter, denoted  DATA_1s and is plotted in the blue-dotted line in Fig. 5a. The Data-driven optimisation seeded with the Luther-condition filter DATA_Luther  is also shown in Fig. 5a (dotted purple line).

We repeat this experiment where we both require the derived filters to transmit at least 20% of the light and also that they belong to the span  the first 6- or first 8-terms in a Cosine basis expansion. When the recovered filter is constrained to lie in the span of a 6-dimensional Cosine basis, the recovered DATA_1s and DATA_Luther - for the two initialisation conditions -  filters are shown in Fig. 5b (respectively blue and purple dotted lines). See Fig. 5c for the filters calculated using an 8-dimensional Cosine basis. 

The red lines shown in Figs. 5b and 5c are the filters optimised by our sampling Algorithm (using Algorithm 3 to find a set of filters to seed algorithm 2 and then choosing the best one that has the best overall performance). The experiment for deriving these filters are described in the next subsection.

By examining Figs. 5a, 5b and 5c it is evident that - when no basis, a 6- or an 8-dimensional Cosine basis are used - that changing the initialisation condition results in a different filter being found.



 §.§ Sampled Optimisation


Using Algorithm 3, let us run a sampled optimisation. That is for a given Cosine basis we calculate a set of candidate solutions, the sample set F. Here
we populate F with  20,000 uniformly and randomly generated filters where the angular threshold between any two filters in the set is 1 degree, θ = 1, see step 7 in Algorithm 3. Each filter in F transmits at least 20% of the light (and F
is populated by filters described as linear combination of a 6- or 8-dimensional Cosine basis).

Each filter in F is used to initialise the Data-driven algorithm. That is we find 20,000 optimised filters. The color measurement performance of each filter in this solution set can be calculated. Then we simply choose the filter that delivers the best overall measurement performance. In Figs. 5b and 5c we show the best  sample-optimised filters (red lines)  which respectively lie in the span of a 6- and 8-dimensional Cosine basis (and transmit at least 20% of the light). Here 'best' is defined to be the filter that results in the smallest mean Δ E^*_ab performance.

Table <ref> reports the Δ E_ab^* color error performance 1995 reflectances and 102 lights <cit.> for the Canon 5D Mark II sensitivities. The row NAT reports the baseline color correction results when a per illuminant based linear correction matrix is applied while no filter is used (note that row 7 in Tables I and row  1 in Table II are the same).

In Table II we report the correction performances in 3 tranches. Rows 2 and 3 correspond to the two filters without using Cosine basis as shown in Fig. 5a. Here we find the best filters with only the 20% minimum transmittance bound. 
Rows, 4,5 and 6 report the performance when the 3 filters shown in Fig. 5b are used, where the filter is additionally constrained to be in the span of the 6-dimensional Cosine basis. 
Finally, when the filter is constrained to belong to an 8-dimensional Cosine basis, the 3 derived filters lead to the error statistics shown in rows 7 through 9.

Table I reported the color measurement performance of the filters found using an unconstrained optimisation. Table II reports the color measurement results that are found when filters are constrained to have a bounded transmittance (here at least 20% of the light) and be smooth.
Let us consider the bounded transmittance first. 
Comparing row 9 of Table I to row 3 of Table II we see that adding a lower transmittance bound returns a filter that delivers poorer measurement performance (but still much better compared with the native camera response). 
Additionally, requiring that our filters smooth on top of the minimum also yields relatively poorer performance compared to the unconstrained filter optimisation.

However, with either the 6- and 8-dimensional Cosine basis constraint we can find the best filter by seeding Algorithm 2 with many possible filter initialisations (and then choosing the best filter overall). Here, we find that comparable performance is possible. Compare rows 6 and 9 of Table II to row 9 of Table I. It is remarkable how well a constrained filter can work: the performance is ever so slightly worse than the unconstrained optimisation. But, the filter is much smoother and more likely to be able to be manufactured.






 §.§ Sampling vs Optimisation


It is worth reflecting on our sample-based optimisation. Clearly, that sampling makes such a difference to the performance that our optimisation can deliver (for filtered color measurement) teaches us that the minimisation at hand has many local minima. By sampling we are effectively allowing our minimiser (Algorithm 2) to find many solutions and then we have the latitude to choose the (closer to) global minimum. Given we seed our optimisation with 20,000 filters we might wonder whether we
 need to actually carry out the Data-driven optimisation.
 
In answering this question, first we remark that it is well known that as the dimension of a space increases it is more sparse. On the Cartesian plane if we have more than 360 vectors (anchored at the origin) then the closest angular distance to at least one vector's nearest neighbours must be less than 1 degree. In 3-dimensions we can have thousands of vectors where every vector is more than 1 degree from its nearest neighbour. 

For our 20,000 member sample  set F we calculated the average angular distance for each element to its nearest neighbour in the set. When F is calculated subject to the 6-dimensional Cosine basis constraint, the average nearest-neighbour distance was 2.6 degrees (with a maximum of 7) and for the 8-dimensional case it was 4.6 degrees (with a maximum of 10). Running the optimisation, Algorithm 2, with each element in F, we effectively refine the initial guess. And, the refinement (difference between the starting and endpoint filter) is on the same order as the average nearest-neighbour distance. 

Significantly, running the optimisation - carrying out the refinement - results in a significant performance increment compared to using the only sample filters. That is we cannot use the sampling strategy alone to find the best optimised filter. The importance of the refinement step will increase as a greater number of basis functions are used in the optimisation.




§ CONCLUSION

In this paper, we developed two algorithms that design transmittance filters which, when placed in front of a camera, make the camera more colorimetric. Our first algorithm is driven by the camera sensitivities themselves. It is well known that a camera that has sensitivities that are a linear transform from the XYZ color matching functions - the camera meets the Luther condition - can be used to measure color without error. Our first algorithm finds the filter that best satisfies the Luther condition. A second algorithm that tackles color correction for a given set of measured lights and surfaces, which we call Data-driven filter optimisation, is also developed. Both Luther- and Data-driven Filters provide a step change in how well a camera can measure color.

Our default optimisation - though compellingly simple to formulate - deliver filters which are not smooth (difficult to make) and may also transmit very little light. Our optimisations are reformulated to incorporate both smoothness and a lower bound on how much light must be transmitted. Initially, when these constraints are adopted, the solved-for filters work less well. However, experiments demonstrated that our optimisations were highly dependent on the initialisation parameters, specifically the seed filter (initial guess) that drives the filter evolution. A simple sampling strategy - i.e. severally running the optimisation for a set of judiciously chosen seed filters - allows us to mitigate this problem. Significantly a smooth filter that transmits more than 20% of the light across the visible spectrum delivers almost as good performance as a very non-transmittive and non-smooth filter (found via the unconstrained optimisations).

[Implementation]

For  both Algorithms 1 and 2 presented in Section <ref> the filter and the color correction matrices can be found using simple least-squares regression. To remind the reader, given  A and B - m× n matrices of rank n where m ≥ n, then the least-squares regression M - an n× n matrix, mapping A to B (AM≈ B)
can be found in closed-form using the Moore-Penrose inverse <cit.>:

    M=[A^TA]^-1A^TB=A^+B




 §.§ Algorithm 1: ALS for the Luther-condition Optimisation


In step 4 of the algorithm, the optimal filter is found by finding scalars that maps each row of Q^i-1 to the corresponding row of X. The best scalar α mapping the vector  v=[Q^i-1_j]^T to w=[ X_j]^T (for the jth row of the data matrices) can be written in closed form using the Moore-Penrose inverse: 
α=v^Tw/v^Tv. 

Similarly, in step 5, the Moore-Penrose inverse can be used for finding M. Denoting Q=diag(f^i)Q^i-1 then M^i= Q^+ X= 
[ Q^T Q]^-1 Q^T X.



 §.§ Algorithm 2: ALS for the Data-driven Optimisation


In step <ref>, each M^i_j can be solved directly using the Moore-Penrose inverse. Denoting Q=C_j^T Q^i-1_j then M^i_j= Q^+C_k^T X=[ Q^T Q]^-1 Q^TC^T_k X.

In step <ref> of the Data-driven optimisation, the filter f is embedded in the equation and so we cannot solve for it directly as we could for the Luther-condition case.


To solve for the filter it is useful to vectorise the minimisation. We recapitulate the minimisation statement of step <ref>:

    min_f∑_j=1^cnt∥(C^T_j diag(f) Q_jM_j - C^T_k  X)∥_F^2


The meaning of ∥ ∥^2_F is the Frobenius norm squared, i.e. the sum of all the argument terms squared. This Frobenius norm is generally applied to matrices (as here) but can equally be applied to vectors (where the vec() operator stacks the columns of a matrix on top of each other):

    min_f∑_j=1^cnt∥vec(C^T_j diag(f) Q_jM_j) - vec(C^T_k  X)∥_F^2


Now let us rewrite the diagonal filter matrix as a summation of each value in the diagonal, f_i, multiplied with a single entry matrix D_i as diag(f) = ∑_i = 1^31 f_i D_i. Here, D_i is a 31 × 31 matrix with a single non-zero entry at D(i,i) = 1. 
By substituting this new filter representation into the first term of the minimisation in Equation (<ref>), we obtain

    vec(C^T_jdiag(f)Q_jM_j) = ∑_i=1^31f_i   vec( C^T_jD_iQ_jM_j)


Now let us denote a matrix V_j = [v_1,j v_2,j ⋯ v_31,j] where its column represents v_i,j=vec(C^T_jD_iQ_jM_j), Equation (<ref>) can be expressed more compactly as

    vec(C^T_jdiag(f)Q_jM_j) = V_jf

where f = [f_1   f_2  ...  f_31]^T. Note that if  C_j is a 31 × n matrix then C^T_jD_iQ_jM_j is an n × 3 matrix (where 3 denotes the number of color channels) and thus v_i,j is 3n × 1 which makes matrix V_j have size of 3n × 31. 

Now we stack all V_j, j = 1,2,..,cnt  matrices  (under cnt different lighting conditions) on top of each other making an (3n*cnt)× 31 matrix, V. Similarly we stack all cnt targeting XYZs on top of each other denoted as w = vec(C_k^T X) which has the size of (3n*cnt)× 1. We remind the reader that k might equal j. Or, k might denote a single privileged illuminant such as CIE D65. 

Now the minimisation in Equation (<ref>) can be equivalently rewritten as:


    min_f∥ Vf - w∥^2_F

The best f can be found in closed form using the Moore-Penrose inverse:

    f= (V)^+w=[V^TV]^-1V^Tw.





 §.§ Filter Constraints

Equation (<ref>) solves for the 31-component f in one step. Suppose we write  f=Bc. We constrain the filter to be describable by a linear basis (B is 31× m where 1≤ m ≤ 31).  Additionally, the filter is restrained by a minimum f_min and maximum f_max bounds on the transmittance.
Then to solve for the filter we find the coefficient vector c that minimises:

    min_c∥ VBc - w∥^2_F
         s.t.  f_min≤Bc≤ f_max

Equation (<ref>) where there is a quadratic objective function and linear inequality constraints can be solved using quadratic programming <cit.>. 




§ ACKNOWLEDGMENT

This work was supported by EPSRC under Grant EP/S028730. The authors would also like to thank Dr. Javier Vazquez-Corral for his insightful comments. 



IEEEtran



