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

关于可能使用的另一个函数的任何想法或开发leas-square拟合的方法可以让我控制订单？

Answer 1

Griddata使用样条拟合.三阶样条与三阶多项式不同(相反,它是每个点的不同的三阶多项式).

如果您只想将2D,3阶多项式拟合到数据中,请执行以下操作,以使用所有数据点估计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()

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