function [M,LMLumWeights] = ComputeDKL_M(bg,T_cones,T_Y)
% [M,LMLumWeights] = ComputeDKL_M(bg,T_cones,T_Y)
% 
% Compute the matrix that converts between incremental cone
% coordinates and DKL space.  The order of
% the coordinates in the DKL column vectors is (Lum, RG, S)
%
% The code follows that published by Brainard
% as an appendix to Human Color Vision by Kaiser
% and Boynton, but has been generalized to work
% for passed cone fundamentals and luminosity 
% function.
%
% These should be passed in standard Psychtoolbox
% form (LMS in rows of T_cones; luminosity function
% as single row vector T_Y).
%
% The argument bg is the LMS cone coordinates of the
% background that defines the space.
%
% See DKLDemo for proper use of this function.  Also
% DKLToConeInc and ConeIncToDKL.
%
% 8/30/96   dhb  Pulled it out.
% 4/9/05    dhb  Allow passing of cones and luminance to be used.
% 11/17/05  dhb  Require passing of cones and luminance.
%           dhb  Fixed definition of M_raw to handle arbitrary L,M scaling.
% 10/5/12   dhb  Comment specifying coordinate system convention.  Supress extraneous printout.
% 04/13/17  dhb  Return weights that give luminance from sum of L and M cone excitations.
% 08/20/21  dhb  Added check of a direct method of computing M, provided by Supi Ray (sray@iisc.ac.in).
%                The check is commented out at the end of this routine, but
%                does indeed give the same answer in the limited cases I've
%                checked. At the heart of the derivation is analytic
%                inversion of M_raw.  I imagine that Ray would be happy to
%                provide the derivation if asked.

% If cones and luminance are passed, find how L and
% M cone incrments sum to best approximate change in
% luminance.
if (nargin == 3)
	T_LM = T_cones(1:2,:);
	LMLumWeights = T_LM'\T_Y';
else
    fprintf('ComputeDKL_M now requires explicit specification\n');
    fprintf('of cone fundamentals and luminosity function\n');
    fprintf('See DKLDemo\n');
    error('');   
end

% Set M_raw as in equation A.4.9.
% This is found by plugging the background
% values into equation A.4.8.  Different
% backgrounds produce different matrices.
% The Matlab notation below just 
% fills the desired 3-by-3 matrix.
%
% Note that A.4.8 in the Brainard chapter contains
% a typo: the row 1 col 3 entry of the matrix should
% be 0, not 1.  Also, at the top of page 571, there is
% an erroneous negative sign in front of the term
% for W_S-Lum,S.
%
% Finally, A.4.8 as given in the chatper assumes
% that Lum = L + M.  The formula below generalizes
% to arbitrary scaling.
M_raw = [ LMLumWeights(1) LMLumWeights(2) 0 ; ...
			1 -bg(1)/bg(2) 0 ; ...
			-LMLumWeights(1) -LMLumWeights(2) (LMLumWeights(1)*bg(1)+LMLumWeights(2)*bg(2))/bg(3) ];

% Compute the inverse of M for
% equation A.4.10.  The Matlab inv() function
% computes the matrix inverse of its argument.
M_raw_inv = inv(M_raw);

% Find the three isolating stimuli as
% the columns of M_inv_raw.  The Matlab
% notation X(:,i) extracts the i-th column
% of the matrix X.
isochrom_raw = M_raw_inv(:,1);
rgisolum_raw = M_raw_inv(:,2);
sisolum_raw = M_raw_inv(:,3);

% Find the pooled cone contrast of each
% of these.  The Matlab norm() function returns
% the vector length of its argument.  The Matlab
% ./ operation represents entry-by-entry division.
isochrom_raw_pooled = norm(isochrom_raw ./ bg);
rgisolum_raw_pooled = norm(rgisolum_raw ./ bg);
sisolum_raw_pooled = norm(sisolum_raw ./ bg);

% Scale each mechanism isolating
% modulation by its pooled contrast to obtain
% mechanism isolating modulations that have
% unit length.
isochrom_unit = isochrom_raw / isochrom_raw_pooled;
rgisolum_unit = rgisolum_raw / rgisolum_raw_pooled;
sisolum_unit = sisolum_raw / sisolum_raw_pooled;

% Compute the values of the normalizing
% constants by plugging the unit isolating stimuli
% into A.4.9 and seeing what we get.  Each vector
% should have only one non-zero entry.  The size
% of the entry is the response of the unscaled
% mechanism to the stimulus that should give unit
% response.
lum_resp_raw = M_raw*isochrom_unit;
l_minus_m_resp_raw = M_raw*rgisolum_unit;
s_minus_lum_resp_raw = M_raw*sisolum_unit;
					 
% We need to rescale the rows of M_raw
% so that we get unit response.  This means
% multiplying each row of M_raw by a constant.
% The easiest way to accomplish the multiplication
% is to form a diagonal matrix with the desired
% scalars on the diagonal.  These scalars are just
% the multiplicative inverses of the non-zero
% entries of the vectors obtained in the previous
% step.  The resulting matrix M provides the
% entries of A.4.11.  The three _resp vectors
% computed should be the three unit vectors
% (and they are).
D_rescale = [1/lum_resp_raw(1) 0 0 ; ...
						 0 1/l_minus_m_resp_raw(2) 0 ; ...
						 0 0 1/s_minus_lum_resp_raw(3) ];				 
M = D_rescale*M_raw;
lum_resp = M*isochrom_unit;
l_minus_m_resp = M*rgisolum_unit;
s_minus_lum_resp = M*sisolum_unit;

% Compute the inverse of M to obtain
% the matrix in equation A.4.12.
M_inv = inv(M);

% Alternate method of computing M that does not
% require inversion of M_raw. This method was developed by 
% Supi Ray.  It gives the same answer (M_alt) for M as the
% code above, to numerical precision.  It would be possible
% to do out the matrix multiplications below to obtain 
% direct expression for each entry of M_alt, but I have 
% not done so. Examinging the final expression might provide
% additional intution about the factors influencing the normalized
% matrix M.
%{
W = diag([LMLumWeights(1) LMLumWeights(2) 1]);
bgNorm = W*bg;
B1 = [1 1 0 ; 1 -bgNorm(1)/bgNorm(2) 0 ; -1 -1 (bgNorm(1) + bgNorm(2))/bgNorm(3)];
B2 = B1*W;
D_rescale_alt = diag([sqrt(3) LMLumWeights(1)*sqrt( 1+(bgNorm(2)/bgNorm(1))^2 ) 1])/(bgNorm(1)+bgNorm(2));
M_raw_alt = [ LMLumWeights(1) LMLumWeights(2) 0 ; ...
			1 -(LMLumWeights(1)*bg(1))/(LMLumWeights(2)*bg(2)) 0 ; ...
			-LMLumWeights(1) -LMLumWeights(2) (LMLumWeights(1)*bg(1)+LMLumWeights(2)*bg(2))/bg(3) ];
M_alt = D_rescale_alt*M_raw;
if ( any(abs(M(:)-M_alt(:)) > 1e-10) )
    error('Two ways of computing M do not agree');
end
%}