%Chapter 7,  Fig 7.7
%Eric Dubois, updated 2018-12-13
%Draw the cone of physically realizable colors in the XYZ coordinate
%system.
% 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',',');

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

%Plot the boundary of the cone of physically realizable colors
%
%Draw a radial line for about 20 wavelengths about equally spaced on the
%shark fin of monochromatic lights. Have these run out to about x_y+z=2 in
%10 equal steps
index_equal = lam_equal-359;
xmesh = zeros(11,21);
ymesh = zeros(11,21);
zmesh = zeros(11,21);
for i=1:11
    alpha = .2*(i-1);
    xmesh(i,:)= alpha*xyzc(index_equal,1);
    ymesh(i,:)= alpha*xyzc(index_equal,2);
    zmesh(i,:)= alpha*xyzc(index_equal,3);
end
figure;

mesh(xmesh,zmesh,ymesh)
view(65,38);
axis square
colormap([0 0 0]);    
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 2.5];
plot3(zaxis_X, zaxis_Z, zaxis_Y,'k')
text(1.2,0,0,'C_X'); text(0,2.5,0,'C_Z');text(0,0,2,'C_Y')
xlabel('X'), ylabel('Z'), zlabel('Y')
%Add the XYZ primaries
[SPX,SPZ,SPY] = sphere(10);
sc = 0.03;
SPXs = sc*SPX+1; SPYs = sc*SPY; SPZs=sc*SPZ;
surf(SPXs,SPZs,SPYs)
text(1,0-.3,-.1,'[X]');
SPXs = sc*SPX; SPYs = sc*SPY; SPZs=sc*SPZ+1;
surf(SPXs,SPZs,SPYs)
text(0,-0.5,1,'[Y]');
SPXs = sc*SPX; SPYs = sc*SPY+1; SPZs=sc*SPZ;
surf(SPXs,SPZs,SPYs)
%text(0,1,-.1,'[Z]')
shading interp
colormap(gray)
grid off
axis off
set(gcf,'Color',[1 1 1]);