Hi Sander,
ooops, how embarrassing! Thanks for spotting!
EDIT:
Here comes a hopefully correct approach: I am using my 9 points of SpherePoints[9]
around the origin {0,0,0} and assume the 9-dice will be the center cell of a 3D Voronoi mesh. To my great luck @Chip Hurst has provided a nice algorithm for this, so all credits should definitely go to him!
pts = Append[SpherePoints[9], {0, 0, 0}];
pad[\[Delta]_][{min_, max_}] := {min,
max} + \[Delta] (max - min) {-1, 1};
VoronoiCells[pts_] /;
MatrixQ[pts, NumericQ] && 2 <= Last[Dimensions[pts]] <= 3 :=
Block[{bds, dm, conn, adj, lc, pc, cpts, hpts, hns, hp, vcells},
bds = pad[.1] /@ MinMax /@ Transpose[pts];
dm = DelaunayMesh[pts];
conn = dm["ConnectivityMatrix"[0, 1]];
adj = conn.Transpose[conn];
lc = conn["MatrixColumns"];
pc = adj["MatrixColumns"];
cpts = MeshCoordinates[dm];
vcells =
Table[hpts = PropertyValue[{dm, {1, lc[[i]]}}, MeshCellCentroid];
hns =
Transpose[
Transpose[cpts[[DeleteCases[pc[[i]], i]]]] - cpts[[i]]];
hp = MapThread[HalfSpace, {hns, hpts}];
BoundaryDiscretizeGraphics[#, PlotRange -> bds] & /@ hp, {i,
MeshCellCount[dm, 0]}];
AssociationThread[cpts, RegionIntersection @@@ vcells]];
vc = VoronoiCells[pts];
This way we end up with:

The whole thing looks like (again using code from Chip Hurst):

Regards -- Henrik