## Copyright (C) 2012   Ming Liu   <u5221472@gmail.com>
##
## This program is free software; you can redistribute it and/or modify
## it under the terms of the GNU General Public License as published by
## the Free Software Foundation; either version 3 of the License, or
## (at your option) any later version.
##
## This program is distributed in the hope that it will be useful,
## but WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
## GNU General Public License for more details.
##
## You should have received a copy of the GNU General Public License
## along with this program; if not, see <http://www.gnu.org/licenses/>.

function Ilab = rgb2lab(I)

% RGB -> XYZ
I = double(I);
R = I(:,:,1)/255;
G = I(:,:,2)/255;
B = I(:,:,3)/255;

R(R > 0.03928 ) = ((R(R > 0.03928 ) + 0.055 ) ./ 1.055 ) .^ 2.4;
G(G > 0.03928 ) = ((G(G > 0.03928 ) + 0.055 ) ./ 1.055 ) .^ 2.4;
B(B > 0.03928 ) = ((B(B > 0.03928 ) + 0.055 ) ./ 1.055 ) .^ 2.4;

R(R <= 0.03928 ) = R(R <= 0.03928 ) ./ 12.92;
G(G <= 0.03928 ) = G(G <= 0.03928 ) ./ 12.92;
B(B <= 0.03928 ) = B(B <= 0.03928 ) ./ 12.92;

X = R .* 0.4124 + G .* 0.3576 + B .* 0.1805;
Y = R .* 0.2126 + G .* 0.7152 + B .* 0.0722;
Z = R .* 0.0193 + G .* 0.1192 + B .* 0.9505;

%XYZ -> LAB
xyz = [0.9504, 1.0000, 1.0889]; %X/Xn Y/Yn Z/Zn
X = X./xyz(1);
Y = Y./xyz(2);
Z = Z./xyz(3);

%Calculate L
L = Y;
L(L > 0.008856) = 116 * L(L > 0.008856).^(1/3) - 16;
L(L<= 0.008856) = 903.3 * L(L<= 0.008856) ;

X(X > 0.008856) = X(X > 0.008856) .^ (1/3);
Y(Y > 0.008856) = Y(Y > 0.008856) .^ (1/3);
Z(Z > 0.008856) = Z(Z > 0.008856) .^ (1/3);

X(X<= 0.008856) = 7.787 * X(X <= 0.008856) + 4/29;
Y(Y<= 0.008856) = 7.787 * Y(Y <= 0.008856) + 4/29;
Z(Z<= 0.008856) = 7.787 * Z(Z <= 0.008856) + 4/29;

a = 500 * (X - Y) + 128; a = uint8(a);
b = 200 * (Y - Z) + 128; b = uint8(b);
L = L*255/100; L = uint8(L);

Ilab = uint8(I);
Ilab(:,:,1) = L;
Ilab(:,:,2) = a;
Ilab(:,:,3) = b;

end


%simple colorspace transformation, rgb to lab
%usage: I = imread('test.png');
%Ilab = rgb2lab(I);
%imshow(Ilab);
