1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
|
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
%}
|
