function [a,b,c] = parab(f);
% [a,b,c] = parab(f);
[m n] = size(f);
a = zeros(m,n);
a = zeros(m,n);
a = zeros(m,n);
con = [6,4,6; 4,6,6; 6,6,9];
for i1 = 2:m-1
  for j1 = 2:n-1
    params(1) = 0; params(2) = 0; params(3) = 0;
    x1 = 0;
    x2 = 0;
    x3 = 0;
    for i2 = 1:3
      i3 = i2-2;
      for j2 = 1:3
        j3 = j2-2;
        fval = f(i1+i3,j1+j3);
        x1 = x1 + (j3^2)*fval;
        x2 = x2 + (i3^2)*fval;
        x3 = x3 + fval;
      end
    end
    params = con\[x1,x2,x3]';
  a(i1,j1) = params(1);
  b(i1,j1) = params(2);
  c(i1,j1) = params(3);
  end
end
