function val = bicubicInterpolationMulti(x, y, data)
%B. Vandeportaele
%https://en.wikipedia.org/wiki/Bicubic_interpolation
%https://fr.wikipedia.org/wiki/Interpolation_bicubique
%A is a k.l.m matrix, m being the number of channels
nbchannels=size(data,3);
%preallocate the return value
val=zeros(nbchannels,1);
%Control points
x0 = floor(x);
y0 = floor(y);
x_1 = x0-1;
x1 = x0+1;
x2 = x0+2;
y_1 = y0-1;
y1 = y0+1;
y2 = y0+2;

for n=1:nbchannels
    % Function evaluated at control points
    f_1_1 = data(y_1,x_1,n);
    f_10 = data(y0,x_1,n);
    f_11 = data(y1,x_1,n);
    f_12 = data(y2,x_1,n);
    f0_1 = data(y_1,x0,n);
    f00 = data(y0,x0,n);
    f01 = data(y1, x0,n);
    f02 = data(y2, x0,n);
    f1_1 = data(y_1, x1,n);
    f10 = data(y0, x1,n);
    f11 = data(y1, x1,n);
    f12 = data(y2, x1,n);
    f2_1 = data(y_1, x2,n);
    f20 = data(y0, x2,n);
    f21 = data(y1, x2,n);
    f22 = data(y2, x2,n);
    
    % Derivatives in x evaluated at control points
    fx0_1 = (f1_1 - f_1_1)/2;
    fx00 = (f10 - f_10)/2;
    fx01 = (f11 - f_11)/2;
    fx02 = (f12 - f_12)/2;
    fx1_1 = (f2_1 - f0_1)/2;
    fx10 = (f20 - f00)/2;
    fx11 = (f21 - f01)/2;
    fx12 = (f22 - f02)/2;
    
    % Derivatives in y evaluated at control points
    fy00 = (f01 - f0_1)/2;
    fy10 = (f11 - f1_1)/2;
    fy01 = (f02 - f00)/2;
    fy11 = (f12 - f10)/2;
    
    % Derivatives in x and y evaluated at control points
    fxy00 = (fx01 - fx0_1)/2;
    fxy10 = (fx11 - fx1_1)/2;
    fxy01 = (fx02 - fx00)/2;
    fxy11 = (fx12 - fx10)/2;
    
    % Interpolation
    Px=[1,(x-x0),(x-x0)^2,(x-x0)^3];
    %don't know why but sometime the generated Px is a column vector
    Px=reshape(Px,1,4);
    Py=[1;(y-y0);(y-y0)^2;(y-y0)^3];
    
    A=[ 1  0  0  0; ... Coefficients aij
        0  0  1  0; ...
        -3  3 -2 -1; ...
        2 -2  1  1];
    F=[ f00  f01  fy00  fy01; ...
        f10  f11  fy10  fy11; ...
        fx00 fx01 fxy00 fxy01; ...
        fx10 fx11 fxy10 fxy11];
    val(n)=Px*A*F*A'*Py;
end