使用 Scipy 在凸包上找到一个点的投影 [英] Find the projection of a point on the convex hull with Scipy

问题描述

从一组点来看，我得到了带有 scipy.spatial 的凸包，或者带有 Delaunay 或 ConvexHull(来自qhull 库).现在我想得到这个凸包外的一个点到船体上的投影(即船体上与外面的点距离最小的点).

这是我目前的代码:

from scipy.spatial import Delaunay, ConvexHull将 numpy 导入为 nphu = np.random.rand(10, 2) ## 得到船体的点集pt = np.array([1.1, 0.5]) ## 外面的一个点pt2 = np.array([0.4, 0.4]) ##里面的一个点hull = ConvexHull(hu) ## 只得到凸包#hull2 = Delaunay(hu) ## 或获取完整的 Delaunay 三角剖分导入 matplotlib.pyplot 作为 pltplt.plot(hu[:,0], hu[:,1], "ro") ## 绘制所有点#plt.triplot(hu[:,0], hu[:,1], hull2.simplices.copy()) ## 绘制 Delaunay 三角剖分## 绘制凸包对于 hull.simplices 中的单纯形:plt.plot(hu[simplex,0], hu[simplex,1], "ro-")## 绘制凸包内外的点plt.plot(pt[0], pt[1], "bs")plt.plot(pt2[0], pt2[1], "bs")plt.show()

使用图片可能更容易，我想从凸包外的蓝点获得绿色的 x 和 y 坐标.这个例子是二维的，但我也需要在更高的维度上应用它.感谢您的帮助.

这里解决了问题，但我在实现它时遇到了麻烦:

From a set of points, I'm getting the convex hull with scipy.spatial, either with Delaunay or ConvexHull (from the qhull library). Now I would like to get the projection of a point outside this convex hull onto the hull (i.e. the point on the hull that is the smallest distance from the point outside).

This is the code I have so far:

from scipy.spatial import Delaunay, ConvexHull
import numpy as np

hu = np.random.rand(10, 2) ## the set of points to get the hull from
pt = np.array([1.1, 0.5]) ## a point outside
pt2 = np.array([0.4, 0.4]) ## a point inside
hull = ConvexHull(hu) ## get only the convex hull
#hull2 = Delaunay(hu) ## or get the full Delaunay triangulation

import matplotlib.pyplot as plt
plt.plot(hu[:,0], hu[:,1], "ro") ## plot all points
#plt.triplot(hu[:,0], hu[:,1], hull2.simplices.copy()) ## plot the Delaunay triangulation
## Plot the convexhull
for simplex in hull.simplices:
    plt.plot(hu[simplex,0], hu[simplex,1], "ro-")

## Plot the points inside and outside the convex hull 
plt.plot(pt[0], pt[1], "bs")
plt.plot(pt2[0], pt2[1], "bs")
plt.show()

With a picture it might be easier, I would like to obtain the x and y coordinates in green from the blue point outside the convex hull. The example is 2d but I would need to apply it in higher dimension as well. Thanks for the help.

EDIT: The problem is addressed here but I have trouble implementing it: https://mathoverflow.net/questions/118088/projection-of-a-point-to-a-convex-hull-in-d-dimensions

解决方案

I am answering to myself. As 0Tech pointed out, ConvexHull.equations gives you the plane equations for each plane (in 2d --- a line therefore) with the form : [A, B, C]. The plane is therefore defined by

A*x + B*y + C = 0

Projecting the point P=(x0, y0) on the plane is explained here. You want a point on the vector parallel to the plane vector (A, B) and passing through the point to project P, this line is parameterised by t:

P_proj = (x, y) = (x0 + A*t, y0 + B*t)

You then want your point to be on the plane and uses the full plane equation to do that:

A*(x0 + A*t) + B*(y0 + B*t) + C = 0
=> t=-(C + A*x0 + B*y0)/(A**2+B**2)

In (clumsy) python, it gives for any dimension:

from scipy.spatial import Delaunay, ConvexHull
import numpy as np

hu = np.random.rand(10, 2) ## the set of points to get the hull from
pt = np.array([1.1, 0.5]) ## a point outside
pt2 = np.array([0.4, 0.4]) ## a point inside
hull = ConvexHull(hu) ## get only the convex hull
#hull2 = Delaunay(hu) ## or get the full Delaunay triangulation

import matplotlib.pyplot as plt
plt.plot(hu[:,0], hu[:,1], "ro") ## plot all points
#plt.triplot(hu[:,0], hu[:,1], hull2.simplices.copy()) ## plot the Delaunay triangulation
## Plot the convexhull
for simplex in hull.simplices:
    plt.plot(hu[simplex,0], hu[simplex,1], "ro-")

## Plot the points inside and outside the convex hull 
plt.plot(pt[0], pt[1], "bs")
plt.plot(pt2[0], pt2[1], "bs")

for eq in hull.equations:
    t = -(eq[-1] + np.dot(eq[:-1], pt))/(np.sum(eq[:-1]**2))
    pt_proj = pt + eq[:-1]*t
    plt.plot(pt_proj[0], pt_proj[1], "gD-.")
plt.show()

---

Browsing stackoverflow, led me to another solution, that has the advantage of using segments instead of the lines, so the projection on one of the segment always lie on the segment:

def min_distance(pt1, pt2, p):
    """ return the projection of point p (and the distance) on the closest edge formed by the two points pt1 and pt2"""
    l = np.sum((pt2-pt1)**2) ## compute the squared distance between the 2 vertices
    t = np.max([0., np.min([1., np.dot(p-pt1, pt2-pt1) /l])]) # I let the answer of question 849211 explains this
    proj = pt1 + t*(pt2-pt1) ## project the point
    return proj, np.sum((proj-p)**2) ## return the projection and the point

Then we can browse each vertices and project the point:

for i in range(len(hull.vertices)):
    pt_proj, d = min_distance(hu[hull.vertices[i]], hu[hull.vertices[(i+1)%len(hull.vertices)]], pt)
    plt.plot([pt[0], pt_proj[0]], [pt[1], pt_proj[1]], "c<:")

And the picture, with the projection of the blue point on the right on each plane (line) in green for the first method and cyan for the second method:

这篇关于使用 Scipy 在凸包上找到一个点的投影的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！