mechanical-compiler/src/mechcomp/geom/rounding.py

"""
The BOSL2 operations that ``sb-geom`` and ``sb-join`` call but Shapely has no
equivalent for: corner rounding, and the path cleanup that follows it.

Ported from BOSL2 at commit ``92d697c2856de2fed93a33e858068589cefc2898``, which
is the commit the frozen oracle records. Read from source rather than from
recollection; the call chain is
``round_corners -> _circlecorner -> arc -> segs``.

**Why this is transcribed rather than reimplemented.** A circular arc becomes a
finite number of straight segments, and the segment count determines the
enclosed area. ``SECTION_AREA_MM2`` is compared against the oracle at 1e-3 mm2,
and the extruded solid is these segments -- so the discretisation is not a
quality setting that a smoother modern approach could improve on. It is the
definition of the surface. Any adaptive subdivision or tolerance-based
flattening would produce a better curve and a failed build.

The generators set ``$fn = facets`` with ``facets = 48``, and all 123 oracle
cases run at that default. With ``$fn`` positive, ``segs()`` ignores the radius
entirely, so segment counts depend only on swept angle. That is reproduced here
via the module-level default; it is a parameter rather than a constant because
the generators expose it as one.

``round_corners`` **raises** when the requested roundovers do not fit the path,
matching BOSL2, which asserts rather than clamping. Callers are expected to
derive safe radii up front -- that is exactly what ``sb_corner_radii`` is for.
Silently clamping here would let a build succeed where the reference aborted.
"""

from __future__ import annotations

import math
from typing import List, Optional, Sequence, Tuple

Point = Tuple[float, float]
Path = Sequence[Point]

# BOSL2's private library epsilon (math.scad).
EPSILON = 1e-9

# OpenSCAD facet count set by the generators. Overridable, but every frozen
# oracle case was produced at this value.
DEFAULT_FN = 48


class RoundoverTooLarge(ValueError):
    """
    Raised where BOSL2 asserts "Roundovers are too big for the path."

    A distinct type because this is a caller error -- radii that were never
    derived against the path -- and not a rejected profile.
    """


# ----------------------------------------------------------------------------
# Comparison helpers (comparisons.scad)
# ----------------------------------------------------------------------------

def approx(a: float, b: float, eps: float = EPSILON) -> bool:
    return abs(a - b) <= eps


def approx_pt(a: Point, b: Point, eps: float = EPSILON) -> bool:
    """Componentwise, as BOSL2's list branch of ``approx``."""
    return abs(a[0] - b[0]) <= eps and abs(a[1] - b[1]) <= eps


def deduplicate(path: Path, closed: bool = False,
                eps: float = EPSILON) -> List[Point]:
    """
    Drop each point that equals the next one.

    Closed paths compare the last point against the first, which is how a
    roundover that fully consumes a segment gets collapsed.
    """
    pts = list(path)
    n = len(pts)
    if n == 0:
        return []
    end = n if closed else n - 1
    return [pts[i] for i in range(n)
            if i == end or not approx_pt(pts[i], pts[(i + 1) % n], eps)]


def _dist2line(d: Point, n: Point) -> float:
    dot = d[0] * n[0] + d[1] * n[1]
    return math.hypot(d[0] - dot * n[0], d[1] - dot * n[1])


def is_collinear(pts: Sequence[Point], eps: float = EPSILON) -> bool:
    """
    BOSL2's ``is_collinear`` via ``_noncollinear_triple``.

    The tolerance is *relative*: the furthest point from the first defines the
    chord, and collinearity holds when every point lies within ``eps`` times
    that chord length of the line. An absolute tolerance would behave
    differently at the scales this library works at.
    """
    if len(pts) < 3:
        return True
    pa = pts[0]
    b = max(range(len(pts)), key=lambda i: math.hypot(pts[i][0] - pa[0],
                                                      pts[i][1] - pa[1]))
    pb = pts[b]
    nrm = math.hypot(pb[0] - pa[0], pb[1] - pa[1])
    if nrm <= eps:
        return True
    n = ((pb[0] - pa[0]) / nrm, (pb[1] - pa[1]) / nrm)
    distlist = [_dist2line((p[0] - pa[0], p[1] - pa[1]), n) for p in pts]
    return max(distlist) < eps * nrm


