function rho=sRGB2refl(T,sRGB)
% Solves min sum(z_i+1 - z_i)^2 s.t. T exp(z) = rgb, z <=0, where
% z=log(reflectance), using Matlab's constrained optimization function fmincon.

% T is 3x36 matrix converting reflectance to D65-referenced 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, in the range (0-1],

% 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/).

% handle special case of (0,0,0)
if all(sRGB==0)
    rho=0.0001*ones(36,1);
    return
end

% compute target linear rgb values
sRGB=sRGB(:)/255; % convert to 0-1 column vector

rgb=zeros(size(sRGB));
% 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

% do optimization
options=optimset('Algorithm','SQP','MaxFunEvals',100000,'MaxIter',1000);
z_soln=fmincon(@sRGB2refl_obj,zeros(36,1),[],[],[],[],[],zeros(36,1),...
    @sRGB2refl_con,options);
rho=exp(z_soln);

% nested functions for objective function and constraints

    function f=sRGB2refl_obj(z)
        % minimize norm of differences of successive log reflectances
        diff=z(1:end-1)-z(2:end);
        f=norm(diff);
    end

    function [c,ceq]=sRGB2refl_con(z)
        c=[];
        % require rgb to match
        current_rho=exp(z);
        current_rgb=T*current_rho;
        ceq=current_rgb-rgb;
    end

end