KNearestConcaveHull.py 21 KB
Newer Older
xuebingbing's avatar
xuebingbing committed
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428 429 430 431 432 433 434 435 436 437 438 439 440 441 442 443 444 445 446 447 448 449 450 451 452 453 454 455 456 457 458 459 460 461 462 463 464 465 466 467 468 469 470 471 472 473 474 475 476 477 478 479 480 481 482 483 484 485 486 487 488 489 490 491 492 493 494 495 496 497 498 499 500 501 502 503 504 505 506 507 508 509 510 511 512 513 514 515 516 517 518 519 520 521 522 523 524 525 526 527 528 529 530 531 532 533 534 535 536 537 538 539 540 541 542 543 544 545 546 547
# -*- coding: utf-8 -*-
"""
/***************************************************************************
    KNearestConcaveHull.py
    ----------------------
    Date                 : November 2014
    Copyright            : (C) 2014 by Detlev Neumann
                           Dr. Neumann Consulting - Geospatial Services
    Email                : dneumann@geospatial-services.de
 ***************************************************************************/

/***************************************************************************
 *                                                                         *
 *   This program is free software; you can redistribute it and/or modify  *
 *   it under the terms of the GNU General Public License as published by  *
 *   the Free Software Foundation; either version 2 of the License, or     *
 *   (at your option) any later version.                                   *
 *                                                                         *
 ***************************************************************************/
"""

__author__ = 'Detlev Neumann'
__date__ = 'November 2014'
__copyright__ = '(C) 2014, Detlev Neumann'

import os.path
import math

from qgis.PyQt.QtGui import QIcon
from qgis.PyQt.QtCore import QVariant

from qgis.core import (QgsApplication,
                       QgsExpression,
                       QgsFeature,
                       QgsFeatureRequest,
                       QgsFeatureSink,
                       QgsField,
                       QgsFields,
                       QgsGeometry,
                       QgsProcessing,
                       QgsProcessingException,
                       QgsProcessingParameterFeatureSink,
                       QgsProcessingParameterFeatureSource,
                       QgsProcessingParameterField,
                       QgsProcessingParameterNumber,
                       QgsPoint,
                       QgsPointXY,
                       QgsWkbTypes)
from processing.algs.qgis.QgisAlgorithm import QgisAlgorithm


