File talk:CIE CRI TCS under FL4.svg

Source

function result=CRI(SPD,CCT,L)

[LAMBDA, XFCN, YFCN, ZFCN]=colorMatchFcn('1931_full');

TCS=textread('CRI_TCS.csv', '', 'delimiter', ',', 'emptyvalue', NaN);
TCS_spd=interp1(TCS(:,1),TCS(:,2:end),LAMBDA);

if CCT<5000             % use a black body
    c2=1.4387752;       % second radiation constant
    ref_spd=LAMBDA.^(-5)./(exp(c2*1e+7./(CCT*LAMBDA))-1);
else
    ref_spd=CIE_daylight(CCT,LAMBDA);
end

[ref_x,ref_y,ref_Y]=chromaticity(ref_spd,LAMBDA);
ref_spd=ref_spd*100/ref_Y; % normalize "ref_spd" to Y=100;
ref_XYZ=applycform([ref_x ref_y 100],makecform('xyl2xyz'));
ref_uvl=applycform(ref_XYZ,makecform('xyz2uvl'));
ref_TCS=repmat(ref_spd',1,15).*TCS_spd;
[ref_X,ref_Y,ref_Z]=tristimulus(ref_TCS,LAMBDA);
ref_tcs_uvl=applycform([ref_X' ref_Y' ref_Z'],makecform('xyz2uvl'));

test_spd=interp1(L,SPD,LAMBDA);
test_spd(isnan(test_spd))=0; % zero out the invalid ones
[test_x,test_y,test_Y]=chromaticity(SPD,L); % no point in using test_spd
test_XYZ=applycform([test_x test_y test_Y],makecform('xyl2xyz'));
test_uvl=applycform(test_XYZ,makecform('xyz2uvl'));
test_spd=test_spd*100/test_Y; % normalize "test_spd" to Y=100
test_TCS=repmat(test_spd',1,15).*TCS_spd;
[test_X,test_Y,test_Z]=tristimulus(test_TCS,LAMBDA);
test_tcs_uvl=applycform([test_X' test_Y' test_Z'],makecform('xyz2uvl'));

cr=(4-ref_uvl(1)-10*ref_uvl(2))./ref_uvl(2);
ct=(4-test_uvl(1)-10*test_uvl(2))./test_uvl(2);
dr=(1.708*ref_uvl(2)-1.481*ref_uvl(1)+0.404)./ref_uvl(2);
dt=(1.708*test_uvl(2)-1.481*test_uvl(1)+0.404)./test_uvl(2);

corr_uvl=test_tcs_uvl;
% this section shows the tristimulus values before chromatic adaptation
% corr_uvw(:,3)=25*corr_uvl(:,3).^(1/3)-17;
% corr_uvw(:,1)=13*corr_uvw(:,3).*(corr_uvl(:,1)-repmat(test_uvl(1),15,1));
% corr_uvw(:,2)=13*corr_uvw(:,3).*(corr_uvl(:,2)-repmat(test_uvl(2),15,1));
for i=1:15
    u=test_tcs_uvl(i,1);
    v=test_tcs_uvl(i,2);
    c=(4-u-10*v)/v;
    d=(1.708*v-1.481*u+0.404)/v;
    corr_uvl(i,1)=(10.872+0.404*c*cr/ct-4*d*dr/dt)/(16.518+1.481*c*cr/ct-d*dr/dt);
    corr_uvl(i,2)=5.52/(16.518+1.481*c*cr/ct-d*dr/dt);
end

corr_uvw=corr_uvl;
corr_uvw(:,3)=25*corr_uvl(:,3).^(1/3)-17;
corr_uvw(:,1)=13*corr_uvw(:,3).*(corr_uvw(:,1)-repmat(test_uvl(1),15,1));
corr_uvw(:,2)=13*corr_uvw(:,3).*(corr_uvw(:,2)-repmat(test_uvl(2),15,1));

ref_tcs_uvw=corr_uvw;
ref_tcs_uvw(:,3)=25*ref_tcs_uvl(:,3).^(1/3)-17;
ref_tcs_uvw(:,1)=13*ref_tcs_uvw(:,3).*(ref_tcs_uvl(:,1)-repmat(ref_uvl(1),15,1));
ref_tcs_uvw(:,2)=13*ref_tcs_uvw(:,3).*(ref_tcs_uvl(:,2)-repmat(ref_uvl(2),15,1));

result=100-4.6*sqrt(sum((corr_uvw-ref_tcs_uvw)'.^2));

% The next section simply draws the CAT diagram and is optional
N=1000; % size of the image
a=0.19; b=0.31; c=0.32; d=.37;          % CRI TCS
x=linspace(a,b,N*min(1,(b-a)/(d-c)));   % right number of pixels
y=linspace(c,d,N*min(1,(d-c)/(b-a)));   % right number of pixels
[u,v,Y]=meshgrid(x,y,50);
xyz=applycform(cat(3,u,v,Y),makecform('uvl2xyz'));
M=max(xyz,[],3);
norm_XYZ=xyz./cat(3,M,M,M);
mean_XYZ=norm_XYZ*0.5./repmat(mean(norm_XYZ,3), [1 1 3]);
srgb=applycform(mean_XYZ,makecform('xyz2srgb'));
image(x,y,srgb);
axis image;
set(gca,'YDir','normal');   % right way up
hold;

plot(ref_uvl(1),ref_uvl(2),'x',test_uvl(1),test_uvl(2),'x');
line([test_tcs_uvl(1:8,1)' ; corr_uvl(1:8,1)'],...
[ test_tcs_uvl(1:8,2)' ; corr_uvl(1:8,2)'],...
'Color','k','LineWidth',1);
text(ref_uvl(1),ref_uvl(2),'2940K')
text(test_uvl(1),test_uvl(2),'FL4')
plot(ref_tcs_uvl(1:8,1),ref_tcs_uvl(1:8,2),'o')
for i=1:8
    text(ref_tcs_uvl(i,1),ref_tcs_uvl(i,2),strcat('TCS',i+48))
end