function [exists, error, height, radius] = q2par(f, K, c, b, a)
% [exists, error, height, radius] = q2par(f, K, c, b, a)
% Compute error  = SUM(c*m^2+b*m+a-f(n+m))^2  m= -k to k
%         height = b^2/(4*c) - b^2/(2*c) + a
%         radius = -sqrt(b^2-4ac)/(2*c)
[m n] = size(f);
k = round((K-1)/2);
for indexi = 1:m
  for indexj = 1:n
    error(indexi,indexj) = 0; height(indexi,indexj) = 0; 
    radius(indexi,indexj) = 0; exists(indexi,indexj) = 0;
    if c(indexi,indexj) < 0.0
      x0 = (-b(indexi,indexj)/(2.0*c(indexi,indexj)));
      if abs(x0) <= 0.8
        exists(indexi,indexj) = 255.0;
        e = 0.0;
        for l = -k:k
          e = e + (c(indexi,indexj)*l*l + b(indexi,indexj)*l +  ...
              a(indexi,indexj) - f(indexi+l,indexj))^2;
        end
        error(indexi,indexj) = e;
        height(indexi,indexj) = ...
            (-b(indexi,indexj)^2)/(4*c(indexi,indexj)) + a(indexi,indexj);
        radius(indexi,indexj) = -sqrt(b(indexi,indexj)^2 - ...
            4*a(indexi,indexj)*c(indexi,indexj))/(2*c(indexi,indexj));
        if height(indexi,indexj) > 255
          error(indexi,indexj) = 0; height(indexi,indexj) = 0; 
          radius(indexi,indexj) = 0; exists(indexi,indexj) = 0;
        end
        if radius(indexi,indexj) > 50
          error(indexi,indexj) = 0; height(indexi,indexj) = 0; 
          radius(indexi,indexj) = 0; exists(indexi,indexj) = 0;
        end				
      end
    end
  end
end
