From 02b4fa1ca7497c3c3cc80f6dbc169820eb481fbd Mon Sep 17 00:00:00 2001 From: Mark Tolmacs Date: Mon, 19 May 2025 16:55:51 +0200 Subject: [PATCH] Fixing mid point calculation precision --- packages/element/src/linearElementEditor.ts | 4 +- .../tests/linearElementEditor.test.tsx | 36 +- packages/excalidraw/tests/helpers/ui.ts | 31 +- packages/math/src/constants.ts | 55 +++ packages/math/src/curve.ts | 323 +++++++++++------- packages/math/src/index.ts | 1 + packages/math/tests/curve.test.ts | 44 +++ 7 files changed, 332 insertions(+), 162 deletions(-) create mode 100644 packages/math/src/constants.ts diff --git a/packages/element/src/linearElementEditor.ts b/packages/element/src/linearElementEditor.ts index 0a8500d09..4ba051807 100644 --- a/packages/element/src/linearElementEditor.ts +++ b/packages/element/src/linearElementEditor.ts @@ -10,7 +10,7 @@ import { isCurve, isLineSegment, curveLength, - curveHalfPoint, + curvePointAtLength, } from "@excalidraw/math"; import { getCurvePathOps } from "@excalidraw/utils/shape"; @@ -707,7 +707,7 @@ export class LinearElementEditor { switch (true) { case isCurve(shape): - return curveHalfPoint(shape as Curve); + return curvePointAtLength(shape as Curve, 0.5); case isLineSegment(shape): return pointCenter(shape[0] as GlobalPoint, shape[1] as GlobalPoint); } diff --git a/packages/element/tests/linearElementEditor.test.tsx b/packages/element/tests/linearElementEditor.test.tsx index f4ba1fbe8..5a5114176 100644 --- a/packages/element/tests/linearElementEditor.test.tsx +++ b/packages/element/tests/linearElementEditor.test.tsx @@ -346,12 +346,12 @@ describe("Test Linear Elements", () => { expect(midPointsWithRoundEdge).toMatchInlineSnapshot(` [ [ - "56.87500", - "48.12500", + "54.27552", + "46.16120", ], [ - "79.37500", - "48.12500", + "76.95494", + "44.56052", ], ] `); @@ -411,12 +411,12 @@ describe("Test Linear Elements", () => { expect(newMidPoints).toMatchInlineSnapshot(` [ [ - "106.87500", - "68.12500", + "104.27552", + "66.16120", ], [ - "129.37500", - "68.12500", + "126.95494", + "64.56052", ], ] `); @@ -727,12 +727,12 @@ describe("Test Linear Elements", () => { expect(newMidPoints).toMatchInlineSnapshot(` [ [ - "31.87500", - "23.12500", + "29.28349", + "20.91105", ], [ - "82.50000", - "51.25000", + "78.86048", + "46.12277", ], ] `); @@ -816,12 +816,12 @@ describe("Test Linear Elements", () => { expect(newMidPoints).toMatchInlineSnapshot(` [ [ - "56.87500", - "48.12500", + "54.27552", + "46.16120", ], [ - "79.37500", - "48.12500", + "76.95494", + "44.56052", ], ] `); @@ -983,8 +983,8 @@ describe("Test Linear Elements", () => { ); expect(position).toMatchInlineSnapshot(` { - "x": "90.23300", - "y": "81.68631", + "x": "86.17305", + "y": "76.11251", } `); }); diff --git a/packages/excalidraw/tests/helpers/ui.ts b/packages/excalidraw/tests/helpers/ui.ts index abfadf331..c22d8b467 100644 --- a/packages/excalidraw/tests/helpers/ui.ts +++ b/packages/excalidraw/tests/helpers/ui.ts @@ -383,25 +383,22 @@ const proxy = ( the proxy */ get(): typeof element; } => { - return new Proxy( - {}, - { - get(target, prop) { - const currentElement = h.elements.find( - ({ id }) => id === element.id, - ) as any; - if (prop === "get") { - if (currentElement.hasOwnProperty("get")) { - throw new Error( - "trying to get `get` test property, but ExcalidrawElement seems to define its own", - ); - } - return () => currentElement; + return new Proxy(element, { + get(target, prop) { + const currentElement = h.elements.find( + ({ id }) => id === element.id, + ) as any; + if (prop === "get") { + if (currentElement.hasOwnProperty("get")) { + throw new Error( + "trying to get `get` test property, but ExcalidrawElement seems to define its own", + ); } - return currentElement[prop]; - }, + return () => currentElement; + } + return currentElement[prop]; }, - ) as any; + }) as any; }; /** Tools that can be used to draw shapes */ diff --git a/packages/math/src/constants.ts b/packages/math/src/constants.ts new file mode 100644 index 000000000..6ba5eb8f7 --- /dev/null +++ b/packages/math/src/constants.ts @@ -0,0 +1,55 @@ +// Legendre-Gauss abscissae (x values) and weights for n=24 +// Refeerence: https://pomax.github.io/bezierinfo/legendre-gauss.html +export const LegendreGaussN24TValues = [ + -0.0640568928626056260850430826247450385909, + 0.0640568928626056260850430826247450385909, + -0.1911188674736163091586398207570696318404, + 0.1911188674736163091586398207570696318404, + -0.3150426796961633743867932913198102407864, + 0.3150426796961633743867932913198102407864, + -0.4337935076260451384870842319133497124524, + 0.4337935076260451384870842319133497124524, + -0.5454214713888395356583756172183723700107, + 0.5454214713888395356583756172183723700107, + -0.6480936519369755692524957869107476266696, + 0.6480936519369755692524957869107476266696, + -0.7401241915785543642438281030999784255232, + 0.7401241915785543642438281030999784255232, + -0.8200019859739029219539498726697452080761, + 0.8200019859739029219539498726697452080761, + -0.8864155270044010342131543419821967550873, + 0.8864155270044010342131543419821967550873, + -0.9382745520027327585236490017087214496548, + 0.9382745520027327585236490017087214496548, + -0.9747285559713094981983919930081690617411, + 0.9747285559713094981983919930081690617411, + -0.9951872199970213601799974097007368118745, + 0.9951872199970213601799974097007368118745, +]; + +export const LegendreGaussN24CValues = [ + 0.1279381953467521569740561652246953718517, + 0.1279381953467521569740561652246953718517, + 0.1258374563468282961213753825111836887264, + 0.1258374563468282961213753825111836887264, + 0.121670472927803391204463153476262425607, + 0.121670472927803391204463153476262425607, + 0.1155056680537256013533444839067835598622, + 0.1155056680537256013533444839067835598622, + 0.1074442701159656347825773424466062227946, + 0.1074442701159656347825773424466062227946, + 0.0976186521041138882698806644642471544279, + 0.0976186521041138882698806644642471544279, + 0.086190161531953275917185202983742667185, + 0.086190161531953275917185202983742667185, + 0.0733464814110803057340336152531165181193, + 0.0733464814110803057340336152531165181193, + 0.0592985849154367807463677585001085845412, + 0.0592985849154367807463677585001085845412, + 0.0442774388174198061686027482113382288593, + 0.0442774388174198061686027482113382288593, + 0.0285313886289336631813078159518782864491, + 0.0285313886289336631813078159518782864491, + 0.0123412297999871995468056670700372915759, + 0.0123412297999871995468056670700372915759, +]; diff --git a/packages/math/src/curve.ts b/packages/math/src/curve.ts index 4909276e7..7664c70f6 100644 --- a/packages/math/src/curve.ts +++ b/packages/math/src/curve.ts @@ -1,16 +1,7 @@ -import { invariant } from "@excalidraw/common"; - -import type { Bounds } from "@excalidraw/element"; - -import { - isPoint, - pointCenter, - pointDistance, - pointFrom, - pointFromVector, -} from "./point"; +import { isPoint, pointDistance, pointFrom, pointFromVector } from "./point"; import { rectangle, rectangleIntersectLineSegment } from "./rectangle"; import { vector, vectorNormal, vectorNormalize, vectorScale } from "./vector"; +import { LegendreGaussN24CValues, LegendreGaussN24TValues } from "./constants"; import type { Curve, GlobalPoint, LineSegment, LocalPoint } from "./types"; @@ -31,81 +22,6 @@ export function curve( return [a, b, c, d] as Curve; } -function gradient( - f: (t: number, s: number) => number, - t0: number, - s0: number, - delta: number = 1e-6, -): number[] { - return [ - (f(t0 + delta, s0) - f(t0 - delta, s0)) / (2 * delta), - (f(t0, s0 + delta) - f(t0, s0 - delta)) / (2 * delta), - ]; -} - -function solve( - f: (t: number, s: number) => [number, number], - t0: number, - s0: number, - tolerance: number = 1e-3, - iterLimit: number = 10, -): number[] | null { - let error = Infinity; - let iter = 0; - - while (error >= tolerance) { - if (iter >= iterLimit) { - return null; - } - - const y0 = f(t0, s0); - const jacobian = [ - gradient((t, s) => f(t, s)[0], t0, s0), - gradient((t, s) => f(t, s)[1], t0, s0), - ]; - const b = [[-y0[0]], [-y0[1]]]; - const det = - jacobian[0][0] * jacobian[1][1] - jacobian[0][1] * jacobian[1][0]; - - if (det === 0) { - return null; - } - - const iJ = [ - [jacobian[1][1] / det, -jacobian[0][1] / det], - [-jacobian[1][0] / det, jacobian[0][0] / det], - ]; - const h = [ - [iJ[0][0] * b[0][0] + iJ[0][1] * b[1][0]], - [iJ[1][0] * b[0][0] + iJ[1][1] * b[1][0]], - ]; - - t0 = t0 + h[0][0]; - s0 = s0 + h[1][0]; - - const [tErr, sErr] = f(t0, s0); - error = Math.max(Math.abs(tErr), Math.abs(sErr)); - iter += 1; - } - - return [t0, s0]; -} - -export const bezierEquation = ( - c: Curve, - t: number, -) => - pointFrom( - (1 - t) ** 3 * c[0][0] + - 3 * (1 - t) ** 2 * t * c[1][0] + - 3 * (1 - t) * t ** 2 * c[2][0] + - t ** 3 * c[3][0], - (1 - t) ** 3 * c[0][1] + - 3 * (1 - t) ** 2 * t * c[1][1] + - 3 * (1 - t) * t ** 2 * c[2][1] + - t ** 3 * c[3][1], - ); - /** * Computes the intersection between a cubic spline and a line segment. */ @@ -113,12 +29,19 @@ export function curveIntersectLineSegment< Point extends GlobalPoint | LocalPoint, >(c: Curve, l: LineSegment): Point[] { // Optimize by doing a cheap bounding box check first - const bounds = curveBounds(c); + const [p0, p1, p2, p3] = c; + if ( rectangleIntersectLineSegment( rectangle( - pointFrom(bounds[0], bounds[1]), - pointFrom(bounds[2], bounds[3]), + pointFrom( + Math.min(p0[0], p1[0], p2[0], p3[0]), + Math.min(p0[1], p1[1], p2[1], p3[1]), + ), + pointFrom( + Math.max(p0[0], p1[0], p2[0], p3[0]), + Math.max(p0[1], p1[1], p2[1], p3[1]), + ), ), l, ).length === 0 @@ -303,15 +226,6 @@ export function curveTangent( ); } -function curveBounds( - c: Curve, -): Bounds { - const [P0, P1, P2, P3] = c; - const x = [P0[0], P1[0], P2[0], P3[0]]; - const y = [P0[1], P1[1], P2[1], P3[1]]; - return [Math.min(...x), Math.min(...y), Math.max(...x), Math.max(...y)]; -} - export function curveCatmullRomQuadraticApproxPoints( points: GlobalPoint[], tension = 0.5, @@ -417,42 +331,201 @@ export function offsetPointsForQuadraticBezier( return offsetPoints; } +/** + * Implementation based on Legendre-Gauss quadrature for more accurate arc + * length calculation. + * + * Reference: https://pomax.github.io/bezierinfo/#arclength + * + * @param c The curve to calculate the length of + * @returns The approximated length of the curve + */ export function curveLength

( c: Curve

, ): number { - // Use numerical integration to approximate the curve length - const steps = 50; - let length = 0; + const z2 = 0.5; + let sum = 0; - // Calculate length by summing segments - let prevPoint = bezierEquation(c, 0); - for (let i = 1; i <= steps; i++) { - const t = i / steps; - const currentPoint = bezierEquation(c, t); - length += pointDistance(prevPoint, currentPoint); - prevPoint = currentPoint; + for (let i = 0; i < 24; i++) { + const t = z2 * LegendreGaussN24TValues[i] + z2; + const derivativeVector = curveTangent(c, t); + const magnitude = Math.sqrt( + derivativeVector[0] * derivativeVector[0] + + derivativeVector[1] * derivativeVector[1], + ); + sum += LegendreGaussN24CValues[i] * magnitude; } - return length; + return z2 * sum; } -export function curveHalfPoint

( +/** + * Calculates the curve length from t=0 to t=parameter using the same + * Legendre-Gauss quadrature method used in curveLength + * + * @param c The curve to calculate the partial length for + * @param t The parameter value (0 to 1) to calculate length up to + * @returns The length of the curve from beginning to parameter t + */ +export function curveLengthAtParameter

( c: Curve

, -): P { - const steps = 50; - const halfLength = curveLength(c) / 2; - - let length = 0; - let prevPoint = bezierEquation(c, 0); - for (let i = 1; i <= steps; i++) { - const t = i / steps; - const currentPoint = bezierEquation(c, t); - length += pointDistance(prevPoint, currentPoint); - if (length > halfLength) { - return pointCenter(prevPoint, currentPoint); - } - prevPoint = currentPoint; + t: number, +): number { + if (t <= 0) { + return 0; + } + if (t >= 1) { + return curveLength(c); } - invariant(false, "No half point found (this should not happen)"); + // Scale and shift the integration interval from [0,t] to [-1,1] + // which is what the Legendre-Gauss quadrature expects + const z1 = t / 2; + const z2 = t / 2; + + let sum = 0; + + for (let i = 0; i < 24; i++) { + const parameter = z1 * LegendreGaussN24TValues[i] + z2; + const derivativeVector = curveTangent(c, parameter); + const magnitude = Math.sqrt( + derivativeVector[0] * derivativeVector[0] + + derivativeVector[1] * derivativeVector[1], + ); + sum += LegendreGaussN24CValues[i] * magnitude; + } + + return z1 * sum; // Scale the result back to the original interval +} + +/** + * Calculates the point at a specific percentage of a curve's total length + * using binary search for improved efficiency and accuracy. + * + * @param c The curve to calculate point on + * @param percent A value between 0 and 1 representing the percentage of the curve's length + * @returns The point at the specified percentage of curve length + */ +export function curvePointAtLength

( + c: Curve

, + percent: number, +): P { + if (percent <= 0) { + return bezierEquation(c, 0); + } + + if (percent >= 1) { + return bezierEquation(c, 1); + } + + const totalLength = curveLength(c); + const targetLength = totalLength * percent; + + // Binary search to find parameter t where length at t equals target length + let tMin = 0; + let tMax = 1; + let t = percent; // Start with a reasonable guess (t = percent) + let currentLength = 0; + + // Tolerance for length comparison and iteration limit to avoid infinite loops + const tolerance = totalLength * 0.0001; // 0.01% of total length + const maxIterations = 20; + + for (let iteration = 0; iteration < maxIterations; iteration++) { + currentLength = curveLengthAtParameter(c, t); + const error = Math.abs(currentLength - targetLength); + + // If we're close enough, return the point + if (error < tolerance) { + break; + } + + // Update search range + if (currentLength < targetLength) { + tMin = t; + } else { + tMax = t; + } + + // Update parameter using bisection method + t = (tMin + tMax) / 2; + } + + return bezierEquation(c, t); +} + +function bezierEquation( + c: Curve, + t: number, +) { + return pointFrom( + (1 - t) ** 3 * c[0][0] + + 3 * (1 - t) ** 2 * t * c[1][0] + + 3 * (1 - t) * t ** 2 * c[2][0] + + t ** 3 * c[3][0], + (1 - t) ** 3 * c[0][1] + + 3 * (1 - t) ** 2 * t * c[1][1] + + 3 * (1 - t) * t ** 2 * c[2][1] + + t ** 3 * c[3][1], + ); +} + +function gradient( + f: (t: number, s: number) => number, + t0: number, + s0: number, + delta: number = 1e-6, +): number[] { + return [ + (f(t0 + delta, s0) - f(t0 - delta, s0)) / (2 * delta), + (f(t0, s0 + delta) - f(t0, s0 - delta)) / (2 * delta), + ]; +} + +function solve( + f: (t: number, s: number) => [number, number], + t0: number, + s0: number, + tolerance: number = 1e-3, + iterLimit: number = 10, +): number[] | null { + let error = Infinity; + let iter = 0; + + while (error >= tolerance) { + if (iter >= iterLimit) { + return null; + } + + const y0 = f(t0, s0); + const jacobian = [ + gradient((t, s) => f(t, s)[0], t0, s0), + gradient((t, s) => f(t, s)[1], t0, s0), + ]; + const b = [[-y0[0]], [-y0[1]]]; + const det = + jacobian[0][0] * jacobian[1][1] - jacobian[0][1] * jacobian[1][0]; + + if (det === 0) { + return null; + } + + const iJ = [ + [jacobian[1][1] / det, -jacobian[0][1] / det], + [-jacobian[1][0] / det, jacobian[0][0] / det], + ]; + const h = [ + [iJ[0][0] * b[0][0] + iJ[0][1] * b[1][0]], + [iJ[1][0] * b[0][0] + iJ[1][1] * b[1][0]], + ]; + + t0 = t0 + h[0][0]; + s0 = s0 + h[1][0]; + + const [tErr, sErr] = f(t0, s0); + error = Math.max(Math.abs(tErr), Math.abs(sErr)); + iter += 1; + } + + return [t0, s0]; } diff --git a/packages/math/src/index.ts b/packages/math/src/index.ts index d00ab469d..d1ebf536c 100644 --- a/packages/math/src/index.ts +++ b/packages/math/src/index.ts @@ -1,4 +1,5 @@ export * from "./angle"; +export * from "./constants"; export * from "./curve"; export * from "./line"; export * from "./point"; diff --git a/packages/math/tests/curve.test.ts b/packages/math/tests/curve.test.ts index 739562096..ccdb5194b 100644 --- a/packages/math/tests/curve.test.ts +++ b/packages/math/tests/curve.test.ts @@ -4,6 +4,9 @@ import { curve, curveClosestPoint, curveIntersectLineSegment, + curveLength, + curveLengthAtParameter, + curvePointAtLength, curvePointDistance, } from "../src/curve"; import { pointFrom } from "../src/point"; @@ -99,4 +102,45 @@ describe("Math curve", () => { expect(curvePointDistance(c, p)).toBeCloseTo(6.695873043213627); }); }); + + describe("length", () => { + it("can be determined", () => { + const c = curve( + pointFrom(-50, -50), + pointFrom(10, -50), + pointFrom(10, 50), + pointFrom(50, 50), + ); + + expect(curveLength(c)).toBeCloseTo(150.0, 0); + }); + }); + + describe("point at given parameter", () => { + it("can be determined", () => { + const c = curve( + pointFrom(-50, -50), + pointFrom(10, -50), + pointFrom(10, 50), + pointFrom(50, 50), + ); + + expect(curveLengthAtParameter(c, 0.5)).toBeCloseTo(80.83); + }); + }); + + describe("point at given length", () => { + it("can be determined", () => { + const c = curve( + pointFrom(-50, -50), + pointFrom(10, -50), + pointFrom(10, 50), + pointFrom(50, 50), + ); + + expect(curvePointAtLength(c, 0.5)).toEqual([ + 4.802938740176614, -5.301185927237384, + ]); + }); + }); });