class KNearestConcaveHull(QgisAlgorithm):
    KNEIGHBORS = 'KNEIGHBORS'
    INPUT = 'INPUT'
    OUTPUT = 'OUTPUT'
    FIELD = 'FIELD'

    def name(self):
        return 'knearestconcavehull'

    def displayName(self):
        return self.tr('Concave hull (k-nearest neighbor)')

    def shortDescription(self):
        return self.tr('Creates a concave hull using the k-nearest neighbor algorithm.')

    def icon(self):
        return QgsApplication.getThemeIcon("/algorithms/mAlgorithmConcaveHull.svg")

    def svgIconPath(self):
        return QgsApplication.iconPath("/algorithms/mAlgorithmConcaveHull.svg")

    def group(self):
        return self.tr('Vector geometry')

    def groupId(self):
        return 'vectorgeometry'

    def __init__(self):
        super().__init__()

    def initAlgorithm(self, config=None):
        self.addParameter(QgsProcessingParameterFeatureSource(self.INPUT,
                                                              self.tr('Input layer')))
        self.addParameter(QgsProcessingParameterNumber(self.KNEIGHBORS,
                                                       self.tr('Number of neighboring points to consider (a lower number is more concave, a higher number is smoother)'),
                                                       QgsProcessingParameterNumber.Integer,
                                                       defaultValue=3, minValue=3))
        self.addParameter(QgsProcessingParameterField(self.FIELD,
                                                      self.tr('Field (set if creating concave hulls by class)'),
                                                      parentLayerParameterName=self.INPUT, optional=True))
        self.addParameter(QgsProcessingParameterFeatureSink(self.OUTPUT, self.tr('Concave hull'),
                                                            QgsProcessing.TypeVectorPolygon))

    def processAlgorithm(self, parameters, context, feedback):
        # Get variables from dialog
        source = self.parameterAsSource(parameters, self.INPUT, context)
        if source is None:
            raise QgsProcessingException(self.invalidSourceError(parameters, self.INPUT))

        field_name = self.parameterAsString(parameters, self.FIELD, context)
        kneighbors = self.parameterAsInt(parameters, self.KNEIGHBORS, context)
        use_field = bool(field_name)

        field_index = -1

        fields = QgsFields()
        fields.append(QgsField('id', QVariant.Int, '', 20))

        current = 0

        # Get properties of the field the grouping is based on
        if use_field:
            field_index = source.fields().lookupField(field_name)
            if field_index >= 0:
                fields.append(source.fields()[field_index]) # Add a field with the name of the grouping field

                # Initialize writer
                (sink, dest_id) = self.parameterAsSink(parameters, self.OUTPUT, context,
                                                       fields, QgsWkbTypes.Polygon, source.sourceCrs())
                if sink is None:
                    raise QgsProcessingException(self.invalidSinkError(parameters, self.OUTPUT))

                success = False
                fid = 0

                # Get unique values of grouping field
                unique_values = source.uniqueValues(field_index)
                total = 100.0 / float(source.featureCount() * len(unique_values))

                for unique in unique_values:
                    points = []
                    filter = QgsExpression.createFieldEqualityExpression(field_name, unique)
                    request = QgsFeatureRequest().setFilterExpression(filter)
                    request.setSubsetOfAttributes([])
                    # Get features with the grouping attribute equal to the current grouping value
                    features = source.getFeatures(request)
                    for in_feature in features:
                        if feedback.isCanceled():
                            break
                        # Add points or vertices of more complex geometry
                        points.extend(extract_points(in_feature.geometry()))
                        current += 1
                        feedback.setProgress(int(current * total))
                    # A minimum of 3 points is necessary to proceed
                    if len(points) >= 3:
                        out_feature = QgsFeature()
                        the_hull = concave_hull(points, kneighbors)
                        if the_hull:
                            vertex = [QgsPointXY(point[0], point[1]) for point in the_hull]
                            poly = QgsGeometry().fromPolygonXY([vertex])

                            out_feature.setGeometry(poly)
                            # Give the polygon the same attribute as the point grouping attribute
                            out_feature.setAttributes([fid, unique])
                            sink.addFeature(out_feature, QgsFeatureSink.FastInsert)
                            success = True  # at least one polygon created
                    fid += 1
                if not success:
                    raise QgsProcessingException('No hulls could be created. Most likely there were not at least three unique points in any of the groups.')
            else:
                # Field parameter provided but can't read from it
                raise QgsProcessingException('Unable to find grouping field')

        else:
            # Not grouped by field
            # Initialize writer
            (sink, dest_id) = self.parameterAsSink(parameters, self.OUTPUT, context,
                                                   fields, QgsWkbTypes.Polygon, source.sourceCrs())
            if sink is None:
                raise QgsProcessingException(self.invalidSinkError(parameters, self.OUTPUT))

            points = []
            request = QgsFeatureRequest()
            request.setSubsetOfAttributes([])
            features = source.getFeatures(request) # Get all features
            total = 100.0 / source.featureCount() if source.featureCount() else 0
            for in_feature in features:
                if feedback.isCanceled():
                    break
                # Add points or vertices of more complex geometry
                points.extend(extract_points(in_feature.geometry()))
                current += 1
                feedback.setProgress(int(current * total))

            # A minimum of 3 points is necessary to proceed
            if len(points) >= 3:
                out_feature = QgsFeature()
                the_hull = concave_hull(points, kneighbors)
                if the_hull:
                    vertex = [QgsPointXY(point[0], point[1]) for point in the_hull]
                    poly = QgsGeometry().fromPolygonXY([vertex])

                    out_feature.setGeometry(poly)
                    out_feature.setAttributes([0])
                    sink.addFeature(out_feature, QgsFeatureSink.FastInsert)
                else:
                    # the_hull returns None only when there are less than three points after cleaning
                    raise QgsProcessingException('At least three unique points are required to create a concave hull.')
            else:
                raise QgsProcessingException('At least three points are required to create a concave hull.')

        return {self.OUTPUT: dest_id}


