% lms2rgb.m
%
% Takes a vector of LMS cone excitations, calculates CIE XYZ, and converts those to
% approximate sRGB (0) or Sony Trinitron (1) 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.
% 
% 12/14/2001 rdb * from rgb2lms
%		 * NOTE: not complete, not tested yet!	
%
function [RGB] = lms2rgb(LMS, sRGBorTrinitron) 

maxLUM = 1.0;      % maximum L+M

% first calculate conversion matrix:

if sRGBorTrinitron == 0
% 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];
end

if sRGBorTrinitron == 1
% Sony Trinitron phosphor XYZ coordinates
rXYZ = [0.390503 0.209930 0.012619];
gXYZ = [0.31950 0.676605 0.078995];
bXYZ = [0.254196 0.113449 0.733276];
end

% primary 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 = inv(XYZLMS)*LMS';

RGBXYZmat = [[ rXYZ(1) gXYZ(1) bXYZ(1)]
             [ rXYZ(2) gXYZ(2) bXYZ(2)] 
             [ rXYZ(3) gXYZ(3) bXYZ(3)]];

RGBlin = RGBXYZmat*XYZ;

% add nonlinear (gamma) encoding:

   % for standard Sony Trinitrons:
if sRGBorTrinitron == 0
   RGBnorm = RGBlin.^(1./2.5);
end

   % for sRGB:
if sRGBorTrinitron == 1
   for n=1:3
      if RGBlin(n) <= 0.0031308
         RGBnorm(n) = RGBlin(n).*12.92;
      else
         RGBnorm(n) = 1.055.*RGBlin(n).^(1./2.4)-0.055;
      end
   end
end

RGB = round(255.*RGBnorm);