% N = rgb2lab(M) Colour transformation from RGB to CIE L*a*b*
%
% Reference:
%
%   William K. Pratt: "Digital Image Processing", 3rd ed. pp. 64ff
%
% Ranges for L*a*b* values:
%   L*: [        0, 100.0397]
%   a*: [-137.8664,  96.1144]
%   b*: [ -78.4155, 142.8521]
%
% In:
%   M: Nx3 array containing rows of RGB values in [0,1].
% Out:
%   N: Nx3 array containing rows of CIE L*a*b* values.

% Copyright (c) 2005 by Miguel A. Carreira-Perpinan

function N = rgb2lab(M)

% CIE XYZ
[X Y Z x y z] = rgb2xyz(M);

% CIE L*a*b*
L = f1(Y/y); a = 500*(f2(X/x)-f2(Y/y)); b = 200*(f2(X/x)-f2(Z/z));

N = [L a b];


function y = f1(x)

z1 = find(x>0.008856); z2 = setdiff(1:length(x),z1);
y = x; y(z1) = 25*(100*x(z1)).^(1/3)-16; y(z2) = 903.3*x(z2);


function y = f2(x)

z1 = find(x>0.008856); z2 = setdiff(1:length(x),z1);
y = x; y(z1) = x(z1).^(1/3); y(z2) = 7.787*x(z2)+0.1379;