path: root/lib/stitches/running_stitch.py

# Authors: see git history
#
# Copyright (c) 2010 Authors
# Licensed under the GNU GPL version 3.0 or later.  See the file LICENSE for details.

import math
from math import tau
import typing
from copy import copy

import numpy as np
from shapely import geometry as shgeo
from ..utils import prng
from ..utils.geometry import Point
from ..utils.threading import check_stop_flag

""" Utility functions to produce running stitches. """


def split_segment_even_n(a, b, segments: int, jitter_sigma: float = 0.0, random_seed=None) -> typing.List[shgeo.Point]:
    if segments <= 1:
        return []
    line = shgeo.LineString((a, b))

    splits = np.array(range(1, segments)) / segments
    if random_seed is not None:
        jitters = (prng.n_uniform_floats(len(splits), random_seed) * 2) - 1
        splits = splits + jitters * (jitter_sigma / segments)

    # sort the splits in case a bad roll transposes any of them
    return [line.interpolate(x, normalized=True) for x in sorted(splits)]


def split_segment_even_dist(a, b, max_length: float, jitter_sigma: float = 0.0, random_seed=None) -> typing.List[shgeo.Point]:
    distance = shgeo.Point(a).distance(shgeo.Point(b))
    segments = math.ceil(distance / max_length)
    return split_segment_even_n(a, b, segments, jitter_sigma, random_seed)


def split_segment_random_phase(a, b, length: float, length_sigma: float, random_seed: str) -> typing.List[shgeo.Point]:
    line = shgeo.LineString([a, b])
    progress = length * prng.uniform_floats(random_seed, "phase")[0]
    splits = [progress]
    distance = line.length
    if progress >= distance:
        return []
    for x in prng.iter_uniform_floats(random_seed):
        progress += length * (1 + length_sigma * (x - 0.5) * 2)
        if progress >= distance:
            break
        splits.append(progress)
    return [line.interpolate(x, normalized=False) for x in splits]


class AngleInterval():
    # Modular interval containing either the entire circle or less than half of it
    # partially based on https://fgiesen.wordpress.com/2015/09/24/intervals-in-modular-arithmetic/

    def __init__(self, a: float, b: float, all: bool = False):
        self.all = all
        self.a = a
        self.b = b

    @staticmethod
    def all():
        return AngleInterval(0, math.tau, True)

    @staticmethod
    def fromBall(p: Point, epsilon: float):
        d = p.length()
        if d <= epsilon:
            return AngleInterval.all()
        center = p.angle()
        delta = math.asin(epsilon / d)
        return AngleInterval(center - delta, center + delta)

    @staticmethod
    def fromSegment(a: Point, b: Point):
        angleA = a.angle()
        angleB = b.angle()
        diff = (angleB - angleA) % tau
        if diff == 0 or diff == math.pi:
            return None
        elif diff < math.pi:
            return AngleInterval(angleA - 1e-6, angleB + 1e-6)
            # slightly larger than normal to avoid rounding error when this method is used in cutSegment
        else:
            return AngleInterval(angleB - 1e-6, angleA + 1e-6)

    def containsAngle(self, angle: float):
        if self.all:
            return True
        return (angle - self.a) % tau <= (self.b - self.a) % tau

    def containsPoint(self, p: Point):
        return self.containsAngle(math.atan2(p.y, p.x))

    def intersect(self, other):
        # assume that each interval contains less than half the circle (or all of it)
        if other is None:
            return None
        elif self.all:
            return other
        elif other.all:
            return self
        elif self.containsAngle(other.a):
            if other.containsAngle(self.b):
                return AngleInterval(other.a, self.b)
            else:
                return AngleInterval(other.a, other.b)
        elif other.containsAngle(self.a):
            if self.containsAngle(other.b):
                return AngleInterval(self.a, other.b)
            else:
                return AngleInterval(self.a, self.b)
        else:
            return None

    def cutSegment(self, origin: Point, a: Point, b: Point):
        if self.all:
            return None
        segArc = AngleInterval.fromSegment(a - origin, b - origin)
        if segArc is None:
            return a  # b is exactly behind origin from a
        if segArc.containsAngle(self.a):
            return cut_segment_with_angle(origin, self.a, a, b)
        elif segArc.containsAngle(self.b):
            return cut_segment_with_angle(origin, self.b, a, b)
        else:
            return None


def cut_segment_with_angle(origin: Point, angle: float, a: Point, b: Point) -> Point:
    # Assumes the crossing is inside the segment
    p = a - origin
    d = b - a
    c = Point(math.cos(angle), math.sin(angle))
    t = (p.y*c.x - p.x*c.y) / (d.x*c.y - d.y*c.x)
    if t < -0.000001 or t > 1.000001:
        raise Exception("cut_segment_with_angle returned a parameter of {0} with points {1} {2} and cut line {3} ".format(t, p, b-origin, c))
    return a + d*t


