invert=solve('u=(0.4661*x+0.1593*y)/(y-0.15735*x+0.2424)','v=0.6581*y/(y-0.15735*x+0.2424)')
% you should get the following
% x=4848.*(-1593.*v+6581.*u)/(-98233171.*v+20710407.*u+61348082.)
% y=22596528.*v/(-98233171.*v+20710407.*u+61348082.)
N=1000; % size of the image
a=0; b=1.0; c=0; d=0.6;
xaxis=linspace(a,b,N*min((b-a)/(d-c),1));% area of interest
yaxis=linspace(c,d,N*min((d-c)/(b-a),1));% area of interest
[u,v,Y]=meshgrid(xaxis,yaxis,50);
x=4848.*(-1593.*v+6581.*u)./(-98233171.*v+20710407.*u+61348082.);
y=22596528.*v./(-98233171.*v+20710407.*u+61348082.);
xyz=applycform(cat(3,x,y,Y),makecform('xyl2xyz'));
M=max(xyz,[],3);
srgb=applycform(xyz./cat(3,M,M,M),makecform('xyz2srgb'));
image(xaxis,yaxis,srgb);
axis image;
set(gca,'YDir','normal'); % right way up
hold;
[LAMBDA, XFCN, YFCN, ZFCN]=colorMatchFcn('1931_full');
sl=applycform(cat(3,XFCN,YFCN,ZFCN),makecform('xyz2xyl'));
sluv=sl;
sluv(:,:,1)=(0.4661*sl(:,:,1)+0.1593*sl(:,:,2))./(sl(:,:,2)-0.15735*sl(:,:,1)+0.2424);
sluv(:,:,2)=0.6581*sl(:,:,2)./(sl(:,:,2)-0.15735*sl(:,:,1)+0.2424);
sl_loopx=[sluv(:,:,1) sluv(:,1,1)]; % make a closed polygon
sl_loopy=[sluv(:,:,2) sluv(:,1,2)]; % make a closed polygon
set(patch([b a a b sl_loopx b],[d d c c sl_loopy c],[0.5 0.5 0.5]),'EdgeColor','none');
plot(sl_loopx,sl_loopy,'k-'); % k- means solid (-) black (k) line
axis([a b c d]); xlabel('u'); ylabel('v'); shg; grid;
[lambda,T]=meshgrid(LAMBDA,[1e3:100:1.1e4 2.5e4 5e4 1e5 5e5 1e6]);
spd=lambda.^(-5)./(exp(1.4388e+7./(lambda.*T))-1);
XYZbb=spd*[XFCN' YFCN' ZFCN'];
xybb=applycform(XYZbb./(XYZbb(:,2)*[1 1 1]),makecform('xyz2xyl'));
uvbb=xybb;
uvbb(:,1)=(0.4661*xybb(:,1)+0.1593*xybb(:,2))./(xybb(:,2)-0.15735*xybb(:,1)+0.2424);
uvbb(:,2)=0.6581*xybb(:,2)./(xybb(:,2)-0.15735*xybb(:,1)+0.2424);
plot(uvbb(:,1),uvbb(:,2),'Color',[0.5 0.5 0.5],'LineWidth',1.5);
xyz_wp=[ [1 1 1] ; ...
whitepoint('a') ; ...
whitepoint('c') ; ...
whitepoint('d50') ; ...
whitepoint('d55') ; ...
whitepoint('d65') ];
xyl_wp=applycform(xyz_wp,makecform('xyz2xyl'));
uvl_wp=xyl_wp;
uvl_wp(:,1)=(0.4661*xyl_wp(:,1)+0.1593*xyl_wp(:,2))./(xyl_wp(:,2)-0.15735*xyl_wp(:,1)+0.2424);
uvl_wp(:,2)=0.6581*xyl_wp(:,2)./(xyl_wp(:,2)-0.15735*xyl_wp(:,1)+0.2424);
plot(uvl_wp(:,1),uvl_wp(:,2),'xk');
slope=diff(uvbb(:,[1 2]));
for i=1:10:100
xi=uvbb(i,1);
yi=uvbb(i,2);
a=-slope(i,1)/slope(i,2);
l=0.05/sqrt(1+a^2);
% set the line to the color of the black body (poor legibility):
% plot([xi-l xi+l],[yi-a*l yi+a*l],'Color',applycform(XYZbb(i,:)./(XYZbb(i,2)*[1 1 1]),makecform('xyz2srgb')),'LineWidth',2);
plot([xi-l xi+l],[yi-a*l yi+a*l],'Color',[0.5 0.5 0.5]);
end
