path: root/lib/stitches/tangential_fill_stitch_pattern_creator.py

import math
from collections import namedtuple
import networkx as nx
import numpy as np
import trimesh
from depq import DEPQ
from shapely.geometry import Point, LineString, LinearRing, MultiLineString
from shapely.ops import nearest_points

from .running_stitch import running_stitch

from ..debug import debug
from ..stitches import constants
from ..stitches import point_transfer
from ..stitches import sample_linestring
from ..stitch_plan import Stitch
from ..utils.geometry import roll_linear_ring

nearest_neighbor_tuple = namedtuple(
    "nearest_neighbor_tuple",
    [
        "nearest_point_parent",
        "nearest_point_child",
        "proj_distance_parent",
        "child_node",
    ],
)


@debug.time
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

    Output:
    -tuple or None
      - the tuple structure is:
        (nearest point in travel_line, nearest 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:
        if isinstance(portion_within_threshold, MultiLineString):
            portion_within_threshold = portion_within_threshold.geoms[0]

        parent_point = Point(portion_within_threshold.coords[0])
        return nearest_points(parent_point, next_line)


def create_nearest_points_list(
        travel_line, tree, children_list, threshold, threshold_hard, preferred_direction=0):
    """
    Takes a line and calculates the nearest distance along this line to
    enter the childs in children_list
    The method calculates the distances along the line and along the
    reversed line to find the best direction which minimizes the overall
    distance for all childs.
    Input:
    -travel_line: The "parent" line for which the distance should
     be minimized to enter the childs
    -children_list: contains the childs of travel_line which need to be entered
    -threshold: The distance between travel_line and a child needs to be
     below threshold to be a valid point for entering
    -preferred_direction: Put a bias on the desired travel direction along
     travel_line. If equals zero no bias is applied.
     preferred_direction=1 means we prefer the direction of travel_line;
     preferred_direction=-1 means we prefer the opposite direction.
    Output:
    -stitching direction for travel_line
    -list of tuples (one tuple per child). The tuple structure is:
     ((nearest point in travel_line, nearest point in child),
       distance along travel_line, belonging child)
    """

    result_list_in_order = []
    result_list_reversed_order = []

    travel_line_reversed = LinearRing(travel_line.coords[::-1])

    weight_in_order = 0
    weight_reversed_order = 0
    for child in children_list:
        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
            )
            assert result is not None
        proj = travel_line.project(result[0])
        weight_in_order += proj
        result_list_in_order.append(
            nearest_neighbor_tuple(
                nearest_point_parent=result[0],
                nearest_point_child=result[1],
                proj_distance_parent=proj,
                child_node=child,
            )
        )

        result = get_nearest_points_closer_than_thresh(
            travel_line_reversed, 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_reversed, tree.nodes[child].val, threshold_hard
            )
            assert result is not None
        proj = travel_line_reversed.project(result[0])
        weight_reversed_order += proj
        result_list_reversed_order.append(
            nearest_neighbor_tuple(
                nearest_point_parent=result[0],
                nearest_point_child=result[1],
                proj_distance_parent=proj,
                child_node=child,
            )
        )

    if preferred_direction == 1:
        # Reduce weight_in_order to make in order stitching more preferred
        weight_in_order = min(
            weight_in_order / 2, max(0, weight_in_order - 10 * threshold)
        )
        if weight_in_order == weight_reversed_order:
            return 1, result_list_in_order
    elif preferred_direction == -1:
        # Reduce weight_reversed_order to make reversed
        # stitching more preferred
        weight_reversed_order = min(
            weight_reversed_order /
            2, max(0, weight_reversed_order - 10 * threshold)
        )
        if weight_in_order == weight_reversed_order:
            return (-1, result_list_reversed_order)

    if weight_in_order < weight_reversed_order:
        return (1, result_list_in_order)
    else:
        return (-1, result_list_reversed_order)


def calculate_replacing_middle_point(line_segment, abs_offset, max_stitch_distance):
    """
    Takes a line segment (consisting of 3 points!)
    and calculates a new middle point if the line_segment is
    straight enough to be resampled by points max_stitch_distance apart FROM THE END OF line_segment.
    Returns None if the middle point is not needed.
    """
    angles = sample_linestring.calculate_line_angles(line_segment)
    if angles[1] < abs_offset * constants.limiting_angle_straight:
        if line_segment.length < max_stitch_distance:
            return None
        else:
            return line_segment.interpolate(line_segment.length - max_stitch_distance).coords[0]
    else:
        return line_segment.coords[1]


