MATLAB/Octave and Python Source Code for Reflectance Reconstruction and Chromatic Adaptation

Published January 31, 2020 by Scott Allen Burns; last updated Sept 13, 2020

Change Log

1/31/2020 Original publication.
3/26/2020 Added discussion on reflectances > 1 and a note on Python packages.
9/13/2020 Added note that the CAT functions call the method_2 functions.

Introduction

Over the years, I’ve published several pages on this website concerning color science applications. This page is dedicated to collecting the source code developed for the various topics. Both MATLAB (and the free open-source version Octave) and Python source code versions are provided.

The topics covered are listed below (click on item to jump to topic):

REFLECTANCE RECONSTRUCTION

Reflectance Reconstruction Method 2 (or LLSS: Least Log Slope Squared): This method reconstructs a reflectance curve from tristimulus values (or other three-element color specifier, such as RGB). It guarantees that the reflectance values will be strictly positive, making it suitable for applications such as subtractive color mixture simulation. It can sometimes return reflectance values greater than 1 (which can be useful in some applications). The target tristimulus values must fall within the spectral locus (real colors) for it to work (see UTILITY below).
How can we have reflectance values greater than 1?
Normally reflectance values >1 are associated with fluorescent objects. But there is another very useful way to interpret them. They can be used to yield an emissive source instead of an object reflectance. Recall that the color matching functions are weighted to some reference illluminant, like D65. An emissive source (which an be thought of as the result of an illuminant shining on an object) can be produced by multiplying the reflectance values by the illuminant values. This works even if the reflectance values are >1.
Reflectance Reconstruction Method 3 (or LHTSS: Least Hyperbolic Tangent Slope Squared): This method reconstructs a reflectance curve from tristimulus values (or other three-element color specifier, such as RGB). It guarantees that the reflectance values will be strictly within the range 0 to 1, making it suitable for applications requiring true 0-1 reflectance curves. The tristimulus values must fall within the object color solid (object colors) for it to work (see UTILITY below).

What about method 1?

Method 1 (LSS: Least Slope Squared) is not presented here because it is fairly limited in its application. Even though it is very fast (only requiring a single matrix multiplication), it can produce reflectance curves with some negative components, which is physically impossible. Most useful applications require a strictly non-negative reflectance curve.

ADAPTATION

Chromatic Adaptation: This method addresses the question, “Given a set of tristimulus values referenced to a particular illuminant, how will the tristimulus values change when moving to another illuminant?” This is also known as a “chromatic adaptation transform” or CAT.

What about the symmetric CAT version?

The version of chromatic adaptation presented here is the non-symmetric version. The symmetric version guarantees that if a color is transformed from illuminant A to B, then that transformed color, when transformed back from B to A, will match the original color. The version implemented here will have very slight differences that will likely not have any impact in practical applications. The non-symmetric version was chosen because it is simpler to implement and can use the method 2 code for reflectance reconstruction, which is already presented on this web page. The symmetric version would require a modified version of method 2. Reference [3] presents detailed information on how to modify method 2 for the symmetric case.

UTILITY

Test if XYZ is Within the Spectral Locus: This function will determine if a set of tristimulus values falls within the spectral locus, as required by method 2 above.
Test if XYZ is Within the Object Color Solid: This function will determine if a set of tristimulus values falls within the object color solid, as required by method 3 above.

Note: Methods 2 and 3 refer to methods presented in Reference [1]. Methods “LLSS”, and “LHTSS” refer to the corresponding methods (for sRGB) presented in Reference [2]. The chromatic adaptation code is developed in Reference [3]. The utility functions are adapted from Reference [4].

Notation

The following symbols are used. Python generally prefers that lower case letters are used for variables, so that convention is (mostly) followed. XYZ remains capitalized to distinguish it from chromaticity coordinates xyz or xyY.

Reflectance Reconstruction Methods

Two methods are presented here, corresponding to method 2 (LLSS) and method 3 (LHTSS). Two functions are provided for each method. The first function computes some preparatory arrays that depend only on the illuminant. The second function accepts these arrays and a target XYZ and generates the reconstructed reflectance. The preparatory arrays can be reused for a series of target XYZ values referenced to a common illuminant, which is a common use case.

Reflectance Reconstruction Method 2 (or LLSS Least Log Slope Squared)

MATLAB/Octave

Preparatory function source code.
Reflectance reconstruction source code.

Python

Preparatory function source code.
Reflectance reconstruction source code.

