"""
Reflectance Reconstruction from Tristimulus Values (Real Colors)
================================================================
Computes a reconstructed spectral reflectance curve from an
illuminant-referenced triplet of tristimulus values
representing a real color (within the spectral locus).
Version 2020_01_24 Scott Allen Burns
"""
import numpy as np
from scipy import linalg
import warnings
warnings.simplefilter("ignore")
def method_2(d, cmfs_w, XYZ_w, max_iterations=20, tolerance=1e-8):
"""
Computes a reconstructed spectral reflectance curve from
an illuminant-referenced triplet of tristimulus values
representing a real color (within the spectral locus).
Parameters
----------
d: nxn array of finite differencing constants, obtained
from function method_2_prep.py. (n is the number
of wavelength bins (rows) used in the CMFs.)
cmfs_w: nx3 array of illuminant-w-referenced CMFs,
obtained from function method_2_prep.py.
XYZ_w: a three element vector of illum-w-referenced
tristimulus values in the 0 <= Y <= 1 convention range.
max_iterations: upper limit on iterations (defaults to 20).
tolerance: determines when the linear equations for a
stationary point are satisfied (default is 1e-8).
Returns
-------
nx1 vector of spectral reflectance values (or zeros if failure).
Notes
-----
The input parameters d and cmfs_w can be generated by function
method_2_prep.py. For a given illuminant, method_2_prep
needs to be called only once. Multiple calls to method_2
may then follow for various tristimulus values. This is a
typical use case, and avoids repeated computation of d and
cmfs_w.
The XYZ_w supplied to method_2 must fall strictly within the
spectral locus (i.e., be a real color), or the method will fail
and return a vector of zeros.
The reflectance curve generated will always have positive
values, but may have some values >1 (which can be interpreted
as an emissive source divided by the illuminant). To obtain
reflectance curves strictly in the range 0 to 1, representing
valid object reflectances, use method_3.py instead.
method_2 is based on "method 2" of Reference [1]. It also
corresponds to method "LLSS" (Least Log Slope Squared) of
Reference [2].
References
----------
[1] Burns SA. Numerical methods for smoothest reflectance
reconstruction. Color Research & Application, Vol 45,
No 1, 2020, pp 8-21.
[2] Generating Reflectance Curves from sRGB Triplets, 2015
http://scottburns.us/reflectance-curves-from-srgb/#LLSS
Acknowledgements
----------------
Thanks to Adam J. Burns and Mark (kram1032) for assistance
with translating from MATLAB to Python.
"""
n = cmfs_w.shape[0] # number of rows
z = np.zeros(n)
lam = np.zeros(3)
for _ in range(max_iterations):
r = np.exp(z)
dr = np.diag(r)
dra = dr.dot(cmfs_w)
v = dra.dot(lam)
f = np.block([d.dot(z)+v, cmfs_w.T.dot(r)-XYZ_w])
j = np.block([[d+np.diag(v), dra], [dra.T, np.zeros((3, 3))]])
try:
delta = linalg.solve(j, -f, assume_a='sym')
except (linalg.misc.LinAlgError, ValueError):
print('Ill-conditioned or singular linear system detected.')
print('Check to make sure XYZ_w is within the spectral locus.')
return np.zeros(n)
z += delta[0:n]
lam += delta[n:n+3]
if np.all(np.abs(f) < tolerance):
return np.exp(z)
print('Maximum number of iterations reached.')
print('Check to make sure XYZ_w is within the spectral locus.')
return np.zeros(n)