def connect_raster_tree_from_inner_to_outer(tree, node, used_offset, stitch_distance, min_stitch_distance, close_point,
                                            offset_by_half):  # noqa: C901
    """
    Takes the offsetted curves organized as tree, connects and samples them.
    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.
    Input:
    -tree: contains the offsetted curves in a hierachical organized
     data structure.
    -used_offset: used offset when the offsetted curves were generated
    -stitch_distance: maximum allowed distance between two points
     after sampling
    -min_stitch_distance stitches within a row shall be at least min_stitch_distance apart. Stitches connecting
     offsetted paths might be shorter.
    -close_point: defines the beginning point for stitching
     (stitching starts always from the undisplaced curve)
    -offset_by_half: If true the resulting points are interlaced otherwise not.
    Returnvalues:
    -All offsetted curves connected to one line and sampled with points obeying
     stitch_distance and offset_by_half
    -Tag (origin) of each point to analyze why a point was placed
     at this position
    """

    current_node = tree.nodes[node]
    current_coords = current_node.val
    abs_offset = abs(used_offset)
    result_coords = []
    result_coords_origin = []

    start_distance = current_coords.project(close_point)
    # We cut the current path so that its index 0 is closest to close_point
    if start_distance > 0:
        current_coords = roll_linear_ring(current_coords, start_distance)
        current_node.val = current_coords

        if not current_node.transferred_point_priority_deque.is_empty():
            new_DEPQ = DEPQ(iterable=None, maxlen=None)
            for item, priority in current_node.transferred_point_priority_deque:
                new_DEPQ.insert(
                    item,
                    math.fmod(
                        priority - start_distance + current_coords.length,
                        current_coords.length,
                    ),
                )
            current_node.transferred_point_priority_deque = new_DEPQ

    # We try to use always the opposite stitching direction with respect to the
    # parent to avoid crossings when entering and leaving the child
    # LEX: this seems like a lie ^^
    parent_stitching_direction = -1
    if current_node.parent is not None:
        parent_stitching_direction = tree.nodes[current_node.parent].stitching_direction

    # Find the nearest point in current_coords and its children and
    # sort it along the stitching direction
    stitching_direction, nearest_points_list = create_nearest_points_list(
        current_coords,
        tree,
        tree[node],
        constants.offset_factor_for_adjacent_geometry * abs_offset,
        2.05 * abs_offset,
        parent_stitching_direction,
    )
    nearest_points_list.sort(
        reverse=False, key=lambda tup: tup.proj_distance_parent)

    # Have a small offset for the starting and ending to avoid double points
    # at start and end point (since the paths are closed rings)
    if nearest_points_list:
        start_offset = min(
            abs_offset * constants.factor_offset_starting_points,
            nearest_points_list[0].proj_distance_parent,
        )
        end_offset = max(
            current_coords.length
            - abs_offset * constants.factor_offset_starting_points,
            nearest_points_list[-1].proj_distance_parent,
        )
    else:
        start_offset = abs_offset * constants.factor_offset_starting_points
        end_offset = (current_coords.length - abs_offset * constants.factor_offset_starting_points)

    if stitching_direction == 1:
        (own_coords, own_coords_origin) = sample_linestring.raster_line_string_with_priority_points(
            current_coords,
            start_offset,  # We add start_offset to not sample the initial/end
            # point twice (avoid double points for start
            # and end)
            end_offset,
            stitch_distance,
            min_stitch_distance,
            current_node.transferred_point_priority_deque,
            abs_offset,
            offset_by_half,
            False
        )
    else:
        (own_coords, own_coords_origin) = sample_linestring.raster_line_string_with_priority_points(
            current_coords,
            current_coords.length - start_offset,  # We subtract
            # start_offset to not
            # sample the initial/end point
            # twice (avoid double
            # points for start
            # and end)
            current_coords.length - end_offset,
            stitch_distance,
            min_stitch_distance,
            current_node.transferred_point_priority_deque,
            abs_offset,
            offset_by_half,
            False
        )
        current_coords.coords = current_coords.coords[::-1]

    assert len(own_coords) == len(own_coords_origin)

    current_node.stitching_direction = stitching_direction
    current_node.already_rastered = True

    to_transfer_point_list = []
    to_transfer_point_list_origin = []
    for k in range(0, len(own_coords)):
        # TODO: maybe do not take the first and the last
        # since they are ENTER_LEAVING_POINT points for sure
        if (
                not offset_by_half
                and own_coords_origin[k] == sample_linestring.PointSource.EDGE_NEEDED
                or own_coords_origin[k] == sample_linestring.PointSource.FORBIDDEN_POINT):
            continue
        if own_coords_origin[k] == sample_linestring.PointSource.ENTER_LEAVING_POINT:
            continue
        to_transfer_point_list.append(Point(own_coords[k]))
        to_transfer_point_list_origin.append(own_coords_origin[k])

    assert len(to_transfer_point_list) == len(to_transfer_point_list_origin)

    # Next we need to transfer our rastered points to siblings and childs
    # Since the projection is only in ccw direction towards inner we
    # need to use "-used_offset" for stitching_direction==-1
    point_transfer.transfer_points_to_surrounding(
        tree,
        node,
        stitching_direction * used_offset,
        offset_by_half,
        to_transfer_point_list,
        to_transfer_point_list_origin,
        overnext_neighbor=False,
        transfer_forbidden_points=False,
        transfer_to_parent=False,
        transfer_to_sibling=True,
        transfer_to_child=True,
    )

    # We transfer also to the overnext child to get a more straight
    # arrangement of points perpendicular to the stitching lines
    if offset_by_half:
        point_transfer.transfer_points_to_surrounding(
            tree,
            node,
            stitching_direction * used_offset,
            False,
            to_transfer_point_list,
            to_transfer_point_list_origin,
            overnext_neighbor=True,
            transfer_forbidden_points=False,
            transfer_to_parent=False,
            transfer_to_sibling=True,
            transfer_to_child=True,
        )

    if not nearest_points_list:
        # If there is no child (inner geometry) we can simply
        # take our own rastered coords as result
        result_coords = own_coords
        result_coords_origin = own_coords_origin
    else:
        # There are childs so we need to merge their coordinates
        # with our own rastered coords

        # Create a closed ring for the following code
        own_coords.append(own_coords[0])
        own_coords_origin.append(own_coords_origin[0])

        # own_coords does not start with current_coords but has an offset
        # (see call of raster_line_string_with_priority_points)
        total_distance = start_offset

        cur_item = 0
        result_coords = [own_coords[0]]
        result_coords_origin = [own_coords_origin[0]]

        for i in range(1, len(own_coords)):
            next_distance = math.sqrt(
                (own_coords[i][0] - own_coords[i - 1][0]) ** 2
                + (own_coords[i][1] - own_coords[i - 1][1]) ** 2
            )
            while (
                cur_item < len(nearest_points_list)
                and total_distance + next_distance + constants.eps
                > nearest_points_list[cur_item].proj_distance_parent
            ):
                # The current and the next point in own_coords enclose the
                # nearest point tuple between this geometry and child
                # geometry. Hence we need to insert the child geometry points
                # here before the next point of own_coords.
                item = nearest_points_list[cur_item]
                (child_coords, child_coords_origin) = connect_raster_tree_from_inner_to_outer(
                    tree,
                    item.child_node,
                    used_offset,
                    stitch_distance,
                    min_stitch_distance,
                    item.nearest_point_child,
                    offset_by_half,
                )

                # Imagine the nearest point of the child is within a long
                # segment of the parent. Without additonal points
                # on the parent side this would cause noticeable deviations.
                # Hence we add here points shortly before and after
                # the entering of the child to have only minor deviations to
                # the desired shape.
                # Here is the point for the entering:
                if (Point(result_coords[-1]).distance(item.nearest_point_parent) > constants.factor_offset_starting_points * abs_offset):
                    result_coords.append(item.nearest_point_parent.coords[0])
                    result_coords_origin.append(
                        sample_linestring.PointSource.ENTER_LEAVING_POINT
                    )

                # Check whether the number of points of the connecting lines
                # from child to child can be reduced
                if len(child_coords) > 1:
                    point = calculate_replacing_middle_point(
                        LineString(
                            [result_coords[-1], child_coords[0], child_coords[1]]
                        ),
                        abs_offset,
                        stitch_distance,
                    )

                    if point is not None:
                        result_coords.append(point)
                        result_coords_origin.append(child_coords_origin[0])

                    result_coords.extend(child_coords[1:])
                    result_coords_origin.extend(child_coords_origin[1:])
                else:
                    result_coords.extend(child_coords)
                    result_coords_origin.extend(child_coords_origin)

                # And here is the point for the leaving of the child
                # (distance to the own following point should not be too large)
                d = item.nearest_point_parent.distance(Point(own_coords[i]))
                if cur_item < len(nearest_points_list) - 1:
                    d = min(
                        d,
                        abs(nearest_points_list[cur_item + 1].proj_distance_parent - item.proj_distance_parent),
                    )

                if d > constants.factor_offset_starting_points * abs_offset:
                    result_coords.append(
                        current_coords.interpolate(item.proj_distance_parent + 2 * constants.factor_offset_starting_points * abs_offset).coords[0]
                    )
                    result_coords_origin.append(
                        sample_linestring.PointSource.ENTER_LEAVING_POINT
                    )
                    # Check whether this additional point makes the last point
                    # of the child unnecessary
                    point = calculate_replacing_middle_point(
                        LineString(
                            [result_coords[-3], result_coords[-2], result_coords[-1]]
                        ),
                        abs_offset,
                        stitch_distance,
                    )
                    if point is None:
                        result_coords.pop(-2)
                        result_coords_origin.pop(-2)

                cur_item += 1
            if i < len(own_coords) - 1:
                if (Point(result_coords[-1]).distance(Point(own_coords[i])) > abs_offset * constants.factor_offset_remove_points):
                    result_coords.append(own_coords[i])
                    result_coords_origin.append(own_coords_origin[i])

            # Since current_coords and own_coords are rastered differently
            # there accumulate errors regarding the current distance.
            # Since a projection of each point in own_coords would be very
            # time consuming we project only every n-th point which resets
            # the accumulated error every n-th point.
            if i % 20 == 0:
                total_distance = current_coords.project(Point(own_coords[i]))
            else:
                total_distance += next_distance

    assert len(result_coords) == len(result_coords_origin)
    return result_coords, result_coords_origin