Examples

MATLAB/Octave example.
Show example
The following code demonstrates the method 2 code for the case of XYZ both inside and outside the spectral locus, using an equal energy reference illuminant. The CMFs (A) are hard coded.

% Test the function method_2.m

% Version 2020_01_18 Scott Allen Burns

% create the input arguments:
% CIE 1931 2 degree color matching functions
% 380nm – 730nm, 10nm bins
A = [0.001368,0.000039,0.00645;
0.004243,0.00012,0.02005;
0.01431,0.000396,0.06785;
0.04351,0.00121,0.2074;
0.13438,0.004,0.6456;
0.2839,0.0116,1.3856;
0.34828,0.023,1.74706;
0.3362,0.038,1.77211;
0.2908,0.06,1.6692;
0.19536,0.09098,1.28764;
0.09564,0.13902,0.81295;
0.03201,0.20802,0.46518;
0.0049,0.323,0.272;
0.0093,0.503,0.1582;
0.06327,0.71,0.07825;
0.1655,0.862,0.04216;
0.2904,0.954,0.0203;
0.43345,0.99495,0.00875;
0.5945,0.995,0.0039;
0.7621,0.952,0.0021;
0.9163,0.87,0.00165;
1.0263,0.757,0.0011;
1.0622,0.631,0.0008;
1.0026,0.503,0.00034;
0.85445,0.381,0.00019;
0.6424,0.265,0.00005;
0.4479,0.175,0.00002;
0.2835,0.107,0;
0.1649,0.061,0;
0.0874,0.032,0;
0.04677,0.017,0;
0.0227,0.00821,0;
0.011359,0.004102,0;
0.00579,0.002091,0;
0.002899,0.001047,0;
0.00144,0.00052,0];
% number of rows of cmfs
n = size(A,1);
% equal energy illuminant
W = ones(n,1);

% create preparatory arrays d and cmfs_w
[D,Aw] = method_2_prep(A,W);

