6 files changed, 256 insertions, 14 deletions

diff --git a/lib/utils/clamp_path.py b/lib/utils/clamp_path.py
new file mode 100644
index 00000000..e5ef78d8
--- /dev/null
+++ b/lib/utils/clamp_path.py
@@ -0,0 +1,121 @@
+from shapely.geometry import LineString, Point as ShapelyPoint, MultiPolygon
+from shapely.prepared import prep
+from .geometry import Point, ensure_multi_line_string
+
+
+def path_to_segments(path):
+    """Convert a path of Points into a list of segments as LineStrings"""
+    for start, end in zip(path[:-1], path[1:]):
+        yield LineString((start, end))
+
+
+def segments_to_path(segments):
+    """Convert a list of contiguous LineStrings into a list of Points."""
+    coords = [segments[0].coords[0]]
+
+    for segment in segments:
+        coords.extend(segment.coords[1:])
+
+    return [Point(x, y) for x, y in coords]
+
+
+def fix_starting_point(border_pieces):
+    """Reconnect the starting point of a polygon border's pieces.
+
+    When splitting a polygon border with two lines, we want to get two
+    pieces.  However, that's not quite how Shapely works.  The outline
+    of the polygon is a LinearRing that starts and ends at the same place,
+    but Shapely still knows where that starting point is and splits there
+    too.
+
+    We don't want that third piece, so we'll reconnect the segments that
+    touch the starting point.
+    """
+
+    if len(border_pieces) == 3:
+        # Fortunately, Shapely keeps the starting point of the LinearRing
+        # as the starting point of the first segment.  That means it's also
+        # the ending point of the last segment.  Reconnecting is super simple:
+        return [border_pieces[1],
+                LineString(border_pieces[2].coords[:] + border_pieces[0].coords[1:])]
+    else:
+        # We probably cut exactly at the starting point.
+        return border_pieces
+
+
+def adjust_line_end(line, end):
+    """Reverse line if necessary to ensure that it ends near end."""
+
+    line_start = ShapelyPoint(*line.coords[0])
+    line_end = ShapelyPoint(*line.coords[-1])
+
+    if line_end.distance(end) < line_start.distance(end):
+        return line
+    else:
+        return LineString(line.coords[::-1])
+
+
+def find_border(polygon, point):
+    for border in polygon.interiors:
+        if border.intersects(point):
+            return border
+    else:
+        return polygon.exterior
+
+
+def clamp_path_to_polygon(path, polygon):
+    """Constrain a path to a Polygon.
+
+    Description: https://gis.stackexchange.com/questions/428848/clamp-linestring-to-polygon
+    """
+
+    path = LineString(path)
+
+    # This splits the path at the points where it intersects with the polygon
+    # border and returns the pieces in the same order as the original path.
+    split_path = ensure_multi_line_string(path.difference(polygon.boundary))
+
+    # contains() checks can fail without this.
+    buffered_polygon = prep(polygon.buffer(1e-9))
+
+    last_segment_inside = None
+    was_inside = False
+    result = []
+
+    for segment in split_path.geoms:
+        if buffered_polygon.contains(segment):
+            if not was_inside:
+                if last_segment_inside is not None:
+                    # The path crossed out of the polygon, and now it's crossed
+                    # back in.  We need to add a path along the border between
+                    # the exiting and entering points.
+
+                    # First, find the two points.  Buffer them just a bit to
+                    # ensure intersection with the border.
+                    x, y = last_segment_inside.coords[-1]
+                    exit_point = ShapelyPoint(x, y).buffer(0.01, resolution=1)
+                    x, y = segment.coords[0]
+                    entry_point = ShapelyPoint(x, y).buffer(0.01, resolution=1)
+
+                    if not exit_point.intersects(entry_point):
+                        # Now break the border into pieces using those points.
+                        border = find_border(polygon, exit_point)
+                        border_pieces = border.difference(MultiPolygon((entry_point, exit_point))).geoms
+                        border_pieces = fix_starting_point(border_pieces)
+
+                        # Pick the shortest way to get from the exiting to the
+                        # entering point along the border.
+                        shorter = min(border_pieces, key=lambda piece: piece.length)
+
+                        # We don't know which direction the polygon border
+                        # piece should be.  adjust_line_end() will figure
+                        # that out.
+                        result.append(adjust_line_end(shorter, entry_point))
+
+            result.append(segment)
+            was_inside = True
+            last_segment_inside = segment
+        else:
+            was_inside = False
+
+    return segments_to_path(result)
diff --git a/lib/utils/geometry.py b/lib/utils/geometry.py
index 789f8720..8f34c467 100644
--- a/lib/utils/geometry.py
+++ b/lib/utils/geometry.py
@@ -220,6 +220,9 @@ class Point:
     def rotate(self, angle):
         return self.__class__(self.x * math.cos(angle) - self.y * math.sin(angle), self.y * math.cos(angle) + self.x * math.sin(angle))
 