def clean_list(list_of_points):
    """
    Deletes duplicate points in list_of_points
    """
    return list(set(list_of_points))


def find_min_y_point(list_of_points):
    """
    Returns that point of *list_of_points* having minimal y-coordinate

    :param list_of_points: list of tuples
    :return: tuple (x, y)
    """
    min_y_pt = list_of_points[0]
    for point in list_of_points[1:]:
        if point[1] < min_y_pt[1] or (point[1] == min_y_pt[1] and point[0] < min_y_pt[0]):
            min_y_pt = point
    return min_y_pt


def add_point(vector, element):
    """
    Returns vector with the given element append to the right
    """
    vector.append(element)
    return vector


def remove_point(vector, element):
    """
    Returns a copy of vector without the given element
    """
    vector.pop(vector.index(element))
    return vector


def euclidian_distance(point1, point2):
    """
    Returns the euclidian distance of the 2 given points.

    :param point1: tuple (x, y)
    :param point2: tuple (x, y)
    :return: float
    """
    return math.sqrt(math.pow(point1[0] - point2[0], 2) + math.pow(point1[1] - point2[1], 2))


def nearest_points(list_of_points, point, k):
    """
    Returns a list of the indices of the k closest neighbors from list_of_points to the specified point. The measure
    of proximity is the Euclidean distance. Internally, k becomes the minimum between the given value for k and the
    number of points in list_of_points

    :param list_of_points: list of tuples
    :param point: tuple (x, y)
    :param k: integer
    :return: list of k tuples
    """
    # build a list of tuples of distances between point *point* and every point in *list_of_points*, and
    # their respective index of list *list_of_distances*
    list_of_distances = []
    for index in range(len(list_of_points)):
        list_of_distances.append((euclidian_distance(list_of_points[index], point), index))

    # sort distances in ascending order
    list_of_distances.sort()

    # get the k nearest neighbors of point
    nearest_list = []
    for index in range(min(k, len(list_of_points))):
        nearest_list.append((list_of_points[list_of_distances[index][1]]))
    return nearest_list


def angle(from_point, to_point):
    """
    Returns the angle of the directed line segment, going from *from_point* to *to_point*, in radians. The angle is
    positive for segments with upward direction (north), otherwise negative (south). Values ranges from 0 at the
    right (east) to pi at the left side (west).

    :param from_point: tuple (x, y)
    :param to_point: tuple (x, y)
    :return: float
    """
    return math.atan2(to_point[1] - from_point[1], to_point[0] - from_point[0])


def angle_difference(angle1, angle2):
    """
    Calculates the difference between the given angles in clockwise direction as radians.

    :param angle1: float
    :param angle2: float
    :return: float; between 0 and 2*Pi
    """
    if (angle1 > 0 and angle2 >= 0) and angle1 > angle2:
        return abs(angle1 - angle2)
    elif (angle1 >= 0 and angle2 > 0) and angle1 < angle2:
        return 2 * math.pi + angle1 - angle2
    elif (angle1 < 0 and angle2 <= 0) and angle1 < angle2:
        return 2 * math.pi + angle1 + abs(angle2)
    elif (angle1 <= 0 and angle2 < 0) and angle1 > angle2:
        return abs(angle1 - angle2)
    elif angle1 <= 0 < angle2:
        return 2 * math.pi + angle1 - angle2
    elif angle1 >= 0 >= angle2:
        return angle1 + abs(angle2)
    else:
        return 0


