function rho=LHTSS(T,sRGB) % This is the Least Hyperbolic Tangent Slope Squared (LHTSS) 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 ((tanh(z)+1)/2) = rgb, % using Lagrangian formulation and Newton's method. % Reflectance will always be in the open interval (0->1). % 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->1), % Written by Scott Allen Burns, May 2019. % 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/ % 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 % handle special case of (255,255,255) if all(sRGB==255) rho=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 rho=1/2 lambda=zeros(3,1); % starting Lagrange mult maxit=100; % max number of iterations ftol=1.0e-8; % function solution tolerance count=0; % iteration counter % Newton's method iteration while count <= maxit d0 = (tanh(z) + 1)/2; d1 = diag((sech(z).^2)/2); d2 = diag(-sech(z).^2.*tanh(z)); F = [D*z + d1*T'*lambda; T*d0 - rgb]; % 39x1 F vector J = [D + diag(d2*T'*lambda), d1*T'; T*d1, zeros(3)]; % 39x39 J 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)