Skip to content

Commit acbb27c

Browse files
havatvnyalldawson
havatv
authored and
nyalldawson
committedAug 8, 2018
[processing] Update to VoronoiPolygons.py, some fixes for #8002
* Update to VoronoiPolygons.py, fixes #8002 adds comments, ++ * Fixed a spelling an an indentation mistake reported by travis * Found some new cases (vertical line clipping and a mistake in the handling of extreme points) (cherry-picked from 28fa839)
·
final-3_2_3final-3_2_2
1 parent ce3cba1 commit acbb27c

File tree

1 file changed

+233
-59
lines changed

1 file changed

+233
-59
lines changed
 

‎python/plugins/processing/algs/qgis/VoronoiPolygons.py

100644100755
Lines changed: 233 additions & 59 deletions
Original file line numberDiff line numberDiff line change
@@ -17,7 +17,7 @@
1717
***************************************************************************
1818
"""
1919

20-
__author__ = 'Victor Olaya'
20+
__authors__ = 'Victor Olaya, Håvard Tveite'
2121
__date__ = 'August 2012'
2222
__copyright__ = '(C) 2012, Victor Olaya'
2323

@@ -71,15 +71,18 @@ def __init__(self):
7171
super().__init__()
7272

7373
def initAlgorithm(self, config=None):
74-
self.addParameter(QgsProcessingParameterFeatureSource(
75-
self.INPUT, self.tr('Input layer'), [QgsProcessing.TypeVectorPoint]))
74+
self.addParameter(
75+
QgsProcessingParameterFeatureSource(
76+
self.INPUT, self.tr('Input layer'),
77+
[QgsProcessing.TypeVectorPoint]))
7678
self.addParameter(
7779
QgsProcessingParameterNumber(
7880
self.BUFFER, self.tr('Buffer region (% of extent)'),
7981
minValue=0.0, maxValue=9999999999, defaultValue=0.0))
80-
81-
self.addParameter(QgsProcessingParameterFeatureSink(
82-
self.OUTPUT, self.tr('Voronoi polygons'), type=QgsProcessing.TypeVectorPolygon))
82+
self.addParameter(
83+
QgsProcessingParameterFeatureSink(
84+
self.OUTPUT, self.tr('Voronoi polygons'),
85+
type=QgsProcessing.TypeVectorPolygon))
8386

8487
def name(self):
8588
return 'voronoipolygons'
@@ -94,12 +97,13 @@ def processAlgorithm(self, parameters, context, feedback):
9497
self.invalidSourceError(parameters, self.INPUT))
9598

9699
buf = self.parameterAsDouble(parameters, self.BUFFER, context)
97-
(sink, dest_id) = self.parameterAsSink(parameters, self.OUTPUT, context,
98-
source.fields(), QgsWkbTypes.Polygon, source.sourceCrs())
100+
(sink, dest_id) = self.parameterAsSink(parameters, self.OUTPUT,
101+
context, source.fields(),
102+
QgsWkbTypes.Polygon,
103+
source.sourceCrs())
99104
if sink is None:
100105
raise QgsProcessingException(
101106
self.invalidSinkError(parameters, self.OUTPUT))
102-
103107
outFeat = QgsFeature()
104108
extent = source.sourceExtent()
105109
extraX = extent.width() * (buf / 100.0)
@@ -115,7 +119,11 @@ def processAlgorithm(self, parameters, context, feedback):
115119
pts = []
116120
ptDict = {}
117121
ptNdx = -1
118-
122+
# Find the minimum and maximum x and y for the input points
123+
xmin = width
124+
xmax = 0
125+
ymin = height
126+
ymax = 0
119127
features = source.getFeatures()
120128
total = 100.0 / source.featureCount() if source.featureCount() else 0
121129
for current, inFeat in enumerate(features):
@@ -128,93 +136,131 @@ def processAlgorithm(self, parameters, context, feedback):
128136
pts.append((x, y))
129137
ptNdx += 1
130138
ptDict[ptNdx] = inFeat.id()
139+
if x < xmin:
140+
xmin = x
141+
if y < ymin:
142+
ymin = y
143+
if x > xmax:
144+
xmax = x
145+
if y > ymax:
146+
ymax = y
131147
feedback.setProgress(int(current * total))
132-
148+
if xmin == xmax or ymin == ymax:
149+
raise QgsProcessingException('The extent of the input points is '
150+
'not a polygon (all the points are '
151+
'on a vertical or horizontal line) '
152+
'- cannot make a Voronoi diagram!')
153+
xyminmax = [xmin, ymin, xmax, ymax]
133154
if len(pts) < 3:
134155
raise QgsProcessingException(
135156
self.tr('Input file should contain at least 3 points. Choose '
136157
'another file and try again.'))
137-
158+
# Eliminate duplicate points
138159
uniqueSet = set(item for item in pts)
139160
ids = [pts.index(item) for item in uniqueSet]
140-
sl = voronoi.SiteList([voronoi.Site(i[0], i[1], sitenum=j) for (j,
141-
i) in enumerate(uniqueSet)])
161+
sl = voronoi.SiteList([voronoi.Site(i[0], i[1], sitenum=j)
162+
for (j, i) in enumerate(uniqueSet)])
142163
voronoi.voronoi(sl, c)
143-
inFeat = QgsFeature()
144-
145-
current = 0
146164
if len(c.polygons) == 0:
147165
raise QgsProcessingException(
148166
self.tr('There were no polygons created.'))
149167

168+
inFeat = QgsFeature()
169+
current = 0
150170
total = 100.0 / len(c.polygons)
151-
171+
# Clip each of the generated "polygons"
152172
for (site, edges) in list(c.polygons.items()):
153173
if feedback.isCanceled():
154174
break
155-
156175
request = QgsFeatureRequest().setFilterFid(ptDict[ids[site]])
157176
inFeat = next(source.getFeatures(request))
158-
lines = self.clip_voronoi(edges, c, width, height, extent)
159-
160-
geom = QgsGeometry.fromMultiPointXY(lines)
161-
geom = QgsGeometry(geom.convexHull())
177+
boundarypoints = self.clip_voronoi(edges, c, width,
178+
height, extent,
179+
inFeat.geometry().asPoint(),
180+
xyminmax)
181+
ptgeom = QgsGeometry.fromMultiPointXY(boundarypoints)
182+
geom = QgsGeometry(ptgeom.convexHull())
162183
outFeat.setGeometry(geom)
163184
outFeat.setAttributes(inFeat.attributes())
164185
sink.addFeature(outFeat, QgsFeatureSink.FastInsert)
165-
166186
current += 1
167187
feedback.setProgress(int(current * total))
168-
169188
return {self.OUTPUT: dest_id}
170189

171-
def clip_voronoi(self, edges, c, width, height, extent):
190+
def clip_voronoi(self, edges, c, width, height, extent, point, xyminmax):
172191
"""Clip voronoi function based on code written for Inkscape.
173192
Copyright (C) 2010 Alvin Penner, penner@vaxxine.com
193+
Clips one Thiessen polygon (convex polygon) to extent
174194
"""
175195

196+
pt_x = point.x() - extent.xMinimum()
197+
pt_y = point.y() - extent.yMinimum()
198+
(xmin, ymin, xmax, ymax) = xyminmax
199+
176200
def clip_line(x1, y1, x2, y2, w, h):
177201
if x1 < 0 and x2 < 0:
202+
# Completely to the left
178203
return [0, 0, 0, 0]
179204
if x1 > w and x2 > w:
205+
# Completely to the right
206+
return [0, 0, 0, 0]
207+
if y1 < 0 and y2 < 0:
208+
# Completely below
209+
return [0, 0, 0, 0]
210+
if y1 > h and y2 > h:
211+
# Completely above
180212
return [0, 0, 0, 0]
181213
if x1 < 0:
214+
# First point to the left
182215
y1 = (y1 * x2 - y2 * x1) / (x2 - x1)
183216
x1 = 0
184217
if x2 < 0:
218+
# Last point to the left
185219
y2 = (y1 * x2 - y2 * x1) / (x2 - x1)
186220
x2 = 0
187221
if x1 > w:
222+
# First point to the right
188223
y1 = y1 + (w - x1) * (y2 - y1) / (x2 - x1)
189224
x1 = w
190225
if x2 > w:
226+
# Last point to the right
191227
y2 = y1 + (w - x1) * (y2 - y1) / (x2 - x1)
192228
x2 = w
193-
if y1 < 0 and y2 < 0:
194-
return [0, 0, 0, 0]
195-
if y1 > h and y2 > h:
196-
return [0, 0, 0, 0]
197229
if x1 == x2 and y1 == y2:
198230
return [0, 0, 0, 0]
199231
if y1 < 0:
232+
# First point below
200233
x1 = (x1 * y2 - x2 * y1) / (y2 - y1)
201234
y1 = 0
202235
if y2 < 0:
236+
# Second point below
203237
x2 = (x1 * y2 - x2 * y1) / (y2 - y1)
204238
y2 = 0
205239
if y1 > h:
240+
# First point above
206241
x1 = x1 + (h - y1) * (x2 - x1) / (y2 - y1)
207242
y1 = h
208243
if y2 > h:
244+
# Second point above
209245
x2 = x1 + (h - y1) * (x2 - x1) / (y2 - y1)
210246
y2 = h
211247
return [x1, y1, x2, y2]
212248

213-
lines = []
249+
bndpoints = []
214250
hasXMin = False
215251
hasYMin = False
216252
hasXMax = False
217253
hasYMax = False
254+
XMinNumber = 0
255+
XMaxNumber = 0
256+
YMinNumber = 0
257+
YMaxNumber = 0
258+
# The same line may appear twice for collinear input points,
259+
# so have to remember which lines have contributed
260+
XMinLine = -1
261+
XMaxLine = -1
262+
YMinLine = -1
263+
YMaxLine = -1
218264
for edge in edges:
219265
if edge[1] >= 0 and edge[2] >= 0:
220266
# Two vertices
@@ -227,15 +273,18 @@ def clip_line(x1, y1, x2, y2, w, h):
227273
height
228274
)
229275
elif edge[1] >= 0:
230-
# Only one vertex
276+
# Only one (left) vertex
231277
if c.lines[edge[0]][1] == 0:
232278
# Vertical line
233-
xtemp = c.lines[edge[0]][2] / c.lines[edge[0]][0]
234-
if c.vertices[edge[1]][1] > height / 2:
235-
ytemp = height
236-
else:
237-
ytemp = 0
279+
#xtemp = c.lines[edge[0]][2] / c.lines[edge[0]][0]
280+
xtemp = c.vertices[edge[1]][0]
281+
ytemp = 0 - 1
282+
#if c.vertices[edge[1]][1] > height / 2:
283+
# ytemp = height
284+
#else:
285+
# ytemp = 0
238286
else:
287+
# Create an end of the line at the right edge - OK
239288
xtemp = width
240289
ytemp = (c.lines[edge[0]][2] - width *
241290
c.lines[edge[0]][0]) / c.lines[edge[0]][1]
@@ -248,15 +297,18 @@ def clip_line(x1, y1, x2, y2, w, h):
248297
height
249298
)
250299
elif edge[2] >= 0:
251-
# Only one vertex
300+
# Only one (right) vertex
252301
if c.lines[edge[0]][1] == 0:
253302
# Vertical line
254-
xtemp = c.lines[edge[0]][2] / c.lines[edge[0]][0]
255-
if c.vertices[edge[2]][1] > height / 2:
256-
ytemp = height
257-
else:
258-
ytemp = 0.0
303+
#xtemp = c.lines[edge[0]][2] / c.lines[edge[0]][0]
304+
xtemp = c.vertices[edge[2]][0]
305+
ytemp = height + 1
306+
#if c.vertices[edge[2]][1] > height / 2:
307+
# ytemp = height
308+
#else:
309+
# ytemp = 0.0
259310
else:
311+
# End the line at the left edge - OK
260312
xtemp = 0.0
261313
ytemp = c.lines[edge[0]][2] / c.lines[edge[0]][1]
262314
[x1, y1, x2, y2] = clip_line(
@@ -267,31 +319,153 @@ def clip_line(x1, y1, x2, y2, w, h):
267319
width,
268320
height,
269321
)
322+
else:
323+
# No vertex, only a line
324+
if c.lines[edge[0]][1] == 0:
325+
# Vertical line - should not happen
326+
xtemp = c.lines[edge[0]][2] / c.lines[edge[0]][0]
327+
ytemp = 0.0
328+
xend = xtemp
329+
yend = height
330+
else:
331+
# End the line at both edges - ???
332+
xtemp = 0.0
333+
ytemp = c.lines[edge[0]][2] / c.lines[edge[0]][1]
334+
xend = width
335+
yend = (c.lines[edge[0]][2] - width *
336+
c.lines[edge[0]][0]) / c.lines[edge[0]][1]
337+
[x1, y1, x2, y2] = clip_line(
338+
xtemp,
339+
ytemp,
340+
xend,
341+
yend,
342+
width,
343+
height,
344+
)
270345
if x1 or x2 or y1 or y2:
271-
lines.append(QgsPointXY(x1 + extent.xMinimum(),
272-
y1 + extent.yMinimum()))
273-
lines.append(QgsPointXY(x2 + extent.xMinimum(),
274-
y2 + extent.yMinimum()))
346+
bndpoints.append(QgsPointXY(x1 + extent.xMinimum(),
347+
y1 + extent.yMinimum()))
348+
bndpoints.append(QgsPointXY(x2 + extent.xMinimum(),
349+
y2 + extent.yMinimum()))
275350
if 0 in (x1, x2):
276351
hasXMin = True
352+
if XMinLine != edge[0]:
353+
XMinNumber = XMinNumber + 1
354+
XMinLine = edge[0]
277355
if 0 in (y1, y2):
278356
hasYMin = True
357+
if YMinLine != edge[0]:
358+
YMinNumber = YMinNumber + 1
359+
YMinLine = edge[0]
279360
if height in (y1, y2):
280361
hasYMax = True
362+
if YMaxLine != edge[0]:
363+
YMaxNumber = YMaxNumber + 1
364+
YMaxLine = edge[0]
281365
if width in (x1, x2):
282366
hasXMax = True
283-
if hasXMin:
284-
if hasYMax:
285-
lines.append(QgsPointXY(extent.xMinimum(),
286-
height + extent.yMinimum()))
287-
if hasYMin:
288-
lines.append(QgsPointXY(extent.xMinimum(),
367+
if XMaxLine != edge[0]:
368+
XMaxNumber = XMaxNumber + 1
369+
XMaxLine = edge[0]
370+
371+
# Add auxiliary points for corner cases, if necessary (duplicate
372+
# points is not a problem - will be ignored later).
373+
# a) Extreme input points (lowest, leftmost, rightmost, highest)
374+
# A point can be extreme on both axis
375+
if pt_x == xmin: # leftmost point
376+
if XMinNumber == 0:
377+
bndpoints.append(QgsPointXY(extent.xMinimum(),
378+
extent.yMinimum()))
379+
bndpoints.append(QgsPointXY(extent.xMinimum(),
380+
height + extent.yMinimum()))
381+
elif XMinNumber == 1:
382+
if hasYMin:
383+
bndpoints.append(QgsPointXY(extent.xMinimum(),
384+
extent.yMinimum()))
385+
elif hasYMax:
386+
bndpoints.append(QgsPointXY(extent.xMinimum(),
387+
height + extent.yMinimum()))
388+
elif pt_x == xmax: # rightmost point
389+
if XMaxNumber == 0:
390+
bndpoints.append(QgsPointXY(width + extent.xMinimum(),
391+
extent.yMinimum()))
392+
bndpoints.append(QgsPointXY(width + extent.xMinimum(),
393+
height + extent.yMinimum()))
394+
elif XMaxNumber == 1:
395+
if hasYMin:
396+
bndpoints.append(QgsPointXY(width + extent.xMinimum(),
397+
extent.yMinimum()))
398+
elif HasYMax:
399+
bndpoints.append(QgsPointXY(width + extent.xMinimum(),
400+
height + extent.yMinimum()))
401+
if pt_y == ymin: # lowest point
402+
if YMinNumber == 0:
403+
bndpoints.append(QgsPointXY(extent.xMinimum(),
404+
extent.yMinimum()))
405+
bndpoints.append(QgsPointXY(width + extent.xMinimum(),
406+
extent.yMinimum()))
407+
elif YMinNumber == 1:
408+
if hasXMin:
409+
bndpoints.append(QgsPointXY(extent.xMinimum(),
410+
extent.yMinimum()))
411+
elif hasXMax:
412+
bndpoints.append(QgsPointXY(width + extent.xMinimum(),
413+
extent.yMinimum()))
414+
elif pt_y == ymax: # highest point
415+
if YMaxNumber == 0:
416+
bndpoints.append(QgsPointXY(extent.xMinimum(),
417+
height + extent.yMinimum()))
418+
bndpoints.append(QgsPointXY(width + extent.xMinimum(),
419+
height + extent.yMinimum()))
420+
elif YMaxNumber == 1:
421+
if hasXMin:
422+
bndpoints.append(QgsPointXY(extent.xMinimum(),
423+
height + extent.yMinimum()))
424+
elif hasXMax:
425+
bndpoints.append(QgsPointXY(width + extent.xMinimum(),
426+
height + extent.yMinimum()))
427+
# b) Polygon that covers the x or the y extent:
428+
if hasYMin and hasYMax:
429+
if YMaxNumber > 1:
430+
if hasXMin:
431+
bndpoints.append(QgsPointXY(extent.xMinimum(),
432+
extent.yMinimum()))
433+
elif hasXMax:
434+
bndpoints.append(QgsPointXY(width + extent.xMinimum(),
435+
extent.yMinimum()))
436+
elif YMinNumber > 1:
437+
if hasXMin:
438+
bndpoints.append(QgsPointXY(extent.xMinimum(),
439+
height + extent.yMinimum()))
440+
elif hasXMax:
441+
bndpoints.append(QgsPointXY(width + extent.xMinimum(),
442+
height + extent.yMinimum()))
443+
elif hasXMin and hasXMax:
444+
if XMaxNumber > 1:
445+
if hasYMin:
446+
bndpoints.append(QgsPointXY(extent.xMinimum(),
447+
extent.yMinimum()))
448+
elif hasYMax:
449+
bndpoints.append(QgsPointXY(width + extent.xMinimum(),
450+
extent.yMinimum()))
451+
elif XMinNumber > 1:
452+
if hasYMin:
453+
bndpoints.append(QgsPointXY(extent.xMinimum(),
454+
height + extent.yMinimum()))
455+
elif hasYMax:
456+
bndpoints.append(QgsPointXY(width + extent.xMinimum(),
457+
height + extent.yMinimum()))
458+
# c) Simple corners:
459+
if XMinNumber == 1 and YMinNumber == 1 and not hasXMax and not hasYMax:
460+
bndpoints.append(QgsPointXY(extent.xMinimum(),
289461
extent.yMinimum()))
290-
if hasXMax:
291-
if hasYMax:
292-
lines.append(QgsPointXY(width + extent.xMinimum(),
462+
if XMinNumber == 1 and YMaxNumber == 1 and not hasXMax and not hasYMin:
463+
bndpoints.append(QgsPointXY(extent.xMinimum(),
293464
height + extent.yMinimum()))
294-
if hasYMin:
295-
lines.append(QgsPointXY(width + extent.xMinimum(),
465+
if XMaxNumber == 1 and YMinNumber == 1 and not hasXMin and not hasYMax:
466+
bndpoints.append(QgsPointXY(width + extent.xMinimum(),
296467
extent.yMinimum()))
297-
return lines
468+
if XMaxNumber == 1 and YMaxNumber == 1 and not hasXMin and not hasYMin:
469+
bndpoints.append(QgsPointXY(width + extent.xMinimum(),
470+
height + extent.yMinimum()))
471+
return bndpoints

0 commit comments

Comments
 (0)
Please sign in to comment.