Sage Journals: Discover world-class research

Abstract

Multi-primary tunable LED lighting systems can generate a vast number of spectral outputs, with the exact number depending on the resolution of the control signal and number of LED primaries. Computing all combinations to identify optimal spectral power distributions (SPDs) with a specified chromaticity would require a tremendous amount of time and computational power. Here, an enumeration reduction (ER) algorithm is described to reduce the computation time by defining a bounding pyramid in the three-dimensional tristimulus space that maps to the circumscribed square of the target region in the chromaticity diagram. This method enables computing only the necessary amount of tristimulus values (and resulting chromaticity coordinates) for SPDs that fall into a user-defined target chromaticity area, reducing computation time by avoiding the need to generate each combination and calculate its chromaticity coordinates. The results show that the proposed ER algorithm can greatly reduce the time – up to 922 times in a series of tests – to determine a set of metamers compared to a full enumeration computation.

1. Introduction

1.1 Multi-primary lighting systems

Solid-state lighting (SSL) devices have become increasingly prevalent due to their efficiency, compact size, robustness, lifetime, spatial and spectral output control, and ability to be dimmed. SSL devices, such as LEDs, utilize additive colour mixing principles to generate white light, with energy radiated directly from the diode or from phosphors that down convert energy in shorter wavelengths to energy in longer wavelengths (i.e. Stoke’s shift). The most common type of LED system in use today relies on a nominally blue ‘pump’ LED and phosphors that produce nominally green, yellow, orange and red light, which in combination produces white light. These products provide a single, fixed output at full power, such that system performance is straightforward to measure and report. They are cost-effective and can have the highest efficacy of any available light source.

White light can also be generated by mixing the output of multiple LEDs, such as nominally red-, green- and blue-emitting diodes, with each referred to as a primary. Because the output of each primary can be controlled independently, multi-primary LED (mpLED) systems can offer flexibility in spectral output, with a single product capable of delivering a wide range of spectral power distributions (SPDs). This can allow the performance of a lighting system to vary after installation, being optimized over the course of use for colour quality,¹ non-visual stimulation,² consistency of appearance,³ energy efficiency,⁴ daylight and electric lighting balance,⁵ horticultural effectiveness⁶ or other attributes. MpLED systems are also of interest because they offer a higher theoretical maximum luminous efficacy due to their potential elimination of Stoke’s shift losses and ability to deliver a more tailored SPD comprised of specific narrowband components. The U.S. Department of Energy projects mpLED systems to be the most efficacious light source by the 2030s.⁷

In different lighting applications, optimization targets (i.e. photometric, colorimetric and circadian metrics) can vary, and finding the optimal combinations of spectral primaries to meet any given parameter might require an extensive search, resulting in a long processing time. For example, determining the SPDs with chromaticity of 3500 K (threshold radius r_t = 0.0005) among combinations produced from 16 LED primaries at just 20% output increments (0%, 20%, 40%, 60%, 80% and 100%) took 96 core hours (1 wall clock hour on four 24 core nodes) at the Pacific Northwest National Laboratory’s (PNNL) Computational Science Facility in Richland, Washington – PNNL’s Constance cluster is composed of 528 nodes each with 24 core Intel^® Haswell E5-2670 v3 (running at 2.3 GHz) with 64 GB of 2133 MHz ECC memory and connected through fourteen data rate InfiniBand^™ network card. For nominal white light combinations that fell within an area close to Planckian radiation defined by ANSI,⁸ colour rendering metrics, such as ANSI/IES TM-30 fidelity index R_f, gamut index R_g and local chroma shifts R_cs,hi were calculated. Despite the long calculation time, this result is still too coarse for practical use.

The number of possible spectral output states of an mpLED system is equal to the number of control increments raised to the number of primaries. For example, a system with 8-bit control (256 levels) and 5 primaries can generate 256⁵ (1,099,511,627,776) unique conditions. To reduce the computational burden, the search space – defined as all the possible spectral power combinations in an mpLED system – could be limited to a particular region of chromaticity, perhaps even one small enough to contain only visually indistinguishable lights. However, there is currently no method for searching only this limited range, and instead the chromaticity of all combinations must be computed and checked against a criterion. Addressing this problem by utilizing the tristimulus values (TSVs) of the primaries is shown in this work to provide substantial computational time savings. The proposed enumeration reduction (ER) algorithm can help facilitate real-time computations and expedite computing and characterizing the performance of mpLED systems.

While only three primaries are necessary to produce mixtures of light with any chromaticity within their gamut, such systems cannot produce metamers, which are traditionally defined as two or more conditions having the same chromaticity but different SPDs. Metamers are important because they allow other performance attributes to vary while chromaticity is held constant, providing a useful tool for lighting practitioners to optimize performance without visible differences in the emitted light. Besides not delivering metamers, typical three-component mpLED systems do not offer adequate colour rendition for many architectural lighting applications.^9–12 Given the limitations of three primary systems, four and five primary mpLED systems are the most common commercial products today. Nevertheless, products with up to eight primaries are readily available, and lighting systems intended for research may have more than 30 unique primaries. LED systems with several primaries can generate a vast number of combinations that have identical chromaticity, but varying SPDs, hence varying colour and melanopic qualities, as shown in Figure 1.

Figure 1

Demonstration of different spectral power distributions of lights, each represented by a different line type, having identical chromaticity (3500 K on the Planckian locus) that could be generated by a multi-primary LED system

Even with only four primaries, physical measurement of all combinations is not practical. However, estimation of combinations is possible by measuring the output of each primary. Once each primary is characterized, combinations can be determined mathematically to estimate the true performance of the physical combinations. In some cases, one measurement of each primary can be scaled to predict its output at any given output-control (i.e. power or drive current) level, although non-linear relationships between drive current and output characteristics are common. More measurements can allow for interpolation to predict the output of a primary at any given control state (i.e. a coefficient x_i in the problem statement shown in Section 2 can be adjusted to accurately estimate the actual light output from each individual primary). This consideration is not within the scope of this article, which instead focuses on the enumeration aspect, rather than understanding the accuracy of estimations.

