Fitting a quadric function to a set of points

[daviddoria]

I have a set of points and I want to find a "best fit quardric surface" through the points. I did the following:

1) Assume the function is in the form:
$$a x^2 + b y^2 + c xy + d x + e y + f = z$$

2) Make a nx6 matrix of the points put into (1), that is A=:
x1^2 y1^2 x1y1 x1 y1 1
x2^2 y2^2 x2y2 x2 y2 1
...
xn^2 yn^2 xnyn xn yn 1

3) Make an nx1 matrix of the z values of the points b =
z1
z2
...
zn

4) Solve $$M [a b c d e f]^T = b$$

Now I have a function that will plot a quadratic surface given the following a0 - a9:

$$F(x,y,z) = a0*x^2 + a1*y^2 + a2*z^2 + a3*x*y + a4*y*z + a5*x*z + a6*x + a7*y + a8*z + a9$$

clearly a0 = a, a1 = b, a2 = 0, a3 = c, a4 = 0, a5 = 0, a6 = d, a7 = e

But I'm not sure about a8 and a9. I think a9=f, but what do I do with the coefficient of z (a8) in F(x,y,z) ?

Thanks,

Dave

 

[hotvette]

Steps 1-3 are fine, but the equation in step 4 is overdetermined and can't be solved. If by 'best fit' you mean ordinary least squares, the applicable equation for step 4 is:
$$M^TM [a b c d e f]^T = M^Tb$$

The last part with the 2nd equation I don't follow at all.

 

[daviddoria]

1) Hm... isn't it just the standard least squares solution, or
$$Ax = b$$
$$x = A^{-1}b$$

or in this case since A is not invertible
$$x = pinv(A) b$$
where pinv is the pseudo inverse?

2) The second part - I was saying that in order to "see" this fit, I have a function in a toolbox that will plot implicit functions (I guess the zero level set?) in the form F(x,y,z). I was asking how the unknowns that I solved for in steps 1-3 would fit into that F(x,y,z) form.

Thanks for the help so far!

Dave

 

[Zaphos]

The coefficient a8 should be -1. So you're subtracting z from both sides of your original equation. a9 should be f like you said.

 

[daviddoria]

Where did you get $$a8 = -1$$ ?

It seems to work with some point sets, but for others (points sampled from a sphere, for example) the result is wrong (it says all of the coefficients are 0). Is there a reason that would happen?

Thanks,

David

 

[Zaphos]

Where did you get $$a8 = -1$$ ?

It seems to work with some point sets, but for others (points sampled from a sphere, for example) the result is wrong (it says all of the coefficients are 0). Is there a reason that would happen?

Thanks,

David

The a8 = -1 is because a8 is the coefficient of z, and your original equation is f(x,y) = z while your plotting function's equation is f(x,y,z) = 0, so we need to subtract z from both sides of the original to get the equations to match up. That gives f(x,y) - z = 0, so the coefficient of z is -1.

It giving all zero coefficients for a sphere is probably a problem with your fitting method, not with the plotting method. You're fitting f(x,y) to the point set, and that's fundamentally a 'height field' where z is only allowed to have one value for a given x,y. A sphere is not a height field; it will have multiple values for some x,y pairs. If your sphere point set is on the plane z=0, then the best fitting 'height field' to that point set actually is the plane z=0, so all the coefficients your solver solves for should come out to be zero. If you want to solve for a general quadric in which all 9 coefficients are free to vary, you need to use a different fitting method.

 

[daviddoria]

Zaphos,

You have been very helpful, thanks! You are right on - I just thought it was an "error", but in fact it makes a lot of sense that the plane is the best fit to a set of sphere points.

Do you have a recommendation of a similar non-height field type of fitting procedure?

Thanks,

Dave

 

[Zaphos]

Unfortunately that fitting problem is a lot trickier ...

With the height field method, your 'error' was 'distance in z from the point to the surface' which falls out quite naturally from the function f(x,y) = z. But when you have a surface implicitly defined by f(x,y,z) there is no similarly convenient and meaningful error measure. Instead you have two basic options:

(a) Use 'geometric distance,' where your distance is the actual shortest distance from the point to the surface. This is the best approach in terms of accuracy, but unfortunately there isn't afaik a simple expression linear for that shortest geometric distance to a quadric, so you end up with a non-linear problem and typically need some iterative approach to solving it: you guess a set of coefficients, find the points on that surface which are closest to your data points, update your coefficients based on that matching, and repeat until satisfied. You might check "Orthogonal Distance Fitting of Implicit Curves and Surfaces" for a detailed algorithm -- http://www.springerlink.com/content/a8hph1jrfkp9u92h/

(b) Use 'algebraic distance,' where your distance is simply defined by f(x,y,z), since the f() = 0 on the surface and is non-zero away from it. This is much easier -- you could extend your current solution to do it -- but also degenerate, as you can see that you can shrink the error of fit trivially just by scaling down all the coefficients uniformly, even though that doesn't change the actual surface at all. So you need to introduce some additional constraint on how those coefficients are allowed to change, which in turn may introduce singularities ... so some care need be taken with the exact constraints used. In a few special cases, eg for spheres, you can find a constraint that works out just right, and gives you the same result as geometric distance fitting. But I don't think there's quite such a nice constraint for quadrics, so fundamentally it's not going to give the 'best fit' directly.
For more on algebraic surface fitting I like Pratt's paper "Direct Least-Squares Fitting of Algebraic Surfaces" -- http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.18.6178 -- which has the geometric-equivalent solution for spheres and also does suggest a related iterative approach for quadrics.