加速双曲抛物面算法上的最近点 [英] Speeding up a closest point on a hyperbolic paraboloid algorithm
问题描述
我编写了一个 python 脚本，它从查询点 (p) 找到表面上最近点的 UV 坐标.该表面由由逆时针列出的四个已知点 (p0,p1,p2,p3) 组成的四个线性边定义.
(请忽略小红球)
我的方法的问题是它非常慢(以低精度阈值执行 5000 次查询大约需要 10 秒.
我正在寻找一种更好的方法来实现我想要的目标，或者是让我的代码更高效的建议.我唯一的限制是它必须用 python 编写.
将 numpy 导入为 np# 定义常量LARGE_VALUE=99999999.0SMALL_VALUE=0.00000001子样本=10.0def closePointOnLineSegment(a,b,c):''' 给定定义线段的两点 (a,b) 和查询点 (c)返回该段上最近的点，之间的距离查询点和最近点，以及从结果中得出的 u 值'''# 检查 c 是否与 a 或 b 相同ac=c-aAC=np.linalg.norm(ac)如果 AC==0.:返回 c,0.,0.bc=c-bBC=np.linalg.norm(bc)如果 BC==0.:返回 c,0.,1.# 查看段长是否为0ab=b-aAB=np.linalg.norm(ab)如果 AB == 0.:返回 a,0.,0.# 标准化段并进行边缘测试ab=ab/AB测试= np.dot(ac,ab)如果测试<0.:返回 a,AC,0.elif 测试 >AB:返回 b,BC,1.# 返回最近的 xyz 段、距离和 u 值p=(测试*ab)+a返回 p,np.linalg.norm(c-p),(test/AB)def surfaceWalk(e0,e1,p,v0=0.,v1=1.):''' 沿着 2 条边走在表面上，对于每个样本段寻找最近点，递归直到两个采样边缘小于 SMALL_VALUE'''edge0=(e1[0]-e0[0])edge1=(e1[1]-e0[1])len0=np.linalg.norm(edge0*(v1-v0))len1=np.linalg.norm(edge1*(v1-v0))vMin=v0vMax=v1最近的_d=0.最近_u=0.最近_v=0.ii=0.dist=LARGE_VALUE对于范围内的 i(int(SUBSAMPLES)+1):v=v0+((v1-v0)*(i/SUBSAMPLES))a=(edge0*v)+e0[0]b=(edge1*v)+e0[1]closest_p,closest_d,closest_u=closestPointOnLineSegment(a,b,p)如果最近的_d <距离:距离=最近的_d最近_v=vii=i# 如果两个边长都<= SMALL_VALUE，我们在我们的精度阈值内，所以返回结果如果 len0 <= SMALL_VALUE 和 len1 <= SMALL_VALUE:返回最接近的_p，最接近的_d，最接近的_u，最接近的_v# 未达到阈值，将 v0 anf v1 限制设置为最近的_v 的任一侧并继续递归vMin=v0+((v1-v0)*(np.clip((ii-1),0.,SUBSAMPLES)/SUBSAMPLES))vMax=v0+((v1-v0)*(np.clip((ii+1),0.,SUBSAMPLES)/SUBSAMPLES))返回surfaceWalk(e0,e1,p,vMin,vMax)def closePointToPlane(p0,p1,p2,p3,p,debug=True):''' 给定定义四边形曲面的四个点 (p0,p1,p2,3) 和一个查询点 p.找到最近的边缘并开始行走穿过一端到另一端，直到我们找到最近的点'''# 找到最近的边，我们将使用该边开始我们的步行c,d,u,v=surfaceWalk([p0,p1],[p3,p2],p)如果调试:打印最近点:%s"%c打印到点的距离:%s"%d打印 'U 坐标:%s'%u打印 'V 坐标:%s'%v返回 c,d,u,vp0 = np.array([1.15, 0.62, -1.01])p1 = np.array([1.74, 0.86, -0.88])p2 = np.array([1.79, 0.40, -1.46])p3 = np.array([0.91, 0.79, -1.84])p = np.array([1.17, 0.94, -1.52])最近点到平面(p0，p1，p2，p3，p)最近点:[ 1.11588876 0.70474519 -1.52660706]到点的距离:0.241488104197U 坐标:0.164463481066V 坐标:0.681959858995如果表面看起来是双曲抛物面，则可以将其上的点 s 参数化为:
s = p0 + u * (p1 - p0) + v * (p3 - p0) + u * v * (p2 - p3 - p1 + p0)这样做，p0p3 行有等式 u = 0，p1p2 是 u = 1, p0p1 是 v = 0 而 p2p3 是 v = 1.我还没有想出一种方法来为最接近点 p 的点提出解析表达式，但是 scipy.optimize 可以用数字来完成这项工作给你:
将 numpy 导入为 np从 scipy.optimize 导入最小化p0 = np.array([1.15, 0.62, -1.01])p1 = np.array([1.74, 0.86, -0.88])p2 = np.array([1.79, 0.40, -1.46])p3 = np.array([0.91, 0.79, -1.84])p = np.array([1.17, 0.94, -1.52])def fun(x, p0, p1, p2, p3, p):你，v = xs = u*(p1-p0) + v*(p3-p0) + u*v*(p2-p3-p1+p0) + p0返回 np.linalg.norm(p - s)>>>最小化(有趣，(0.5，0.5)，(p0，p1，p2，p3，p))状态:0成功:正确涅夫:9非传染性疾病:36乐趣:0.24148810420527048x: 数组([ 0.16446403, 0.68196253])消息:'优化成功终止.'赫斯:数组([[ 0.38032445, 0.15919791],[ 0.15919791, 0.44908365]])jac:数组([-1.27032399e-06，6.74091280e-06])minimize 的返回是一个对象，你可以通过它的属性访问值:
<小时>
关于如何在没有 scipy 的情况下找到解决方案的一些指示......连接点 s 的向量与通用点 p 是ps.要找出最近的点，您可以采用两种不同的方法来得到相同的结果:
- 计算该向量的长度，
(p-s)**2，取其导数 w.r.t.u和v并将它们等同于零. - 计算与
s处的 hypar 相切的两个向量，可以找到ds/du和ds/dv，并强加它们与ps的内积为零.
如果你解决了这些问题，你最终会得到两个方程，这些方程需要大量的操作才能得出 u 或 v 的三次或四次方程，所以没有精确的解析解，尽管你可以只用 numpy 数值求解.一个更简单的选择是计算这些方程，直到得到这两个方程，其中 a = p1-p0, b = p3-p0, c = p2-p3-p1+p0, s_ = s-p0, p_ = p-p0:
u = (p_ - v*b)*(a + v*c)/(a + v*c)**2v = (p_ - u*a)*(b + u*c)/(b + u*c)**2你不能轻易地想出一个解析解，但你可以希望如果你使用这两个关系来迭代一个试验解，它会收敛.对于您的测试用例，它确实有效:
def solve_uv(x0, p0, p1, p2, p3, p, tol=1e-6, niter=100):a = p1 - p0b = p3 - p0c = p2 - p3 - p1 + p0p_ = p - p0u, v = x0错误 = 无while niter and (error is None or error > tol):硝基 -= 1u_ = np.dot(p_ - v*b, a + v*c)/np.dot(a + v*c, a + v*c)v_ = np.dot(p_ - u*a, b + u*c)/np.dot(b + u*c, b + u*c)错误 = np.linalg.norm([u - u_, v - v_])u, v = u_, v_返回 np.array([u, v]), np.linalg.norm(u*a + v*b +u*v*c + p0 - p)>>>solve_uv([0.5, 0.5], p0, p1, p2, p3, p)(数组([0.16446338，0.68195998])，0.2414881041967159)我认为这不能保证收敛，但对于这种特殊情况，它似乎工作得很好，只需要 15 次迭代即可达到要求的容差.
I wrote a python script which finds the UV coords of the closest point on surface from a query point (p). The surface is defined by four linear edges made from four known points (p0,p1,p2,p3) listed counter clockwise.
(Please ignore the little red ball)
The problem with my approach is that it is very slow (~10 seconds to do 5000 queries with a low precision threshold.
I'm looking for a better approach to achieve what i want, or suggestions to make my code more efficient. My only constraint is that it must be written in python.
import numpy as np
# Define constants
LARGE_VALUE=99999999.0
SMALL_VALUE=0.00000001
SUBSAMPLES=10.0
def closestPointOnLineSegment(a,b,c):
''' Given two points (a,b) defining a line segment and a query point (c)
return the closest point on that segment, the distance between
query and closest points, and a u value derived from the results
'''
# Check if c is same as a or b
ac=c-a
AC=np.linalg.norm(ac)
if AC==0.:
return c,0.,0.
bc=c-b
BC=np.linalg.norm(bc)
if BC==0.:
return c,0.,1.
# See if segment length is 0
ab=b-a
AB=np.linalg.norm(ab)
if AB == 0.:
return a,0.,0.
# Normalize segment and do edge tests
ab=ab/AB
test=np.dot(ac,ab)
if test < 0.:
return a,AC,0.
elif test > AB:
return b,BC,1.
# Return closest xyz on segment, distance, and u value
p=(test*ab)+a
return p,np.linalg.norm(c-p),(test/AB)
def surfaceWalk(e0,e1,p,v0=0.,v1=1.):
''' Walk on the surface along 2 edges, for each sample segment
look for closest point, recurse until the both sampled edges
are smaller than SMALL_VALUE
'''
edge0=(e1[0]-e0[0])
edge1=(e1[1]-e0[1])
len0=np.linalg.norm(edge0*(v1-v0))
len1=np.linalg.norm(edge1*(v1-v0))
vMin=v0
vMax=v1
closest_d=0.
closest_u=0.
closest_v=0.
ii=0.
dist=LARGE_VALUE
for i in range(int(SUBSAMPLES)+1):
v=v0+((v1-v0)*(i/SUBSAMPLES))
a=(edge0*v)+e0[0]
b=(edge1*v)+e0[1]
closest_p,closest_d,closest_u=closestPointOnLineSegment(a,b,p)
if closest_d < dist:
dist=closest_d
closest_v=v
ii=i
# If both edge lengths <= SMALL_VALUE, we're within our precision treshold so return results
if len0 <= SMALL_VALUE and len1 <= SMALL_VALUE:
return closest_p,closest_d,closest_u,closest_v
# Threshold hasn't been met, set v0 anf v1 limits to either side of closest_v and keep recursing
vMin=v0+((v1-v0)*(np.clip((ii-1),0.,SUBSAMPLES)/SUBSAMPLES))
vMax=v0+((v1-v0)*(np.clip((ii+1),0.,SUBSAMPLES)/SUBSAMPLES))
return surfaceWalk(e0,e1,p,vMin,vMax)
def closestPointToPlane(p0,p1,p2,p3,p,debug=True):
''' Given four points defining a quad surface (p0,p1,p2,3) and
a query point p. Find the closest edge and begin walking
across one end to the next until we find the closest point
'''
# Find the closest edge, we'll use that edge to start our walk
c,d,u,v=surfaceWalk([p0,p1],[p3,p2],p)
if debug:
print 'Closest Point: %s'%c
print 'Distance to Point: %s'%d
print 'U Coord: %s'%u
print 'V Coord: %s'%v
return c,d,u,v
p0 = np.array([1.15, 0.62, -1.01])
p1 = np.array([1.74, 0.86, -0.88])
p2 = np.array([1.79, 0.40, -1.46])
p3 = np.array([0.91, 0.79, -1.84])
p = np.array([1.17, 0.94, -1.52])
closestPointToPlane(p0,p1,p2,p3,p)
Closest Point: [ 1.11588876 0.70474519 -1.52660706]
Distance to Point: 0.241488104197
U Coord: 0.164463481066
V Coord: 0.681959858995
If your surface is, as it seems, a hyperbolic paraboloid, you can parametrize a point s on it as:
s = p0 + u * (p1 - p0) + v * (p3 - p0) + u * v * (p2 - p3 - p1 + p0)
Doing things this way, the line p0p3 has equation u = 0, p1p2 is u = 1, p0p1 is v = 0 and p2p3 is v = 1. I haven't been able to figure out a way of coming up with an analytical expression for the closest point to a point p, but scipy.optimize can do the job numerically for you:
import numpy as np
from scipy.optimize import minimize
p0 = np.array([1.15, 0.62, -1.01])
p1 = np.array([1.74, 0.86, -0.88])
p2 = np.array([1.79, 0.40, -1.46])
p3 = np.array([0.91, 0.79, -1.84])
p = np.array([1.17, 0.94, -1.52])
def fun(x, p0, p1, p2, p3, p):
u, v = x
s = u*(p1-p0) + v*(p3-p0) + u*v*(p2-p3-p1+p0) + p0
return np.linalg.norm(p - s)
>>> minimize(fun, (0.5, 0.5), (p0, p1, p2, p3, p))
status: 0
success: True
njev: 9
nfev: 36
fun: 0.24148810420527048
x: array([ 0.16446403, 0.68196253])
message: 'Optimization terminated successfully.'
hess: array([[ 0.38032445, 0.15919791],
[ 0.15919791, 0.44908365]])
jac: array([ -1.27032399e-06, 6.74091280e-06])
The return of minimize is an object, you can access the values through its attributes:
>>> res = minimize(fun, (0.5, 0.5), (p0, p1, p2, p3, p))
>>> res.x # u and v coordinates of the nearest point
array([ 0.16446403, 0.68196253])
>>> res.fun # minimum distance
0.24148810420527048
Some pointers on how to go about finding a solution without scipy... The vector joining the point s parametrized as above with a generic point p is p-s. TO find out the closest point you can go two different ways that give the same result:
- Compute the length of that vector,
(p-s)**2, take its derivatives w.r.t.uandvand equate them to zero. - Compute two vectors tangent to the hypar at
s, which can be found asds/duandds/dv, and impose that their inner product withp-sbe zero.
If you work these out, you'll end up with two equations that would require a lot of manipulation to arrive at something like a third or fourth degree equation for either u or v, so there is no exact analytical solution, although you could solve that numerically with numpy only. An easier option is to work out those equations until you get these two equations, where a = p1-p0, b = p3-p0, c = p2-p3-p1+p0, s_ = s-p0, p_ = p-p0:
u = (p_ - v*b)*(a + v*c) / (a + v*c)**2
v = (p_ - u*a)*(b + u*c) / (b + u*c)**2
You can't come up with an analytical solution for this easily, but you can hope that if you use those two relations to iterate a trial solution, it will converge. For your test case it does work:
def solve_uv(x0, p0, p1, p2, p3, p, tol=1e-6, niter=100):
a = p1 - p0
b = p3 - p0
c = p2 - p3 - p1 + p0
p_ = p - p0
u, v = x0
error = None
while niter and (error is None or error > tol):
niter -= 1
u_ = np.dot(p_ - v*b, a + v*c) / np.dot(a + v*c, a + v*c)
v_ = np.dot(p_ - u*a, b + u*c) / np.dot(b + u*c, b + u*c)
error = np.linalg.norm([u - u_, v - v_])
u, v = u_, v_
return np.array([u, v]), np.linalg.norm(u*a + v*b +u*v*c + p0 - p)
>>> solve_uv([0.5, 0.5], p0, p1, p2, p3, p)
(array([ 0.16446338, 0.68195998]), 0.2414881041967159)
I don't think this is guaranteed to converge, but for this particular case it seems to work pretty fine, and only needs 15 iterations to get to the requested tolerance.
这篇关于加速双曲抛物面算法上的最近点的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
