I want to compute the distance between a point in 3D and a convex set in 3D which is given as the intersection of halfplanes in numpy (python).

The distance between a point and a set is defined as the infimum (minimum) of the distance from the point and any point in the set.

Suppose the point is given as a 1D numpy array eg:

P = np.array([-1,-2,-3], ndmin = 1)

P represents the point (x,y,z) = (-1,-2,-3). And suppose the convex set is given as a 2D numpy array eg:

C = np.array([1,2,3,4], [5,6,7,8] , [9,10,11,12] )

C represents the convex set given by (x 2y 3y >= 4) /\ (5x 6y 7y >= 8) /\ (9x 10y 11y >= 12). How can I compute the distance between P and C? (We can assume that all halfplanes only use the operator >=)

Note: As mentioned by Pierre in the comments, the distance of the point is not just the minimum of the distance of the point and all line segments on the boundary of the halfplanes. We have to consider the faces of the convex set too.

CodePudding user response：

A simple way could be to consider the convex optimization (in this case a quadratic program) of finding the closest point to P, s.t. constraints that Ax b <= 0. Note that there might be more efficient, specialized and accurate ways of doing this (see. e.g. this mathoverflow answer and the links it refers to).

You can solve a convex problem using scipy.optimize.minimize. Your half-spaces define the constraints, and the function to minimize is the distance (or its square) to the point P.

Note that we use the convention Ax b <= 0 for the half-spaces, and encode them [A; b].

import numpy as np
from scipy.optimize import minimize, LinearConstraint

def dist(x, p):
    return np.linalg.norm(x - p)

def find_closest(halfspaces, p):
    return minimize(
        dist, 
        np.zeros(3),
        args=(p,),
        constraints=[LinearConstraint(halfspaces[:, :-1], -np.inf, -halfspaces[:, -1])],
    )

As an example, consider the tetrahedron defined by the points (0,0,0), (2,0,0), (1,2,0), and (1,1,1):

halfspaces = np.array(
    [[0,0,-1,0],  # Ax   b <= 0; each row of A is a vector normal to 
                  # hyperplane and pointing outside C
     [2,1,1,-4],
     [-2,1,1,0],
     [0,-1,1,0]])

Then:

>>> find_closest(halfspaces, [0,0,10])
     fun: 9.110433579144301
     jac: array([ 0.10976422,  0.10976422, -0.98787833])
 message: 'Optimization terminated successfully'
    nfev: 12
     nit: 3
    njev: 3
  status: 0
 success: True
       x: array([1., 1., 1.])

# note: closest to a point (the top point)


>>> find_closest(halfspaces, [1,-0.1,0])
     fun: 0.10000032194980375
     jac: array([2.53758579e-03, 9.99996780e-01, 7.45058060e-08])
 message: 'Optimization terminated successfully'
    nfev: 13
     nit: 3
    njev: 3
  status: 0
 success: True
       x: array([ 1.00025375e 00, -2.22044605e-17,  6.04908606e-17])

# note: closest to an edge (the x-aligned line segment)


>>> find_closest(halfspaces, [1,1,-0.5])
     fun: 0.5000000000000011
     jac: array([0.00000000e 00, 5.96046448e-08, 1.00000000e 00])
 message: 'Optimization terminated successfully'
    nfev: 26
     nit: 6
    njev: 6
  status: 0
 success: True
       x: array([9.99999991e-01, 1.00000002e 00, 4.83315889e-16])

# note: closest to a face (the bottom face)