पायथन 3 डी बहुपद सतह फिट, आदेश निर्भर

Question

मैं वर्तमान में खगोलीय डेटा के साथ काम कर रहा हूं जिसमें मेरे पास धूमकेतु छवियां हैं। कैप्चर (ट्वाइलाइट) के समय के कारण मैं इन छवियों में पृष्ठभूमि आकाश ढाल को हटाना चाहता हूं। ऐसा करने वाला पहला प्रोग्राम मैंने मैप्लोट्लिब के "जीनपूट" (एक्स, वाई) से उपयोगकर्ता चुने गए अंकों को प्रत्येक समन्वय (जेड) के लिए डेटा खींच लिया और फिर साइपी के "ग्रिडाटा" के साथ डेटा को एक नई सरणी में ग्रिड किया।

चूंकि पृष्ठभूमि को केवल थोड़ा अलग माना जाता है, इसलिए मैं इस एक्स (वाई, वाई, जेड) बिंदुओं के सेट के लिए 3 डी कम ऑर्डर बहुपद फिट करना चाहता हूं। हालांकि, "griddata" किसी इनपुट आदेश के लिए अनुमति नहीं देता:

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

एक और समारोह है कि इस्तेमाल किया जा सकता है या विकासशील एक एल ई ऐज़ वर्गों फिट है कि मुझे आदेश पर नियंत्रण कर पाएँगे के लिए एक विधि पर कोई भी विचार?

Answer 1

Griddata एक स्पलीन फिटिंग का उपयोग करता है। एक तीसरी ऑर्डर स्पलीन एक तीसरी ऑर्डर बहुपद के समान नहीं है (इसके बजाय, यह प्रत्येक बिंदु पर एक अलग तीसरा ऑर्डर बहुपद है)।

यदि आप अपने डेटा के लिए 2 डी, तीसरे क्रम बहुपद को फिट करना चाहते हैं, तो अपने डेटा बिंदुओं के सभी का उपयोग करके 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()