पायथन में एक बिंदु और एक बाइबिक स्पलीन सतह के बीच की दूरी के लिए मुझे तेजी से अनुमान कैसे मिल सकता है?

Question

पायथन में एक बिंदु और एक बाइबिक स्पलीन सतह के बीच की दूरी के लिए मुझे तेजी से अनुमान कैसे मिल सकता है? क्या कोई मौजूदा समाधान है जिसे मैं साइपी, न्यूपी, या कुछ अन्य पैकेज में लाभ उठा सकता हूं?

मैं सतह इस के रूप में एक bicubic प्रक्षेप द्वारा परिभाषित किया गया है:

import numpy as np 
import scipy.interpolate 

# Define regular grid surface 
xmin,xmax,ymin,ymax = 25, 125, -50, 50 
x = np.linspace(xmin,xmax, 201) 
y = np.linspace(ymin,ymax, 201) 
xx, yy = np.meshgrid(x, y) 
z_ideal = (xx**2 + yy**2)/400 
z_ideal += z_ideal + np.random.uniform(-0.5, 0.5, z_ideal.shape) 
s_ideal = scipy.interpolate.interp2d(x, y, z_ideal, kind='cubic') 

और मुझे लगता है कि सतह के कुछ मापा अंक मिल गया है:

# Fake some measured points on the surface 
z_measured = z_ideal + np.random.uniform(-0.1, 0.1, z_ideal.shape) 
s_measured = scipy.interpolate.interp2d(x, y, z_measured, kind='cubic') 
p_x = np.random.uniform(xmin,xmax,10000) 
p_y = np.random.uniform(ymin,ymax,10000) 
p_z = s_measured(p_x, p_y) 

मैं निकटतम लगाना चाहते हैं p में प्रत्येक बिंदु पर सतह s_ideal पर इंगित करें। एक सामान्य मामले में जंगली रूप से भिन्न splines के लिए कई समाधान हो सकते हैं, इसलिए मैं समस्या को उन सतहों तक सीमित कर रहा हूं जिन्हें ज़ेड के साथ बिंदु के प्रक्षेपण के आसपास के केवल एक समाधान के लिए जाना जाता है। यह माप या सतह परिभाषा बिंदुओं की एक छोटी संख्या नहीं है, इसलिए मैं 1E-5 पर सटीकता की कीमत पर भी गति को अनुकूलित करना चाहता हूं।

1. उपयोग pt = [p_x, p_y, p_z] प्रारंभिक परीक्षण बिंदु के रूप में, जहां p_z = s_ideal(pt)
2. ढलान की गणना (:

   विधि जो मन में आता एक ढाल वंश दृष्टिकोण का उपयोग और प्रत्येक माप बिंदु p के लिए की तरह कुछ करने के लिए है स्पर्श) वेक्टर m = [ m_x, m_y ]pt

3. में से वेक्टर r की गणना ptp रहे हैं: r = p - pt
012,
4. यदि कोण thetar और m के बीच 90 डिग्री की कुछ सीमा के भीतर है, तो pt अंतिम बिंदु है।
5. अन्यथा, अद्यतन pt के रूप में:

r_len = numpy.linalg.norm(r) 
dx = r_len * m_x 
dy = r_len * m_y 
if theta > 90: 
    pt = [ p_x + dx, p_y + dy ] 
else: 
    pt = [ p_x - dx, p_y - dy ] 

मैं एक विधि -1 डी मामले के लिए एक बहुत ही उच्च सटीकता के लिए तेजी से परिणाम उत्पादन कर सकता का सुझाव दे this मिला, लेकिन यह एक भी आयाम के लिए है और मेरे लिए दो में परिवर्तित होना बहुत कठिन हो सकता है।

Answer 1

हाँ! क्लस्टरिंग के साथ के-मीन्स का उपयोग करना ठीक उसी तरह करेगा। तो s_ideal लक्ष्य होगा, फिर आप p_z पर ट्रेन करेंगे। आखिरकार आपको सेंट्रॉइड मिलेगा जो आपको p में प्रत्येक बिंदु पर सतह s_ideal पर निकटतम बिंदु देगा।

Here एक उदाहरण है, जो आप चाहते हैं उसके काफी करीब है।

Answer 2

प्रश्न तीन आयामी सतह S(x,y,z) और एक और बिंदु x0,y0,z0 के बीच यूक्लिडियन दूरी को कम करने का प्रयास करता है। सतह को आयताकार (x,y) जाल पर परिभाषित किया गया है, जहां z(x,y) = f(x,y) + random_noise(x,y)। "आदर्श" सतह पर शोर की शुरूआत समस्या के लिए काफी जटिलता जोड़ती है, क्योंकि इसे सतह को 2-आयामी तीसरे क्रम के स्पलीन का उपयोग करके इंटरपोलेट करने की आवश्यकता होती है।

