path: root/lib/stitches/StitchPattern.py

from shapely.geometry.polygon import LinearRing, LineString
from shapely.geometry import Polygon, MultiLineString
from shapely.ops import polygonize
from shapely.geometry import MultiPolygon
from anytree import AnyNode, PreOrderIter
from shapely.geometry.polygon import orient
from depq import DEPQ
from enum import IntEnum
from ..stitches import ConnectAndSamplePattern
from ..stitches import constants


def offset_linear_ring(ring, offset, side, resolution, join_style, mitre_limit):
    """
    Solves following problem: When shapely offsets a LinearRing the
    start/end point might be handled wrongly since they
    are only treated as LineString.
    (See e.g. https://i.stack.imgur.com/vVh56.png as a problematic example)
    This method checks first whether the start/end point form a problematic
    edge with respect to the offset side. If it is not a problematic
    edge we can use the normal offset_routine. Otherwise we need to
    perform two offsets:
    -offset the ring
    -offset the start/end point + its two neighbors left and right
    Finally both offsets are merged together to get the correct
    offset of a LinearRing
    """

    coords = ring.coords[:]
    # check whether edge at index 0 is concave or convex. Only for
    # concave edges we need to spend additional effort
    dx_seg1 = dy_seg1 = 0
    if coords[0] != coords[-1]:
        dx_seg1 = coords[0][0] - coords[-1][0]
        dy_seg1 = coords[0][1] - coords[-1][1]
    else:
        dx_seg1 = coords[0][0] - coords[-2][0]
        dy_seg1 = coords[0][1] - coords[-2][1]
    dx_seg2 = coords[1][0] - coords[0][0]
    dy_seg2 = coords[1][1] - coords[0][1]
    # use cross product:
    crossvalue = dx_seg1 * dy_seg2 - dy_seg1 * dx_seg2
    sidesign = 1
    if side == "left":
        sidesign = -1

    # We do not need to take care of the joint n-0 since we
    # offset along a concave edge:
    if sidesign * offset * crossvalue <= 0:
        return ring.parallel_offset(offset, side, resolution, join_style, mitre_limit)

    # We offset along a convex edge so we offset the joint n-0 separately:
    if coords[0] != coords[-1]:
        coords.append(coords[0])
    offset_ring1 = ring.parallel_offset(
        offset, side, resolution, join_style, mitre_limit
    )
    offset_ring2 = LineString((coords[-2], coords[0], coords[1])).parallel_offset(
        offset, side, resolution, join_style, mitre_limit
    )

    # Next we need to merge the results:
    if offset_ring1.geom_type == "LineString":
        return LinearRing(offset_ring2.coords[:] + offset_ring1.coords[1:-1])
    else:
        # We have more than one resulting LineString for offset of
        # the geometry (ring) = offset_ring1.
        # Hence we need to find the LineString which belongs to the
        # offset of element 0 in coords =offset_ring2
        # in order to add offset_ring2 geometry to it:
        result_list = []
        thresh = constants.offset_factor_for_adjacent_geometry * abs(offset)
        for offsets in offset_ring1:
            if (
                abs(offsets.coords[0][0] - coords[0][0]) < thresh
                and abs(offsets.coords[0][1] - coords[0][1]) < thresh
            ):
                result_list.append(
                    LinearRing(offset_ring2.coords[:] + offsets.coords[1:-1])
                )
            else:
                result_list.append(LinearRing(offsets))
        return MultiLineString(result_list)


def take_only_valid_linear_rings(rings):
    """
    Removes all geometries which do not form a "valid" LinearRing
    (meaning a ring which does not form a straight line)
    """
    if rings.geom_type == "MultiLineString":
        new_list = []
        for ring in rings:
            if len(ring.coords) > 3 or (
                len(ring.coords) == 3 and ring.coords[0] != ring.coords[-1]
            ):
                new_list.append(ring)
        if len(new_list) == 1:
            return LinearRing(new_list[0])
        else:
            return MultiLineString(new_list)
    else:
        if len(rings.coords) <= 2:
            return LinearRing()
        elif len(rings.coords) == 3 and rings.coords[0] == rings.coords[-1]:
            return LinearRing()
        else:
            return rings


def make_tree_uniform_ccw(root):
    """
    Since naturally holes have the opposite point ordering than non-holes we
    make all lines within the tree "root" uniform (having all the same
    ordering direction)
    """
    for node in PreOrderIter(root):
        if node.id == "hole":
            node.val.coords = list(node.val.coords)[::-1]


# Used to define which stitching strategy shall be used
class StitchingStrategy(IntEnum):
    CLOSEST_POINT = 0
    INNER_TO_OUTER = 1
    SPIRAL = 2


