*Published April 29, 2015 by Scott Allen Burns, last updated May 25, 2020*

**Change Log**

9/19/2017 Updated symbols to reflect common colormetric usage: instead of , instead of , instead of , instead of .

9/20/2017 Added explicit formula for in LSS method.

9/26/2018 Added link to a fast reflectance reconstruction method.

10/2/2018 Changed matrix to be called , which matches established usage.

10/10/2018 Changed to the more correct non-italic form .

10/18/2018 Added a line-by-line explanation of how to read the Matlab code for ILSS.

11/1/2018 Did some minor editing for clarification of some points.

5/23/2019 Added attribution to earlier work related to the LSS method. Added notice of upcoming publication in the

*Color Research and Application*journal.

6/4/2019 Added new LHTSS (Least Hyperbolic Tangent Slope Squared) method that gives best quality reflectances of all six methods. It appears at the end of the article.

12/25/2019 Added link to short note on implementing the Rec.BT.2020 RGB primaries.

1/9/20 Added notifications in the body pointing out that LHTSS is a better option than ILLSS in general.

1/18/20 Added link to CR&A paper.

3/29/20 Added preface and reference to TLSS.

4/13/20 Updated the preface.

5/11/20 Corrected a typo in the F expression for the LHTSS method. (Thanks Marnix!)

5/15/20 Added links to Excel spreadsheets that implement the LLSS and LHTSS methods.

5/25/20 Added link to the Wiley Online Library for the CR&A journal paper.

##### Preface

**A note to those working in the color science field:** This presentation is aimed at a general audience and will cover much background material with which you are already familiar. A much more concise version of the methods presented below is available in the journal article “Numerical Methods for Smoothest Reflectance Reconstruction,” which has been published in Color Research and Application (Vol 45, No 1, 2020, pp 8-21, DOI: 10.1002/col.22437). The paper PDF appears to be available at the journal’s web page, but if not, please contact me if you’d like for me to send you a copy.

Additional development is available online at the supplementary documentation web page. Note that the three methods presented in the journal paper (called simply methods 1, 2, and 3) correspond to the methods below called LSS, LLSS, and LHTSS, respectively.

*In more recent news, I have been in contact with Glenn Davis, who is in the process of incorporating the LLSS and LHTSS methods into his colorSpec library. He combines the two into one function called TLSS, which stands for Transformed Least Slope Squared. It is written in the R programming language. An input parameter is passed, which effectively selects between the LLSS method (for any color within the spectral locus or for emissive sources) and the LHTSS method (for object colors produced by an illuminant reflecting off of a real object). He also implements some shortcuts to avoid explicit construction of diagonal matrices.*

#### Overview

I present several algorithms for generating a reflectance curve from a specified sRGB triplet, written for a general audience. Although there are an infinite number of reflectance curves that can give rise to the specific color sensation associated with an sRGB triplet, the algorithms presented here are designed to generate reflectance curves that are similar to those found with naturally occurring colored objects. My hypothesis is that the reflectance curve with the least sum of slope squared (or in the continuous case, the integral of the squared first derivative) will do this. After presenting the algorithms, I examine the quality of the computed reflectance curves compared to thousands of documented reflectance curves measured from paints and pigments available commercially or in nature. Being able to generate reflectance curves from three-dimensional color information is useful in computer graphics, particularly when modeling color transformations that are wavelength specific.

#### Introduction

There are many different 3D color space models, such as XYZ, RGB, HSV, L*a*b*, etc., and one thing they all have in common is that they require only three quantities to describe a unique color sensation. This reflects the “trichromatic” nature of human color perception. The space of *color stimuli*, however, is not three dimensional. To specify a unique color stimulus that enters the eye, the power level at every wavelength over the visible range (e.g., 380 nm to 730 nm) must be specified. Numerically, this is accomplished by discretizing the spectrum into narrow wavelength bands (e.g., 10 nm bands), and specifying the total power in each band. In the case of 10 nm bands between 380 and 730 nm, the space of color stimuli is 36 dimensional. As a result, there are many different color stimuli that give rise to the same color sensation (infinitely many, in fact).

For most color-related applications, the three-dimensional representation of color is efficient and appropriate. But it is sometimes necessary to have the full wavelength-based description of a color, for example, when modeling color transformations that are wavelength specific, such as dispersion or scattering of light, or the subtractive mixture of colors, for example, when mixing paints or illuminating colored objects with various illuminants. In fact, this web page was developed in support of another page on this site, concerning how to compute the RGB color produced by subtractive mixture of two RGB colors.

I present several algorithms for converting a three-dimensional color specifier (sRGB) into a wavelength-based color specifier, expressed in the form of a reflectance curve. When quantifying object colors, the reflectance curve describes the fraction of light that is reflected from the object by wavelength, across the visible spectrum. This provides a convenient, dimensionless color specification, a curve that varies between zero and one (although fluorescent objects can have reflectance values >1). The motivating idea behind these algorithms is that the one reflectance curve that has the least sum of slope squared (integral of the first derivative squared, in the continuous case) seems to match reasonably well the reflectance curves measured from real paints and pigments available commercially and in nature. After presenting the algorithms, I compare the computed reflectance curves to thousands of documented reflectance measurements of paints and pigments to demonstrate the quality of the match.

#### sRGB Triplet from a Reflectance Curve

