MATLAB实现CIE色彩计算：CCT、CRI与相关应用

%Appendix B: MATLAB script for calculating measures of light source color: CCT, CRI, GA, and FSI
% Second, calculate Correlated Color Temperature (CCT), Tc.
load ('isoTempLinesNewestFine.mat', 'T', 'ut', 'vt', 'tt'); 
% Find adjacent lines to (us, vs) 
n = length (T); 
index = 0; 
d1 = ((v-vt(1)) - tt(1)*(u-ut(1)))/sqrt(1+tt(1)*tt(1)); 
for i=2:n 
d2 = ((v-vt(i)) - tt(i)*(u-ut(i)))/sqrt(1+tt(i)*tt(i)); 
if (d1/d2 < 0) 
index = i; 
break; 
else 
d1 = d2; 
end 
end 
if index == 0 
Tc = -1; % Not able to calculate CCT, u, v coordinates outside range. 
return 
end 


% Calculate CCT by interpolation between isotemperature lines
Tc = 1/(1/T(index-1)+d1/(d1-d2)*(1/T(index)-1/T(index-1)));

% Third, calculate the Color Rendering Indices (CRI and its 14 indices)
% Calculate Reference Source Spectrum, spdref. 
if (Tc < 5000) 
c1 = 3.7418e-16; 
c2 = 1.4388e-2; 

spdref = c1 * (1e-9*wavelength_spd).^-5 ./ (exp(c2./(Tc.* 1e-9*wavelength_spd)) - 1); 
else 
if (Tc <= 25000) 
load('CIEDaySn','wavelength','S0','S1','S2'); 
if (Tc <= 7000) 
xd = -4.6070e9 / Tc.^3 + 2.9678e6 / Tc.^2 + 0.09911e3 / Tc + 0.244063; 
else 
xd = -2.0064e9 / Tc.^3 + 1.9018e6 / Tc.^2 + 0.24748e3 / Tc + 0.237040; 
end 
yd = -3.000*xd*xd + 2.870*xd - 0.275; 
M1 = (-1.3515 - 1.7703*xd + 5.9114*yd) / (0.0241 + 0.2562*xd - 0.7341*yd); 
M2 = (0.0300 - 31.4424*xd + 30.0717*yd) / (0.0241 + 0.2562*xd - 0.7341*yd); 
spdref = S0 + M1*S1 + M2*S2; 
spdref = interp1(wavelength,spdref,wavelength_spd); 
spdref(isnan(spdref)) = 0.0; 
else 
R = -1; 
return 
else 
end

 
% Load data for the spectral reflectance data of 14 color samples
TCS = load ('Tcs.txt');
TCS = TCS/1000;

 
% Interpolate TCS values from 5 nm to spd nm increments
TCS_1 = zeros (14,length(wavelength_spd));
wavelength_5 = 380:5:750;