Python三维多项式曲面拟合，顺序依赖

我目前正在使用其中有彗星图像的天文数据。由于拍摄时间（黄昏），我想删除这些图像中的背景天空渐变。我开发的第一个程序是从Matplotlib的“ginput”（x，y）中抽取用户选定的点，然后将每个坐标（z）的数据拉到一起，然后用SciPy的“griddata”格式化数据到一个新数组中。由于背景假设只是略有变化，所以我想将一个3d低阶多项式拟合到这组（x，y，z）点上。然而，“的GridData”不允许输入命令：Python三维多项式曲面拟合，顺序依赖

griddata(points,values, (dimension_x,dimension_y), method='nearest/linear/cubic')

上可使用的其他功能或开发一个执法机关方拟合，让我来控制顺序的方法的任何想法？

来源

2011-11-03 Sebastian De Pascuale

Griddata使用样条拟合。三阶样条与三阶多项式不同（相反，它是每个点上不同的三阶多项式）。

如果您只是想为您的数据拟合一个二维的三阶多项式，那么请执行如下操作，使用您的数据点的全部来估计16个系数。

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

def main(): 
    # Generate Data... 
    numdata = 100 
    x = np.random.random(numdata) 
    y = np.random.random(numdata) 
    z = x**2 + y**2 + 3*x**3 + y + np.random.random(numdata) 

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

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

    # Plot 
    plt.imshow(zz, extent=(x.min(), y.max(), x.max(), y.min())) 
    plt.scatter(x, y, c=z) 
    plt.show() 

def polyfit2d(x, y, z, order=3): 
    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()

enter image description here

来源

2011-11-03 15:56:08

这是一个非常优雅的解决问题的方法。我已经使用你的建议代码来适应一个椭圆抛物面，稍微修改我想分享的内容。我感兴趣的是仅以“z = a *（x-x0）** 2 + b *（y-y0）** 2 + c”的形式拟合线性解。我的更改的完整代码可以在[这里]看到（http://www.nublia.com/dev/stackoverflow/stow_polyfit2d.py）。 – regeirk

注意：对于numpy的最新版本，请参阅下面的@ klaus的答案。在我原来的答案'polyvander2d'等时不存在，但他们现在是要走的路。 –

这真的是一个三阶多项式吗？除非我对它有错误的理解，但是它不会有第6项的“X ** 3 * Y ** 3”这个词吗？ – maxymoo

的polyfit2d下实现使用可用numpy的方法numpy.polynomial.polynomial.polyvander2d和numpy.polynomial.polynomial.polyval2d

#!/usr/bin/env python3 

import unittest 


def polyfit2d(x, y, f, deg): 
    from numpy.polynomial import polynomial 
    import numpy as np 
    x = np.asarray(x) 
    y = np.asarray(y) 
    f = np.asarray(f) 
    deg = np.asarray(deg) 
    vander = polynomial.polyvander2d(x, y, deg) 
    vander = vander.reshape((-1,vander.shape[-1])) 
    f = f.reshape((vander.shape[0],)) 
    c = np.linalg.lstsq(vander, f)[0] 
    return c.reshape(deg+1) 

class MyTest(unittest.TestCase): 

    def setUp(self): 
     return self 

    def test_1(self): 
     self._test_fit(
      [-1,2,3], 
      [ 4,5,6], 
      [[1,2,3],[4,5,6],[7,8,9]], 
      [2,2]) 

    def test_2(self): 
     self._test_fit(
      [-1,2], 
      [ 4,5], 
      [[1,2],[4,5]], 
      [1,1]) 

    def test_3(self): 
     self._test_fit(
      [-1,2,3], 
      [ 4,5], 
      [[1,2],[4,5],[7,8]], 
      [2,1]) 

    def test_4(self): 
     self._test_fit(
      [-1,2,3], 
      [ 4,5], 
      [[1,2],[4,5],[0,0]], 
      [2,1]) 

    def test_5(self): 
     self._test_fit(
      [-1,2,3], 
      [ 4,5], 
      [[1,2],[4,5],[0,0]], 
      [1,1]) 

    def _test_fit(self, x, y, c, deg): 
     from numpy.polynomial import polynomial 
     import numpy as np 
     X = np.array(np.meshgrid(x,y)) 
     f = polynomial.polyval2d(X[0], X[1], c) 
     c1 = polyfit2d(X[0], X[1], f, deg) 
     np.testing.assert_allclose(c1, 
           np.asarray(c)[:deg[0]+1,:deg[1]+1], 
           atol=1e-12) 

unittest.main()

来源

2015-08-30 14:25:40

Python三维多项式曲面拟合，顺序依赖

回答

相关问题