The reverse process of computing an sRGB triplet from a reflectance curve is straightforward. It requires two main components: (1) a mathematical model of the “standard observer,” which is an empirical mathematical relationship between color stimuli and color sensation (tristimulus values), and (2) a definition of the sRGB color space that specifies the reference illuminant and the mathematical transformation between tristimulus values and sRGB values.

The linear transformation relating a color stimulus, , to its corresponding tristimulus values, , is

The column vector has three elements, , , and . The matrix has three rows (called the three “color matching functions”) and columns, where is the number of discretized wavelength bands. In this study, all computations are performed with 36 wavelength bands of width 10 nm, running over the range 380 nm to 730 nm. The stimulus vector also has components, representing the total power of all wavelengths within each band. The specific color matching functions I use in this work are the CIE 1931 color matching functions. (Note that I’m using the symbol here; the standard matrix is x 3, and so indicates that it has been transposed.)

The stimulus vector can be constructed as the product of an diagonal illuminant matrix, , and an reflectance vector, . The computation of is usually normalized so the tristimulus value is equal to 1 when is a perfect reflector (contains all 1s). The normalizing factor, , is thus the inner product of the second row of and the illuminant vector, , yielding the alternate form of the tristimulus value equation

The transformation from tristimulus values to sRGB is a two-step process. First, a 33 linear transformation, , is applied to convert to , which is a triplet of “linear RGB” values:

The second step is to apply a “gamma correction” to the linear values, also known as “companding” or applying the “color component transfer function.” This is how it is done: for each , , and component of , let’s generically call it , if , use , otherwise use . This gives sRGB values in the range of 0 to 1. To convert them to the alternate integer range of 0 to 255 (used in most 24-bit color devices), we multiply each by 255 and round to the nearest integer.

**sRGB to rgb**

`sRGB=sRGB/255; % convert 0-255 range to 0-1 range`

for i=1:3

if sRGB(i)<0.04045

rgb(i)=sRGB(i)/12.92;

else

rgb(i)=((sRGB(i)+0.055)/1.055)^2.4;

end

end

The expression relating and above can be simplified by combining the three matrices and the normalizing factor into a single matrix, so that

The formal definition of the sRGB color space uses an illuminant similar to daylight, called D65, as its “reference” illuminant. Here are the specific values for the matrix, the matrix, the D65 vector, and the matrix. The normalizing factor, , has a value of 10.5677. Most of the RGB-related theory I present here I learned from Bruce Lindbloom’s highly informative website.

Now that we have a simple expression for computing sRGB from a reflectance curve, we can use that as the basis of doing the opposite, computing a reflectance curve from an sRGB triplet. In the sections that follow, I will present five different algorithms for doing this. Each has its strengths and weaknesses. Once they are presented, I will then compare them to each other and to reflectance curves found in nature.

#### Linear Least Squares (LLS) Method

Since there are so many more columns of than rows, the linear system is under-determined and gives rise to an dimensional subspace (33-dimensional in our case) of reflectance curves for a single sRGB triplet. There are well-established techniques for solving under-determined linear systems. The most common method goes by various names: the linear least squares method, the pseudo-inverse, the least-squares inverse, the Moore-Penrose inverse, and even the “Fundamental Color Space” fundamental metamer.

Suppose we pose this optimization problem:

This linearly constrained minimization can be solved easily by forming the Lagrangian function

The solution can be found by finding a stationary point of , i.e., setting partial derivatives with respect to and equal to zero:

Solving this system by eliminating gives the LLS solution

Thus, a reflectance curve can be found from a sRGB triplet by simply converting it to linear and multiplying it by a 33 matrix, . Unfortunately, the resulting solution is sometimes not very useful. Consider its application to this reflectance curve, which represents a bright red object color:

The LLS solution contains negative reflectance values, which don’t have physical meaning and limit its usefulness in realistic applications. Computationally, this is a very efficient method. The matrix can be computed in advance, as shown here, and each new sRGB value needs only be multiplied by it to get a reflectance curve.

#### Least Slope Squared (LSS) Method

Note that the standard LLS method minimizes , that is, it finds the solution that is nearest to the origin, or the reflectance curve that oscillates most tightly about the wavelength axis. For purposes of computing reflectance curves, I can’t think of a compelling reason why this should be a useful objective.

It dawned on me that it might be better to try a different objective function. The reflectance curves of most natural colored objects don’t tend to oscillate up and down very much. I came up with the idea to minimize the square of the slope of the reflectance curve, summed over the entire curve. In the continuous case, this would be equivalent to

The square is used because it equally penalizes upward and downward movement of . This objective will favor flatter reflectance curves and avoid curves that have a lot of up and down movement. (Edit: I later learned that this objective function has been previously investigated by C van Trigt in 1990 in a series of papers called “Smoothest Reflectance Functions.”)

Other researchers have developed methods for reconstructing reflectance curves from tristimulus values. In order to reduce the oscillations, they typically introduce basis functions that oscillate very little, such as segments of low-frequency sinusoids, or they “frequency limit” or “band limit” the solution by constraining portions of the Fourier transform of the reflectance curve. To my mind, these approaches seem to ignore that fact that realistic reflectance curves can sometimes exhibit sudden steep changes in reflectance at certain frequencies, which would have relatively large high-frequency Fourier components. These methods would not be able to create such reflectance curves.

My proposed method would be able to create steep changes in reflectance, but only as a last resort when flatter-sloped curves are not able to match the target tristimulus values. The other advantage of the minimum slope squared approach is that it can be expressed as a quadratic objective function subject to linear constraints, which is solvable by standard least-square strategies.

