function in = in_object_color_solid(A, W, XYZ)
% 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_18 Scott Allen Burns

% A is an nx3 array of color matching functions, by columns.
% W is an nx1 illuminant vector, any scale.
% XYZ is a three-element vector of tristimulus values, 0<Y<1.
% "in" is a logical indicating if XYZ is inside the object color solid.

% Notes
% -----
% 1.  If 'XYZ' is directly on the boundary of the object color solid,
%     this function will return false (limited to floating point precision).
% 2.  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.
% 
% 
% 3.  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)).
% 
% 4.  The resolution of cmfs (number of rows) affects the shape of the
%     computed spectral locus. See Section 6 of Reference [1] for
%     further discussion.
% 
% 5.  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

n = size(A,1); % resolution of CMFs
Aw = diag(W)*A/(A(:,2)'*W); % illuminant-referenced CMFs
xyz = XYZ/sum(XYZ(:)); % projective transform to x-y space
Y = XYZ(2); % target Y value
% Build a Y slice of object color solid
count = 0;
for step_width = 1:n
    Q = zeros(3,n);
    rho = zeros(n,1);
    % create step function square wave rho
    rho(1:step_width,1) = ones(step_width,1);
    for i = 1:n
        % compute TSVs of rho
        Q(:,i) = Aw'*rho;
        % shift rho with wraparound
        temp = rho(n);
        rho(2:n) = rho(1:n-1);
        rho(1) = temp;
    end
    % look for Q(2) values that straddle Y
    for i = 1:n-1
        % check if Y differences have a change in sign
        if (Q(2,i) - Y) * (Q(2,i+1) - Y) <= 0
            % interpolate to find where line segment crosses Y plane
            Xval=(Q(1,i+1)-Q(1,i))*(Y-Q(2,i))/(Q(2,i+1)-Q(2,i))+Q(1,i);
            Zval=(Q(3,i+1)-Q(3,i))*(Y-Q(2,i))/(Q(2,i+1)-Q(2,i))+Q(3,i);
            count = count + 1;
            x(count) = Xval/(Xval+Y+Zval);
            y(count) = Y/(Xval+Y+Zval);
        end
    end
end
% Q not in order, so find convex hull to sort xy in CCW order
K = convhull(x, y);
x = x(K);
y = y(K);
% check if xyz is within object color slice
[in_or_on, on] = inpolygon(xyz(1), xyz(2), x, y);
in = in_or_on && ~on;
