Pythonでポイントとバイキュービックスプラインサーフェスの間の距離を素早く推定する方法はありますか？

Question

Pythonでポイントとバイキュービックスプラインサーフェスの間の距離を見積もるにはどうすればよいですか？ SciPy、NumPy、または他のパッケージで活用できる既存のソリューションはありますか？私は、最も近いを見つけたい

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

：

私はこのようバイキュービック補間によって定義された表面持っている：

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

をし、私はその表面の一部の測定点を持っています表面上の点s_ideal各点はpになります。一般的なケースでは、スプラインを大きく変化させるための複数の解法がある可能性があるので、zに沿った点の射影の近傍に1つの解のみを持つことが知られているサーフェスに限定しています。 これは計測点や表面の定義点が非常に少ないため、正確さを犠牲にしても速度を最適化したいと考えています（1E-5）。

頭に浮かぶ方法は、勾配降下アプローチを使用し、各測定点pのためのような何かを行うことです。

1. p_z = s_ideal(pt)
2. は（傾きを計算し、最初のテストポイントとして使用pt = [p_x, p_y, p_z]pt
3. における接線）ベクターm = [ m_x, m_y ]ptからpにベクターrを計算する：r = p - pt
4. rとmの間の角度thetaがある程度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 ] 

私は方法は、1Dケースのための非常に高い精度で高速な結果を生む可能性が示唆thisを見つけたが、それは、単一のディメンションのためだと私が2に変換するにはあまりにも難しいかもしれません。

Answer 1

はい！クラスタリングを使ってK-Meansを使うことはまさにそれを行うでしょう。だからs_idealがターゲットになり、p_zに列車に乗る。最終的には、重心を得て、pの各点に、表面上の最も近い点をs_idealに与えます。

Hereは、あなたが望むものにかなり近い例です。

Answer 2

質問は、3次元表面S(x,y,z)と別の点x0,y0,z0との間のユークリッド距離を最小化しようとしています。表面は、長方形の(x,y)メッシュに定義されており、z(x,y) = f(x,y) + random_noise(x,y)です。 「理想的な」サーフェスにノイズを導入すると、2次元の3次スプラインを使用してサーフェスを補間する必要があるため、問題がかなり複雑になります。

理想的な表面へのノイズの導入が実際に必要な理由はわかりません。理想的な表面が本当に理想的であれば、少なくとも経験的には解析的にではなく、xとyの真の多項式の適合を求めることができます。ランダムノイズが実際の測定値をシミュレートする場合、ノイズが平均化されてゼロになるまで測定を記録するだけで十分です。同様に、信号フィルタリングを使用すると、ノイズを除去して信号の真の動作を明らかにすることができます。

表面上の最も近い点を別の点に見つけるには、距離の式とその派生関数を使用する必要があります。スプラインの基底を使ってサーフェスを記述することが本当にしかできない場合は、reconstructスプライン表現を使用し、その導関数を見つける必要があります。これは自明ではありません。あるいは、細かいメッシュを使用してサーフェスを評価することもできますが、ここではすぐにメモリの問題が発生します。そのため、まず補間が使用されました。我々は、表面がxとyで簡単な式を使用して定義することができることに同意することができれば

はしかし、その後、最小化は些細次のようになります。

最小化の目的のためには、の広場を見て、より便利です距離d^2(x,y)（zは、xとyの関数に過ぎません。D(x,y)の2つの点の間に平方根が存在しません。臨界点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を使用して解くことができます。我々は、推測（ptの表面上への投影）と方程式系のJacobianをrootに渡すだけでよい。

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