def intersect(line1, line2):
    """
    Returns True if the two given line segments intersect each other, and False otherwise.

    :param line1: 2-tuple of tuple (x, y)
    :param line2: 2-tuple of tuple (x, y)
    :return: boolean
    """
    a1 = line1[1][1] - line1[0][1]
    b1 = line1[0][0] - line1[1][0]
    c1 = a1 * line1[0][0] + b1 * line1[0][1]
    a2 = line2[1][1] - line2[0][1]
    b2 = line2[0][0] - line2[1][0]
    c2 = a2 * line2[0][0] + b2 * line2[0][1]
    tmp = (a1 * b2 - a2 * b1)
    if tmp == 0:
        return False
    sx = (c1 * b2 - c2 * b1) / tmp
    if (sx > line1[0][0] and sx > line1[1][0]) or (sx > line2[0][0] and sx > line2[1][0]) or\
            (sx < line1[0][0] and sx < line1[1][0]) or (sx < line2[0][0] and sx < line2[1][0]):
        return False
    sy = (a1 * c2 - a2 * c1) / tmp
    if (sy > line1[0][1] and sy > line1[1][1]) or (sy > line2[0][1] and sy > line2[1][1]) or\
            (sy < line1[0][1] and sy < line1[1][1]) or (sy < line2[0][1] and sy < line2[1][1]):
        return False
    return True


def point_in_polygon_q(point, list_of_points):
    """
    Return True if given point *point* is laying in the polygon described by the vertices *list_of_points*,
    otherwise False

    Based on the "Ray Casting Method" described by Joel Lawhead in this blog article:
    http://geospatialpython.com/2011/01/point-in-polygon.html

    """
    x = point[0]
    y = point[1]
    poly = [(pt[0], pt[1]) for pt in list_of_points]
    n = len(poly)
    inside = False

    p1x, p1y = poly[0]
    for i in range(n + 1):
        p2x, p2y = poly[i % n]
        if y > min(p1y, p2y):
            if y <= max(p1y, p2y):
                if x <= max(p1x, p2x):
                    if p1y != p2y:
                        xints = (y - p1y) * (p2x - p1x) / (p2y - p1y) + p1x
                    if p1x == p2x or x <= xints:
                        inside = not inside
        p1x, p1y = p2x, p2y

    return inside


def extract_points(geom):
    """
    Generate list of QgsPoints from QgsGeometry *geom* ( can be point, line, or polygon )
    Code taken from fTools plugin

    :param geom: an arbitrary geometry feature
    :return: list of points
    """
    multi_geom = QgsGeometry()
    temp_geom = []
    # point geometry
    if geom.type() == 0:
        if geom.isMultipart():
            temp_geom = geom.asMultiPoint()
        else:
            temp_geom.append(geom.asPoint())
    # line geometry
    if geom.type() == 1:
        # if multipart feature explode to single part
        if geom.isMultipart():
            multi_geom = geom.asMultiPolyline()
            for i in multi_geom:
                temp_geom.extend(i)
        else:
            temp_geom = geom.asPolyline()
    # polygon geometry
    elif geom.type() == 2:
        # if multipart feature explode to single part
        if geom.isMultipart():
            multi_geom = geom.asMultiPolygon()
            # now single part polygons
            for i in multi_geom:
                # explode to line segments
                for j in i:
                    temp_geom.extend(j)
        else:
            multi_geom = geom.asPolygon()
            # explode to line segments
            for i in multi_geom:
                temp_geom.extend(i)
    return temp_geom


def sort_by_angle(list_of_points, last_point, last_angle):
    """
    returns the points in list_of_points in descending order of angle to the last segment of the envelope, measured
    in a clockwise direction. Thus, the rightmost of the neighboring points is always selected. The first point of
    this list will be the next point of the envelope.
    """
    def getkey(item):
        return angle_difference(last_angle, angle(last_point, item))

    vertex_list = sorted(list_of_points, key=getkey, reverse=True)
    return vertex_list


