7

I have some points in 3-space and I'd like to fit a quadratic surface through them.

I tried this code 

import itertools
import numpy as np
import matplotlib.pyplot as plt


def main():
    points = [ [ 175697888, -411724928, 0.429621160030365 ], [ 175697888, -411725144, 0.6078286170959473 ], [ 175698072, -411724640, 0.060898926109075546 ], [ 175698008, -411725360, 0.6184252500534058 ], [ 175698248, -411725720, 0.0771455243229866 ], [ 175698448, -411724456, -0.5925689935684204 ], [ 175698432, -411725936, -0.17584866285324097 ], [ 175698608, -411726152, -0.24736160039901733 ], [ 175698840, -411724360, -1.27967369556427 ], [ 175698800, -411726440, -0.21100902557373047 ], [ 175699016, -411726744, -0.12785470485687256 ], [ 175699280, -411724208, -2.472576856613159 ], [ 175699536, -411726688, -0.19858847558498383 ], [ 175699760, -411724104, -3.5765910148620605 ], [ 175699976, -411726504, -0.7432857155799866 ], [ 175700224, -411723960, -4.770215034484863 ], [ 175700368, -411726304, -1.2959377765655518 ], [ 175700688, -411723760, -6.518451690673828 ], [ 175700848, -411726080, -3.02254056930542 ], [ 175701160, -411723744, -7.941056251525879 ], [ 175701112, -411725896, -3.884831428527832 ], [ 175701448, -411723824, -8.661275863647461 ], [ 175701384, -411725720, -5.21607780456543 ], [ 175701704, -411725496, -6.181706428527832 ], [ 175701800, -411724096, -9.490276336669922 ], [ 175702072, -411724344, -10.066594123840332 ], [ 175702216, -411724560, -10.098011016845703 ], [ 175702256, -411724864, -9.619892120361328 ], [ 175702032, -411725160, -6.936516284942627 ] ]

    n = len(points)
    x, y, z = map(np.array, zip(*points))

    plt.figure()
    plt.subplot(1, 1, 1)

    # Fit a 3rd order, 2d polynomial
    m = polyfit2d(x,y,z, order=2)

    # Evaluate it on a grid...
    nx, ny = 100, 100
    xx, yy = np.meshgrid(np.linspace(x.min(), x.max(), nx), np.linspace(y.min(), y.max(), ny))
    zz = polyval2d(xx, yy, m)

    plt.scatter(xx, yy, c=zz, marker=2)
    plt.scatter(x, y, c=z)
    plt.show()

def polyfit2d(x, y, z, order=2):
    ncols = (order + 1)**2
    G = np.zeros((x.size, ncols))
    ij = itertools.product(range(order+1), range(order+1))
    for k, (i,j) in enumerate(ij):
        G[:,k] = x**i * y**j
    m, _, _, _ = np.linalg.lstsq(G, z)
    return m

def polyval2d(x, y, m):
    order = int(np.sqrt(len(m))) - 1
    ij = itertools.product(range(order+1), range(order+1))
    z = np.zeros_like(x)
    for a, (i,j) in zip(m, ij):
        z += a * x**i * y**j
    return z

main()

based on this answer: Python 3D polynomial surface fit, order dependent

But it actually gives the opposite result:

enter image description here

Look at the colour of the points compared to the surface. Any idea what I'm doing wrong?

EDIT: Update the code to remove the imshow showing that isn't the issue.

Community
  • 1
  • 1
