线性代数方法（我不认为这里需要LP）。使用一个矩阵乘法进行超平面half-space测试，然后进行一些max（）和sign（）后处理。
您可以通过在任何半空间测试之前执行直线修剪来变得更聪明，并在任何一个半空间测试失败时将矩阵乘法分区并停止。当一些重要部分的点在单纯形之外测试时，它帮助最大。在最极端的情况下，如果单纯形中不包含任何点（尝试半径=1.1），则非分区算法需要约0.6秒，而50个分区需要约0.01秒。

from time import monotonic

import numpy as np
from numpy.random import default_rng
from scipy.spatial import ConvexHull

def make_homogeneous(test_points: np.ndarray) -> np.ndarray:
    """
    Pre-process an array of (p, ndim) test points into a homogeneous
    transformation matrix of size (ndim+1, p). This only needs to be
    done once for a given point cloud.
    """
    return np.vstack((
        test_points.T,
        np.ones(shape=test_points.shape[0], dtype=test_points.dtype),
    ))

def test_hull(
    hull: ConvexHull, test_homogeneous: np.ndarray,
    n_partitions: int = 1, trim: bool = True,
) -> np.ndarray:
    """
    Vectorized test of whether each test point falls within the simplex.

    :param hull: Hull defining the simplex. Number of dimensions (i.e. hull.equations.shape[1]-1)
                 must be equal to number of dimensions in the test point cloud.
    :param test_homogeneous:
                 Test point cloud, in homogeneous transformation matrix format (from
                 make_homoegeneous()).
    :param n_partitions:
                 Number of inner product partitions. If the number of points falling inside of the
                 simplex is high, partitioning will not help and this should be left as 1 (non-
                 partitioned). If the number of points falling inside of the simplex is low, set
                 this on the order of ~ 1% of the number of hull faces.
    :param trim: Whether to perform a rectilinear trim before any dot products.
    :return: A boolean array with length as the number of test points: true for "inside", false for
             "outside". Values exactly on the simplex boundary are treated as "true" (inside the
             simplex) due to `<= 0` below.
    """

    # m: number of hull faces (to be partitioned)
    # n: 1 + ndim
    # p: number of test points
    m, n = hull.equations.shape
    n, p = test_homogeneous.shape
    partition_size = m // n_partitions

    if trim:
        extents0 = hull.points.min(axis=0)[:, np.newaxis]
        extents1 = hull.points.max(axis=0)[:, np.newaxis]
        inside = (
            (test_homogeneous[:3, :] >= extents0).all(axis=0) &
            (test_homogeneous[:3, :] <= extents1).all(axis=0)
        )
        test_subset = test_homogeneous[:, inside]
        # print(f'Trimmed to {np.count_nonzero(inside)}/{p} points')
    else:
        inside = np.ones(shape=p, dtype=bool)
        test_subset = test_homogeneous

    for i in range(0, m, partition_size):
        hull_subset = hull.equations[i: i + partition_size, :]
        product = hull_subset @ test_subset

        inside_subset = product.max(axis=0) <= 0
        # inside_subset = (product < 0).all(axis=0)  # Equivalent, marginally slower?

        inside[inside] = inside_subset
        if not inside_subset.any():
            break

        test_subset = test_subset[:, inside_subset]

    return inside