The evaluation of mpLED product performance is crucial in terms of market acceptance, usability and quality determination. A reliable and repeatable quantification of mpLED spectrum-based performance (e.g. colour rendition, melanopic stimulation) will enable manufacturers to build better performing products by optimizing the quantity and type of included primaries. There are two existing challenges: (1) determining ways to summarize performance that are useful and (2) efficiently computing the performance of estimated combinations. This article is focused on the latter, but with consideration for preliminary ideas for addressing the former.

1.2 Methods for quantifying estimated mpLED system performance

There are broadly two approaches to understanding the performance of estimated combinations. One is to perform single or multi-objective heuristic optimizations with specified performance criteria. This method produces a single result for each optimization that is used for addressing a specific goal, such as determining an experimental stimulus or, if multiple optimizations are completed, finding the Pareto boundary of two metrics. Computing optimal solutions is relatively fast – although not presently fast enough for real-time control – but requires some skill to program. This approach has been taken in numerous studies of LED system performance, with techniques including linear,⁹ non-linear,¹² stochastic,¹⁰ genetic,¹¹ multi-objective genetic,¹³ swarm¹⁴ and Nelder–Mead simplex¹⁵ algorithms.

Alternatively, brute-force enumeration techniques can be utilized to generate large numbers of SPDs. This approach determines specified performance characteristics for all (or many) combinations that can be created from an mpLED system. This approach gives a broader characterization of the trade-offs between different performance attributes, and through a lookup table has the potential to facilitate real-time performance tuning. Enumeration enables one to identify patterns, detect anomalies and perform statistically sound analysis,¹⁶ and is most useful when the specific aim of optimization is initially unknown or may vary. Carefully considered enumeration techniques can be preferable for understanding and reporting the performance of mpLED systems but can require substantial computation time and data storage. One recent work utilized a combination of heuristic optimizations and enumeration.¹⁷

One initial step to reduce computational burden of enumeration is to characterize system performance with a discrete set of metamers, rather than across the full range of chromaticity – the range of chromaticity can be easily determined and reported separately. This is the starting point of this work. Because the practical output of each primary occurs at discrete levels, this requires consideration of the precision of metamer determination. That is, with discrete control levels, only a small number of true metamers exist, but by defining a small tolerance – taking the form of a radius – for the chromaticity of the metamers (r_t), a much larger set of practical metamers can be created that are still visually equivalent in terms of the chromaticity of the light to a human observer. This region is deemed a target chromaticity area (TCA).

This article proposes a novel method for determining the combinations of primaries that are nominally metamers (i.e. have chromaticity within the TCA). The goal was to eliminate the need to determine the chromaticity of all possible combinations to check whether it falls within the TCA, and thus reduce the computational burden. To our knowledge, there is no alternative algorithm that serves this function; only objective-based optimizations and brute-force enumeration have previously been documented. After describing the proposed method, the article provides example results to demonstrate how it reduces the computation burden.

2. Proposed method

In an mpLED lighting system with n number of primaries, TSVs X, Y and Z are generated by scaling each primary using a coefficient (x_i), as shown in Equations (1)–(3):

X = X_{1} x_{1} + X_{2} x_{2} + \dots . . + X_{n} x_{n}

(1)

Y = Y_{1} x_{1} + Y_{2} x_{2} + \dots + Y_{n} x_{n}

(2)

Z = Z_{1} x_{1} + Z_{2} x_{2} + \dots + Z_{n} x_{n}

(3)

The TSVs are later converted to chromaticity scales and colour spaces. Solving this problem with a brute-force approach requires testing all combinations possible for the coefficients. For example, in a lighting system with 8-bit dimming control, each individual primary can be dimmed between 0 and 255 (x_i = 0, 1, 2, . . ., 255).

At its core, the proposed ER algorithm is a way to predict metamers with calculations that are less burdensome than computing the chromaticity coordinates of every possible combination of primaries comprising a composite SPD. This is possible given the geometric relationship between TSVs and chromaticity coordinates. Chromaticity coordinates provide quantification of a colour stimulus generated by a light source in a two-dimensional diagram. The International Commission on Illumination (CIE) defines the CIE 1976 $(u^{'}, v^{'})$ uniform chromaticity scale diagram coordinates as shown in Equations (4) and (5):

