Matlab Program for CIE

%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;


for i = 1:14
TCS_1(i,:) = interp1(wavelength_5,TCS(i,:),wavelength_spd');
TCS_1(i,isnan(TCS_1(i,:))) = 0.0; % remove NaN from vector.
end


% Calculate u, v chromaticity coordinates of samples under test illuminant, uk, vk and
% reference illuminant, ur, vr.
uki = zeros(1,14);
vki = zeros(1,14);
uri = zeros(1,14);
vri = zeros(1,14);
Yknormal = 100 / Y;
Yk = Y*Yknormal;
uk = 4*X/(X+15*Y+3*Z);
vk = 6*Y/(X+15*Y+3*Z);
X = sum(spdref .* xbar);
Y = sum(spdref .* ybar);
Z = sum(spdref .* zbar);
Yrnormal = 100 / Y;
Yr = Y*Yrnormal;
ur = 4*X/(X+15*Y+3*Z);
vr = 6*Y/(X+15*Y+3*Z);

for i = 1:14
X = sum(spd .* TCS_1(i,:)' .* xbar);
Y = sum(spd .* TCS_1(i,:)' .* ybar);
Z = sum(spd .* TCS_1(i,:)' .* zbar);
Yki(i) = Y*Yknormal;
uki(i) = 4*X/(X+15*Y+3*Z);
vki(i) = 6*Y/(X+15*Y+3*Z);
X = sum(spdref .* TCS_1(i,:)' .* xbar);
Y = sum(spdref .* TCS_1(i,:)' .* ybar);
Z = sum(spdref .* TCS_1(i,:)' .* zbar);
Yri(i) = Y*Yrnormal;
uri(i) = 4*X/(X+15*Y+3*Z);
vri(i) = 6*Y/(X+15*Y+3*Z);
end


% Check tolerance for reference illuminant
DC = sqrt((uk-ur).^2 + (vk-vr).^2);
if DC>0.0054
return
end


% Apply adaptive (p

erceived) color shift.
ck = (4 - uk - 10*vk) / vk;
dk = (1.708*vk + 0.404 - 1.481*uk) / vk;
cr = (4 - ur - 10*vr) / vr;
dr = (1.708*vr + 0.404 - 1.481*ur) / vr;

for i = 1:14
cki = (4 - uki(i) - 10*vki(i)) / vki(i);
dki = (1.708*vki(i) + 0.404 - 1.481*uki(i)) / vki(i);
ukip(i) = (10.872 + 0.404*cr/ck*cki - 4*dr/dk*dki) / (16.518 + 1.481*cr/ck*cki - dr/dk*dki);
vkip(i) = 5.520 / (16.518 + 1.481*cr/ck*cki - dr/dk*dki);
end


% Transformation into 1964 Uniform space coordinates.

for i = 1:14
Wstarr(i) = 25*Yri(i).^.333333 - 17;
Ustarr(i) = 13*Wstarr(i)*(uri(i) - ur);
Vstarr(i) = 13*Wstarr(i)*(vri(i) - vr);
Wstark(i) = 25*Yki(i).^.333333 - 17;
Ustark(i) = 13*Wstark(i)*(ukip(i) - ur); % after applying the adaptive color shift, u'k = ur
Vstark(i) = 13*Wstark(i)*(vkip(i) - vr); % after applying the adaptive color shift, v'k = vr
end


% Determination of resultant color shift, delta E.
deltaE = zeros(1,14);

for i = 1:14
deltaE(i) = sqrt((Ustarr(i) - Ustark(i)).^2 + (Vstarr(i) - Vstark(i)).^2 + (Wstarr(i) - Wstark(i)).^2);
Ra14(i) = 100 - 4.6*deltaE(i);
end
Ra = sum(Ra14(1:8))/8;


% fourth, calculate the gamut area formed by the 8 CIE standard color samples
ukii=[uki(:,1:8),uki(1)];
vkii=1.5*[vki(:,1:8),vki(1)];
Ga=polyarea(ukii,vkii);
% Normalize gamut area to equal energy source
Ga=Ga/0.00728468*100;


% fifth, calculate the FSI (full spectrum index)
% Calculates the Full-spectrum Index

% FSI = fsi(spd,startw,endw,incrementw)
% spd is a single column or row vector and startw, endw and incrementw specify the
% starting, ending and increment of wavelength in nm.

% FSI = fsi(spd)
% spd is a two column matrix with wavelength values in column 1 and spd values in
% column 2. Wavelength values are assumed to be in units of nm.
if length(varargin)==0
[rows columns] = size(spd);
if columns > 2
error('Not column oriented data. Try transposing spd');
end
wavelength_spd = spd(:,1);
spd = spd(:,2);
else
startw = varargin{1}
endw = varargin{2}
incrementw = varargin{3}
wavelength_spd = (startw:incrementw:endw)';
[rows columns] = size(spd);
if columns > 1
error('Detected multiple columns of data. Try transposing spd');
end
end


% Interpolate to wavelength interval of 1nm from 380nm to 730nm
numWave = 351;
t=(380:1:730)';
spd=interp1(wavelength_spd,spd,t,'spline');
spd(isnan(spd)) = 0.0;
spd = spd/sum(spd); % Normalize the relative spd so that the total power equals 1
%Equal energy cumulative spd
EEcum=(1/numWave:1/numWave:1)';
%Calculate FSI

for j=1:numWave
cum = cumsum(spd); % A MatLab function for cumulative sums
sqrDiff = (cum-EEcum).^2;
sumSqrDiff(j)=sum(sqrDiff);
spd=circshift(spd,1);
end
FSI=mean(sumSqrDiff);


% Note: To make FSI more directly comparable to CRI, FSI values in Figure 14 have been converted to a 0-100 scale, with an equal energy spectrum defined as having an FSI value of 100, and all practical light sources

having FSI values lower than 100; a monochromatic light source (e.g., low pressure sodium) has a value of 0.