function [k1,k2,k3,k4,k5,k6,k7,k8,k9,k10] = bicubicm(f,l)
%
% bicubicm: fits bicubic function at each point to an image
%           uses least squares fit
% On input:
%     f: m by n array
%     l: window size for least squares fit at each point
% On output:
%     k1..k10: bicubic coefficients (see Haralick and Shapiro, Computer Vision)
% Call:
%     [k1,k2,k3,k4,k5,k6,k7,k8,k9,k10] = bicubicm(f,l)
%

[m,n] = size(f);
w = (l-1)/2;
for i1 = 1:m
  for j1 = 1:n           % for every point
    f1 = zeros(l*l,1);
    index = 0;
    for r = -w:w
      for c = -w:w            % build equations in l by l window
        index = index+1;
        A(index,1) = 1;
        A(index,2) = r;
        A(index,3) = c;
        A(index,4) = r*r;
        A(index,5) = r*c;
        A(index,6) = c*c;
        A(index,7) = r*r*r;
        A(index,8) = r*r*c;
        A(index,9) = r*c*c;
        A(index,10) = c*c*c;
        ind1 = i1+r;  ind2 = j1+c;
        if (ind1<1) | (ind1>m)
          ind1 = i1;
        end
        if (ind2<1) | (ind2>n)
          ind2 = j1;
        end
        f1(index) = f(ind1,ind2);
      end
    end
    k = A\f1;        % solve system
    k1(i1,j1) = k(1);  k2(i1,j1) = k(2);  k3(i1,j1) = k(3); k4(i1,j1) = k(4);
    k5(i1,j1) = k(5);  k6(i1,j1) = k(6);
    k7(i1,j1) = k(7); k8(i1,j1) = k(8);
    k9(i1,j1) = k(9); k10(i1,j1) = k(10);
  end
end
