# VoronoiDiagram returns empty regions

Dear all,

Given the set of 2 points

P = [[-2687.19000000000, -2088.53000000000], [-2686.81000000000, -2084.19000000000]]


I would like to compute the voronoi diagram of $P$. When doing so i saw that the regions of the voronoi are:

{P(-2687.19000000000, -2088.53000000000): The empty polyhedron in RDF^0,
P(-2686.81000000000, -2084.19000000000): A 1-dimensional polyhedron in RDF^2 defined as the convex hull of 1 vertex and 1 ray}


In my code, I do need that those regions are not empty! Why is this happening and how can i fix it?

P.s. Note that when I make the points a little I get that the voronoi is not empty.

Full Code:

P = [[-2687.19, -2088.53], [-2686.81, -2084.19]]
V = VoronoiDiagram(P)
R = V.regions()

edit retag close merge delete

Sort by » oldest newest most voted

I think I have fixed the bug in ticket 29176.

Hopefully, it will be available from Sage version 9.1+.

In the meantime you can use the following workaround.

from sage.all import *
from sage.geometry.polyhedron.parent import Polyhedra_RDF_cdd
from sage.geometry.polyhedron.backend_cdd import Polyhedron_RDF_cdd
class Polyhedron_RDF_cdd_workaround(Polyhedron_RDF_cdd):
def _init_from_cdd_output(self, cddout):
Polyhedron_RDF_cdd._init_from_cdd_output(self, cddout)
cddout = cddout.splitlines()

def parse_indices(count, cdd_indices, cdd_indices_to_sage_indices=None):
cdd_indices = [int(x) for x in cdd_indices]
if cdd_indices_to_sage_indices is None:
cdd_indices_to_sage_indices = {i: i-1 for i in cdd_indices}
if count < 0:
assert cdd_indices_to_sage_indices is not None, "Did not expect negative counts here"
count = -count
cdd_indices = list(set(cdd_indices_to_sage_indices.keys()) - set(cdd_indices))
assert count in [len(cdd_indices), len(cdd_indices) - 1]
assert count == len(cdd_indices)
return [cdd_indices_to_sage_indices[i] for i in cdd_indices if cdd_indices_to_sage_indices[i] is not None]

def parse_incidence_matrix(intro, data):
if 'incidence_matrix' in self.__dict__:
raise NotImplementedError("can not replace internal representation as this breaks caching")
N = len(self._Hrepresentation)
M = len(self._Vrepresentation)
ret = matrix(ZZ, M, N, 0)
data.pop(0)
data.reverse()
if cdd_vertex not in self._cdd_V_to_sage_V:
self._cdd_V_to_sage_V[cdd_vertex] = None
v = self._cdd_V_to_sage_V[cdd_vertex]
if v is None:
continue
for w in parse_indices(count, adjacencies[3:], self._cdd_H_to_sage_H):
if w is None:
continue
ret[v, w] = 1
ret.set_immutable()
self.incidence_matrix.set_cache(ret)

Polyhedron_RDF_cdd._parse_block(cddout, 'Vertex incidence', parse_incidence_matrix)

class Polyhedra_RDF_cdd_workaround(Polyhedra_RDF_cdd):
Element = Polyhedron_RDF_cdd_workaround

class VoronoiDiagram_workaround(VoronoiDiagram):
def __init__(self, points):
VoronoiDiagram.__init__(self, points)
self._P = {}
if self._n > 0:
e = [([sum(vector(i)[k] ** 2
for k in range(self._d))] +
[(-2) * vector(i)[l] for l in range(self._d)] + )
for i in self._points]
e = [[self._base_ring(i) for i in k] for k in e]
parent = Polyhedra_RDF_cdd_workaround(RDF, self._d + 1, 'cdd')
p = parent(None, [e, []])
for i in range(self._n):
parent = Polyhedra_RDF_cdd_workaround(RDF, self._d, 'cdd')
equ = p.Hrepresentation(i)
pvert = [[u[k] for k in range(self._d)] for u in equ.incident()
if u.is_vertex()]
prays = [[u[k] for k in range(self._d)] for u in equ.incident()
if u.is_ray()]
pline = [[u[k] for k in range(self._d)] for u in equ.incident()
if u.is_line()]
(self._P)[self._points[i]] = parent([pvert, prays, pline], None)


Then

P = [[-2687.19000000000, -2088.53000000000], [-2686.81000000000, -2084.19000000000]]
V = VoronoiDiagram_workaround(P)
R = V.regions()
R


will produce correctly

{P(-2687.19000000000, -2088.53000000000): A 2-dimensional polyhedron in RDF^2 defined as the convex hull of 1 vertex, 1 ray, 1 line,
P(-2686.81000000000, -2084.19000000000): A 2-dimensional polyhedron in RDF^2 defined as the convex hull of 1 vertex, 1 ray, 1 line}

more

This looks like a bug. I have reported it on trac: ticket 29176

Thanks for reporting this!

more