加速双曲抛物面算法的最近点

Question

我写了一个python脚本,它从查询点(p)中找到表面上最近点的UV坐标.表面由四个线性边缘限定,这四个线性边缘由逆时针列出的四个已知点(p0,p1,p2,p3)组成.

(请忽略小红球)

我的方法的问题是它非常慢(使用低精度阈值进行5000次查询约10秒).

我正在寻找一种更好的方法来实现我想要的,或者建议使我的代码更有效率.我唯一的限制是它必须用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

Answer 1

如果您的表面看起来像是一个双曲抛物面,您可以将其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可以用数字方式完成这项工作:

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

返回minimize是一个对象,您可以通过其属性访问值:

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

关于如何在没有scipy的情况下寻找解决方案的一些指示...将s如上所述参数化的点与一般点相连接的矢量p是p-s.要找出最接近的点,您可以采用两种不同的方式来得到相同的结果:

1. 计算该向量的长度(p-s)**2,取其导数wrt u并将v它们等于零.
2. 计算与hypar相切的两个向量s,可以找到ds/du和ds/dv,并强制它们的内积p-s为零.

如果这些工作了,你会拥有两个方程,将需要大量的操作,以在像第三或第四度方程要么到达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)**2
v = (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 - 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)

我不认为这可以保证收敛,但对于这种特殊情况,它似乎工作得非常好,并且只需要15次迭代才能达到所要求的容差.