function in = in_spectral_locus(A, XYZ)
% Returns a logical true/false to indicate if a triplet of
% tristimulus values falls strictly within the spectral locus.

% Version 2020_01_18 Scott Allen Burns

% A is an nx3 array of color matching functions, by columns.
% XYZ is a three-element vector of tristimulus values, any
%     scale; Y > 0.
% in is a logical indicating if XYZ is strictly inside the
%     spectral locus.

% 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 function computes the projective transformation from
%    XYZ space to the x-y plane, so that the test for being within
%    the locus is reduced a dimension to being within a 2-D polygon.
% 
% 3. While the shape of the object color solid depends on the
%    illuminant, the location of the spectral locus in the x-y plane
%    is independent of illuminant.
% 
% 4. 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)).
% 
% 5. The resolution of cmfs (number of rows) affects the shape of the
%    computed spectral locus. See Section 6 of Reference [1] for
%    further discussion.
% 
% 6. This code was translated from the MATLAB code presented in Section
%    3 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

locus = A(:,1:2)./sum(A,2); % get x,y of spectral locus
xyz = XYZ/sum(XYZ(:)); % convert XYZ to xyz
[in_or_on, on] = inpolygon(xyz(1), xyz(2), locus(:,1), locus(:,2));
in = in_or_on && ~on;
