# This file is part of https://github.com/KurtBoehm/svg-path-editor.
#
# This Source Code Form is subject to the terms of the Mozilla Public
# License, v. 2.0. If a copy of the MPL was not distributed with this
# file, You can obtain one at https://mozilla.org/MPL/2.0/.
from __future__ import annotations
from collections.abc import Iterator, Sequence
from dataclasses import dataclass
from decimal import Decimal
from typing import TYPE_CHECKING, Any, overload, override
from .math import (
Boolean,
Expr,
Number,
Precision,
Symbol,
are_equal,
as_bool,
dec_to_rat,
eq,
evalf,
ge,
le,
rat_to_dec,
subs,
)
if TYPE_CHECKING:
import sympy as sp
# ------------------------------------------------------------------------------
# SymPy helpers
# ------------------------------------------------------------------------------
[docs]
def rotation_matrix(phi: Expr) -> Mat2:
r"""
Rotation matrix for angle :math:`φ` in degrees.
Uses
.. math::
R(φ) = \begin{pmatrix}
\cos φ & -\sin φ \\
\sin φ & \cos φ
\end{pmatrix},
where :math:`\cos` and :math:`\sin` are evaluated after converting
:math:`φ` from degrees to radians via :func:`sympy.rad`.
"""
import sympy as sp
rad: sp.Expr = sp.rad(phi)
c, s = sp.cos(rad), sp.sin(rad)
return Mat2(c, -s, s, c)
# ------------------------------------------------------------------------------
# Basic geometric primitives
# ------------------------------------------------------------------------------
[docs]
class Point:
"""2D point with :class:`decimal.Decimal` coordinates."""
def __init__(self, x: Number, y: Number) -> None:
self.x: Decimal = Decimal(x)
self.y: Decimal = Decimal(y)
[docs]
def __iter__(self) -> Iterator[Decimal]:
"""Iterate as ``(x, y)``."""
yield self.x
yield self.y
[docs]
@override
def __eq__(self, other: Any) -> bool:
"""Compare coordinates for equality."""
if isinstance(other, Point):
return self.x == other.x and self.y == other.y
return False
[docs]
@override
def __ne__(self, value: object, /) -> bool:
"""Compare coordinates for inequality."""
return not self == value
@property
def vec2(self) -> Vec2:
"""
Exact conversion to :class:`Vec2`.
Coordinates are converted to SymPy rationals via :func:`dec_to_rat`.
"""
return Vec2.from_point(self)
[docs]
@override
def __str__(self) -> str:
"""Human-readable representation ``(x, y)`` with decimal formatting."""
return f"({self.x:f}, {self.y:f})"
[docs]
@override
def __repr__(self) -> str:
"""Debug representation ``Point(x, y)`` with decimal formatting."""
return f"Point({self.x:f}, {self.y:f})"
@property
def length(self) -> Decimal:
"""Euclidean norm :math:`‖v‖_2 = \\sqrt{x^2 + y^2}`."""
return (self.x * self.x + self.y * self.y).sqrt()
@property
def normalized(self) -> Point:
"""
Unit vector :math:`v / ‖v‖_2`.
The zero vector is returned unchanged.
"""
length = self.length
if length == 0:
return Point(0, 0)
return self / length
# ---- vector arithmetic -------------------------------------------------------
[docs]
def __neg__(self) -> Point:
"""Unary minus :math:`-v`."""
return Point(-self.x, -self.y)
[docs]
def __add__(self, other: Point) -> Point:
"""Vector addition :math:`v + w`."""
return Point(self.x + other.x, self.y + other.y)
[docs]
def __sub__(self, other: Point) -> Point:
"""Vector subtraction :math:`v - w`."""
return Point(self.x - other.x, self.y - other.y)
[docs]
def __mul__(self, other: Number) -> Point:
r"""Scalar multiplication :math:`v ⋅ λ`."""
other = Decimal(other)
return Point(self.x * other, self.y * other)
[docs]
def __truediv__(self, other: Number) -> Point:
"""Scalar division :math:`v / λ`."""
other = Decimal(other)
return Point(self.x / other, self.y / other)
[docs]
@dataclass
class Vec2:
"""
2D vector with SymPy coordinates.
Supports exact arithmetic and simple linear operations.
"""
x: Expr
y: Expr
# ---- construction / conversion -----------------------------------------------
[docs]
@staticmethod
def from_point(p: Point) -> Vec2:
"""
Construct a :class:`Vec2` from a :class:`Point`.
Coordinates are converted to SymPy rationals via :func:`dec_to_rat`.
"""
return Vec2(dec_to_rat(p.x), dec_to_rat(p.y))
@property
def point(self) -> Point:
"""
Convert to numeric :class:`Point`.
Uses :func:`rat_to_dec` to convert SymPy expressions to
:class:`~decimal.Decimal`.
"""
return Point(rat_to_dec(self.x), rat_to_dec(self.y))
# ---- basic protocol ----------------------------------------------------------
[docs]
def __iter__(self) -> Iterator[Expr]:
"""Iterate as ``(x, y)``."""
yield self.x
yield self.y
[docs]
def subs(self, sub: dict[Symbol, Expr], *, n: Precision | None = None) -> Vec2:
"""
Substitute symbols in both coordinates.
:param sub: Substitution dictionary mapping symbols to expressions.
:param n: Optional precision passed through to :func:`subs`.
"""
return Vec2(subs(self.x, sub, n=n), subs(self.y, sub, n=n))
@property
def swapped(self) -> Vec2:
"""Swap coordinates: :math:`(x, y) ↦ (y, x)`."""
return Vec2(self.y, self.x)
[docs]
def evalf(self, *, n: Precision | None = None) -> Vec2:
"""
Evaluate coordinates numerically.
:param n: Optional precision passed to :func:`evalf`.
"""
return Vec2(evalf(self.x, n=n), evalf(self.y, n=n))
[docs]
def simplify(self) -> Vec2:
"""
Simplify both components using :func:`sympy.simplify`.
"""
import sympy as sp
return Vec2(x=sp.simplify(self.x), y=sp.simplify(self.y))
[docs]
@override
def __eq__(self, other: object) -> bool:
"""
Equality comparison between :class:`Vec2` instances,
``False`` when comparing to non-:class:`Vec2`.
"""
import sympy as sp
if isinstance(other, Vec2):
return as_bool(sp.And(eq(self.x, other.x), eq(self.y, other.y)))
return False
[docs]
@override
def __ne__(self, value: object, /) -> bool:
"""
Inequality comparison between :class:`Vec2` instances,
``False`` when comparing to non-:class:`Vec2`.
"""
return not self == value
# ---- elementary geometry -----------------------------------------------------
@property
def length(self) -> Expr:
"""Euclidean norm :math:`‖v‖_2 = \\sqrt{x^2 + y^2}`."""
import sympy as sp
return sp.sqrt(self.x * self.x + self.y * self.y)
@property
def normalized(self) -> Vec2:
"""
Unit vector :math:`v / ‖v‖_2`.
The zero vector is returned unchanged.
"""
import sympy as sp
length = self.length
if are_equal(length, 0):
return Vec2(sp.Integer(0), sp.Integer(0))
return self / length
# ---- vector arithmetic -------------------------------------------------------
[docs]
def __neg__(self) -> Vec2:
"""Unary minus :math:`-v`."""
return Vec2(-self.x, -self.y)
[docs]
def __add__(self, other: Vec2) -> Vec2:
"""Vector addition :math:`v + w`."""
return Vec2(self.x + other.x, self.y + other.y)
[docs]
def __sub__(self, other: Vec2) -> Vec2:
"""Vector subtraction :math:`v - w`."""
return Vec2(self.x - other.x, self.y - other.y)
[docs]
def __mul__(self, other: Expr) -> Vec2:
r"""Scalar multiplication :math:`v ⋅ λ`."""
return Vec2(self.x * other, self.y * other)
[docs]
def __truediv__(self, other: Expr) -> Vec2:
"""Scalar division :math:`v / λ`."""
return Vec2(self.x / other, self.y / other)
@overload
def dot(v1: Point, v2: Point) -> Decimal: ...
@overload
def dot(v1: Vec2, v2: Vec2) -> Expr: ...
[docs]
def dot[V: Point | Vec2](v1: V, v2: V) -> Decimal | Expr:
"""Compute the dot product of two 2D vectors."""
return v1.x * v2.x + v1.y * v2.y
[docs]
@dataclass
class Mat2:
r"""
:math:`2×2` matrix
.. math::
M =
\begin{pmatrix}
a & b \\
c & d
\end{pmatrix}
acting on :class:`Vec2` by standard matrix-vector multiplication.
:ivar a: Entry :math:`a_{11}`.
:ivar b: Entry :math:`a_{12}`.
:ivar c: Entry :math:`a_{21}`.
:ivar d: Entry :math:`a_{22}`.
"""
a: Expr
b: Expr
c: Expr
d: Expr
[docs]
def __matmul__(self, v: Vec2) -> Vec2:
"""
Matrix-vector product ``M @ v``.
.. math::
(x', y') = (a x + b y, \\; c x + d y).
"""
return Vec2(self.a * v.x + self.b * v.y, self.c * v.x + self.d * v.y)
# ------------------------------------------------------------------------------
# Polygon utility
# ------------------------------------------------------------------------------
[docs]
def polygon_signed_area(poly: Sequence[Vec2]) -> Expr:
r"""
Signed area of a simple polygon.
Uses the shoelace formula
.. math::
A = \frac12 \sum_i (x_i y_{i+1} - x_{i+1} y_i),
with positive area for counter-clockwise vertex order.
:param poly: Vertex sequence, implicitly closed.
"""
import sympy as sp
area: sp.Expr = sp.Integer(0)
n = len(poly)
for i in range(n):
x1, y1 = poly[i]
x2, y2 = poly[(i + 1) % n]
area += x1 * y2 - x2 * y1
return area / 2
# ------------------------------------------------------------------------------
# Line segment
# ------------------------------------------------------------------------------
[docs]
class Line:
r"""
Line segment from :attr:`p` to :attr:`q`.
Parametric form
.. math::
L(t) = p + (q - p)\,t, \quad t \in \mathbb{R}.
"""
def __init__(self, p: Vec2, q: Vec2) -> None:
self.p: Vec2 = p
self.q: Vec2 = q
@property
def delta(self) -> Vec2:
"""Direction vector :math:`q - p` of the segment."""
return Vec2(x=self.q.x - self.p.x, y=self.q.y - self.p.y)
@property
def length(self) -> Expr:
"""Euclidean segment length :math:`‖q - p‖_2`."""
return self.delta.length
[docs]
def inward_normal(self, is_ccw: bool) -> Vec2:
r"""
Unit inward normal.
For a CCW-oriented boundary, the inward normal is obtained by rotating
the edge direction :math:`(Δx, Δy)` 90° clockwise; for a CW boundary,
by rotating 90° counter-clockwise:
.. math::
n = \begin{cases}
(Δy, -Δx) & \text{if CCW} \\
(-Δy, Δx) & \text{if CW}
\end{cases}
The resulting vector is then normalized.
:param is_ccw: ``True`` if the enclosing polygon is CCW oriented.
"""
dx, dy = self.delta
n = Vec2(dy, -dx) if is_ccw else Vec2(-dy, dx)
return n.normalized
[docs]
def offset(self, *, d: Expr, is_ccw: bool, n: Precision | None = None) -> Line:
r"""
Offset the line along its inward normal by distance ``d``.
Both endpoints are translated by :math:`d ⋅ n_{\mathit{in}}`.
:param d: Signed offset distance.
:param is_ccw: Orientation of the surrounding boundary.
:param n: Optional precision used in :meth:`Vec2.evalf`.
"""
normal = self.inward_normal(is_ccw=is_ccw) * d
normal = normal.evalf(n=n)
return Line(self.p + normal, self.q + normal)
[docs]
def __call__(self, t: Expr) -> Vec2:
r"""
Evaluate the parametric line :math:`L(t) = p + (q - p)\,t`.
No restriction is imposed on :math:`t`; points with :math:`t \notin [0, 1]`
lie on the infinite supporting line but outside the segment.
"""
return self.p + (self.q - self.p) * t
[docs]
@override
def __str__(self) -> str:
"""Human-readable representation ``(p, q)``."""
return f"({self.p}, {self.q})"
[docs]
@override
def __repr__(self) -> str:
"""Debug representation ``Line(p, q)``."""
return f"Line({self.p!r}, {self.q!r})"
# ------------------------------------------------------------------------------
# Elliptical arc
# ------------------------------------------------------------------------------
[docs]
@dataclass
class ParametricEllipticalArc:
r"""
Elliptical arc in parametric form.
The underlying full ellipse is
.. math::
E(θ) &= R(φ) ⋅
\begin{pmatrix}
r_x \cos θ \\
r_y \sin θ
\end{pmatrix}
+
\begin{pmatrix}
c_x \\
c_y
\end{pmatrix}
where :math:`θ` and :math:`φ` are in degrees and :math:`(c_x, c_y)` is the center.
This arc covers the interval :math:`[θ_0, θ_0 + Δθ]` modulo :math:`360°`.
:ivar c: Center :math:`(c_x, c_y)`.
:ivar r: Radii :math:`(r_x, r_y)`.
:ivar theta0: Start angle :math:`θ_0` in degrees.
:ivar dtheta: Sweep :math:`Δθ` in degrees (signed).
:ivar phi: Rotation angle :math:`φ` in degrees.
"""
c: Vec2
r: Vec2
theta0: Expr
dtheta: Expr
phi: Expr
# ---- conversion --------------------------------------------------------------
@property
def theta1(self) -> Expr:
"""
End angle of the arc.
:return: :math:`θ_1 = θ_0 + Δθ` in degrees.
"""
return self.theta0 + self.dtheta
[docs]
def locally_convex(self, *, is_ccw: bool) -> bool:
"""
Test if the arc is locally convex with respect to the boundary orientation.
:return: ``True`` iff the interior lies on the convex side of the arc
for a boundary with orientation ``is_ccw``.
"""
return as_bool(self.dtheta < 0) == is_ccw
[docs]
def offset(
self,
*,
d: Expr,
is_ccw: bool,
n: Precision | None = None,
) -> ParametricEllipticalArc:
"""
Offset the arc by changing its radii.
* :math:`d > 0`: move inward
* :math:`d < 0`: move outward
“Inward” is defined with respect to the polygon orientation:
for a CCW boundary, the interior is to the *left* of the path,
for a CW boundary to the *right*.
:param d: Signed offset distance applied to both radii.
:param is_ccw: Orientation of the surrounding boundary.
:param n: Unused; kept for API symmetry with :meth:`Line.offset`.
"""
radial_delta = -d if self.locally_convex(is_ccw=is_ccw) else d
return ParametricEllipticalArc(
c=self.c,
r=Vec2(self.r.x + radial_delta, self.r.y + radial_delta),
theta0=self.theta0,
dtheta=self.dtheta,
phi=self.phi,
)
# ---- angle range condition ---------------------------------------------------
@overload
def angle_condition(
self, theta: Symbol, *, n: Precision | None = None
) -> Boolean: ...
@overload
def angle_condition(
self, theta: Expr, *, n: Precision | None = None
) -> Boolean: ...
[docs]
def angle_condition(self, theta: Expr, *, n: Precision | None = None) -> Boolean:
"""
Test whether ``theta`` (in degrees) lies on this arc, modulo 360°.
Works for positive and negative :math:`Δθ` and wrap-around intervals.
:param theta: Angle to test, interpreted in degrees.
:param n: Optional precision passed to :func:`evalf`, :func:`ge`,
and :func:`le`.
"""
import sympy as sp
t0, t1 = evalf(self.theta0, n=n) % 360, evalf(self.theta1, n=n) % 360
dtheta = evalf(self.dtheta, n=n)
theta = evalf(theta, n=n) % 360
lo, hi = (t0, t1) if as_bool(ge(dtheta, sp.S.Zero, n=n)) else (t1, t0)
if as_bool(le(lo, hi, n=n)):
return sp.And(le(lo, theta, n=n), le(theta, hi, n=n))
return sp.Or(le(lo, theta, n=n), le(theta, hi, n=n))
# ---- evaluation and differential geometry -----------------------------------
[docs]
def point_tangent(
self,
theta: Expr,
n: Precision | None = None,
) -> tuple[Vec2, Vec2]:
r"""
Point and tangent at parameter ``theta`` (in degrees).
Returns :math:`(p(θ), ±p'(θ))`, where the derivative is w.r.t. :math:`θ`
in degrees and not normalized. The sign of the derivative is chosen so that
the tangent at :math:`θ_0` points along the arc and that at :math:`θ_1`
points away from the arc.
:param n: Optional precision used in :func:`evalf` and sign checks.
"""
import sympy as sp
rphi = sp.rad(self.phi)
cos_phi, sin_phi = sp.cos(rphi), sp.sin(rphi)
rtheta = sp.rad(theta)
cos_theta, sin_theta = sp.cos(rtheta), sp.sin(rtheta)
c, r = self.c, self.r
# Position
x = c.x + r.x * cos_theta * cos_phi - r.y * sin_theta * sin_phi
y = c.y + r.x * cos_theta * sin_phi + r.y * sin_theta * cos_phi
# Derivative w.r.t. θ (chain rule via d/dθ cos(rad(θ)), sin(rad(θ)))
dxdt = -r.x * sin_theta * cos_phi - r.y * cos_theta * sin_phi
dydt = -r.x * sin_theta * sin_phi + r.y * cos_theta * cos_phi
if as_bool(evalf(self.dtheta, n=n) < 0):
dxdt, dydt = -dxdt, -dydt
return Vec2(x, y).evalf(n=n), Vec2(dxdt, dydt).evalf(n=n)
# ---- transform / implicit form -----------------------------------------------
[docs]
def implicit(self, x: Expr, y: Expr) -> Expr:
r"""
Implicit ellipse equation at ``(x, y)``.
Returns
.. math::
F(x, y) = u^2 + v^2 - 1,
where :math:`(u, v)` is the image of :math:`(x, y)` under the inverse transform
to the unit circle. Points on the ellipse satisfy :math:`F(x, y) = 0`.
"""
uv = self.transform(Vec2(x, y), inverse=True)
return uv.x**2 + uv.y**2 - 1