def orient_linear_ring(ring):
    if not ring.is_ccw:
        return LinearRing(reversed(ring.coords))
    else:
        return ring


def reorder_linear_ring(ring, start):
    # TODO: actually use start?
    start_index = np.argmin(np.linalg.norm(ring, 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.
    """

    ring1 = orient_linear_ring(ring1)
    ring2 = orient_linear_ring(ring2)

    # 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(ring1, count=num_points)
    ring2_resampled = trimesh.path.traversal.resample_path(ring2, 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)

    # TODO: remove when rastering is cheaper
    return result.simplify(constants.simplification_threshold, False)


def connect_raster_tree_spiral(tree, used_offset, stitch_distance, min_stitch_distance, close_point, offset_by_half):  # noqa: C901
    """
    Takes the offsetted curves organized as tree, connects and samples them as a spiral.
    It expects that each node in the tree has max. one child
    Input:
    -tree: contains the offsetted curves in a hierarchical organized
     data structure.
    -used_offset: used offset when the offsetted curves were generated
    -stitch_distance: maximum allowed distance between two points
     after sampling
    -min_stitch_distance stitches within a row shall be at least min_stitch_distance apart. Stitches connecting
     offsetted paths might be shorter.
    -close_point: defines the beginning point for stitching
     (stitching starts always from the undisplaced curve)
    -offset_by_half: If true the resulting points are interlaced otherwise not.
    Return values:
    -All offsetted curves connected to one spiral and sampled with
     points obeying stitch_distance and offset_by_half
    -Tag (origin) of each point to analyze why a point was
     placed at this position
    """

    if not tree['root']:  # if node has no children
        stitches = [Stitch(*point) for point in tree.nodes['root'].val.coords]
        return running_stitch(stitches, stitch_distance)

    starting_point = close_point.coords[0]
    path = []
    for node in nx.dfs_preorder_nodes(tree, 'root'):
        if tree[node]:
            ring1 = tree.nodes[node].val
            child = list(tree.successors(node))[0]
            ring2 = tree.nodes[child].val

            spiral_part = interpolate_linear_rings(ring1, ring2, stitch_distance, starting_point)
            path.extend(spiral_part.coords)

    path = [Stitch(*stitch) for stitch in path]

    return running_stitch(path, stitch_distance), None