nickponline
  • 25,354
  • 32
  • 99
  • 167
  • Looks like your coordinate system is flipped from bottom-left to top-right – OneCricketeer Jan 27 '16 at 19:26
  • 1
    The `extent` kwarg to `imshow` should be `[xmin, xmax, ymax, ymin]`, not [`xmin, xmax, ymin, ymax]`, unless you change the `origin` parameter. I had a typo w.r.t. `extent` in the linked answer, but despite the x/y mixup, note that the second half is `max, min`, not `min, max`. – Joe Kington Jan 27 '16 at 19:27
  • What is `polyval2d`? – Alex Alifimoff Jan 27 '16 at 19:32
  • Also, I have a mistake somewhere in the original answer. Your example will work if you take out the `y = -y` and use the `extent=[xmax, xmin, ymax, ymin]` that's in the original answer... However, that shouldn't be correct. – Joe Kington Jan 27 '16 at 19:36
  • @JoeKington I don't think so, I got rid of the imshow totally and the y=-y and just used another scatter plt.scatter(xx, yy, c=zz, marker=1) but I get this: https://www.dropbox.com/s/s2gn53u6cxl4kir/Screenshot%202016-01-27%2011.40.05.png?dl=0 – nickponline Jan 27 '16 at 19:40
  • @JoeKingston When I use a order-1 polynomial seems to work as intended though: https://www.dropbox.com/s/roapwsafox41sj3/Screenshot%202016-01-27%2011.40.38.png?dl=0 – nickponline Jan 27 '16 at 19:40
  • Seems that order 1, 2, 3 and 4 regressions give totally different fits. https://www.dropbox.com/s/o12s53oarmp2lsi/Screenshot%202016-01-27%2011.46.44.png?dl=0 the high ones may just be over fitting but the order 2 is opposite to the order 1. – nickponline Jan 27 '16 at 19:47
  • The root problem is that it isn't fitting things at all... Something strange is happening, probably because your x & y coordinates are much larger than your z-coordinates. Note that the range of `zz` is ~1e13, while the range of the input `z` is ~10. – Joe Kington Jan 27 '16 at 20:09

2 Answers2

8

There seems to be a problem with the floating point accuracy. I played with your code a bit and change the range of x and y made the least square solution work. Doing

x, y = x - x[0], y - y[0]

solved the accuracy issue. You can try:

import itertools
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
# from matplotlib import cbook
from matplotlib import cm
from matplotlib.colors import LightSource


def poly_matrix(x, y, order=2):
    """ generate Matrix use with lstsq """
    ncols = (order + 1)**2
    G = np.zeros((x.size, ncols))
    ij = itertools.product(range(order+1), range(order+1))
    for k, (i, j) in enumerate(ij):
        G[:, k] = x**i * y**j
    return G


points = np.array([[175697888, -411724928, 0.429621160030365],
                   [175697888, -411725144, 0.6078286170959473],
                   [175698072, -411724640, 0.060898926109075546],
                   [175698008, -411725360, 0.6184252500534058],
                   [175698248, -411725720, 0.0771455243229866],
                   [175698448, -411724456, -0.5925689935684204],
                   [175698432, -411725936, -0.17584866285324097],
                   [175698608, -411726152, -0.24736160039901733],
                   [175698840, -411724360, -1.27967369556427],
                   [175698800, -411726440, -0.21100902557373047],
                   [175699016, -411726744, -0.12785470485687256],
                   [175699280, -411724208, -2.472576856613159],
                   [175699536, -411726688, -0.19858847558498383],
                   [175699760, -411724104, -3.5765910148620605],
                   [175699976, -411726504, -0.7432857155799866],
                   [175700224, -411723960, -4.770215034484863],
                   [175700368, -411726304, -1.2959377765655518],
                   [175700688, -411723760, -6.518451690673828],
                   [175700848, -411726080, -3.02254056930542],
                   [175701160, -411723744, -7.941056251525879],
                   [175701112, -411725896, -3.884831428527832],
                   [175701448, -411723824, -8.661275863647461],
                   [175701384, -411725720, -5.21607780456543],
                   [175701704, -411725496, -6.181706428527832],
                   [175701800, -411724096, -9.490276336669922],
                   [175702072, -411724344, -10.066594123840332],
                   [175702216, -411724560, -10.098011016845703],
                   [175702256, -411724864, -9.619892120361328],
                   [175702032, -411725160, -6.936516284942627]])

ordr = 2  # order of polynomial
x, y, z = points.T
x, y = x - x[0], y - y[0]  # this improves accuracy

# make Matrix:
G = poly_matrix(x, y, ordr)
# Solve for np.dot(G, m) = z:
m = np.linalg.lstsq(G, z)[0]


# Evaluate it on a grid...
nx, ny = 30, 30
xx, yy = np.meshgrid(np.linspace(x.min(), x.max(), nx),
                     np.linspace(y.min(), y.max(), ny))
GG = poly_matrix(xx.ravel(), yy.ravel(), ordr)
zz = np.reshape(np.dot(GG, m), xx.shape)

# Plotting (see http://matplotlib.org/examples/mplot3d/custom_shaded_3d_surface.html):
fg, ax = plt.subplots(subplot_kw=dict(projection='3d'))
ls = LightSource(270, 45)
rgb = ls.shade(zz, cmap=cm.gist_earth, vert_exag=0.1, blend_mode='soft')
surf = ax.plot_surface(xx, yy, zz, rstride=1, cstride=1, facecolors=rgb,
                       linewidth=0, antialiased=False, shade=False)
ax.plot3D(x, y, z, "o")

fg.canvas.draw()
plt.show()

which gives 3DResult Plot]

To evaluate the quality of you fit read the documentation for np.linalg.lstsq(). The rank should be the size of your result vector and the residual divided by the number of data points gives the average error (distance between point and plane).

Dietrich
  • 5,241
  • 3
  • 24
  • 36
0

The expression on the right in

G[:,k] = x**i * y**j

is silently overflowing. x and y numpy arrays of integers, and therefore the result x**i * y**j is also an integer array. But when i and j are both 2, some of the products overflow, and wrap around to negative values.

Change your code to make x and y floating point. For example,

    points = np.array([ [ 175697888, -411724928, 0.429621160030365 ], [ 175697888, -411725144, 0.6078286170959473 ], [ 175698072, -411724640, 0.060898926109075546 ], [ 175698008, -411725360, 0.6184252500534058 ], [ 175698248, -411725720, 0.0771455243229866 ], [ 175698448, -411724456, -0.5925689935684204 ], [ 175698432, -411725936, -0.17584866285324097 ], [ 175698608, -411726152, -0.24736160039901733 ], [ 175698840, -411724360, -1.27967369556427 ], [ 175698800, -411726440, -0.21100902557373047 ], [ 175699016, -411726744, -0.12785470485687256 ], [ 175699280, -411724208, -2.472576856613159 ], [ 175699536, -411726688, -0.19858847558498383 ], [ 175699760, -411724104, -3.5765910148620605 ], [ 175699976, -411726504, -0.7432857155799866 ], [ 175700224, -411723960, -4.770215034484863 ], [ 175700368, -411726304, -1.2959377765655518 ], [ 175700688, -411723760, -6.518451690673828 ], [ 175700848, -411726080, -3.02254056930542 ], [ 175701160, -411723744, -7.941056251525879 ], [ 175701112, -411725896, -3.884831428527832 ], [ 175701448, -411723824, -8.661275863647461 ], [ 175701384, -411725720, -5.21607780456543 ], [ 175701704, -411725496, -6.181706428527832 ], [ 175701800, -411724096, -9.490276336669922 ], [ 175702072, -411724344, -10.066594123840332 ], [ 175702216, -411724560, -10.098011016845703 ], [ 175702256, -411724864, -9.619892120361328 ], [ 175702032, -411725160, -6.936516284942627 ] ])
    x, y, z = points.T
Warren Weckesser
  • 110,654
  • 19
  • 194
  • 214