Fitting "quadratic" surface to data points in 3D

Question

I am trying to fit a quadratic plane to a cloud of data points in python. My plane function is of the form

   f(x,y,z) = a*x**2 + b*y**2 + c*x*y + d*x + e*y + f - z

Currently, my data points do not have errors associated with them, however, some errors can be assumed if necessary. Following suggestions from here, I work out the vertical distance from a point p0=(x0,y0,z0) (which correspond to my data points) to a point on the plane p=(x,y,z) following this method. I then end up with

def vertical_distance(params,p0):
    *** snip ***
    nominator = f + a*x**2 + b*y**2 + c*x*y - x0*(2*a*x-c*y-d) - y0*(2*b*y-c*x-e) + z0
    denominator = sqrt((2*a*x+c*y+d)**2 + (2*b*y+c*x+e)**2 + (-1)**2)
    return nominator/denominator

Ultimately, I think it is the vertical_distance function that I need to minimise. I can happily feed a list of starting parameters (params) and the array of data points to it in two dimensions, however, I am unsure on how to achieve this in 3D. ODR pack seems to only allow data containing x,y or two dimensions. Furthermore, how do I implement the points on the plane (p) into the minimising routine? I guess that during the fit operations the points vary according to the parameter optimisation and thus the exact equation of the plane at that very moment.

Those links are for linear planes. – meowgoesthedog Mar 01 '18 at 17:04 — meowgoesthedog, Mar 01 '18 at 17:04

score 0 · Accepted Answer · answered Mar 01 '18 at 18:13

I guess that «quadratic surface» would be a more correct term than «plane».

And the problem is to fit z = ax^2 + by^2 + cxy + dx + ey + f to given set of points P.

To do that via optimization you need to formulate residual function (for instance, vertical distance).

For each 3D points p from P residual is

|p_2 – ap_0^2 + bp_1^2 + c*p_0*p_1 + dp_0 + ep_1 + f|

You need to minimize all residuals, i.e. sum of square of them, variating parameters a…f.

The following code technically should solve above problem. But fitting the problem is multi-extremal and such routine may fail to find right set of parameters without good starting point or globalization of search.

import numpy
import scipy.optimize

P = numpy.random.rand(3,10) # given point set

def quadratic(x,y, a, b, c, d, e, f): 
  #fit quadratic surface
  return a*x**2 + b*y**2 + c*x*y + d*x + e*y + f

def residual(params, points):
  #total residual
  residuals = [
    p[2] - quadratic(p[0], p[1],
    params[0], params[1], params[2], params[3], params[4], params[5]) for p in points]

  return numpy.linalg.norm(residuals)

result = scipy.optimize.minimize(residual, 
                                 (1, 1, 0, 0, 0, 0),#starting point
                                 args=P)

Fitting "quadratic" surface to data points in 3D

1 Answers1