Ajuste de superficie polinomial 3D de Python, dependiente del pedido

Question

Actualmente estoy trabajando con datos astronómicos entre los que tengo imágenes de cometas. Me gustaría eliminar el gradiente de fondo del cielo en estas imágenes debido al tiempo de captura (crepúsculo). El primer programa que desarrollé para hacerlo tomó los puntos seleccionados por el usuario de la "ginput" de Matplotlib (x, y) extrajo los datos para cada coordenada (z) y luego cuadriculó los datos en una nueva matriz con "griddata" de SciPy.

Como se supone que el fondo varía ligeramente, me gustaría ajustar un polinomio de orden bajo de 3D a este conjunto de puntos (x, y, z). Sin embargo, el "gridData" no permite la entrada de una orden:

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

¿Alguna idea sobre otra función que se puede utilizar o un método para desarrollar un LEA cuadrados encajan que me permitirá controlar el orden?

Answer 1

Griddata utiliza una conexión spline. Una spline de 3er orden no es lo mismo que un polinomio de 3er orden (en cambio, es un polinomio diferente de 3er orden en cada punto).

Si solo desea ajustar un polinomio 2D de 3er orden a sus datos, haga lo siguiente para estimar los 16 coeficientes usando todos de sus puntos de datos.

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

La siguiente aplicación de polyfit2d utiliza los métodos disponibles numpy numpy.polynomial.polynomial.polyvander2d y 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()