% ================================================
% *** FUNCTION sprague
% ***
% *** Interpolates an n by w matrix of spectra
% *** where n is the number of spectra
% *** f is an interpolation factor
% *** e.g. if f=2 the sampling rate is doubled
% ===================================================
function [p] = sprague(spectra,f)

if (f<2 | ((f-floor(f))>0))
    disp('invalid f value - premature termination');
    p = 0;
    return;
end
% set the parameters
c = [884 -1960 3033 -2648 1080 -180; ...
508 -540 488 -367 144 -24; ...
-24 144 -367 488 -540 508; ...
-180 1080 -2648 3033 -1960 884];
[numSpectra lengthSpectra] = size(spectra);
for i=1:numSpectra
    % select a spectrum
    r = spectra(i,:);
    % add the extra start and end points
    k = c(1,:);
    p1 = (k*(r(1:6))')/209;
    k = c(2,:);
    p2 = (k*(r(1:6))')/209;
    k = c(3,:);
    p3 = (k*(r(11:16))')/209;
    k = c(4,:);
    p4 = (k*(r(11:16))')/209;
    r = [p1 p2 r p3 p4];
    N = lengthSpectra+4;
    p = zeros(numSpectra,f*(N-5)+1);
    xx = linspace(1/f,1-1/f,f-1);
    for j=3:N-3
        a0 = r(j);
        a1 = (2*r(j-2)-16*r(j-1)+16*r(j+1)-2*r(j+2))/24;
        a2 = (-r(j-2)+16*r(j-1)-30*r(j)+16*r(j+1)-r(j+2))/24;
        a3 = (-9*r(j-2)+39*r(j-1)-70*r(j)+66*r(j+1)-33*r(j+2)+7*r(j+3))/24;
        a4 = (13*r(j-2)-64*r(j-1)+126*r(j)-124*r(j+1)+61*r(j+2)-12*r(j+3))/24;
        a5 = (-5*r(j-2)+25*r(j-1)-50*r(j)+50*r(j+1)-25*r(j+2)+5*r(j+3))/24;
        y = a0+a1*xx+a2*xx.^2+a3*xx.^3+a4*xx.^4+a5*xx.^5;
        index = j-2;
        p(i,(index-1)*f+1) = r(j);
        p(i,(index-1)*f+1+1:(index-1)*f+1+f-1)=y;
    end
    p(i,f*(N-5)+1)=r(N-2);
end

end
% ===================================================
% *** END FUNCTION sprague
% ===================================================

% ===================================================
% *** FUNCTION r2xyz
% ***
% *** function [xyz] = r2xyz(p, startw, endw, obs)
% *** computes XYZ from reflectance p
% *** p is an n by w matrix of n spectra
% *** e.g. set obs to 'd65_64 for D65 and 1964
% *** the startw and endw variables denote first and
% *** last wavelengths (e.g. 400 and 700) for p which
% *** must be 10-nm data in the range 360-780
% ===================================================
function [xyz] = r2xyz(p, startlam, endlam, obs)
if ((endw>780) | (startw<360) | (rem(endw,10)~=0) | (rem(start,10)~=0))
    disp('startw and endw must be divisible by 10')
    disp('wavelength range must be 360-780 or less');
    return;
end
load weights.mat
% weights.mat contains the tables of weights
if strcmp('a_64',obs)
    cie = a_64;
elseif strcmp('a_31', obs)
    cie = a_31;
elseif strcmp('c_64', obs)
    cie = c_64;
elseif strcmp('c_31', obs)
    cie = c_31;
elseif strcmp('d50_64', obs)
    cie = d50_64;
elseif strcmp('d50_31', obs)
    cie = d50_31;
elseif strcmp('d55_64', obs)
    cie = d55_64;
elseif strcmp('d55_31', obs)
    cie = d55_31;
elseif strcmp('d65_64', obs)
    cie = d65_64;
elseif strcmp('d65_31', obs)
    cie = d65_31;
elseif strcmp('d75_64', obs)
    cie = d75_64;
elseif strcmp('d75_31', obs)
    cie = d75_31;
elseif strcmp('f2_64', obs)
    cie = f2_64;
elseif strcmp('f2_31', obs)
    cie = f2_31;
elseif strcmp('f7_64', obs)
    cie = f7_64;
elseif strcmp('f7_31', obs)
    cie = f7_31;
elseif strcmp('f11_64', obs)
    cie = f11_64;
elseif strcmp('f11_31', obs)
    cie = f11_31;
elseif
    disp('unknown option obs');
    disp('use d65_64 for D65 and 1964 observer');
    return;
end
% check dimensions of p
dim = size(p);
N = ((endw-startw)/10 + 1);
if (length(p) ~= N)
    disp('dimensions of p inconsistent');
    return;
end
% deal with possible truncation of reflectance
i = (startw - 360)/10 + 1;
if (i>1)
    cie(i,:) = cie(i,:) + sum(cie(1:i-1,:));
end
j = i + N - 1;
if (j<43)
    cie(j,:) = cie(j,:) + sum(cie(j+1:43,:));
end
cie = cie(i:j,:);
% the main calculation
xyz = (p*cie)*100/sum(cie(:,2));

end
% ===================================================
% *** END FUNCTION r2xyz
% ===================================================

% ===================================================
% *** FUNCTION cieplot
% ***
% *** function = cieplot()
% *** colour representation of chromaticity diagram
% ===================================================
function [] = cieplot()
% load spectral locus xy values at 1-nm intervals
load locus.mat
plot(locus(:,1),locus(:,2),'k','LineWidth',2)
grid on; hold on
axis([0.0 0.75 0.0 0.85])
xlabel('x')
ylabel('y')
% plot the non-spectral locus
plot([locus(1,1) locus(end,1)], [locus(1,2)
locus(end,2)],'k','LineWidth',2)
% chromaticity coordinates of spectrum locus
x = [0.175596 0.172787 0.170806 0.170085 0.160343 ...
0.146958 0.139149 0.133536 0.126688 0.115830 ...
0.109616 0.099146 0.091310 0.078130 0.068717 ...
0.054675 0.040763 0.027497 0.016270 0.008169 ...
0.004876 0.003983 0.003859 0.004646 0.007988 ...
0.013870 0.022244 0.027273 0.032820 0.038851 ...
0.045327 0.052175 0.059323 0.066713 0.074299 ...
0.089937 0.114155 0.138695 0.154714 0.192865 ...
0.229607 0.265760 0.301588 0.337346 0.373083 ...
0.408717 0.444043 0.478755 0.512467 0.544767 ...
0.575132 0.602914 0.627018 0.648215 0.665746 ...
0.680061 0.691487 0.700589 0.707901 0.714015 ...
0.719017 0.723016 0.734674 ]';
y = [0.005295 0.004800 0.005472 0.005976 0.014496 ...
0.026643 0.035211 0.042704 0.053441 0.073601 ...
0.086866 0.112037 0.132737 0.170464 0.200773 ...
0.254155 0.317049 0.387997 0.463035 0.538504 ...
0.587196 0.610526 0.654897 0.675970 0.715407 ...
0.750246 0.779682 0.792153 0.802971 0.812059 ...
0.819430 0.825200 0.829460 0.832306 0.833833 ...
0.833316 0.826231 0.814796 0.805884 0.781648 ...
0.754347 0.724342 0.692326 0.658867 0.624470 ...
0.589626 0.554734 0.520222 0.486611 0.454454 ...
0.424252 0.396516 0.372510 0.351413 0.334028 ...
0.319765 0.308359 0.299317 0.292044 0.285945 ...
0.280951 0.276964 0.265326 ]';
N = length(x);
i = 1;
e = 1/3;
steps = 25;
xy4rgb = zeros(N*steps*4,5,'double');
for w = 1:N % wavelength
    w2 = mod(w,N)+1;
    a1 = atan2(y(w)-e,x(w)-e); % start angle
    a2 = atan2(y(w2)-e,x(w2)-e); % end angle
    r1 = ((x(w)-e)^2 + (y(w)- e)^2)^0.5;% start radius
    r2 = ((x(w2)-e)^2 + (y(w2)-e)^2)^0.5; % end radius
    for c = 1:steps % colourfulness
        % patch polygon
        xyz(1,1) = e+r1*cos(a1)*c/steps;
        xyz(1,2) = e+r1*sin(a1)*c/steps;
        xyz(1,3) = 1 - xyz(1,1) - xyz(1,2);
        xyz(2,1) = e+r1*cos(a1)*(c-1)/steps;
        xyz(2,2) = e+r1*sin(a1)*(c-1)/steps;
        xyz(2,3) = 1 - xyz(2,1) - xyz(2,2);
        xyz(3,1) = e+r2*cos(a2)*(c-1)/steps;
        xyz(3,2) = e+r2*sin(a2)*(c-1)/steps;
        xyz(3,3) = 1 - xyz(3,1) - xyz(3,2);
        xyz(4,1) = e+r2*cos(a2)*c/steps;
        xyz(4,2) = e+r2*sin(a2)*c/steps;
        xyz(4,3) = 1 - xyz(4,1) - xyz(4,2);
        % compute sRGB for vertices
        rgb = xyz2srgb(xyz');
        % store the results
        xy4rgb(i:i+3,1:2) = xyz(:,1:2);
        xy4rgb(i:i+3,3:5) = rgb';
        i = i + 4;
    end
end
[rows cols] = size(xy4rgb);
f = [1 2 3 4];
v = zeros(4,3,'double');
for i = 1:4:rows
    v(:,1:2) = xy4rgb(i:i+3,1:2);
    patch('Vertices',v, 'Faces',f,'EdgeColor','none', ...
    'FaceVertexCData',xy4rgb(i:i+3,3:5), 'FaceColor','interp')
end

end
% ===================================================
% *** END FUNCTION cieplot
% ===================================================

% ===================================================
% *** FUNCTION xyz2srgb
% ***
% *** function [RGB] = xyz2srgb(XYZ)
% *** computes 8-bit sRGB from XYZ
% *** XYZ is n by 3 and in the range 0-1
% *** see also srgb2xyz
function [RGB] = xyz2srgb(XYZ)
if (size(XYZ,2)~=3)
    disp('XYZ must be n by 3'); return;
end
M = [3.2406 -1.5372 -0.4986; -0.9689 1.8758 0.0415; 0.0557 -0.2040 1.0570];
RGB = (M*XYZ')';
RGB(RGB<0) = 0;
RGB(RGB>1) = 1;
DACS = zeros(size(XYZ));
index = (RGB<=0.0031308);
DACS = DACS + (index).*(12.92*RGB);
DACS = DACS + (1-index).*(1.055*RGB.^(1/2.4)-0.055);
RGB=round(DACS*255);
RGB(RGB<0) = 0;

RGB(RGB>255) = 255;

end
% ===================================================
% *** END FUNCTION xyz2srgb
% ===================================================

% ===================================================
% *** FUNCTION srgb2xyz
% ***
% *** function [XYZ] = srgb2xyz(RGB)
% *** computes XYZ from 8-bit RGB
% *** RGB is n by 3 and in the range 0-255
% *** XYZ is returned in the range 0-1
% *** see also xyz2srgb
function [XYZ] = srgb2xyz(RGB)
if (size(RGB,2)~=3)
    disp('RGB must be n by 3');
    return;
end

XYZ = zeros(size(RGB));
M = [0.4124 0.3576 0.1805; 0.2126 0.7152 0.0722; 0.0193 0.1192 0.9505];
DACS=RGB/255;
RGB = zeros(size(RGB));
index = (DACS<=0.04045);
RGB = RGB + (index).*(DACS/12.92);
RGB = RGB + (1-index).*((DACS+0.055)/1.055).^2.4;
XYZ = (M*RGB')';

end

% ===================================================
% *** FUNCTION xyz2lab
% ***
% *** function [lab] = xyz2lab(xyz, obs, xyzw)
% *** computes LAB from XYZ
% *** xyz is an n by 3 matrix
% *** e.g. set obs to 'd65_64 for D65 and 1964
% *** set obs to 'user' to use optional argument
% *** xyzw as the white point
% ===================================================
function [lab] = xyz2lab(xyz,obs,xyzw)

if (size(xyz,2)~=3)
    disp('xyz must be n by 3');
    return;
end
lab = zeros(size(xyz,1),size(xyz,2));

if strcmp('a_64',obs)
    white=[111.144 100.00 35.200];
elseif strcmp('a_31', obs)
    white=[109.850 100.00 35.585];
elseif strcmp('c_64', obs)
    white=[97.285 100.00 116.145];
elseif strcmp('c_31', obs)
    white=[98.074 100.00 118.232];
elseif strcmp('d50_64', obs)
    white=[96.720 100.00 81.427];
elseif strcmp('d50_31', obs)
    white=[96.422 100.00 82.521];
elseif strcmp('d55_64', obs)
    white=[95.799 100.00 90.926];
elseif strcmp('d55_31', obs)
    white=[95.682 100.00 92.149];
elseif strcmp('d65_64', obs)
    white=[94.811 100.00 107.304];
elseif strcmp('d65_31', obs)
    white=[95.047 100.00 108.883];
elseif strcmp('d75_64', obs)
    white=[94.416 100.00 120.641];
elseif strcmp('d75_31', obs)
    white=[94.072 100.00 122.638];
elseif strcmp('f2_64', obs)
    white=[103.279 100.00 69.027];
elseif strcmp('f2_31', obs)
    white=[99.186 100.00 67.393];
elseif strcmp('f7_64', obs)
    white=[95.792 100.00 107.686];
elseif strcmp('f7_31', obs)
    white=[95.041 100.00 108.747];
elseif strcmp('f11_64', obs)
    white=[103.863 100.00 65.607];
elseif strcmp('f11_31', obs)
    white=[100.962 100.00 64.350];
elseif strcmp('user', obs)
    white=xyzw;
else
    disp(''unknown option obs');
    disp(''use d65_64 for D65 and 1964 observer');
    return;
end

lab = zeros(size(xyz,1),3);

fx = zeros(size(xyz,1),3);
for i=1:3
    index = (xyz(:,i)/white(i) > (6/29)^3);
    fx(:,i) = fx(:,i) + index.*(xyz(:,i)/white(i)).^(1/3);
    fx(:,i) = fx(:,i) + (1-index).*((841/108)*(xyz(:,i)/white(i)) + 4/29);
end

lab(:,1)=116*fx(:,2)-16;
lab(:,2) = 500*(fx(:,1)-fx(:,2));
lab(:,3) = 200*(fx(:,2)-fx(:,3));

end
% ====================================================
% *** END FUNCTION xyz2lab
% ====================================================