Ajuste de superficie polinomial Python 3D, depende del orden

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

Dado que se supone que el fondo varía solo ligeramente, me gustaría ajustar un polinomio de bajo orden 3D a este conjunto de puntos (x, y, z). Sin embargo, el “griddata” no permite un orden de entrada:

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

¿Alguna idea sobre otra función que pueda usarse o un método para desarrollar un ajuste de leas-square que me permita controlar el pedido?

Related of "Ajuste de superficie polinomial Python 3D, depende del orden"

Griddata utiliza un ajuste spline. Una spline de tercer orden no es lo mismo que un polinomio de tercer orden (en cambio, es un polinomio de tercer orden diferente en cada punto).

Si solo desea ajustar un polinomio 2D de 3er orden a sus datos, haga algo como lo siguiente para estimar los 16 coeficientes usando todos 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() 

introduzca la descripción de la imagen aquí

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

Si alguien está buscando un polinomio de un orden específico (en lugar de los polinomios en los que la potencia más alta es igual a la order , puede hacer este ajuste a la polyfit y polyval la respuesta aceptada polyval :

en lugar de:

 ij = itertools.product(range(order+1), range(order+1)) 

el cual, para el order=2 da [(0, 0), (0, 1), (0, 2), (1, 0), (1, 1), (1, 2), (2, 0), (2, 1), (2, 2)] (también conocido como polinomio de hasta 4º grado), puede usar

 def xy_powers(order): powers = itertools.product(range(order + 1), range(order + 1)) return [tup for tup in powers if sum(tup) <= order] 

Esto devuelve [(0, 0), (0, 1), (0, 2), (1, 0), (1, 1), (2, 0)] para order=2

De acuerdo con el principio de los cuadrados mínimos , e imita el estilo de Kington, mientras mueves el argumento m al argumento m_1 y el argumento m_2.

 import numpy as np import matplotlib.pyplot as plt import itertools # w = (Phi^T Phi)^{-1} Phi^T t # where Phi_{k, j + i (m_2 + 1)} = x_k^i y_k^j, # t_k = z_k, # i = 0, 1, ..., m_1, # j = 0, 1, ..., m_2, # k = 0, 1, ..., n - 1 def polyfit2d(x, y, z, m_1, m_2): # Generate Phi by setting Phi as x^iy^j nrows = x.size ncols = (m_1 + 1) * (m_2 + 1) Phi = np.zeros((nrows, ncols)) ij = itertools.product(range(m_1 + 1), range(m_2 + 1)) for h, (i, j) in enumerate(ij): Phi[:, h] = x ** i * y ** j # Generate t by setting t as Z t = z # Generate w by solving (Phi^T Phi) w = Phi^T t w = np.linalg.solve(Phi.T.dot(Phi), (Phi.T.dot(t))) return w # t' = Phi' w # where Phi'_{k, j + i (m_2 + 1)} = x'_k^i y'_k^j # t'_k = z'_k, # i = 0, 1, ..., m_1, # j = 0, 1, ..., m_2, # k = 0, 1, ..., n' - 1 def polyval2d(x_, y_, w, m_1, m_2): # Generate Phi' by setting Phi' as x'^i y'^j nrows = x_.size ncols = (m_1 + 1) * (m_2 + 1) Phi_ = np.zeros((nrows, ncols)) ij = itertools.product(range(m_1 + 1), range(m_2 + 1)) for h, (i, j) in enumerate(ij): Phi_[:, h] = x_ ** i * y_ ** j # Generate t' by setting t' as Phi' w t_ = Phi_.dot(w) # Generate z_ by setting z_ as t_ z_ = t_ return z_ if __name__ == "__main__": # Generate x, y, z n = 100 x = np.random.random(n) y = np.random.random(n) z = x ** 2 + y ** 2 + 3 * x ** 3 + y + np.random.random(n) # Generate w w = polyfit2d(x, y, z, m_1=3, m_2=2) # Generate x', y', z' n_ = 1000 x_, y_ = np.meshgrid(np.linspace(x.min(), x.max(), n_), np.linspace(y.min(), y.max(), n_)) z_ = np.zeros((n_, n_)) for i in range(n_): z_[i, :] = polyval2d(x_[i, :], y_[i, :], w, m_1=3, m_2=2) # Plot plt.imshow(z_, extent=(x_.min(), y_.max(), x_.max(), y_.min())) plt.scatter(x, y, c=z) plt.show() 

introduzca la descripción de la imagen aquí