Python 3D多项式曲面拟合，依赖于顺序。

Question

我目前处理天文数据，其中包括彗星图像。由于拍摄时间（黄昏），我希望去除这些图像中的背景天空梯度。我开发了第一个程序来实现这个目标，该程序使用Matplotlib的"ginput"获取用户选择的点（x、y），然后提取每个坐标（z）的数据，并使用SciPy的"griddata"将数据网格化到新的数组中。

由于假定背景只略微变化，因此我想将3D低阶多项式拟合到这组（x、y、z）点上。但是，"griddata"不允许输入顺序：

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

有没有其他可用的函数或开发最小二乘拟合并允许我控制顺序的方法？

- Sebastian De Pascuale

---

你能解释一下“输入顺序”是什么意思吗？我不明白它如何适用于你的问题。谢谢。 - Owen

Answer 1

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()

- Joe Kington

---

这是一个非常优雅的解决问题的方案。我用你建议的代码进行了微小修改以适应一个椭圆抛物面的拟合。我只对形式为z = a*(x-x0)**2 + b*(y-y0)**2 + c的线性解决方案感兴趣。我修改后的完整代码可以在这里看到：http://www.nublia.com/dev/stackoverflow/stow_polyfit2d.py。 - regeirk

注意：对于最近版本的numpy，请参见下面@klaus的答案。在我最初回答时，polyvander2d等函数并不存在，但现在它们是首选的方式。 - Joe Kington

这真的是一个三次多项式吗？除非我理解错了，不然它不会有一个阶数为6的 X ** 3 * Y ** 3 项吗？ - maxymoo

你是正确的 @maxymoo。我猜作者的意思是每个变量的最高幂次为3。 - Tal J. Levy

regeirk的代码链接已失效。 - Steve Boege

Answer 2

下面的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()

- klaus se

---

deg=np.array([deg,deg]) 而不是（如果原始的 deg 是整数？） - jtlz2

你能否扩展 polyfit2d 函数，使其像 numpy.polynomial.polynomial.polyfit 一样接受一个权重数组？ - Owen

@Owen 添加权重将会像这样：c = np.linalg.lstsq(vander * weights, f * weights)[0] - klaus se

Answer 3

根据最小二乘法的原理，并模仿Kington的风格，将参数m移动到参数m_1和参数m_2。

import numpy as np
import matplotlib.pyplot as plt

import itertools


# w = (Phi^T Phi)^{-1} Phi^T t
# where Phi_{k, j + i (m_2 + 1)} = x_k^i y_k^j,
#       t_k = z_k,
#           i = 0, 1, ..., m_1,
#           j = 0, 1, ..., m_2,
#           k = 0, 1, ..., n - 1
def polyfit2d(x, y, z, m_1, m_2):
    # Generate Phi by setting Phi as x^i y^j
    nrows = x.size
    ncols = (m_1 + 1) * (m_2 + 1)
    Phi = np.zeros((nrows, ncols))
    ij = itertools.product(range(m_1 + 1), range(m_2 + 1))
    for h, (i, j) in enumerate(ij):
        Phi[:, h] = x ** i * y ** j
    # Generate t by setting t as Z
    t = z
    # Generate w by solving (Phi^T Phi) w = Phi^T t
    w = np.linalg.solve(Phi.T.dot(Phi), (Phi.T.dot(t)))
    return w


# t' = Phi' w
# where Phi'_{k, j + i (m_2 + 1)} = x'_k^i y'_k^j
#       t'_k = z'_k,
#           i = 0, 1, ..., m_1,
#           j = 0, 1, ..., m_2,
#           k = 0, 1, ..., n' - 1
def polyval2d(x_, y_, w, m_1, m_2):
    # Generate Phi' by setting Phi' as x'^i y'^j
    nrows = x_.size
    ncols = (m_1 + 1) * (m_2 + 1)
    Phi_ = np.zeros((nrows, ncols))
    ij = itertools.product(range(m_1 + 1), range(m_2 + 1))
    for h, (i, j) in enumerate(ij):
        Phi_[:, h] = x_ ** i * y_ ** j
    # Generate t' by setting t' as Phi' w
    t_ = Phi_.dot(w)
    # Generate z_ by setting z_ as t_
    z_ = t_
    return z_


if __name__ == "__main__":
    # Generate x, y, z
    n = 100
    x = np.random.random(n)
    y = np.random.random(n)
    z = x ** 2 + y ** 2 + 3 * x ** 3 + y + np.random.random(n)

    # Generate w
    w = polyfit2d(x, y, z, m_1=3, m_2=2)

    # Generate x', y', z'
    n_ = 1000
    x_, y_ = np.meshgrid(np.linspace(x.min(), x.max(), n_),
                         np.linspace(y.min(), y.max(), n_))
    z_ = np.zeros((n_, n_))
    for i in range(n_):
        z_[i, :] = polyval2d(x_[i, :], y_[i, :], w, m_1=3, m_2=2)

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

- Togtokh Khorchin

Answer 4

如果有人想要拟合特定阶数的多项式（而不是最高次幂等于order的多项式），您可以对已接受答案中的polyfit和polyval进行以下调整：

改为：

ij = itertools.product(range(order+1), range(order+1))

对于order=2，它会返回[(0, 0), (0, 1), (0, 2), (1, 0), (1, 1), (1, 2), (2, 0), (2, 1), (2, 2)]（即多项式的4次方），您可以使用

def xy_powers(order):
    powers = itertools.product(range(order + 1), range(order + 1))
    return [tup for tup in powers if sum(tup) <= order]

对于order=2，这将返回[(0, 0), (0, 1), (0, 2), (1, 0), (1, 1), (2, 0)]

- Scott Staniewicz