function [out] = csi1D (in1, in2)
% [out] = csi1D (in1, in2)
%
% 1-D Cubic-Spline Interpolation
% Written by Matt Golub for Psych 221, Winter Quarter, 2008.
%
% INPUTS:  in1 -  A 1D set of CFA output data to be interpolated
%          in2  - A 1D set of CFA output data to be used for hue 
%                 information when interpolating in1.
%Both inputs must be row vectors of the same length.  Length must be
%greater than or equal to 7.
%
% OUTPUTS: out - The result of cubic-spline interpolation over in1.

% This matrix is used to solve the 12 simultaneous equations
% needed to construct the cubic-splines.
M = [0 0 0 1 0 0 0 0 0 0 0 0;
     8 4 2 1 0 0 0 0 0 0 0 0;
     0 0 0 0 8 4 2 1 0 0 0 0;
     0 0 0 0 64 16 4 1 0 0 0 0;
     0 0 0 0 0 0 0 0 64 16 4 1;
     0 0 0 0 0 0 0 0 216 36 6 1;
     12 4 1 0 -12 -4 -1 0 0 0 0 0;
     0 0 0 0 48 8 1 0 -48 -8 -1 0;
     12 2 0 0 -12 -2 0 0 0 0 0 0; 
     0 0 0 0 24 2 0 0 -24 -2 0 0;
     -1 -1 -1 -1 27 9 3 1 0 0 0 0;
     0 0 0 0 -27 -9 -3 -1 125 25 5 1];

Minv=inv(M);

N = numel(in1);
if N~=numel(in2) || N<7
    error('Invalid input vector(s)');
end

%Due to the shift-variant nature of this algorithm, we must
%ensure that green data is available at the bein1ning of each
%1x7 interpolation vector. We will determine here if we must offset
%the inputs to ensure green data in the first position.
%(We want in1 to look like [G 0 G 0...], not [0 G 0 G...]
if (in1(1)+in1(3)+in1(5)+in1(7)==0)
       offset = true;
       in1=[in1(2) in1];
       in2 = [0 in2];
else offset = false;
end

%Each interpolation is over 7 contiguous pixels
%Both inputs must be padded with nearest neighbor value
%so that each input has an odd length of at least 7.
if (N<7)
    padLen = 7-N;
else padLen = -(mod(N,2)-1);
end
if in1(end)==0
    pad1 = in1(end-1);
    pad2 = 0;
else
    pad1 = 0;
    pad2 = in2(end-1);
end
in1 = [in1 pad1*ones(1,padLen)];
in2 = [in2 pad2*ones(1,padLen)];

newLen = numel(in1);
numInterps = 1+floor((newLen-7)/2);



out=in1;

for(n=1:numInterps)
    start=2*(n-1);
    
    V = [in1(start+1) in1(start+3) in1(start+3) in1(start+5) in1(start+5) in1(start+7) 0 0 0 0 in2(start+4)-in2(start+2) in2(start+6)-in2(start+4)]';
    C = Minv*V;
    
    %Full equations (aren't computed to improve efficiency)
    %Only compute where we will use interpolated value
    %x=0:6;
    %s0 = C(1).*x.^3+C(2).*x.^2+C(3).*x+C(4);
    %s1 = C(5).*x.^3+C(6).*x.^2+C(7).*x+C(8);
    %s2 = C(9).*x.^3+C(10).*x.^2+C(11).*x+C(12);

    if (n==1)  %Border case
        out(start+2) = C(1)*1^3+C(2)*1^2+C(3)*1+C(4); % = s0 evaluated at x=1
    elseif (n==numInterps) %Border case
        out(start+6) = C(9)*5^3+C(10)*5^2+C(11)*5+C(12); % = s2 evaluated at x=5
    end
    
    out(start+4) = C(5)*3^3+C(6)*3^2+C(7)*3+C(8); % = s1 evaluated at x=3
    
  %  if (var(in1(start+1:2:start+7))>.2)
  %  debugPlot;
  %  end;
    
end

% Correct for the offsets made earlier
if (offset)
    out = [out(2:end-(padLen))];
else out=out(1:(end-padLen));
end

% Limit saturated values
negatives = find(out<0);
out(negatives) = 0;
saturates = find(out>1);
out(saturates) = 1;
    
end