An accurate method for computing correlated color temperature

: For the correlated color temperature (CCT) of a light source to be estimated, a nonlinear optimization problem must be solved. In all previous methods available to compute CCT, the objective function has only been approximated, and their predictions have achieved limited accuracy. For example, different unacceptable CCT values have been predicted for light sources located on the same isotemperature line. In this paper, we propose to compute CCT using the Newton method, which requires the first and second derivatives of the objective function. Following the current recommendation by the International Commission on Illumination (CIE) for the computation of tristimulus values (summations at 1 nm steps from 360 nm to 830 nm), the objective function and its first and second derivatives are explicitly given and used in our computations. Comprehensive tests demonstrate that the proposed method, together with an initial estimation of CCT using Robertson ’s method [J. Opt. Soc. Am. 58, 1528-1535 (1968)], gives highly accurate predictions below 0.0012 K for light sources with CCTs ranging from 500 K to 10 6 K.


Introduction
The spectral power distribution (SPD) of a Planckian radiator [1], denoted by M(, T), can be completely determined from its absolute temperature T in Kelvin (K)using the expression: (1) where constant c1 = 3.741832×10 -16 Wm 2 and c2 = 1.4388×10 -2 m K. The color temperature T of a light source is the color temperature of a Planckian radiator which emits radiation of the same chromaticity as the light source. When a light source has a chromaticity which is not the same as the one of the Planckian radiator, the term "correlated color temperature" (CCT) is used, defined as the temperature of the Planckian radiator whose chromaticity is the one nearest to that of the light source in a uniform chromaticity scale (UCS) diagram [1][2][3]. Specifically, the UCS diagram recommended by the International Commission on Illumination (CIE) in 1971 for CCT computation is the CIE 1960 UCS diagram [1,2], where the chromaticity coordinates are given by : 4 / ( 2 12 3), 6 / ( 2 12 3), and x, y are the chromaticity coordinates of the light source in the CIE 1931 coordinate system [4,5] defined in terms of its tristimulus values (TSVs) X, Y, and Z as follows: / ( ), / ( ).
x X X Y Z y Y X Y Z Since it is a single number, CCT is simpler to communicate than is the SPD, which led the lighting industry to accept CCT as a shorthand means of reporting the color appearance of "white" light emitted from electric light sources. CCT values are intended by the lighting industry to give a general indication of the apparent "warmth" or "coolness" of the light emitted by a light source. Also, scientists [6][7][8][9] found that CCT can be used to determine the relative SPD of daylight illuminants. Furthermore, CCT is also used to evaluate the colorrendering index (CRI) [10]. For many lighting applications, in offices, hotels, galleries, the textile industry, etc., the CCT and CRI are the two most important factors governing the choice of appropriate light sources.
The computation of the CCT based on the given SPD of a light source/illuminant has a long history. Let u, v be the chromaticity coordinates of the light source/illuminant, and u(T), v(T) be the chromaticity coordinates of the Planckian radiator with color temperature T. Then, the square of the distance between (u, v) and (u(T), v(T)) is the function f(T): It follows from Eqs. (1)(2)(3)(4)(5) that the computation of the CCT is a problem of non-linear optimization. Therefore, there is no analytical expression for computing CCT, and scientists and engineers have long been searching for simple and practical methods to compute CCT in industrial applications.
In the early stage, the so-called isotemperature lines were used to estimate the CCT from graphs. Judd [11] appears to have been the first gave a table of CIE x, y chromaticity coordinates of the intersections of a series of lines of constant CCT (isotemperature lines) with the Planckian radiator locus. The separation between isotemperature lines was represented in terms of the unit called "micro-reciprocal degree" or "reciprocal mega-Kelvin" [8, page 224] ("rd" or "mired") Mrd ( = 10 6 /T) rather than the correlated color temperature T, in order to achieve greater uniformity in plots in the CIE 1960 UCS diagram.
The isotemperature lines were by definition perpendicular to the Planckian radiator locus when plotted in the CIE 1960 UCS diagram. In 1963, Kelly [12] published graphs with isotemperature lines in the x, y, and u, v chromaticity diagrams to estimate the CCT. In 1968, Robertson [13] gave an explicit interpolation formula for computing CCT based on the isotemperature lines. The speed and accuracy of Robertson's method depends on how many isotemperature lines are chosen. That is, the more lines used, the greater the accuracy, but also the more time spent for computation.
Another type of estimation of the CCT of a light source is based on the explicit approximation of the chromaticity coordinates x(T), y(T) or u(T), v(T) of the Planckian radiator, for less computational cost. Then the bisection method [14] can be used for minimizing the objective function defined by Eq. (4). In 1977-78, Schanda and Dányi [15] and Schanda et al. [16] proposed an approximation for x(T), y(T) by using up to a 7 th -order polynomial of the micro reciprocal degree for color temperatures in the range from 1700 K to 50000 K. In 1984, Krystek [17] offered an approximation for u(T) and v(T) using a rational Chebyshev approximation for the color temperature in the range from 1000 K to 15000 K.
In 2000, Gardner [18] proposed the use of the Newton method [14] to compute the CCT from the chromaticity coordinates u, v of a light source. The Newton method is an iterative method for finding the roots of a differentiable function f'(T) (i.e. solutions to the equation f'(T) = 0 in our case), by the use of the equation where f'(Tj) and f''(Tj) are the first and second derivatives of the function f(T) defined by Eq.
(4). It should be noted that Gardner [18] did not use the exact first and second derivatives, but rather approximations based on fitted 7 th -order polynomial relationships between v(T) and u(T), and between the inverse color temperature (1/T) and u(T). Thus, Gardner's iterative method depends on the degree of approximation of the formulae used, and cannot provide an exact solution, as intended in the current paper. A third type of method to estimate the CCT is based on established functional relationships between the CCT and the chromaticity coordinates x, y or u, v of the light source, as intended, for example, in the method proposed in 1987 by Qiu [19]. Qiu's formula for estimating the CCT, Te, can be expressed by: For the exact function q, please consult the original paper [19]. Qiu's method was recommended for CCTs ranging from 2500 K to 10000 K.
Later, in 1992 McCamy [20,21] proposed a simpler formula for estimating Te from x, y chromaticity coordinates, for CCTs ranging from 2000 K to 12500 K: 3 2 , e T a n bn cn d     n x x y y x y      (9) In 1999, Hernández-Andrés et al. [22] developed another explicit formula valid for estimating a wider range of CCTs up to 10 6 K. This formula is given by the equation: where n is defined as in previous Eq. (9), but using two different epicenter coordinates xe and ye within two different ranges of CCTs, and the values of the coefficients Ai (i=0,1,2,3) and tj (j=1,2,3) can be found in the original paper [22].  Table 1 shows the CCT ranges recommended by the different computational methods proposed in the literature. As we can see, most methods cover the range only from 2000 K to 10000 K, which is the most important one in industrial applications.
While CCT is a classical topic in color science, it can be said that the ideal method to compute CCT currently remains an open question. Thus, in draft #11 of future CIE Publication 15, CIE TC 1-85 [23] has added a new note stating that "computed correlated colour temperature is highly dependent on specifications of the computing process…. An agreement about all details of the computing method is necessary to avoid undesirable discrepancies between computed correlated colour temperatures." As will be shown in next sections, the method proposed in the current paper can be useful to achieve accurate computed values of CCT in comparison with methods previously reported in the literature.