Consider this optimization formulation:

where is the number of discrete wavelength bands (36 in this study). This optimization can be solved by solving the system of linear equations arising from the Lagrangian stationary conditions

where is a 3636 tridiagonal matrix

Since and do not depend on , the matrix can be inverted ahead of time, instead of each time an sRGB value is processed. Defining

we have

or

where is the upper-right 363 portion of the matrix. Alternatively, the matrix inversion leading to can be done explicitly, yielding

This 36×3 matrix is shown here.

Since computing is a simple matter of matrix multiplication, the LSS method is just as computationally efficient as the LLS method, and tends to give much better reflectance curves, as I’ll demonstrate later.

Here is a Matlab program for the LSS (Least Slope Squared) method. It also works in the open source free alternative to Matlab, called Octave.

#### Least Log Slope Squared (LLSS) Method

It is generally not a good idea to allow reflectance curves with negative values. Not only is this physically meaningless, but it also can cause problems down the road when the reflectance curves are used in other computations. For example, when modeling subtractive color mixture, it may be necessary to require the reflectance curves to be strictly positive.

One way to modify the LSS method to keep reflectances positive is to operate the algorithm in the space of the logarithm of reflectance values, . I call this the **Least Log Slope Squared (LLSS) method**:

This new optimization is not as easy to solve as the previous one. Nevertheless, the Lagrangian formulation can still be used, giving rise to a system of 39 *nonlinear* equations and 39 unknowns:

where is the same 3636 tridiagonal matrix presented earlier.

Newton’s method solves this system of equations with ease, typically in just a few iterations. Forming the Jacobian matrix,

the change in the variables with each Newton iteration is found by solving the linear system

Here is a Matlab program for the LLSS (Least Log Slope Squared) method. It also works in the free Octave software. I added a check for the special case of sRGB = (0,0,0), which simply returns . This is necessary since the log formulation is not able to create a reflectance of exactly zero (nor is that desirable in some applications). It can come very close to zero as approaches , but it is numerically better to handle this one case specially. I chose the value of 0.0001 because it is the largest power of ten that translates back to an integer sRGB triplet of (0,0,0).

If you prefer not to work in a programming environment, here is an Excel spreadsheet that implements the LLSS method. Just follow the instructions at the top of the spreadsheet. It uses the “Solver” add-in within Excel, which you may need to load. Here are instructions for doing that.

The LLSS method requires substantially more computational effort than the previous two methods. Each iteration of Newton’s method requires the solution of 39 linear equations in 39 unknowns.

This Matlab program has also been tested in Octave and was found to work fine.

#### Iterative Least Log Slope Squared (ILLSS) Method

The LLSS method above can return reflectance curves with values >1. Although this is physically meaningful phenomenon (fluorescent objects can exhibit this), it may be desirable in some applications to have the entire reflectance curve between 0 and 1. It dawned on me that I might be able to modify the Lagrangian formulation to cap the reflectance values at 1. The main obstacle to doing this is that the use of inequality constraints in the Lagrangian approach greatly complicates the solution process, requiring the solution of the “KKT conditions,” and in particular, the myriad “complementary slackness” conditions. If only there were some way to know which reflectance values need to be constrained at 1, then these could be treated by a set of equality constraints and no KKT solution would be necessary.

That led me to investigate the nature of the LLSS reflectance curves with values >1. I ran the LLSS routine on every value of sRGB by intervals of five, that is, sRGB = (0,0,0), (0,0,5), (0,0,10), …, (255,255,250), (255,255,255). In every one of those 140,608 cases, the algorithm found a solution in less than a dozen or so iterations (usually just a handful), and 38,445 (27.3%) of them had reflectance values >1.

Of the 38,445 solutions with values >1, 36,032 of them had a single contiguous region of reflectance values >1. The remaining 2,413 had two regions, always located at both ends of the visible spectrum. Since the distribution of values >1 is so well defined, I started thinking of an algorithm that would iteratively force the reflectance to 1. It would start by running the LLSS method. If any of the reflectance values ended up >1, I would add a single equality constraint forcing the reflectance at the maximum of each contiguous >1 region to equal 1, solve that optimization, then force the adjacent values that were still >1 to 1, optimize again, and repeat until all values were 1.

That was getting to be an algorithmic headache to implement, so I tried a simpler approach, as follows. First, run the LLSS method. If any reflectance values end up >1, constrain ALL of them to equal 1, and re-solve the optimization. This will usually cause some more values adjacent to the old contiguous regions to become >1, so constrain them in addition to the previous ones. Re-solve the optimization. Repeat the last two steps until all values are 1. Here is an animation of this process, which I call the **Iterative Least Log Slope Squared** (ILLSS) process, applied to sRGB = (75, 255, 255):

To express the ILLSS algorithm mathematically, let’s begin with the LLSS optimization statement and add the additional equality constraints:

“FixedSet” is the set of reflectance indices that are constrained to equal 1, or equivalently, the set of indices constrained to equal zero (since ). Initially, FixedSet is set to be the empty set. Each time the optimization is repeated, the values 0 have their indices added to this set.

We can define a matrix that summarizes the fixed set, for example:

This example indicates that there are two reflectance values being constrained (because has two rows), and the third and fifth reflectance values are the particular ones being constrained.

