# 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).
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.

• 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

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)

##### 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 () 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 matrix is replaced with a matrix relating reflectance to RGB. Reference [2] demonstrates this for the sRGB color space, where is replaced with a matrix called .

Another possibile color space is that of cone sensitivity functions (cone fundamentals). The cone fundamentals would be used in place of the 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, , instead of , and pass the illuminant white point tristimlulus values, , instead of , and method 2 will generate a suitable illuminant curve, . If the white point has =1, then the illuminant will also be appropriately scaled so that an illuminant-referenced 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 .

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 . 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 , 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, , and the object reflectance, . (In matrix algebra terms, it would be the element-by-element product of the two.) When we use a reflectance reconstruction algorithm to generate , we also specify a reference illuminant, let’s call it . Thus the stimulus produced by the colored object characterized by , under the action of the reference illuminant, would be . When has components , 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)

##### Examples
• MATLAB/Octave example.
Show example
The following code demonstrates the method 3 code for the case of XYZ both inside and outside the object color solid, using an equal energy reference illuminant. The CMFs (A) are hard coded.

% Test the function method_3.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_3_prep(A,W);

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

disp('Test outside object color solid:');
XYZw = [0.2, 0.7, 0.5];
% reflectance reconstruction
rho = method_3(D,Aw,XYZw);
disp(rho');

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

Test inside object color solid:
Columns 1 through 10
0.0582 0.0582 0.0583 0.0584 0.0589 0.0606 0.0648 0.0732 0.0874 0.1106
Columns 11 through 20
0.1461 0.1964 0.2602 0.3305 0.3927 0.4317 0.4406 0.4200 0.3758 0.3186
Columns 21 through 30
0.2610 0.2123 0.1760 0.1512 0.1351 0.1253 0.1194 0.1162 0.1144 0.1135
Columns 31 through 36
0.1130 0.1128 0.1127 0.1127 0.1126 0.1126

Test outside object color solid:
Maximum number of iterations reached.
Check to make sure XYZw is within the object color solid.
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 3 code for the case of XYZ both inside and outside the object color solid, using an equal energy reference illuminant. The CMFs (cmfs) are hard coded.

# Test the function method_3.py

# Version 2020_01_24 Scott Allen Burns

import numpy as np
from method_3_prep import method_3_prep
from method_3 import method_3

# 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_3_prep(cmfs, w)

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

print('Test outside object color solid')
XYZ_w = np.array([0.2, 0.7, 0.5])
# reflectance reconstruction
rho = method_3(d, cmfs_w, XYZ_w)
print(rho)

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

Test inside object color solid:
[0.05819499 0.05820581 0.05825025 0.05840874 0.05891865 0.06055314
0.06481247 0.0731528 0.08741221 0.11059585 0.14614887 0.19637775
0.26022948 0.33048433 0.3926965 0.43165695 0.44060517 0.41999403
0.37576787 0.3185753 0.26098453 0.21232254 0.17602657 0.15115906
0.13511843 0.12527071 0.11944401 0.11616137 0.11439869 0.11349234
0.11303233 0.11280935 0.11270106 0.11265007 0.11262823 0.11262098]

Test outside object color solid
Ill-conditioned or singular linear system detected.
Check to make sure XYZ_w is within the object color solid.
[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

Three functions are provided. The first two functions compute some preparatory arrays that are passed to the third function. The first two functions differ in how the illuminant white points are specified, and either one can be used. The “XYZ” version uses white point tristimulus values; the “xy” version uses white point chromaticity coordinates. The preparatory arrays they create can be reused for a series of source XYZ values, which is a common use case. Please note that the Chromatic Adaptation functions make use of the “method_2” functions presented above.

##### MATLAB/Octave
• Preparatory function given illuminant white point XYZ source code.
• Preparatory function given illuminant white point xy source code.
##### Python
• Preparatory function given illuminant white point XYZ source code.
• Preparatory function given illuminant white point xy source code.
##### Examples
• MATLAB/Octave example.
Show example
The following code demonstrates the adaptation code for the case of XYZ = (0.2, 0.3, 0.1) going from illuminant A to illuminant D65, for both the XYZ and xy versions of the prep code. The CMFs (A) are hard coded.

% 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];
% source tristimulus values
XYZs = [0.2, 0.3, 0.1];
degree = 1;

% test the XYZ version
% source illuminant white point (illum A)
XYZwps = [1.09847, 1, 0.35582];
% destination illuminant white point (D65)
XYZwpd = [0.95047, 1, 1.08883];
disp(['Source XYZs = (',num2str(XYZs),') adapted from']);
disp(['white point XYZs = (',num2str(XYZwps),') to']);
disp(['white point XYZd = (',num2str(XYZwpd),') is']);
disp(['XYZd = (',num2str(XYZd'),')']);
disp(' ');

% test the xy version
% source illuminant white point (illum A)
xy_wps = [0.44757, 0.40745];
% destination illuminant white point (D65)
xy_wpd = [0.31271, 0.32902];
disp(['Source XYZs = (',num2str(XYZs),') adapted from']);
disp(['white point xy = (',num2str(xy_wps),') to']);
disp(['white point xy = (',num2str(xy_wpd),') is']);
disp(['XYZd = (',num2str(XYZd'),')']);
disp(' ');

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

Source XYZs = (0.2 0.3 0.1) adapted from
white point XYZs = (1.0985 1 0.35582) to
white point XYZd = (0.95047 1 1.0888) is
XYZd = (0.17074 0.3 0.24253)

Source XYZs = (0.2 0.3 0.1) adapted from
white point xy = (0.44757 0.40745) to
white point xy = (0.31271 0.32902) is
XYZd = (0.17073 0.3 0.24253)

• Python example.
Show example
The following code demonstrates the adaptation code for the case of XYZ = (0.2, 0.3, 0.1) going from illuminant A to illuminant D65, for both the XYZ and xy versions of the prep code. The CMFs (cmfs) are hard coded.

# Version 2020_01_24 Scott Allen Burns

import numpy as np

# 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]])

# source tristimulus values
XYZ_s = np.array([0.2, 0.3, 0.1])
degree = 1

# test the XYZ version
# source illuminant white point (illum A)
XYZ_wps = np.array([1.09847, 1, 0.35582])
# destination illuminant white point (D65)
XYZ_wpd = np.array([0.95047, 1, 1.08883])
d, cmfs_s, cmfs_d = adaptation_prep_XYZ(cmfs, XYZ_wps, XYZ_wpd, degree)
XYZ_d = adaptation(d, cmfs_s, cmfs_d, XYZ_s)
print("Source XYZ_s = " + str(XYZ_s) + " adapted from")
print("white point " + str(XYZ_wps) + " to")
print("white point " + str(XYZ_wpd) + " is")
print("XYZ_d = " + str(XYZ_d))
print("(degree of adaptation = " + str(degree) + ")")
print(" ")

# test the xy version
# source illuminant white point (illum A)
xy_wps = np.array([0.44757, 0.40745])
# destination illuminant white point (D65)
xy_wpd = np.array([0.31271, 0.32902])
d, cmfs_s, cmfs_d = adaptation_prep_xy(cmfs, xy_wps, xy_wpd, degree)
XYZ_d = adaptation(d, cmfs_s, cmfs_d, XYZ_s)
print("Source XYZ_s = " + str(XYZ_s) + " adapted from")
print("white point " + str(xy_wps) + " to")
print("white point " + str(xy_wpd) + " is")
print("XYZ_d = " + str(XYZ_d))
print("(degree of adaptation = " + str(degree) + ")")
print(" ")

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

Source XYZ_s = [0.2 0.3 0.1] adapted from
white point [1.09847 1. 0.35582] to
white point [0.95047 1. 1.08883] is
XYZ_d = [0.17074153 0.3 0.24252612] (degree of adaptation = 1)

Source XYZ_s = [0.2 0.3 0.1] adapted from
white point [0.44757 0.40745] to
white point [0.31271 0.32902] is
XYZ_d = [0.17073486 0.3 0.24253345] (degree of adaptation = 1)

Process finished with exit code 0

### Utilities

These functions determine whether a target XYZ triplet is suitable for one of the reflectance reconstruction methods above.

#### Test if XYZ is within Spectral Locus

This function returns true if the XYZ falls strictly within the spectral locus, confirming that it is suitable for method 2.

##### Examples
• MATLAB/Octave example.
Show example
The following code demonstrates the in_spectral_locus code for three cases: an XYZ within the object color solid, an XYZ within the spectral locus but outside the object color solid, and an XYZ outside the spectral locus. The CMFs (A) are hard coded.

% Test the function in_spectral_locus.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];

disp('Test XYZ inside object color solid:');
XYZ = [0.2; 0.3; 0.1];
result = in_spectral_locus(A, XYZ);
disp(result);
disp(' ');

disp('Test XYZ inside spectral locus but outside object color solid:');
XYZ = [0.2; 0.7; 0.5];
result = in_spectral_locus(A, XYZ);
disp(result);
disp(' ');

disp('Test XYZ outside spectral locus:');
XYZ = [0.8; 0.2; 0.5];
result = in_spectral_locus(A, XYZ);
disp(result);

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

Test XYZ inside object color solid:
1

Test XYZ inside spectral locus but outside object color solid:
1

Test XYZ outside spectral locus:
0

• Python example.
Show example
The following code demonstrates the in_spectral_locus code for three cases: an XYZ within the object color solid, an XYZ within the spectral locus but outside the object color solid, and an XYZ outside the spectral locus. The CMFs (cmfs) are hard coded.

# Test the function in_spectral_locus.py

# Version 2020_01_24 Scott Allen Burns

import numpy as np
from in_spectral_locus import in_spectral_locus

# 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]])

print('Test XYZ inside object color solid:')
XYZ = np.array([0.2, 0.3, 0.1])
result = in_spectral_locus(cmfs, XYZ)
print(result)

print('Test XYZ inside spectral locus but outside object color solid:')
XYZ = np.array([0.2, 0.7, 0.5])
result = in_spectral_locus(cmfs, XYZ)
print(result)

print('Test XYZ outside spectral locus:')
XYZ = np.array([0.8, 0.2, 0.5])
result = in_spectral_locus(cmfs, XYZ)
print(result)

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

Test XYZ inside object color solid:
True
Test XYZ inside spectral locus but outside object color solid:
True
Test XYZ outside spectral locus:
False

Process finished with exit code 0

#### Test if XYZ is within Object Color Solid

This function returns true if the XYZ falls strictly within the object color solid, confirming that it is suitable for method 3.

##### Examples
• MATLAB/Octave example.
Show example
The following code demonstrates the in_object_color_solid code for three cases: an XYZ within the object color solid, an XYZ within the spectral locus but outside the object color solid, and an XYZ outside the spectral locus. The CMFs (A) are hard coded.

% Test the function in_object_color_solid.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);

disp('Test XYZ inside object color solid:');
XYZ = [0.2; 0.3; 0.1];
result = in_object_color_solid(A, W, XYZ);
disp(result);
disp(' ');

disp('Test XYZ inside spectral locus but outside object color solid:');
XYZ = [0.2; 0.7; 0.5];
result = in_object_color_solid(A, W, XYZ);
disp(result);
disp(' ');

disp('Test XYZ outside spectral locus:');
XYZ = [0.8; 0.2; 0.5];
result = in_object_color_solid(A, W, XYZ);
disp(result);

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

Test XYZ inside object color solid:
1

Test XYZ inside spectral locus but outside object color solid:
0

Test XYZ outside spectral locus:
0

• Python example.
Show example
The following code demonstrates the in_object_color_solid code for three cases: an XYZ within the object color solid, an XYZ within the spectral locus but outside the object color solid, and an XYZ outside the spectral locus. The CMFs (cmfs) are hard coded.

# Test the function in_object_color_solid.py

# Version 2020_01_24 Scott Allen Burns

import numpy as np
from in_object_color_solid import in_object_color_solid

# 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)

print('Test XYZ inside object color solid:')
XYZ = np.array([0.2, 0.3, 0.1])
result = in_object_color_solid(cmfs, w, XYZ)
print(result)

print('Test XYZ inside spectral locus but outside object color solid:')
XYZ = np.array([0.2, 0.7, 0.5])
result = in_object_color_solid(cmfs, w, XYZ)
print(result)

print('Test XYZ outside spectral locus:')
XYZ = np.array([0.8, 0.2, 0.5])
result = in_object_color_solid(cmfs, w, XYZ)
print(result)

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

Test XYZ inside object color solid:
True
Test XYZ inside spectral locus but outside object color solid:
False
Test XYZ outside spectral locus:
False

Process finished with exit code 0

### A note on Python packages

The Python SciPy package used in the code above is dependent on the NumPy package. It is important to install the “Numpy+MKL” version of NumPy, which is linked to the Intel Math Kernel Library and includes required DLLs in the numpy.DLLs directory, and uses an optimized BLAS (Basic Linear Algebra Subprograms) library. The SciPy functions for linear algebra will solve MUCH faster this way.

### References

[1] Burns, Scott A., Numerical Methods for Smoothest Reflectance Reconstruction, Color Research & Application, Vol 45, No 1, 2020, pp 8-21.

[2] Burns, Scott A., Generating Reflectance Curves from sRGB Triplets, 2015 http://scottburns.us/reflectance-curves-from-srgb

[3] Burns, Scott A., Chromatic Adaptation Transform by Spectral Reconstruction, Color Research & Application, Vol 44, No 5, 2019, pp 682-693.

[4] Burns, Scott A., Supplementary Documentation: Numerical Methods for Smoothest Reflectance Reconstruction, http://scottburns.us/suppl-docs-num-meth-smoothest-refl-recon/

[5] Meng, J., Simon, F., Hanika, J., and Dachsbacher, C., Physically Meaningful Rendering using Tristimulus Colours, Computer Graphics Forum, Vol 34, No 4, 2015, pp 31-40.

### Acknowledgements

Special thanks to Adam J Burns for orientation to the Python development environment and python programming. Also thanks to Mark (kram1032) for assistance with MATLAB->Python translation.