Predictions of methods currently available for CCT computation
We tested the performance of all the methods mentioned in Table 1 by using sets of samples (hypothetical light sources) with known CCTs. Specifically, we considered sets of 5 samples placed on different isotemperature lines, which by definition must have the same CCT for a perfect computational method. For any selected color temperature T from 500 K to 25000 K at 1 K steps, the u(T) and v(T) of the Planckian radiator were computed from TSVs (see Eqs. (5), (3), (2) and summations from 360 nm to 830 nm at 1 nm steps), and the associated isotemperature lines were determined. The slopes of the isotemperature lines were determined from u'(T)/v'(T) [13] (with u'(T) and v'(T) defined by Eq. (15) below) and linear equations of the isotemperature lines were also defined by the slopes. Then, we computed the next five equi-spaced samples (Sample #1-#5) in each isotemperature line using the linear equation: Sample #3 was placed exactly on the Planckian radiator locus in the u, v space; Samples #1 and #2 were placed above the Planckian locus at distances of 5×10 -2 and 2.5×10 -2 units in u, v space, respectively; Samples #4 and #5 were placed below the Planckian locus at distances of 2.5×10 -2 and 5×10 -2 units in u, v space, respectively . The distance of 5×10 -2 units to the Planckian locus in the u, v space was chosen because it is the limit currently established by CIE [9,23] in order to apply the concept of CCT. The maximum absolute difference between the predicted and exact CCTs of the five aforementioned samples divided by the exact value of the CCT of the isotemperature line will be designated as |T/T|MAX and it can be considered to be the maximum relative prediction error for the isotemperature line. Figure 1 shows the maximum relative prediction errors (|T/T|MAX) for each of the seven methods indicated in Table 1 Figure 1 shows that no method gives perfect predictions over the entire CCT range tested. In our computations for the JS method, the AR method was used for the initial CCT estimation, considering 30 intervals (i.e. 31 isotemperature lines) between 500 K and 25000 K. From Fig. 1, the three best methods appear to be JS, AR, and XQ (in that order). Bearing in mind that the predictions of the AR method were the initial guess for the JS method, we expected the JS method to perform better than the AR method.  Table 1 for CCTs ranging from 2000 K to 10000 K.  Figure 1 also shows that some methods can be very good for specific CCTs. For example, the HA method is very accurate for CCTs around 7400 K. This depends on the training data employed for deriving the formulae. Figure 2 shows the performance of all seven methods for CCTs ranging from 10000 K to 25000 K (note that the y-axes scale in Figs. 1 and 2 are different). In this CCT range the rank of methods from the best to the worst is: JS, AR, XQ, HA, JM, JG, and MK. Therefore, it can be said that JS, AR, and XQ are the three best methods (in that order) to compute CCTs from 2000 K to 25000 K.  Table 1 for CCTs ranging from 10000 K to 25000 K.  Table 1 for CCTs ranging from 25000 K to 50000 K. Figure 3 shows the performance of the AR, JS, and HA methods for CCTs ranging from 25000 K to 50000 K (for this range of high CCTs the remaining four methods shown above in Figs. 1 and 2 cannot be applied, as indicated in Table 1). For this test, 30 intervals (or 31 isotemperature lines) between 500 K and 50000 K were used for the AR method, which is called AR30. The results from the AR30 method were used as initial estimations of the JS method. For comparison, the results from the AR method with 60 intervals (or 61 isotemperature lines) between 500 K and 50000 K were found as well, which are called AR60. It can be seen that the results from AR60 are better than those from AR30, as might be expected. However, results from AR60 and JS are roughly the same.  Table 1 to predict five test samples ('Sample #1 (S#1)' to 'Sample #5 (S#5)' in the order from top downward in u, v space) placed on the isotemperature line with T=8000 K.
All the tested methods often gave different CCT predictions for the five samples we considered on each isotemperature line, as shown in Fig. 4 for the isotemperature line with T=8000 K (similar results were found for other lines). The heights of the bars in Fig. 4 represent relative prediction errors, |T/T|. Samples #4 and #5 are located below the Planckian locus at distances of 2.5×10 -2 and 5×10 -2 units in u, v space, respectively. Similarly, Samples #2 and #1 are located above the Planckian locus at distances of 2.5×10 -2 and 5×10 -2 units in u, v space, respectively. Thus, in Fig. 4 the same color bars represent the prediction errors by the different methods for the same test sample. In general, for all methods the predictions for the five samples located on an isotemperature line were different. Specifically, the prediction errors increased with the distance to the Planckian radiator locus, and in most cases the maximum prediction errors were found for Sample #5, which is the one farthest below the Planckian radiator locus. On the other hand, as expected, the best (but non-null) prediction errors were found for the Sample #3, located on the Planckian radiator locus. This was expected, because each method uses its own approximation formulae, which depend on the development data.

The proposed method
In this paper, the Newton method defined in Eq. (6) is proposed for computing the CCT. The problem with available methods is that either the function f(T) defined by Eq.(4) itself or its first and second derivatives with the color temperature T have not been accurately provided/used. Here we have used the exact formulae for f(T), f'(T), and f''(T), as well as the current CIE recommendation for computation of tristimulus values (TSVs) using wavelengths from 360 nm to 830 nm at 1 nm steps [9,23]. Coordinates u, v of the test light source (or illuminant) must also be computed from accurate TSVs using CIE recommendation (summations from 360 nm to 830 nm at 1 nm steps). Some methods have been recently proposed [24,25] for accurate TSVs computations in the case available spectral data are not as required by CIE. For the function f(T), the exact u(T) and v(T) are needed, which depend on the TSVs X(T), Y(T) and Z(T) (see Eqs. where 0= 360 nm, n= 830 nm, j+1-j = 1 nm, and the function M is given in Eq. (1). The scaling factor k in Eq. (11) is not important, since it will be canceled when the chromaticity coordinates u(T) and v(T) are computed via x(T) and y(T) (see Eqs. (3) and (2)). Next, from Eq. (4) it follows that: and 2 2 ( ) 2( ( )) 2( ( )) ( ) 2( ( )) 2( Thus, we have and therefore, From Eq. (1) we also have: 2 2 Go to Step 1 else Do not converge within an allowed number of iterations and stop. end if Note that the method proposed above needs an initial guess of the CCT. Bearing in mind the results in the previous sections showing that AR method performs the second best overall; we suggest that the initial guess T0 can be determined using this method.
The MATLAB software to compute the CCT of a light source using any of the methods mentioned in this paper is available from its first author.

Performance of the proposed method
To test the performance of the proposed method, we first generated test samples in the same way as above. For each color temperature T, five samples were generated from the associated isotemperature line as mentioned in Section 2. The color temperature T varied from 500 K to 10 6 K at intervals of 1 K, meaning that we used 4,997,505 samples (999,501 isotemperature lines × 5 samples per line). The testing was run on a PC (CPU: Intel i5-4430 with 3.0 GHz, RAM: 8.0 GB and Windows 10: 64 bits) using MATLAB software. The proposed method is compared with Roberson's method with 30/60 intervals (or 31/61 isotemperature lines), and HA method, respectively, since they are the only two methods which are also valid for this wide CCT range (see Table 1). The results are listed in Table 2, where, for samples on each isotemperature line, we provide the maximum prediction error |∆T|MAX and the maximum relative prediction error |∆T/T|MAX for each method. In the last column of Table 2 the CPU time in seconds per sample is given for each method. It should be remembered that the initial guesses for the method proposed in the current paper were provided by the AR30 method, and the time employed to perform such initial guesses have been included in the values shown in the first row of Table 2.  Table 2, we see that the maximum of |∆T|MAX for the proposed method was 0.0012 K, meaning almost no errors at all. In fact, for CCT less than 5000 K, the maximum of |∆T|MAX for the proposed method was found to be less than 1.5×10 7 K. However, when the CCT increased, the maximum of |∆T|MAX also increased, as reflected in Fig. 5. Note that for each CCT in horizontal axis of Fig. 5, there is only one corresponding value |∆T|MAX in the vertical axis. It looks there are many points in a particular vertical line because there are too many vertical lines overlapping (at 1 K steps). Similar performance phenomena were found for |∆T/T|MAX. While for the AR method with 31 and 61 isotemperature lines, the maximum values of |∆T|MAX were 72570 K and 40739 K, respectively, for the HA method, the maximum value of |∆T|MAX was 814590 K. Thus, we can conclude that the method proposed in the current paper is the best one available, providing almost exact CCT values. In any case, it should be considered that our proposed method provides highly satisfactory results because it is based on an initial guess based on the AR method, which was one of the best methods available in the literature (see Section 2). However, for the proposed method, the total CPU time per sample was 0.0162 seconds, which is about 10-fold longer than the time used in HA and AR methods.
In the above test, any possible light source/illuminant around Planckian radiator within a distance of ±5×10 -2 units in u, v space ranging from 500 K to 10 6 K was enumerated and tested at 1 K steps. It can be seen that for a real light source, which does not deviate from the Planckian locus far than ±5×10 -2 units in u, v space, the proposed algorithm should work well. Fig. 5. Maximum absolute prediction errors (|T|MAX) for the proposed method for CCTs ranging from 500 K to 10 6 K. The AR method with 31 isotemperature lines provided the initial guess for the proposed method. Figure 6 shows the maximum number of iterations for convergence in the proposed method using a fixed tolerance of 0.001 K. For CCTs between 500 K and 10 6 K, it takes 9 iterations for convergence in the worst case. About 5 iterations are necessary for convergence with CCTs lower than 15000 K, and 7 iterations for CCTs lower than 50000 K, using the aforementioned tolerance limit. It is important to note that the method proposed in this paper uses the Newton method, which converges locally. This means that the convergence of the proposed method needs a good initial guess. It was found that the proposed method works well for CCTs ranging from 500 K to 10 6 K with the initial guess provided by the AR method. It was also found that, using CCT Maximum number of iterations the initial guess provided by the XQ method, the proposed method works well for CCTs less than 25000 K, but it may not converge for CCTs larger than 25000 K.
It is also worth mentioning that CIE [9,23] recommended the 1 nm summation should be used for computing TSVs in Eq. (11). That is why the 1 nm summations were used here for anything related to tristimulus integration. If sampling rate is changed from 1 nm intervals to 5 nm, 10 nm or 20 nm intervals, the proposed method may converge faster, but using sampling rates larger than 1 nm intervals in Eqs. (11), (21) and (22), the functions in Eqs. (4), (6), (12), and (13) are only approximated rather than exact. Hence, the accuracy of the method proposed in the current paper may decrease when TSVs are computed using summations at steps larger than 1 nm. Specifically, the proposed method has been tested with CCTs ranging from 2000 K to 10000 K at 1 K steps using sampling rates of 1 nm, 2 nm, 5 nm and 10 nm, and the maximum absolute errors |∆T|MAX found were 1.6×10 -9 , 0.0708, 0.3303, and 22.4599 K, respectively. Perhaps for most practical purposes we may use up to 5 nm sampling rate in the proposed method for CCTs ranging from 2000 K to 10000 K, but anyway, we still recommend using the 1 nm sampling rate for two main reasons: the 1 nm sampling rate is the recommendation made by the CIE [9,23], and it gives most accurate predictions for a wider CCT range.

Conclusions
In this paper, we consider CCT prediction for a given light source or illuminant. In the literature, all available methods give CCT predictions with some lack of accuracy, which increases with the distance of the light source to the Planckian radiator locus. In an initial stage, simple methods were proposed for approximating the CCT, such as plots provided by Kelly's isotemperature lines [12], Robertson's interpolation method [13] based on the isotemperature lines, or explicit formulae from Qiu [19], McCamy [20,21], and Hernández-Andrés et al. [22]. With the increase in computational power, iterative methods were proposed for a better prediction of CCT, such as those suggested by Schanda and Dányi [15], Krystek [17], and Gardner [18]. However, these methods used approximated functions f, or approximated first and second derivatives, to estimated u(T) and v(T), rather than exact u(T) and v(T) functions. In this paper, the Newton method is proposed for computing the CCT. This method finds color temperature T satisfying ( ) 0, fT   (23) using the iterative formula defined in Eq. (6). Examining Eq. (6), we note that it needs the first and second derivatives of the function f. In this paper, the function f, and its first and second derivatives are accurately defined and computed following the current CIE recommendations [9] for computation of TSVs using 1nm summations. Thus, the proposed method gives CCT predictions with an accuracy of below 0.0012 K, so long as the distance between the chromaticity of a test light source and the Planckian radiator spectral locus is not greater than 5.0×10  in the CIE 1960 UCS diagram, as required by the CIE [9,23]. The performance test shows that the proposed method converges with no more than nine iterations within a tolerance of 0.001 K for CCTs ranging from 500 K to 10 6 K. Thus, the proposed method can also be considered to be a direct method.