disp('Test inside spectral locus:');
% target w-referenced tristimulus values
XYZw = [0.2; 0.3; 0.1];
% reflectance reconstruction
rho = method_2(D,Aw,XYZw);
disp(rho');
disp(' ');

disp('Test outside spectral locus:');
XYZw = [0.8; 0.2; 0.5];
% reflectance reconstruction
rho = method_2(D,Aw,XYZw);
disp(rho');

When this code is executed, the following output is produced:

Test inside spectral locus:
Columns 1 through 10
0.0635 0.0635 0.0636 0.0637 0.0642 0.0656 0.0693 0.0765 0.0885 0.1077
Columns 11 through 20
0.1368 0.1791 0.2372 0.3110 0.3904 0.4502 0.4658 0.4333 0.3717 0.3049
Columns 21 through 30
0.2475 0.2039 0.1733 0.1528 0.1396 0.1314 0.1266 0.1239 0.1224 0.1217
Columns 31 through 36
0.1213 0.1211 0.1210 0.1210 0.1210 0.1210

Test outside spectral locus:
Maximum number of iterations reached.
Check to make sure XYZw is within the spectral locus.
Columns 1 through 17
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
Columns 18 through 34
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
Columns 35 through 36
0 0
Python example.
Show example
The following code demonstrates the method 2 code for the case of XYZ both inside and outside the spectral locus, using an equal energy reference illuminant. The CMFs (cmfs) are hard coded.

# Test the function method_2.py

# Version 2020_01_24 Scott Allen Burns

import numpy as np
from method_2_prep import method_2_prep
from method_2 import method_2

# create the input arguments:
# CIE 1931 2 degree color matching functions
# 380nm – 730nm, 10nm bins
cmfs = np.array([[0.001368, 0.000039, 0.00645],
[0.004243, 0.00012, 0.02005],
[0.01431, 0.000396, 0.06785],
[0.04351, 0.00121, 0.2074],
[0.13438, 0.004, 0.6456],
[0.2839, 0.0116, 1.3856],
[0.34828, 0.023, 1.74706],
[0.3362, 0.038, 1.77211],
[0.2908, 0.06, 1.6692],
[0.19536, 0.09098, 1.28764],
[0.09564, 0.13902, 0.81295],
[0.03201, 0.20802, 0.46518],
[0.0049, 0.323, 0.272],
[0.0093, 0.503, 0.1582],
[0.06327, 0.71, 0.07825],
[0.1655, 0.862, 0.04216],
[0.2904, 0.954, 0.0203],
[0.43345, 0.99495, 0.00875],
[0.5945, 0.995, 0.0039],
[0.7621, 0.952, 0.0021],
[0.9163, 0.87, 0.00165],
[1.0263, 0.757, 0.0011],
[1.0622, 0.631, 0.0008],
[1.0026, 0.503, 0.00034],
[0.85445, 0.381, 0.00019],
[0.6424, 0.265, 0.00005],
[0.4479, 0.175, 0.00002],
[0.2835, 0.107, 0],
[0.1649, 0.061, 0],
[0.0874, 0.032, 0],
[0.04677, 0.017, 0],
[0.0227, 0.00821, 0],
[0.011359, 0.004102, 0],
[0.00579, 0.002091, 0],
[0.002899, 0.001047, 0],
[0.00144, 0.00052, 0]])

# number of rows of cmfs
n = cmfs.shape[0] # equal energy illuminant
w = np.ones(n)

# create preparatory arrays d and cmfs_w
d, cmfs_w = method_2_prep(cmfs, w)

print('Test inside spectral locus:')
# target w-referenced tristimulus values
XYZ_w = np.array([0.2, 0.3, 0.1])
# reflectance reconstruction
rho = method_2(d, cmfs_w, XYZ_w)
print(rho)
print(' ')

print('Test outside spectral locus')
XYZ_w = np.array([0.8, 0.2, 0.5])
# reflectance reconstruction
rho = method_2(d, cmfs_w, XYZ_w)
print(rho)

When this code is executed, the following output is produced:

Test inside spectral locus:
[0.06353106 0.06354059 0.06357975 0.06371935 0.06416815 0.06560276
0.0693139 0.07647965 0.08849564 0.10765595 0.13683705 0.1791222
0.23723889 0.31103304 0.39039904 0.45021201 0.46584055 0.43330252
0.37167925 0.3048644 0.24746618 0.20394604 0.17329319 0.15275626
0.13956546 0.13144754 0.12662724 0.12390385 0.12243879 0.1216847
0.1213018 0.12111615 0.12102598 0.12098351 0.12096532 0.12095929]

Test outside spectral locus:
Maximum number of iterations reached.
Check to make sure XYZ_w is within the spectral locus.
[0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.
0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
Process finished with exit code 0

Applications

Application to other color spaces/visual systems/imaging systems.
Expand
The reflectance reconstruction methods are very versatile in their application. To review, the functions above accept tristimulus values referenced to some illuminant along with a matrix ( $A_w$ ) that relates reflectance to tristimulus values (referenced to the same illuminant). Other color spaces can be treated. For example, if working the the RGB color space, a reconstructed reflectance can be generated from a specified RGB triplet as long as the $A_w$ matrix is replaced with a matrix relating reflectance to RGB. Reference [2] demonstrates this for the sRGB color space, where $A_w$ is replaced with a matrix called $T$ .

Another possibile color space is that of cone sensitivity functions (cone fundamentals). The cone fundamentals would be used in place of the $A_w$ matrix and cone signals would be specified instead of XYZ values. If the cone fundamentals are a simple 3×3 linear transformation of the color matching functions (also referred to as satisfying the Luther-Ives condition), then the reconstructed reflectance obtained from a specific set of cone signals will be identical to a reflectance produced from XYZ values related to the cone signals through the same 3×3 transformation.

Even non-human vision systems (e.g., digital cameras) can be treated with these functions. By supplying the three channel sensitivities as CMFs and the three camera signals in place of XYZ, a smoothest reflectance curve can be generated that represents the smoothest reflectance that will give rise to the target camera signals.

Cameras with more than three sensors (multispectral cameras) can also be treated. The code would need some modification, since there will be more than three Lagrange multiplies needed to solve the optimization. But the modifications would be fairly straightforward to implement.
Illuminant reconstruction from white point values.
Expand
In some applications, it is useful to be able to generate an illuminant spectral power distribution (SPD) that has some specified white point. This is done, for example, in the chromatic adaptation method presented below. It is simple to do.

Just pass the standard (non-illuminant-referenced) color matching functions, $A$ , instead of $A_w$ , and pass the illuminant white point tristimlulus values, ${XYZ}^{wp}$ , instead of ${XYZ}_w$ , and method 2 will generate a suitable illuminant curve, $W$ . If the white point has $Y$ =1, then the illuminant will also be appropriately scaled so that an illuminant-referenced $A_w$ can be computed simply as diag(W)*A.
Comparison to Meng’s implementation for smoothest reconstruction.
Expand
Meng et al (Reference [5]) have published a Python implementation of the least slope squared algorithm for reflectance reconstruction as part of the Colour open-source Python package for color science. Their implementation minimizes the sum of squares of reflectance differences (slope), without using the log or tanh transformation, optimized using a general-purpose optimization utility. Without the transformation of variables, the solution will sometimes contain negative components. To prevent this, Meng et al instruct the optimizer to apply explicit lower bounds to keep reflectance $\ge 0$ .

I was curious to learn how much computational overhead is required to use a general-purpose optimization package in comparison to the approach used in methods 2 and 3. I extracted Meng’s code from the Colour package and stripped away any Colour-specific interactions that would only increase its execution time. It is available here.

The following function tests the Meng code (and the method 2 code), and then performs a very crude execution time comparison by running each 1000 times.

# Test meng2015.py

# Version 2020_01_24 Scott Allen Burns

import numpy as np
from method_2_prep import method_2_prep
from method_2 import method_2
from meng2015 import XYZ_to_sd_Meng2015
import time

# create the input arguments:
# CIE 1931 2 degree color matching functions
# 380nm – 730nm, 10nm bins
cmfs = np.array([[0.001368, 0.000039, 0.00645],
[0.004243, 0.00012, 0.02005],
[0.01431, 0.000396, 0.06785],
[0.04351, 0.00121, 0.2074],
[0.13438, 0.004, 0.6456],
[0.2839, 0.0116, 1.3856],
[0.34828, 0.023, 1.74706],
[0.3362, 0.038, 1.77211],
[0.2908, 0.06, 1.6692],
[0.19536, 0.09098, 1.28764],
[0.09564, 0.13902, 0.81295],
[0.03201, 0.20802, 0.46518],
[0.0049, 0.323, 0.272],
[0.0093, 0.503, 0.1582],
[0.06327, 0.71, 0.07825],
[0.1655, 0.862, 0.04216],
[0.2904, 0.954, 0.0203],
[0.43345, 0.99495, 0.00875],
[0.5945, 0.995, 0.0039],
[0.7621, 0.952, 0.0021],
[0.9163, 0.87, 0.00165],
[1.0263, 0.757, 0.0011],
[1.0622, 0.631, 0.0008],
[1.0026, 0.503, 0.00034],
[0.85445, 0.381, 0.00019],
[0.6424, 0.265, 0.00005],
[0.4479, 0.175, 0.00002],
[0.2835, 0.107, 0],
[0.1649, 0.061, 0],
[0.0874, 0.032, 0],
[0.04677, 0.017, 0],
[0.0227, 0.00821, 0],
[0.011359, 0.004102, 0],
[0.00579, 0.002091, 0],
[0.002899, 0.001047, 0],
[0.00144, 0.00052, 0]])

# number of rows of cmfs
n = cmfs.shape[0] # equal energy illuminant
w = np.ones(n)

# create preparatory arrays d and cmfs_w
d, cmfs_w = method_2_prep(cmfs, w)

print('method 2 test inside spectral locus:')
XYZ_w = np.array([0.2, 0.3, 0.1])
rho = method_2(d, cmfs_w, XYZ_w)
print(rho)
print('check:')
print(cmfs_w.T.dot(rho))
print(' ')

print('Meng test inside spectral locus:')
XYZ_w = np.array([0.2, 0.3, 0.1])
rho = XYZ_to_sd_Meng2015(XYZ_w, cmfs_w)
print(rho)
print('check:')
print(cmfs_w.T.dot(rho))
print(' ')

print('method 2 test outside spectral locus:')
XYZ_w = np.array([0.8, 0.2, 0.5])
rho = method_2(d, cmfs_w, XYZ_w)
print(rho)
print('check:')
print(cmfs_w.T.dot(rho))
print(' ')

print('Meng test outside spectral locus:')
XYZ_w = np.array([0.8, 0.2, 0.5])
rho = XYZ_to_sd_Meng2015(XYZ_w, cmfs_w)
print(rho)
print('check:')
print(cmfs_w.T.dot(rho))
print(' ')

# timing study
XYZ_w = np.array([0.2, 0.3, 0.1])
# do method 2
start = time.time()
for _ in range(1000):
rho = method_2(d, cmfs_w, XYZ_w)
end = time.time()
print('elapsed time method 2 (sec)')
print(end – start)
# do meng
start = time.time()
for _ in range(1000):
rho = XYZ_to_sd_Meng2015(XYZ_w, cmfs_w)
end = time.time()
print('elapsed time Meng (sec)')
print(end – start)

The output produced by this function is shown below. After each reflectance is reconstructed, the solution is checked by recomputing the tristimulus values produced by the reflectance solution.

method 2 test inside spectral locus:
[0.06353106 0.06354059 0.06357975 0.06371935 0.06416815 0.06560276
0.0693139 0.07647965 0.08849564 0.10765595 0.13683705 0.1791222
0.23723889 0.31103304 0.39039904 0.45021201 0.46584055 0.43330252
0.37167925 0.3048644 0.24746618 0.20394604 0.17329319 0.15275626
0.13956546 0.13144754 0.12662724 0.12390385 0.12243879 0.1216847
0.1213018 0.12111615 0.12102598 0.12098351 0.12096532 0.12095929] check:
[0.2 0.3 0.1]

Meng test inside spectral locus:
[0.01168061 0.01174062 0.01191253 0.01255591 0.01464514 0.02114432
0.0369829 0.06417914 0.10224493 0.14950563 0.2022331 0.25601156
0.30697802 0.35142211 0.38503908 0.40395434 0.40661952 0.39336736
0.3661289 0.32808212 0.28353645 0.23747087 0.19466611 0.15872304
0.13134669 0.11236025 0.10017805 0.09295798 0.08896425 0.08687811
0.08581701 0.08529703 0.08503886 0.08491505 0.08486197 0.08484061] check:
[0.2 0.3 0.1]
method 2 test outside spectral locus:
Maximum number of iterations reached.
Check to make sure XYZ_w is within the spectral locus.
[0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.
0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.] check:
[0. 0. 0.]
Meng test outside spectral locus:
Optimization failed.
[0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.
0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.] check:
[0. 0. 0.]
elapsed time method 2 (sec)
1.0691401958465576
elapsed time Meng (sec)
30.164045095443726

Process finished with exit code 0

It is evident that the general-purpose optimization function requires a considerable amount of additional computational effort. Repeated tests show that method 2 is about 25-30 times faster.

The timing results are highly dependent upon the specific BLAS (Basic Linear Algebra Subprograms) that is installed with the NumPy package. The default BLAS appears to be about 25 times slower than the fastest BLAS versions available. Here is a discussion on this issue. My installation of NumPy uses the fast MKL version, which can be checked with using np.__config__.show():

import numpy as np
np.__config__.show()
blas_mkl_info:
libraries = [‘mkl_lapack95_lp64’, ‘mkl_blas95_lp64’, ‘mkl_rt’]

Without the fast MKL BLAS, I found that the method 2 speedup was only about 5 times that of the Meng function. I would recommend using the MKL BLAS version if possible.
How to interpret non-fluorescent “reflectances” greater than 1.
Expand
Method 3, presented below, generates reflectance curves that are always strictly within the range 0 to 1, representing true (non-fluorescent) object reflectances. Method 2 can sometimes yield reflectances $>1$ . So it might seem natural to ask why method 2 would ever be preferred over method 3.

One reason is that the domain of applicability of method 2 is larger than that of method 3. True object colors reside in the object color solid, which is a subset of the domain of all “real” colors, which is bounded the spectral envelope. The application above on finding an illuminant SPD from a specified white point takes advantage of the expanded domain of method 2. There is no reason to limit an illuminant SPD to be $<1$ , which method 3 would require.

The question remains, “how should I interpret a reflectance curve that has components greater than 1?” The answer is to treat the “reflectance” as a means of specifying an emissive source, rather than as an object reflectance.

The stimulus that is emitted from a colored object under the action of an illuminant is the product of the illuminant, $W$ , and the object reflectance, $\rho$ . (In matrix algebra terms, it would be the element-by-element product of the two.) When we use a reflectance reconstruction algorithm to generate $\rho$ , we also specify a reference illuminant, let’s call it $W_{ref}$ . Thus the stimulus produced by the colored object characterized by $\rho$ , under the action of the reference illuminant, would be $\rho \times W_{ref}$ . When $\rho$ has components $>1$ , we should ignore the notion that the stimulus arises from an illuminant acting on an object reflectance, and instead, simply consider the product of the two to be a spectral power distribution produced by an emissive source.

Reflectance Reconstruction Method 3 (or LHTSS Least Hyperbolic Tangent Slope Squared)