def path_merge_collinear(path: Path, closed: bool = True,
                         eps: float = EPSILON) -> List[Point]:
    """
    Remove vertices that lie on the line between their neighbours.

    Exact butt joints and zero-radius fillets leave collinear vertices. They are
    harmless in 2D but leave zero-area triangles the tessellator cannot resolve,
    so a section that measures perfectly can still fail to extrude.
    """
    pts = deduplicate(path, closed=closed, eps=eps)
    n = len(pts)
    if n <= 2:
        return pts
    out: List[Point] = []
    if not closed:
        out.append(pts[0])
        rng = range(1, n - 1)
    else:
        rng = range(n)
    for i in rng:
        triple = (pts[(i - 1) % n], pts[i], pts[(i + 1) % n])
        if not is_collinear(triple, eps=eps):
            out.append(triple[1])
    if not closed:
        out.append(pts[-1])
    return out


# ----------------------------------------------------------------------------
# Segmentation (utility.scad: segs)
# ----------------------------------------------------------------------------

def segs(r: float, angle: Optional[float] = None,
         fn: int = DEFAULT_FN, fa: float = 12.0, fs: float = 2.0) -> int:
    """
    Number of sides OpenSCAD gives a circle, or an arc of ``angle`` degrees.

    The ``2e-15`` subtraction is BOSL2's, guarding an angle that is fractionally
    over its true value through rounding. It is reproduced because dropping it
    can add a segment at exactly 360-divisible angles.
    """
    if angle is not None:
        return math.ceil(segs(r, None, fn, fa, fs) * abs(angle) / 360.0 - 2e-15)
    if fn > 0:
        return fn if fn > 3 else 3
    rr = r if math.isfinite(r) else 0.0
    return math.ceil(max(5.0, min(360.0 / fa, abs(rr) * 2.0 * math.pi / fs)))


# ----------------------------------------------------------------------------
# Arc through two points about a centre (drawing.scad: arc)
# ----------------------------------------------------------------------------

def _vector_angle3(a: Point, b: Point, c: Point) -> float:
    """Angle at ``b``, degrees, in [0, 180]."""
    ux, uy = a[0] - b[0], a[1] - b[1]
    vx, vy = c[0] - b[0], c[1] - b[1]
    nu, nv = math.hypot(ux, uy), math.hypot(vx, vy)
    if nu == 0.0 or nv == 0.0:
        return 0.0
    cosv = (ux * vx + uy * vy) / (nu * nv)
    return math.degrees(math.acos(max(-1.0, min(1.0, cosv))))


def arc(n: int, cp: Point, points: Tuple[Point, Point]) -> List[Point]:
    """
    ``n`` points along the arc from ``points[0]`` to ``points[1]`` about ``cp``.

    Sweep direction follows the sign of the 2D cross product, taking the short
    way round -- BOSL2's ``long``/``cw``/``ccw`` flags are never passed by the
    call sites this port needs, so the minor arc is always the one drawn.
    """
    start, end = points
    angle = _vector_angle3(start, cp, end)
    v1 = (start[0] - cp[0], start[1] - cp[1])
    v2 = (end[0] - cp[0], end[1] - cp[1])
    prelim = v1[0] * v2[1] - v1[1] * v2[0]
    direction = 1.0 if prelim > 0 else (-1.0 if prelim < 0 else 1.0)
    r = math.hypot(v1[0], v1[1])
    final_angle = direction * angle
    sa = math.degrees(math.atan2(v1[1], v1[0]))
    out: List[Point] = []
    for i in range(n):
        theta = sa + i * final_angle / (n - 1)
        out.append((r * math.cos(math.radians(theta)) + cp[0],
                    r * math.sin(math.radians(theta)) + cp[1]))
    return out


