function [rho_R, rho_G, rho_B]=RGB_components_optim(T) % Solves min sum(c={r,g,b})(sum(i=1:35)(zc_i+1 - zc_i)^2) % s.t. T exp(zc) = s, c = {r,g,b}, s = {[1;0;0], [0;1;0], [0;0;1]} % rho_R + rho_G + rho_B < 1, % where zc=ln(rho_c), % using Matlab's constrained optimization function fmincon. % T is 3x36 matrix converting reflectance to D65-referenced linear rgb: T*rho=rgb. % Each rho is a 36x1 vector of reconstructed reflectance values, in the range (0-1]. % For more information, see: % http://scottburns.us/fast-rgb-to-spectrum-conversion-for-reflectances/ % Written by Scott Allen Burns, September 2018. % Licensed under a Creative Commons Attribution-ShareAlike 4.0 International % License (http://creativecommons.org/licenses/by-sa/4.0/). options=optimoptions(@fmincon,'Algorithm','sqp', 'MaxFunEvals',500000, 'MaxIter',3000,... 'ConstraintTolerance',1.0e-8, 'OptimalityTolerance',1.0e-8); z_soln=fmincon(@my_obj, zeros(108,1), [],[],[],[],[],[], @my_con, options); rho_R=exp(z_soln(1:36)); rho_G=exp(z_soln(37:72)); rho_B=exp(z_soln(73:108)); % nested functions for objective function and constraints function f=my_obj(z) % minimize norm of differences of successive log reflectances diff=z(1:35)-z(2:36); f=norm(diff); diff=z(37:71)-z(38:72); f=f+norm(diff); diff=z(73:107)-z(74:108); f=f+norm(diff); end function [c,ceq]=my_con(z) % require rgb of components to match ceq=zeros(9,1); current_rhoR=exp(z(1:36)); current_rgbR=T*current_rhoR; ceq(1:3)=current_rgbR-[1;0;0]; current_rhoG=exp(z(37:72)); current_rgbG=T*current_rhoG; ceq(4:6)=current_rgbG-[0;1;0]; current_rhoB=exp(z(73:108)); current_rgbB=T*current_rhoB; ceq(7:9)=current_rgbB-[0;0;1]; % upper limit on sum of components c = current_rhoR + current_rhoG + current_rhoB - 1; end end