+    def scale(self, x_scale, y_scale):
+        return self.__class__(self.x * x_scale, self.y * y_scale)
+
     def as_int(self):
         return self.__class__(int(round(self.x)), int(round(self.y)))
 
@@ -238,3 +241,7 @@ class Point:
 
 def line_string_to_point_list(line_string):
     return [Point(*point) for point in line_string.coords]
+
+
+def coordinate_list_to_point_list(coordinate_list):
+    return [Point.from_tuple(coords) for coords in coordinate_list]
diff --git a/lib/utils/list.py b/lib/utils/list.py
new file mode 100644
index 00000000..efa3969e
--- /dev/null
+++ b/lib/utils/list.py
@@ -0,0 +1,23 @@
+import random
+
+
+def _uniform_rng():
+    while True:
+        yield random.uniform(0, 1)
+
+
+_rng = _uniform_rng()
+
+
+def poprandom(sequence, rng=_rng):
+    index = int(round(next(rng) * (len(sequence) - 1)))
+    item = sequence[index]
+
+    # It's O(1) to pop the last item, and O(n) to pop any other item. So we'll
+    # always pop the last item and put it in the slot vacated by the item we're
+    # popping.
+    last_item = sequence.pop()
+    if index < len(sequence):
+        sequence[index] = last_item
+
+    return item
diff --git a/lib/utils/prng.py b/lib/utils/prng.py
index 2ec037c6..33102205 100644
--- a/lib/utils/prng.py
+++ b/lib/utils/prng.py
@@ -3,7 +3,7 @@ from math import ceil
 from itertools import count, chain
 import numpy as np
 
-# Framework for reproducable pseudo-random number generation.
+# Framework for reproducible pseudo-random number generation.
 
 # Unlike python's random module (which uses a stateful generator based on global variables),
 # a counter-mode PRNG like uniformFloats can be used to generate multiple, independent random streams
@@ -13,7 +13,7 @@ import numpy as np
 # Using multiple counters for n-dimentional random streams is also possible and is useful for grid-like structures.
 
 
-def joinArgs(*args):
+def join_args(*args):
     # Stringifies parameters into a slash-separated string for use in hash keys.
     # Idempotent and associative.
     return "/".join([str(x) for x in args])
@@ -22,37 +22,37 @@ def joinArgs(*args):
 MAX_UNIFORM_INT = 2 ** 32 - 1
 
 
-def uniformInts(*args):
+def uniform_ints(*args):
     # Single pseudo-random drawing determined by the joined parameters.
     # To get a longer sequence of random numbers, call this loop with a counter as one of the parameters.
     # Returns 8 uniformly random uint32.
 
-    s = joinArgs(*args)
+    s = join_args(*args)
     # blake2s is python's fastest hash algorithm for small inputs and is designed to be usable as a PRNG.
     h = blake2s(s.encode()).hexdigest()
     nums = []
     for i in range(0, 64, 8):
-        nums.append(int(h[i:i+8], 16))
+        nums.append(int(h[i:i + 8], 16))
     return np.array(nums)
 
 
-def uniformFloats(*args):
+def uniform_floats(*args):
     # Single pseudo-random drawing determined by the joined parameters.
     # To get a longer sequence of random numbers, call this loop with a counter as one of the parameters.
     # Returns an array of 8 floats in the range [0,1]
-    return uniformInts(*args) / MAX_UNIFORM_INT
+    return uniform_ints(*args) / MAX_UNIFORM_INT
 
 
-def nUniformFloats(n: int, *args):
+def n_uniform_floats(n: int, *args):
     # returns a fixed number (which may exceed 8) of floats in the range [0,1]
-    seed = joinArgs(*args)
-    nBlocks = ceil(n/8)
-    blocks = [uniformFloats(seed, x) for x in range(nBlocks)]
+    seed = join_args(*args)
+    nBlocks = ceil(n / 8)
+    blocks = [uniform_floats(seed, x) for x in range(nBlocks)]
     return np.concatenate(blocks)[0:n]
 
 
