diff options
Diffstat (limited to 'lib/stitches/contour_fill.py')
| -rw-r--r-- | lib/stitches/contour_fill.py | 551 |
1 files changed, 551 insertions, 0 deletions
diff --git a/lib/stitches/contour_fill.py b/lib/stitches/contour_fill.py new file mode 100644 index 00000000..c42cc6f2 --- /dev/null +++ b/lib/stitches/contour_fill.py @@ -0,0 +1,551 @@ +from collections import namedtuple +from itertools import chain + +import networkx as nx +import numpy as np +import trimesh +from shapely.geometry import GeometryCollection, MultiPolygon, Polygon, LineString, Point +from shapely.geometry.polygon import orient +from shapely.ops import nearest_points +from shapely.ops import polygonize + +from .running_stitch import running_stitch +from ..stitch_plan import Stitch +from ..utils import DotDict +from ..utils.geometry import cut, reverse_line_string, roll_linear_ring +from ..utils.geometry import ensure_geometry_collection, ensure_multi_polygon + + +class Tree(nx.DiGraph): + # This lets us do tree.nodes['somenode'].parent instead of the default + # tree.nodes['somenode']['parent']. + node_attr_dict_factory = DotDict + + def __init__(self, *args, **kwargs): + self.__node_num = 0 + super().__init__(**kwargs) + + def generate_node_name(self): + node = self.__node_num + self.__node_num += 1 + + return node + + +nearest_neighbor_tuple = namedtuple( + "nearest_neighbor_tuple", + [ + "nearest_point_parent", + "nearest_point_child", + "proj_distance_parent", + "child_node", + ], +) + + +def _offset_linear_ring(ring, offset, resolution, join_style, mitre_limit): + result = Polygon(ring).buffer(-offset, resolution, cap_style=2, join_style=join_style, mitre_limit=mitre_limit, single_sided=True) + result = ensure_multi_polygon(result) + rings = GeometryCollection([poly.exterior for poly in result.geoms]) + rings = rings.simplify(0.01, False) + + return _take_only_valid_linear_rings(rings) + + +def _take_only_valid_linear_rings(rings): + """ + Removes all geometries which do not form a "valid" LinearRing. + + A "valid" ring is one that does not form a straight line. + """ + + valid_rings = [] + + for ring in ensure_geometry_collection(rings).geoms: + if len(ring.coords) > 3 or (len(ring.coords) == 3 and ring.coords[0] != ring.coords[-1]): + valid_rings.append(ring) + + return GeometryCollection(valid_rings) + + +def _orient_linear_ring(ring, clockwise=True): + # Unfortunately for us, Inkscape SVGs have an inverted Y coordinate. + # Normally we don't have to care about that, but in this very specific + # case, the meaning of is_ccw is flipped. It actually tests whether + # a ring is clockwise. That makes this logic super-confusing. + if ring.is_ccw != clockwise: + return reverse_line_string(ring) + else: + return ring + + +def _orient_tree(tree, clockwise=True): + """ + Orient all linear rings in the tree. + + Since naturally holes have the opposite point ordering than non-holes we + make all lines within the tree uniform (having all the same ordering + direction) + """ + + for node in tree.nodes.values(): + node.val = _orient_linear_ring(node.val, clockwise) + + +def offset_polygon(polygon, offset, join_style, clockwise): + """ + Convert a polygon to a tree of isocontours. + + An isocontour is an offset version of the polygon's boundary. For example, + the isocontours of a circle are a set of concentric circles inside the + circle. + + This function takes a polygon (which may have holes) as input and creates + isocontours until the polygon is filled completely. The isocontours are + returned as a Tree, with a parent-child relationship indicating that the + parent isocontour contains the child isocontour. + + Arguments: + polygon - The shapely Polygon which may have holes + offset - The spacing between isocontours + join_style - Join style used when offsetting the Polygon border to create + isocontours. Can be round, mitered or bevel, as defined by + shapely: + https://shapely.readthedocs.io/en/stable/manual.html#shapely.geometry.JOIN_STYLE + clockwise - If True, isocontour points are in clockwise order; if False, counter-clockwise. + + Return Value: + Tree - see above + """ + + ordered_polygon = orient(polygon, -1) + tree = Tree() + tree.add_node('root', type='node', parent=None, val=ordered_polygon.exterior) + active_polygons = ['root'] + active_holes = [[]] + + for hole in ordered_polygon.interiors: + hole_node = tree.generate_node_name() + tree.add_node(hole_node, type="hole", val=hole) + active_holes[0].append(hole_node) + + while len(active_polygons) > 0: + current_poly = active_polygons.pop() + current_holes = active_holes.pop() + + outer, inners = _offset_polygon_and_holes(tree, current_poly, current_holes, offset, join_style) + polygons = _match_polygons_and_holes(outer, inners) + + for polygon in polygons.geoms: + new_polygon, new_holes = _convert_polygon_to_nodes(tree, polygon, parent_polygon=current_poly, child_holes=current_holes) + + if new_polygon is not None: + active_polygons.append(new_polygon) + active_holes.append(new_holes) + + for previous_hole in current_holes: + # If the previous holes are not + # contained in the new holes they + # have been merged with the + # outer polygon + if not tree.nodes[previous_hole].parent: + tree.nodes[previous_hole].parent = current_poly + tree.add_edge(current_poly, previous_hole) + + _orient_tree(tree, clockwise) + return tree + + +def _offset_polygon_and_holes(tree, poly, holes, offset, join_style): + outer = _offset_linear_ring( + tree.nodes[poly].val, + offset, + resolution=5, + join_style=join_style, + mitre_limit=10, + ) + + inners = [] + for hole in holes: + inner = _offset_linear_ring( + tree.nodes[hole].val, + -offset, # take negative offset for holes + resolution=5, + join_style=join_style, + mitre_limit=10, + ) + if not inner.is_empty: + inners.append(Polygon(inner.geoms[0])) + + return outer, inners + + +def _match_polygons_and_holes(outer, inners): + result = MultiPolygon(polygonize(outer.geoms)) + if len(inners) > 0: + result = ensure_geometry_collection(result.difference(MultiPolygon(inners))) + + return result + + +def _convert_polygon_to_nodes(tree, polygon, parent_polygon, child_holes): + polygon = orient(polygon, -1) + + if polygon.area < 0.1: + return None, None + + valid_rings = _take_only_valid_linear_rings(polygon.exterior) + + try: + exterior = valid_rings.geoms[0] + except IndexError: + return None, None + + node = tree.generate_node_name() + tree.add_node(node, type='node', parent=parent_polygon, val=exterior) + tree.add_edge(parent_polygon, node) + + hole_nodes = [] + for hole in polygon.interiors: + hole_node = tree.generate_node_name() + tree.add_node(hole_node, type="hole", val=hole) + for previous_hole in child_holes: + if Polygon(hole).contains(Polygon(tree.nodes[previous_hole].val)): + tree.nodes[previous_hole].parent = hole_node + tree.add_edge(hole_node, previous_hole) + hole_nodes.append(hole_node) + + return node, hole_nodes + + +def _get_nearest_points_closer_than_thresh(travel_line, next_line, threshold): + """ + Find the first point along travel_line that is within threshold of next_line. + + Input: + travel_line - The "parent" line for which the distance should + be minimized to enter next_line + next_line - contains the next_line which need to be entered + threshold - The distance between travel_line and next_line needs + to below threshold to be a valid point for entering + + Return value: + tuple or None + - the tuple structure is: + (point in travel_line, point in next_line) + - None is returned if there is no point that satisfies the threshold. + """ + + # We'll buffer next_line and find the intersection with travel_line. + # Then we'll return the very first point in the intersection, + # matched with a corresponding point on next_line. Fortunately for + # us, intersection of a Polygon with a LineString yields pieces of + # the LineString in the same order as the input LineString. + threshold_area = next_line.buffer(threshold) + portion_within_threshold = travel_line.intersection(threshold_area) + + if portion_within_threshold.is_empty: + return None + else: + # Projecting with 0 lets us avoid distinguishing between LineString and + # MultiLineString. + parent_point = Point(portion_within_threshold.interpolate(0)) + return nearest_points(parent_point, next_line) + + +def _create_nearest_points_list(travel_line, tree, children, threshold, threshold_hard): + """Determine the best place to enter each of parent's children + + Arguments: + travel_line - The "parent" line for which the distance should + be minimized to enter each child + children - children of travel_line that need to be entered + threshold - The distance between travel_line and a child should + to be below threshold to be a valid point for entering + threshold_hard - As a last resort, we can accept an entry point + that is this far way + + Return value: + list of nearest_neighbor_tuple - indicating where to enter each + respective child + """ + + children_nearest_points = [] + + for child in children: + result = _get_nearest_points_closer_than_thresh(travel_line, tree.nodes[child].val, threshold) + if result is None: + # where holes meet outer borders a distance + # up to 2 * used offset can arise + result = _get_nearest_points_closer_than_thresh(travel_line, tree.nodes[child].val, threshold_hard) + + proj = travel_line.project(result[0]) + children_nearest_points.append( + nearest_neighbor_tuple( + nearest_point_parent=result[0], + nearest_point_child=result[1], + proj_distance_parent=proj, + child_node=child, + ) + ) + + return children_nearest_points + + +def _find_path_inner_to_outer(tree, node, offset, starting_point, avoid_self_crossing, forward=True): + """Find a stitch path for this ring and its children. + + Strategy: A connection from parent to child is made as fast as possible to + reach the innermost child as fast as possible in order to stitch afterwards + from inner to outer. + + This function calls itself recursively to find a stitch path for each child + (and its children). + + Arguments: + tree - a Tree of isocontours (as returned by offset_polygon) + offset - offset that was passed to offset_polygon + starting_point - starting point for stitching + avoid_self_crossing - if True, tries to generate a path that does not + cross itself. + forward - if True, this ring will be stitched in its natural direction + (used internally by avoid_self_crossing) + + Return value: + LineString -- the stitching path + """ + + current_node = tree.nodes[node] + current_ring = current_node.val + + if not forward and avoid_self_crossing: + current_ring = reverse_line_string(current_ring) + + # reorder the coordinates of this ring so that it starts with + # a point nearest the starting_point + start_distance = current_ring.project(starting_point) + current_ring = roll_linear_ring(current_ring, start_distance) + current_node.val = current_ring + + # Find where along this ring to connect to each child. + nearest_points_list = _create_nearest_points_list( + current_ring, + tree, + tree[node], + threshold=1.5 * offset, + threshold_hard=2.05 * offset + ) + nearest_points_list.sort(key=lambda tup: tup.proj_distance_parent) + + result_coords = [] + if not nearest_points_list: + # We have no children, so we're at the center of a spiral. Reversing + # the innermost ring gives a nicer visual appearance. + if not avoid_self_crossing: + current_ring = reverse_line_string(current_ring) + else: + # This is a recursive algorithm. We'll stitch along this ring, pausing + # to jump to each child ring in turn and sew it before continuing on + # this ring. We'll end back where we started. + + result_coords.append(current_ring.coords[0]) + distance_so_far = 0 + for child_connection in nearest_points_list: + # Cut this ring into pieces before and after where this child will connect. + before, after = cut(current_ring, child_connection.proj_distance_parent - distance_so_far) + distance_so_far = child_connection.proj_distance_parent + + # Stitch the part leading up to this child. + if before is not None: + result_coords.extend(before.coords) + + # Stitch this child. The child will start and end in the same + # place, which should be close to our current location. + child_path = _find_path_inner_to_outer( + tree, + child_connection.child_node, + offset, + child_connection.nearest_point_child, + avoid_self_crossing, + not forward + ) + result_coords.extend(child_path.coords) + + # Skip ahead a little bit on this ring before resuming. This + # gives a nice spiral pattern, where we spiral out from the + # innermost child. + if after is not None: + skip, after = cut(after, offset) + distance_so_far += offset + + current_ring = after + + if current_ring is not None: + # skip a little at the end so we don't end exactly where we started. + remaining_length = current_ring.length + if remaining_length > offset: + current_ring, skip = cut(current_ring, current_ring.length - offset) + + result_coords.extend(current_ring.coords) + + return LineString(result_coords) + + +def inner_to_outer(tree, offset, stitch_length, starting_point, avoid_self_crossing): + """Fill a shape with spirals, from innermost to outermost.""" + + stitch_path = _find_path_inner_to_outer(tree, 'root', offset, starting_point, avoid_self_crossing) + points = [Stitch(*point) for point in stitch_path.coords] + stitches = running_stitch(points, stitch_length) + + return stitches + + +def _reorder_linear_ring(ring, start): + distances = ring - start + start_index = np.argmin(np.linalg.norm(distances, axis=1)) + return np.roll(ring, -start_index, axis=0) + + +def _interpolate_linear_rings(ring1, ring2, max_stitch_length, start=None): + """ + Interpolate between two LinearRings + + Creates a path from start_point on ring1 and around the rings, ending at a + nearby point on ring2. The path will smoothly transition from ring1 to + ring2 as it travels around the rings. + + Inspired by interpolate() from https://github.com/mikedh/pocketing/blob/master/pocketing/polygons.py + + Arguments: + ring1 -- LinearRing start point will lie on + ring2 -- LinearRing end point will lie on + max_stitch_length -- maximum stitch length (used to calculate resampling accuracy) + start -- Point on ring1 to start at, as a tuple + + Return value: Path interpolated between two LinearRings, as a LineString. + """ + + # Resample the two LinearRings so that they are the same number of points + # long. Then take the corresponding points in each ring and interpolate + # between them, gradually going more toward ring2. + # + # This is a little less accurate than the method in interpolate(), but several + # orders of magnitude faster because we're not building and querying a KDTree. + + num_points = int(20 * ring1.length / max_stitch_length) + ring1_resampled = trimesh.path.traversal.resample_path(np.array(ring1.coords), count=num_points) + ring2_resampled = trimesh.path.traversal.resample_path(np.array(ring2.coords), count=num_points) + + if start is not None: + ring1_resampled = _reorder_linear_ring(ring1_resampled, start) + ring2_resampled = _reorder_linear_ring(ring2_resampled, start) + + weights = np.linspace(0.0, 1.0, num_points).reshape((-1, 1)) + points = (ring1_resampled * (1.0 - weights)) + (ring2_resampled * weights) + result = LineString(points) + + return result.simplify(0.1, False) + + +def _check_and_prepare_tree_for_valid_spiral(tree): + """Check whether spiral fill is possible, and tweak if necessary. + + Takes a tree consisting of isocontours. If a parent has more than one child + we cannot create a spiral. However, to make the routine more robust, we + allow more than one child if only one of the children has own children. The + other children are removed in this routine then. If the routine returns true, + the tree will have been cleaned up from unwanted children. + + If even with these weaker constraints, a spiral is not possible, False is + returned. + """ + + def process_node(node): + children = set(tree[node]) + + if len(children) == 0: + return True + elif len(children) == 1: + child = children.pop() + return process_node(child) + else: + children_with_children = {child for child in children if tree[child]} + if len(children_with_children) > 1: + # Node has multiple children with children, so a perfect spiral is not possible. + # This False value will be returned all the way up the stack. + return False + elif len(children_with_children) == 1: + children_without_children = children - children_with_children + child = children_with_children.pop() + tree.remove_nodes_from(children_without_children) + return process_node(child) + else: + # None of the children has its own children, so we'll just take the longest. + longest = max(children, key=lambda child: tree[child]['val'].length) + shorter_children = children - {longest} + tree.remove_nodes_from(shorter_children) + return process_node(longest) + + return process_node('root') + + +def single_spiral(tree, stitch_length, starting_point): + """Fill a shape with a single spiral going from outside to center.""" + return _spiral_fill(tree, stitch_length, starting_point, _make_spiral) + + +def double_spiral(tree, stitch_length, starting_point): + """Fill a shape with a double spiral going from outside to center and back to outside. """ + return _spiral_fill(tree, stitch_length, starting_point, _make_fermat_spiral) + + +def _spiral_fill(tree, stitch_length, close_point, spiral_maker): + starting_point = close_point.coords[0] + + rings = _get_spiral_rings(tree) + path = spiral_maker(rings, stitch_length, starting_point) + path = [Stitch(*stitch) for stitch in path] + + return running_stitch(path, stitch_length) + + +def _get_spiral_rings(tree): + rings = [] + + node = 'root' + while True: + rings.append(tree.nodes[node].val) + + children = tree[node] + if len(children) == 0: + break + elif len(children) == 1: + node = list(children)[0] + else: + # We can only really fill a shape with a single spiral if each + # parent has only one child. We'll do our best though, because + # that is probably more helpful to the user than just refusing + # entirely. We'll pick the child that's closest to the center. + parent_center = rings[-1].centroid + node = min(children, key=lambda child: parent_center.distance(tree.nodes[child].val.centroid)) + + return rings + + +def _make_fermat_spiral(rings, stitch_length, starting_point): + forward = _make_spiral(rings[::2], stitch_length, starting_point) + back = _make_spiral(rings[1::2], stitch_length, starting_point) + back.reverse() + + return chain(forward, back) + + +def _make_spiral(rings, stitch_length, starting_point): + path = [] + + for ring1, ring2 in zip(rings[:-1], rings[1:]): + spiral_part = _interpolate_linear_rings(ring1, ring2, stitch_length, starting_point) + path.extend(spiral_part.coords) + + return path |