function I_interp = get_intensity_interp(I, x, mode)

% I_interp = get_intensity_interp(I, x, mode)
%
% DESC:
% Computes the intesity values at the points x using nearest neighbour or
% bilinear interpolation. Note that the coordinates x do not need to be
% interger values.
%
% AUTHOR
% Marco Zuliani - zuliani@ece.ucsb.edu
%
% VERSION:
% 1.1.1
%
% INPUT:
% I             = single channel image
% x             = point coordinates
% mode          = interpolation mode
%                 0 -> nearest neighbour
%                 1 -> bilinear (default)
%                 2 -> cubic (not tested yet, don't use)
%
% OUTPUT:
% I_interp      = interpolated values of intensities

% HISTORY
% 1.0.0         - ??/??/04 - initial version
% 1.1.0         - ??/??/04 - improved clarity of the code and speed
% 1.1.1         - 11/02/07 - minor improvements


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if size(I, 3) > 1
    error('I has to be a single channel image')
end;

N = size(x, 2);
I_interp = zeros(1, N);

if nargin == 2
    mode = 1;
end;

if mode == 2
    offsets = [-1 0 1 2];
    F = zeros(4, 4);
end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Interpolate
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for h = 1:N

    xx = floor(x(:, h));

    switch mode

        % nearest neighbour
        case 0

            I_interp(h) = double(I(xx(1), xx(2)));

        % bilinear
        case 1

            delta = x(:, h) - xx;
            delta_m = 1 - delta;

            F_00 = double(I(xx(1), xx(2)));
            F_01 = double(I(xx(1), xx(2)+1));
            F_10 = double(I(xx(1)+1, xx(2)));
            F_11 = double(I(xx(1)+1, xx(2)+1));

            F_A = F_00*delta_m(1)+F_10*delta(1);
            F_B = F_01*delta_m(1)+F_11*delta(1);

            I_interp(h) = F_A*delta_m(2)+F_B*delta(2);

        case 2

            for j = 1:4
                for i = 1:4
                    F(i, j) = double(I(xx(1)+offsets(i), xx(2)+offsets(j)));
                end;
            end;

            dx = x(1, h) - xx(1);
            dy = x(2, h) - xx(2);

            F_0 = F(1);
            F_1 = F(2);
            F_2 = F(3);
            F_3 = F(4);
            F_4 = F(5);
            F_5 = F(6);
            F_6 = F(7);
            F_7 = F(8);
            F_8 = F(9);
            F_9 = F(10);
            F_10 = F(11);
            F_11 = F(12);
            F_12 = F(13);
            F_13 = F(14);
            F_14 = F(15);
            F_15 = F(16);
            
            % this has been automatically generated...
            % looks ugly and should e manually simplified
            I_interp(h) = ...
                1/6*F_0*(-2/3*dx^3+(dx+1)^3-2/3*(dx+2)^3+1/6*(dx+3)^3)*(1-dy)^3+...
                1/6*F_1*(dx^3-2/3*(dx+1)^3+1/6*(dx+2)^3)*(1-dy)^3+...
                1/6*F_2*(-2/3*dx^3+1/6*(dx+1)^3)*(1-dy)^3+...
                1/36*F_3*dx^3*(1-dy)^3+...
                F_4*(-2/3*dx^3+(dx+1)^3-2/3*(dx+2)^3+1/6*(dx+3)^3)*(-2/3*(1-dy)^3+1/6*(-dy+2)^3)+...
                F_5*(dx^3-2/3*(dx+1)^3+1/6*(dx+2)^3)*(-2/3*(1-dy)^3+1/6*(-dy+2)^3)+...
                F_6*(-2/3*dx^3+1/6*(dx+1)^3)*(-2/3*(1-dy)^3+1/6*(-dy+2)^3)+...
                1/6*F_7*dx^3*(-2/3*(1-dy)^3+1/6*(-dy+2)^3)+...
                F_8*(-2/3*dx^3+(dx+1)^3-2/3*(dx+2)^3+1/6*(dx+3)^3)*((1-dy)^3-2/3*(-dy+2)^3+1/6*(3-dy)^3)+...
                F_9*(dx^3-2/3*(dx+1)^3+1/6*(dx+2)^3)*((1-dy)^3-2/3*(-dy+2)^3+1/6*(3-dy)^3)+...
                F_10*(-2/3*dx^3+1/6*(dx+1)^3)*((1-dy)^3-2/3*(-dy+2)^3+1/6*(3-dy)^3)+...
                1/6*F_11*dx^3*((1-dy)^3-2/3*(-dy+2)^3+1/6*(3-dy)^3)+...
                F_12*(-2/3*dx^3+(dx+1)^3-2/3*(dx+2)^3+1/6*(dx+3)^3)*(-2/3*(1-dy)^3+(-dy+2)^3-2/3*(3-dy)^3+1/6*(-dy+4)^3)+...
                F_13*(dx^3-2/3*(dx+1)^3+1/6*(dx+2)^3)*(-2/3*(1-dy)^3+(-dy+2)^3-2/3*(3-dy)^3+1/6*(-dy+4)^3)+...
                F_14*(-2/3*dx^3+1/6*(dx+1)^3)*(-2/3*(1-dy)^3+(-dy+2)^3-2/3*(3-dy)^3+1/6*(-dy+4)^3)+...
                1/6*F_15*dx^3*(-2/3*(1-dy)^3+(-dy+2)^3-2/3*(3-dy)^3+1/6*(-dy+4)^3);

    end;

end;

return