-def iterUniformFloats(*args):
+def iter_uniform_floats(*args):
     # returns an infinite sequence of floats in the range [0,1]
-    seed = joinArgs(*args)
-    blocks = map(lambda x: list(uniformFloats(seed, x)), count(0))
+    seed = join_args(*args)
+    blocks = map(lambda x: list(uniform_floats(seed, x)), count(0))
     return chain.from_iterable(blocks)
diff --git a/lib/utils/smoothing.py b/lib/utils/smoothing.py
new file mode 100644
index 00000000..9d43a9f1
--- /dev/null
+++ b/lib/utils/smoothing.py
@@ -0,0 +1,84 @@
+import numpy as np
+from scipy.interpolate import splprep, splev
+
+from .geometry import Point, coordinate_list_to_point_list
+from ..stitches.running_stitch import running_stitch
+from ..debug import debug
+
+
+def _remove_duplicate_coordinates(coords_array):
+    """Remove consecutive duplicate points from an array.
+
+    Arguments:
+        coords_array -- numpy.array
+
+    Returns:
+        a numpy.array of coordinates, minus consecutive duplicates
+    """
+
+    differences = np.diff(coords_array, axis=0)
+    zero_differences = np.isclose(differences, 0)
+    keepers = np.r_[True, np.any(zero_differences == False, axis=1)]  # noqa: E712
+
+    return coords_array[keepers]
+
+
+@debug.time
+def smooth_path(path, smoothness=1.0):
+    """Smooth a path of coordinates.
+
+    Arguments:
+        path -- an iterable of coordinate tuples or Points
+        smoothness -- float, how much smoothing to apply.  Bigger numbers
+            smooth more.
+
+    Returns:
+        A list of Points.
+    """
+    from ..debug import debug
+
+    if smoothness == 0:
+        # s of exactly zero seems to indicate a default level of smoothing
+        # in splprep, so we'll just exit instead.
+        return path
+
+    # Smoothing seems to look nicer if the line segments in the path are mostly
+    # similar in length.  If we have some especially long segments, then the
+    # smoothed path sometimes diverges more from the original path as the
+    # spline curve struggles to fit the path.  This can be especially bad at
+    # the start and end.
+    #
+    # Fortunately, we can convert the path to segments that are mostly the same
+    # length by using the running stitch algorithm.
+    path = running_stitch(coordinate_list_to_point_list(path), 5 * smoothness, smoothness / 2)
+
+    # splprep blows up on duplicated consecutive points with "Invalid inputs"
+    coords = _remove_duplicate_coordinates(np.array(path))
+    num_points = len(coords)
+
+    if num_points <= 3:
+        # splprep throws an error unless num_points > k
+        return path
+
+    # s is explained in this issue: https://github.com/scipy/scipy/issues/11916
+    # the smoothness parameter limits how much the smoothed path can deviate
+    # from the original path.  The standard deviation of the distance between
+    # the smoothed path and the original path is equal to the smoothness.
+    # In practical terms, if smoothness is 1mm, then the smoothed path can be
+    # up to 1mm away from the original path.
+    s = num_points * (smoothness ** 2)
+
+    # .T transposes the array (for some reason splprep expects
+    # [[x1, x2, ...], [y1, y2, ...]]
+    tck, fp, ier, msg = splprep(coords.T, s=s, k=3, nest=-1, full_output=1)
+    if ier > 0:
+        debug.log(f"error {ier} smoothing path: {msg}")
+        return path
+
+    # Evaluate the spline curve at many points along its length to produce the
+    # smoothed point list.  2 * num_points seems to be a good number, but it
+    # does produce a lot of points.
+    smoothed_x_values, smoothed_y_values = splev(np.linspace(0, 1, int(num_points * 2)), tck[0])
+    coords = np.array([smoothed_x_values, smoothed_y_values]).T
+
+    return [Point(x, y) for x, y in coords]
diff --git a/lib/utils/string.py b/lib/utils/string.py
index cb852ce3..e9204076 100644
--- a/lib/utils/string.py
+++ b/lib/utils/string.py
@@ -8,3 +8,10 @@ def string_to_floats(string, delimiter=","):
 
     floats = string.split(delimiter)
     return [float(num) for num in floats]
+
+
+def remove_suffix(string, suffix):
+    if string.endswith(suffix):
+        return string[:-len(suffix)]
+    else:
+        return string