"""
Determine if XYZ is within the Object Color Solid
================================================================
Returns a logical true/false to indicate if a triplet of
tristimulus values falls strictly within the object color solid
for a given illuminant.
Version 2020_01_24 Scott Allen Burns
"""
import numpy as np
from scipy.spatial import ConvexHull
from shapely.geometry import Point
from shapely.geometry.polygon import Polygon
def in_object_color_solid(cmfs, w, XYZ):
"""
Checks to see if tristimulus values 'XYZ' are strictly
within the object color solid for a given illuminant, 'w',
returning true of false. If 'XYZ' is on the boundary of the
object color solid, this function will return false (limited
to floating point precision).
Parameters
----------
cmfs: nx3 array of color matching functions (not referenced to
an illuminant).
w: an n-element vector of illuminant relative spectral power
magnitudes, scaled arbitrarily.
XYZ: a three element vector of tristimulus values in the
0 <= Y <= 1 convention range.
Returns
-------
logical true or false.
Notes
-----
1. The object color solid boundary is constructed here by a brute-
force approach, considering all possible “optimal” colors. First,
a 0-1 step function is constructed of varying width (from one
wavelength bin to all n bins), and then each step function is
shifted along the wavelength axis (with wraparound) to produce
n shifted versions. This produces n^2 reflectance functions and
all possible optimal colors that can be represented with n
wavelength bins. As each step function is shifted, we check if the
Y tristimulus values corresponding to adjacent step functions
straddle the target Y value. If so, linear interpolation is used
to estimate the chromaticity coords of the intermediate color at
the intersection with the target Y plane. These chromaticity coords
are collected, and a convex hull function is used to sort the object
color solid slice in CCW order, producing a 2D polygon. Finally, the
chromaticity coords of 'XYZ' are tested to determine if they fall
within the Y-slice polygon.
2. Be sure to pass an XYZ triplet that sums to a number far enough
from zero to prevent overflow when computing the chromaticity
coordinates (xy = XYZ[0:2]/np.sum(XYZ)).
3. The resolution of cmfs (number of rows) affects the shape of the
computed spectral locus. See Section 6 of Reference [1] for
further discussion.
4. This code was translated from the MATLAB code presented in Section
4 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] Supplementary Documentation: Numerical Methods for Smoothest
Reflectance Reconstruction, 2019
http://scottburns.us/suppl-docs-num-meth-smoothest-refl-recon/#Sec3
Acknowledgements
----------------
Thanks to Adam J. Burns and Mark (kram1032) for assistance
with translating from MATLAB to Python.
"""
# number of wavelength bins (rows)
n = cmfs.shape[0]
# form illuminant-w-referenced cmfs
w1 = np.squeeze(w)
w_norm = w1/(w1.dot(cmfs[:, 1]))
cmfs_w = np.diag(w_norm).dot(cmfs)
# get chromaticity coords
xy = XYZ[0:2]/np.sum(XYZ)
# target Y slice
Y = XYZ[1]
x = np.array([])
y = np.array([])
# Build a Y slice of object color solid
for step_width in range(1, n+1): # 1 to n
q = np.zeros((n, 3))
rho = np.zeros(n)
# create step function square wave rho
rho[0:step_width] = np.ones(step_width)
for i in range(n):
# compute TSVs of rho
q[i, :] = cmfs_w.T.dot(rho)
# shift rho right with wraparound
temp = rho[n-1]
rho[1:n] = rho[0:n-1]
rho[0] = temp
# look for Q(1,:) values that straddle Y
for i in range(n-1): # 0 to n-2
# check if Y differences have a change in sign
if (q[i, 1] - Y) * (q[i+1, 1] - Y) <= 0:
# interpolate to find where line segment crosses Y plane
xval = (q[i+1,0]-q[i,0])*(Y-q[i,1])/(q[i+1,1]-q[i,1])+q[i,0]
zval = (q[i+1,2]-q[i,2])*(Y-q[i,1])/(q[i+1,1]-q[i,1])+q[i,2]
x = np.append(x, xval/(xval+Y+zval))
y = np.append(y, Y/(xval+Y+zval))
# Q not in order, so find convex hull to sort xy in CCW order
qxy = np.array([x, y]).T
hull = ConvexHull(qxy)
poly = Polygon(hull.points)
p = Point(xy)
# check if point is in polygon
return poly.contains(p)