From e027bcaa8f7b6156ea9bc355e7e288949900ba50 Mon Sep 17 00:00:00 2001 From: Daan Bloembergen Date: Wed, 8 Feb 2023 19:01:44 +0100 Subject: [PATCH] Add alpha shape implementation (#60) * Alpha shape implementation. Fixes #50 * Fix bug where an empty BGT reference file messes up the labelling. * Bump version to 0.2. --- README.md | 2 +- environment.yml | 2 +- requirements.txt | 2 +- setup.cfg | 4 +- src/upcp/fusion/ahn_fuser.py | 26 ++- src/upcp/fusion/building_fuser.py | 2 +- src/upcp/fusion/car_fuser.py | 2 +- src/upcp/fusion/pole_fuser.py | 2 +- src/upcp/fusion/road_fuser.py | 2 +- src/upcp/fusion/street_furniture_fuser.py | 2 +- src/upcp/utils/alpha_shape_utils.py | 192 ++++++++++++++++++++++ src/upcp/utils/math_utils.py | 3 +- 12 files changed, 222 insertions(+), 19 deletions(-) create mode 100644 src/upcp/utils/alpha_shape_utils.py diff --git a/README.md b/README.md index 8fcf863..a862c80 100644 --- a/README.md +++ b/README.md @@ -73,7 +73,7 @@ This code has been tested with `Python >= 3.8` on `Linux` and `MacOS`, and shoul ```bash # Install the latest release as Wheel - python -m pip install https://github.com/Amsterdam-AI-Team/Urban_PointCloud_Processing/releases/download/v0.1/upcp-0.1-py3-none-any.whl + python -m pip install https://github.com/Amsterdam-AI-Team/Urban_PointCloud_Processing/releases/download/v0.2/upcp-0.2-py3-none-any.whl # Alternatively, install the latest version from source python -m pip install git+https://github.com/Amsterdam-AI-Team/Urban_PointCloud_Processing.git#egg=upcp diff --git a/environment.yml b/environment.yml index 7a74c93..0c8e9a8 100644 --- a/environment.yml +++ b/environment.yml @@ -15,7 +15,7 @@ dependencies: - owslib>=0.27 - pandas>=1.3 - pyntcloud>=0.3.1 - - shapely>=1.7 + - shapely>=1.8 - scikit-learn>=0.24 - scipy>=1.6 - tifffile>=2021.6 diff --git a/requirements.txt b/requirements.txt index 20e97d1..927ce8c 100644 --- a/requirements.txt +++ b/requirements.txt @@ -10,7 +10,7 @@ open3d~=0.16.0 owslib>=0.27 pandas>=1.3 pyntcloud>=0.3.1 -shapely>=1.7 +shapely>=1.8 scikit-learn>=0.24 scipy>=1.6 tifffile>=2021.6 diff --git a/setup.cfg b/setup.cfg index 643bb3f..e045acf 100644 --- a/setup.cfg +++ b/setup.cfg @@ -1,6 +1,6 @@ [metadata] name = upcp -version = 0.11 +version = 0.2 author = Daan Bloembergen author_email = daanbl.dev@gmail.com description = Repository for automatic classification and labeling of Urban PointClouds using data fusion and region growing techniques. @@ -33,7 +33,7 @@ install_requires = owslib>=0.27 pandas>=1.3 pyntcloud>=0.3.1 - shapely>=1.7 + shapely>=1.8 scikit-learn>=0.24 scipy>=1.6 tifffile>=2021.6 diff --git a/src/upcp/fusion/ahn_fuser.py b/src/upcp/fusion/ahn_fuser.py index cdda8e7..b194aaf 100644 --- a/src/upcp/fusion/ahn_fuser.py +++ b/src/upcp/fusion/ahn_fuser.py @@ -10,6 +10,7 @@ from ..labels import Labels from ..utils import clip_utils from ..utils import math_utils +from ..utils.alpha_shape_utils import alpha_shape logger = logging.getLogger(__name__) @@ -66,6 +67,10 @@ def _set_defaults(self, params): params['min_comp_size'] = 50 if 'buffer' not in params: params['buffer'] = 0.05 + if 'use_concave' not in params: + params['use_concave'] = True + if 'concave_alpha' not in params: + params['concave_alpha'] = 2.0 return params def _refine_layer(self, points, points_z, @@ -87,13 +92,20 @@ def _refine_layer(self, points, points_z, for cc in cc_labels: # select points that belong to the cluster cc_mask = (point_components == cc) - # Compute convex hull and add a small buffer - poly = (math_utils - .convex_hull_poly(points[label_ids[cc_mask], 0:2]) - .buffer(self.params['buffer'])) - # Select ground points within poly - poly_mask = clip_utils.poly_clip(points[ground_mask], poly) - mask[ground_ids[poly_mask]] = True + if not self.params['use_concave']: + # Compute convex hull + polys = [math_utils.convex_hull_poly( + points[label_ids[cc_mask], 0:2])] + else: + # Compute concave hull + polys = alpha_shape(points[label_ids[cc_mask], 0:2], + alpha=self.params['concave_alpha']) + # Select ground points within poly(s) + for poly in polys: + poly_mask = clip_utils.poly_clip( + points[ground_mask], + poly.buffer(self.params['buffer'])) + mask[ground_ids[poly_mask]] = True return mask def _refine_ground(self, points, points_z, ground_mask, diff --git a/src/upcp/fusion/building_fuser.py b/src/upcp/fusion/building_fuser.py index 86c0a2e..877daf2 100644 --- a/src/upcp/fusion/building_fuser.py +++ b/src/upcp/fusion/building_fuser.py @@ -73,7 +73,7 @@ def get_labels(self, points, labels, mask, tilecode): merge=True) if len(building_polygons) == 0: logger.debug('No buildings found for tile, skipping.') - return label_mask + return labels if mask is None: mask = np.ones((len(points),), dtype=bool) diff --git a/src/upcp/fusion/car_fuser.py b/src/upcp/fusion/car_fuser.py index 4d6b1f9..fde8571 100644 --- a/src/upcp/fusion/car_fuser.py +++ b/src/upcp/fusion/car_fuser.py @@ -114,7 +114,7 @@ def get_labels(self, points, labels, mask, tilecode): tilecode, exclude_types=self.exclude_types, merge=False) if len(road_polygons) == 0: logger.debug('No road parts found for tile, skipping.') - return label_mask + return labels # Get the interpolated ground points of the tile ground_z = self.ahn_reader.interpolate( diff --git a/src/upcp/fusion/pole_fuser.py b/src/upcp/fusion/pole_fuser.py index f95d9b7..c154b5f 100644 --- a/src/upcp/fusion/pole_fuser.py +++ b/src/upcp/fusion/pole_fuser.py @@ -254,7 +254,7 @@ def get_labels(self, points, labels, mask, tilecode): if len(bgt_points) == 0: logger.debug(f'No {self.bgt_type} objects found in tile, ' + ' skipping.') - return label_mask + return labels ahn_tile = self.ahn_reader.filter_tile(tilecode) fast_z = FastGridInterpolator(ahn_tile['x'], ahn_tile['y'], diff --git a/src/upcp/fusion/road_fuser.py b/src/upcp/fusion/road_fuser.py index 4b73a56..73353d5 100644 --- a/src/upcp/fusion/road_fuser.py +++ b/src/upcp/fusion/road_fuser.py @@ -75,7 +75,7 @@ def get_labels(self, points, labels, mask, tilecode): merge=True) if len(road_polygons) == 0: logger.debug('No road parts found in tile, skipping.') - return label_mask + return labels # Already labelled ground points can be labelled as road. mask = labels == Labels.GROUND diff --git a/src/upcp/fusion/street_furniture_fuser.py b/src/upcp/fusion/street_furniture_fuser.py index 3bd1e49..1a25420 100644 --- a/src/upcp/fusion/street_furniture_fuser.py +++ b/src/upcp/fusion/street_furniture_fuser.py @@ -117,7 +117,7 @@ def get_labels(self, points, labels, mask, tilecode): if len(bgt_points) == 0: logger.debug(f'No {self.bgt_type} objects found in tile, ' + ' skipping.') - return label_mask + return labels # Get the interpolated ground points of the tile ground_z = self.ahn_reader.interpolate( diff --git a/src/upcp/utils/alpha_shape_utils.py b/src/upcp/utils/alpha_shape_utils.py new file mode 100644 index 0000000..227d512 --- /dev/null +++ b/src/upcp/utils/alpha_shape_utils.py @@ -0,0 +1,192 @@ +# Urban_PointCloud_Processing by Amsterdam Intelligence, GPL-3.0 license + +"""This module provides methods to generate a concave hull (alpha shape).""" + +import numpy as np +import shapely.geometry as sg +from scipy.spatial import Delaunay + + +def alpha_shape(points, alpha=1.): + """ + Return the concave polygon (alpha shape) of a set of points. Depending on + the value of alpha, multiple polygons may be returned, and individual + polygons may contain holes (interior rings). + + Parameters + ---------- + points : np.array of shape (n,2) + Points of which concave hull should be returned. + alpha : float (default: 1.0) + Alpha value, determines "concaveness". A value of 0 equals the convex + hull. + + The derivation is as follows: to generate the concave hull, + "circumcircles" are fitted to the points by triangulation. Whenever the + radius of such a circle is larger than 1/alpha, a whole is generated. + + Returns + ------- + A list of shapely.Polygon objects representing the concave hull. + """ + edges = get_alpha_shape_edges(points, alpha=alpha, only_outer=True) + concave_polys = generate_poly_from_edges(edges, points) + return concave_polys + + +# https://stackoverflow.com/a/50159452 +# CC BY-SA 4.0 +def get_alpha_shape_edges(points, alpha, only_outer=True): + """ + Compute the alpha shape (concave hull) of a set of points. + :param points: np.array of shape (n,2) points. + :param alpha: alpha value. + :param only_outer: boolean value to specify if we keep only the outer + border or also inner edges. + :return: set of (i,j) pairs representing edges of the alpha-shape. (i,j) + are the indices in the points array. + """ + assert points.shape[0] > 3, "Need at least four points" + + def add_edge(edges, i, j): + """ + Add a line between the i-th and j-th points, + if not in the list already + """ + if (i, j) in edges or (j, i) in edges: + # already added + assert (j, i) in edges, ( + "Can't go twice over same directed edge right?") + if only_outer: + # if both neighboring triangles are in shape, it is not a + # boundary edge + edges.remove((j, i)) + return + edges.add((i, j)) + + tri = Delaunay(points) + edges = set() + # Loop over triangles: + # ia, ib, ic = indices of corner points of the triangle + for ia, ib, ic in tri.simplices: + pa = points[ia] + pb = points[ib] + pc = points[ic] + # Computing radius of triangle circumcircle + # www.mathalino.com/reviewer/derivation-of-formulas/derivation-of-formula-for-radius-of-circumcircle + a = np.sqrt((pa[0] - pb[0]) ** 2 + (pa[1] - pb[1]) ** 2) + b = np.sqrt((pb[0] - pc[0]) ** 2 + (pb[1] - pc[1]) ** 2) + c = np.sqrt((pc[0] - pa[0]) ** 2 + (pc[1] - pa[1]) ** 2) + s = (a + b + c) / 2.0 + area = np.sqrt(s * (s - a) * (s - b) * (s - c)) + if area == 0: + # TODO something more clever? + continue + circum_r = a * b * c / (4.0 * area) + if circum_r < (1 / alpha): + add_edge(edges, ia, ib) + add_edge(edges, ib, ic) + add_edge(edges, ic, ia) + return edges + + +# https://stackoverflow.com/a/50714300 +# CC BY-SA 4.0 +def find_edges_with(i, edge_set): + i_first = [j for (x, j) in edge_set if x == i] + i_second = [j for (j, x) in edge_set if x == i] + return i_first, i_second + + +# https://stackoverflow.com/a/50714300 +# CC BY-SA 4.0 +def stitch_boundaries(edges): + edge_set = edges.copy() + boundary_lst = [] + while len(edge_set) > 0: + boundary = [] + edge0 = edge_set.pop() + boundary.append(edge0) + last_edge = edge0 + while len(edge_set) > 0: + i, j = last_edge + j_first, j_second = find_edges_with(j, edge_set) + if j_first: + edge_set.remove((j, j_first[0])) + edge_with_j = (j, j_first[0]) + boundary.append(edge_with_j) + last_edge = edge_with_j + elif j_second: + edge_set.remove((j_second[0], j)) + edge_with_j = (j, j_second[0]) # flip edge rep + boundary.append(edge_with_j) + last_edge = edge_with_j + + if edge0[0] == last_edge[1]: + break + + # boundary_lst.append(boundary) + boundary_lst.extend(split_loops(boundary)) + return boundary_lst + + +# Own work +def split_loops(boundary): + pts = [j for i, j in boundary] + u, c = np.unique(pts, return_counts=True) + dups = u[c == 2] # TODO: edge cases + if len(dups) == 0: + return [boundary] + loops = [] + for dup in dups: + locs = np.where(pts == dup)[0] + if len(locs) != 2: + continue # TODO: loops within loops + loop = pts[locs[0]:locs[1]] + loop.append(dup) + # The loop may potentially have loops itself. Add some recursion. + loops.append([(loop[i], loop[i+1]) for i in range(len(loop)-1)]) + new_pts = pts[:locs[0]] + new_pts.extend(pts[locs[1]:]) + pts = new_pts + pts.append(pts[0]) + boundaries = [[(pts[i], pts[i+1]) for i in range(len(pts)-1)]] + boundaries.extend(loops) + return boundaries + + +# Own work +def boundary_to_poly(boundary, points): + xs = [] + ys = [] + for i, j in boundary: + xs.append(points[i, 0]) + ys.append(points[i, 1]) + xs.append(points[j, 0]) + ys.append(points[j, 1]) + return sg.Polygon([[x, y] for x, y in zip(xs, ys)]) + + +# Own work +def generate_poly_from_edges(edges, points): + def get_poly_with_hole(polys): + biggest = np.argmax([p.area for p in polys]) + outer = polys.pop(biggest) + inners = [] + for idx, poly in enumerate(polys): + if outer.contains(poly): + inners.append(idx) + outer = outer - poly + for index in sorted(inners, reverse=True): + del polys[index] + if type(outer) == sg.MultiPolygon: + return outer.geoms + else: + return [outer] + + boundary_lst = stitch_boundaries(edges) + polys = [boundary_to_poly(b, points) for b in boundary_lst] + outers = [] + while len(polys) > 0: + outers.extend(get_poly_with_hole(polys)) + return outers diff --git a/src/upcp/utils/math_utils.py b/src/upcp/utils/math_utils.py index 9d1d8bc..ebcfd58 100644 --- a/src/upcp/utils/math_utils.py +++ b/src/upcp/utils/math_utils.py @@ -59,8 +59,7 @@ def compute_bounding_box(points): def convex_hull_poly(points): """Return convex hull as a shapely Polygon.""" - hull = points[ConvexHull(points).vertices] - return Polygon(np.vstack((hull, hull[0]))) + return Polygon(points[ConvexHull(points, qhull_options='QJ').vertices]) def minimum_bounding_rectangle(points):