% COLOR12 Varous plots for color calculations
%
%
% by Chuck DiMarzio
%    Northeastern Univesity
%    December 2008
%
%
%  !! This file may be copied, used, or modified for educational and
%  !! research purposes provided that this header information is not
%  !! removed or altered, and provided that the book is cited in
%  !! publications, as DiMarzio, Charles A., Optics for Engineers,
%  !! CRC Press, Boca Raton, FL, 2011.
%  !! http://www.crcpress.com
%  !! Other distribution is prohibited without permission.
%
%
% This part from xyztest.m
% Tests xyz.m and shows how close it is to
% accepted data.
% by Chuck DiMarzio, Northeastern University
% 2001
%
set(0,'DefaultAxesColorOrder',[0 0 0]);
w=(300:5:800);
[xbar,ybar,zbar]=xyz(w); % returns functional approximations
load xyzdata;% a(:,1) is wavelength, a(:,2) is xbar, b(:,2) is ybar, c(:,2)
             % is zbar
fig1=figure;
plot(w,ybar,'-',b(:,1)*1000,b(:,2),'.'); % ybar for photometry section
xlabel('Wavelength, nm.');
ylabel('ybar');
grid on;

fig2=figure;     % xyz curves - compare functional approximation to data
plot(w,xbar,'-',w,ybar,'-',w,zbar,'-',...
   a(:,1)*1000,a(:,2),'.',...
   b(:,1)*1000,b(:,2),'.',...
   c(:,1)*1000,c(:,2),'.');
xlabel('Wavelength, nm.');
ylabel('xbar,ybar,zbar');
text(480,1.4,'z');
text(600,1.2,'x');
text(440,0.45,'x');
text(550,1.2,'y');
grid on;


fig3=figure;     % xyz for lasers
lambdal=[457.9,530.9,647.1]; % Argon Krypton Krypton
[xl,yl,zl]=xyz(lambdal);
plot(lambdal,xl,'x',...
     lambdal,yl,'v',...
     lambdal,zl,'<',...
     w,xbar,'-',w,ybar,'-',w,zbar,'-',...
     lambdal(1)*[1,1],[0,2],'-',...
     lambdal(2)*[1,1],[0,2],'-',...
     lambdal(3)*[1,1],[0,2],'-')
xlabel('Wavelength, nm.');
ylabel('xbar,ybar,zbar');
text(480,1.4,'z');
text(600,1.2,'x');
text(440,0.45,'x');
text(550,1.2,'y');
legend('xbar','ybar','zbar');
grid on;