def cut_segment_with_circle(origin: Point, r: float, a: Point, b: Point) -> Point:
    # assumes that a is inside the circle and b is outside
    p = a - origin
    d = b - a
    # inner products
    p2 = p * p
    d2 = d * d
    r2 = r * r
    pd = p * d
    # r2 = p2 + 2*pd*t + d2*t*t, quadratic formula
    t = (math.sqrt(pd*pd + r2*d2 - p2*d2) - pd) / d2
    if t < -0.000001 or t > 1.000001:
        raise Exception("cut_segment_with_circle returned a parameter of {0}".format(t))
    return a + d*t


def take_stitch(start: Point, points: typing.Sequence[Point], idx: int, stitch_length: float, tolerance: float):
    # Based on a single step of the Zhao-Saalfeld curve simplification algorithm.
    # https://cartogis.org/docs/proceedings/archive/auto-carto-13/pdf/linear-time-sleeve-fitting-polyline-simplification-algorithms.pdf
    # Adds early termination condition based on stitch length.
    if idx >= len(points):
        return None, None

    sleeve = AngleInterval.all()
    last = start
    for i in range(idx, len(points)):
        p = points[i]
        if sleeve.containsPoint(p - start):
            if start.distance(p) < stitch_length:
                sleeve = sleeve.intersect(AngleInterval.fromBall(p - start, tolerance))
                last = p
                continue
            else:
                cut = cut_segment_with_circle(start, stitch_length, last, p)
                return cut, i
        else:
            cut = sleeve.cutSegment(start, last, p)
            if start.distance(cut) > stitch_length:
                cut = cut_segment_with_circle(start, stitch_length, last, p)
            return cut, i
    return points[-1], None


def stitch_curve_evenly(points: typing.Sequence[Point], stitch_length: float, tolerance: float):
    # Will split a straight line into even-length stitches while still handling curves correctly.
    # Includes end point but not start point.
    if len(points) < 2:
        return []
    distLeft = [0] * len(points)
    for i in reversed(range(0, len(points) - 1)):
        distLeft[i] = distLeft[i + 1] + points[i].distance(points[i+1])

    i = 1
    last = points[0]
    stitches = []
    while i is not None and i < len(points):
        check_stop_flag()

        d = last.distance(points[i]) + distLeft[i]
        if d == 0:
            return stitches
        stitch_len = d / math.ceil(d / stitch_length) + 0.000001  # correction for rounding error

        stitch, newidx = take_stitch(last, points, i, stitch_len, tolerance)
        i = newidx
        if stitch is not None:
            stitches.append(stitch)
            last = stitch
    return stitches


def path_to_curves(points: typing.List[Point], min_len: float):
    # split a path at obvious corner points so that they get stitched exactly
    # min_len controls the minimum length after splitting for which it won't split again,
    # which is used to avoid creating large numbers of corner points when encouintering micro-messes.
    if len(points) < 3:
        return [points]
    curves = []

    last = 0
    last_seg = points[1] - points[0]
    seg_len = last_seg.length()
    for i in range(1, len(points) - 1):
        # vectors of the last and next segments
        a = last_seg
        b = points[i + 1] - points[i]
        aabb = (a * a) * (b * b)
        abab = (a * b) * abs(a * b)

        # Test if the turn angle from vectors a to b is more than 45 degrees.
        # Optimized version of checking if cos(angle(a,b)) <= sqrt(0.5) and is defined
        if aabb > 0 and abab <= 0.5 * aabb:
            if seg_len >= min_len:
                curves.append(points[last: i + 1])
                last = i
            seg_len = 0

        if b * b > 0:
            last_seg = b
        seg_len += b.length()

    curves.append(points[last:])
    return curves


def running_stitch(points, stitch_length, tolerance):
    # Turn a continuous path into a running stitch.
    stitches = [points[0]]
    for curve in path_to_curves(points, 2 * tolerance):
        # segments longer than twice the tollerance will usually be forced by it, so set that as the minimum for corner detection
        stitches.extend(stitch_curve_evenly(curve, stitch_length, tolerance))
    return stitches


def bean_stitch(stitches, repeats):
    """Generate bean stitch from a set of stitches.

    "Bean" stitch is made by backtracking each stitch to make it heavier.  A
    simple bean stitch would be two stitches forward, one stitch back, two
    stitches forward, etc.  This would result in each stitch being tripled.

    We'll say that the above counts as 1 repeat.  Backtracking each stitch
    repeatedly will result in a heavier bean stitch.  There will always be
    an odd number of threads piled up for each stitch.

    Repeats is a list of a repeated pattern e.g. [0, 1, 3] doesn't repeat the first stitch,
    goes back and forth on the second stitch, goes goes 3 times back and forth on the third stitch,
    and starts the pattern again by not repeating the fourth stitch, etc.
    """

    if len(stitches) < 2:
        return stitches

    repeat_list_length = len(repeats)
    repeat_list_pos = 0

    new_stitches = [stitches[0]]

    for stitch in stitches:
        new_stitches.append(stitch)

        for i in range(repeats[repeat_list_pos]):
            new_stitches.extend(copy(new_stitches[-2:]))

        repeat_list_pos += 1
        if repeat_list_pos == repeat_list_length:
            repeat_list_pos = 0

    return new_stitches