function W = cubic_kernel(width);

% to be called in REDUCE

N = (width-1)/2;

% returns a matrix descirbing the optimal 2-D cubic kernel
% the width of the filter is width

%first step : compute the 1D optimal linear filter
w = zeros(1,width);

%bunch of constants
a1 = -0.32274;
a3 = 0.08709;
t = atan(0.032932/0.080621);
D = 0.42253;
E = 0.22387;
F = 0.27174;

%computation of the vector x
l = round(N/2-eps);
k = [0:l];
% cf loosy description of x; assume k>=0
x = D*a1.^k+(E*cos(k*t)+F*sin(k*t)).*a3.^k;
% now, get the mirrored image
for i=2:l+1
  X(l+i) = x(i);
  X(l-i+2) = x(i);
end
X(l+1)=x(1);

% now get the definition of w
% do w for k>=0 first and then get the mirrored image
if (mod(N,2)==0)
  w(N+1:2:width) = X(l+1:2*l+1);
  w(width-1) = (3*X(2*l+1)+6*X(2*l)-X(2*l-1))/8;
  for i=width-3:-2:N+1
    w(i) = (-X((i-N-4)/2+l+1)+9*X((i-N-2)/2+l+1)+9*X((i-N)/2+l+1)-X((i-N+2)/2+l+1))/16; 
  end
  
else
  w(N+1:2:width-1) = X(l+1:2*l);
  for i=width-2:-2:N+1
    w(i) = (-X((i-N-4)/2+l+1)+9*X((i-N-2)/2+l+1)+9*X((i-N)/2+l+1)-X((i-N+2)/2+l+1))/16; 
  end
end

% mirror image
for i=1:N
  w(N+1-i) = w(N+1+i);
end


W = w'*w;

%normalisation
W = W/sum(sum(W));
