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
|
function ls = LMSToMacBoyn(LMS,T_cones,T_lum)
% ls = LMSToMacBoyn(LMS,[T_cones,T_lum])
%
% Compute MacLeod-Boynton chromaticity from cone exciation coordinates.
% This is L/Lum and S/Lum, with appropriate normalization as described
% below. Recommended usage yields MB chromacities according to CIE 170-2:2015.
%
% ** Recommended Usage: Compute LMS with respect to some specified T_cones,
% and pass both T_cones as well as the cooresponding T_lum (the photopic
% luminosity function.) T_lum should be a linear combination of the L and M
% cone fundamentals.
%
% This routine will then scale passed L and M values so that they sum to
% the best linear approximation of luminance and then normalize the L
% excitation by luminance (as computed by the linear combination) to obtain
% the l chromaticity. It will normalize the S excitaiton by luminance to
% obtain s chromaticity, with an overall scaling so that the maximum value
% of this chromaticity is 1 taken over the visible spectrum.
%
% Note that the s cone scaling can vary a bit depending on the wavelength
% sampling of the passed T_cones and T_lum, since the max is taken over
% these. If you use the T_cones_ss2/T_cones_ss10 and T_CIE_Y2/T_CIE_Y10
% files provided in PTB, the default sampling is at 1 nm and this is fine.
% If you use subsampled wavelength spacing, the computation of the s cone
% scaling will begin to deviate from the standard. But so will your
% computation of LMS values, so this isn't really an issue specific to this
% routine.
%
% When you use the CIE cone fundamentals and corresponding luminance
% functions, this procedure yields the MacLeod-Boynton chromaticity
% diagrams as specified in CIE 170-2:2015.
%
% ** Legacy Usage: Just pass LMS values. In this case, we assume that the
% passed LMS values were computed with respect to the Smith-Pokorny
% fundamentals normalized to a peak of 1 and Judd-Vos luminance (more or
% less). That is, this usage assumes LMS was computed using the
% fundamentals stored in PTB's T_cones_sp. This is old usage and preserved
% for backwards compatibility, but the three argument usage as described
% above is preferred for clarity. Moreover, in this case, the s
% chromaticity is not further normalized. This leads to S chromaticities
% considerably larger than those obtained with the new usage.
%
% ** A Backwards Incompatibility. The scaling for s chromaticity to match
% CIE 170-2:2015 was introduced in Janurary 2019 and is not backwards
% compatible with previous behavior when T_cones and T_lum are passed.
% Preserving such compatibility did not seem important enough relative to
% the gains of having this work as now specified in the CIE standard.
%
% 10/30/97 dhb Wrote it.
% 7/9/02 dhb T_cones_sp -> T_cones on line 20. Thanks to Eiji Kimura.
% 1/23/19 dhb Scale s chromaticity value to be consistent with CIE
% 170-2:2015, when T_cones and T_lum are passed. This is
% not backwards compatible with previous scaling, but it
% seems good to match the standard. Thanks to Danny Garside
% for pointing out the scaling specified in the 2015
% standard.
% Examples:
%{
% Recreate the spectrum locus and equal energy white shown in Figure 8.2
% of CIE 170-2:2015. Also performs a regression check.
load T_cones_ss2
load T_CIE_Y2
lsSpectrumLocus = LMSToMacBoyn(T_cones_ss2,T_cones_ss2,T_CIE_Y2);
% Compute the sum of the ls values in the spectrum locus, and compare
% to the value that this example computed in February 2019, entered
% here to four places as 412.2608. This comparison provides a
% check that this routine still works the way it did when we put in the
% check.
check = round(sum(lsSpectrumLocus(:)),4);
if (abs(check-412.2608) > 1e-4)
error('No longer get same check value as we used to');
end
% Compute ls for equal energy white
LMSEEWhite = sum(T_cones_ss2,2);
lsEEWhite = LMSToMacBoyn(LMSEEWhite,T_cones_ss2,T_CIE_Y2);
% Plot
figure; hold on;
plot(lsSpectrumLocus(1,:)',lsSpectrumLocus(2,:)','r','LineWidth',3);
plot(lsEEWhite(1),lsEEWhite(2),'bs','MarkerFaceColor','b','MarkerSize',12);
xlim([0.4 1]); ylim([0,1]);
xlabel('l chromaticity'); ylabel('s chromaticity');
title('CIE 170-2:2015 Figure 8.2');
%}
%{
% Demonstrate invariance of ls after scaling of cones and luminance, as
% long as LMS valued are computed with respect to passed cones and
% luminance.
load T_cones_ss2
load T_CIE_Y2
% Spectrum locus ls chromaticity with no scaling.
lsSpectrumLocus = LMSToMacBoyn(T_cones_ss2,T_cones_ss2,T_CIE_Y2);
% Do some scaling and recompute spectrum locus.
% The choices of 1.8, 3, 0.05 as scaling are
% just three numbers I made up. You can use any three numbers and it
% should still work, modulo any numerical issues with very big or
% very small numbers.
T_CIE_Y2_scaled = 1.8*T_CIE_Y2;
T_cones_ss2_scaled = T_cones_ss2;
T_cones_ss2_scaled(1,:) = 3*T_cones_ss2(1,:);
T_cones_ss2_scaled(3,:) = 0.05*T_cones_ss2(3,:);
lsSpectrumLocusScaled = LMSToMacBoyn(T_cones_ss2_scaled,T_cones_ss2_scaled,T_CIE_Y2_scaled);
% Make sure the difference is very small numerically.
diff = max(abs(lsSpectrumLocus(:)-lsSpectrumLocusScaled(:)));
fprintf('Difference in spectrum locus chromaticity after scaling is %0.5f (should be small)\n',diff);
% Plot the locus computed both ways to compare visually.
figure; hold on;
plot(lsSpectrumLocus(1,:)',lsSpectrumLocus(2,:)','r','LineWidth',3);
plot(lsSpectrumLocusScaled(1,:)',lsSpectrumLocusScaled(2,:)','g','LineWidth',2);
xlim([0.4 1]); ylim([0,1]);
xlabel('l chromaticity'); ylabel('s chromaticity');
%}
% Scale LMS so that L+M = luminance and S cone value corresponds to a
% fundamental with a max of 1.
if (nargin == 1)
LMS = diag([0.6373 0.3924 1]')*LMS;
elseif (nargin == 3 )
factorsLM = (T_cones(1:2,:)'\T_lum');
factorS = 1/max(T_cones(3,:)./(factorsLM(1)*T_cones(1,:) + factorsLM(2)*T_cones(2,:)));
LMS = diag([factorsLM ; factorS])*LMS;
else
error('Number of input arguments should be either 1 or 3');
end
% Compute ls coordinates from LMS
n = size(LMS,2);
ls = zeros(2,n);
denom = [1 1 0]*LMS;
ls = LMS([1 3],:) ./ ([1 1]'*denom);
|
