function xyz = lab2xyz(lab,white)
% 
%  xyz = lab2xyz(lab,white)
%
%    Return the xyz values from a set of La*b* values.
% 
% ARUGMENTS:
%    lab is a 3 x N matrix of Lab values
%    white is 3 vector specifying the white point
%
% 
% 01.19.98 BW.  Modified by inserting new tests (find, isempty)
% for Matlab 5.0 syntax.  The old strategy produced errors under 5.0

% DEBUGGING
% lab = grayLAB;
% white = whiteXYZ;

if length(white) ~= 3,
  error('White point is not a three-vector');
elseif size(lab,1) ~= 3,
  error('lab argument must have three rows');
end

Lstar = lab(1,:);
astar = lab(2,:);
bstar = lab(3,:);

Xn = white(1);
Yn = white(2);
Zn = white(3);


%
%	Allocate space
%
fX = zeros(1,size(lab,2));
fY = zeros(1,size(lab,2));
fZ = zeros(1,size(lab,2));

Xratio = zeros(1,size(lab,2));
Yratio = zeros(1,size(lab,2));
Zratio = zeros(1,size(lab,2));


% Usual formula for Lstar
Y = Yn .*  ( (Lstar + 16.0) ./ 116.0) .^ 3.0;

% but if Y/Yn is too small, use other formula
l = (Y/Yn) < 0.008856;
l = find(l > 0) ;
if ~isempty(l)
  Y(l) =  (Lstar(l) / 903.3) * Yn;
end
%
%	Now, we estimate X/Xn, Y/Yn, and Z/Zn
%
Yratio = (Y/Yn);
l = find(Yratio > .008856);
if ~isempty(l)
  fY(l) = Yratio(l).^(1/3);
end

l = find(Yratio <= .008856);
if ~isempty(l)
  fY(l) = 7.787 * Yratio(l) + 16/116;
end


% Inverse range for Lab ratios
fX = (astar/500.0) + fY;

l = find(fX > 0.206893);
if ~isempty(l)
  Xratio(l) = fX(l) .^ 3;
end

l = find(fX <= 0.206893);
if ~isempty(l)
  Xratio(l) = (fX(l) - (16/116)) / 7.787;
end

% Inverse range for Lab ratios
fZ = fY - (bstar/200.0);

l = find(fZ > 0.206893);
if ~isempty(l)
  Zratio(l) = fZ(l) .^ 3;
end
l = find(fZ <= 0.206893);
if ~isempty(l)
  Zratio(l) = (fZ(l) - (16/116)) / 7.787;
end

%	Finally, we deduce X, Y and Z
X = Xratio * Xn;
Z = Zratio * Zn;

xyz = [X;Y;Z];