function rho=ILSS(B11,B12,sRGB) % This is the Iterative Least Slope Squared (ILSS) 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(rho_i+1 - rho_i)^2 s.t. T rho = rgb, K1 rho = 1, K0 rho = 0, % using Lagrangian formulation and iteration to keep all rho (0-1]. % B11 is upper-left 36x36 part of inv([D,T';T,zeros(3)]) % B12 is upper-right 36x3 part of inv([D,T';T,zeros(3)]) % 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/26/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/ rho=ones(36,1)/2; % initialize output to 0.5 rhomin=0.00001; % smallest refl value % handle special case of (255,255,255) if all(sRGB==255) rho=ones(36,1); return end % handle special case of (0,0,0) if all(sRGB==0) rho=rhomin*ones(36,1); return end % 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 R=B12*rgb; % iteration to get all refl 0-1 maxit=10; % max iterations count=0; % counter for iteration while ( (any(rho>1) || any(rho=1; fixed_upper=find(fixed_upper_logical); num_upper=length(fixed_upper); K1=zeros(num_upper,36); for i=1:num_upper K1(i,fixed_upper(i))=1; end % create K0 matrix for fixed refl at rhomin fixed_lower_logical = rho<=rhomin; fixed_lower=find(fixed_lower_logical); num_lower=length(fixed_lower); K0=zeros(num_lower,36); for i=1:num_lower K0(i,fixed_lower(i))=1; end % set up linear system K=[K1;K0]; C=B11*K'/(K*B11*K'); % M*K'*inv(K*M*K') rho=R-C*(K*R-[ones(num_upper,1);rhomin*ones(num_lower,1)]); rho(fixed_upper_logical)=1; % eliminate FP noise rho(fixed_lower_logical)=rhomin; % eliminate FP noise count=count+1; end if count>=maxit disp(['No solution found after ',num2str(maxit),' iterations.']) end