Color Processing

clear all;
close all;

[filename, pathname] = uigetfile({'*.jpg;*.tif;*.png;*.gif;*.bmp','All Image Files';...
          '*.*','All Files' },'mytitle',...
          'C:\Work\myfile.jpg');

img = imread(filename);
[ligne colonne dimension]=size(img);
img=double(img);
R=img(:,:,1);
G=img(:,:,2);
B=img(:,:,3);
 % XYZ to Lab
Rn=R/255;
Gn=G/255;
Bn=B/255;

X = 0.412453*Rn + 0.357580*Gn + 0.180423*Bn;
Y = 0.212671*Rn + 0.715160*Gn + 0.072169*Bn;
Z = 0.019334*Rn + 0.119193*Gn + 0.950227*Bn;
%illuminant D50
% Xref = 0.964296;
% Yref = 1;
% Zref = 0.825105;
%illuminant D65
Xref = 0.950456;
Yref = 1;
Zref = 1.088754;
X=X./Xref;
Y=Y./Yref;
Z=Z./Zref;



% seuil = 0.008856;
% Xp= (X>seuil)  .* (X.^(1/3)) + (X<=seuil)  .* (7.787 .* X + 16/116);
% Yp= (Y>seuil)  .* (Y.^(1/3)) + (Y<=seuil)  .* (7.787 .* Y + 16/116);
% Zp= (Z>seuil)  .* (Z.^(1/3)) + (Z<=seuil)  .* (7.787 .* Z + 16/116);
seuil = (6/29)^3;
f=(1/3) * (29/6)^2;
Xp= (X>seuil)  .* (X.^(1/3)) + (X<=seuil)  .* ( f.* X + 4/29);
Yp= (Y>seuil)  .* (Y.^(1/3)) + (Y<=seuil)  .* ( f.* Y + 4/29);
Zp= (Z>seuil)  .* (Z.^(1/3)) + (Z<=seuil)  .* ( f.* Z + 4/29);

L= 116.* Yp - 16;
a= 500.*(Xp-Yp);
b= 200.*(Yp-Zp);

Ypp= (L + 16) ./ 116;
Xpp= (a./500)+Ypp;
Zpp= Ypp-(b./200);

% seuil = 0.008856;
% Xr =  (Xpp.^3>seuil) .* Xpp.^3 + (Xpp.^3<=seuil) .* ((Xpp - 16/116)/7.787) ;  
% Yr =  (Ypp.^3>seuil) .* Ypp.^3 + (Ypp.^3<=seuil) .* ((Ypp - 16/116)/7.787) ;
% Zr =  (Zpp.^3>seuil) .* Zpp.^3 + (Zpp.^3<=seuil) .* ((Zpp - 16/116)/7.787) ;
seuil = 6/29;
%f= 3 * (6/29)^2;
Xr =  (Xpp>seuil) .* Xpp.^3 + (Xpp<=seuil) .* ((Xpp - 4/29 )/f) ;  
Yr =  (Ypp>seuil) .* Ypp.^3 + (Ypp<=seuil) .* ((Ypp - 4/29 )/f) ;
Zr =  (Zpp>seuil) .* Zpp.^3 + (Zpp<=seuil) .* ((Zpp - 4/29 )/f) ;

Xr=Xr*Xref;
Yr=Yr*Yref;
Zr=Zr*Zref;
 
Rr =  3.240479*Xr - 1.537150*Yr - 0.498535*Zr;
Gr = -0.969256*Xr + 1.875992*Yr + 0.041556*Zr;
Br =  0.055648*Xr - 0.204043*Yr + 1.057311*Zr;

Rr=Rr*255;
Gr=Gr*255;
Br=Br*255;
figure(  'Name','RGB --> XYZ --> Lab --> XYZ -->RGB',...
            'NumberTitle','off',...
            'color',[0.3137 0.3137 0.5098]);
vide=zeros(ligne,colonne);
subplot(4,4,1)
imshow(uint8(cat(3,R,G,B)))
title('Image origine')
subplot(4,4,5)
imshow(uint8(cat(3,R,vide,vide)))
title('Image R origine')
subplot(4,4,9)
imshow(uint8(cat(3,vide,G,vide)))
title('Image G origine')
subplot(4,4,13)
imshow(uint8(cat(3,vide,vide,B)))
title('Image B origine')
subplot(4,4,2)
imagesc(X,[0 1])
colormap('gray')
axis image
title('X')
subplot(4,4,6)
imagesc(Y,[0 1])
colormap('gray')
axis image
title('Y')
subplot(4,4,10)
imagesc(Z,[0 1])
colormap('gray')
axis image
title('Z')
% subplot(4,4,14)
% imshow(uint8(cat(3,vide,vide,Bp)*255))
% title('Image B linéarisé gamma correction')
subplot(4,4,3)
imagesc(Xr,[0 1])
colormap('gray')
axis image
title('X reconstuit')
subplot(4,4,7)
imagesc(Yr,[0 1])
colormap('gray')
axis image
title('Y reconstruit')
subplot(4,4,11)
imagesc(Zr,[0 1])
colormap('gray')
axis image
title('Z reconstruit')
% subplot(4,4,15)
% imshow(uint8(cat(3,vide,vide,Bpp)*255))
% title('Image B reconstuite linéarisé gamma correction')
 subplot(4,4,4)
imshow(uint8(cat(3,Rr,Gr,Br)))
title('Image reconstuite')
 subplot(4,4,8)
imshow(uint8(cat(3,Rr,vide,vide)))
title('Image R reconstuite')
subplot(4,4,12)
imshow(uint8(cat(3,vide,Gr,vide)))
title('Image G reconstuite')
subplot(4,4,16)
imshow(uint8(cat(3,vide,vide,Br)))
title('Image B reconstuite')