The Lagrangian formulation now has additional Lagrange multipliers, called , one for each of the constrained reflectances. The system of nonlinear equations produced by finding a stationary point of the Lagrangian (setting partial derivatives of the Lagrangian with respect to each set of variables (, , and ) equal to zero) is

where is the same 3636 tridiagonal matrix presented earlier. As before, we solve this nonlinear system with Newton’s method. Forming the Jacobian matrix,

the change in the variables with each Newton iteration is found by solving the linear system

Here is a Matlab program that performs the ILLSS (Iterative Least Log Slope Squared) optimization. I included a check for the two special cases of = (0,0,0) or (255,255,255), which simply return = (0.0001, 0.0001, …, 0.0001) or (1, 1, …, 1). The additional special case of (255,255,255) is needed because numerical issues arise if the matrix grows to 3636, as it would in that second special case. This program works in Octave as well.

**Optimization Packages**

I am aware of two commercial optimization packages that work. First, Excel has an Add-In called the Solver. I believe it comes with every version of Excel, but you may have to enable it to get it to appear in the Data tab/menu. One of its solving methods is called “GRG Nonlinear” and it works very well on this problem. Here is an Excel spreadsheet with the optimization set up. Just follow the instructions on the sheet to get a physically valid reflectance curve from a specified sRGB triplet.

The second optimization package I recommend is the Optimization Toolbox in Matlab. I don’t know whether or not the free alternative to Matlab (called Octave) has similar optimization capabilities. Here is a Matlab function that does the inequality-constrained optimization.

I consider both of these optimization packages to be “heavy artillery” for solving this relatively simple optimization problem, in comparison to the tight and efficient ILLSS solution method.

#### Iterative Least Slope Squared (ILSS) Method

For completeness, I thought it would be a good idea to add one more algorithm. Recall the ILLSS method modifies the LLSS method to cap reflectances >1. Similarly, the ILSS method will modify the LSS method to cap values both >1 and <0. The ILSS may reduce computational effort in comparison to the ILLSS method since the inner loop of the ILLSS method requires an iterative Newton’s method solution, whereas there would be no inner loop needed with the ILSS method; it is simply the solution of a linear system of equations. Here is the ILSS formulation:

is the set of reflectance indices that are constrained to equal 1, and is the set of indices that are constrained to equal (the smallest allowable reflectance, typically 0.00001). Initially, both fixed sets are the empty set. Each time the optimization is repeated, the values have their indices added to and those have their indices added to .

We define two matrices that summarize the fixed sets, for example:

This example indicates that there are two reflectance values being constrained to equal 1 (because has two rows), and the third and fifth reflectance values are the particular ones being constrained. There are three reflectance values being constrained to (because has three rows), and the first, sixth, and fourth are the particular ones being constrained. The order in which the rows appear in these matrices is not important.

At each iteration of the ILSS method, this linear system is solved:

where and are Lagrange multipliers associated with the and sets, respectively.

This system can be solved for , yielding the expression

where

and are the upper-left 3636 and upper-right 363 parts of , respectively. They can be computed ahead of time. Note that only an matrix needs to be inverted at each iteration, where is the number of values being held fixed, typically zero or a small number. When is zero, the ILSS method simplifies to the LSS method. Here is a Matlab program that performs the ILSS (Iterative Least Slope Squared) method. It also works in Octave. Matlab code can be dense and difficult to read, so I’ve prepared a line-by-line commentary for understanding the code.

#### Comparison of Methods

*(Update 1/9/20: The newly-developed LHTSS method presented below is a better choice than ILLSS when requiring the reflectance curves to fall strictly between 0 and 1. It matches the Munsell colors better and takes about 60% of the computational effort of ILLSS. See the 6/4/2019 update below for more information.)*

I ran each of the five algorithms with every sRGB value (by fives) and have summarized the results in the table below. The total number of runs for each solution was 140,608.

Name |
ExecutionTime (sec) |
Max |
Min |
Num>1 |
Num<0 |
MaxIter. |
MeanIter. |
ComputationalEffort |

LLS, Linear Least Squares |
2.7 | 1.36 | -0.28 | 26,317 | 50,337 | n/a | n/a | Matrix mult only |

LSS, Least Slope Squared |
2.7 | 1.17 | -0.17 | 9,316 | 48,164 | n/a | n/a | Matrix mult only |

ILSS, Iterative Least Slope Squared |
25.7 | 1 | 0 | 0 | 0 | 5 | 1.49 | Mult soln of linear eqns** |

LLSS, Least Log Slope Squared |
322. | 3.09 | 0 | 38,445 | 0 | 16 | 6.77 | Mult soln of 39 linear eqns |

ILLSS, Iterative Least Log Slope Squared |
525. | 1 | 0 | 0 | 0 | 5* | 1.41* | Mult soln of (39+) linear eqns** |

* This is outer loop iteration count. The inner loop has iteration count similar to previous line.

** The quantity is the number of reflectance values that end up being constrained at either 1 or 0.

**Explanation of Columns**

- Execution Time: Real-time duration to compute 140,608 reflectance curves of all sRGB values (in intervals of five), on a relatively slow Thinkpad X61 tablet. Relative times are more important than absolute times.
- Max : The maximum reflectance value of all computed curves.
- Min : The minimum reflectance value of all computed curves. Zero is used to represent some small specified lower bound, typically 0.0001.
- Num >1: The number of reflectance curves with maxima above 1.
- Num <0: The number of reflectance curves with minima below 0.
- Max Iter.: The largest number of iterations required for any reflectance curve (for iterative methods).
- Mean Iter.: The mean value of all iteration counts (for iterative methods).
- Computational Effort: Comments on the type of computation required for the method.

