%Chapter 7,  Fig 7.8
%Eric Dubois, updated 2018-12-13
%This program generates all quantities related to the CIE 1931 RGB and XYZ 
%primaries to be used in Chapter 7.
%data files obtained from CVRL website
%Data files should be placed in a data directory on the MATLAB path.

clear all, close all;
%read the CIE 1931 XYZ color matching functions.
%they are defined from 360 nm to 830 nm in steps of 1 nm.
xyz = dlmread('ciexyz31_1.txt',',');
lamxyz = xyz(:,1);

%
%compute chromaticities of the monochromatic lights
npts=size(xyz);
S = xyz(:,2:4) * [1 1 1]';
xyzc = [xyz(:,2)./S xyz(:,3)./S xyz(:,4)./S];

%Plotting the color cone
% Determine uniform distribution of wavelengths on the spectral locus
% to help plot the color cone
%Determine the cumulative distance along the spectral locus in the xy 
%chromaticity diagram as a function of wavelength
%interpolate again to .2 nm
total_length(1) = 0;
i=1;
for lambda = 360:829
  i=i+1;
  total_length(i)=total_length(i-1) + sqrt((xyzc(i,1)-xyzc(i-1,1))^2+(xyzc(i,2)-xyzc(i-1,2))^2);
end
total_length = 10*total_length/total_length(i);
%find wavelengths of 20 equally spaced values
total_length_fun =@(x) interp1((360:830),total_length,x,'spline');
lam_equal(1) = 400;
for i = 2:20
    total_length_fun =@(x) interp1((360:830),total_length,x,'spline')-(i-1)*.5;
    lam_equal(i) = round(fzero(total_length_fun,500));
end
lam_equal(21)=700;


%
%Now add the cone of Sony CRT RGB realizable colors
%plot the cone again, as a surface
%Increase the number of wavelengths to 41 to get a smoother plot
%and also the number of radials
%find wavelengths of 41 equally spaced values
total_length_fun =@(x) interp1((360:830),total_length,x,'spline');
lam_equal2(1) = 400;
for i = 2:40
    total_length_fun =@(x) interp1((360:830),total_length,x,'spline')-(i-1)*.25;
    lam_equal2(i) = round(fzero(total_length_fun,500));
end
lam_equal2(41)=700;
%Plot the boundary of the cone of physically realizable colors
%
%Draw a radial line for about 41 wavelengths about equally spaced on the
%shark fin of monochromatic lights. Have these run out to about x_y+z=2 in
%40 equal steps
index_equal2 = lam_equal2-359;
xmesh2 = zeros(41,41);
ymesh2 = zeros(41,41);
zmesh2 = zeros(41,41);
for i=1:41
    alpha = .05*(i-1);
    xmesh2(i,:)= alpha*xyzc(index_equal2,1);
    ymesh2(i,:)= alpha*xyzc(index_equal2,2);
    zmesh2(i,:)= alpha*xyzc(index_equal2,3);
end
figure;

surfl(xmesh2,zmesh2,ymesh2)
view(65,38);
axis square
shading interp
colormap copper;  
hold on
xaxis_X = [0 1.2]; xaxis_Y = [0 0]; xaxis_Z = [0 0];
plot3(xaxis_X, xaxis_Z, xaxis_Y,'k')
yaxis_X = [0 0]; yaxis_Y = [0 2]; yaxis_Z = [0 0];
plot3(yaxis_X, yaxis_Z, yaxis_Y,'k')
zaxis_X = [0 0]; zaxis_Y = [0 0]; zaxis_Z = [0 4];
plot3(zaxis_X, zaxis_Z, zaxis_Y,'k')
text(1.2,0,0,'C_X'); text(0,4,0,'C_Z');text(0,0,2,'C_Y')
xlabel('X'), ylabel('Z'), zlabel('Y')
grid off
axis off
load 'A_RGBCRTtoXYZ';
%Add three edges  -- don't start right at zero and label RGB primaries
E=3.5; E0=0.8 %extent of edges
raxis_X = [E0*A(1,1) E*A(1,1)]; raxis_Y = [E0*A(2,1) E*A(2,1)]; raxis_Z = [E0*A(3,1) E*A(3,1)];
plot3(raxis_X, raxis_Z, raxis_Y,'r')
gaxis_X = [E0*A(1,2) E*A(1,2)]; gaxis_Y = [E0*A(2,2) E*A(2,2)]; gaxis_Z = [E0*A(3,2) E*A(3,2)];
plot3(gaxis_X, gaxis_Z, gaxis_Y,'g')
baxis_X = [E0*A(1,3) E*A(1,3)]; baxis_Y = [E0*A(2,3) E*A(2,3)]; baxis_Z = [E0*A(3,3) E*A(3,3)];
plot3(baxis_X, baxis_Z, baxis_Y,'b')
text(E*A(1,1),E*A(3,1),E*A(2,1)-.15,'[R_C]');
text(E*A(1,2),E*A(3,2),E*A(2,2)+.1,'[G_C]');
text(E*A(1,3),E*A(3,3)+.1,E*A(2,3),'[B_C]');
%Add NC contours on each face
NC = 10;
for con = 1:NC;
r = E*con/NC;    
%RB face    
    rbcon_X = [r*A(1,1) r*A(1,3)]; rbcon_Y = [r*A(2,1) r*A(2,3)];rbcon_Z = [r*A(3,1) r*A(3,3)];
    plot3(rbcon_X, rbcon_Z, rbcon_Y,'k');
