function rho=LLSS_XYZ(A, W, XYZ) % This is the Least Log Slope Squared (LLSS) algorithm for generating % a "reasonable" reflectance curve, rho, from a given XYZ color triplet. % The reflectance spans the wavelength range 380-730 nm in 10 nm increments. % It solves min sum(z_i+1 - z_i)^2 s.t. T exp(z) = XYZ, where % z=log(rho), using a Lagrangian formulation and Newton's method. % A is the 3x36 matrix of color matching functions (380-730 nm in 10nm intervals). % W is a 36x1 illuminant vector spanning 380nm to 730 nm in 10 nm intervals. % XYZ is a 3 element vector of target tristimulus values in 0-1 range. % rho is a 36x1 vector of reflectance values over wavelengths 380-730 nm, % Written by Scott Allen Burns, 4/11/2015. % Licensed under a Creative Commons Attribution-ShareAlike 4.0 International % License (http://creativecommons.org/licenses/by-sa/4.0/). % For more information, see http://scottburns.us/reflectance-curves-from-srgb/ % NOTE: Not all input XYZ values will work. Some XYZ cannot be represented by % real-world reflectance curves, and the method will fail. % initialize output to zeros rho=zeros(36,1); % convert XYZ to column vector XYZ=XYZ(:); % Create CMFs weighted and normalized for illuminant W. % NOTE: to speed up this routine, compute T ahead of time and pass it in as % an input argument. T=A*diag(W)/(A(2,:)*W); % handle special case of (0,0,0) if all(XYZ<0.0001) rho=0.0001*ones(36,1); return end % 36x36 difference matrix having 4 on main diagonal and -2 on off diagonals, % except first and last main diagonal are 2. (Note: to speed up this routine, % pass D in as an input argument instead of forming it here.) D=full(gallery('tridiag',36,-2,4,-2)); D(1,1)=2; D(36,36)=2; % initialize z=zeros(36,1); % starting point all zeros lambda=zeros(3,1); % starting Lagrange mult maxit=100; % max number of iterations ftol=1.0e-8; % function solution tolerance deltatol=1.0e-8; % change in oper pt tolerance count=0; % iteration counter % Newton's method iteration while count <= maxit r=exp(z); v=-diag(r)*T'*lambda; % 36x1 m1=-T*r; % 3x1 m2=-T*diag(r); % 3x36 F=[D*z+v;m1+XYZ]; % 39x1 function vector J=[D+diag(v),m2';m2,zeros(3)]; % 39x39 Jacobian matrix if rcond(J) < 1.e-10 % The main source of failure of this method will be when an XYZ triplet % cannot be represented by a real-world object reflectance. disp('Singular matrix encountered.') disp('XYZ possibly out of "real-world" object color range.') return end delta=J\(-F); % solve Newton system of equations J*delta = -F z=z+delta(1:36); % update z lambda=lambda+delta(37:39); % update lambda if all(abs(F)