Clearly, the methods that don’t need to solve linear systems of equations are much faster. The reflectance curves they create, however, may not be very realistic.

#### Comparison of Reflectance Curves

In this section, I compare the computed reflectance curves against two large sets of measured reflectance data. My goal is to assess how “realistic” the computed reflectance curves are in comparison to reflectances measured from real colored objects.

**Munsell Color Samples**

The first set comes from the Munsell Book of Colors, specifically the 2007 glossy version. The Munsell system organizes colors by hue, chroma (like saturation), and value (like lightness). In 2012, Paul Centore measured 1485 different Munsell samples and published the results here. I computed the sRGB values of each sample, and used those quantities to compute reflectance curves for each of the five methods described above. Of the 1485 samples, 189 of them are outside the sRGB gamut (have values <0 or >255) and were not examined in this study, leaving 1296 in-gamut samples. An Excel file of the 1485 reflectance curves and sRGB values is available here.

When comparing reflectance curves, some parts of the curve are more important than others. As we approach both ends of the visible spectrum, adjacent to the UV and IR ranges, the human eye’s sensitivity to these wavelengths rapidly decreases. This phenomenon is described by the “luminous efficiency” curve shown in the adjacent figure. As I was computing the reflectance curves, I noticed that there are often large discrepancies between the computed and measured reflectance curves at the ends of the visible spectrum. These large differences have relatively little impact on color perception, and would also have little consequence on operations performed on the computed reflectance curves, such as when modeling subtractive color mixture. Consequently, I decided to downplay these end differences when comparing the curves.

To assess how well each computed reflectance matches the measured reflectance curve, I first subtracted one from the other and took absolute values of the difference. Then I multiplied the differences by the luminous efficiency curve. Finally, I summed all of the values to give a single measure of how well the reflectance curves match, or a “reflectance match measure” ():

Lower values of represent a better match, with zero being a perfect match of the two reflectance curves. Keep in mind that regardless of the value of , the sRGB triplets of the computed and measured curves always match exactly, as would the perceived color evoked by the two reflectance curves.

I examined the maximum and mean values of for each of the five methods when applied to all 1296 Munsell samples:

Name |
Max |
Mean |

LLS, Linear Least Squares |
2.46 | 0.88 |

LSS, Least Slope Squared |
1.11 | 0.17 |

ILSS, Iterative Least Slope Squared |
1.04 | 0.16 |

LLSS, Least Log Slope Squared |
0.92 | 0.15 |

ILLSS, Iterative Least Log Slope Squared |
0.86 | 0.15 |

The Linear Least Squares method is clearly much worse than the others. Considering that the Least Slope Squared method (LSS) requires the same computational effort as LLS, but with much better results, I decided to present the comparison figures that follow for the last four methods only: LSS, ILSS, LLSS, and ILLSS.

The following figures compare the four methods. Each gray square represents one Munsell sample, and the shade of gray indicates the corresponding value. The shade of gray is linearly mapped so that it is white for = 0 and black for = 1.11, the maximum value for all four methods.

Overall, the log version of the algorithms do a little better. The worst matches with non-log-based LSS/ILSS take place with the highly chromatic reds and purples, and sometimes with the bright yellows. With log-based LLSS/ILLSS, the worst matches tend to be in the bright yellows, and sometimes in the chromatic reds and oranges. In both pairs, the iterative version clips the reflectance curves between 0 and 1, and usually reduces the mismatch somewhat, but not by very much.

Here are some examples of LSS and ILSS with high s:

This is typical of the mismatch with chromatic reds and purples. This one (Hue=5RP, Value=6, Chroma=12) has an of 1.04. The curve for LSS is not visible because it is directly behind the ILSS curve. There was no clipping needed, so the two curves are the same.

One weakness of the LSS method is that it tends to drop into the negative region for chromatic red colors. In these cases, the clipping provided by the ILSS method greatly improves the match. The following example has an LSS of 1.01 and an improved ILSS of 0.40:

The following is a typical mismatch for LSS/ILSS in the yellow region:

Generally, when LSS/ILSS has a bad match, it is because it oscillates too little.

Here are some examples of log-based LLSS/ILLSS with high s:

This example has an of 0.88, mainly because it takes advantage of the wavelengths of light that give rise to a yellow sensation, in contrast to how the pigments in this sample get the same yellow sensation from a broader mix of wavelengths.

On bright red and orange Munsell samples, the LLSS/ILLSS curves tend to shoot high on the red side:

Even though there is considerable clipping with the ILLSS version, most of it takes place in the long wavelengths, which does not help the score as much. This is evident in how little ILLSS needs to adjust the mid-wavelength range to make up for the huge difference in the long wavelengths, while maintaining the same sRGB values.

In summary, the non-log versions (LSS/ILSS) tend to undershoot peaks and the log versions (LLSS/ILLSS) tend to overshoot them. Overall, the log versions give a somewhat better match, but at the expense of considerably more computation.

**Commercial Paints and Pigments**

The second large dataset of reflectance curves comes from Zsolt Kovacs-Vajna’s RS2color webpage. He has a database of reflectance curves for many sets of commercial paints and pigments, which can be obtained by emailing a request to him. They are grouped into the following families:

