function rho=LLSS(T,sRGB) % This is the Least Log Slope Squared (LLSS) algorithm for generating % a "reasonable" reflectance curve from a given sRGB color triplet. % The reflectance spans the wavelength range 380-730 nm in 10 nm increments. % Solves min sum(z_i+1 - z_i)^2 s.t. T exp(z) = rgb, where % z=log(reflectance), using Lagrangian formulation and Newton's method. % Allows reflectance values >1 to be in solution. % T is 3x36 matrix converting reflectance to D65-weighted linear rgb, % sRGB is a 3 element vector of target D65 referenced sRGB values in 0-255 range, % rho is a 36x1 vector of reconstructed reflectance values, all > 0, % For more information, see http://www.scottburns.us/subtractive-color-mixture/ % Written by Scott Allen Burns, March 2015. % Licensed under a Creative Commons Attribution-ShareAlike 4.0 International % License (http://creativecommons.org/licenses/by-sa/4.0/). % initialize outputs to zeros rho=zeros(36,1); % handle special case of (0,0,0) if all(sRGB==0) rho=0.0001*ones(36,1); return end % 36x36 difference matrix for Jacobian % having 4 on main diagonal and -2 on off diagonals, % except first and last main diagonal are 2. D=full(gallery('tridiag',36,-2,4,-2)); D(1,1)=2; D(36,36)=2; % compute target linear rgb values sRGB=sRGB(:)/255; % convert to 0-1 rgb=zeros(3,1); % remove gamma correction to get linear rgb for i=1:3 if sRGB(i)<0.04045 rgb(i)=sRGB(i)/12.92; else rgb(i)=((sRGB(i)+0.055)/1.055)^2.4; end end % 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+rgb]; % 39x1 function vector J=[D+diag(v),m2';m2,zeros(3)]; % 39x39 Jacobian matrix 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)