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

어떤 아이디어?

Answer 1

Griddata는 스플라인 피팅을 사용합니다. 3 차 스플라인은 3 차 다항식과 동일하지 않습니다 (대신 모든 점에서 다른 3 차 다항식입니다).

데이터에 2D, 3 차 다항식을 맞추려면 다음과 같이하십시오. all 데이터 포인트를 사용하여 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()