def _circlecorner(points: Tuple[Point, Point, Point],
                  parm: Tuple[float, float], fn: int = DEFAULT_FN) -> List[Point]:
    """
    One rounded corner: ``parm`` is ``(d, r)``, the tangent setback and radius.

    A straight vertex -- half-angle 90 degrees -- degenerates to the two tangent
    points with no arc between them.
    """
    prev_p, here, next_p = points
    angle = _vector_angle3(prev_p, here, next_p) / 2.0
    d, r = parm

    pux, puy = prev_p[0] - here[0], prev_p[1] - here[1]
    nux, nuy = next_p[0] - here[0], next_p[1] - here[1]
    pn = math.hypot(pux, puy)
    nn = math.hypot(nux, nuy)
    prev_u = (pux / pn, puy / pn)
    next_u = (nux / nn, nuy / nn)

    start = (here[0] + prev_u[0] * d, here[1] + prev_u[1] * d)
    end = (here[0] + next_u[0] * d, here[1] + next_u[1] * d)

    if approx(angle, 90.0):
        return [start, end]

    bx, by = prev_u[0] + next_u[0], prev_u[1] + next_u[1]
    bn = math.hypot(bx, by)
    scale = r / math.sin(math.radians(angle))
    center = (scale * bx / bn + here[0], scale * by / bn + here[1])

    # The segment count is a ceiling on a value that is frequently an exact
    # integer -- a 60 degree half-angle at $fn=48 gives exactly 8. Floating
    # point delivers that as 8.000000000000004 or 7.999999999999998 depending
    # on how the corner was reached, and the ceiling then differs by a whole
    # segment, changing the enclosed area by about a thousandth of a square
    # millimetre.
    #
    # This is reproduced unguarded because the reference is unguarded. Rounding
    # first was tried and is wrong: it fixes the Y profile and breaks Three-Fin,
    # because BOSL2's own output is asymmetric in exactly this way and the
    # oracle records that asymmetry. See FAILURES.md F-034.
    raw = (90.0 - angle) / 180.0 * segs(r, None, fn)
    n = max(3, math.ceil(raw))
    return arc(n, center, (start, end))


# ----------------------------------------------------------------------------
# round_corners (rounding.scad), method="circle", measure="radius"
# ----------------------------------------------------------------------------

def round_corners(path: Path, radius, closed: bool = True,
                  fn: int = DEFAULT_FN) -> List[Point]:
    """
    Round each corner of ``path`` to its own radius.

    ``radius`` is a scalar or one value per vertex. Zero leaves a vertex
    untouched, which is how ``sb_fillet_concave`` rounds only reflex corners
    while keeping every convex corner bit-exact.

    Raises ``RoundoverTooLarge`` when the setbacks overrun an edge, matching
    BOSL2's assertion. The message carries the same scale factors BOSL2 reports,
    since those say directly how much too large the request was.
    """
    pts = list(path)
    n = len(pts)
    if n < 3:
        raise ValueError("Path has length %d. Length must be 3 or more." % n)

    parm = [float(radius)] * n if isinstance(radius, (int, float)) \
        else [float(x) for x in radius]
    if len(parm) != n:
        raise ValueError("radius list length %d does not match path length %d"
                         % (len(parm), n))

    dk: List[Tuple[float, ...]] = []
    for i in range(n):
        bit = (pts[(i - 1) % n], pts[i], pts[(i + 1) % n])
        degenerate = approx_pt(bit[0], bit[1]) or approx_pt(bit[1], bit[2])
        angle = None if degenerate else _vector_angle3(*bit) / 2.0

        if not closed and (i == 0 or i == n - 1):
            dk.append((0.0,))
            continue
        if parm[i] == 0:
            dk.append((0.0,))
            continue
        if angle is None:
            raise ValueError("Repeated point in path at index %d with nonzero "
                             "rounding" % i)
        if approx(angle, 0.0):
            raise ValueError("Path turns back on itself at index %d with "
                             "nonzero rounding" % i)
        dk.append((parm[i] / math.tan(math.radians(angle)), parm[i]))

    lengths = [math.hypot(pts[i % n][0] - pts[(i - 1) % n][0],
                          pts[i % n][1] - pts[(i - 1) % n][1])
               for i in range(n + 1)]

    scalefactors: List[float] = []
    for i in range(n):
        if not (closed or (i != 0 and i != n - 1)):
            continue
        back = dk[(i - 1) % n][0] + dk[i][0]
        fwd = dk[i][0] + dk[(i + 1) % n][0]
        scalefactors.append(min(
            math.inf if back == 0 else lengths[i] / back,
            math.inf if fwd == 0 else lengths[i + 1] / fwd,
        ))

    if scalefactors and min(scalefactors) < 1.0:
        raise RoundoverTooLarge(
            "Roundovers are too big for the path. If you multiply them by this "
            "vector they should fit: %r" % (scalefactors,))

    out: List[Point] = []
    for i in range(n):
        if dk[i][0] == 0:
            out.append(pts[i])
        else:
            bit = (pts[(i - 1) % n], pts[i], pts[(i + 1) % n])
            out.extend(_circlecorner(bit, (dk[i][0], dk[i][1]), fn))
    return deduplicate(out, closed=False)