%RG face
    rgcon_X = [r*A(1,1) r*A(1,2)]; rgcon_Y = [r*A(2,1) r*A(2,2)];rgcon_Z = [r*A(3,1) r*A(3,2)];
    plot3(rgcon_X, rgcon_Z, rgcon_Y,'k');
%BG face
    bgcon_X = [r*A(1,3) r*A(1,2)]; bgcon_Y = [r*A(2,3) r*A(2,2)];bgcon_Z = [r*A(3,3) r*A(3,2)];
    plot3(bgcon_X, bgcon_Z, bgcon_Y,'k');    
end

set(gcf,'Color',[1 1 1]);

%Plot RGB cube in XYZ axes


figure;

plot3([0 A(1,1)], [0 A(3,1)],[0 A(2,1)], 'k')
view(65,38);
axis square
shading interp
colormap copper;  
hold on
plot3([0 A(1,2)],  [0 A(3,2)],[0 A(2,2)], 'k')
plot3([0 A(1,3)], [0 A(3,3)],[0 A(2,3)],  'k')
plot3([A(1,1) A(1,1)+A(1,2)],  [A(3,1) A(3,1)+A(3,2)],[A(2,1) A(2,1)+A(2,2)], 'k')
plot3([A(1,1) A(1,1)+A(1,3)],  [A(3,1) A(3,1)+A(3,3)],[A(2,1) A(2,1)+A(2,3)], 'k')
plot3([A(1,2) A(1,2)+A(1,3)],  [A(3,2) A(3,2)+A(3,3)], [A(2,2) A(2,2)+A(2,3)],'k')
plot3([A(1,2) A(1,2)+A(1,1)], [A(3,2) A(3,2)+A(3,1)], [A(2,2) A(2,2)+A(2,1)], 'k')
plot3([A(1,3) A(1,3)+A(1,1)],  [A(3,3) A(3,3)+A(3,1)],[A(2,3) A(2,3)+A(2,1)], 'k')
plot3([A(1,3) A(1,3)+A(1,2)], [A(3,3) A(3,3)+A(3,2)], [A(2,3) A(2,3)+A(2,2)], 'k')
plot3([A(1,1)+A(1,2) A(1,1)+A(1,2)+A(1,3)],  [A(3,1)+A(3,2) A(3,1)+A(3,2)+A(3,3)],[A(2,1)+A(2,2) A(2,1)+A(2,2)+A(2,3)], 'k')
plot3([A(1,1)+A(1,3) A(1,1)+A(1,2)+A(1,3)],  [A(3,1)+A(3,3) A(3,1)+A(3,2)+A(3,3)],[A(2,1)+A(2,3) A(2,1)+A(2,2)+A(2,3)], 'k')
plot3([A(1,2)+A(1,3) A(1,1)+A(1,2)+A(1,3)],  [A(3,2)+A(3,3) A(3,1)+A(3,2)+A(3,3)],[A(2,2)+A(2,3) A(2,1)+A(2,2)+A(2,3)], 'k')
xaxis_X = [0 .5]; xaxis_Y = [0 0]; xaxis_Z = [0 0];
plot3(xaxis_X, xaxis_Z, xaxis_Y,'k')
yaxis_X = [0 0]; yaxis_Y = [0 .7]; yaxis_Z = [0 0];
plot3(yaxis_X, yaxis_Z, yaxis_Y,'k')
zaxis_X = [0 0]; zaxis_Y = [0 0]; zaxis_Z = [0 2.6];
plot3(zaxis_X, zaxis_Z, zaxis_Y,'k')
text(.5,0,0.05,'C_X'); text(0,2.7,0,'C_Z');text(0,0.05,.7,'C_Y')
E=1.1;
text(E*A(1,1)+.07,E*A(3,1)-.1,E*A(2,1),'[R_C]');
text(E*A(1,2),E*A(3,2)-.05,E*A(2,2)+.05,'[G_C]');
text(E*A(1,3),E*A(3,3)+.1,E*A(2,3),'[B_C]');
text( A(1,1)+A(1,2)+A(1,3), A(3,1)+A(3,2)+A(3,3)+.05,A(2,1)+A(2,2)+A(2,3),'[W]');
grid off
axis off
set(gcf,'Color',[1 1 1]);