function [pa,pb,pc,pd,pe] = parabg2(f,Ky,Kx,coefm,ridges)
% [pa,pb,pc,pd,pe] = parabg2(f,Ky,Kx,coefm,ridges)
[m n] = size(f);
pa = zeros(m,n); pb = zeros(m,n); pc = zeros(m,n);
pd = zeros(m,n); pe = zeros(m,n);
kx = (Kx-1)/2;
ky = (Ky-1)/2;
for i1 = 1+ky:m-ky
  for j1 = 1+kx:n-kx
   if ridges(i1,j1) > 0
    sumy2 = 0; sumx2 = 0; sumx = 0; sumy = 0; sumf = 0;
    for i2 = -ky:ky
      y2 = i2*i2;
      for j2 = -kx:kx
        x2 = j2*j2;
        fval = f(i1+i2,j1+j2);
        sumf = sumf+fval;
        sumy2 = sumy2 + y2*fval;
        sumx2 = sumx2 + x2*fval;
        sumx = sumx + j2*fval;
        sumy = sumy + i2*fval;
      end
    end
    params = coefm\[sumx2, sumy2, sumx, sumy, sumf]';
    pa(i1,j1) = params(1);
    pb(i1,j1) = params(2);
    pc(i1,j1) = params(3);
    pd(i1,j1) = params(4);
    if abs(pa(i1,j1))>eps
      pc(i1,j1) = -pc(i1,j1)/(2*pa(i1,j1));
    end
    if abs(pb(i1,j1))>eps
      pd(i1,j1) = -pd(i1,j1)/(2*pb(i1,j1));
    end
    pe(i1,j1) = params(5);
   end
  end
end
