ASKSAGE: Sage Q&A Forum - Individual question feedhttps://ask.sagemath.org/questions/Q&A Forum for SageenCopyright Sage, 2010. Some rights reserved under creative commons license.Thu, 19 Mar 2020 10:46:15 -0500VoronoiDiagram returns empty regionshttps://ask.sagemath.org/question/49749/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?
Please advise and thanks in advance.
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()Sun, 02 Feb 2020 10:43:31 -0600https://ask.sagemath.org/question/49749/voronoidiagram-returns-empty-regions/Answer by Jonathan Kliem for <p>Dear all,</p>
<p>Given the set of 2 points</p>
<pre><code>P = [[-2687.19000000000, -2088.53000000000], [-2686.81000000000, -2084.19000000000]]
</code></pre>
<p>I would like to compute the voronoi diagram of $P$. When doing so i saw that the regions of the voronoi are:</p>
<pre><code>{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}
</code></pre>
<p>In my code, I do need that those regions are not empty! Why is this happening and how can i fix it?</p>
<p>Please advise and thanks in advance.</p>
<p>P.s. Note that when I make the points a little I get that the voronoi is not empty.</p>
<p>Full Code:</p>
<pre><code>P = [[-2687.19, -2088.53], [-2686.81, -2084.19]]
V = VoronoiDiagram(P)
R = V.regions()
</code></pre>
https://ask.sagemath.org/question/49749/voronoidiagram-returns-empty-regions/?answer=50303#post-id-50303I think I have fixed the bug in [ticket 29176](https://trac.sagemath.org/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()
for adjacencies in data:
assert adjacencies[2] == ':', "Not a line of adjacency data"
cdd_vertex = int(adjacencies[0])
count = int(adjacencies[1])
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)] + [1])
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}Thu, 19 Mar 2020 10:46:15 -0500https://ask.sagemath.org/question/49749/voronoidiagram-returns-empty-regions/?answer=50303#post-id-50303Answer by jipilab for <p>Dear all,</p>
<p>Given the set of 2 points</p>
<pre><code>P = [[-2687.19000000000, -2088.53000000000], [-2686.81000000000, -2084.19000000000]]
</code></pre>
<p>I would like to compute the voronoi diagram of $P$. When doing so i saw that the regions of the voronoi are:</p>
<pre><code>{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}
</code></pre>
<p>In my code, I do need that those regions are not empty! Why is this happening and how can i fix it?</p>
<p>Please advise and thanks in advance.</p>
<p>P.s. Note that when I make the points a little I get that the voronoi is not empty.</p>
<p>Full Code:</p>
<pre><code>P = [[-2687.19, -2088.53], [-2686.81, -2084.19]]
V = VoronoiDiagram(P)
R = V.regions()
</code></pre>
https://ask.sagemath.org/question/49749/voronoidiagram-returns-empty-regions/?answer=49869#post-id-49869This looks like a bug. I have reported it on trac: [ticket 29176](https://trac.sagemath.org/ticket/29176)
Thanks for reporting this!Mon, 10 Feb 2020 06:47:43 -0600https://ask.sagemath.org/question/49749/voronoidiagram-returns-empty-regions/?answer=49869#post-id-49869