Open adamkrawczyk opened 4 years ago
Hi, I found this package searching for an alternative to the using of Voronoi diagram method from the not-yet-released-version 1.8 of the Shapely package which is based on the Geos package.
With that package, using a list of points like this:
(20.1273, 18.7303), (26.5107, 18.7303), (20.1273, 23.8437), (26.5107, 23.8437)
[Those are the vertexes of a rectangular box.]
I had an Exception which made the algorithm not usable for me. Posting on the Issues page of Shapely, a contributor replied that the same could be seen on PostGIS (using Geos too) so the issue is upstream, in the Geos algorithm for Voronoi diagram.
So I searched for an alternative and found this package Yatoom/voronoi, with which I tried the same series of points. Got the error you posted above.
Only that in this case I only had to make a couple of changes to make it work. BTW, I've tried with 64 points in a grid and it worked after making the simple changes I've mentioned.
I think that this is an ugly hack but it may point the developer in the right direction.
The changes I made were: In voronoi.algorithm.handle_circle_event() method: just after "Updating breakpoint" I test if the updated attribute is None and if None I return False. `
self.beach_line, updated, removed, left, right = self.update_breakpoints(
self.beach_line, self.sweep_line, arc_node, predecessor, successor)
if updated is None:
return None
If the method reach it's end I return True.
and in voronoi.algorithm.handle_circle_event() method:
When handling circle events, when handling the event,
instead of
self.handle_circle_event(event, verbose=verbose)
I replaced it with:
`
Handle the event
res = self.handle_circle_event(event, verbose=verbose)
if res is None: continue`
PS: Sorry for the formatting, I'm not used to ...
Best regards, Marius
LE: Well, while the Matplotlib plot that the create_diagram() method display looks correct for 64 points, in reality with the modification above, after plotting the actual resulting Voronoi polygons in another plot, some of the polygons are fused together. So the modifications are not adequate. https://ibb.co/2npRRy3
Thanks @mars0001, there seems to be an issue with this case indeed. I will investigate this further.
I'm putting this here so that we have a reproducible script.
import numpy as np
from voronoi import Voronoi, Polygon
points = [(20.1273, 18.7303), (26.5107, 18.7303), (20.1273, 23.8437), (26.5107, 23.8437)]
polygon = Polygon([
(15., 15.),
(15., 30.),
(30., 30.),
(30., 15.),
])
v = Voronoi(polygon)
v.create_diagram(points=points, vis_steps=True, verbose=False, vis_before_clipping=True, vis_result=True, vis_tree=False)
This issue was caused by rounding errors, which in turn caused the circle events to be executed in the wrong order. Thanks for your help. If you find anything else that doesn't work, please let me know! :smile:
Hi Jeroen, ( @Yatoom )
Thank you for your effort!
Yet, using the latest version, the voronoi diagram fail with the same error from the first post for the following points:
points = [(0.3308, 0.204), (10.1432, 0.204), (19.9556, 0.204), (29.768, 0.204), (0.3308, 7.942), (10.1432, 7.942), (19.9556, 7.942), (29.768, 7.942), (0.3308, 15.6801), (10.1432, 15.6801), (19.9556, 15.6801), (29.768, 15.6801), (0.3308, 23.4181), (10.1432, 23.4181), (19.9556, 23.4181), (29.768, 23.4181)]
I think it has to do with the fact that the source points are on a square grid.
Best regards, Marius
Thanks Marius/@mars0001! In that case let me reopen this issue.
Those grids are indeed a bit tricky as they require special/careful handling when multiple events happen on the same line (see https://github.com/Yatoom/voronoi/issues/2#issuecomment-586625894 for more information on this)
Using your points, I'll try to solve this one next.
Hi @mars0001, it turned out to be rounding errors again. I think it should be fixed now :+1:
Hi Jeroen, ( @Yatoom ),
Now it works for a 4 by 4 grid but not for a higher order grid (5x4, 5x5 and higher).
Points list for which it does not work:
5x4: [[0.3308, 0.204], [7.6901, 0.204], [15.0494, 0.204], [22.4087, 0.204], [29.768, 0.204], [0.3308, 7.942], [7.6901, 7.942], [15.0494, 7.942], [22.4087, 7.942], [29.768, 7.942], [0.3308, 15.6801], [7.6901, 15.6801], [15.0494, 15.6801], [22.4087, 15.6801], [29.768, 15.6801], [0.3308, 23.4181], [7.6901, 23.4181], [15.0494, 23.4181], [22.4087, 23.4181], [29.768, 23.4181]]
5x5: [[0.3308, 0.204], [7.6901, 0.204], [15.0494, 0.204], [22.4087, 0.204], [29.768, 0.204], [0.3308, 6.0075], [7.6901, 6.0075], [15.0494, 6.0075], [22.4087, 6.0075], [29.768, 6.0075], [0.3308, 11.811], [7.6901, 11.811], [15.0494, 11.811], [22.4087, 11.811], [29.768, 11.811], [0.3308, 17.6146], [7.6901, 17.6146], [15.0494, 17.6146], [22.4087, 17.6146], [29.768, 17.6146], [0.3308, 23.4181], [7.6901, 23.4181], [15.0494, 23.4181], [22.4087, 23.4181], [29.768, 23.4181]]
6x6: [[0.3308, 0.204], [6.2182, 0.204], [12.1057, 0.204], [17.9931, 0.204], [23.8806, 0.204], [29.768, 0.204], [0.3308, 4.8468], [6.2182, 4.8468], [12.1057, 4.8468], [17.9931, 4.8468], [23.8806, 4.8468], [29.768, 4.8468], [0.3308, 9.4896], [6.2182, 9.4896], [12.1057, 9.4896], [17.9931, 9.4896], [23.8806, 9.4896], [29.768, 9.4896], [0.3308, 14.1325], [6.2182, 14.1325], [12.1057, 14.1325], [17.9931, 14.1325], [23.8806, 14.1325], [29.768, 14.1325], [0.3308, 18.7753], [6.2182, 18.7753], [12.1057, 18.7753], [17.9931, 18.7753], [23.8806, 18.7753], [29.768, 18.7753], [0.3308, 23.4181], [6.2182, 23.4181], [12.1057, 23.4181], [17.9931, 23.4181], [23.8806, 23.4181], [29.768, 23.4181]]
6x7: [[0.3308, 0.204], [6.2182, 0.204], [12.1057, 0.204], [17.9931, 0.204], [23.8806, 0.204], [29.768, 0.204], [0.3308, 4.073], [6.2182, 4.073], [12.1057, 4.073], [17.9931, 4.073], [23.8806, 4.073], [29.768, 4.073], [0.3308, 7.942], [6.2182, 7.942], [12.1057, 7.942], [17.9931, 7.942], [23.8806, 7.942], [29.768, 7.942], [0.3308, 11.811], [6.2182, 11.811], [12.1057, 11.811], [17.9931, 11.811], [23.8806, 11.811], [29.768, 11.811], [0.3308, 15.6801], [6.2182, 15.6801], [12.1057, 15.6801], [17.9931, 15.6801], [23.8806, 15.6801], [29.768, 15.6801], [0.3308, 19.5491], [6.2182, 19.5491], [12.1057, 19.5491], [17.9931, 19.5491], [23.8806, 19.5491], [29.768, 19.5491], [0.3308, 23.4181], [6.2182, 23.4181], [12.1057, 23.4181], [17.9931, 23.4181], [23.8806, 23.4181], [29.768, 23.4181]]
7x7: [[0.3308, 0.204], [5.237, 0.204], [10.1432, 0.204], [15.0494, 0.204], [19.9556, 0.204], [24.8618, 0.204], [29.768, 0.204], [0.3308, 4.073], [5.237, 4.073], [10.1432, 4.073], [15.0494, 4.073], [19.9556, 4.073], [24.8618, 4.073], [29.768, 4.073], [0.3308, 7.942], [5.237, 7.942], [10.1432, 7.942], [15.0494, 7.942], [19.9556, 7.942], [24.8618, 7.942], [29.768, 7.942], [0.3308, 11.811], [5.237, 11.811], [10.1432, 11.811], [15.0494, 11.811], [19.9556, 11.811], [24.8618, 11.811], [29.768, 11.811], [0.3308, 15.6801], [5.237, 15.6801], [10.1432, 15.6801], [15.0494, 15.6801], [19.9556, 15.6801], [24.8618, 15.6801], [29.768, 15.6801], [0.3308, 19.5491], [5.237, 19.5491], [10.1432, 19.5491], [15.0494, 19.5491], [19.9556, 19.5491], [24.8618, 19.5491], [29.768, 19.5491], [0.3308, 23.4181], [5.237, 23.4181], [10.1432, 23.4181], [15.0494, 23.4181], [19.9556, 23.4181], [24.8618, 23.4181], [29.768, 23.4181]]
8x8: [[0.3308, 0.204], [4.5361, 0.204], [8.7414, 0.204], [12.9467, 0.204], [17.1521, 0.204], [21.3574, 0.204], [25.5627, 0.204], [29.768, 0.204], [0.3308, 3.5203], [4.5361, 3.5203], [8.7414, 3.5203], [12.9467, 3.5203], [17.1521, 3.5203], [21.3574, 3.5203], [25.5627, 3.5203], [29.768, 3.5203], [0.3308, 6.8366], [4.5361, 6.8366], [8.7414, 6.8366], [12.9467, 6.8366], [17.1521, 6.8366], [21.3574, 6.8366], [25.5627, 6.8366], [29.768, 6.8366], [0.3308, 10.1529], [4.5361, 10.1529], [8.7414, 10.1529], [12.9467, 10.1529], [17.1521, 10.1529], [21.3574, 10.1529], [25.5627, 10.1529], [29.768, 10.1529], [0.3308, 13.4692], [4.5361, 13.4692], [8.7414, 13.4692], [12.9467, 13.4692], [17.1521, 13.4692], [21.3574, 13.4692], [25.5627, 13.4692], [29.768, 13.4692], [0.3308, 16.7855], [4.5361, 16.7855], [8.7414, 16.7855], [12.9467, 16.7855], [17.1521, 16.7855], [21.3574, 16.7855], [25.5627, 16.7855], [29.768, 16.7855], [0.3308, 20.1018], [4.5361, 20.1018], [8.7414, 20.1018], [12.9467, 20.1018], [17.1521, 20.1018], [21.3574, 20.1018], [25.5627, 20.1018], [29.768, 20.1018], [0.3308, 23.4181], [4.5361, 23.4181], [8.7414, 23.4181], [12.9467, 23.4181], [17.1521, 23.4181], [21.3574, 23.4181], [25.5627, 23.4181], [29.768, 23.4181]]
And so on... I hope those can help in pinpointing the bugs.
Thanks, Marius
Hi @mars0001, thanks for these grids! Those really help.
There was another place where rounding errors where causing trouble. This time it was when querying the binary tree to find the breakpoints that needed to be updated. I used a workaround for this: instead of returning None
when the matching break-point could not be found, it will look into the other branches of the three as well. I think it shouldn't impact performance too much, because normally, it should just be able to find the right breakpoint without the workaround.
Let me know if it works for you :smile:
Hi Jeroen (@Yatoom ),
It seems that the latest fix solved the issue with point on grid. Here is a picture with 20x20 points (400 points in a gid). https://ibb.co/sK9j5KD
But I did found that by adding manual points (is done by clicking with the mouse) and perhaps adding overlapped points (points with same or very close coordinates) I get this result. https://ibb.co/pJFwtqn
The black points are actually the clicked points that are used in the algorithm. Each seed point should have its own polygon but that is not the case ...
I get the Shapely polygons using something like this:
def generate_voronoi_geometry(self, pts):
env = self.solid_geo.envelope # an outline of a complex geometry
fact = 1 if self.units == 'MM' else 0.039
env = env.buffer(fact) # create a polygon from the outline
# voronoi_poly is the Polygon from Voronoi package so it does not get confused with the Shapely Polygon
env_poly = voronoi_poly(list(env.exterior.coords)
new_pts = [[pt.x, pt.y] for pt in pts]
# Initialize the algorithm
v = Voronoi(env_poly)
# calculate the Voronoi diagram
try:
v.create_diagram(new_pts)
except Exception as e:
print("CNCJobObject.generate_voronoi_geometry_2() --> %s" % str(e))
# try again with slightly changed coordinates (random) perhaps will make it work
new_pts_2 = []
for pt_index in range(len(new_pts)):
new_pts_2.append([
new_pts[pt_index][0] + random.random() * 1e-06,
new_pts[pt_index][1] + random.random() * 1e-06
])
try:
v.create_diagram(new_pts_2)
except Exception:
print("Didn't work.")
return
new_voronoi = []
for p in v.points:
p_coords = [(coord.x, coord.y) for coord in p.get_coordinates()]
new_pol = Polygon(p_coords)
new_voronoi.append(new_pol)
new_voronoi = MultiPolygon(new_voronoi)
return new_voronoi
Should I open a new issue for this? Thank you, Marius
LE: in my case I can get away with storing manual points in a set (therefore eliminating duplicates) or test if it is already added, but maybe fixing it may make your package more foolproof.
Hi @mars0001, could you confirm whether this issue only pops up if there are duplicates?
Hi Jeroen @Yatoom ,
I did made sure that the duplicated points are not added at all and that fixed the issue.
Best regards, Marius
When point number is too large (around 30) Algorithm stops workin, and i can see this error
After limiting the number of points graph is computed properly