%Cargar las funciones de igualación CIEXYZ-1931 (ciexyzn.mat)
[XYZ_l,Yxw,Mxx]=loadsys;
%Cargamos los datos del monitor
[tm,a,g]=loadmon(Mxx);

%Definimos un conjunto de colores
XYZ=[61  105   30;91   87  156;81   39  144]
[n_col,sat,XYZp]=tri2val(XYZ,Yxw,tm,a,g,8);
im=(size(n_col,1)+1)*ones(3,2*size(n_col,1)+1);
im(2,2:2:(2*size(n_col,1)+1))=1:size(n_col,1);
colormap([n_col;0 0 0]);image(im);axis('off');

%Representamos estos colores, y el blanco CIEXYZ, en el
%diagrama cromático:

lambda=XYZ_l([41 101 111:5:251 261:20:301],1);
figure;colordgm(XYZ_l([41 101 111:5:251 261:20:301],2:4),1,XYZ_l,Yxw,'show_numb',1)
ch=get(gca,'Children');
for i=1:size(lambda)
   set(ch(i),'string',num2str(lambda(length(lambda)-i+1)))
end
hold on;colordgm([1 1 1;XYZ],1,XYZ_l,Yxw,'show_numb',1)

%El dominante de un color C es aquel color espectral o de la recta 
%de los púrpuras, que, mezclado con blanco en las proporciones adecuadas,
%iguala la apariencia de C
%
%Para obtener gráficamente el dominante, trazamos la recta que va desde el
%blanco al color C, hasta cortar al locus espectral o a la recta de los púrpuras y
%medimos las coordenadas cromáticas del dominante

xyY=tri2coor([1 1 1;XYZ],Yxw);
m=(xyY(2:end,2)-xyY(1,2))./(xyY(2:end,1)-xyY(1,1));
n=xyY(2:end,2)-m.*xyY(2:end,1);
x=linspace(0,0.73);
hold on
title('Selecciona el dominante');
for i=1:size(XYZ,1)
   h(i)=plot(x,m(i)*x+n(i),'r-');
   [xD(i),yD(i)]=ginput(1);
   set(h(i),'Color',[0.8 0.8 0.8])
   plot(xD(i),yD(i),'o');
   tx(i)=text(xD(i)+0.01,yD(i),['(' num2str(xD(i),3) ',' num2str(yD(i),3) ')'])
   set(tx,'Color',[1 0 0]);
end

%Conocidas las coordenadas del dominante, podemos calcular su longitud de onda o la
%longitud de onda del complementario.

xyD=[xD' yD'];
xy_l=tri2coor(XYZ_l(:,2:4),Yxw);

%Para ello, determinamos a qué intervalo de longitudes de onda pertenece la dominante
%y calculamos la longitud de onda dominante por interpolación.
%Si el color pertenece a la región púrpura, no encontraremos su dominante mediante este
%proceso.

lD=zeros(size(xyY,1)-1,1);
for j=1:size(XYZ,1)
 for i=1:5:size(XYZ_l(:,1))-5
    dentro=inpolygon(xyY(j+1,1),xyY(j+1,2),[xy_l(i,1) 1/3 xy_l(i+5,1)],[xy_l(i,2) 1/3  xy_l(i+5,2)]);
    if dentro>0 
       lD(j)=interp1(atan2(xy_l(i:i+5,2),xy_l(i:i+5,1)),XYZ_l(i:i+5,1),atan2(xyD(j,2),xyD(j,1)));
       break
    end
 end
end
noespurp=find(lD>0)'
for i=noespurp
   set(tx(i),'String',[num2str(lD(i),4) ' nm']);
end

%Si la dominante es púrpura, calculamos el color complementario, esto es, el estímulo
%espectral que sumado al color C iguala al blanco.
%
espurp=find(lD==0)';
for i=espurp
   set(h(i),'Color',[1 0 0])
   title('Selecciona el complementario');
   [xC(i),yC(i)]=ginput(1);
   xyC(i,:)=[xC(i) yC(i)];
   set(h(i),'Color',[0.8 0.8 0.8])
   plot(xC(i),yC(i),'o');
   tx(i)=text(xC(i)+0.01,yC(i),['(' num2str(xC(i),3) ',' num2str(yC(i),3) ')'])
   set(tx(i),'Color',[0 0 1]);
