From 302f0fe068a819ff85fad5923b8e7be521f1af0b Mon Sep 17 00:00:00 2001 From: Fabian Wunsch Date: Mon, 11 Apr 2022 16:28:30 +0200 Subject: [PATCH] Changed codebase and minor bugfixes Changed the codebase to https://github.com/Twista/python-polylabel/ and rounded results --- tools/calculate_center.py | 256 +++++++++++++++++++++----------------- 1 file changed, 139 insertions(+), 117 deletions(-) diff --git a/tools/calculate_center.py b/tools/calculate_center.py index cafa041c..ff1bcd8e 100644 --- a/tools/calculate_center.py +++ b/tools/calculate_center.py @@ -1,131 +1,35 @@ """ -From https://github.com/stefda/polylabel, +From https://github.com/Twista/python-polylabel/, which is in turn implemented from https://github.com/mapbox/polylabel """ -from queue import Queue -import math +from math import sqrt +import time -SQRT2 = math.sqrt(2) +# Python3 +from queue import PriorityQueue +from math import inf -def polylabel(polygon, precision=0.5, debug=False): - """ - Computes the pole of inaccessibility coordinate in [x, y] format. - """ - - # find the bounding box of the outer ring - minX = maxX = polygon[0][0] - minY = maxY = polygon[0][1] - for p in polygon[1:]: - if p[0] < minX: minX = p[0] - if p[1] < minY: minY = p[1] - if p[0] > maxX: maxX = p[0] - if p[1] > maxY: maxY = p[1] - - width = maxX - minX - height = maxY - minY - cellSize = min(width, height) - h = cellSize / 2.0 - - # priority queue of cells in order of their "potential" (max distance to polygon) - cellQueue = Queue(0) - - if cellSize == 0: return [minX, minY]; - - # cover polygon with initial cells - for x in range(math.floor(minX), math.ceil(maxX - 1), round(cellSize)): - for y in range(math.floor(minY), math.ceil(maxY - 1), round(cellSize)): - cellQueue.put(Cell(x + h, y + h, h, polygon)) - - # take centroid as the first best guess - bestCell = getCentroidCell(polygon) - - # special case for rectangular polygons - bboxCell = Cell(minX + width / 2, minY + height / 2, 0, polygon) - if bboxCell.d > bestCell.d: bestCell = bboxCell - - numProbes = cellQueue.qsize() - - while not cellQueue.empty(): - # pick the most promising cell from the queue - cell = cellQueue.get() - - # update the best cell if we found a better one - if cell.d > bestCell.d: - bestCell = cell - if debug: print('found best %d after %d probes' % (round(1e4 * cell.d) / 1e4, numProbes)) - - # do not drill down further if there's no chance of a better solution - if cell.max - bestCell.d <= precision: continue - - # split the cell into four cells - h = cell.h / 2 - cellQueue.put(Cell(cell.x - h, cell.y - h, h, polygon)) - cellQueue.put(Cell(cell.x + h, cell.y - h, h, polygon)) - cellQueue.put(Cell(cell.x - h, cell.y + h, h, polygon)) - cellQueue.put(Cell(cell.x + h, cell.y + h, h, polygon)) - numProbes += 4 - - if debug: - print('num probes: %d' % numProbes) - print('best distance: %d' % bestCell.d) - - return [bestCell.x, bestCell.y] - - -def compareMax(a, b): - return b.max - a.max - - -class Cell: - def __init__(self, x, y, h, polygon): - self.x = x # cell center x - self.y = y # cell center y - self.h = h # half the cell size - self.d = pointToPolygonDist(x, y, polygon) # distance from cell center to polygon - self.max = self.d + self.h * SQRT2 # max distance to polygon within a cell - - -def pointToPolygonDist(x, y, polygon): - """ - Computes the signed distance from point to polygon outline (negative if point is outside). - """ +def _point_to_polygon_distance(x, y, polygon): inside = False - minDistSq = float('inf') + min_dist_sq = inf - for a, b in zip(polygon, rotate(polygon)): - if ((a[1] > y) != (b[1] > y)) and (x < (b[0] - a[0]) * (y - a[1]) / (b[1] - a[1]) + a[0]): + b = polygon[-1] + for a in polygon: + if ((a[1] > y) != (b[1] > y) and + (x < (b[0] - a[0]) * (y - a[1]) / (b[1] - a[1]) + a[0])): inside = not inside - minDistSq = min(minDistSq, getSegDistSq(x, y, a, b)) + min_dist_sq = min(min_dist_sq, _get_seg_dist_sq(x, y, a, b)) + b = a - return math.sqrt(minDistSq) * (1 if inside else -1) + result = sqrt(min_dist_sq) + if not inside: + return -result + return result -def getCentroidCell(polygon): - """ - Gets the polygon centroid. - """ - area = 0 - x = 0 - y = 0 - - for a, b in zip(polygon, rotate(polygon)): - f = a[0] * b[1] - b[0] * a[1] - x += (a[0] + b[0]) * f - y += (a[1] + b[1]) * f - area += f * 3 - - if area == 0: - return Cell(polygon[0][0], polygon[0][1], 0, polygon) - - return Cell(x / area, y / area, 0, polygon) - - -def getSegDistSq(px, py, a, b): - """ - Gets the squared distance from a point to a segment. - """ +def _get_seg_dist_sq(px, py, a, b): x = a[0] y = a[1] dx = b[0] - x @@ -133,9 +37,11 @@ def getSegDistSq(px, py, a, b): if dx != 0 or dy != 0: t = ((px - x) * dx + (py - y) * dy) / (dx * dx + dy * dy) + if t > 1: x = b[0] y = b[1] + elif t > 0: x += dx * t y += dy * t @@ -146,5 +52,121 @@ def getSegDistSq(px, py, a, b): return dx * dx + dy * dy -def rotate(l): - return l[-1:] + l[:-1] +class Cell(object): + def __init__(self, x, y, h, polygon): + self.h = h + self.y = y + self.x = x + self.d = _point_to_polygon_distance(x, y, polygon) + self.max = self.d + self.h * sqrt(2) + + def __lt__(self, other): + return self.max < other.max + + def __lte__(self, other): + return self.max <= other.max + + def __gt__(self, other): + return self.max > other.max + + def __gte__(self, other): + return self.max >= other.max + + def __eq__(self, other): + return self.max == other.max + + +def _get_centroid_cell(polygon): + area = 0 + x = 0 + y = 0 + previous = polygon[-1] + for current in polygon: + f = current[0] * previous[1] - previous[0] * current[1] + x += (current[0] + previous[0]) * f + y += (current[1] + previous[1]) * f + area += f * 3 + previous =current + if area == 0: + return Cell(polygon[0][0], polygon[0][1], 0, polygon) + return Cell(x / area, y / area, 0, polygon) + + +def polylabel(polygon, precision=0.5, debug=False, with_distance=False): + # find bounding box + first_item = polygon[0] + min_x = first_item[0] + min_y = first_item[1] + max_x = first_item[0] + max_y = first_item[1] + for p in polygon: + if p[0] < min_x: + min_x = p[0] + if p[1] < min_y: + min_y = p[1] + if p[0] > max_x: + max_x = p[0] + if p[1] > max_y: + max_y = p[1] + + width = max_x - min_x + height = max_y - min_y + cell_size = min(width, height) + h = cell_size / 2.0 + + cell_queue = PriorityQueue() + + if cell_size == 0: + if with_distance: + return [min_x, min_y], None + else: + return [min_x, min_y] + + # cover polygon with initial cells + x = min_x + while x < max_x: + y = min_y + while y < max_y: + c = Cell(x + h, y + h, h, polygon) + y += cell_size + cell_queue.put((-c.max, time.time(), c)) + x += cell_size + + best_cell = _get_centroid_cell(polygon) + + bbox_cell = Cell(min_x + width / 2, min_y + height / 2, 0, polygon) + if bbox_cell.d > best_cell.d: + best_cell = bbox_cell + + num_of_probes = cell_queue.qsize() + while not cell_queue.empty(): + _, __, cell = cell_queue.get() + + if cell.d > best_cell.d: + best_cell = cell + + if debug: + print('found best {} after {} probes'.format( + round(1e4 * cell.d) / 1e4, num_of_probes)) + + if cell.max - best_cell.d <= precision: + continue + + h = cell.h / 2 + c = Cell(cell.x - h, cell.y - h, h, polygon) + cell_queue.put((-c.max, time.time(), c)) + c = Cell(cell.x + h, cell.y - h, h, polygon) + cell_queue.put((-c.max, time.time(), c)) + c = Cell(cell.x - h, cell.y + h, h, polygon) + cell_queue.put((-c.max, time.time(), c)) + c = Cell(cell.x + h, cell.y + h, h, polygon) + cell_queue.put((-c.max, time.time(), c)) + num_of_probes += 4 + + if debug: + print('num probes: {}'.format(num_of_probes)) + print('best distance: {}'.format(best_cell.d)) + if with_distance: + return [best_cell.x, best_cell.y], best_cell.d + else: + return [best_cell.x, best_cell.y]