diff --git a/CHANGELOG.md b/CHANGELOG.md index 7c8b047..e45226d 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -7,6 +7,21 @@ semantic versioning. ## Unreleased ### Added +- **`tensorcad.nurbs`: piecewise (break-aware) degree elevation** (`NURBSCurve.split_at_discontinuities`, `NURBSCurve.elevate_degree_piecewise`, + `docs/nurbs_piecewise_degree_elevation.md`, `examples/nurbs_piecewise_elevation/`, + 17 tests). Degree elevation already refused a curve with an interior knot of + multiplicity `degree + 1` -- a break, not a `C0` join -- and told the caller to + split the curve first, but nothing in the package could split it: a composite + profile (polyline plus arc, a duplicated seam, two concatenated patches) could not + reach the elevated representation at all. `split_at_discontinuities` returns the + maximal `C0`-or-better pieces as a pure slice of the control polygon (nothing + inserted, nothing re-solved, so the split is bit-exact and keeps every interior + knot and rational weight; a smooth curve returns itself), and + `elevate_degree_piecewise(times)` composes it with the exact homogeneous elevation. + Evidence on the exact rational quarter circle plus straight segment: radius + `2 m + 8.9e-16 m`, total length `pi + 2 sqrt 2` to `8.9e-16 m`, point deviation + `<= 6.3e-16 m`, the 45 deg kink and the two jumps of a three-piece fixture + preserved, control-point growth minimal. - **`tensorcad.girder`: hull-girder elastic and plastic section properties** (`tensorcad/girder.py`, `docs/girder_section.md`, `examples/girder_section/`, `tools/girder_section_doc.py`, 25 tests). Geometry could produce a hull and the hydrostatics could say what it displaces, but nothing diff --git a/docs/nurbs_degree_elevation.md b/docs/nurbs_degree_elevation.md index f7a3c02..58412c0 100644 --- a/docs/nurbs_degree_elevation.md +++ b/docs/nurbs_degree_elevation.md @@ -132,7 +132,11 @@ with the degree.* documentation already states. A curve with an interior knot of multiplicity `degree + 1` is genuinely discontinuous (its two sides share no control point), so the single-curve collocation matrix is singular; the method refuses with a named - remedy ("split it at that knot first") rather than returning a wrong curve. + remedy rather than returning a wrong curve. That remedy now exists as real code: + `NURBSCurve.split_at_discontinuities` and `NURBSCurve.elevate_degree_piecewise` + (`docs/nurbs_piecewise_degree_elevation.md`) split a broken curve into its exact + pieces and elevate each -- the elevated result is a tuple, because a break stays a + break at every degree. 6. **Cost.** The implementation solves a dense `count x count` system with one elimination per call. For the curves in this library (tens of control points) that is microseconds; a production kernel would exploit the band explicitly. The diff --git a/docs/nurbs_piecewise_degree_elevation.md b/docs/nurbs_piecewise_degree_elevation.md new file mode 100644 index 0000000..e35b83b --- /dev/null +++ b/docs/nurbs_piecewise_degree_elevation.md @@ -0,0 +1,164 @@ +# Piecewise (break-aware) degree elevation + +## 1. Case introduction + +`NURBSCurve.elevate_degree` (`docs/nurbs_degree_elevation.md`) raises the degree of a +curve exactly, but it refuses a curve with an interior knot of multiplicity +`degree + 1` -- a genuine **break**, where the two sides share no control point and the +collocation matrix is singular by construction. The refusal was honest but it named a +remedy that did not exist in the package: "split it at that knot first" had no +implementation, so a user holding a broken composite -- the natural output of a +polyline-plus-arc profile, a duplicated seam, or an assembly where two patches were +concatenated -- had no way to reach the elevated representation at all. + +This increment lands the missing half: `NURBSCurve.split_at_discontinuities()` and +`NURBSCurve.elevate_degree_piecewise(times)`. The reference object is a degree-2 curve +that is a real composite: + +| item | value | +|---|---| +| piece A | exact rational quadratic quarter circle, radius `R = 2.0 m`, control points `(2,0), (2,2), (0,2)`, weights `1, sqrt(2)/2, 1` | +| piece B | straight segment `(0,2) -> (-2,0)` as a degenerate quadratic, weights `1, 1, 1` | +| knot vector | `0,0,0, 1,1,1, 2,2,2` (multiplicity 3 = `degree + 1` at `u = 1`: a break) | +| control points / degree | 6 / 2 | +| analytic arc length | `pi R / 2 + 2 sqrt 2 = 5.970019778335983 m` | +| analytic kink at the join | 45 deg between the one-sided tangents | + +Second fixture, for jumps: a degree-3 curve with **two** breaks whose sides do not even +share a point (`knots 0,0,0,0, 0.4, 1^4, 1.5^4, 2,2,2,2`), which is what a collection of +unconnected patches looks like when someone concatenates their control polygons. + +## 2. Method + +The insight that makes the split free: at a break the multiplicity is already +`degree + 1`, so the curve *is* the direct sum of its pieces. No knot insertion, no +re-fitting, no solve is needed. + +1. Find the break parameters: interior knots whose multiplicity exceeds `degree`. +2. For the interval between consecutive bounds (domain start, breaks, domain end), + the piece's knot vector is `[a]*(degree+1) + (original knots strictly inside) + + [b]*(degree+1)`, and its control points are the contiguous window + `points[s_a - degree : ...]` where `s_a` is the index of the **last** occurrence of + `a` in the original knot vector. The window length follows from + `count = len(piece_knots) - degree - 1`, which is exactly one point per basis + function that is non-zero on the interval. + +Consequences, all of them testable rather than asserted: + +* **Zero arithmetic.** Because the pieces are slices of the original polygon, the + geometry is bit-for-bit preserved: the deviation between the whole curve and its + pieces is `0.0`, not `1e-16`, and rational weights travel untouched, so the arc piece + is still the exact circle. +* **Nothing is inserted.** Interior knots *inside* a piece (a `C0`/`C1`/`C2` join that + is not a break) keep their multiplicity, so the split does not degrade continuity or + grow the representation. A curve with no break returns `(self,)` -- the identity. +* **Parameters are preserved.** The i-th piece lives on the original knot interval, so + no reparameterisation happens and a deviation can be measured at the same `u`. + +A third, very practical fixture is checked in the tests: a chine profile stored as one +degree-1 polygon with duplicated vertices (`0,0, 1,1, 2,2, 3,3` knots) splits into its +three segments, whose lengths add to exactly `6.0 m` (`2 + 2 + 2`), and one elevation +puts each middle control point exactly on the analytic segment midpoint -- the degree-2 +form of a straight segment is the original endpoints plus the midpoint. + +`elevate_degree_piecewise(times)` is then just +`tuple(piece.elevate_degree(times) for piece in split_at_discontinuities())`. Each piece +is `C0`-or-better, so the existing exact algorithm (ISO 10303-42 / Piegl & Tiller A5.9, +homogeneous collocation at Greville abscissae) applies unchanged, and the result is a +tuple rather than a curve because a break stays a break at every degree -- no single +spline represents it. + +The whole broken curve is still point-evaluable (that is how it is visualised), but it is +not differentiable across a break: the derivative net divides by a zero-length knot +interval, so `arc_length`, `curvature` and `tangent` raise on the composite and need the +pieces. That asymmetry is the practical argument for the split. + +Reproduce everything below with + +```bash +PYTHONPATH=. python examples/nurbs_piecewise_elevation/run_example.py --out artifacts/m15428 +python -m pytest tests/test_nurbs_piecewise_degree_elevation.py -q +``` + +## 3. Results + +Elevation of the composite, `times = 0..3` (generated into +`artifacts/m15428/table.md` and `metrics.json`; not hand-edited): + +| times | degree | control points | max |elev - original| [m] | radius error [m] | arc length [m] | line length [m] | total - analytic [m] | kink angle [deg] | +|---|---|---|---|---|---|---|---|---| +| 0 | [2, 2] | [3, 3] | 0.000e+00 | 4.441e-16 | 3.141592653589793 | 2.828427124746191 | +1.776e-15 | 45.000000000000 | +| 1 | [3, 3] | [4, 4] | 6.280e-16 | 8.882e-16 | 3.141592653589794 | 2.828427124746191 | +1.776e-15 | 45.000000000000 | +| 2 | [4, 4] | [5, 5] | 4.965e-16 | 4.441e-16 | 3.141592653589793 | 2.828427124746191 | +1.776e-15 | 45.000000000000 | +| 3 | [5, 5] | [6, 6] | 4.965e-16 | 8.882e-16 | 3.141592653589793 | 2.828427124746191 | +8.882e-16 | 45.000000000000 | + +Analytic references: arc `pi = 3.141592653589793`, line `2 sqrt 2 = 2.8284271247461903`, +total `5.970019778335983`. Control-point growth is the minimum +`times * (#spans of the piece)` (1 per single-span piece here), and the radius error is +the circle anchor: the raised arc is still exactly the circle of radius 2, because the +elevation solve runs in homogeneous coordinates. + +Second fixture, the jumps (a break stays a break, elevations do not close it): + +| break | jump before [m] | jump after x2 [m] | +|---|---|---| +| 1.0 | 1.414213562373095 | 1.414213562373095 | +| 1.5 | 3.162277660168380 | 3.162277660168380 | + +The two jumps are `sqrt 2` and `sqrt 10`; the first piece keeps its interior `C2` knot +(`0.4`, multiplicity 1) verbatim through the split, and grows by exactly 2 control +points (two spans) under one elevation. + +![geometry of the composite and its elevated pieces; point deviation at machine precision](nurbs_piecewise_degree_elevation/figure.png) + +Left: the composite sampled piece by piece -- original (degree 2), one elevation +(degree 3) and three elevations (degree 5) are visually indistinguishable, which is the +claim. Right: `|elevated(u) - original(u)|` per piece on a log axis, flat at +`1e-16 .. 1e-15` on both sides of the break at `u = 1`. + +## 4. Notes, limits and things that bit + +* **Multiple spans per piece without refinement.** The first draft inserted knots to + make every interior join `C0` before slicing. That was unnecessary (the window + identity holds for any multiplicity) and it silently degraded `C1`/`C2` joins and grew + the representation, so it was removed; the final code inserts nothing. +* **Evaluation on a break is one-sided.** `find_span` returns the last span, so + `evaluate(1.0)` on the composite returns piece B's start. Sampling a deviation + *across* a break therefore shows the jump (1.41 m / 3.16 m in the second fixture), + which is geometry, not error. The tests and the example sample strictly between the + knots. +* **Rational exactness must be checked on observables.** `weights` equality, `|P(u)| = R` + and the analytic arc length are checked separately, because a wrong elevation could + keep the endpoints and still bulge. +* **Scope.** Curves only. The tensor-product counterpart has the same hole per + direction (`NURBSSurface.elevate_degree` routes each row/column through the same + helper), but a surface has to be split at an isoparametric break in one direction and + at an arbitrary parameter in the other, which is a separate increment; trimming and + stitching remain outside this module. +* **Cost.** The split is `O(#knots)` slicing; the elevation cost per piece is unchanged + (dense elimination, microseconds at these sizes). + +## 5. Error judgement against references + +| quantity | measured | reference | error | within 3 %? | +|---|---:|---:|---:|---| +| arc piece radius, raised x1 / x3 | 2.0 m + 8.9e-16 m | 2.0 m (analytic circle) | 4.4e-16 relative | **yes** | +| arc piece length, raised x1 / x3 | 3.141592653589793 m | `pi` m | < 1e-15 relative | **yes** | +| straight piece length | 2.828427124746191 m | `2 sqrt 2 = 2.8284271247461903` m | 8.9e-16 absolute | **yes** | +| total composite length after x3 | 5.970019778335985 m | `pi + 2 sqrt 2 = 5.970019778335983` m | 8.9e-16 absolute | **yes** | +| point deviation, original vs each piece | 0.0 m | 0 (pure slice) | exact | **yes** | +| point deviation, elevated vs original | <= 6.3e-16 m | 0 (identity) | absolute 6.3e-16 m | **yes** | +| kink angle after x3 | 45.000000000000 deg | 45 deg (bisector of the axes) | < 1e-12 deg | **yes** | +| jumps of the second fixture after x2 | `sqrt 2`, `sqrt 10` | same values before elevation | < 1e-14 absolute | **yes** | + +Every reference here is an exact analytic identity or an exact structural invariant, +not an engineering approximation, so the meaningful statement is "at double-precision +round-off": all errors are at least twelve orders of magnitude below the 3 % threshold. +No hydrodynamics, no structural suitability: this is a representation change, and the +only claims made are geometric. + +Regression evidence: `artifacts/m15428/red_new_test_before_implementation.txt` +(15 failed / 1 passed before the methods existed) and +`artifacts/m15428/green_new_test_after_implementation.txt` (67 passed across +`test_nurbs.py`, `test_nurbs_degree_elevation.py`, +`test_nurbs_surface_degree_elevation.py`, `test_nurbs_piecewise_degree_elevation.py`). \ No newline at end of file diff --git a/docs/nurbs_piecewise_degree_elevation/figure.png b/docs/nurbs_piecewise_degree_elevation/figure.png new file mode 100644 index 0000000..f43625a Binary files /dev/null and b/docs/nurbs_piecewise_degree_elevation/figure.png differ diff --git a/examples/nurbs_piecewise_elevation/run_example.py b/examples/nurbs_piecewise_elevation/run_example.py new file mode 100644 index 0000000..0d977b8 --- /dev/null +++ b/examples/nurbs_piecewise_elevation/run_example.py @@ -0,0 +1,246 @@ +"""Piecewise (break-aware) degree elevation: evidence generator for the docs. + +The object under test is a degree-2 curve with a genuine break: an exact rational +quarter circle of radius 2 followed by a straight segment, joined with multiplicity +``degree + 1 = 3``. The script measures, for ``times = 0..3``: + +* how many pieces the split produces, their degrees and control-point counts; +* the maximum point deviation between each elevated piece and the original curve + (sampled between the knots, never on a break, where evaluation is one-sided); +* the deviation of the raised arc from the analytic radius (``|P(u)| - R``); +* the arc length of the pieces against the analytic ``pi R / 2 + 2 sqrt 2``; +* the kink angle between the two one-sided tangents at the break; +* for a three-piece curve whose joins do not even share a point, the size of each + jump before and after elevation (a break stays a break). + +Outputs (``--out``): ``metrics.json``, ``table.md`` (the table quoted in +``docs/nurbs_piecewise_degree_elevation.md``) and ``figure.png``. + +Run from the repository root:: + + python examples/nurbs_piecewise_elevation/run_example.py --out artifacts/m15428 +""" +from __future__ import annotations + +import argparse +import json +import math +from itertools import pairwise +from pathlib import Path + +import matplotlib + +matplotlib.use("Agg") + +import matplotlib.pyplot as plt # backend must be set before this import + +from tensorcad.nurbs import NURBSCurve, NURBSError + +RADIUS = 2.0 +ARC_LENGTH = math.pi * RADIUS / 2.0 +LINE_LENGTH = math.hypot(2.0, 2.0) +TOTAL_LENGTH = ARC_LENGTH + LINE_LENGTH + + +def composite_curve() -> NURBSCurve: + """Exact rational quarter circle (R = 2) + straight segment, joined at a break.""" + return NURBSCurve( + degree=2, + knots=(0.0, 0.0, 0.0, 1.0, 1.0, 1.0, 2.0, 2.0, 2.0), + control_points=( + (2.0, 0.0, 0.0), + (2.0, 2.0, 0.0), + (0.0, 2.0, 0.0), + (0.0, 2.0, 0.0), + (-1.0, 1.0, 0.0), + (-2.0, 0.0, 0.0), + ), + weights=(1.0, math.cos(math.pi / 4.0), 1.0, 1.0, 1.0, 1.0), + ) + + +def jump_curve() -> NURBSCurve: + """Degree 3, three pieces, two breaks that do not share a point.""" + knots = (0.0, 0.0, 0.0, 0.0, 0.4, 1.0, 1.0, 1.0, 1.0, 1.5, 1.5, 1.5, 1.5, 2.0, 2.0, 2.0, 2.0) + points = tuple((float(index), float((index * index) % 5), 0.0) for index in range(len(knots) - 4)) + return NURBSCurve(degree=3, knots=knots, control_points=points) + + +def breaks_of(curve: NURBSCurve) -> list[float]: + low, high = curve.domain + values = sorted({value for value in curve.knots if low < value < high}) + return [value for value in values if sum(1 for knot in curve.knots if knot == value) > curve.degree] + + +def sample(segment: NURBSCurve, per_span: int = 61) -> list[tuple[float, tuple[float, float, float]]]: + low, high = segment.domain + distinct = sorted({value for value in segment.knots if low <= value <= high}) + parameters = [ + left + (right - left) * index / per_span + for left, right in pairwise(distinct) + for index in range(per_span + 1) + ] + return [(u, segment.evaluate(u)) for u in parameters] + + +def max_deviation(curve: NURBSCurve, piece: NURBSCurve) -> float: + worst = 0.0 + for u, point in sample(piece): + reference = curve.evaluate(u) + worst = max(worst, math.dist(point, reference)) + return worst + + +def elevation_report(times: int) -> dict[str, object]: + curve = composite_curve() + pieces = curve.split_at_discontinuities() + raised = curve.elevate_degree_piecewise(times) if times else pieces + deviations = [max_deviation(curve, piece) for piece in raised] + arc = raised[0] + radius_error = max( + abs(math.hypot(*arc.evaluate(index / 200.0)[:2]) - RADIUS) for index in range(201) + ) + lengths = [piece.arc_length() for piece in raised] + left, right = raised[0].tangent(1.0), raised[1].tangent(1.0) + cosine = max(-1.0, min(1.0, sum(a * b for a, b in zip(left, right, strict=True)))) + return { + "times": times, + "pieces": len(raised), + "degrees": [piece.degree for piece in raised], + "control_points": [len(piece.control_points) for piece in raised], + "max_deviation": max(deviations), + "radius_error": radius_error, + "arc_length": lengths[0], + "line_length": lengths[1], + "total_length": sum(lengths), + "total_length_error": sum(lengths) - TOTAL_LENGTH, + "kink_angle_deg": math.degrees(math.acos(cosine)), + } + + +def jump_report() -> list[dict[str, float]]: + curve = jump_curve() + pieces = curve.split_at_discontinuities() + raised = curve.elevate_degree_piecewise(2) + rows = [] + for index in range(len(pieces) - 1): + before = math.dist( + pieces[index].evaluate(pieces[index].domain[1]), + pieces[index + 1].evaluate(pieces[index + 1].domain[0]), + ) + after = math.dist( + raised[index].evaluate(raised[index].domain[1]), + raised[index + 1].evaluate(raised[index + 1].domain[0]), + ) + rows.append({"break": pieces[index].domain[1], "jump_before": before, "jump_after": after}) + return rows + + +def figure(out: Path, reports: list[dict[str, object]]) -> None: + curve = composite_curve() + pieces = curve.split_at_discontinuities() + fig, axes = plt.subplots(1, 2, figsize=(11.0, 4.6)) + + ax = axes[0] + designs = [ + ("original (break at u=1)", pieces, "-"), + ("elevated x1 (degree 3)", curve.elevate_degree_piecewise(1), "--"), + ("elevated x3 (degree 5)", curve.elevate_degree_piecewise(3), ":"), + ] + for label, segments, style in designs: + for segment in segments: + points = sample(segment, per_span=24) + ax.plot([p[0] for _, p in points], [p[1] for _, p in points], style, lw=1.6, label=label) + label = None # one legend entry per design + ax.plot([-2, 2], [0, 0], "k.", ms=3, alpha=0.25) + ax.set_aspect("equal") + ax.set_xlabel("x [m]") + ax.set_ylabel("y [m]") + ax.set_title("quarter circle R=2 + straight segment, split at the break") + ax.legend(loc="lower left", fontsize=7) + ax.grid(alpha=0.3) + + ax = axes[1] + parameters = [index / 200.0 for index in range(201)] + for report in reports: + raised = curve.elevate_degree_piecewise(int(report["times"])) if report["times"] else pieces + deviations = [] + for u in parameters: + segment = raised[0] if u <= 1.0 else raised[1] + deviations.append(max(math.dist(segment.evaluate(u), curve.evaluate(u)), 1e-17)) + ax.semilogy(parameters, deviations, lw=1.4, label=f"degree {report['degrees'][0]} (x{report['times']})") + ax.axvline(1.0, color="k", lw=0.8, ls="--") + ax.set_xlabel("parameter u") + ax.set_ylabel("|elevated - original| [m]") + ax.set_title("point deviation (machine precision), break at u=1") + ax.legend(fontsize=7) + ax.grid(alpha=0.3, which="both") + + fig.tight_layout() + fig.savefig(out / "figure.png", dpi=160) + plt.close(fig) + + +def main() -> int: + parser = argparse.ArgumentParser(description=__doc__) + parser.add_argument("--out", type=Path, default=Path("artifacts/m15428")) + arguments = parser.parse_args() + arguments.out.mkdir(parents=True, exist_ok=True) + + reports = [elevation_report(times) for times in range(4)] + rows = jump_report() + + try: + composite_curve().elevate_degree(1) + refusal = "NOT RAISED" + except NURBSError as error: + refusal = str(error) + + payload = { + "case": { + "degree": 2, + "pieces": 2, + "break_multiplicity": 3, + "radius": RADIUS, + "analytic_arc_length": ARC_LENGTH, + "analytic_line_length": LINE_LENGTH, + "analytic_total_length": TOTAL_LENGTH, + }, + "elevation": reports, + "jumps": rows, + "single_curve_refusal": refusal, + } + (arguments.out / "metrics.json").write_text(json.dumps(payload, indent=2) + "\n") + + header = ( + "| times | degree | control points | max |elev - original| [m] | radius error [m] " + "| arc length [m] | line length [m] | total - analytic [m] | kink angle [deg] |\n" + "|---|---|---|---|---|---|---|---|---|\n" + ) + body = "".join( + f"| {report['times']} | {report['degrees']} | {report['control_points']} " + f"| {report['max_deviation']:.3e} | {report['radius_error']:.3e} " + f"| {report['arc_length']:.15f} | {report['line_length']:.15f} " + f"| {report['total_length_error']:+.3e} | {report['kink_angle_deg']:.12f} |\n" + for report in reports + ) + jump_table = ( + "\n| break | jump before [m] | jump after x2 [m] |\n|---|---|---|\n" + + "".join( + f"| {row['break']} | {row['jump_before']:.15f} | {row['jump_after']:.15f} |\n" + for row in rows + ) + ) + (arguments.out / "table.md").write_text(header + body + jump_table) + + figure(arguments.out, reports) + + print(header + body + jump_table) + print(f"analytic total length = {TOTAL_LENGTH:.15f} m") + print(f"single-curve elevation refused with: {refusal}") + print(f"artifacts written to {arguments.out}") + return 0 + + +if __name__ == "__main__": + raise SystemExit(main()) diff --git a/tensorcad/nurbs.py b/tensorcad/nurbs.py index 9528f3b..221d0ac 100644 --- a/tensorcad/nurbs.py +++ b/tensorcad/nurbs.py @@ -68,6 +68,7 @@ import math from collections.abc import Sequence from dataclasses import dataclass +from itertools import pairwise from pathlib import Path from typing import Any @@ -389,8 +390,10 @@ def _elevate_homogeneous( raise NURBSError( f"degree elevation needs a C0-or-better {label}; knot " f"{disconnected[0]!r} has multiplicity {multiplicity[disconnected[0]]} " - f"> degree {degree}, so the {label} is discontinuous there -- split it " - "at that knot first" + f"> degree {degree}, so the {label} is discontinuous there -- for a curve " + "use NURBSCurve.split_at_discontinuities(), or " + "NURBSCurve.elevate_degree_piecewise(times) which splits and elevates in " + "one step" ) elevated_knots: list[float] = [low] * (target + 1) @@ -765,6 +768,74 @@ def elevate_degree(self, times: int = 1) -> NURBSCurve: weights=weights, ) + # -- piecewise (break-aware) operations --------------------------------- + def split_at_discontinuities(self) -> tuple[NURBSCurve, ...]: + """The maximal ``C0``-or-better pieces, split where the curve is broken. + + An interior knot of multiplicity ``degree + 1`` is not a ``C0`` join: the two + sides share no control point, so a single-curve collocation matrix is singular + there and :meth:`elevate_degree` (rightly) refuses. This method is the + splitting step that refusal asks for; :meth:`elevate_degree_piecewise` + composes the two. + + The pieces are exact, not refitted: because a break already carries + multiplicity ``degree + 1``, the curve is the direct sum of its pieces -- the + i-th piece's control points are a contiguous window of the original polygon and + its knots the corresponding sub-range of the original knot vector. Nothing is + inserted, removed or re-solved, so every interior knot (and therefore every + ``C0``/``C1`` join *inside* a piece) is preserved verbatim, parameters are + unchanged (the i-th piece is parameterised on the knot interval it came from), + and rational weights travel untouched. A curve with no break returns + ``(self,)``: the identity, not a copy. + + Only true breaks (multiplicity ``degree + 1``) split the curve. Note that the + *whole* broken curve is point-evaluable but not differentiable across a break + (the derivative net divides by a zero-length knot interval), which is exactly + why ``arc_length`` and ``curvature`` need the pieces this method returns. + """ + degree = self.degree + low, high = self.domain + interior = sorted({value for value in self.knots if low < value < high}) + breaks = [value for value in interior if _multiplicity(self.knots, value) > degree] + if not breaks: + return (self,) + + knots = self.knots + bounds = [low, *breaks, high] + pieces: list[NURBSCurve] = [] + for left, right in pairwise(bounds): + last_left = max(index for index, value in enumerate(knots) if value == left) + piece_knots = [left] * (degree + 1) + piece_knots.extend(value for value in knots if left < value < right) + piece_knots.extend([right] * (degree + 1)) + count = len(piece_knots) - degree - 1 + start = last_left - degree + pieces.append( + NURBSCurve( + degree=degree, + knots=tuple(piece_knots), + control_points=self.control_points[start : start + count], + weights=self.weights[start : start + count], + ) + ) + return tuple(pieces) + + def elevate_degree_piecewise(self, times: int = 1) -> tuple[NURBSCurve, ...]: + """Elevate every continuous piece of a curve that is allowed to have breaks. + + The piecewise counterpart of :meth:`elevate_degree`: split at every + multiplicity ``degree + 1`` knot (:meth:`split_at_discontinuities`) and elevate + each piece. The result is a tuple of segments, because a break stays a break at + every degree -- no single spline of degree ``degree + times`` represents a + broken curve. Every segment is exact, so the union traces the original curve + point-for-point while each segment is now usable by operations that need one + connected control polygon (patch-degree alignment against a neighbour, + curvature, collocation). + """ + if not isinstance(times, int) or isinstance(times, bool) or times < 0: + raise NURBSError(f"times must be a non-negative integer, got {times!r}") + return tuple(piece.elevate_degree(times) for piece in self.split_at_discontinuities()) + # -- continuity --------------------------------------------------------- def continuity(self, order: int = 2) -> dict[str, Any]: records = curve_derivative_continuity(self, order) diff --git a/tests/test_nurbs_piecewise_degree_elevation.py b/tests/test_nurbs_piecewise_degree_elevation.py new file mode 100644 index 0000000..263f00e --- /dev/null +++ b/tests/test_nurbs_piecewise_degree_elevation.py @@ -0,0 +1,291 @@ +"""Piecewise (break-aware) degree elevation: splitting, exactness, analytic anchors. + +The reference object is a degree-2 curve with a real break in it: an exact rational +quarter circle of radius 2 (control polygon ``(2,0), (2,2), (0,2)``, weights +``1, sqrt(2)/2, 1``) followed by a straight segment to ``(-2, 0)``. The join carries +multiplicity ``degree + 1 = 3``, so it is not a ``C0`` join but a break: the two sides +share no control point, the tangent turns by 45 degrees across it, and every +single-curve operation that differentiates (``arc_length``, ``curvature``, +``elevate_degree``) is singular there. + +Analytic anchors (no reference implementation involved): + +* the arc piece is a circle of radius 2, so ``|P(u)| == 2`` for every ``u``; +* the arc piece has length ``pi R / 2 = pi`` and the straight piece length ``2 sqrt 2``; +* the straight piece is the linear map ``u -> (-2u, 2 - 2u)`` on its own parameter + interval, derived here by hand from the Bernstein form; +* the kink angle between the two one-sided tangents is 45 degrees. +""" +from __future__ import annotations + +import math + +import pytest + +from tensorcad.nurbs import NURBSCurve, NURBSError, circle_curve + +RADIUS = 2.0 +ARC_LENGTH = math.pi * RADIUS / 2.0 # pi +LINE_LENGTH = math.hypot(2.0, 2.0) # 2 sqrt 2 +TOTAL_LENGTH = ARC_LENGTH + LINE_LENGTH # pi + 2 sqrt 2 = 5.970019778335983 + + +def composite_curve() -> NURBSCurve: + """Quarter circle followed by a straight segment, joined at a genuine break.""" + half_angle_weight = math.cos(math.pi / 4.0) # sqrt(2)/2: exact rational quarter circle + return NURBSCurve( + degree=2, + knots=(0.0, 0.0, 0.0, 1.0, 1.0, 1.0, 2.0, 2.0, 2.0), + control_points=( + (2.0, 0.0, 0.0), + (2.0, 2.0, 0.0), + (0.0, 2.0, 0.0), + (0.0, 2.0, 0.0), + (-1.0, 1.0, 0.0), + (-2.0, 0.0, 0.0), + ), + weights=(1.0, half_angle_weight, 1.0, 1.0, 1.0, 1.0), + ) + + +def jump_curve() -> NURBSCurve: + """Degree 3, three pieces, two breaks: the joins do not even share a point.""" + knots = (0.0, 0.0, 0.0, 0.0, 0.4, 1.0, 1.0, 1.0, 1.0, 1.5, 1.5, 1.5, 1.5, 2.0, 2.0, 2.0, 2.0) + points = tuple((float(index), float((index * index) % 5), 0.0) for index in range(len(knots) - 4)) + return NURBSCurve(degree=3, knots=knots, control_points=points) + + +def interior_parameters(curve: NURBSCurve, samples: int = 97) -> list[float]: + """Parameters inside the domain, never on a break (evaluation there is one-sided).""" + low, high = curve.domain + breaks = { + value + for value in curve.knots + if low < value < high + and sum(1 for knot in curve.knots if knot == value) > curve.degree + } + return [ + low + (high - low) * (index + 0.5) / samples + for index in range(samples) + if all(abs(low + (high - low) * (index + 0.5) / samples - b) > 1e-12 for b in breaks) + ] + + +# --------------------------------------------------------------------------- +# splitting +# --------------------------------------------------------------------------- +def test_split_returns_one_piece_per_continuous_piece() -> None: + pieces = composite_curve().split_at_discontinuities() + assert len(pieces) == 2 + assert [piece.degree for piece in pieces] == [2, 2] + assert [piece.domain for piece in pieces] == [(0.0, 1.0), (1.0, 2.0)] + # a pure slice: three control points per single-span quadratic piece + assert [len(piece.control_points) for piece in pieces] == [3, 3] + # weights travel untouched, so the arc piece is still the exact rational arc + assert pieces[0].weights == (1.0, math.cos(math.pi / 4.0), 1.0) + assert pieces[1].weights == (1.0, 1.0, 1.0) + assert pieces[0].rational and not pieces[1].rational + + +def test_split_pieces_reproduce_every_interior_point() -> None: + curve = composite_curve() + pieces = curve.split_at_discontinuities() + deviation = 0.0 + for piece in pieces: + for u in interior_parameters(curve): + if piece.domain[0] <= u <= piece.domain[1]: + deviation = max(deviation, math.dist(curve.evaluate(u), piece.evaluate(u))) + assert deviation < 1e-14 + + +def test_split_is_a_pure_slice_on_a_three_piece_curve() -> None: + """Interior knots *inside* a piece are reproduced verbatim -- nothing is inserted.""" + curve = jump_curve() + pieces = curve.split_at_discontinuities() + assert len(pieces) == 3 + assert pieces[0].knots == (0.0, 0.0, 0.0, 0.0, 0.4, 1.0, 1.0, 1.0, 1.0) + assert pieces[1].knots == (1.0, 1.0, 1.0, 1.0, 1.5, 1.5, 1.5, 1.5) + assert pieces[2].knots == (1.5, 1.5, 1.5, 1.5, 2.0, 2.0, 2.0, 2.0) + assert [len(piece.control_points) for piece in pieces] == [5, 4, 4] + # the C2 interior knot of the first piece survives the split + assert pieces[0].continuity()["passed"] is True + deviation = max( + math.dist(curve.evaluate(u), piece.evaluate(u)) + for piece in pieces + for u in [ + piece.domain[0] + (piece.domain[1] - piece.domain[0]) * (index + 0.5) / 61.0 + for index in range(61) + ] + ) + assert deviation == 0.0 + + +def test_a_polyline_with_duplicated_vertices_splits_into_its_segments() -> None: + """The realistic case: a chine profile stored as one degree-1 polygon.""" + curve = NURBSCurve( + degree=1, + knots=(0.0, 0.0, 1.0, 1.0, 2.0, 2.0, 3.0, 3.0), + control_points=( + (0.0, 0.0, 0.0), + (2.0, 0.0, 0.0), + (2.0, 0.0, 0.0), + (2.0, 2.0, 0.0), + (2.0, 2.0, 0.0), + (0.0, 2.0, 0.0), + ), + ) + pieces = curve.split_at_discontinuities() + assert [piece.domain for piece in pieces] == [(0.0, 1.0), (1.0, 2.0), (2.0, 3.0)] + assert sum(piece.arc_length() for piece in pieces) == pytest.approx(6.0, abs=1e-14) + elevated = curve.elevate_degree_piecewise(1) + assert [piece.degree for piece in elevated] == [2, 2, 2] + # elevating a straight degree-1 segment puts its middle control point exactly on + # the analytic midpoint of that segment + expected_first = ((0.0, 0.0, 0.0), (1.0, 0.0, 0.0), (2.0, 0.0, 0.0)) + expected_last = ((2.0, 2.0, 0.0), (1.0, 2.0, 0.0), (0.0, 2.0, 0.0)) + for raised, expected in ((elevated[0], expected_first), (elevated[-1], expected_last)): + for point, target in zip(raised.control_points, expected, strict=True): + assert math.dist(point, target) < 1e-15 + + +def test_a_smooth_curve_is_returned_as_itself() -> None: + circle = circle_curve(RADIUS, arcs=3) + pieces = circle.split_at_discontinuities() + assert len(pieces) == 1 + assert pieces[0] is circle # identity, not a copy + + +def test_the_broken_curve_cannot_be_differentiated_across_the_break() -> None: + """Why the split exists: the whole curve's derivative net has a zero-length interval.""" + with pytest.raises(NURBSError, match="zero-length knot interval"): + composite_curve().arc_length() + + +def test_elevation_of_a_broken_curve_still_refuses_and_names_the_remedy() -> None: + with pytest.raises(NURBSError, match="discontinuous") as error: + composite_curve().elevate_degree(1) + assert "elevate_degree_piecewise" in str(error.value) + + +# --------------------------------------------------------------------------- +# analytic anchors on the pieces (before and after elevation) +# --------------------------------------------------------------------------- +def test_the_arc_piece_is_the_analytic_circle() -> None: + arc = composite_curve().split_at_discontinuities()[0] + radii = [math.hypot(*arc.evaluate(u)[:2]) for u in (index / 200.0 for index in range(201))] + assert max(abs(radius - RADIUS) for radius in radii) < 1e-14 + assert arc.evaluate(0.0) == pytest.approx((RADIUS, 0.0, 0.0), abs=1e-15) + assert arc.evaluate(1.0) == pytest.approx((0.0, RADIUS, 0.0), abs=1e-15) + # the 45-degree point of the rational quadratic is (sqrt 2, sqrt 2) + assert arc.evaluate(0.5) == pytest.approx((math.sqrt(2.0), math.sqrt(2.0), 0.0), abs=1e-15) + + +def test_the_straight_piece_is_the_analytic_line() -> None: + """Independent reference: the Bernstein form of this piece is ``u -> (-2u, 2 - 2u)``.""" + line = composite_curve().split_at_discontinuities()[1] + assert line.domain == (1.0, 2.0) # parameters stay on the original ones + worst = 0.0 + for index in range(101): + u = index / 100.0 + expected = (-2.0 * u, 2.0 - 2.0 * u, 0.0) # local parameter u = knot - 1 + worst = max(worst, math.dist(line.evaluate(1.0 + u), expected)) + assert worst < 1e-15 + + +def test_piecewise_arc_length_matches_the_analytic_length() -> None: + pieces = composite_curve().split_at_discontinuities() + lengths = [piece.arc_length() for piece in pieces] + assert lengths[0] == pytest.approx(ARC_LENGTH, abs=1e-13) + assert lengths[1] == pytest.approx(LINE_LENGTH, abs=1e-13) + assert sum(lengths) == pytest.approx(TOTAL_LENGTH, abs=1e-13) + + +def test_the_kink_angle_is_45_degrees_and_survives_elevation() -> None: + pieces = composite_curve().split_at_discontinuities() + elevated = composite_curve().elevate_degree_piecewise(3) + for left, right in ( + (pieces[0].tangent(1.0), pieces[1].tangent(1.0)), + (elevated[0].tangent(1.0), elevated[1].tangent(1.0)), + ): + cosine = max(-1.0, min(1.0, sum(a * b for a, b in zip(left, right, strict=True)))) + assert math.degrees(math.acos(cosine)) == pytest.approx(45.0, abs=1e-12) + + +# --------------------------------------------------------------------------- +# piecewise elevation +# --------------------------------------------------------------------------- +def test_piecewise_elevation_reproduces_the_geometry_exactly() -> None: + curve = composite_curve() + pieces = curve.split_at_discontinuities() + elevated = curve.elevate_degree_piecewise(2) + assert [piece.degree for piece in elevated] == [4, 4] + deviation = 0.0 + for original, raised in zip(pieces, elevated, strict=True): + low, high = original.domain + for index in range(201): + u = low + (high - low) * index / 200.0 + deviation = max(deviation, math.dist(original.evaluate(u), raised.evaluate(u))) + assert deviation < 1e-14 + # the raised arc is still exactly the rational circle: elevation is exact for + # rational geometry because the solve happens in homogeneous coordinates + raised_arc = elevated[0] + radii = [abs(math.hypot(*raised_arc.evaluate(index / 200.0)[:2]) - RADIUS) for index in range(201)] + assert max(radii) < 1e-13 + + +def test_piecewise_elevation_keeps_the_arc_length() -> None: + elevated = composite_curve().elevate_degree_piecewise(3) + assert sum(piece.arc_length() for piece in elevated) == pytest.approx(TOTAL_LENGTH, abs=1e-12) + + +def test_growth_is_the_minimum_times_spans_per_piece() -> None: + """Each piece grows by exactly ``times * (#spans of the piece)`` control points.""" + curve = jump_curve() + pieces = curve.split_at_discontinuities() + elevated = curve.elevate_degree_piecewise(1) + for original, raised in zip(pieces, elevated, strict=True): + low, high = original.domain + spans = len({value for value in original.knots if low <= value <= high}) - 1 + assert len(raised.control_points) - len(original.control_points) == spans + assert raised.degree == original.degree + 1 + + +def test_piecewise_elevation_does_not_close_the_breaks() -> None: + """A break stays a break: the two sides of each jump stay exactly where they were.""" + curve = jump_curve() + pieces = curve.split_at_discontinuities() + elevated = curve.elevate_degree_piecewise(2) + for index in (0, 1): + left, right = pieces[index], pieces[index + 1] + raised_left, raised_right = elevated[index], elevated[index + 1] + before = math.dist(left.evaluate(left.domain[1]), right.evaluate(right.domain[0])) + after = math.dist( + raised_left.evaluate(raised_left.domain[1]), + raised_right.evaluate(raised_right.domain[0]), + ) + assert after == pytest.approx(before, abs=1e-14) + assert before > 1.0 # it really is a jump, not a shared point + + +def test_piecewise_elevation_of_a_smooth_curve_equals_direct_elevation() -> None: + circle = circle_curve(RADIUS, arcs=3) + direct = circle.elevate_degree(2) + piecewise = circle.elevate_degree_piecewise(2) + assert len(piecewise) == 1 + assert piecewise[0].degree == direct.degree + deviation = max( + math.dist(direct.evaluate(index / 300.0), piecewise[0].evaluate(index / 300.0)) + for index in range(301) + ) + assert deviation < 1e-13 + + +def test_zero_times_is_the_identity_and_negative_times_is_rejected() -> None: + curve = composite_curve() + pieces = curve.split_at_discontinuities() + unchanged = curve.elevate_degree_piecewise(0) + assert [piece.degree for piece in unchanged] == [2, 2] + assert [piece.knots for piece in unchanged] == [piece.knots for piece in pieces] + for bad in (-1, 1.5, True): + with pytest.raises(NURBSError, match="non-negative integer"): + curve.elevate_degree_piecewise(bad)