यह समझा नहीं गया है कि आदर्श सतह पर शोर का परिचय वास्तव में क्यों आवश्यक है। यदि आदर्श सतह वास्तव में आदर्श थी, तो यह अच्छी तरह से समझना चाहिए कि x और y में एक वास्तविक बहुपद फिट निर्धारित किया जा सकता है, अगर विश्लेषणात्मक रूप से कम से कम अनुभवजन्य नहीं है। यदि यादृच्छिक शोर एक वास्तविक माप को अनुकरण करना था, तो केवल एक को माप को पर्याप्त समय तक रिकॉर्ड करने की आवश्यकता होती है जब तक कि शोर औसत से औसत न हो जाए। इसी तरह, सिग्नल फ़िल्टरिंग का उपयोग शोर को खत्म करने और सिग्नल के सच्चे व्यवहार को प्रकट करने में मदद कर सकता है।

सतह पर निकटतम बिंदु को दूसरी बिंदु पर खोजने के लिए, दूरी समीकरण और इसके डेरिवेटिव का उपयोग किया जाना चाहिए। अगर सतह को वास्तव में केवल splines के आधार का उपयोग करके वर्णित किया जा सकता है, तो एक को reconstruct स्पलीन प्रतिनिधित्व होना चाहिए और इसके डेरिवेटिव को ढूंढना चाहिए, जो गैर-तुच्छ है। वैकल्पिक रूप से, सतह को ठीक जाल का उपयोग करके मूल्यांकन किया जा सकता है, लेकिन यहां, एक मेमोरी के मुद्दों में जल्दी से चलता है, यही कारण है कि इंटरपोलेशन का इस्तेमाल पहली जगह में किया जाता था।

हालांकि, अगर हम इस बात से सहमत कर सकते हैं x और y में एक सरल अभिव्यक्ति का उपयोग करते हुए कि सतह परिभाषित किया जा सकता है, तो न्यूनतम तुच्छ हो जाता है:

न्यूनीकरण के प्रयोजनों के लिए, इसके बारे में वर्ग को देखने के लिए और अधिक सुविधाजनक है दूरी d^2(x,y) (zx और y का एक फ़ंक्शन है) दो बिंदुओं, D(x,y) के बीच, क्योंकि यह वर्ग रूट को समाप्त करता है। D(x,y) के महत्वपूर्ण बिंदुओं को ढूंढने के लिए, हम इसके आंशिक डेरिवेटिव w.r.t x और y लेते हैं और अपनी जड़ों को = 0: d/dx D(x,y) = f1(x,y) = 0 और d/dy D(x,y) = f2(x,y)=0 सेट करके ढूंढते हैं। यह समीकरणों की एक गैर-रैखिक प्रणाली है, जिसके लिए हम scipy.optimize.root का उपयोग करके हल कर सकते हैं। हमें केवल root अनुमान लगाने की आवश्यकता है (सतह पर ब्याज के pt का प्रक्षेपण) और समीकरणों की प्रणाली के Jacobian।

import numpy as np 
import scipy.interpolate 
import scipy.optimize 

# Define regular grid surface 
xmin,xmax,ymin,ymax = 25, 125, -50, 50 
x = np.linspace(xmin,xmax, 201) 
y = np.linspace(ymin,ymax, 201) 
xx, yy = np.meshgrid(x, y) 
z_ideal = (xx**2 + yy**2)/400 

# Fake some measured points on the surface 
z_measured = z_ideal + np.random.uniform(-0.1, 0.1, z_ideal.shape) 
s_measured = scipy.interpolate.interp2d(x, y, z_measured, kind='cubic') 
p_x = np.random.uniform(xmin,xmax,10000) 
p_y = np.random.uniform(ymin,ymax,10000) 

# z_ideal function 
def z(x): 
    return (x[0] ** 2 + x[1] ** 2)/400 

# returns the system of equations 
def f(x,pt): 
    x0,y0,z0 = pt 
    f1 = 2*(x[0] - x0) + (z(x)-z0)*x[0]/100 
    f2 = 2*(x[1] - y0) + (z(x)-z0)*x[1]/100 
    return [f1,f2] 

# returns Jacobian of the system of equations 
def jac(x, pt): 
    x0,y0,z0 = pt 
    return [[2*x[0]+1/100*(1/400*(z(x)+2*x[0]**2))-z0, x[0]*x[1]/2e4], 
    [2*x[1]+1/100*(1/400*(z(x)+2*x[1]**2))-z0, x[0]*x[1]/2e4]] 

def minimize_distance(pt): 
    guess = [pt[0],pt[1]] 
    return scipy.optimize.root(f,guess,jac=jac, args=pt) 

# select a random point from the measured data 
x0,y0 = p_x[30], p_y[30] 
z0 = float(s_measured(x0,y0)) 

minimize_distance([x0,y0,z0]) 

आउटपुट:

fjac: array([[-0.99419141, -0.1076264 ], 
     [ 0.1076264 , -0.99419141]]) 
    fun: array([ -1.05033229e-08, -2.63163477e-07]) 
message: 'The solution converged.' 
    nfev: 19 
    njev: 2 
    qtf: array([ 2.80642738e-07, 2.13792093e-06]) 
     r: array([-2.63044477, -0.48260582, -2.33011149]) 
    status: 1 
success: True 
     x: array([ 110.6726472 , 39.28642206])