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);