Python 3D多項式の表面適合、次数依存

Question

私は現在、私が彗星の画像を持っている天文データを扱っています。私はキャプチャ（夕暮れ）の時間のために、これらの画像の背景の空の勾配を削除したいと思います。私がこれを行うために開発した最初のプログラムでは、Matplotlibの "ginput"（x、y）からユーザーが選択したポイントが各座標（z）のデータを取得し、SciPyの "griddata"を使って新しい配列にグリッドされました。

バックグラウンドはわずかにしか変化しないと仮定されているので、この（x、y、z）点のセットに3d低次多項式を当てはめたいと思います。しかし、「griddataは、」入力順序を許可しない：

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

のLEA-正方形を開発するために使用することができる他の関数やメソッド上の任意のアイデアそれは私が順序を制御することができます収まりますか？

Answer 1

Griddataはスプラインフィッティングを使用します。 3次スプラインは3次多項式と同じではありません（代わりに、すべての点で異なる3次多項式です）。

2次、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()