使用Scipy找到点在凸包上的投影

Question

我正在使用scipy.spatial中的Delaunay或ConvexHull（来自qhull库）从一组点中获取凸包。现在，我想获取一个在凸包外部的点到凸包上的投影（即来自点外最小距离的凸包上的点）。

这是我目前的代码：

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

有一张图片可能会更容易，我想获取凸包外的蓝点中绿色坐标的x和y值。这个例子是2D的，但我也需要在更高的维度中应用它。 感谢您的帮助。

编辑：这个问题在这里得到了解决，但我在实现时遇到了麻烦：https://mathoverflow.net/questions/118088/projection-of-a-point-to-a-convex-hull-in-d-dimensions

- Jean

---

您可以从ConvexHull.equations获取凸包的超平面，因此您可以计算点与每个超平面之间的距离。然后，您可以使用这些距离找到最近的超平面。但是，您应该按照自己的方式_定义_坐标，因为您只有法向量。 - 0Tech

Answer 1

我在回答自己的问题。 正如0Tech所指出的那样，ConvexHull.equations会为每个平面（在2D中是一条直线）提供平面方程，形式为：[A, B, C]。因此，该平面由以下公式定义：

A*x + B*y + C = 0

在平面上投影点P =（x0，y0）的方法在这里解释。您希望在与平面向量（A，B）平行且通过要投影的点P的向量上找到一个点，该线由t参数化：

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

你需要使你的点位于平面上，并使用完整的平面方程来实现这一点：

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

在（笨拙的）Python 中，它适用于任何维度：

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

---

在浏览stackoverflow时，我发现了另一种解决方案，它使用线段而不是直线，因此在一个线段上的投影总是位于该线段上，具有这样的优势：

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

然后，我们可以遍历每个顶点并投影该点：

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<:")

以下是第一种方法中每个平面（线）上蓝色点的投影以绿色显示，第二种方法中以青色显示：

- Jean

---

第一种解决方案适用于任何维度，但是投影不一定在凸包上。 - Jean

请查看答案http://math.stackexchange.com/questions/2146961/find-a-point-on-the-convex-hull-of-a-given-set-of-points-which-is-closest-to-a-g/2153051#2153051，其中提供了使用带约束的二次优化求解的方法，但对于高维空间可能太慢。 - Jean

Answer 2

因为似乎这里（或任何地方）还没有一个好的答案，所以我遵循了这篇帖子（和其他一些帖子），使用二次规划解决了这个问题：

Original Answer翻译成“最初的回答”。

$$ \begin{align}\text{minimize} & \quad \frac{1}{2} x^{T} \mathbf{I} x - z^{T} x\\ \text{subject to} & \quad C x \leq b\\ \end{align} $$

其中$C$是法线方程，$b$是偏移量。关键区别在于我将点投影到凸包内部，而不一定是投影到凸包表面。也就是说，在凸包内部的点将保持不变，而在凸包外部的点将被投影到凸包中最靠近它的点上（该点总是在凸包表面上）。通过去除这个等式约束，这是一个相对简单的二次规划问题，我使用quadprog软件包来解决。

可能有一种理论上更快的方法来实现，但这种方法足够快速、简单和稳健。

import numpy as np
from scipy.spatial import ConvexHull
from quadprog import solve_qp

def proj2hull(z, equations):
    """
    Project `z` to the convex hull defined by the
    hyperplane equations of the facets

    Arguments
        z: array, shape (ndim,)
        equations: array shape (nfacets, ndim + 1)

    Returns
        x: array, shape (ndim,)
    """
    G = np.eye(len(z), dtype=float)
    a = np.array(z, dtype=float)
    C = np.array(-equations[:, :-1], dtype=float)
    b = np.array(equations[:, -1], dtype=float)
    x, f, xu, itr, lag, act = solve_qp(G, a, C.T, b, meq=0, factorized=True)
    return x

一个简单的例子：

一个简单的例子：

X = np.random.normal(size=(1000, 5))
z = np.random.normal(scale=2, size=(5))
hull = ConvexHull(X)
y = proj2hull(z, hull.equations)

最初的回答：抱歉，Latex格式似乎没有被正确显示。

- stevens_mom