% N = rgb2luv(M) Colour transformation from RGB to CIE L*u*v*
%
% Reference:
%
%   William K. Pratt: "Digital Image Processing", 3rd ed. pp. 64ff
%
% Ranges for L*u*v* values:
%   L*: [        0, 100.0397]
%   u*: [-132.0772, 220.9554]
%   v*: [-139.1071, 121.4934]
%
% In:
%   M: Nx3 array containing rows of RGB values in [0,1].
% Out:
%   N: Nx3 array containing rows of CIE L*u*v* values.

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

function N = rgb2luv(M)

% Deal with black (RGB = [0 0 0]) separately to avoid division-by-zero
tmp = find(M(:,1)==0 & M(:,2)==0 & M(:,3)==0); M(tmp,:) = 1;

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

% CIE L*u*v*
L = f1(Y/y);
u = 4*X./(X+15*Y+3*Z); u0 = 4*x./(x+15*y+3*z); u = 13*L.*(u-u0);
v = 9*Y./(X+15*Y+3*Z); v0 = 9*y./(x+15*y+3*z); v = 13*L.*(v-v0);

L(tmp) = 0; u(tmp) = 0; v(tmp) = 0;	% Should L be -16 for black??

N = [L u v];


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);