end
for j=espurp
   xxyyY(j,:)=tri2coor([1 1 1]+coor2tri([xyC(j,:) 1],Yxw),Yxw);%Obtenemos un color entre el blanco y el complementario
   for i=1:5:size(XYZ_l(:,1))-5
    dentro=inpolygon(xxyyY(j,1),xxyyY(j,2),[xy_l(i,1) 1/3 xy_l(i+5,1)],[xy_l(i,2) 1/3  xy_l(i+5,2)]);
    if dentro>0 
       lD(j)=-interp1(atan2(xy_l(i:i+5,2),xy_l(i:i+5,1)),XYZ_l(i:i+5,1),atan2(xyC(j,2),xyC(j,1)));
       break
    end
 end
 set(tx(j),'String',[num2str(lD(j),3) ' nm'])
end

%Calculamos la pureza de excitación, pe=d(W,C)/d(W,D)
dWC=sqrt(sum(((xyY(2:end,1:2)-xyY(ones(size(xyY,1)-1,1),1:2)).^2)'));
dWD=sqrt(sum(((xyD-xyY(ones(size(xyY,1)-1,1),1:2)).^2)'));
pe=dWC./dWD

%Calculamos la pureza colorimetrica pc=YD/YW
%Para ello, resolvemos la ecuación siguiente:
%TC=TD+TW=YD*tD/(Yw*tD')+YW*tW/(Yw*tW')
%Esto es un sistema lineal de tres ecuaciones con 2 incógnita
%Calculamos la matriz del sistema de ecuaciones, para cada color
xyzD=[xyD 1-sum(xyD')'];
xyzW=[xyY(1,1:2) 1-sum(xyY(1,1:2))];
for i=1:size(XYZ,1)
   M(:,:,i)=[xyzD(i,:)'/(Yxw*xyzD(i,:)') xyzW'/(Yxw*xyzW')];
end

%inspeccionamos las matrices M(:,:,i) hasta encontrar un determinante
%2x2 no nulo y resolvemos el sistema de ecuaciones
com=[1 2;1 3;2 3];
no=[3 2 1];
for i=1:size(XYZ,1)
   for cual=1:3
      if abs(det(M(com(cual,:),:,i)))>0
         YDW(:,i)=inv(M(com(cual,:),:,i))*XYZ(i,com(cual,:))';
         cual_no(i)=no(cual);
         break;
      end
   end
end
%Comprobamos que el resultado es compatible con la ecuación que no hemos usado
for i=1:size(XYZ,1)
   [XYZ(i,cual_no(i)) M(cual_no(i),:,i)*YDW(:,i)]
   [XYZ(i,cual_no(i))-M(cual_no(i),:,i)*YDW(:,i)]
end   
   
%Calculamos la pureza colorimétrica
pc=YDW(1,:)'./xyY(2:end,3)

%Comprobamos la relación entre ambos tipos de pureza
xyz=[xyY(2:end,1:2) 1-sum(xyY(2:end,1:2)')']
for i=1:size(XYZ)
   [pc(i)  pe(i)*(xyzD(i,:)*Yxw')/(xyz(i,:)*Yxw')]
end

%Si resolvemos el problema analítica y no gráficamente

lpeY=coor2lp(xyY(2:end,:),1,XYZ_l,Yxw)
lpcY=coor2lp(xyY(2:end,:),2,XYZ_l,Yxw)

%Comparamos los resultados obtenidos en CIEXYZ con los que se obtienen 
%en CIERGB

%Cargamos las funciones de igualación de RBG (ciergb_31)
[RGB_l,Yrw,Mrx]=loadsys;

%Transformamos los colores del problema a este espacio

RGB=newbasis(XYZ,inv(Mrx));

%Calculamos las coordenadas (r,g) cromáticas de estos colores y su
%luminancia

rgY=tri2coor(RGB,Yrw);

%Calculamos dominante, pureza y luminancia
lpeY_r=coor2lp(rgY,1,RGB_l,Yrw)
lpcY_r=coor2lp(rgY,2,RGB_l,Yrw)


%Al comparar con las obtenidas en CIEXYZ, vemos que la dominante, 
%la pureza colorimétrica y la luminancia no cambian al cambiar de espacio,
%pero la pureza de excitación sí. ¿Por qué?
%
%La dominante y la pureza colorimétrica no cambian si los dos espacios
%a) son una transformación lineal el uno del otro y b) tienen el mismo blanco
%de referencia.

