% rgb2lms.m
%
% Takes a vector of sRGB numbers, calculates CIE XYZ, and converts those to
% approximate LMS coordinates.
%
% sRGB is a standard RGB space used by Microsoft OS's, Hewlett-Packard imaging devices,
% and on the World-Wide Web. All images produced by HP scanners and digital cameras are
% encoded in the sRGB space. For more information, see:
%
% 	http://www.srgb.com/
%
% * Conversion to CIE(1931) XYZ is done using the sRGB standard (nonlinearity as
%   described on the above web pages, sRGB XYZ coordinates from an ICM profile
%   downloaded there).
% * Conversion to LMS cone excitation coordinates is done using the method described
%   in Y. Nakano's appendix in COLOR VISION by Kaiser & Boynton.
% 
% 1/6/2001 rdb   * Wrote it: convertRGB.m.
% 12/14/2001 rdb * Changed yourXYZ calculation to simple weighting of phosphor XYZ coordinates
%                * Corrected error in Yrb calculation b=S/(L+M), not M/(L+M)
%		 * Turned it into function rgb2lms
%		 * Changed to use sRGB as input space 
%		 * NOTE: not tested yet!	
%
function [LMS] = rgb2lms(sRGB) 

% normalize:
sRGBnorm = sRGB./255;

% remove sRGB gamma nonlinearity:
for n=1:3
   if sRGBnorm(n) <= 0.04045
      RGB(n) = sRGBnorm(n)./12.92;
   else
      RGB(n) = ((sRGBnorm(n)+0.055)./1.055).^2.4;
   end
end

% sRGB primaries (ITU-R BT.709 reference primaries) 
rXYZ = [0.436066 0.222488 0.013916];
gXYZ = [0.385147 0.716873 0.097076];
bXYZ = [0.143066 0.060608 0.714096];

% phosphor chromaticity coordinates
xr = rXYZ(1) ./ (rXYZ(1) + rXYZ(2) + rXYZ(3));
yr = rXYZ(2) ./ (rXYZ(1) + rXYZ(2) + rXYZ(3));
zr = rXYZ(3) ./ (rXYZ(1) + rXYZ(2) + rXYZ(3));
xg = gXYZ(1) ./ (gXYZ(1) + gXYZ(2) + gXYZ(3));
yg = gXYZ(2) ./ (gXYZ(1) + gXYZ(2) + gXYZ(3));
zg = gXYZ(3) ./ (gXYZ(1) + gXYZ(2) + gXYZ(3));
xb = bXYZ(1) ./ (bXYZ(1) + bXYZ(2) + bXYZ(3));
yb = bXYZ(2) ./ (bXYZ(1) + bXYZ(2) + bXYZ(3));
zb = bXYZ(3) ./ (bXYZ(1) + bXYZ(2) + bXYZ(3));

% from Kaiser & Boynton:
RGBXYZ = [[xr./yr xg./yg xb./yb]       % phosphor coordinate (RGB) to XYZ conversion matrix
          [1      1      1     ]
          [zr./yr zg./yg zb./yb]];
k = [[0.2535 0       0      ]          % arbitrary scaling
     [0      -0.4000 0      ]
     [0      0       0.01327]];
SP75 = [[0.7465 1.4000  0.1748]        % Smith and Pokorny (1975) copunctual points for dichromats
        [0.2535 -0.4000 0.0000]
        [0.0000 0.0000  0.8252]];
XYZLMS = k * inv(SP75);                % XYZ to LMS conversion matrix


XYZ = RGB(1).*rXYZ + RGB(2).*gXYZ + RGB(3).*bXYZ;
Yxy = [XYZ(2) XYZ(1)./(XYZ(1)+XYZ(2)+XYZ(3)) XYZ(2)./(XYZ(1)+XYZ(2)+XYZ(3))];
LMS = XYZLMS * XYZ';
Yrb = [LMS(1)+LMS(2) LMS(1)./(LMS(1)+LMS(2)) LMS(3)./(LMS(1)+LMS(2))];