Computes power triangulation in N dimensions

power-triangulation.py

import numpy
from scipy.spatial import ConvexHull


def is_oriented_simplex(S):
    M = numpy.concatenate([S, numpy.ones((S.shape[0], 1))], axis = 1)
    return numpy.linalg.det(M) > 0


def get_power_triangulation(S, R):    
    def permut(A):
        return numpy.concatenate([A[:1], numpy.flip(A[1:])])

    # Compute the lifted weighted points
    S_norm = numpy.sum(S ** 2, axis = 1) - R ** 2
    S_lifted = numpy.concatenate([S, S_norm[:,None]], axis = 1)

    # Compute the convex hull of the lifted weighted points
    hull = ConvexHull(S_lifted)
    
    # Extract the Delaunay triangulation from the lower hull
    simplex_list = numpy.stack([simplex if is_oriented_simplex(S[simplex]) else permut(simplex) for simplex, eq in zip(hull.simplices, hull.equations) if eq[-2] <= 0], axis = 0)

    # Job done
    return simplex_list