Ticket #13517: trac_#13517_voronoi_first_try.patch

File trac_#13517_voronoi_first_try.patch, 3.0 KB (added by moritz, 9 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 b  
    1717
    1818
    1919import toric_plotter
     20
     21from 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
    - +  
     1"""
     2Voronoi diagram
     3
     4This module provides the function :func: `Voronoi_diagram` for computing the
     5Voronoi diagram of a finite list of points in \RR^d.
     6"""
     7
     8#*****************************************************************************
     9#       Copyright (C) 2012 Moritz Firsching <moritz@math.fu-berlin.de>
     10#
     11#  Distributed under the terms of the GNU General Public License (GPL)
     12#
     13#                  http://www.gnu.org/licenses/
     14#*****************************************************************************
     15
     16def Voronoi_diagram(points,**kwds):
     17    r"""
     18    Compute the Voronoi diagram of a list of points
     19
     20    INPUT:
     21
     22    - ``points`` -- list, ``points[i]`` should be a point in \RR^d
     23   
     24    OUTPUT:
     25
     26    - a list of polyhedra that are the voronoi cells
     27
     28    .. NOTE::
     29
     30        preliminary version.
     31
     32    ALGORITHM:
     33
     34    We use hyperplanes tangent to the paraboloid one dimension higher to
     35    get a convex polyhedron and then project back to one dimension lower
     36
     37    See for example [M2002]
     38
     39    REFERENCES:
     40
     41    ..  [M2002]
     42        Jiri Matousek,
     43        "Lectures on Discrete Geometry", Springer, Ch.5.7, p.118.
     44    AUTHORS:
     45
     46    - Moritz Firsching (2012-09-21)
     47   
     48    EXAMPLES:
     49    Get the Voronoi diagram for some points in \RR^2 and plot it:
     50
     51    d=2;   #dimension, plotting works well only in dimension 2
     52    n=7;   #number of points
     53    s=[[random() for k in range(d)] for i in range(n)]
     54 
     55    P=voronoi(s)
     56           
     57    S=line([])
     58    for i,j in enumerate(s):
     59        S+=(P[i]).render_solid(color=rainbow(n)[i], zorder=1)
     60        S+=point(j, color=rainbow(n)[i], pointsize=10,zorder=3)
     61        S+=point(vector(j), color='black',pointsize=20,zorder=2)
     62    show(S, xmax=1,xmin=0,ymax=1,ymin=0)
     63
     64    """
     65    from sage.geometry.polyhedron.constructor import Polyhedron
     66    from sage.all import RDF
     67
     68
     69    n=len(points)
     70    d=len(points[0])
     71    P=[]
     72    e=[([sum(i[k]**2 for k in range(d))]+[(-2)*i[l] for l in range(d)]+[1]) for i in points]
     73    p=Polyhedron(ieqs = e, base_ring=RDF)
     74    for i in range(len(points)):
     75        equ=p.Hrepresentation(i)
     76        pvert=[[u[k] for k in range(d)] for u in equ.incident() if u.is_vertex()]
     77        prays=[[u[k] for k in range(d)] for u in equ.incident() if u.is_ray()]
     78        pline=[[u[k] for k in range(d)] for u in equ.incident() if u.is_line()]
     79        P.append(Polyhedron(vertices=pvert, lines=pline, rays=prays, base_ring=RDF))
     80    print P
     81    return P
     82   
     83