Re: Voronoi Volume calculation
- To: mathgroup at smc.vnet.net
- Subject: [mg50590] Re: Voronoi Volume calculation
- From: "Steve Luttrell" <steve_usenet at _removemefirst_luttrell.org.uk>
- Date: Sat, 11 Sep 2004 06:44:51 -0400 (EDT)
- References: <chp91s$jln$1@smc.vnet.net>
- Sender: owner-wri-mathgroup at wolfram.com
Here is a possible way of doing it based on using the Boole function in the
Calculus`Integration` package. It computes the volume of a Voronoi cell very
quickly in 2 dimensions, in a few tens of seconds in 3 dimensions, but took
too long in 4 dimensions for me to get any answer out (I use a 2GHz Athlon
PC running Windows XP). I presume the method will work if you are patient.
Input the package.
<< Calculus`Integration`
Define a Eucliean distance function.
distsq[x1_, x2_] := #.# &[x1 - x2];
Generate a set of points in the region [-1,1] X [-1,1]. Here I generate
rationals rather than reals.
denom = 10;
numpoints = 10;
points = Table[Random[Integer, {-denom, denom}]/denom, {numpoints}, {2}];
Use the origin as the reference point whose Voronoi cell's volume you want
to compute.
Compute the distance of an arbitrary point {x,y} from the origin.
d0sq = distsq[{x, y}, {0, 0}];
Compute a list of inequality constraints that ensure that {x,y} lies closer
to the origin than to any of the above points. If you have the Delaunay
triangulation available then (for efficiency) you can use only those points
whose Voronoi cells share a common side with the cell at the origin;
Mathematica can compute a 2D Delaunay triangulation but I can't find a
higher dimensional version, so I have NOT made use of Delaunay triangulation
information here.
boundslist = Map[Simplify[distsq[{x, y}, #] - d0sq] >= 0 &, points];
Convert this list of constraints into a single Boolean constraint.
constraint = Apply[And, boundslist];
Finally compute the 2-dimensional Voronoi cell volume. I have assumed that
the Voronoi cell around the origin as computed using the above list of
points is contained within the region [-1,1] X [-1,1]; and if this condition
happens to be false (e.g. not enough points or the points are arranged in an
unlucky configuration) then the Voronoi cell is truncated by its
intersection with the boundary of [-1,1] X [-1,1].
Integrate[Boole[constraint], {x, -1, 1}, {y, -1, 1}]
That deals with the 2-dimensional case, and evaluates the result very
quickly.
Now gather all the code together in the 3-dimensional case to obtain the
following code:
<< Calculus`Integration`
distsq[x1_, x2_] := #.# &[x1 - x2];
denom = 10;
numpoints = 10;
points = Table[Random[Integer, {-denom, denom}]/denom, {numpoints}, {3}];
d0sq = distsq[{x, y, z}, {0, 0, 0}];
boundslist=Map[Simplify[distsq[{x,y,z},#]-d0sq] >= 0 &,points];
constraint = Apply[And, boundslist];
Integrate[Boole[constraint], {x, -1, 1}, {y, -1, 1}, {z, -1, 1}]
The result is a rational number with an impressively large number of digits
in its numerator and denominator.
The 4-dimensional case is analogous, but needs more time to run than my
patience allowed.
Good kuck!
Steve Luttrell
"Ran" <efeldesh at zahav.net.il> wrote in message
news:chp91s$jln$1 at smc.vnet.net...
> Hi everyone,
>
> I was trying to find the "volume" of a voronoi shape for a set of
> points in a 4-D space.
> what is the best way to do it in Mathematica?
>
> couldn't find any remarks on volume of these shapes at the wolfarm
> site
>
> thanks ahead,
> - Ran :)
>