u^{'} = \frac{4 X}{X + 15 Y + 3 Z}

(4)

v^{'} = \frac{9 Y}{X + 15 Y + 3 Z}

(5)

where X, Y and Z are TSVs, which are derived from the sensitivity of the human visual system. For more information, see CIE 15:2018.¹⁸

The proposed algorithm avoids determining the chromaticity coordinates for every possible spectral output by first partitioning the TSVs of each individual LED primary into three groups based on the dominant coordinate $X, Y o r Z$ . The values within each of those groups are then sorted by their directional distance from planes forming a geometric shape, called the limiting pyramid. The sorted lists are divided into fixed width bins and the additive property of directional distances is used to limit the ranges of combinations from each of the coordinate-dominant groups that possibly map into the TCA. The TCA can be any region in a chromaticity diagram.

The definition of the limiting pyramid and directional distances is provided in the following sections. First, a simpler notion of a limiting quadrilateral will be presented in two dimensions in the $X, Y$ plane of TSV space and the same concept will be extended to the limiting pyramid in three dimensions.

2.1 Geometric analysis of the chromaticity diagram

2.1.1 The limiting quadrilateral

The desired $(u^{'}, v^{'})$ limits for a test SPD are given in Equations (6) and (7):

{u^{'}}_{a} \leq u^{'} \leq {u^{'}}_{b}

(6)

{v^{'}}_{a} \leq v^{'} \leq {v^{'}}_{b}

(7)

where ${u^{'}}_{a} = {u^{'}}_{t} - r_{t}$ , ${u^{'}}_{b} = {u^{'}}_{t} + r_{t}$ , ${v^{'}}_{a} = {v^{'}}_{t} - r_{t}$ and ${v^{'}}_{b} = {v^{'}}_{t} + r_{t}$ . The ${u^{'}}_{a}$ can be rewritten as Equation (8):

{u^{'}}_{a} \leq \frac{4 X}{X + 15 Y + 3 Z}

(8)

and since the TSVs X, Y and Z are all positive, this can become Equation (9):

(X + 15 Y + 3 Z) {u^{'}}_{a} \leq 4 X

(9)

Subtracting $(X + 3 Z) {u^{'}}_{a}$ from both sides of Equation (9) and dividing by $15 {u^{'}}_{a}$ results in Equation (10):

Y \leq (\frac{4 - {u^{'}}_{a}}{15 {u^{'}}_{a}}) X - \frac{Z}{5}

(10)

Similarly, Equation (11) can be derived from the $u^{'} \leq {u^{'}}_{b}$ inequality:

Y \geq (\frac{4 - {u^{'}}_{b}}{15 {u^{'}}_{b}}) X - \frac{Z}{5}

(11)

and the inequalities for bounds on v′ provides the following two relationships, provided in Equations (12) and (13):

Y \leq (\frac{{v^{'}}_{b}}{9 - 15 {v^{'}}_{b}}) X + (\frac{{v^{'}}_{b}}{3 - 5 {v^{'}}_{b}}) Z

(12)

Y \geq (\frac{{v^{'}}_{a}}{9 - 15 {v^{'}}_{a}}) X + (\frac{{v^{'}}_{a}}{3 - 5 {v^{'}}_{a}}) Z

(13)

The inequalities in Equations (12) and (13) are elucidated in the X-Y plane for some fixed value of Z. For a fixed value of Z, the lines in the X-Y plane can be defined using Equations (10) through (13) where the limiting lines can be referred to as ${u^{'}}_{a}, {u^{'}}_{b}, {v^{'}}_{b} and {v^{'}}_{a}$ , respectively. In short, the collective terms for these lines are the X-Y limiting lines.

Thus, for a fixed value of Z, values of X and Y that will result in chromaticity within the TCA must lie below the ${u^{'}}_{a}$ line and above the ${u^{'}}_{b}$ line as well as lie below the ${v^{'}}_{b}$ line and above the ${v^{'}}_{a}$ line, and hence must lie in the limiting quadrilateral defined by intersections of the four X-Y limiting lines, as shown in Figure 2. It should be noted that the slopes of the X-Y limiting lines are independent of Z, and only their intercepts change with Z. When Z increases, the limiting quadrilateral will grow slightly larger, maintain the same orientation, and move upward and to the right. Because all the X-Y limiting line slopes are positive, the minimum X and Y coordinates of the limiting quadrilateral occur at the intersection of the u_a and v_a lines with coordinates given by (X_min, Y_min). The maximum X and Y coordinates of the limiting quadrilateral are at the intersection of the u_b and v_b lines, (X_max, Y_max). Allowing m_a be the slope of the line u_a, m_b be the slope of the line u_b, p_a be the slope of the line v_a, and p_b be the slope of the line v_b, the equations of the X-Y limiting lines can be rewritten as shown in Equations (14) through (22). For line $u_{a}$ :

Y = m_{a} X - 0.2 Z

(14)

Figure 2

The limiting quadrilateral showing the boundaries of the tristimulus combinations selected for further analysis

where

m_{a} = \frac{4 - {u^{'}}_{a}}{15 {u^{'}}_{a}}

(15)

For line $u_{b}$ :

Y = m_{b} X - 0.2 Z

(16)

where

m_{b} = \frac{4 - {u^{'}}_{b}}{15 {u^{'}}_{b}}

(17)

For line $v_{a}$ :

Y = p_{a} X + (\frac{{v^{'}}_{a}}{3 - 5 {v^{'}}_{a}}) Z

(18)

where

p_{a} = \frac{{v^{'}}_{a}}{9 - 15 {v^{'}}_{a}}

(19)

For line $v_{b}$ :

Y = p_{b} X + (\frac{{v^{'}}_{b}}{3 - 5 {v^{'}}_{b}}) Z

(20)

where

p_{b} = \frac{{v^{'}}_{b}}{9 - 15 {v^{'}}_{b}}

(21)

The coordinates of the extreme corners of the limiting quadrilateral are $(X_{\min}, Y_{\min})$ at the intersection of lines $u_{a}$ and $v_{a}$ , and $(X_{\max}, Y_{\max})$ at the intersection of lines $u_{b}$ and $v_{b}$ . These are given by Equations (22) through (25):

X_{\min} = (\frac{0.2 + 3 p_{a}}{m_{a} - p_{a}}) Z

(22)

Y_{\min} = p_{a} * X_{\min} + 3 p_{a} Z

(23)

X_{\max} = (\frac{0.2 + 3 p_{b}}{m_{b} - p_{b}}) Z

(24)

Y_{\max} = p_{b} X_{\max} + 3 p_{b} Z

(25)

Once the Z component of the tristimulus triple is fixed, there are fairly narrow restrictions on the X and Y components for the triple to map to the TCA. It is also apparent that for a given set of TSVs, there is a maximum amount of X and Y components available (with all primaries at maximum output). The combination of the formulae for X_min and Y_min yields the limit for the maximum amount of Z component that can be balanced by the available X and Y components to position the combinations within the TCA. The sum of the TSV X components for the included primaries can be denoted by X_tot and the sum of the TSV Y components can be denoted by Y_tot. Since X_min ⩽ X_tot and Y_min ⩽ Y_tot, the geometric boundaries of the quadrilateral can be described as shown in Equations (26) and (27)

Z_{\max} \leq (\frac{m_{a} - p_{a}}{0.2 + 3 p_{a}}) X_{tot}

(26)

Z_{\max} \leq (\frac{m_{a} - p_{a}}{p_{a} (0.2 + 3 m_{a})}) Y_{tot}

(27)

The smallest of these two bounds is Z_max.

2.1.2 The limiting pyramid

From the equations for the slopes of the $X - Y$ limiting lines (Equations (12), (14), (16) and (18)), it is evident that slopes are independent of $Z$ , depending only on the TCA location and its radius (that define the four values ${u^{'}}_{a}$ , $u_{b}$ , ${v^{'}}_{a}$ and ${v^{'}}_{b}$ ). Thus, as the value of $Z$ shifts from 0 to $Z_{m a x}$ , the limiting quadrilateral sweeps out an inverted pyramid with apex at the origin and its base is the limiting quadrilateral in the $Z = Z_{m a x}$ plane. This structure is termed the limiting pyramid and is depicted in Figure 3 by the shaded blue region and dashed blue and red lines. The crux of the new algorithm depicted here is that combinations of the TSVs that will map to the TCA must lie in the limiting pyramid.

Figure 3

The limiting pyramid shows the three-dimensional boundary for the tristimulus values of the combined spectral power distributions

2.2 Directional distance

Directional distance between a point and a line or a plane is simply the Euclidian distance between them with a sign attached indicating on which side of the line or plane the point lies. This notion is first illustrated with an example in two dimensions and extended to three dimensions in Section 2.2.2.

2.2.1 Directional distance in two dimensions

Figure 4 shows a graph of the line $y = 2 x$ , for $x > 0$ and the point (2,4) on the line, and two points off the line, (0,5) and (4,3). It is useful in discussing the Euclidian distance between a point and a line to refer to the standard form of an equation for a line, provided as Equation (28), in variables $x$ and $y$ :

a x + b y + c = 0

(28)

Figure 4

Directional distance, the distance between a point and y = 2x line, is used to sort out the outcomes on the limiting quadrilateral

This line can be denoted as $line (a, b, c)$ . Equation (28) can be multiplied by any constant without changing the position of the line. If it is the case that $a^{2} + b^{2} + c^{2} = 1$ , Equation (28) will be the normalized standard form of the line. The Euclidian distance from point $(x_{1,} y_{1})$ to $line (a, b, c)$ , denoted $d ((x_{1,} y_{1}), l i n e (a, b, c))$ , is shown in Equation (29):

d ((x_{1}, y_{1}), line(a, b, c)) = \frac{| a x_{1} + b y_{1} + c |}{\sqrt{a^{2} + b^{2} + c^{2}}}

(29)

A normalized standard form of the line is shown in Equation (30):

d ((x_{1}, y_{1}), line (a, b, c)) = | a x_{1} + b y_{1} + c |

(30)

The directional distance, which is denoted $D ((x_{1}, y_{1}), l i n e (a, b, c))$ , is found by removing the absolute value from the Euclidian distance calculation, as shown in Equation (31):

D ((x_{1}, y_{1}), line (a, b, c)) = a x_{1} + b y_{1} + c

(31)

In Figure 3, the normalized standard form of $y = 2 x$ is given as Equation (32):

\frac{- 2 \sqrt{5}}{5} x + \frac{\sqrt{5}}{5} y + 0 = 0

(32)

The directional distance from point (0,5) to $line (- 2 \sqrt{5} / 5, \sqrt{5} / 5, 0)$ is as provided in Equation (33):

\begin{matrix} D ((0, 5), line (\frac{- 2 \sqrt{5}}{5}, \frac{\sqrt{5}}{5}, 0)) = (\frac{- 2 \sqrt{5}}{5} \times 0) \\ + (\frac{\sqrt{5}}{5} \times 5) + 0 \\ = \sqrt{5} \end{matrix}

(33)

The directional distance from point (4,3) to $line (- 2 \sqrt{5} / 5, \sqrt{5} / 5, 0)$ is given in Equation (34):

\begin{matrix} D ((4, 3), line (\frac{- 2 \sqrt{5}}{5}, \frac{\sqrt{5}}{5}, 0)) = (\frac{- 2 \sqrt{5}}{5} \times 4) \\ + (\frac{\sqrt{5}}{5} \times 3) + 0 \\ = - \sqrt{5} \end{matrix}

(34)

2.2.2 Directional distance in three dimensions

Directional distance in three dimensions is directly analogous to the two-dimensional definition. The standard form for the equation of a plane in variables $x, y and z$ is given in Equation (35):

a x + b y + c z + d = 0

(35)

A particular plane specified in Equation (32) can be denoted as $plane (a, b, c, d)$ . If Equation (32) has been normalized so that $\sqrt{a^{2} + b^{2} + c^{2} + d^{2}} = 1$ , the directional distance from point $(x_{1}, y_{1}, z_{1})$ to $plane (a, b, c, d), D ((x_{1}, y_{1}, z_{1}), plane (a, b, c, d))$ is given in Equation (36):

D ((x_{1}, y_{1}, z_{1}), plane (a, b, c, d)) = a x_{1} + b y_{1} + c z_{1} + d

(36)

The planes of interest in the enumeration process will be the planes swept out by the $v_{b}$ and $u_{b}$ $X, Y$ limiting lines, which will be referred to as the v_b plane and u_b plane, respectively.

2.2.3 Additive property of directional distances

An important property of directional distances that they are additive. When forming combinations of two TSVs, the directional distance of their sum from a given plane is the sum of their directional distances from that plane. Equation (37) follows directly from Equation (36):

\begin{array}{l} D ((x_{1} + x_{2}, y_{1} + y_{2}, z_{1} + z_{2}), plane (a, b, c, d)) = D ((x_{1}, y_{1}, z_{1}), plane (a, b, c, d)) + \\ D ((x_{2}, y_{2}, z_{2}), plane (a, b, c, d)) \end{array}

(37)

2.3 ER algorithm

A rough outline of the process is given below (steps 1–9), and subsequent subsections will elaborate on the individual steps.

Determine the normalized coefficients of the v_b plane and u_b plane.

Sort the primaries into three groups based on the dominant coordinate of their TSV (X_p1, Y_p1, Z_p1), (X_p2, Y_p2, Z_p2), . . ., (X_pn, Y_pn, Z_pn) for n number of primaries.

Within each group, form all possible output and primary SPD combinations (X_x,_com, Y_x,_comp, Z_x,_com), (X_y,_com, Y_y,_com, Z_y,_com), (X_z,_com, Y_z,_com, Z_z,_com).

Sort the Z-dominant combinations by their TSV Z component (X_z,_com, Y_z,_com, Z_z,_com) and store the sorted list of TSV combinations in a file removing those combinations with Z component greater than $Z_{\max}$ .

Sort the Y-dominant combinations by their directional distance from the v_b plane (X_y,_com, Y_y,_com, Z_y,_com) and store the sorted list of TSVs. Also bin the sorted list into a number of fixed width bins, determining the index in the sorted list of the first element in each bin. These bin indices, the bin width, the minimum and maximum directional distance from the v_b plane, and the minimum and maximum directional distance from the u_b plane are stored.

Similarly sort the X-dominant combinations by their directional distance from the u_b plane (X_x,_com, Y_x,_com, Z_x,_com), storing the sorted list along with the two directional distances for each combination and the minimum and maximum directional distances from each of the two limiting planes. Also bin this sorted list into fixed width bins, storing the bin width and bin indices in the sorted list.

For each Z-dominant combination, compute its directional distance from the $u_{b}$ and $v_{b}$ planes. Use these and the bounds on the directional distances of the X-dominant combinations to limit the bins of Y-dominant combinations that need to be examined.

For each Z-dominant plus Y-dominant combination, denoted as a $Z - Y$ combination, determine the directional distances from the $u_{b}$ and $v_{b}$ planes (by adding their individual directional distances). Use this and the bin width and bin indices for the X-dominant combinations to restrict which X-dominant combinations to test with this $Z - Y$ combination.

For each of these $Z - Y$ combination plus X-dominant combinations, denoted a $Z - Y - X$ combination, use Equations (4) and (5) to determine the $u^{'}$ and $v^{'}$ values for the combined TSV and test to see whether it lies within the TCA.

2.3.1 The $v_{b}$ and $u_{b}$ plane equations

The equations of the $v_{b}$ and $u_{b}$ planes may be determined from the coordinates of any three non-linear points that lie on the plane with combinations of determinants from groups of the coordinates of the three points. A lower–upper factorization¹⁹ with total pivoting was used to evaluate the determinants, which is more numerically stable compared to other algorithms, such as Cramer’s Rule.²⁰

The points for the $v_{b}$ plane are the origin and the points $(X_{\max}, Y_{\max}, Z)$ and $(X_{d}, Y_{d}, Z)$ in Figure 1 for some value of $Z$ . Although the value of $Z$ is not critical, it is important to choose a relatively small Z for numerical stability of the determinant computation. In this enumeration process, the value $Z = 1$ is determined to be sufficient. The equation scaling is chosen so that the directional distance from $(X_{c}, Y_{c}, Z)$ to the $v_{b}$ plane is negative.

The points used to determine the equation of the $u_{b}$ plane are the origin and the points $(X_{\max}, Y_{\max}, Z)$ and $(X_{c}, Y_{c}, Z)$ . The equation scaling is chosen so that the directional distance from $(X_{d}, Y_{d}, Z)$ to the $u_{b}$ plane is positive.

2.3.2 The coordinate-dominant regions of TSV space

All possible combinations of a set of TSVs at all discrete output-control levels lie in a box with lower left corner at the origin (all output levels are zero) and upper right corner $(X_{tot}, Y_{tot}, Z_{tot})$ as shown in Figure 5. The distribution of TSVs for these various combinations is non-uniform through this limiting box. As the TSVs are grouped by their dominant coordinate partition, the box is partitioned into three regions depicted in Figures 6 –8. The diagonal planes forming the triangular faces of the depicted regions are the $x = y, x = z and y = z$ planes.

Figure 5

All tristimulus value combinations range showing the boundaries of X_tot, Y_tot and Z_tot in a box

Figure 6

X-dominant tristimulus value combinations range

Figure 7

Y-dominant tristimulus value combinations range

Figure 8

Z-dominant tristimulus value combinations range

2.3.3 The sorted Z-dominant TSVs

For combinations of the Z-dominant TSVs, three coordinates in three eight-byte binary floating-point numbers (doubles) are stored, as the directional distances are computed just once in step 7 of the ER method. Any combination that has a $Z$ component greater than $Z_{\max}$ is not stored. The sorting is done with a merge sort. If $n_{z}$ is the number of Z-dominant TSVs that is available and the number of output-control levels per TSV is $m$ , then the total number of possible Z-dominant combinations, $N_{z}$ , is ${(m + 1)}^{n_{Z}}$ as a output-control level of 0 for each TSV is allowed. That yields a maximum size for the Z-dominant combinations of $24 N_{z}$ bytes.

2.3.3 The sorted Y-dominant TSVs

The Y-dominant TSVs and their directional distances from the $u_{b}$ and $v_{b}$ planes are computed and stored. Hence, five eight-byte numbers for each Y-dominant combination must be stored. This list is sorted based on the directional distance of the TSV from the $v_{b}$ plane. As the directional distances from the $u_{b}$ and $v_{b}$ planes are computed, the minimum and maximum of these distances in $Y_{u_{\min}}$ , $Y_{u_{\max}}$ , $Y_{v_{\min}}$ and $Y_{v_{\max}}$ for use in step 8 of the ER algorithm are recorded. After sorting the Y-dominant TSVs, they are binned into equal width bins, storing the index in the sorted list of the smallest element in each bin in an array which is termed an array of bin pointers.

If there are $n_{y}$ Y-dominant TSVs, then the total number of possible Y-dominant combinations, $N_{Y}$ , is ${(m + 1)}^{n_{y}}$ . The number of bin pointers (B_Y) is chosen to be the smaller of 4096 and the largest power of 2 less than or equal to $N_{Y} / 4$ . The bin width, $W_{Y}$ , is $(Y_{v_{\max}} - Y_{v_{\min}}) / (B_{Y} - 1)$ . The number of bins would be B_Y − 1 so that the last index points to one beyond the last Y-dominant combination and successive bin pointers delimit the combinations within a bin. The total storage for the sorted Y-dominant TSVs is $40 N_{Y} + 8 B_{Y} + 44$ with the last 44 bytes used to store the minimums and maximums of the directional distances and the number of bins and the bin width.

2.3.4 The sorted X-dominant TSVs

The X-dominant TSVs and their directional distances from the $u_{b}$ and $v_{b}$ planes are computed and stored. Hence, five eight-byte numbers for each X-dominant combination are stored. This list is sorted based on the directional distance of the TSV from the $u_{b}$ plane. As the directional distances from the $u_{b}$ and $v_{b}$ planes are computed, the minimum and maximum of these distances in $X_{u_{\min}}$ , $X_{u_{\max}}$ , $X_{v_{\min}}$ and $X_{v_{\max}}$ for use in step 8 of the ER algorithm are recorded. After sorting the X-dominant TSVs, they are binned into equal width bins, storing the index in the sorted list of the smallest element in each bin in an array which is termed an array of bin pointers.

If there are $n_{x}$ X-dominant TSVs, the total number of possible X-dominant combinations, $N_{X}$ , is ${(m + 1)}^{n_{x}}$ . The number of bin pointers (B_X) is chosen to be the smaller than 4096 and the largest output of 2 less than or equal to $N_{X} / 4$ . The bin width, $W_{X}$ , is $(X_{u_{\max}} - X_{u_{\min}}) / (B_{X} - 1)$ . The number of bins is B_X − 1 so that the last index points to one beyond the last X-dominant combination and successive bin pointers delimit the combinations within a bin. The total storage for the sorted X-dominant TSVs is $40 N_{X} + 8 B_{X} + 44$ with the last 44 bytes used to store the minimums and maximums of the directional distances and the number of bins and the bin width.

2.3.5 Combining dominant-coordinate group combinations that could land in the TCA

In steps 7–9 of the of the ER algorithm forming all the $N_{Z} \times N_{Y} \times N_{X}$ combinations of the TSVs should be avoided because a TSV combination can only map to the TCA if it lies inside the limiting pyramid. This means that a TSV combination lies between the $v_{a}$ and $v_{b}$ planes and between the $u_{a}$ and $u_{b}$ planes. The previous discussion is limited to computing and storing the distances from the $u_{b}$ and $v_{b}$ planes due to the approximation of the $u_{a}$ and $v_{a}$ planes by planes parallel to the $u_{b}$ and $v_{b}$ planes. Therefore, the limiting pyramid becomes a limiting parallelepiped that contains the limiting pyramid. While this does increase the volume of this limiting geometric figure, it halves the number of directional distances needed to compute and store for the coordinate-dominant combinations formed in steps 4, 5 and 6 of the ER algorithm. The distance between these approximating planes is the maximum distance between the corresponding planes in the range from $Z = 0$ to $Z = Z_{\max}$ .

If the maximum distance between the $v_{a}$ and $v_{b}$ planes in that range are denoted ∆_v given by Equation (38) where $D ()$ refers to the directional distance of Section 2.2.2. Note that ∆_v is negative.

Δ_{v} = D ((X_{c}, Y_{c}, Z_{\max}), v_{b} - plane)

(38)

Similarly, the maximum distance between the $u_{a}$ and $u_{b}$ planes can be defined as ∆_u given by Equation (39). ∆_u is positive in our definition of the equation of the $u_{b}$ plane.

Δ_{u} = D ((X_{d}, Y_{d}, Z_{\max}), u_{b} plane)

(39)

Thus, a tristimulus combination that maps to the TCA must have its directional distance from the $v_{b}$ plane in the closed interval [∆_v, 0] and its directional distance from the $u_{b}$ plane must be in the closed interval [0, ∆_u]. For a given Z-dominant TSV, if its directional distance from the $v_{b}$ plane, say $σ_{v}$ , are computed, it is adequate to analyse only Y-dominant combinations whose distance from the $v_{b}$ plane, say $ρ_{v}$ , satisfies the inequality shown in Equation (40):

Δ_{v} - X_{v_{\max}} \leq σ_{v} + ρ_{v} \leq - X_{v_{\min}}

(40)

The $X_{v_{\min}}$ and $X_{v_{\max}}$ in Equation (40) are included to allow for the possible range that the yet to be added X-dominant TSV combination might shift the overall directional distance to the $v_{b}$ plane. Subtracting $σ_{v}$ from all parts of Equation (37) shows that $ρ_{v}$ must lie in the interval shown in Equation (41):

Δ_{v} - X_{v_{\max}} - σ_{v} \leq ρ_{v} \leq - X_{v_{\min}} - σ_{v}

(41)

As the Y-dominant TSVs are sorted by their directional distance from the $v_{b}$ plane, and those directional distances are in the interval $[Y_{v_{\min}}, Y_{v_{\max}}]$ , and have been assigned into $B_{Y} - 1$ bins of width $W_{Y}$ as discussed in Section 2.3.4, it is possible to determine the range of bins of the Y-dominant TSVs that hold combinations that when combined with a given Z-dominant TSV and a yet to be determined X-dominant TSV could map to the TCA. The first bin index, $b_{0}$ is given by Equation (42):

b_{0} = \max (0, floor (\frac{Δ_{v} - X_{v_{\max}} - σ_{v} - Y_{v_{\min}}}{W_{Y}}))

(42)

where $floor (x)$ is the largest integer less than or equal to $x$ .

The last bin index, $b_{1}$ , is given by Equation (43), where $ceil (x)$ is the smallest integer greater than or equal to $x$ .

b_{1} = \min (B_{Y} - 1, ceil (\frac{- X_{v_{\min}} - σ_{v} - Y_{v_{\min}}}{W_{Y}}))

(43)

In a similar fashion, for a given Z-dominant TSV combination paired with each Y-dominant TSV combination in bins $b_{0}$ through $b_{1}$ , their combined directional distance from the $u_{b}$ plane, say $σ_{u}$ , is determined and it is necessary to limit the search of X-dominant TSV combinations with directional distance from the $u_{b}$ plane, $ρ_{u}$ , to be such that $σ_{u} + ρ_{u}$ lies in the interval [0, ∆_u]. This yields the interval shown in Equation (44) for $ρ_{u}$ .

- σ_{u} \leq ρ_{u} \leq Δ_{u} - σ_{u}

(44)

Analogously, and with the definitions in Section 2.3.5, it is evident that the delimiting indices for the X-dominant TSV combinations bins, $c_{0}$ and $c_{1}$ , are given in Equations (45) and (46):

c_{0} = \max (0, floor (\frac{- σ_{u} - X_{u_{\min}}}{W_{X}}))

(45)

c_{1} = \min (B_{X} - 1, ceil (\frac{Δ_{u} - σ_{u} - X_{u_{\min}}}{W_{X}}))

(46)

2.3.6 Parallel implementation of the ER algorithm

Steps 1 through 6 in the ER algorithm take very little time and are run serially in a program called gen_combos to generate the sorted coordinate-dominant TSV groups with their corresponding directional distances from the $u_{b}$ and $v_{b}$ planes and their bin pointers. Steps 7–9 can be executed in parallel in a message passing interface-based C program.²¹ An individual task is all the comparisons for a particular Z-dominant combination. They can be round robin scheduled using the order of the Z-dominant TSVs generated in gen_combos. Before starting the tasks, the program reads in and distributes the three coordinate-dominant TSV databases, storing them in a shared memory segment on each node to facilitate keeping only one copy on each of the shared memory nodes participating in the parallel job.

3 Performance evaluation

A test was performed to capture the benefits of the ER algorithm by comparing the calculation time needed to generate a metamer set for the proposed method to that of a full enumeration. In all, 12 theoretical mpLED systems were established by calculating primaries as Gaussian distributions. The 12 mpLED systems included three systems with four primaries, three systems with five primaries, three systems with six primaries and three systems with seven primaries. The data comprising these systems are available as supplemental material. A TCA was centred on the Planckian locus with a colour temperature of 3500 K, chromaticity coordinates $({u^{'}}_{t}, v_{t}) = (0.2357, 0.5113)$ , and radius r_t = 0.0005.

For each mpLED system, a full enumeration was calculated using the Constance computer previously described, with the chromaticity of each spectral combination of LED primaries checked for compliance with the TCA criterion of 3500 K on Planckian locus with a $Δ {u^{'}}_{t}, v_{t}$ tolerance of 0.0005. The previously linked code includes a parameter choice to complete this full calculation. The ER algorithm was then used to generate the same metamer set. The accuracy of the algorithm was confirmed by comparing the details of the metamer sets; no errors were found (i.e. the ER algorithm did not miss any combinations that do fall within the TCA). Table 1 summarizes the results, including time reduction provided by the ER algorithm.

Table 1

Results of wall clock time comparisons for the ER algorithm and full calculation of all combinations

Number of primaries	Combinations in TCA	Nodes	Cores/node	ER compute time (seconds)	Full compute time (seconds)	ER time multiple
4	1.3880E+03	1	24	0.9	1.0	1.1
4	7.3700E+02	1	24	1.3	1.1	0.8
4	3.1430E+03	1	24	1.3	1.2	0.9
5	2.1110E+06	1	24	1.2	10.1	8.7
5	7.0098E+05	1	24	1.5	10.0	6.7
5	6.2186E+05	1	24	1.2	9.1	7.6
6	2.8914E+08	1	24	45.3	985.5	21.8
6	1.4931E+08	1	24	26.3	779.0	29.6
6	6.1680E+06	1	24	2.2	757.6	344.4
7	9.7919E+09	1	24	1586.3	83816.6	52.8
7	1.8840E+10	1	24	2632.7	76475.1	29.0
7	1.3174E+09	1	24	82.9	76407.7	921.7

ER: enumeration reduction; TCA: target chromaticity area.

As shown in Table 1, the ER algorithm reduced the time to find the set of metamers (i.e. combinations with a chromaticity within the TCA) by up to 99.89%. However, the reduction in time was highly variable, depending on both the number of primaries and the number of combinations within the TCA. The number of combinations in the TCA computed on one core, not provided, is also a factor. Likewise, the full computation varies slightly for systems with the same number of primaries due to differences in parallel computing.

For four primary systems, the ER algorithm offered no benefit because the time needed to determine and sort the groups of TSVs offset any reduction in chromaticity coordinate calculations. Likewise, there were only modest benefits for five primary systems. As the number of primaries increased, the benefit of the ER algorithm also increased, and the algorithm is anticipated to be essential for understanding the performance of research-grade mpLED systems with eight or more primaries. Table 1 is limited to a maximum of seven primaries because the cost and calculation time associated with a full enumeration of combinations for eight primary systems became unjustifiable. Running on 30 nodes, the ER algorithm for three computed eight primary mpLED systems took between 130 and 12,026 seconds. The calculation time for the baseline full enumeration of an eight primary mpLED systems could be expected to take days, even on a high-performance computer.

The ER algorithm showed more benefit for mpLED systems with fewer combinations falling within the TCA for a given number of primaries because more chromaticity coordinate calculations could be avoided. It is important to remember that the ER algorithm does not result in chromaticity coordinates itself, but a set of combinations for which chromaticity is calculated. The set is slightly larger than what is in the TCA, so some combinations identified by the ER algorithm are ultimately not kept in the metamer set.

Two possible improvements to the ER algorithm are apparent with a possibility of additional overhead. The first is to sort the Z-dominant TSVs by their directional distance from the $v_{b}$ plane instead of their $Z$ component as currently proposed. This might improve the spatial location of access to the Y-dominant bins over time. A second possible improvement would be to sort and bin the Y-dominant TSVs by their directional distance from the $u_{b}$ plane within each bin sorted by distance from the $v_{b}$ plane. Similarly, it is also possible to sort and bin the X-dominant TSVs by their directional distance from the $v_{b}$ plane within each bin sorted by distance from the $u_{b}$ plane.

4 Conclusions

While tunable, mpLED systems can emit a vast number of spectral outputs, calculating and describing these lighting conditions requires a tremendous amount of computational time and data storage. Among a broader set of methods being developed to improve the efficiency and impact of mpLED system assessments, a new method, denoted the ER algorithm, was developed to predict the combinations of LED primaries that will fall within a specified area of a chromaticity diagram. Using this method, chromaticity coordinates for only a small subset of the total spectral outputs that can be generated by a system need to be calculated to identify a set of nominal metamers (i.e. having chromaticity coordinates within the specified area). For a set of 12 theoretical mpLED systems, the ER algorithm reduced the time to determine a set of nominal metamers by up to 99.89%. More benefit was found for systems with more primaries, and with a lower percentage of all possible combinations having chromaticity coordinates within the area of interest, which varies based on the specific primaries included in the system. The ER algorithm offers substantial benefit for researchers trying to understand the performance of products and for engineers trying to quantify and market product performance.

Footnotes

Declaration of conflicting interests

The authors declared no potential conflicts of interest with respect to the research, authorship and/or publication of this article.

Funding

The authors disclosed receipt of the following financial support for the research, authorship and/or publication of this article: This work was supported by the U.S. Department of Energy’s Lighting R&D Program, part of the Building Technologies Office within the Office of Energy Efficiency and Renewable Energy (EERE).

ORCID iDs

D Durmus

M Royer

Replication of results

Details of the optimization algorithm have been provided in this manuscript to enable others to replicate the results of the proposed method. Due to the large file size of the data collected using this algorithm, the raw calculation data have been omitted, except for one example made available via GitHub. Code for the ER algorithm is available from GitHub: . Example data are also provided.

References

Sisto

Gauvin

Accurate chromatic control and color rendering optimization in LED lighting systems using junction temperature feedback. Thirteenth International Conference on Solid State Lighting, San Diego, CA, USA, 18 September 2014, 9190: 919002.

Sánchez-Cano

Aporta

Optimization of lighting projects including photopic and circadian criteria: a simplified action protocol. Applied Sciences 2020; 10: 8068.

Durmus

Davis

Object color naturalness and attractiveness with spectrally optimized illumination. Optics Express 2017; 25: 12839–12850.

Aldrich

Zhao

Paradiso

Energy efficient control of polychromatic solid state lighting using a sensor network. In Tenth International Conference on Solid State Lighting, San Diego, CA, USA, 18 August 2010, 7784: 17–31.

Hertog

Llenas

Carreras

Optimizing indoor illumination quality and energy efficiency using a spectrally tunable lighting system to augment natural daylight. Optics Express 2015; 23: A1564–A1574.

Huang

Laser and LED hybrid plant lighting system design based on the particle swarm algorithm. Applied Sciences 2020; 10: 7588.

Pattison

Hansen

Bardsley

Elliott

Lee

Pattison

, et al. 2019 Lighting R&D Opportunities. Santa Barbara, CA: U.S. DOE Office of Energy Efficiency and Renewable Energy (EERE), Energy Efficiency Office, 2020.

National Electrical Manufacturers Association. American National Standard for Electric Lamps-Specifications for the Chromaticity of Solid-State Lighting Products. Rosslyn, VA: NEMA, 2017.

Ohno

Spectral design considerations for white LED color rendering. Optical Engineering 2005; 44: 111302.

10.

Žukauskas

Vaicekauskas

Ivanauskas

Gaska

Shur

MS.

Optimization of white polychromatic semiconductor lamps. Applied Physics Letters 2002; 80: 234–236.

11.

Smet

Ryckaert

Pointer

Deconinck

Hanselaer

Optimal colour quality of LED clusters based on memory colours. Optics Express 2011; 19: 6903–6912.

12.

Royer

Evaluating tradeoffs between energy efficiency and color rendition. OSA Continuum 2019; 2: 2308–2327.

13.

Durmus

Abdalla

Duis

Davis

Spectral optimization to minimize light absorbed by artwork. Leukos 2020; 16: 45–54.

14.

Luo

Chen

Liu

ST.

Color-tunable light emitting diodes based on quantum dot suspension. Applied Optics 2015; 54: 2845–2850.

15.

Vázquez

Fernández-Balbuena

Canabal

Muro

Durmus

Davis

, et al. Energy optimization of a light projection system for buildings that virtually restores artworks. Digital Applications in Archaeology and Cultural Heritage 2020; 16: e00128.

16.

Bradley

Fayyad

Mangasarian

OL.

Mathematical programming for data mining: formulations and challenges. INFORMS Journal on Computing 1999; 11: 217–238.

17.

Zandi

Stefani

Herzog

Schlangen

Trinh

Khanh

TQ.

Optimising metameric spectra for integrative lighting to modulate the circadian system without affecting visual appearance. Scientific Reports 2021; 11: 1–4.

18.

International Commission on Illumination. Colorimetry. CIE 15:2018. Vienna, Austria: CIE, 2015.

19.

Schwarzenberg-Czerny

On matrix factorization and efficient least squares solution. Astronomy and Astrophysics Supplement Series 1995; 110: 405.

20.

Trefethen

Bau

III . Numerical Linear Algebra, 1st edition. Philadelphia, PA: SIAM: Society for Industrial and Applied Mathematics, 1997.

21.

Pacheco

Parallel Programming with MPI. San Francisco, CA: Morgan Kaufmann, 1997.

Enumeration reduction algorithm for the characterization of multi-primary LED systems

Abstract

1. Introduction

1.1 Multi-primary lighting systems

1.2 Methods for quantifying estimated mpLED system performance

2. Proposed method

2.1 Geometric analysis of the chromaticity diagram

2.1.1 The limiting quadrilateral

2.1.2 The limiting pyramid

2.2 Directional distance

2.2.1 Directional distance in two dimensions

2.2.2 Directional distance in three dimensions

2.2.3 Additive property of directional distances

2.3 ER algorithm

2.3.1 The v b and u b plane equations

2.3.2 The coordinate-dominant regions of TSV space

2.3.3 The sorted Z-dominant TSVs

2.3.3 The sorted Y-dominant TSVs

2.3.4 The sorted X-dominant TSVs

2.3.5 Combining dominant-coordinate group combinations that could land in the TCA

2.3.6 Parallel implementation of the ER algorithm

3 Performance evaluation

4 Conclusions

Footnotes

Declaration of conflicting interests

Funding

ORCID iDs

Replication of results

References

2.3.1 The $v_{b}$ and $u_{b}$ plane equations