# Intersecting cones for weighted Voronoi diagram - DiscretizeRegion was unable to discretize the region

I want to make a diagram (and also a Graph of the vertex + edge config) of the orthogonal projection of the intersections of a set of conical surfaces in order to create weighted Voronoi diagrams as in these pictures:

additively weighted Voronoi diagram (heights vary)

multiplicatively weighted Voronoi diagram (radii vary)

I tried to do this using Mathematica's geometric regions tools but they're not working well for me. For example, here is some code:

(* Define a cone centred at x,y,z with radius r and height h \
by subtracting a disk from the boundary of a solid cone *)

cone[x_, y_, z_, r_, h_] :=
RegionDifference[RegionBoundary[Cone[{{x, y, z}, {x, y, z + h}}, r]],
TransformedRegion[Ball[{x, y, z}, r],
Function[{p}, {p[[1]], p[[2]], 0}]]]

(* Create a random set of cones sitting on the xy-plane*)

n = 3;

SeedRandom[123]

points = RandomReal[{-1, 1}, {n, 2}];
heights = RandomReal[{1/2, 1}, n];

cones = cone[#[[1, 1]], #[[1, 2]], 0, #[[2]], #[[3]]] & /@

(* Check discretize region for one cone*)

DiscretizeRegion[cones[[1]]]

(* Discretize union of intersections of cones *)

DiscretizeRegion[
RegionUnion[RegionIntersection /@ Subsets[cones, {2}]]]


For the particular random seed 123, DiscretizeRegion[cones[[1]]] outputs

However it fails for most other random seeds. And then the next line has failed for every seed I tried (and does so quite slowly).

What am I doing wrong and is there a better way to do this?

Thanks a lot.

EDIT: I hope to get this working nicely for about $n=60$.

Mathematica seems to prefer implicit regions, you can get the cones easily enough:

cone[xc_, yc_, r_, h_] := ImplicitRegion[(x - xc)^2 + (y - yc)^2 == r^2 (h - z)^2/h^2 &&
z >= 0 && z <= h, {x, y, z}];
n = 3; SeedRandom[123];
points = RandomReal[{-1, 1}, {n, 2}];
heights = RandomReal[{1/2, 1}, n];
cones = cone[#[[1, 1]], #[[1, 2]], #[[2]], #[[3]]] & /@
DiscretizeRegion[#, Method -> "RegionPlot3D"] & /@ cones


You can even get some kind of intersection with

intersect = BooleanRegion[BooleanCountingFunction[{2, n}, n], cones];
DiscretizeRegion[intersect, Method -> "Semialgebraic"]


though not without artefacts

If you just want to plot Voronoi diagram, better do it directly with:

cone[xc_, yc_, r_, h_] :=
h - h/r Sqrt[(#[[1]] - xc)^2 + (#[[2]] - yc)^2] &;
Plot3D[Max@Through@cones@{x, y}, {x, -1, 1}, {y, -1, 1}]


DensityPlot[Ordering[Through@cones@{x, y}, -1], {x, -1, 1}, {y, -1, 1}, PlotPoints -> 100]


• Thanks. This method doesn't seem to be effective though since it is giving an error for $n>4$ and is quite slow anyway. I'd like to have this working for larger, say $n=60$. Maybe using the Geometric Regions toolset is not the best way to about this? – Nasos Evangelou-Oost Oct 12 '15 at 1:12
• Maybe a lower-level approach could be better? For example, plotting the cone arrangement with Plot3D and looking at the view !from the top the intersection diagram is visible from the lighting effect within the 3D graphics and the image is produced instantly even for a very large number of cones. – Nasos Evangelou-Oost Oct 12 '15 at 1:20
• @AthanasiosEvangelou, you asked about Regions, if you just want to plot it, you can just plot it, I added some ways to do it – panda-34 Oct 12 '15 at 13:46