1. apaFerrario_PenColor

2. Chroma_AtelierInteractive

3. ETAC_EFX500b

4. ETAC_EFX500t

5. GamblinConservationColors

6. Golden_HB

7. Golden_OpenAcrylics

8. GretagMacbethMini

9. Holbein_Aeroflash

10. Holbein_DuoP

11. KremerHistorical

12. Liquitex_HB

13. Maimeri_Acqua

14. Maimeri_Brera

15. Maimeri_Polycolor

16. MunsellGlossy5B

17. MunsellGlossy5G

18. MunsellGlossy5P

19. MunsellGlossy5R

20. MunsellGlossy5Y

21. MunsellGlossyN

22. pigments

23. supports

24. Talens_Ecoline

25. Talens_Rembrandt

26. Talens_VanGoghH2Oil

27. whites

28. WinsorNewton_ArtAcryl

29. WinsorNewton_Artisan

30. WinsorNewton_Finity

31. WinsorNewtonHandbook

In total, there are 1493 different samples in these 31 groups. I discarded the 275 of them that were out of the sRGB gamut, and computed reflectance curves from the sRGB values of the remaining 1218 samples. I again computed values to compare the computed curves to the actual measured curves. The following animated GIFs show the values (color coded according to the colormap on the right) plotted in sRGB space.

Above, values (by color) for the LSS method plotted in sRGB space.

Above, values (by color) for the ILSS method plotted in sRGB space.

Above, values (by color) for the LLSS method plotted in sRGB space.

Above, values (by color) for the ILLSS method plotted in sRGB space.

It is apparent from these GIFs that the non-log-based methods have a fairly large region of mismatch in the red region (large R, small G and B). The log-based versions have a much smaller region of mismatch in the yellow region (large R and G, small B). The iterative version of both methods improve the matches in both cases. The relative sizes of the mismatch regions can be seen more clearly in a projection onto the R-G plane:

Above, values (by color) for the LSS method projected onto the sRGB R-G plane.

Above, values (by color) for the ILSS method projected onto the sRGB R-G plane.

Above, values (by color) for the LLSS method projected onto the sRGB R-G plane.

Above, values (by color) for the ILLSS method projected onto the sRGB R-G plane.

#### Conclusions

*(Update 1/9/20: The newly-developed LHTSS method presented below is a better choice than ILLSS when requiring the reflectance curves to fall strictly between 0 and 1. It matches the Munsell colors better and takes about 60% of the computational effort of ILLSS. See the 6/4/2019 update below for more information.)*

In summary, the method that best matches paint and pigment colors found commercially and in nature is the ILLSS (Iterative Least Log Slope Squared) method. It suffers, however, from very large computational requirements. If efficiency is more important, then the ILSS (Iterative Least Slope Squared) method is the preferred one. A third alternative that is midway between the match quality and computing effort is LLSS (Least Log Slope Squared), but be aware that some reflectance curves can end up with values >1. I would not recommend the LSS (Least Slope Squared) method, despite its spectacular computational efficiency, because it can give reflectance curves with physically meaningless negative values. The following table summarizes the three suggested methods:

Algorithm Name |
Computational Effort |
Comments |
Link to Matlab Code |

ILSS (Iterative Least Slope Squared) |
Relatively little. | Very fast, but tends to undershoot reflectance curve peaks, especially for bright red and purple colors. Always returns reflectance values in the range 0-1. | link |

LLSS (Least Log Slope Squared) |
About 12 times that of ILSS. | Better quality matches overall, but tends to overshoot peaks in the yellow region. Some reflectance values can be >1, especially for bright red colors. | link |

ILLSS (Iterative Least Log Slope Squared) |
About 20 times that of ILSS. | Best quality matches. Tends to overshoot peaks in the yellow region. Always returns reflectance values in the range 0-1. | link |

__________________________________

**Update 9/26/2018:** I’ve added a new page that presents an even faster way to generate reflectance curves from sRGB values. It is as fast as the fastest methods presented above, guarantees reflectance values in the 0-1 range, but produces reflectance curves that are not as smooth or realistic as the slower, higher quality reconstruction methods. Nevertheless, it permits good quality subtractive mixture with little computational overhead. (This new method is the result of the idea being suggested to me by Brien Dieterle, who asked me to help find a way to implement the idea.)

__________________________________

**Update 10/2/2018:** I have recently added a new page that shows how to apply these methods to the color space called Recommendation ITU-R BT.2020-2 (10/2015), or Rec. 2020 for short.

________________________________________

**Update 5/23/2019:** I’ve recently applied the reflectance reconstruction ideas presented here to the subject of chromatic adaptation, and have written a paper called “Chromatic Adaptation Transform by Spectral Reconstruction,” which was published in the journal *Color Research and Application* (citation: Burns SA. Chromatic adaptation transform by spectral reconstruction. Color Res Appl. 2019;44(5):682-693). Here is the full text of the published paper courtesy of the Wiley Content Sharing initiative.

________________________________________

**Update 6/4/2019:** I’ve just developed a new method that gives even better natural-reflectance-matching results than the ILLSS method above, with considerably less computational effort required (comparable to that of LLSS). Here it is:

#### Least Hyperbolic Tangent Slope Squared (LHTSS) Method