def concave_hull(points_list, k):
    """
    Calculates a valid concave hull polygon containing all given points. The algorithm searches for that
    point in the neighborhood of k nearest neighbors which maximizes the rotation angle in clockwise direction
    without intersecting any previous line segments.

    This is an implementation of the algorithm described by Adriano Moreira and Maribel Yasmina Santos:
    CONCAVE HULL: A neighborhood_k-NEAREST NEIGHBORS APPROACH FOR THE COMPUTATION OF THE REGION OCCUPIED BY A SET OF POINTS.
    GRAPP 2007 - International Conference on Computer Graphics Theory and Applications; pp 61-68.

    :param points_list: list of tuples (x, y)
    :param k: integer
    :return: list of tuples (x, y)
    """
    # return an empty list if not enough points are given
    if k > len(points_list):
        k = len(points_list)

    # the number of nearest neighbors k must be greater than or equal to 3
    kk = max(k, 3)

    # delete duplicate points
    point_set = clean_list(points_list)

    # if point_set has less then 3 points no polygon can be created and an empty list will be returned
    if len(point_set) < 3:
        return None

    # if point_set has 3 points then these are already vertices of the hull. Append the first point to
    # close the hull polygon
    if len(point_set) == 3:
        return add_point(point_set, point_set[0])

    # make sure that k neighbors can be found
    kk = min(kk, len(point_set))

    # start with the point having the smallest y-coordinate (most southern point)
    first_point = find_min_y_point(point_set)

    # add this points as the first vertex of the hull
    hull = [first_point]

    # make the first vertex of the hull to the current point
    current_point = first_point

    # remove the point from the point_set, to prevent him being among the nearest points
    point_set = remove_point(point_set, first_point)
    previous_angle = math.pi

    # step counts the number of segments
    step = 2

    # as long as point_set is not empty or search is returning to the starting point
    while (current_point != first_point) or (step == 2) and (len(point_set) > 0):

        # after 3 iterations add the first point to point_set again, otherwise a hull cannot be closed
        if step == 5:
            point_set = add_point(point_set, first_point)

        # search the k nearest neighbors of the current point
        k_nearest_points = nearest_points(point_set, current_point, kk)

        # sort the candidates (neighbors) in descending order of right-hand turn. This way the algorithm progresses
        # in clockwise direction through as many points as possible
        c_points = sort_by_angle(k_nearest_points, current_point, previous_angle)

        its = True
        i = -1

        # search for the nearest point to which the connecting line does not intersect any existing segment
        while its is True and (i < len(c_points) - 1):
            i += 1
            if c_points[i] == first_point:
                last_point = 1
            else:
                last_point = 0
            j = 2
            its = False

            while its is False and (j < len(hull) - last_point):
                its = intersect((hull[step - 2], c_points[i]), (hull[step - 2 - j], hull[step - 1 - j]))
                j += 1

        # there is no candidate to which the connecting line does not intersect any existing segment, so the
        # for the next candidate fails. The algorithm starts again with an increased number of neighbors
        if its is True:
            return concave_hull(points_list, kk + 1)

        # the first point which complies with the requirements is added to the hull and gets the current point
        current_point = c_points[i]
        hull = add_point(hull, current_point)

        # calculate the angle between the last vertex and his precursor, that is the last segment of the hull
        # in reversed direction
        previous_angle = angle(hull[step - 1], hull[step - 2])

        # remove current_point from point_set
        point_set = remove_point(point_set, current_point)

        # increment counter
        step += 1

    all_inside = True
    i = len(point_set) - 1

    # check if all points are within the created polygon
    while (all_inside is True) and (i >= 0):
        all_inside = point_in_polygon_q(point_set[i], hull)
        i -= 1

    # since at least one point is out of the computed polygon, try again with a higher number of neighbors
    if all_inside is False:
        return concave_hull(points_list, kk + 1)

    # a valid hull has been constructed
    return hull