def offset_poly(
        poly, offset, join_style, stitch_distance, offset_by_half, strategy, starting_point):
    """
    Takes a polygon (which can have holes) as input and creates offsetted
    versions until the polygon is filled with these smaller offsets.
    These created geometries are afterwards connected to each other and
    resampled with a maximum stitch_distance.
    The return value is a LineString which should cover the full polygon.
    Input:
    -poly: The shapely polygon which can have holes
    -offset: The used offset for the curves
    -join_style: Join style for the offset - can be round, mitered or bevel
     (https://shapely.readthedocs.io/en/stable/manual.html#shapely.geometry.JOIN_STYLE)
     For examples look at
     https://shapely.readthedocs.io/en/stable/_images/parallel_offset.png
    -stitch_distance maximum allowed stitch distance between two points
    -offset_by_half: True if the points shall be interlaced
    -strategy: According to StitchingStrategy enum class you can select between
     different strategies for the connection between parent and childs
    -starting_point: Defines the starting point for the stitching
    Output:
    -List of point coordinate tuples
    -Tag (origin) of each point to analyze why a point was placed
     at this position
    """
    ordered_poly = orient(poly, -1)
    ordered_poly = ordered_poly.simplify(
        constants.simplification_threshold, False)
    root = AnyNode(
        id="node",
        val=ordered_poly.exterior,
        already_rastered=False,
        transferred_point_priority_deque=DEPQ(iterable=None, maxlen=None),
    )
    active_polys = [root]
    active_holes = [[]]

    for holes in ordered_poly.interiors:
        active_holes[0].append(
            AnyNode(
                id="hole",
                val=holes,
                already_rastered=False,
                transferred_point_priority_deque=DEPQ(
                    iterable=None, maxlen=None),
            )
        )

    while len(active_polys) > 0:
        current_poly = active_polys.pop()
        current_holes = active_holes.pop()
        poly_inners = []

        outer = offset_linear_ring(
            current_poly.val,
            offset,
            "left",
            resolution=5,
            join_style=join_style,
            mitre_limit=10,
        )
        outer = outer.simplify(constants.simplification_threshold, False)
        outer = take_only_valid_linear_rings(outer)

        for j in range(len(current_holes)):
            inner = offset_linear_ring(
                current_holes[j].val,
                offset,
                "left",
                resolution=5,
                join_style=join_style,
                mitre_limit=10,
            )
            inner = inner.simplify(constants.simplification_threshold, False)
            inner = take_only_valid_linear_rings(inner)
            if not inner.is_empty:
                poly_inners.append(Polygon(inner))
        if not outer.is_empty:
            if len(poly_inners) == 0:
                if outer.geom_type == "LineString":
                    result = Polygon(outer)
                else:
                    result = MultiPolygon(polygonize(outer))
            else:
                if outer.geom_type == "LineString":
                    result = Polygon(outer).difference(
                        MultiPolygon(poly_inners))
                else:
                    result = MultiPolygon(outer).difference(
                        MultiPolygon(poly_inners))

            if not result.is_empty and result.area > offset * offset / 10:
                result_list = []
                if result.geom_type == "Polygon":
                    result_list = [result]
                else:
                    result_list = list(result)

                for polygon in result_list:
                    polygon = orient(polygon, -1)

                    if polygon.area < offset * offset / 10:
                        continue

                    polygon = polygon.simplify(
                        constants.simplification_threshold, False
                    )
                    poly_coords = polygon.exterior
                    poly_coords = take_only_valid_linear_rings(poly_coords)
                    if poly_coords.is_empty:
                        continue

                    node = AnyNode(
                        id="node",
                        parent=current_poly,
                        val=poly_coords,
                        already_rastered=False,
                        transferred_point_priority_deque=DEPQ(
                            iterable=None, maxlen=None
                        ),
                    )
                    active_polys.append(node)
                    hole_node_list = []
                    for hole in polygon.interiors:
                        hole_node = AnyNode(
                            id="hole",
                            val=hole,
                            already_rastered=False,
                            transferred_point_priority_deque=DEPQ(
                                iterable=None, maxlen=None
                            ),
                        )
                        for previous_hole in current_holes:
                            if Polygon(hole).contains(Polygon(previous_hole.val)):
                                previous_hole.parent = hole_node
                        hole_node_list.append(hole_node)
                    active_holes.append(hole_node_list)
        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 previous_hole.parent is None:
                previous_hole.parent = current_poly

    # DebuggingMethods.drawPoly(root, 'r-')

    make_tree_uniform_ccw(root)
    # print(RenderTree(root))
    if strategy == StitchingStrategy.CLOSEST_POINT:
        (
            connected_line,
            connected_line_origin,
        ) = ConnectAndSamplePattern.connect_raster_tree_nearest_neighbor(
            root, offset, stitch_distance, starting_point, offset_by_half
        )
    elif strategy == StitchingStrategy.INNER_TO_OUTER:
        (
            connected_line,
            connected_line_origin,
        ) = ConnectAndSamplePattern.connect_raster_tree_from_inner_to_outer(
            root, offset, stitch_distance, starting_point, offset_by_half
        )
    else:
        raise ValueError("Invalid stitching stratety!")

    return connected_line, connected_line_origin