function Ljg = XYZToLjg(XYZ)
% Ljg = XYZToLjg(XYZ)
%
% Convert XYZ (10 degree) to OSA Ljg coordinates.  Formulae
% derived from MacAdam (1974, JOSA, 64, pp. 1691-1702).  Note that
% MacAdam's formulae are in error for the case Y0 less than 30.  The
% correct formulae may be deduced by examing Semelroth's (JOSA, 1970, 16,
% 1685-1689) Eqn. 3, from which MacAdam's incomplete Eqn. 1 is derived.
%
% Note that the above problem is propogated into the formulae
% reported in Wyszecki and Stiles, 2cd edition.  In addition, W+S
% reverse the formulae for G and J and leave out the final
% transformation from scriptL to L.
%
% Finally the formualae published in Brainard (2003, Color appearance and
% color difference specification. In The Science of Color, 2cd edition,
% S. K. Shevell (ed.), Optical Society of America, Washington D.C., 191-216),
% which correct for the above errors, introduce a new mistake: The coefficient
% on B^1/3 for the j coordinate in Eq. 5-2 should be -9.7, rather than the
% the +9.7 that is published.
%
% The output of this routine was verified against the tabulated
% values in W+S, Table I(6.6.4).  These are republished from
% MacAdam (1978, JOSA, 68, 121-130).  See OSAUCSTest.
%
% 3/27/01  dhb  Wrote it.
% 7/14/10  dhb  Added comment that the formulae in my chapter have a typo.
%               The code here is and was correct.

% Define XYZToRGB matrix.
M_XYZToRGB = [0.799 0.4194 -0.1648 ; 
						 -0.4493 1.3265 0.0927 ;
						 -0.1149 0.3394 0.7170];
RGB = M_XYZToRGB*XYZ;
RGB3 = RGB.^(1/3);

% Compute xyY from XYZ
xyY = XYZToxyY(XYZ);

% Compute Y0
x = xyY(1,:);
y = xyY(2,:);
Y = xyY(3,:);
Y0 = Y.* ...
				  (4.4934*(x.^2)+4.3034*(y.^2)-4.276*(x.*y) ...
			    -1.3744*x - 2.5643*y + 1.8103);

% Compute scriptL.  Note that MacAdam does not correctly
% handle the case of Y0 < 30.
scriptL = zeros(size(Y0));
index = find(Y0 > 30);
if (~isempty(index))
	scriptL(index) = 5.9 * ((Y0(index).^(1/3))-(2/3)+0.042*((abs(Y0(index)-30)).^(1/3)));
end
index = find(Y0 <= 30);
if (~isempty(index))
	scriptL(index) = 5.9 * ((Y0(index).^(1/3))-(2/3)-0.042*((abs(Y0(index)-30)).^(1/3)));
end

% Compute C.  Use version that depends on scriptL, as I'm not sure
% the alternate version is correct for Y0 < 30 (I didn't check).
C = scriptL./(5.9*((Y0.^(1/3))-(2/3)));

% Compute L,g,j.
Ljg = zeros(size(XYZ));
Ljg(1,:) = (scriptL-14.4)/sqrt(2);
Ljg(2,:) = C.*([1.7 8 -9.7]*RGB3);
Ljg(3,:) = C.*([-13.7 17.7 -4]*RGB3);