Recall that the Least Log Slope Squared (LLSS) method generates reflectance curves that are strictly positive, without requiring explicit bounds in the optimization statement. The lower bound is handled implicitly in the logarithmic transformation. This section presents a newly developed approach to generate reflectances that are strictly within the 0-1 range, also not needing explicit bounds. Why is this useful? The methods that use explicit bounds (presented above: ILSS and ILLSS) tend to have abrupt discontinuities in slope when a bounding constraint is engaged. This makes the reflectance curves seem somewhat unnatural. The LHTSS method described in this section makes reflectance curves that gently approach the upper and/or lower bounds, and match the reflectance curves of commercial paints and pigments even better than any of the previous methods above.

The **Least Hyperbolic Tangent Slope Squared (LHTSS) method** uses a transformation defined by

to keep all reflectance values between 0 and 1. Below is a plot of the hyperbolic tangent function:

Note how the curve approaches +1/-1 as grows large in either the positive or negative directions. By adding 1 to tanh and then dividing by 2, the upper and lower bounds shift to 1 and 0. This causes the reflectances to remain strictly between 0 and 1 for any value of :

The least slope optimization statement for the Hyperbolic Tangent Slope Squared method is:

Noting that the first derivative of is , the stationary conditions of the Lagrangian function associated with this optimization comprise a system of 39 nonlinear equations and 39 unknowns:

where is the 3636 tridiagonal matrix of finite differencing constants presented earlier.

Newton’s method solves this system of equations with ease, typically in a fairly small number iterations. Noting that the first derivative of is , the Jacobian matrix of first partial derivatives of is

The change in the variables with each Newton iteration is found by solving the linear system

At each iteration, the values of and are updated using and .

Here is a Matlab program for the LHTSS (Least Hyperbolic Tangent Slope Squared) method. It also works in the free Octave software. It handles the two cases of sRGB = (0,0,0) and (255,255,255) specially, since the reflectance curves associated with them fall outside the strict interval (0->1) imposed by the hyperbolic tangent transform.

If you prefer not to work in a programming environment, here is an Excel spreadsheet that implements the LHTSS method. Just follow the instructions at the top of the spreadsheet. It uses the “Solver” add-in within Excel, which you may need to load. Here are instructions for doing that.

An example of how this hyperbolic tangent method (LHTSS) compares to the logarithmic method (LLSS) and the regular least slope squared method (LSS), when applied to the Munsell color 7.5R 5/16, is shown in the figure below.

This figure demonstrates an instance where LSS has negative elements, LLSS has elements >1, and LHTSS corrects both of these deficiencies. It is quite striking how well the LHTSS reconstruction matches the measured reflectance curve of this particular Munsell chip.

To check just how well the curves produced by LHTSS compare to the reflectance curves measured from 1296 Munsell chip colors, the “reflectance match measure, RMM” was computed and added to the previous table:

Name |
Max |
Mean |

LLS, Linear Least Squares |
2.46 | 0.88 |

LSS, Least Slope Squared |
1.11 | 0.17 |

ILSS, Iterative Least Slope Squared |
1.04 | 0.16 |

LLSS, Least Log Slope Squared |
0.92 | 0.15 |

ILLSS, Iterative Least Log Slope Squared |
0.86 | 0.15 |

LHTSS, Least Hyperbolic Tangent Slope Squared |
0.84 | 0.14 |

It is evident that LHTSS performs the best of all methods with regard to the comparison to Munsell reflectance curves.

The computational effort required for the LHTSS method is comparable to that of LLSS. The computational effort study performed above, using 140,608 sRGB triplets (every value possible in intervals of five), showed that the LLSS method required 6.77 iterations on average. The LHTSS method requires an average of 5.66 iterations, so its average run times are somewhat smaller (approx 5%).

________________________________________

**Update 1/18/20:** These methods (including the new LHTSS method) have been published in a paper called “Numerical Methods for Smoothest Reflectance Reconstruction,” in the journal *Color Research and Application* (citation: Burns SA. Numerical methods for smoothest reflectance reconstruction. Color Res Appl. 2020;45(1):8-21).

________________________________________

#### Acknowledgments

This work has been a tremendous learning experience for me, and I want to thank several people who graciously posted high quality material, upon which most of my work has been based: Bruce Lindbloom, Bruce MacEvoy, Zsolt Kovacs-Vajna, David Briggs, and Paul Centore. If you are interested in learning more about color theory, these five links are excellent places to start!

________________________________________

**Generating Reflectance Curves from sRGB Triplets**by Scott Allen Burns is licensed under a Creative Commons Attribution-ShareAlike 4.0 International License.

________________________________________

It would be awesome if you could provide a bit more detail on using other colorspaces and light sources. I can’t quite figure out how I would modify the program to use Adobe RGB, for instance, which uses a normal 2.4 gamma instead of the funky sRGB. Ultimately it should be possible to use an ICC profile to get the RGB working space. Then you could paint in a space that exactly matches your calibrated/profiled monitor.

I do mention that at the end of the companion article on subtractive mixture. It’s really pretty simple. This website:

http://www.brucelindbloom.com/index.html?Eqn_RGB_XYZ_Matrix.html

outlines the changes you would make to the M matrix for AdobeRGB and other RGB types. This would then be used to compute a new T matrix. See the “more details” drop-downs on that page for, uh, more details!

You’re right, it’s all there. I was looking for the w normalizing factor on Bruce’s site, totally missed your description of how that is calculated. Sorry about that!

