
function [QMh] = ciecam_qm(xyz, xyzw, LA, c, Nc, FLL, F, Yb, D);

%
%  CIECAM_QM	function [QMh] = ciecam_qm(xyz, xyzw, LA, 
%					   c, Nc, FLL, F, 
%					   Yb, D)
%
%	calculates CIECAM97s QMh coordinates of an XYZ color
%
%  inputs:
%	
%	xyz	color in XYZ colorspace
%	xyzw	referencey white in XYZ colorspace
%	LA	absolute adapting luminance
%	c	impact of surround		(see table)
%	Nc	chromatic induction factor	(see table)
%	FLL	lightness contrast factor	(see table)
%	F	degree of adaptation		(see table)
%	Yb	background relative luminance (in the 
%		  same units of xyzw)
%	D 	=1 to discount the illuminant
%
%  outputs:
%
%	QMh	QMh color vector
%
%
%  ref: Condition               c       Nc      FLL     F
%       Avg., Stim.>4deg.       0.69    1.0     0.0     1.0
%       Avg.                    0.69    1.0     1.0     1.0
%       Dim                     0.59    1.1     1.0     0.9
%       Dark                    0.525   0.8     1.0     0.9
%       Cut Sheet               0.41    0.8     1.0     0.9
%
%  (adapted from IDL code by Mark Fairchild)
%  Krishna S. Nayak 3/13/98
%  Leland Stanford Junior University
%


MB =    [       0.8951,         0.2664,        -0.1614 ;  ...
               -0.7502,         1.7135,         0.0367 ;  ...
                0.0389,        -0.0685,         1.0296 ];

MH =    [       0.38971,        0.68898,       -0.07868 ; ...
               -0.22981,        1.18340,        0.04641 ; ...
                0.00,           0.00,           1.00    ];

rgb = MB * (xyz/xyz(2));
rgbw = MB * (xyzw/xyzw(2));

if (D == 1.0) 
  D = F - (F/(1.+2.*LA^0.25+(LA^2.)/300.));
end
p = rgbw(3)^0.0834;

rgbc(1) = (D*(1.0/rgbw(1))+1.0-D)*rgb(1);
rgbc(2) = (D*(1.0/rgbw(2))+1.0-D)*rgb(2);
rgbc(3) = (rgb(3)/abs(rgb(3)))*((D*(1.0/(rgbw(3)^p))+1.0-D)*abs(rgb(3))^p);

rgbcw(1) = (D*(1.0/rgbw(1))+1.0-D)*rgbw(1);
rgbcw(2) = (D*(1.0/rgbw(2))+1.0-D)*rgbw(2);
rgbcw(3) = (rgbw(3)/abs(rgbw(3)))*((D*(1.0/(rgbw(3)^p))+1.0-D)*abs(rgbw(3))^p);

MBinv = inv(MB);

rgbprime = MH * (MBinv * (rgbc'*xyz(2)));

rgbprimew = MH * (MBinv * (rgbcw'*xyzw(2)));

k = 1.0/(5.0*LA +1);
FL = 0.2*(k^4.0)*5.0*LA+0.1*((1-k^4.0)^2.0)*(5.0*LA)^(1./3.);
n = Yb/xyzw(2);
Nbb = 0.725*(1.0/n)^0.2;
Ncb = 0.725*(1.0/n)^0.2;
z = 1.0+FLL*n^0.5;

rgba = (40.0*(FL*rgbprime/100.).^0.73) ./ 	...
       ((FL*rgbprime/100).^0.73+2.0)     + 1.0;

rgbaw = (40.0*(FL*rgbprimew/100.).^0.73) ./	...
       	((FL*rgbprimew/100).^0.73+2.0)   + 1.0;

a = rgba(1)-12.0*rgba(2)/11.0+rgba(3)/11.0;
b = (1.0/9.0)*(rgba(1)+rgba(2)-2.0*rgba(3));

h = (180.0/pi)*atan2(b,a);

if (h < 0.0) 
  h=360.0+h;
end

JCh(3) = h;
QMh(3) = h;

if (h <= 20.14) 
  e1 = 0.8565; e2 = 0.8; h1 = 0.0; h2 = 20.14;
elseif (h <= 90.0) 
  e1 = 0.8; e2 = 0.7; h1 = 20.14; h2 = 90.0;
elseif (h <= 164.25) 
  e1 = 0.7; e2 = 1.0; h1 = 90.0; h2 = 164.25;
elseif (h <= 237.53) 
  e1 = 1.0; e2 = 1.2; h1 = 164.25; h2 = 237.53;
else
  e1 = 1.2; e2 = 0.8565; h1 = 237.53; h2 = 360.0;
end

e = e1 +(e2-e1)*(h-h1)/(h2-h1);

Astim = (2*rgba(1)+rgba(2)+(1.0/20.0)*rgba(3)-2.05)*Nbb;
Awhite = (2*rgbaw(1)+rgbaw(2)+(1.0/20.0)*rgbaw(3)-2.05)*Nbb;

JCh(1) = 100.0*(Astim/Awhite)^(c*z);
QMh(1) = (1.24/c)*((JCh(1)/100.0)^0.67)*(Awhite+3.0)^0.9;

s = (50.0*sqrt(a^2.0+b^2.0)*100.0*e*(10.0/13.0)*Nc*Ncb) / ...
      (rgba(1)+rgba(2)+(21.0/20.0)*rgba(3));

JCh(2) = 2.44*(s^0.69)*((JCh(1)/100.0)^(0.67*n))*(1.64-0.29^n);
QMh(2) = JCh(2)*FL^0.15;

