repair/project: join polygons directly instead of alphashape

Author: Robert Sachunsky

CHANGELOG.md

@@ -4,6 +4,10 @@ Versioned according to [Semantic Versioning](http://semver.org/).
 
 ## Unreleased
 
+Fixed:
+
+ * repair/project: join polygons directly instead of alphashape
+
 ## [0.1.20] - 2022-04-27
 
 Fixed:

ocrd_segment/project.py

@@ -1,11 +1,10 @@
 from __future__ import absolute_import
 
 import os.path
-from itertools import chain
+import itertools
 import numpy as np
-from shapely.geometry import Polygon
-from shapely.ops import unary_union
-import alphashape
+from shapely.geometry import Polygon, LineString
+from shapely.ops import unary_union, nearest_points
 
 from ocrd import Processor
 from ocrd_utils import (
@@ -144,39 +143,38 @@ def _process_segment(self, segment, constituents, page_id):
             else:
                 segment.set_Coords(coords)
 
-def join_polygons(polygons, loc='', scale=20):
-    """construct concave hull (alpha shape) from input polygons"""
-    LOG = getLogger('processor.ProjectHull')
+def pairwise(iterable):
+    a, b = itertools.tee(iterable)
+    next(b, None)
+    return zip(a, b)
+
+def join_polygons(polygons, scale=20):
+    """construct concave hull (alpha shape) from input polygons by connecting their pairwise nearest points"""
     # ensure input polygons are simply typed
-    polygons = list(chain.from_iterable([
+    polygons = list(itertools.chain.from_iterable([
         poly.geoms if poly.type in ['MultiPolygon', 'GeometryCollection']
         else [poly]
         for poly in polygons]))
-    if len(polygons) == 1:
+    npoly = len(polygons)
+    if npoly == 1:
         return polygons[0]
-    # get equidistant list of points along hull
-    # (otherwise alphashape will jump across the interior)
-    points = [poly.exterior.interpolate(dist).coords[0] # .xy
-              for poly in polygons
-              for dist in np.arange(0, poly.length, min(scale / 2, poly.length / 4))]
-    #alpha = alphashape.optimizealpha(points) # too slow
-    alpha = 0.01
-    jointp = alphashape.alphashape(points, alpha)
-    tries = 0
-    # from descartes import PolygonPatch
-    # import matplotlib.pyplot as plt
-    while jointp.is_empty or jointp.area == 0.0 or jointp.type in ['MultiPolygon', 'GeometryCollection'] or len(jointp.interiors):
-        # plt.figure()
-        # plt.gca().scatter(*zip(*points))
-        # for geom in jointp.geoms:
-        #     plt.gca().add_patch(PolygonPatch(geom, alpha=0.2))
-        # plt.show()
-        alpha *= 0.7
-        tries += 1
-        if tries > 10:
-            LOG.warning("cannot find alpha for concave hull on '%s'", loc)
-            alpha = 0
-        jointp = alphashape.alphashape(points, alpha)
+    # find min-dist path through all polygons (travelling salesman)
+    pairs = itertools.combinations(range(npoly), 2)
+    paths = list(itertools.permutations(range(npoly)))
+    dists = np.eye(npoly, dtype=float)
+    for i, j in pairs:
+        dists[i, j] = polygons[i].distance(polygons[j])
+        dists[j, i] = dists[i, j]
+    dists = [sum(dists[i, j] for i, j in pairwise(path))
+             for path in paths]
+    path = paths[min(enumerate(dists), key=lambda x: x[1])[0]]
+    polygons = [polygons[i] for i in path]
+    # iteratively join to next nearest neighbour
+    jointp = polygons[0]
+    for thisp, nextp in pairwise(polygons):
+        nearest = nearest_points(jointp, nextp)
+        bridgep = LineString(nearest).buffer(max(1, scale/5), resolution=1)
+        jointp = unary_union([jointp, bridgep, nextp])
     if jointp.minimum_clearance < 1.0:
         # follow-up calculations will necessarily be integer;
         # so anticipate rounding here and then ensure validity
@@ -226,8 +224,7 @@ def make_intersection(poly1, poly2):
         interp = unary_union([geom for geom in interp.geoms if geom.area > 0])
     if interp.type == 'MultiPolygon':
         # homogeneous result: construct convex hull to connect
-        # FIXME: construct concave hull / alpha shape
-        interp = interp.convex_hull
+        interp = join_polygons(interp.geoms)
     if interp.minimum_clearance < 1.0:
         # follow-up calculations will necessarily be integer;
         # so anticipate rounding here and then ensure validity
@@ -236,15 +233,16 @@ def make_intersection(poly1, poly2):
     return interp
 
 def make_valid(polygon):
-    for split in range(1, len(polygon.exterior.coords)-1):
+    points = list(polygon.exterior.coords)
+    for split in range(1, len(points)):
         if polygon.is_valid or polygon.simplify(polygon.area).is_valid:
             break
         # simplification may not be possible (at all) due to ordering
         # in that case, try another starting point
-        polygon = Polygon(polygon.exterior.coords[-split:]+polygon.exterior.coords[:-split])
-    for tolerance in range(1, int(polygon.area)):
+        polygon = Polygon(points[-split:]+points[:-split])
+    for tolerance in range(int(polygon.area)):
         if polygon.is_valid:
             break
         # simplification may require a larger tolerance
-        polygon = polygon.simplify(tolerance)
+        polygon = polygon.simplify(tolerance + 1)
     return polygon

ocrd_segment/repair.py

@@ -162,7 +162,7 @@ def process(self):
                             try:
                                 ensure_consistent(element)
                             except Exception as e:
-                                LOG.error(str(e))
+                                LOG.error(str(e)) # exc_info=e
                                 report.add_error(error)
                                 continue
                         else:
@@ -556,7 +556,7 @@ def shrink_regions(page_image, page_coords, page, page_id, padding=0):
         # pick contour and convert to absolute:
         region_polygon = join_polygons([make_valid(Polygon(contour[:, 0, ::]))
                                         for contour in contours
-                                        if len(contour) >= 3], loc=region.id, scale=scale)
+                                        if len(contour) >= 3], scale=scale)
         if padding:
             region_polygon = region_polygon.buffer(padding)
         region_polygon = coordinates_for_segment(region_polygon.exterior.coords[:-1], page_image, page_coords)

requirements.txt

@@ -13,4 +13,3 @@ opencv-python-headless
 # Is not allowed on PyPI:
 #pycocotools@git+https://github.com/bertsky/pycocotools#subdirectory=PythonAPI
 ocrd-fork-pycocotools
-alphashape