%
% Compute color vectors, x,y,z,
cbluelaser=[xl(1),yl(1),zl(1)];
cgreenlaser=[xl(2),yl(2),zl(2)];
credlaser=[xl(3),yl(3),zl(3)];
disp('Red, Green, Blue Lasers in Columns, rows are x,y,z');
lasermatrix=[credlaser',cgreenlaser',cbluelaser']
% and normalize for later
cbluelaser=cbluelaser/sum(cbluelaser);
cgreenlaser=cgreenlaser/sum(cgreenlaser);
credlaser=credlaser/sum(credlaser);

% Black body three temperatures
tlow=2000;ttungsten=3000;tqh=3500;
dlambdab=1;lambdab=(300:dlambdab:800); % nanometers
mlow=bbspec(tlow,lambdab/1000); % Low  % Convert wavelength to microns
mt=bbspec(ttungsten,lambdab/1000);   % t for Tungsten.
mqh=bbspec(tqh,lambdab/1000);   % qh for Quartz Halogen
[xbar,ybar,zbar]=xyz(lambdab);
xlow=xbar.*mlow;Xlow=sum(xlow)*dlambdab;
ylow=ybar.*mlow;Ylow=sum(ylow)*dlambdab;
zlow=zbar.*mlow;Zlow=sum(zlow)*dlambdab;
tristimlow=[Xlow,Ylow,Zlow]';
fig4=figure;subplot(3,1,1);
semilogy(lambdab,xlow,'-',lambdab,ylow,'-.',lambdab,zlow,'--',...
         lambdab,mlow,':');
axis([300,800,1e3,1e5]);xlabel('\lambda, Wavelength, nm');
ylabel('xs,ys,zs');
clow=[Xlow,Ylow,Zlow]'
xt=xbar.*mt;Xt=sum(xt)*dlambdab;
yt=ybar.*mt;Yt=sum(yt)*dlambdab;
zt=zbar.*mt;Zt=sum(zt)*dlambdab;
tristimt=[Xt,Yt,Zt]';
subplot(3,1,2);
semilogy(lambdab,xt,'-',lambdab,yt,'-.',lambdab,zt,'--',...
         lambdab,mt,':');
axis([300,800,1e5,1e7]);xlabel('\lambda, Wavelength, nm');
ylabel('xs,ys,zs');
ct=[Xt,Yt,Zt]'
xqh=xbar.*mqh;Xqh=sum(xqh)*dlambdab;
yqh=ybar.*mqh;Yqh=sum(yqh)*dlambdab;
zqh=zbar.*mqh;Zqh=sum(zqh)*dlambdab;
tristimqh=[Xqh,Yqh,Zqh]';
subplot(3,1,3);
semilogy(lambdab,xqh,'-',lambdab,yqh,'-.',lambdab,zqh,'--',...
         lambdab,mqh,':');
axis([300,800,1e5,1e7]);xlabel('\lambda, Wavelength, nm');
ylabel('xs,ys,zs');

cqh=[Xqh,Yqh,Zqh]'

% Normalize the hot objects
xplow=Xlow/(Xlow+Ylow+Zlow);
xpt=Xt/(Xt+Yt+Zt);
xpqh=Xqh/(Xqh+Yqh+Zqh);
yplow=Ylow/(Xlow+Ylow+Zlow);
ypt=Yt/(Xt+Yt+Zt);
ypqh=Yqh/(Xqh+Yqh+Zqh);

% Find the colorspace of the rgb colors used.
matrix_xyz2rgb=xyz2rgb(0);
%matrix_rgb2xyz=inv(matrix_xyz2rgb);
matrix_rgb2xyz=eye(size(matrix_xyz2rgb))/matrix_xyz2rgb;
xyzfromrgb=matrix_rgb2xyz*(eye(3));
xyzfromrgb=xyzfromrgb./repmat(sum(xyzfromrgb,1),3,1); % Check normalization

% Get file of xy data vs temperature
load hotcolors;

%  Find the boundary of the xy plot, to plot all color data on xy
denom=(a(:,2)+b(:,2)+c(:,2));
xbound=a(:,2)./denom;
ybound=b(:,2)./denom;
fig5=figure;plot(xplow,yplow,'*',...
            xpt,ypt,'+',...
            xpqh,ypqh,'o',...
            cbluelaser(1),cbluelaser(2),'<',...
            cgreenlaser(1),cgreenlaser(2),'^',...
            credlaser(1),credlaser(2),'>',...
            xyzfromrgb(1,2),xyzfromrgb(2,2),'s',...
            xbound,ybound,'-',...
            xhot,yhot,'--',...
            [cbluelaser(1),credlaser(1),cgreenlaser(1),cbluelaser(1)],...
            [cbluelaser(2),credlaser(2),cgreenlaser(2),cbluelaser(2)]);
axis([-0.1,1,-0.1,1]);
xlabel('x');ylabel('y');
axis equal;
hold on;
plot(xyzfromrgb(1,[1:3,1]),xyzfromrgb(2,[1:3,1]),'-')
legend([num2str(tlow),'K'],[num2str(ttungsten),'K'],...
       [num2str(tqh),'K'],...
       ['\lambda=',num2str(lambdal(1)),' nm'],...
       ['\lambda=',num2str(lambdal(2)),' nm'],...
       ['\lambda=',num2str(lambdal(3)),' nm']);
hold off;

% Work in progress to generate the chromaticity coordinates
[x,y]=meshgrid(0:0.001:1,0:0.001:1);z=1-x-y;
axyz(:,:,1)=x;axyz(:,:,2)=y;axyz(:,:,3)=z;
rgb=xyz2rgb(axyz);

c=rgb/max(max(max(rgb)));
c(c<0)=0;
csave=c;
c=c./repmat(max(c,[],3),[1,1,3]);
% Generate unmask which is 1 outside the boundary and 0 inside
%  This is done by checking the cross producct of a vector from each
%  boundary point to the next with a vector from that boundary point to
%  a point in the image, and unmasking if negative.
unmask=(zeros(size(x))==1);
nstart=7;nstop=64;
for n=nstart:nstop;
x1=xbound(n);y1=ybound(n);x2=xbound(n+1);y2=ybound(n+1);
cp=(x-x1)*(y2-y1)-(y-y1)*(x2-x1);
test=cp<0;
unmask=unmask|test;
end;
x1=xbound(nstop);y1=ybound(nstop);x2=xbound(nstart);y2=ybound(nstart);
cp=(x-x1)*(y2-y1)-(y-y1)*(x2-x1);
test=cp<0;
unmask=unmask|test;
% Whiten everything that is outside the boundary
for n=1:3;
c(:,:,n)=c(:,:,n).*(1-unmask)+unmask;
end;
fig6=figure;imagesc(x(1,:),y(:,1),c);axis xy;axis equal;axis image;
hold on;
plot(xbound,ybound,'w');
xlabel('x');ylabel('y');
hold off;

%  The big color plot
fig7=figure;imagesc(x(1,:),y(:,1),c);axis xy;axis equal;axis image;
hold on;
plot(xhot,yhot,'w-',...
     xyzfromrgb(1,[1:3,1]),xyzfromrgb(2,[1:3,1]),'w-')
h=plot(xplow,yplow,'k*',...
     xpt,ypt,'k+',...
     xpqh,ypqh,'ko',...
          cbluelaser(1),cbluelaser(2),'k<',...
     cgreenlaser(1),cgreenlaser(2),'k^',...
     credlaser(1),credlaser(2),'k>',...
     xyzfromrgb(1,2),xyzfromrgb(2,2),'ks',...
     xbound,ybound,'k:',...
     [cbluelaser(1),credlaser(1),cgreenlaser(1),cbluelaser(1)],...
     [cbluelaser(2),credlaser(2),cgreenlaser(2),cbluelaser(2)],'k--');
xlabel('x');ylabel('y');
legend(h,...
       [num2str(tlow),'K'],[num2str(ttungsten),'K'],...
       [num2str(tqh),'K'],...
       ['\lambda=',num2str(lambdal(1)),' nm'],...
       ['\lambda=',num2str(lambdal(2)),' nm'],...
       ['\lambda=',num2str(lambdal(3)),' nm'],...
       'RGB')

%  Computing the rgb power for 4000 Lumens of white light with the lasers
%  described above.
%rgb_for_white=inv(lasermatrix)*[1,1,1]'/3*3*4000/683
rgb_for_white=eye(size(lasermatrix))/lasermatrix*[1,1,1]'/3*3*4000/683

% Computing the xyz for a 488-nm laser
[X488,Y488,Z488]=xyz(488)
%rgb488=inv(lasermatrix)*[X488,Y488,Z488]'
rgb488=eye(size(lasermatrix))/lasermatrix*[X488,Y488,Z488]'

rgb2xyzwikip=[0.49,0.31,0.20;0.17697,0.81240,0.01063;0.00,0.01,0.99]/ ...
    0.17697;
% http://en.wikipedia.org/wiki/CIE_1931_color_space
% Fairman H.S., Brill M.H., Hemmendinger H. (February 1997). "How the CIE
% 1931 Color-Matching Functions Were Derived from the
% Wright-Guild Data". Color Research and Application 22
% (1):
%
% doi:10.1002/(SICI)1520-6378(199702)22:1<11::AID-COL4>3.0.CO;2-7.  and
% Fairman H.S., Brill M.H., Hemmendinger H. (August 1998). "Erratum: How
% the CIE 1931 Color-Matching Functions Were Derived from the
% Wright--Guild Data". Color Research and Application 23
% (4): 259.
Red, Green, Blue Lasers in Columns, rows are x,y,z
lasermatrix =
    0.3369    0.1709    0.2890
    0.1344    0.8313    0.0313
    0.0000    0.0347    1.6962
clow =
   1.0e+06 *
    2.7170
    2.1592
    0.3101
ct =
   1.0e+08 *
    1.5087
    1.3908
    0.5480
cqh =
   1.0e+08 *
    4.8600
    4.6534
    2.4384
Warning: Ignoring extra legend entries. 
rgb_for_white =
   11.9840
    4.9821
    3.3507
X488 =
    0.0605
Y488 =
    0.2112
Z488 =
    0.5630
rgb488 =
   -0.2429
    0.2811
    0.3261