Find the projection of a point on the convex hull with Scipy

Question

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

Answer 1

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:

Answer 2

Since there still doesn't seem to be a good answer here (or anywhere), I followed this post (among others) and solved the problem using quadratic programming:

$$ \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} $$

Where $C$ are the normal equations and $b$ is the offsets. The key difference here is that I am projecting points to within the convex hull, not necessarily onto it. i.e. a point within the convex hull will remain unchanged and points outside the convex hull will be projected to the closest point in the convex hull (which will always be on the surface of it). By removing this equality constraint, it's a relatively simple quadratic programming problem which I solve using the quadprog package.

There might be a theoretically faster way to do it but this way is fast enough, simple, and robust:

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

A simple example:

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

(Edit: Sorry that it doesn't look like the Latex is formatting)