function [pa,pb,pc] = parabg(f,Kx,Ky,coefm)
% [pa,pb,pc] = parabg(f,Kx,Ky,coefm)
[m n] = size(f);
pa = zeros(m,n); pb = zeros(m,n); pc = zeros(m,n);
kx = (Kx-1)/2;
ky = (Ky-1)/2;
onek = ones(Ky,Kx);
sumf = conv2(f,onek);
for i1 = 1+ky:m-ky
  for j1 = 1+kx:n-kx
    sumy2 = 0; sumx2 = 0;
    for i2 = -ky:ky
      y2 = i2*i2;
      for j2 = -kx:kx
        x2 = j2*j2;
        fval = f(i1+i2,j1+j2);
        sumy2 = sumy2 + y2*fval;
        sumx2 = sumx2 + x2*fval;
      end
    end
    params = coefm\[sumf(i1,j1), sumy2, sumx2]';
    pa(i1,j1) = params(1);
    pb(i1,j1) = params(2);
    pc(i1,j1) = params(3);
  end
end
