# Ticket #13517: trac_#13517_voronoi_first_try.patch

File trac_#13517_voronoi_first_try.patch, 3.0 KB (added by moritz, 10 years ago)

a first try

• ## sage/geometry/all.py

```# HG changeset patch
# User Moritz Firsching <moritz@math.fu-berlin.de>
# Date 1348242724 -7200
# Node ID 5f891da9d858306458820757c4c55b9a059e5bd8
# Parent  d06cf4b2215d37d3a87a58f65ac53234502dd471
Trac #13517 a first try

diff --git a/sage/geometry/all.py b/sage/geometry/all.py```
 a import toric_plotter from voronoi_diagram import *
• ## new file sage/geometry/voronoi_diagram.py

```diff --git a/sage/geometry/voronoi_diagram.py b/sage/geometry/voronoi_diagram.py
new file mode 100644```
 - """ Voronoi diagram This module provides the function :func: `Voronoi_diagram` for computing the Voronoi diagram of a finite list of points in \RR^d. """ #***************************************************************************** #       Copyright (C) 2012 Moritz Firsching # #  Distributed under the terms of the GNU General Public License (GPL) # #                  http://www.gnu.org/licenses/ #***************************************************************************** def Voronoi_diagram(points,**kwds): r""" Compute the Voronoi diagram of a list of points INPUT: - ``points`` -- list, ``points[i]`` should be a point in \RR^d OUTPUT: - a list of polyhedra that are the voronoi cells .. NOTE:: preliminary version. ALGORITHM: We use hyperplanes tangent to the paraboloid one dimension higher to get a convex polyhedron and then project back to one dimension lower See for example [M2002] REFERENCES: ..  [M2002] Jiri Matousek, "Lectures on Discrete Geometry", Springer, Ch.5.7, p.118. AUTHORS: - Moritz Firsching (2012-09-21) EXAMPLES: Get the Voronoi diagram for some points in \RR^2 and plot it: d=2;   #dimension, plotting works well only in dimension 2 n=7;   #number of points s=[[random() for k in range(d)] for i in range(n)] P=voronoi(s) S=line([]) for i,j in enumerate(s): S+=(P[i]).render_solid(color=rainbow(n)[i], zorder=1) S+=point(j, color=rainbow(n)[i], pointsize=10,zorder=3) S+=point(vector(j), color='black',pointsize=20,zorder=2) show(S, xmax=1,xmin=0,ymax=1,ymin=0) """ from sage.geometry.polyhedron.constructor import Polyhedron from sage.all import RDF n=len(points) d=len(points[0]) P=[] e=[([sum(i[k]**2 for k in range(d))]+[(-2)*i[l] for l in range(d)]+[1]) for i in points] p=Polyhedron(ieqs = e, base_ring=RDF) for i in range(len(points)): equ=p.Hrepresentation(i) pvert=[[u[k] for k in range(d)] for u in equ.incident() if u.is_vertex()] prays=[[u[k] for k in range(d)] for u in equ.incident() if u.is_ray()] pline=[[u[k] for k in range(d)] for u in equ.incident() if u.is_line()] P.append(Polyhedron(vertices=pvert, lines=pline, rays=prays, base_ring=RDF)) print P return P