I’m having trouble wrapping my brain around alpha channel. Do you know if there is a fundamental issue with dealing with transparent paints using this technique? It seems that when I ignore the alpha channel blend with pure alpha areas on the canvas, then it’s behaving like I’m blending with white paint. If I premulitply the RGB triplet with the alpha value, then it’s like mixing with a black or darker paint. I’m wondering if the whole concept of alpha channel only works for light/additive model?

I really don’t know much about alpha channel blending. Out of curiosity I found this page that might answer some questions:

https://en.wikipedia.org/wiki/Alpha_compositing

Darn, but thanks. I’ve been looking at that and a few other sites and haven’t figured it out yet, but I’m almost sure it will need to use a different technique. If any color is totally or partially transparent, maybe I can use those opacity values within the WGM mixing. So, if I mix two paints at a 1:1 ratio, but paint A is 100% opaque and paint B is 50% opaque, then the WGM ratio should be modified to be 3:1 and the overall alpha level of 75% would be applied to the resulting RGB. Maybe. . . I’ll give it a shot and let you know.

I meant, 2:1, not 3:1…

I think I figured it out, seems to work really well. Basically the idea is that if a paint is 50% opaque, we can treat it as 50% less paint and adjust the ratio. I’m sure that’s a crude approximation to actual diluted paint, but it gets the job done maybe. Here’s a spreadsheet with the formula I came up with:

https://docs.google.com/spreadsheets/d/1xOPfDzAA7kzHRaxPX05GFRky19mXdQFLD5CDDtCcioU/edit?usp=sharing

I’m sure it can be simplified. I just had to add some checks to avoid div/0

Hi Scott. This recent publication seems to confirm the method I arrived at for handling alpha. Basically treating alpha as the ratio for the weighted geometric mean, scaled to sum to 1.0:

https://hal.archives-ouvertes.fr/hal-01700606/document

Also, I’ve made some more progress that makes your algorithm much more convenient. Instead of calculating a new SPD for each individual color, I’ve just generated 3 spectral primary curves for R, G and B. Then whenever I mix two colors I upscale the RGB to the spectral curves, sum those 3 curves into one unique spectral curve for each color. Then I mix as you’ve demonstrated with the weighted geometric mean and convert back to RGB. It seems to work really well! I made a thread here where you can download the program too:

https://community.mypaint.org/t/smudge-tweaks-testers-needed/1202/

Incredible work on your part. Really impressive. I’m working on implementing your ILLSS algorithm in JavaScript inside an upcoming color library, so it can be used by web developers.

Thanks for your interest! You might find this latest work to be of interest as well: http://scottburns.us/fast-rgb-to-spectrum-conversion-for-reflectances/

Hi,

For the LHTSS method their seems to be a typo either in the matlab or webpage.

In matlab d0 = (tanh(z)+1/2) (which is rho I believe)

This d0 is then used in the F vector but in on the web page it shows that to be diag(sech(z)^2/2. I believe the webpage is correct

Hi Marnix,

Wow, thanks for catching the inconsistency! The mistake is in the webpage; the Matlab code is right. The second term in the F vector should be T (tanh(z)+1)/2 – rgb. This is simply a restatement of the constraint expression in the optimization. I’ve updated the webpage with the correction. Many thanks for catching the problem!

Scott

Great. A question if I may. If I want to start from XYZ values instead of sRGB then I just remove the M matrix from the T and continue with that “new” T and of course use XYZ values instead of linearized RGB values further in the math.

SP the “new” T would be equal to A’ diag(W)/w

Yes, using XYZ instaed of sRGB is quite straightforward. I’ve emailed you the journal paper that grew out of this webpage as an attachment (Burns SA. Numerical methods for smoothest reflectance reconstruction. Color Res Appl. 2020;45(1):8-21). Please feel free to contact me for further discussion if you wish.

Scott thanks for sharing your computations in a spreadsheet…I was messing around with the LHTSS calculations and I wanted to adapt them to D50 illuminant but I could not understand where or how you derived the CMF coefficients for xbar ybar and zbar at the bottom…is this from a chart or calculation??

Hi Todd, Yes, the conversion from RGB to reflectance can be done with Illuminant D50, but be aware that you will not be working with sRGB anymore, as it is specifically referenced to D65. To start with a D50-referenced RGB, you would have to change not only the illuminant-referenced color matching functions (“CMFs” in cells B49 to AK51 in the spreadsheet), but also the rgb-to-XYZ 3×3 matrix (cells J9 to L11 in the spreadsheet) and the XYZ-to-rgb 3×3 matrix (cells J14 to L16). To change the CMFs matrix to be D50 referenced, you would multiply the un-referenced CMFs (http://scottburns.us/wp-content/uploads/2015/03/CIE-1931-2-deg-CMFs-10-nm.txt) by a square matrix with the D50 illuminant on the main diagonal (called diag(W) above) and then divide by the scalar product of the second row of A’ and D50 (called w above). The rgb-to-XYZ and XYZ-to-rgb matrices would be computed according to Bruce Lindbloom’s webpage using D50 instead of D65 (http://www.brucelindbloom.com/index.html?Eqn_RGB_XYZ_Matrix.html). Regarding the “companding” or “gamma correction” that is part of the rgb-to-sRGB standard (http://www.brucelindbloom.com/index.html?Eqn_RGB_to_XYZ.html), you would need to research whether or not there is a similar conversion standard in the D50-referenced color space that you are using.