Phù hợp với bề mặt đa thức Python 3D, thứ tự phụ thuộc

Question

Tôi hiện đang làm việc với dữ liệu thiên văn trong đó tôi có hình ảnh sao chổi. Tôi muốn loại bỏ gradient nền trong những hình ảnh này do thời gian chụp (chạng vạng). Chương trình đầu tiên tôi phát triển để làm như vậy đã lấy điểm người dùng được lựa chọn từ "ginput" của Matplotlib (x, y) kéo dữ liệu cho mỗi tọa độ (z) và sau đó gridded dữ liệu trong một mảng mới với SciPy "griddata".

Vì nền được giả định chỉ thay đổi một chút, tôi muốn phù hợp với đa thức bậc thấp 3d cho tập hợp các điểm (x, y, z) này. Tuy nhiên, "griddata" không cho phép cho một trật tự đầu vào:

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

Bất kỳ ý tưởng về một chức năng mà có thể được sử dụng hoặc một phương pháp để phát triển một leas-ô vuông phù hợp mà sẽ cho phép tôi để kiểm soát trật tự?

Answer 1

Griddata sử dụng khớp nối spline. Một spline thứ 3 không giống với đa thức bậc 3 (thay vào đó, nó là một đa thức bậc 3 khác nhau ở mọi điểm).

Nếu bạn chỉ muốn phù hợp với đa thức bậc 2, thứ 2 cho dữ liệu của mình, hãy thực hiện như sau để ước tính 16 hệ số sử dụng tất cả điểm dữ liệu của bạn.

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

Việc thực hiện sau đây của polyfit2d sử dụng các phương pháp NumPy sẵn numpy.polynomial.polynomial.polyvander2d và 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()