1# -*- coding: utf-8 -*-
2
3"""
4***************************************************************************
5    VoronoiPolygons.py
6    ---------------------
7    Date                 : August 2012
8    Copyright            : (C) 2012 by Victor Olaya
9    Email                : volayaf at gmail dot com
10***************************************************************************
11*                                                                         *
12*   This program is free software; you can redistribute it and/or modify  *
13*   it under the terms of the GNU General Public License as published by  *
14*   the Free Software Foundation; either version 2 of the License, or     *
15*   (at your option) any later version.                                   *
16*                                                                         *
17***************************************************************************
18"""
19
20__authors__ = 'Victor Olaya, Håvard Tveite'
21__date__ = 'August 2012'
22__copyright__ = '(C) 2012, Victor Olaya'
23
24import os
25
26from qgis.PyQt.QtGui import QIcon
27
28from qgis.core import (QgsApplication,
29                       QgsFeatureRequest,
30                       QgsFeatureSink,
31                       QgsFeature,
32                       QgsGeometry,
33                       QgsPointXY,
34                       QgsWkbTypes,
35                       QgsProcessing,
36                       QgsProcessingException,
37                       QgsProcessingParameterFeatureSource,
38                       QgsProcessingParameterFeatureSink,
39                       QgsProcessingParameterNumber)
40
41from processing.algs.qgis.QgisAlgorithm import QgisAlgorithm
42
43from . import voronoi
44
45pluginPath = os.path.split(os.path.split(os.path.dirname(__file__))[0])[0]
46
47
48class VoronoiPolygons(QgisAlgorithm):
49    INPUT = 'INPUT'
50    BUFFER = 'BUFFER'
51    OUTPUT = 'OUTPUT'
52
53    def icon(self):
54        return QgsApplication.getThemeIcon("/algorithms/mAlgorithmVoronoi.svg")
55
56    def svgIconPath(self):
57        return QgsApplication.iconPath("/algorithms/mAlgorithmVoronoi.svg")
58
59    def group(self):
60        return self.tr('Vector geometry')
61
62    def groupId(self):
63        return 'vectorgeometry'
64
65    def __init__(self):
66        super().__init__()
67
68    def initAlgorithm(self, config=None):
69        self.addParameter(
70            QgsProcessingParameterFeatureSource(
71                self.INPUT, self.tr('Input layer'),
72                [QgsProcessing.TypeVectorPoint]))
73        self.addParameter(
74            QgsProcessingParameterNumber(
75                self.BUFFER, self.tr('Buffer region (% of extent)'),
76                minValue=0.0, defaultValue=0.0))
77        self.addParameter(
78            QgsProcessingParameterFeatureSink(
79                self.OUTPUT, self.tr('Voronoi polygons'),
80                type=QgsProcessing.TypeVectorPolygon))
81
82    def name(self):
83        return 'voronoipolygons'
84
85    def displayName(self):
86        return self.tr('Voronoi polygons')
87
88    def tags(self):
89        return self.tr('thiessen').split(',')
90
91    def processAlgorithm(self, parameters, context, feedback):
92        source = self.parameterAsSource(parameters, self.INPUT, context)
93        if source is None:
94            raise QgsProcessingException(
95                self.invalidSourceError(parameters, self.INPUT))
96
97        buf = self.parameterAsDouble(parameters, self.BUFFER, context)
98        (sink, dest_id) = self.parameterAsSink(parameters, self.OUTPUT,
99                                               context, source.fields(),
100                                               QgsWkbTypes.Polygon,
101                                               source.sourceCrs())
102        if sink is None:
103            raise QgsProcessingException(
104                self.invalidSinkError(parameters, self.OUTPUT))
105        outFeat = QgsFeature()
106        extent = source.sourceExtent()
107        extraX = extent.width() * (buf / 100.0)
108        # Adjust the extent
109        extent.setXMinimum(extent.xMinimum() - extraX)
110        extent.setXMaximum(extent.xMaximum() + extraX)
111        extraY = extent.height() * (buf / 100.0)
112        extent.setYMinimum(extent.yMinimum() - extraY)
113        extent.setYMaximum(extent.yMaximum() + extraY)
114        height = extent.height()
115        width = extent.width()
116        c = voronoi.Context()
117        pts = []
118        ptDict = {}
119        ptNdx = -1
120        # Find the minimum and maximum x and y for the input points
121        xmin = width
122        xmax = 0
123        ymin = height
124        ymax = 0
125        features = source.getFeatures()
126        total = 100.0 / source.featureCount() if source.featureCount() else 0
127        for current, inFeat in enumerate(features):
128            if feedback.isCanceled():
129                break
130            geom = inFeat.geometry()
131            point = geom.asPoint()
132            x = point.x() - extent.xMinimum()
133            y = point.y() - extent.yMinimum()
134            pts.append((x, y))
135            ptNdx += 1
136            ptDict[ptNdx] = inFeat.id()
137            if x < xmin:
138                xmin = x
139            if y < ymin:
140                ymin = y
141            if x > xmax:
142                xmax = x
143            if y > ymax:
144                ymax = y
145            feedback.setProgress(int(current * total))
146        if xmin == xmax or ymin == ymax:
147            raise QgsProcessingException('The extent of the input points is '
148                                         'not a polygon (all the points are '
149                                         'on a vertical or horizontal line) '
150                                         '- cannot make a Voronoi diagram!')
151        xyminmax = [xmin, ymin, xmax, ymax]
152        if len(pts) < 3:
153            raise QgsProcessingException(
154                self.tr('Input file should contain at least 3 points. Choose '
155                        'another file and try again.'))
156        # Eliminate duplicate points
157        uniqueSet = set(item for item in pts)
158        ids = [pts.index(item) for item in uniqueSet]
159        sl = voronoi.SiteList([voronoi.Site(i[0], i[1], sitenum=j)
160                               for (j, i) in enumerate(uniqueSet)])
161        voronoi.voronoi(sl, c)
162        if len(c.polygons) == 0:
163            raise QgsProcessingException(
164                self.tr('There were no polygons created.'))
165
166        inFeat = QgsFeature()
167        current = 0
168        total = 100.0 / len(c.polygons)
169        # Clip each of the generated "polygons"
170        for (site, edges) in list(c.polygons.items()):
171            if feedback.isCanceled():
172                break
173            request = QgsFeatureRequest().setFilterFid(ptDict[ids[site]])
174            inFeat = next(source.getFeatures(request))
175            boundarypoints = self.clip_voronoi(edges, c, width,
176                                               height, extent,
177                                               inFeat.geometry().asPoint(),
178                                               xyminmax)
179            ptgeom = QgsGeometry.fromMultiPointXY(boundarypoints)
180            geom = QgsGeometry(ptgeom.convexHull())
181            outFeat.setGeometry(geom)
182            outFeat.setAttributes(inFeat.attributes())
183            sink.addFeature(outFeat, QgsFeatureSink.FastInsert)
184            current += 1
185            feedback.setProgress(int(current * total))
186        return {self.OUTPUT: dest_id}
187
188    def clip_voronoi(self, edges, c, width, height, extent, point, xyminmax):
189        """Clip voronoi function based on code written for Inkscape.
190        Copyright (C) 2010 Alvin Penner, penner@vaxxine.com
191        Clips one Thiessen polygon (convex polygon) to extent
192        """
193
194        pt_x = point.x() - extent.xMinimum()
195        pt_y = point.y() - extent.yMinimum()
196        (xmin, ymin, xmax, ymax) = xyminmax
197
198        def isclose(a, b, rel_tol=1e-9, abs_tol=0.0):
199            return abs(a - b) <= max(rel_tol * max(abs(a), abs(b)), abs_tol)
200
201        def clip_line(x1, y1, x2, y2, w, h):
202            if x1 < 0 and x2 < 0:
203                # Completely to the left
204                return [0, 0, 0, 0]
205            if x1 > w and x2 > w:
206                # Completely to the right
207                return [0, 0, 0, 0]
208            if y1 < 0 and y2 < 0:
209                # Completely below
210                return [0, 0, 0, 0]
211            if y1 > h and y2 > h:
212                # Completely above
213                return [0, 0, 0, 0]
214            # Clip on the left envelope boundary
215            if x1 < 0:
216                # First point to the left
217                y1 = (y1 * x2 - y2 * x1) / (x2 - x1)
218                x1 = 0
219            if x2 < 0:
220                # Last point to the left
221                y2 = (y1 * x2 - y2 * x1) / (x2 - x1)
222                x2 = 0
223            # Clip on the right envelope boundary
224            if x1 > w:
225                # First point to the right
226                y1 = y1 + (w - x1) * (y2 - y1) / (x2 - x1)
227                x1 = w
228            if x2 > w:
229                # Last point to the right
230                y2 = y1 + (w - x1) * (y2 - y1) / (x2 - x1)
231                x2 = w
232            if isclose(x1, x2) and isclose(y1, y2):
233                return [0, 0, 0, 0]
234            # Clip on the bottom envelope boundary
235            if y1 < 0:
236                # First point below
237                x1 = (x1 * y2 - x2 * y1) / (y2 - y1)
238                y1 = 0
239            if y2 < 0:
240                # Second point below
241                x2 = (x1 * y2 - x2 * y1) / (y2 - y1)
242                y2 = 0
243            # Clip on the top envelope boundary
244            if y1 > h:
245                # First point above
246                x1 = x1 + (h - y1) * (x2 - x1) / (y2 - y1)
247                y1 = h
248            if y2 > h:
249                # Second point above
250                x2 = x1 + (h - y1) * (x2 - x1) / (y2 - y1)
251                y2 = h
252            if isclose(x1, x2) and isclose(y1, y2):
253                return [0, 0, 0, 0]
254            return [x1, y1, x2, y2]
255
256        bndpoints = []
257        hasXMin = False
258        hasYMin = False
259        hasXMax = False
260        hasYMax = False
261        XMinNumber = 0
262        XMaxNumber = 0
263        YMinNumber = 0
264        YMaxNumber = 0
265        # The same line may appear twice for collinear input points,
266        # so have to remember which lines have contributed
267        XMinLine = -1
268        XMaxLine = -1
269        YMinLine = -1
270        YMaxLine = -1
271        for edge in edges:
272            if edge[1] >= 0 and edge[2] >= 0:
273                # Two vertices
274                [x1, y1, x2, y2] = clip_line(
275                    c.vertices[edge[1]][0],
276                    c.vertices[edge[1]][1],
277                    c.vertices[edge[2]][0],
278                    c.vertices[edge[2]][1],
279                    width,
280                    height
281                )
282            elif edge[1] >= 0:
283                # Only one (left) vertex
284                if c.lines[edge[0]][1] == 0:
285                    # Vertical line
286                    # xtemp = c.lines[edge[0]][2] / c.lines[edge[0]][0]
287                    xtemp = c.vertices[edge[1]][0]
288                    ytemp = 0 - 1
289                    # if c.vertices[edge[1]][1] > height / 2:
290                    #    ytemp = height
291                    # else:
292                    #    ytemp = 0
293                else:
294                    # Create an end of the line at the right edge - OK
295                    xtemp = width
296                    ytemp = (c.lines[edge[0]][2] - width *
297                             c.lines[edge[0]][0]) / c.lines[edge[0]][1]
298                [x1, y1, x2, y2] = clip_line(
299                    c.vertices[edge[1]][0],
300                    c.vertices[edge[1]][1],
301                    xtemp,
302                    ytemp,
303                    width,
304                    height
305                )
306            elif edge[2] >= 0:
307                # Only one (right) vertex
308                if c.lines[edge[0]][1] == 0:
309                    # Vertical line
310                    # xtemp = c.lines[edge[0]][2] / c.lines[edge[0]][0]
311                    xtemp = c.vertices[edge[2]][0]
312                    ytemp = height + 1
313                    # if c.vertices[edge[2]][1] > height / 2:
314                    #    ytemp = height
315                    # else:
316                    #    ytemp = 0.0
317                else:
318                    # End the line at the left edge - OK
319                    xtemp = 0.0
320                    ytemp = c.lines[edge[0]][2] / c.lines[edge[0]][1]
321                [x1, y1, x2, y2] = clip_line(
322                    xtemp,
323                    ytemp,
324                    c.vertices[edge[2]][0],
325                    c.vertices[edge[2]][1],
326                    width,
327                    height,
328                )
329            else:
330                # No vertex, only a line
331                if c.lines[edge[0]][1] == 0:
332                    # Vertical line - should not happen
333                    xtemp = c.lines[edge[0]][2] / c.lines[edge[0]][0]
334                    ytemp = 0.0
335                    xend = xtemp
336                    yend = height
337                else:
338                    # End the line at both edges - ???
339                    xtemp = 0.0
340                    ytemp = c.lines[edge[0]][2] / c.lines[edge[0]][1]
341                    xend = width
342                    yend = (c.lines[edge[0]][2] - width *
343                            c.lines[edge[0]][0]) / c.lines[edge[0]][1]
344                [x1, y1, x2, y2] = clip_line(
345                    xtemp,
346                    ytemp,
347                    xend,
348                    yend,
349                    width,
350                    height,
351                )
352            if x1 or x2 or y1 or y2:
353                bndpoints.append(QgsPointXY(x1 + extent.xMinimum(),
354                                            y1 + extent.yMinimum()))
355                bndpoints.append(QgsPointXY(x2 + extent.xMinimum(),
356                                            y2 + extent.yMinimum()))
357                if 0 in (x1, x2):
358                    hasXMin = True
359                    if XMinLine != edge[0]:
360                        XMinNumber = XMinNumber + 1
361                        XMinLine = edge[0]
362                if 0 in (y1, y2):
363                    hasYMin = True
364                    if YMinLine != edge[0]:
365                        YMinNumber = YMinNumber + 1
366                        YMinLine = edge[0]
367                if height in (y1, y2):
368                    hasYMax = True
369                    if YMaxLine != edge[0]:
370                        YMaxNumber = YMaxNumber + 1
371                        YMaxLine = edge[0]
372                if width in (x1, x2):
373                    hasXMax = True
374                    if XMaxLine != edge[0]:
375                        XMaxNumber = XMaxNumber + 1
376                        XMaxLine = edge[0]
377
378        # Add auxiliary points for corner cases, if necessary (duplicate
379        # points is not a problem - will be ignored later).
380        # a) Extreme input points (lowest, leftmost, rightmost, highest)
381        #    A point can be extreme on both axis
382        if pt_x == xmin:  # leftmost point
383            if XMinNumber == 0:
384                bndpoints.append(QgsPointXY(extent.xMinimum(),
385                                            extent.yMinimum()))
386                bndpoints.append(QgsPointXY(extent.xMinimum(),
387                                            height + extent.yMinimum()))
388            elif XMinNumber == 1:
389                if hasYMin:
390                    bndpoints.append(QgsPointXY(extent.xMinimum(),
391                                                extent.yMinimum()))
392                elif hasYMax:
393                    bndpoints.append(QgsPointXY(extent.xMinimum(),
394                                                height + extent.yMinimum()))
395        elif pt_x == xmax:  # rightmost point
396            if XMaxNumber == 0:
397                bndpoints.append(QgsPointXY(width + extent.xMinimum(),
398                                            extent.yMinimum()))
399                bndpoints.append(QgsPointXY(width + extent.xMinimum(),
400                                            height + extent.yMinimum()))
401            elif XMaxNumber == 1:
402                if hasYMin:
403                    bndpoints.append(QgsPointXY(width + extent.xMinimum(),
404                                                extent.yMinimum()))
405                elif hasYMax:
406                    bndpoints.append(QgsPointXY(width + extent.xMinimum(),
407                                                height + extent.yMinimum()))
408        if pt_y == ymin:  # lowest point
409            if YMinNumber == 0:
410                bndpoints.append(QgsPointXY(extent.xMinimum(),
411                                            extent.yMinimum()))
412                bndpoints.append(QgsPointXY(width + extent.xMinimum(),
413                                            extent.yMinimum()))
414            elif YMinNumber == 1:
415                if hasXMin:
416                    bndpoints.append(QgsPointXY(extent.xMinimum(),
417                                                extent.yMinimum()))
418                elif hasXMax:
419                    bndpoints.append(QgsPointXY(width + extent.xMinimum(),
420                                                extent.yMinimum()))
421        elif pt_y == ymax:  # highest point
422            if YMaxNumber == 0:
423                bndpoints.append(QgsPointXY(extent.xMinimum(),
424                                            height + extent.yMinimum()))
425                bndpoints.append(QgsPointXY(width + extent.xMinimum(),
426                                            height + extent.yMinimum()))
427            elif YMaxNumber == 1:
428                if hasXMin:
429                    bndpoints.append(QgsPointXY(extent.xMinimum(),
430                                                height + extent.yMinimum()))
431                elif hasXMax:
432                    bndpoints.append(QgsPointXY(width + extent.xMinimum(),
433                                                height + extent.yMinimum()))
434        # b) Polygon that covers the x or the y extent:
435        if hasYMin and hasYMax:
436            if YMaxNumber > 1:
437                if hasXMin:
438                    bndpoints.append(QgsPointXY(extent.xMinimum(),
439                                                extent.yMinimum()))
440                elif hasXMax:
441                    bndpoints.append(QgsPointXY(width + extent.xMinimum(),
442                                                extent.yMinimum()))
443            elif YMinNumber > 1:
444                if hasXMin:
445                    bndpoints.append(QgsPointXY(extent.xMinimum(),
446                                                height + extent.yMinimum()))
447                elif hasXMax:
448                    bndpoints.append(QgsPointXY(width + extent.xMinimum(),
449                                                height + extent.yMinimum()))
450        elif hasXMin and hasXMax:
451            if XMaxNumber > 1:
452                if hasYMin:
453                    bndpoints.append(QgsPointXY(extent.xMinimum(),
454                                                extent.yMinimum()))
455                elif hasYMax:
456                    bndpoints.append(QgsPointXY(width + extent.xMinimum(),
457                                                extent.yMinimum()))
458            elif XMinNumber > 1:
459                if hasYMin:
460                    bndpoints.append(QgsPointXY(extent.xMinimum(),
461                                                height + extent.yMinimum()))
462                elif hasYMax:
463                    bndpoints.append(QgsPointXY(width + extent.xMinimum(),
464                                                height + extent.yMinimum()))
465        # c) Simple corners:
466        if XMinNumber == 1 and YMinNumber == 1 and not hasXMax and not hasYMax:
467            bndpoints.append(QgsPointXY(extent.xMinimum(),
468                                        extent.yMinimum()))
469        if XMinNumber == 1 and YMaxNumber == 1 and not hasXMax and not hasYMin:
470            bndpoints.append(QgsPointXY(extent.xMinimum(),
471                                        height + extent.yMinimum()))
472        if XMaxNumber == 1 and YMinNumber == 1 and not hasXMin and not hasYMax:
473            bndpoints.append(QgsPointXY(width + extent.xMinimum(),
474                                        extent.yMinimum()))
475        if XMaxNumber == 1 and YMaxNumber == 1 and not hasXMin and not hasYMin:
476            bndpoints.append(QgsPointXY(width + extent.xMinimum(),
477                                        height + extent.yMinimum()))
478        return bndpoints
479