orisys-sensor-design/blender/scripts/seam_gap.py

"""Seam margin as a post-process on generated UV polygons (stdlib only).

Why not remove fragments
------------------------
The first attempt simply dropped every fragment that came within the margin of the seam.
That produces a sawtooth edge, because a fragment can only be kept or dropped whole, so the
band boundary snaps to fragment boundaries and its width varies from zero to a full fragment.
The band is also the one place where the width actually matters: it exists to absorb two-pass
registration error, so an uncontrolled width means an uncontrolled overlap risk.

What this does instead
----------------------
The margin is a set of half-spaces cut into the *surface* a fixed distance back from each seam
edge. Within one triangle the surface-to-UV map is affine, so a world-space half-space maps to
a UV-space half-space exactly:

    (p - A) . d >= inset,   p = P0 + M (u - U0)
    =>  c . u >= k,   c = M^T d,   k = inset - (P0 - A) . d + c . U0

where ``M`` is the triangle's 3x2 Jacobian built from its world and UV edges. Clipping the
already-generated UV polygon against those half-planes therefore gives a straight edge at the
requested physical distance, whatever the fragment size and however distorted the layout is in
UV. Nothing upstream changes, which is what makes it a post-process.
"""

from __future__ import annotations

from typing import Iterable, Sequence

# Below this the UV triangle or the projected constraint is degenerate and is skipped rather
# than producing a garbage half-plane.
_EPS = 1e-18


def _sub(a, b):
    return (a[0] - b[0], a[1] - b[1], a[2] - b[2])


def _dot(a, b):
    return a[0] * b[0] + a[1] * b[1] + a[2] * b[2]


def _cross(a, b):
    return (
        a[1] * b[2] - a[2] * b[1],
        a[2] * b[0] - a[0] * b[2],
        a[0] * b[1] - a[1] * b[0],
    )


def _length(a):
    return (_dot(a, a)) ** 0.5


def _normalize(a):
    length = _length(a)
    if length <= _EPS:
        return None
    return (a[0] / length, a[1] / length, a[2] / length)


def triangle_gap_planes(
    world_triangle: Sequence[Sequence[float]],
    uv_triangle: Sequence[Sequence[float]],
    seam_edges: Iterable[Sequence[Sequence[float]]],
    inset_mm: float,
) -> list[tuple[tuple[float, float], float]]:
    """Half-planes ``(c, k)`` in UV, each imposing ``c . u >= k`` for this triangle.

    ``seam_edges`` holds the world-space ``(a, b)`` pairs that are seam edges *of this
    triangle*. Returns an empty list when the triangle is not on the seam, when its UV
    triangle is degenerate, or when a seam edge is parallel to the triangle normal (which
    cannot happen for an edge of the triangle itself, but is guarded anyway).
    """
    if len(world_triangle) < 3 or len(uv_triangle) < 3:
        return []
    p0, p1, p2 = (tuple(map(float, p)) for p in world_triangle[:3])
    u0, u1, u2 = (tuple(map(float, u)) for u in uv_triangle[:3])

    du1 = (u1[0] - u0[0], u1[1] - u0[1])
    du2 = (u2[0] - u0[0], u2[1] - u0[1])
    det = du1[0] * du2[1] - du1[1] * du2[0]
    if abs(det) <= _EPS:
        return []

    dp1 = _sub(p1, p0)
    dp2 = _sub(p2, p0)
    normal = _normalize(_cross(dp1, dp2))
    if normal is None:
        return []
    centroid = tuple((p0[k] + p1[k] + p2[k]) / 3.0 for k in range(3))

    planes: list[tuple[tuple[float, float], float]] = []
    for edge in seam_edges:
        a, b = (tuple(map(float, point)) for point in edge[:2])
        tangent = _normalize(_sub(b, a))
        if tangent is None:
            continue
        direction = _normalize(_cross(tangent, normal))
        if direction is None:
            continue
        # Orient the cut so the triangle's own interior is the kept side.
        if _dot(_sub(centroid, a), direction) < 0.0:
            direction = (-direction[0], -direction[1], -direction[2])
        w0 = _dot(dp1, direction)
        w1 = _dot(dp2, direction)
        # c = (J^-1)^T w with J = [[du1.x, du2.x], [du1.y, du2.y]].
        cx = (du2[1] * w0 - du1[1] * w1) / det
        cy = (-du2[0] * w0 + du1[0] * w1) / det
        c = (cx, cy)
        if abs(cx) <= _EPS and abs(cy) <= _EPS:
            continue
        k = inset_mm - _dot(_sub(p0, a), direction) + (cx * u0[0] + cy * u0[1])
        planes.append((c, k))
    return planes


def _lerp(a, b, t):
    return (a[0] + (b[0] - a[0]) * t, a[1] + (b[1] - a[1]) * t)


def clip_polygon(
    points: Sequence[Sequence[float]],
    planes: Iterable[tuple[tuple[float, float], float]],
) -> list[tuple[float, float]]:
    """Sutherland-Hodgman clip of a UV polygon against ``c . u >= k`` half-planes.

    The polygon stays convex if it starts convex, which every emitted fragment is (fragments
    are triangles), so a single pass per plane is exact.
    """
    out = [(float(p[0]), float(p[1])) for p in points]
    for c, k in planes:
        if len(out) < 3:
            return []
        cx, cy = c
        result: list[tuple[float, float]] = []
        previous = out[-1]
        d_previous = cx * previous[0] + cy * previous[1] - k
        for current in out:
            d_current = cx * current[0] + cy * current[1] - k
            if d_current >= 0.0:
                if d_previous < 0.0:
                    result.append(_lerp(previous, current, d_previous / (d_previous - d_current)))
                result.append(current)
            elif d_previous >= 0.0:
                result.append(_lerp(previous, current, d_previous / (d_previous - d_current)))
            previous, d_previous = current, d_current
        out = result
    return out if len(out) >= 3 else []


def polygon_area(points: Sequence[Sequence[float]]) -> float:
    """Absolute shoelace area, for discarding slivers."""

    def _area(seq):
        total = 0.0
        count = len(seq)
        for index in range(count):
            x0, y0 = seq[index]
            x1, y1 = seq[(index + 1) % count]
            total += x0 * y1 - x1 * y0
        return abs(total) * 0.5

    return _area(list(points)) if len(points) >= 3 else 0.0


__all__ = [
    "clip_polygon",
    "polygon_area",
    "triangle_gap_planes",
]