"""Classes associated with general curves."""
from __future__ import annotations
from abc import ABC, abstractmethod
import numpy as np
from buffalo_core.numeric import as_float_array, as_float_scalar
from buffalo_core.typing import (
FloatArray,
FloatInput,
FloatScalar,
)
from scipy.integrate import quad
from scipy.optimize import root_scalar
from buffalo_wings.type_aliases import CurveBreakpointSides
from .runtime_common import (
QUAD_ABS_TOLERANCE,
QUAD_LIMIT,
QUAD_REL_TOLERANCE,
ROOT_ABS_TOLERANCE,
ROOT_MAX_ITERATION,
)
_BREAKPOINT_PARAMETER_STEP = 1e-7
[docs]
class Curve(ABC):
"""
Base class for 1-d curves.
Curves can be interrogated based on their specific parameterization and by
their arc-length parameterization. The specific parametrization variable is
``u``, and the arc-length parameterization variable is ``s`` and is a
measure of the distance from the start of the curve (at the minimum ``u``).
Arc-length queries are more expensive because the mapping from surface
distance to the native parameter is not available in closed form for
general curve shapes.
Breakpoints are the native-parameter locations where one-sided derivative
information matters.
Interior breakpoints are reserved for locations where first- or
higher-derivative behavior may differ on the two sides.
Endpoints are always included as one-sided boundary markers.
The ordinary derivative evaluators, such as :meth:`xy_u`,
:meth:`xy_uu`, :meth:`xy_s`, and :meth:`xy_ss`, use the ``minus``-side
value when the query lands exactly on a reported breakpoint.
The paired ``*_breakpoint`` methods expose both one-sided values
explicitly.
Subclasses should override the breakpoint methods whenever exact
one-sided values are available and should rely on the generic sampled
fallback only when no exact representation is available.
"""
#
# Native curve geometry
#
[docs]
@abstractmethod
def xy_from_u(self, u: FloatInput) -> tuple[FloatArray, FloatArray]:
"""
Calculate the coordinates of geometry at parameter location.
Parameters
----------
u : FloatInput
Parameter for desired locations.
Returns
-------
tuple[FloatArray, FloatArray]
Tuple ``(x, y)`` of ``float64`` arrays matching the normalized
shape of ``u``.
"""
[docs]
@abstractmethod
def xy_u(self, u: FloatInput) -> tuple[FloatArray, FloatArray]:
"""
Calculate rates of change of the coordinates at parameter location.
Parameters
----------
u : FloatInput
Parameter for desired locations.
Returns
-------
tuple[FloatArray, FloatArray]
Tuple ``(dx/du, dy/du)`` of ``float64`` arrays matching the
normalized shape of ``u``.
Notes
-----
If ``u`` matches one of :meth:`breakpoints` exactly, this method
returns the ``minus``-side derivative limit.
"""
[docs]
@abstractmethod
def xy_uu(self, u: FloatInput) -> tuple[FloatArray, FloatArray]:
"""
Calculate second derivative of the coordinates at parameter location.
Parameters
----------
u : FloatInput
Parameter for desired locations.
Returns
-------
tuple[FloatArray, FloatArray]
Tuple ``(d^2x/du^2, d^2y/du^2)`` of ``float64`` arrays matching the
normalized shape of ``u``.
Notes
-----
If ``u`` matches one of :meth:`breakpoints` exactly, this method
returns the ``minus``-side derivative limit.
"""
#
# Derived curve geometry
#
[docs]
def tangent(self, u: FloatInput) -> tuple[FloatArray, FloatArray]:
"""
Calculate the unit tangent at parameter location.
Parameters
----------
u : FloatInput
Parameter for desired locations.
Returns
-------
tuple[FloatArray, FloatArray]
Tuple ``(t_x, t_y)`` of ``float64`` arrays matching the normalized
shape of ``u``.
"""
tx, ty = self.xy_u(u)
temp = np.sqrt(tx**2 + ty**2)
tx = as_float_array(np.divide(tx, temp))
ty = as_float_array(np.divide(ty, temp))
return tx, ty
[docs]
def normal(self, u: FloatInput) -> tuple[FloatArray, FloatArray]:
"""
Calculate the unit normal at parameter location.
Parameters
----------
u : FloatInput
Parameter for desired locations.
Returns
-------
tuple[FloatArray, FloatArray]
Tuple ``(n_x, n_y)`` of ``float64`` arrays matching the normalized
shape of ``u``.
"""
sx, sy = self.tangent(u)
nx = -sy
ny = sx
return nx, ny
[docs]
def k(self, u: FloatInput) -> FloatArray:
"""
Calculate the curvature at parameter location.
Parameters
----------
u : FloatInput
Parameter for desired locations.
Returns
-------
FloatArray
Curvature of surface matching the normalized shape of ``u``.
"""
xu, yu = self.xy_u(u)
xuu, yuu = self.xy_uu(u)
return (xu * yuu - yu * xuu) / (xu**2 + yu**2) ** (3 / 2)
@property
def length(self) -> FloatScalar:
"""
Return the total curve length.
Returns
-------
FloatScalar
Arc length measured from the minimum to maximum native
parameter.
"""
breakpoints = self._breakpoints()
return as_float_scalar(self.arc_length(breakpoints[0], breakpoints[-1]))
#
# Breakpoint interface
#
[docs]
def xy_u_breakpoint(
self,
*,
index: int,
) -> CurveBreakpointSides:
"""
Return both sides of first derivatives at one breakpoint.
Parameters
----------
index : int
Index into :meth:`breakpoints`.
Returns
-------
CurveBreakpointSides
``((x_u_minus, y_u_minus), (x_u_plus, y_u_plus))``.
Notes
-----
Endpoint breakpoints return the same boundary value for both entries.
This method is the exact-breakpoint contract that pairs with
:meth:`xy_u`.
Subclasses should override it whenever they can return exact
one-sided derivative values.
The generic implementation evaluates nearby one-sided parameter
samples and therefore serves only as an approximation fallback.
"""
u_minus, u_plus = self.breakpoint_parameter_limits(index=index)
minus_derivative = self.xy_u(u_minus)
plus_derivative = self.xy_u(u_plus)
minus_values = (
as_float_scalar(minus_derivative[0]),
as_float_scalar(minus_derivative[1]),
)
plus_values = (
as_float_scalar(plus_derivative[0]),
as_float_scalar(plus_derivative[1]),
)
return minus_values, plus_values
[docs]
def xy_uu_breakpoint(
self,
*,
index: int,
) -> CurveBreakpointSides:
"""
Return one-sided second derivatives at one breakpoint.
Parameters
----------
index : int
Index into :meth:`breakpoints`.
Returns
-------
CurveBreakpointSides
``((x_uu_minus, y_uu_minus), (x_uu_plus, y_uu_plus))``.
Notes
-----
Endpoint breakpoints return the same boundary value for both entries.
This method is the exact-breakpoint contract that pairs with
:meth:`xy_uu`.
Subclasses should override it whenever they can return exact
one-sided second-derivative values.
The generic implementation evaluates nearby one-sided parameter
samples and therefore serves only as an approximation fallback.
"""
u_minus, u_plus = self.breakpoint_parameter_limits(index=index)
minus_values = self.xy_uu(u_minus)
plus_values = self.xy_uu(u_plus)
return (
as_float_scalar(minus_values[0]),
as_float_scalar(minus_values[1]),
), (
as_float_scalar(plus_values[0]),
as_float_scalar(plus_values[1]),
)
[docs]
@abstractmethod
def breakpoints(self) -> list[FloatScalar]:
"""
Return the sorted locations of any breakpoints in the curve.
The resulting list must be in ascending parameter order and contain any
parametric locations where one-sided derivative information may be
needed, such as slope, curvature, or higher-derivative changes, as
well as the end points for the curve (if they exist).
Endpoints are included as boundary markers even though they are only
one-sided breakpoints.
Interior breakpoints are the locations where two-sided derivative
information may differ.
Returns
-------
list[FloatScalar]
Parametric coordinates of any breakpoints.
"""
[docs]
def breakpoint_parameter_limits(
self,
*,
index: int,
) -> tuple[FloatScalar, FloatScalar]:
"""
Return parameter limits for one breakpoint.
Notes
-----
Endpoint breakpoints return the exact boundary parameter.
Interior breakpoints return nearby one-sided parameters chosen within
the neighboring breakpoint interval for the current generic
breakpoint-side implementation.
These limits exist to support the sampled fallback in the generic
``*_breakpoint`` methods and should not be treated as the primary
source of truth when a subclass can provide exact one-sided values
directly.
"""
breakpoints = self._breakpoints()
u_breakpoint = breakpoints[index]
if index == 0:
return u_breakpoint, u_breakpoint
if index == len(breakpoints) - 1:
return u_breakpoint, u_breakpoint
u_prev = breakpoints[index - 1]
u_next = breakpoints[index + 1]
minus_step = as_float_scalar(
np.minimum(
_BREAKPOINT_PARAMETER_STEP,
0.5 * (u_breakpoint - u_prev),
)
)
plus_step = as_float_scalar(
np.minimum(
_BREAKPOINT_PARAMETER_STEP,
0.5 * (u_next - u_breakpoint),
)
)
u_minus = as_float_scalar(u_breakpoint - minus_step)
u_plus = as_float_scalar(u_breakpoint + plus_step)
return (
u_minus,
u_plus,
)
#
# Parameter validation
#
def _validate_u(self, u: FloatInput) -> FloatArray:
"""
Validate and normalize native curve parameter values.
Parameters
----------
u : FloatInput
Native curve parameter values.
Returns
-------
FloatArray
Normalized ``float64`` array of native parameter values.
Raises
------
ValueError
If any value lies outside the curve breakpoint domain.
"""
u_array = as_float_array(u)
breakpoints = self._breakpoints()
u_min = breakpoints[0]
u_max = breakpoints[-1]
if (u_array < u_min - ROOT_ABS_TOLERANCE).any() or (
u_array > u_max + ROOT_ABS_TOLERANCE
).any():
msg = (
"Invalid curve parameter provided. "
f"Valid range is {u_min:.6g} <= u <= {u_max:.6g}."
)
raise ValueError(msg)
return u_array
def _validate_s(self, s: FloatInput) -> FloatArray:
"""
Validate and normalize arc-length parameter values.
Parameters
----------
s : FloatInput
Arc lengths measured from the minimum native parameter.
Returns
-------
FloatArray
Normalized ``float64`` array of arc-length values.
Raises
------
ValueError
If any value lies outside the curve arc-length domain.
"""
s_array = as_float_array(s)
total_length = self.length
if (s_array > total_length + ROOT_ABS_TOLERANCE).any() or (
s_array < -ROOT_ABS_TOLERANCE
).any():
msg = (
"Invalid arc length provided. "
f"Valid range is 0 <= s <= {total_length:.6g}."
)
raise ValueError(msg)
return s_array
#
# Arc-length interface
#
[docs]
def arc_length(self, u_s: FloatScalar, u_e: FloatInput) -> FloatArray:
"""
Calculate the arc-length distance between two points on surface.
Parameters
----------
u_s : FloatScalar
Start point of distance calculation.
u_e : FloatInput
End point of distance calculation.
Returns
-------
FloatArray
Distance from start point to end point.
"""
def fun(u: FloatScalar) -> FloatScalar:
xu, yu = self.xy_u(u)
return as_float_scalar(np.sqrt(xu**2 + yu**2))
u_begin = as_float_scalar(u_s)
u_end_array = as_float_array(u_e)
arc_length_array = np.empty_like(u_end_array)
flat_u_end = u_end_array.ravel()
flat_arc_length = arc_length_array.ravel()
for index, u_end in enumerate(flat_u_end):
segment_ends = self._arc_length_segment_ends(
u_begin=u_begin,
u_end=as_float_scalar(u_end),
)
segment_start = u_begin
arc_length_total = 0.0
for t_end in segment_ends:
arc_length_total += quad(
fun,
segment_start,
t_end,
epsabs=QUAD_ABS_TOLERANCE,
epsrel=QUAD_REL_TOLERANCE,
limit=QUAD_LIMIT,
)[0]
segment_start = as_float_scalar(t_end)
flat_arc_length[index] = arc_length_total
return arc_length_array
[docs]
def u_from_s(self, s: FloatInput) -> FloatArray:
"""
Calculate the parametric value for arc-length provided.
Parameters
----------
s : FloatInput
Arc-length location of point.
Returns
-------
FloatArray
Parametric value matching the normalized shape of ``s``.
Raises
------
ValueError
When arc-length provided lies outside the curve arc-length range.
"""
s_array = self._validate_s(s)
u_array = np.empty_like(s_array)
breakpoints = self._breakpoints()
u_min = breakpoints[0]
u_max = breakpoints[-1]
total_length = self.length
def solve(s_target: FloatScalar) -> FloatScalar:
if np.abs(s_target) < ROOT_ABS_TOLERANCE:
return u_min
if np.abs(s_target - total_length) < ROOT_ABS_TOLERANCE:
return u_max
def fun(u: FloatScalar, s: FloatScalar) -> FloatScalar:
return as_float_scalar(self.arc_length(u_min, u)) - s
root = root_scalar(
fun,
args=(s_target,),
bracket=(u_min, u_max),
xtol=ROOT_ABS_TOLERANCE,
rtol=ROOT_ABS_TOLERANCE,
maxiter=ROOT_MAX_ITERATION,
)
return as_float_scalar(root.root)
flat_s = s_array.ravel()
flat_u = u_array.ravel()
for index, s_value in enumerate(flat_s):
s_target = as_float_scalar(s_value)
flat_u[index] = solve(s_target)
return u_array
[docs]
def xy_from_s(self, s: FloatInput) -> tuple[FloatArray, FloatArray]:
"""
Return curve coordinates at arc-length locations.
Parameters
----------
s : FloatInput
Arc length location of point.
Returns
-------
tuple[FloatArray, FloatArray]
``(x, y)`` coordinates matching the normalized shape of ``s``.
"""
u = self.u_from_s(s)
return self.xy_from_u(u)
#
# Arc-length shared helpers
#
def _arc_length_segment_ends(
self,
*,
u_begin: FloatScalar,
u_end: FloatScalar,
) -> list[FloatScalar]:
"""
Return the breakpoint-partitioned segment ends for one arc-length run.
Parameters
----------
u_begin : FloatScalar
Native starting parameter.
u_end : FloatScalar
Native ending parameter.
Returns
-------
list[FloatScalar]
Ordered sequence of segment end parameters that partitions the
interval from ``u_begin`` to ``u_end`` at every interior
breakpoint. The starting value is omitted because callers already
track the current segment start explicitly.
"""
breakpoints = self._breakpoints()
if u_begin <= u_end:
interior_breakpoints = [
value for value in breakpoints if u_begin < value < u_end
]
return [*interior_breakpoints, u_end]
interior_breakpoints = [
value for value in breakpoints if u_end < value < u_begin
]
return [*reversed(interior_breakpoints), u_end]
[docs]
def xy_s(self, s: FloatInput) -> tuple[FloatArray, FloatArray]:
"""
Calculate first derivatives at arc-length location.
Parameters
----------
s : FloatInput
Arc length location of point.
Returns
-------
tuple[FloatArray, FloatArray]
``(dx/ds, dy/ds)`` coordinates matching the normalized shape of
``s``.
Notes
-----
If ``s`` matches one of :meth:`arc_length_breakpoints` exactly, this
method returns the ``minus``-side derivative limit.
Subclasses should override :meth:`xy_s_breakpoint` when exact
one-sided breakpoint derivatives are available analytically.
"""
s_array = as_float_array(s)
x_s = as_float_array(np.empty_like(s_array))
y_s = as_float_array(np.empty_like(s_array))
analytic_mask = np.ones_like(s_array, dtype=bool)
breakpoints = as_float_array(self.arc_length_breakpoints())
flat_s = s_array.ravel()
flat_x_s = x_s.ravel()
flat_y_s = y_s.ravel()
flat_mask = analytic_mask.ravel()
for index, value in enumerate(flat_s):
breakpoint_index = self._arc_length_breakpoint_index(
as_float_scalar(value),
breakpoints=breakpoints,
)
if breakpoint_index is not None:
minus, _ = self.xy_s_breakpoint(index=breakpoint_index)
flat_x_s[index] = minus[0]
flat_y_s[index] = minus[1]
flat_mask[index] = False
if np.any(analytic_mask):
t = self.u_from_s(s_array[analytic_mask])
x_s[analytic_mask], y_s[analytic_mask] = self.tangent(t)
return x_s, y_s
[docs]
def xy_ss(self, s: FloatInput) -> tuple[FloatArray, FloatArray]:
"""
Calculate second derivatives at arc-length location.
Parameters
----------
s : FloatInput
Arc length location of point.
Returns
-------
tuple[FloatArray, FloatArray]
``(d^2x/ds^2, d^2y/ds^2)`` coordinates matching the normalized
shape of ``s``.
Notes
-----
If ``s`` matches one of :meth:`arc_length_breakpoints` exactly, this
method returns the ``minus``-side derivative limit.
Subclasses should override :meth:`xy_ss_breakpoint` when exact
one-sided breakpoint second derivatives are available analytically.
"""
s_array = as_float_array(s)
x_ss = as_float_array(np.empty_like(s_array))
y_ss = as_float_array(np.empty_like(s_array))
analytic_mask = np.ones_like(s_array, dtype=bool)
breakpoints = as_float_array(self.arc_length_breakpoints())
flat_s = s_array.ravel()
flat_x_ss = x_ss.ravel()
flat_y_ss = y_ss.ravel()
flat_mask = analytic_mask.ravel()
for index, value in enumerate(flat_s):
breakpoint_index = self._arc_length_breakpoint_index(
as_float_scalar(value),
breakpoints=breakpoints,
)
if breakpoint_index is not None:
minus, _ = self.xy_ss_breakpoint(index=breakpoint_index)
flat_x_ss[index] = minus[0]
flat_y_ss[index] = minus[1]
flat_mask[index] = False
if np.any(analytic_mask):
t = self.u_from_s(s_array[analytic_mask])
x_ss[analytic_mask], y_ss[analytic_mask] = self._xy_ss_from_t(t)
return x_ss, y_ss
[docs]
def xy_s_breakpoint(
self,
*,
index: int,
) -> CurveBreakpointSides:
"""
Return both sides of first derivatives at a breakpoint.
Parameters
----------
index : int
Index into :meth:`arc_length_breakpoints`.
Returns
-------
CurveBreakpointSides
``((x_s_minus, y_s_minus), (x_s_plus, y_s_plus))``.
Notes
-----
Endpoint breakpoints return the same boundary value for both entries.
This method is the exact-breakpoint contract that pairs with
:meth:`xy_s`.
Subclasses should override it whenever exact one-sided arc-length
derivatives are available.
The generic implementation evaluates nearby one-sided native-
parameter samples and therefore serves only as an approximation
fallback.
"""
t_minus, t_plus = self.breakpoint_parameter_limits(index=index)
minus_derivative = self.tangent(t_minus)
plus_derivative = self.tangent(t_plus)
minus_values = (
as_float_scalar(minus_derivative[0]),
as_float_scalar(minus_derivative[1]),
)
plus_values = (
as_float_scalar(plus_derivative[0]),
as_float_scalar(plus_derivative[1]),
)
return minus_values, plus_values
[docs]
def xy_ss_breakpoint(
self,
*,
index: int,
) -> CurveBreakpointSides:
"""
Return one-sided second derivatives at a breakpoint.
Parameters
----------
index : int
Index into :meth:`arc_length_breakpoints`.
Returns
-------
CurveBreakpointSides
``((x_ss_minus, y_ss_minus), (x_ss_plus, y_ss_plus))``.
Notes
-----
Endpoint breakpoints return the same boundary value for both entries.
This method is the exact-breakpoint contract that pairs with
:meth:`xy_ss`.
Subclasses should override it whenever exact one-sided arc-length
second derivatives are available.
The generic implementation evaluates nearby one-sided native-
parameter samples and therefore serves only as an approximation
fallback.
"""
t_minus, t_plus = self.breakpoint_parameter_limits(index=index)
minus_values = self._xy_ss_from_t(t_minus)
plus_values = self._xy_ss_from_t(t_plus)
return (
as_float_scalar(minus_values[0]),
as_float_scalar(minus_values[1]),
), (
as_float_scalar(plus_values[0]),
as_float_scalar(plus_values[1]),
)
def _xy_ss_from_t(self, t: FloatInput) -> tuple[FloatArray, FloatArray]:
"""
Return arc-length second derivatives from native-parameter data.
Parameters
----------
t : FloatInput
Native curve parameter locations.
Returns
-------
tuple[FloatArray, FloatArray]
``(x_ss, y_ss)`` values obtained from the unit normal and scalar
curvature at ``t``.
"""
nx, ny = self.normal(t)
curvature = self.k(t)
return (
as_float_array(curvature * nx),
as_float_array(curvature * ny),
)
[docs]
def arc_length_breakpoints(self) -> list[FloatScalar]:
"""
Return the breakpoint locations in arc-length coordinates.
Returns
-------
list[FloatScalar]
Arc-length coordinates measured from the minimum native
parameter.
Notes
-----
These values include the two curve endpoints as boundary markers.
Interior breakpoints correspond to the native-parameter interior
breakpoints returned by :meth:`breakpoints`.
"""
u_min = self._breakpoints()[0]
return [
as_float_scalar(self.arc_length(u_min, value))
for value in self._breakpoints()
]
def _breakpoints(self) -> list[FloatScalar]:
"""
Return validated breakpoint values in strictly ascending order.
Raises
------
ValueError
If :meth:`breakpoints` omits the endpoints or does not return
strictly increasing values.
"""
breakpoint_values = self.breakpoints()
if not breakpoint_values:
msg = (
f"{type(self).__name__}.breakpoints() must include the "
"curve endpoints."
)
raise ValueError(msg)
breakpoint_array = as_float_array(breakpoint_values)
if np.any(np.diff(breakpoint_array) <= 0.0):
msg = (
f"{type(self).__name__}.breakpoints() must return values "
"in strictly ascending order."
)
raise ValueError(msg)
return [as_float_scalar(value) for value in breakpoint_array]
@staticmethod
def _arc_length_breakpoint_index(
value: FloatScalar,
*,
breakpoints: FloatArray,
) -> int | None:
"""
Return the matching arc-length breakpoint index within tolerance.
Parameters
----------
value : FloatScalar
Arc-length coordinate to classify.
breakpoints : FloatArray
Precomputed arc-length breakpoint coordinates in ascending order.
Returns
-------
int | None
Index of the matching breakpoint when ``value`` lies within the
root-tolerance comparison band, otherwise ``None``.
"""
matches = np.flatnonzero(
np.isclose(
breakpoints,
value,
atol=ROOT_ABS_TOLERANCE,
rtol=0.0,
)
)
if matches.size == 0:
return None
return int(matches[0])