function [k1, k2, k3, k4, k5, k6, k7, k8, k9, k10] = bicubic(a)
%
% bicubic: generate bicubic fit polynomial to image data
%          uses mask convolution
% On input:
%     a: m by n array
% On output:
%     k1..k10: bicubic coefficients (see Haralick and Shapiro, Computer Vision)
% Call:
%       [k1 k2 k3 k4 k5 k6 k7 k8 k9 k10] = bicubic(a)
%

[m,n] = size(a);
a1 = zeros(m+4,n+4);
a1(3:m+2,3:n+2) = a;

k1m = (1/175)* [-13  2  7  2 -13; ...
                  2 17 22 17   2; ...
	 			  7 22 27 22   7; ...
				  2 17 22 17   2; ...
			    -13  2  7  2 -13];

k2m = (1/420)* [ 31  -5 -17  -5  31; ...
                -44 -62 -68 -62 -44; ...
				  0   0   0   0   0; ...
				 44  62  68  62  44; ...
			    -31   5  17   5 -31];

k3m = (1/420)* [ 31 -44 0 44 -31; ...
                 -5 -62 0 62   5; ...
				-17 -68 0 68  17; ...
				 -5 -62 0 62   5; ...
			     31 -44 0 44 -31];

k4m = (1/70)* [ 2  2  2  2  2; ...
               -1 -1 -1 -1 -1; ...
			   -2 -2 -2 -2 -2; ...
			   -1 -1 -1 -1 -1; ...
			    2  2  2  2  2];

k5m = (1/100)* [ 4  2 0 -2 -4; ...
                 2  1 0 -1 -2; ...
				 0  0 0  0  0; ...
				-2 -1 0  1  2; ...
			    -4 -2 0  2  4];

k6m = (1/420)* [2 -1 -2 -1 2; ...
                2 -1 -2 -1 2; ...
				2 -1 -2 -1 2; ...
				2 -1 -2 -1 2; ...
			    2 -1 -2 -1 2];

k7m = (1/60)* [-1 -1 -1 -1 -1; ...
               2  2  2  2  2; ...
			   0  0  0  0  0; ...
			  -2 -2 -2 -2 -2; ...
			   1  1  1  1  1];

k8m = (1/140)* [-4 -2 0  2  4; ...
                 2  1 0 -1 -2; ...
				 4  2 0 -2 -4; ...
				 2  1 0 -1 -2; ...
			    -4 -2 0  2  4];

k9m = (1/140)* [-4  2  4  2 -4; ...
                -2  1  2  1 -2; ...
				 0  0  0  0  0; ...
				 2 -1 -2 -1  2; ...
			     4 -2 -4 -2  4];

k10m = (1/60)* [-1 2 0 -2 1; ...
                -1 2 0 -2 1; ...
				-1 2 0 -2 1; ...
				-1 2 0 -2 1; ...
			    -1 2 0 -2 1];

mask = zeros(m,n);
mask(3:m-2,3:n-2) = ones(m-4,n-4);

k1 = zeros(m,n);
for i1 = 3:m+2
  for j1 = 3:n+2
    k1(i1,j1) = sum(sum(k1m.*a1(i1-2:i1+2,j1-2:j1+2)));
  end
end
k1 = k1(3:m+2,3:n+2).*mask;

k2 = zeros(m,n);
for i1 = 3:m+2
  for j1 = 3:n+2
    k2(i1,j1) = sum(sum(k2m.*a1(i1-2:i1+2,j1-2:j1+2)));
  end
end
k2 = k2(3:m+2,3:n+2).*mask;

k3 = zeros(m,n);
for i1 = 3:m+2
  for j1 = 3:n+2
    k3(i1,j1) = sum(sum(k3m.*a1(i1-2:i1+2,j1-2:j1+2)));
  end
end
k3 = k3(3:m+2,3:n+2).*mask;

k4 = zeros(m,n);
for i1 = 3:m+2
  for j1 = 3:n+2
    k4(i1,j1) = sum(sum(k4m.*a1(i1-2:i1+2,j1-2:j1+2)));
  end
end
k4 = k4(3:m+2,3:n+2).*mask;

k5 = zeros(m,n);
for i1 = 3:m+2
  for j1 = 3:n+2
    k5(i1,j1) = sum(sum(k5m.*a1(i1-2:i1+2,j1-2:j1+2)));
  end
end
k5 = k5(3:m+2,3:n+2).*mask;

k6 = zeros(m,n);
for i1 = 3:m+2
  for j1 = 3:n+2
    k6(i1,j1) = sum(sum(k6m.*a1(i1-2:i1+2,j1-2:j1+2)));
  end
end
k6 = k6(3:m+2,3:n+2).*mask;

k7 = zeros(m,n);
for i1 = 3:m+2
  for j1 = 3:n+2
    k7(i1,j1) = sum(sum(k7m.*a1(i1-2:i1+2,j1-2:j1+2)));
  end
end
k7 = k7(3:m+2,3:n+2).*mask;

k8 = zeros(m,n);
for i1 = 3:m+2
  for j1 = 3:n+2
    k8(i1,j1) = sum(sum(k8m.*a1(i1-2:i1+2,j1-2:j1+2)));
  end
end
k8 = k8(3:m+2,3:n+2).*mask;

k9 = zeros(m,n);
for i1 = 3:m+2
  for j1 = 3:n+2
    k9(i1,j1) = sum(sum(k9m.*a1(i1-2:i1+2,j1-2:j1+2)));
  end
end
k9 = k9(3:m+2,3:n+2).*mask;

k10 = zeros(m,n);
for i1 = 3:m+2
  for j1 = 3:n+2
    k10(i1,j1) = sum(sum(k10m.*a1(i1-2:i1+2,j1-2:j1+2)));
  end
end
k10 = k10(3:m+2,3:n+2).*mask;
