Làm cách nào để có được ước tính nhanh về khoảng cách giữa một điểm và bề mặt đường tròn bicubic trong Python?

Question

Làm cách nào để có được ước tính nhanh về khoảng cách giữa một điểm và bề mặt đường tròn bicubic bằng Python? Có một giải pháp hiện có mà tôi có thể tận dụng trong SciPy, NumPy, hoặc một số gói khác?

Tôi đã có mặt xác định bởi một suy bicubic như thế này:

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

và tôi đã có một số điểm đo của bề mặt:

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

Tôi muốn tìm gần nhất chỉ trên bề mặt s_ideal đến từng điểm trong p. Một trường hợp tổng quát có thể có nhiều giải pháp cho splines cực kỳ khác nhau, vì vậy tôi giới hạn vấn đề với các bề mặt được biết là chỉ có một giải pháp trong vùng lân cận của phép chiếu của điểm dọc theo z. Đây không phải là một số lượng nhỏ các điểm đo lường hoặc độ nét bề mặt, vì vậy tôi muốn tối ưu hóa tốc độ ngay cả với chi phí chính xác có thể là 1E-5.

Phương pháp mà nói đến cái tâm là sử dụng một phương pháp gradient descent và làm điều gì đó giống như cho mỗi điểm đo p:

1. Sử dụng pt = [p_x, p_y, p_z] như là điểm thử nghiệm ban đầu, nơi p_z = s_ideal(pt)
2. Tính dốc (tangent) vector m = [ m_x, m_y ] tại pt
3. Tính vector rpt-p: r = p - pt
012.
4. Nếu góc theta giữa r và m nằm trong ngưỡng tối đa là 90 độ, thì pt là điểm cuối cùng.
5. Nếu không, cập nhật pt như:

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 ] 

tôi thấy this gợi ý một phương pháp có thể tạo ra kết quả nhanh chóng đến một độ chính xác rất cao đối với trường hợp 1D, nhưng nó cho một chiều hướng duy nhất và có thể là quá khó cho tôi để chuyển đổi thành hai.

Answer 1

Có! Sử dụng K-Means với clustering sẽ làm chính xác điều đó. Vì vậy, s_ideal sẽ là mục tiêu, sau đó bạn đào tạo trên p_z. Cuối cùng, bạn sẽ nhận được các centroids sẽ cho bạn điểm gần nhất trên bề mặt s_ideal đến từng điểm trong p.

Here là một ví dụ, nó khá gần với những gì bạn muốn.

Answer 2

Câu hỏi tìm cách giảm thiểu khoảng cách Euclidian giữa bề mặt ba chiều S(x,y,z) và một điểm khác x0,y0,z0. Bề mặt được xác định trên lưới (x,y) hình chữ nhật, trong đó z(x,y) = f(x,y) + random_noise(x,y). Việc đưa tiếng ồn vào bề mặt "lý tưởng" làm tăng thêm sự phức tạp đáng kể cho vấn đề, vì nó đòi hỏi bề mặt phải được nội suy bằng cách sử dụng một spline thứ ba 2 chiều.

Không hiểu tại sao việc đưa tiếng ồn vào bề mặt lý tưởng thực sự là cần thiết. Nếu bề mặt lý tưởng thực sự lý tưởng, nó phải hiểu rõ rằng một sự đa thức thực sự trong x và y có thể được xác định, nếu không phân tích, ít nhất là theo kinh nghiệm. Nếu tiếng ồn ngẫu nhiên là để mô phỏng một phép đo thực tế, thì người ta chỉ cần ghi lại số đo đủ cho đến khi nhiễu được tính trung bình bằng không. Tương tự, việc sử dụng lọc tín hiệu có thể giúp loại bỏ nhiễu và tiết lộ hành vi thực sự của tín hiệu.

Để tìm điểm gần nhất trên bề mặt đến điểm khác, phương trình khoảng cách và các dẫn xuất của nó phải được sử dụng. Nếu bề mặt thực sự chỉ có thể được mô tả bằng cách sử dụng cơ sở của splines, sau đó một phải reconstruct đại diện spline và tìm thấy các dẫn xuất của nó, đó là không tầm thường. Ngoài ra, bề mặt có thể được đánh giá bằng cách sử dụng một lưới tốt, nhưng ở đây, một trong những nhanh chóng chạy vào các vấn đề bộ nhớ, đó là lý do tại sao nội suy được sử dụng ở nơi đầu tiên.

Tuy nhiên, nếu chúng ta có thể đồng ý rằng bề mặt có thể được xác định bằng cụm đơn giản trong x và y, sau đó giảm thiểu trở nên tầm thường:

Đối với mục đích của giảm thiểu, nó là thuận tiện hơn để nhìn vào vuông khoảng cách d^2(x,y) (z chỉ là một chức năng của x và y) giữa hai điểm, D(x,y), vì nó loại bỏ căn bậc hai. Để tìm các điểm quan trọng của D(x,y), chúng tôi lấy các dẫn xuất một phần của nó w.r.t x và y và tìm nguồn gốc của chúng bằng cách đặt = 0: d/dx D(x,y) = f1(x,y) = 0 và d/dy D(x,y) = f2(x,y)=0. Đây là hệ phương trình phi tuyến tính, mà chúng ta có thể giải quyết bằng cách sử dụng scipy.optimize.root. Chúng ta chỉ cần vượt qua root một dự đoán (chiếu của pt quan tâm lên bề mặt) và Jacobian của hệ phương trình.

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

Output:

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