def cube_test() -> None:
    """
    Unit test for a cube-shaped hull (2 hull facets per cube side, for 12 facets total)
    """
    hull = ConvexHull([
        [-1, -1, -1],
        [ 1, -1, -1],
        [-1,  1, -1],
        [ 1,  1, -1],
        [-1, -1,  1],
        [ 1, -1,  1],
        [-1,  1,  1],
        [ 1,  1,  1],
    ])
    in_points = np.array([[ 0. ,  0. ,  0. ],
                          [ 0.9,  0. ,  0.2],
                          [-0.5, -0.2, -0.3],
                          [-0.9,  0.4,  0.6],
                          [ 0.1,  0.1,  0.1]])
    bound_points = np.array([[ 1. ,  1. ,  1. ],
                             [-1. , -1. , -1. ],
                             [ 0.5,  0. ,  1. ],
                             [ 1. ,  1. ,  0. ],
                             [ 0. ,  0. ,  1. ]])
    out_points = np.array([[ 2. ,  0. ,  0. ],
                           [ 1. ,  1. ,  1.5],
                           [ 0. ,  0. , -2. ],
                           [ 0. ,  1.1,  0. ],
                           [-1.1,  0. ,  1.2]])
    assert np.all(test_hull(hull, make_homogeneous(in_points)))
    assert np.all(test_hull(hull, make_homogeneous(bound_points)))
    assert not np.any(test_hull(hull, make_homogeneous(out_points)))

    assert np.all(test_hull(hull, make_homogeneous(in_points), n_partitions=4))
    assert np.all(test_hull(hull, make_homogeneous(bound_points), n_partitions=4))
    assert not np.any(test_hull(hull, make_homogeneous(out_points), n_partitions=4))

def random_hemisphere(
    rand, n: int, radius: float = 1,
    centre: tuple[float, float, float] = (0, 0, 0),
    theta_limit=np.pi/2,
) -> np.ndarray:
    """
    Generate a 3D hemisphere with randomly-distributed vertices. The "cut" face of the hemisphere is
    not guaranteed to be exact when processed as a convex hull.
    """
    phi = rand.uniform(low=0, high=2*np.pi, size=n)
    theta = rand.uniform(low=0, high=theta_limit, size=n)
    cx = np.sin(theta)*np.cos(phi)
    cy = np.sin(theta)*np.sin(phi)
    cz = np.cos(theta)
    return radius*np.stack((cx, cy, cz), axis=1) + centre

def hemisphere_test() -> None:
    """
    Unit test for a hemisphere-shaped ("bowl") convex hull.
    """
    rand = default_rng(seed=0)

    centre = 5, 10, 12  # not the barycentre: only the centre of the "cut" face

    hull = ConvexHull(random_hemisphere(rand, n=100, radius=2, centre=centre))

    test_in = make_homogeneous(random_hemisphere(rand, n=150, radius=1.5, centre=centre, theta_limit=np.pi*0.4))
    indicators = test_hull(hull, test_in)
    assert np.all(indicators)

    test_close = make_homogeneous(random_hemisphere(rand, n=150, radius=1.975, centre=centre, theta_limit=np.pi*0.45))
    indicators = test_hull(hull, test_close)
    mean = indicators.mean()
    assert 0.48 <= mean <= 0.52  # 0.5067 for this seed

    test_out = make_homogeneous(random_hemisphere(rand, n=150, radius=2.1, centre=centre))
    indicators = test_hull(hull, test_out)
    assert not np.any(indicators)

def time_test() -> None:
    """
    Output timings for non-partitioned and partitioned configurations of a hemisphere hull test
    """
    rand = default_rng(seed=0)
    n = 10_000
    hull_pts = random_hemisphere(rand, n=n)
    test_pts = random_hemisphere(rand, n=n, radius=0.9999)

    t0 = monotonic()
    homogeneous = make_homogeneous(test_pts)
    t1 = monotonic()
    hull = ConvexHull(hull_pts)
    t2 = monotonic()
    print('n =', n)
    print(f'make_homogeneous: {t1-t0:.4f} s')
    print(f'ConvexHull:       {t2-t1:.4f} s')

    t3 = monotonic()
    indicators = test_hull(hull, homogeneous)
    t4 = monotonic()
    print(f'test_hull(part=   1): {t4-t3:.4f} s, {np.count_nonzero(indicators)} inside')

    for part in (5, 10, 20, 50, 100, 200, 500, 1_000):
        t5 = monotonic()
        indicators = test_hull(hull, homogeneous, n_partitions=part)
        t6 = monotonic()
        print(f'test_hull(part={part:4}): {t6-t5:.4f} s, {np.count_nonzero(indicators)} inside')

if __name__ == '__main__':
    cube_test()
    hemisphere_test()
    time_test()

个字符