This is slightly old, but I just wanted to mention one way to find Delaunay triangulations in arbitrary dimensions N. As explained in Wikipedia, if you know how to compute a convex hull in N+1 dimensions (which is supported by e.g. Polyhedra.jl) you can easily get Delaunay triangulations of N-dimensional points ps... by appending |ps[i]|^2 to each point ps[i], finding the convex hull, and then projecting back to N dimensions by removing the last coordinate in the vertices of the resulting hull mesh.
EDIT: I’m not so sure you can use Polyhedra.jl to decompose a convex hull into its faces in more than 3 dimensions, but you can certainly do it in 3D