function rho=ILLSS(T,sRGB) % This is the Iterative Least Log Slope Squared (ILLSS) 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. % It solves min sum(z_i+1 - z_i)^2 s.t. T exp(z) = rgb, K z = 0, where % z=log(reflectance), using Lagrangian approach and Newton's method. % Clips values >1 and repeats optimization until all reflectance <=1. % T is 3x36 matrix converting reflectance to linear rgb over the range 380-730 nm, % sRGB is a 3 element vector of target D65 referenced sRGB values in 0-255 range, % rho is a 36x1 vector of reflectance values (0->1] over wavelengths 380-730 nm, % Written by Scott Allen Burns, 4/11/15. % Licensed under a Creative Commons Attribution-ShareAlike 4.0 International % License (http://creativecommons.org/licenses/by-sa/4.0/). % For more information, see http://www.scottburns.us/subtractive-color-mixture/ % initialize output 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 % handle special case of (255,255,255) if all(sRGB==255) rho=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. 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 column vector 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 % outer iteration to get all refl <=1 maxouter=10; outer_count=0; % counter for outer iteration while (any(rho>1) && outer_count<=maxouter) || all(rho==0) % create K matrix for fixed refl constraints fixed_refl=find(rho>=1)'; numfixed=length(fixed_refl); K=zeros(numfixed,36); for i=1:numfixed K(i,fixed_refl(i))=1; end % initialize z=zeros(36,1); % starting point all zeros lambda=zeros(3,1); % starting point for lambda mu=zeros(numfixed,1); % starting point for mu maxit=50; % 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+K'*mu;m1+rgb;K*z]; % function vector J=[D+diag(v),[m2',K'];[m2;K],zeros(numfixed+3)]; % 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 mu=mu+delta(40:end); if all(abs(F)=maxit disp(['No inner loop solution found after ',num2str(maxit),' iterations.']) end outer_count=outer_count+1; end if outer_count