From 9d36d1a09da43c53b6452cda27f74f5b4a0080ed Mon Sep 17 00:00:00 2001 From: "Saeed.Torabi" <114527638+saeeedtorabi@users.noreply.github.com> Date: Fri, 14 Nov 2025 15:58:45 -0500 Subject: [PATCH 01/15] Improved CurveCurveIntersectXY --- .../internalContexts/CurveCurveIntersectXY.ts | 30 ++++++++--- .../test/curve/CurveCurveIntersectXY.test.ts | 52 +++++++++++++------ 2 files changed, 61 insertions(+), 21 deletions(-) diff --git a/core/geometry/src/curve/internalContexts/CurveCurveIntersectXY.ts b/core/geometry/src/curve/internalContexts/CurveCurveIntersectXY.ts index e6c7007e8bc..842c7514d3c 100644 --- a/core/geometry/src/curve/internalContexts/CurveCurveIntersectXY.ts +++ b/core/geometry/src/curve/internalContexts/CurveCurveIntersectXY.ts @@ -40,6 +40,7 @@ import { Loop } from "../Loop"; import { Path } from "../Path"; import { ProxyCurve } from "../ProxyCurve"; import { TransitionSpiral3d } from "../spiral/TransitionSpiral3d"; +import { CurveCurve } from "../CurveCurve"; // cspell:word XYRR /** @@ -1057,7 +1058,8 @@ export class CurveCurveIntersectXY extends RecurseToCurvesGeometryHandler { const spiralApproximation = LineString3d.create(); this._geometryB.emitStrokes(spiralApproximation); const numPreviousResults = this._results.length; - this.computeSegmentLineString(segmentA, this._extendA0, this._extendA1, spiralApproximation, false, false, false); + // this.computeSegmentLineString(segmentA, this._extendA0, this._extendA1, spiralApproximation, false, false, false); + this._results.push(...CurveCurve.closeApproachProjectedXYPairs(segmentA, spiralApproximation, Geometry.smallMetricDistance)); const numberOfNewResults = this._results.length - numPreviousResults; this.refineSpiralResultsByNewton(segmentA, this._geometryB, numberOfNewResults); } else if (this._geometryB instanceof CurveCollection) { @@ -1088,7 +1090,8 @@ export class CurveCurveIntersectXY extends RecurseToCurvesGeometryHandler { const spiralApproximation = LineString3d.create(); this._geometryB.emitStrokes(spiralApproximation); const numPreviousResults = this._results.length; - this.computeLineStringLineString(lsA, this._extendA0, this._extendA1, spiralApproximation, false, false, false); + // this.computeLineStringLineString(lsA, this._extendA0, this._extendA1, spiralApproximation, false, false, false); + this._results.push(...CurveCurve.closeApproachProjectedXYPairs(lsA, spiralApproximation, Geometry.smallMetricDistance)); const numberOfNewResults = this._results.length - numPreviousResults; this.refineSpiralResultsByNewton(lsA, this._geometryB, numberOfNewResults); } else if (this._geometryB instanceof CurveCollection) { @@ -1119,7 +1122,8 @@ export class CurveCurveIntersectXY extends RecurseToCurvesGeometryHandler { const spiralApproximation = LineString3d.create(); this._geometryB.emitStrokes(spiralApproximation); const numPreviousResults = this._results.length; - this.computeArcLineString(arc0, this._extendA0, this._extendA1, spiralApproximation, false, false, false); + // this.computeArcLineString(arc0, this._extendA0, this._extendA1, spiralApproximation, false, false, false); + this._results.push(...CurveCurve.closeApproachProjectedXYPairs(arc0, spiralApproximation, Geometry.smallMetricDistance)); const numberOfNewResults = this._results.length - numPreviousResults; this.refineSpiralResultsByNewton(arc0, this._geometryB, numberOfNewResults); } else if (this._geometryB instanceof CurveCollection) { @@ -1150,7 +1154,8 @@ export class CurveCurveIntersectXY extends RecurseToCurvesGeometryHandler { const spiralApproximation = LineString3d.create(); this._geometryB.emitStrokes(spiralApproximation); const numPreviousResults = this._results.length; - this.dispatchLineStringBSplineCurve(spiralApproximation, false, false, curve, this._extendA0, this._extendA1, true); + // this.dispatchLineStringBSplineCurve(spiralApproximation, false, false, curve, this._extendA0, this._extendA1, true); + this._results.push(...CurveCurve.closeApproachProjectedXYPairs(curve, spiralApproximation, Geometry.smallMetricDistance)); const numberOfNewResults = this._results.length - numPreviousResults; this.refineSpiralResultsByNewton(curve, this._geometryB, numberOfNewResults); } else if (this._geometryB instanceof CurveCollection) { @@ -1177,6 +1182,10 @@ export class CurveCurveIntersectXY extends RecurseToCurvesGeometryHandler { for (const detail of resultsToBeRefined) { let spiralFraction = reversed ? detail.detailA.fraction : detail.detailB.fraction; let otherFraction = reversed ? detail.detailB.fraction : detail.detailA.fraction; + if (spiral.fractionToPoint(spiralFraction).isAlmostEqualXY(otherCurve.fractionToPoint(otherFraction))) { // already an accurate intersection + this.recordPointWithLocalFractions(otherFraction, otherCurve, 0, 1, spiralFraction, spiral, 0, 1, reversed); + continue; + } const xyMatchingFunction = new CurveCurveIntersectionXYRRToRRD(spiral, otherCurve); const newtonSearcher = new Newton2dUnboundedWithDerivative(xyMatchingFunction); newtonSearcher.setUV(spiralFraction, otherFraction); @@ -1190,11 +1199,20 @@ export class CurveCurveIntersectXY extends RecurseToCurvesGeometryHandler { } /** Double dispatch handler for strongly typed spiral curve. */ public override handleTransitionSpiral(spiral: TransitionSpiral3d): any { - if (this._geometryB instanceof CurvePrimitive) { // this also handles CurveChainWithDistanceIndex + if (this._geometryB instanceof CurveChainWithDistanceIndex) { + this.dispatchCurveChainWithDistanceIndex(spiral, this.handleTransitionSpiral.bind(this)); + } else if (this._geometryB instanceof CurvePrimitive) { const spiralApproximation = LineString3d.create(); spiral.emitStrokes(spiralApproximation); + const geomBApproximation = LineString3d.create(); + if (this._geometryB instanceof TransitionSpiral3d) + this._geometryB.emitStrokes(geomBApproximation); const numPreviousResults = this._results.length; - this.handleLineString3d(spiralApproximation); + // this.handleLineString3d(spiralApproximation); + if (this._geometryB instanceof TransitionSpiral3d) + this._results.push(...CurveCurve.closeApproachProjectedXYPairs(spiralApproximation, geomBApproximation, Geometry.smallMetricDistance)); + else + this._results.push(...CurveCurve.closeApproachProjectedXYPairs(spiralApproximation, this._geometryB, Geometry.smallMetricDistance)); const numberOfNewResults = this._results.length - numPreviousResults; this.refineSpiralResultsByNewton(this._geometryB, spiral, numberOfNewResults, true); } else if (this._geometryB instanceof CurveCollection) { diff --git a/core/geometry/src/test/curve/CurveCurveIntersectXY.test.ts b/core/geometry/src/test/curve/CurveCurveIntersectXY.test.ts index d671d45bf93..c346b902a1f 100644 --- a/core/geometry/src/test/curve/CurveCurveIntersectXY.test.ts +++ b/core/geometry/src/test/curve/CurveCurveIntersectXY.test.ts @@ -2503,11 +2503,13 @@ describe("CurveCurveIntersectXY", () => { } // curve primitives const lineSegment0 = LineSegment3d.create(Point3d.create(70, 30), Point3d.create(70, -30)); + const lineSegment1 = LineSegment3d.create(Point3d.create(-20, 0), Point3d.create(100, 0)); const lineString0 = LineString3d.create( Point3d.create(20, 20), Point3d.create(40, 20), Point3d.create(70, -40), Point3d.create(90, -40), Point3d.create(70, 20), Point3d.create(130, 20), ); const arc0 = Arc3d.createXY(Point3d.create(50, 10), 30); + const arc1 = Arc3d.createXY(Point3d.create(0, -30), 30); const bspline0 = BSplineCurve3d.createUniformKnots( [ Point3d.create(70, 50, 0), @@ -2520,28 +2522,30 @@ describe("CurveCurveIntersectXY", () => { )!; // curve collection (path-loop), curve chain, and bag of curves const lineString1 = LineString3d.create(Point3d.create(-36.33, 64.45), Point3d.create(8.34, 120.78), Point3d.create(76.58, -7)); - const arc1 = Arc3d.create( + const arc2 = Arc3d.create( Point3d.create(80, 40), Vector3d.create(10, 0), Vector3d.create(0, 50), AngleSweep.createStartEndDegrees(250, 90), ); - const arc2 = Arc3d.create( + const arc3 = Arc3d.create( Point3d.create(70, -40), Vector3d.create(20, 0), Vector3d.create(0, 20), AngleSweep.createStartEndDegrees(0, -180), ); const lineString2 = LineString3d.create(Point3d.create(80, 90), Point3d.create(90, -100), Point3d.create(40, -100)); const lineString3 = LineString3d.create(Point3d.create(50, -40), Point3d.create(0, -40), Point3d.create(0, 0)); - const lineSegment1 = LineSegment3d.create(Point3d.create(40, -100), Point3d.create(40, 30)); - const lineSegment2 = LineSegment3d.create(Point3d.create(40, 30), Point3d.create(70, -70)); - const lineSegment3 = LineSegment3d.create(Point3d.create(60, -50), Point3d.create(90, -40)); - const path0 = Path.create(arc1, lineString2, lineSegment1, lineSegment2, integratedSpirals[1]); - const path1 = Path.create(lineSegment3, arc2, lineString3, integratedSpirals[0]); - const loop = Path.create(lineString1, arc1, lineString2, lineSegment1, lineSegment2, integratedSpirals[1]); + const lineSegment2 = LineSegment3d.create(Point3d.create(40, -100), Point3d.create(40, 30)); + const lineSegment3 = LineSegment3d.create(Point3d.create(40, 30), Point3d.create(70, -70)); + const lineSegment4 = LineSegment3d.create(Point3d.create(60, -50), Point3d.create(90, -40)); + const path0 = Path.create(arc2, lineString2, lineSegment2, lineSegment3, integratedSpirals[1]); + const path1 = Path.create(lineSegment4, arc3, lineString3, integratedSpirals[0]); + const loop = Path.create(lineString1, arc2, lineString2, lineSegment2, lineSegment3, integratedSpirals[1]); const curveChain0 = CurveChainWithDistanceIndex.createCapture(path0); const curveChain1 = CurveChainWithDistanceIndex.createCapture(path1); const bagOfCurves = BagOfCurves.create(path0, arc0, lineString0); const curves: AnyCurve[] = [ + lineSegment1, lineSegment0, lineString0, arc0, + arc1, bspline0, // add rotated and non-planar spirals integratedSpirals[1], @@ -2581,14 +2585,14 @@ describe("CurveCurveIntersectXY", () => { bagOfCurves, ]; const numExpectedIntersections = [ - 1, 3, 2, 1, // curve primitives other than spirals + 1, 1, 3, 2, 1, 1, // curve primitives other than spirals 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, // rotated and non-planar integrated spirals 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, // 1, 1, 1, 1, 1, 1, 1, 1 // rotated and non-planar direct spirals 5, 6, 5, 10, // path, loop, curve chain, and bag of curves ]; ck.testCoordinate(curves.length, numExpectedIntersections.length, "matching arrays"); // spiral vs all curves - const test = (spiral: TransitionSpiral3d) => { + const test0 = (spiral: TransitionSpiral3d) => { for (let j = 0; j < curves.length; j++) { const curve = curves[j]; const numExpectedIntersection = numExpectedIntersections[j]; @@ -2600,15 +2604,15 @@ describe("CurveCurveIntersectXY", () => { } for (let i = 0; i < integratedSpirals.length; i++) // skip rotated and non-planar integrated spirals if (i % 3 === 0) - test(integratedSpirals[i]); + test0(integratedSpirals[i]); dx += 250; for (let i = 0; i < directSpirals.length; i++) // skip rotated and non-planar direct spirals if (i % 3 === 0) - test(directSpirals[i]); - // curve chain/collection vs curve chain/collection + test0(directSpirals[i]); dx = 0; - dy = 6000; - const numExpected = 12; + dy = 6400; + let numExpected = 12; + // curve chain/collection vs curve chain/collection visualizeAndTestSpiralIntersection(ck, allGeometry, curveChain0, curveChain1, numExpected, dx, dy); dy += 200; visualizeAndTestSpiralIntersection(ck, allGeometry, path0, path1, numExpected, dx, dy); @@ -2616,6 +2620,24 @@ describe("CurveCurveIntersectXY", () => { visualizeAndTestSpiralIntersection(ck, allGeometry, curveChain0, path1, numExpected, dx, dy); dy += 200; visualizeAndTestSpiralIntersection(ck, allGeometry, curveChain1, path0, numExpected, dx, dy); + // tangency at the interior of the spiral + dy += 200; + numExpected = 1; + const test1 = (spiral: TransitionSpiral3d) => { + const ray = spiral.fractionToPointAndDerivative(0.5); + const ls = LineString3d.create( + ray.origin.plusScaled(ray.direction.normalize()!, 50), ray.origin.plusScaled(ray.direction.normalize()!, -50) + ); + visualizeAndTestSpiralIntersection(ck, allGeometry, spiral, ls, numExpected, dx, dy); + dx += 200; + } + for (let i = 0; i < integratedSpirals.length; i++) // skip rotated and non-planar integrated spirals + if (i % 3 === 0) + test1(integratedSpirals[i]); + dx += 250; + for (let i = 0; i < directSpirals.length; i++) // skip rotated and non-planar direct spirals + if (i % 3 === 0) + test1(directSpirals[i]); GeometryCoreTestIO.saveGeometry(allGeometry, "CurveCurveIntersectXY", "SpiralIntersection"); expect(ck.getNumErrors()).toBe(0); From 7f0db82704aa49272fd3e51bb314992d6e68db4c Mon Sep 17 00:00:00 2001 From: "Saeed.Torabi" <114527638+saeeedtorabi@users.noreply.github.com> Date: Mon, 17 Nov 2025 14:58:15 -0500 Subject: [PATCH 02/15] Added support for extend --- .../internalContexts/CurveCurveIntersectXY.ts | 19 ++++++--- .../test/curve/CurveCurveIntersectXY.test.ts | 41 +++++++++++++++---- 2 files changed, 46 insertions(+), 14 deletions(-) diff --git a/core/geometry/src/curve/internalContexts/CurveCurveIntersectXY.ts b/core/geometry/src/curve/internalContexts/CurveCurveIntersectXY.ts index 842c7514d3c..c2b2c2e03e2 100644 --- a/core/geometry/src/curve/internalContexts/CurveCurveIntersectXY.ts +++ b/core/geometry/src/curve/internalContexts/CurveCurveIntersectXY.ts @@ -1058,8 +1058,9 @@ export class CurveCurveIntersectXY extends RecurseToCurvesGeometryHandler { const spiralApproximation = LineString3d.create(); this._geometryB.emitStrokes(spiralApproximation); const numPreviousResults = this._results.length; - // this.computeSegmentLineString(segmentA, this._extendA0, this._extendA1, spiralApproximation, false, false, false); this._results.push(...CurveCurve.closeApproachProjectedXYPairs(segmentA, spiralApproximation, Geometry.smallMetricDistance)); + if (this._extendA0 || this._extendA1) + this.computeSegmentLineString(segmentA, this._extendA0, this._extendA1, spiralApproximation, false, false, false); const numberOfNewResults = this._results.length - numPreviousResults; this.refineSpiralResultsByNewton(segmentA, this._geometryB, numberOfNewResults); } else if (this._geometryB instanceof CurveCollection) { @@ -1090,8 +1091,9 @@ export class CurveCurveIntersectXY extends RecurseToCurvesGeometryHandler { const spiralApproximation = LineString3d.create(); this._geometryB.emitStrokes(spiralApproximation); const numPreviousResults = this._results.length; - // this.computeLineStringLineString(lsA, this._extendA0, this._extendA1, spiralApproximation, false, false, false); this._results.push(...CurveCurve.closeApproachProjectedXYPairs(lsA, spiralApproximation, Geometry.smallMetricDistance)); + if (this._extendA0 || this._extendA1) + this.computeLineStringLineString(lsA, this._extendA0, this._extendA1, spiralApproximation, false, false, false); const numberOfNewResults = this._results.length - numPreviousResults; this.refineSpiralResultsByNewton(lsA, this._geometryB, numberOfNewResults); } else if (this._geometryB instanceof CurveCollection) { @@ -1122,8 +1124,9 @@ export class CurveCurveIntersectXY extends RecurseToCurvesGeometryHandler { const spiralApproximation = LineString3d.create(); this._geometryB.emitStrokes(spiralApproximation); const numPreviousResults = this._results.length; - // this.computeArcLineString(arc0, this._extendA0, this._extendA1, spiralApproximation, false, false, false); this._results.push(...CurveCurve.closeApproachProjectedXYPairs(arc0, spiralApproximation, Geometry.smallMetricDistance)); + if (this._extendA0 || this._extendA1) + this.computeArcLineString(arc0, this._extendA0, this._extendA1, spiralApproximation, false, false, false); const numberOfNewResults = this._results.length - numPreviousResults; this.refineSpiralResultsByNewton(arc0, this._geometryB, numberOfNewResults); } else if (this._geometryB instanceof CurveCollection) { @@ -1154,8 +1157,9 @@ export class CurveCurveIntersectXY extends RecurseToCurvesGeometryHandler { const spiralApproximation = LineString3d.create(); this._geometryB.emitStrokes(spiralApproximation); const numPreviousResults = this._results.length; - // this.dispatchLineStringBSplineCurve(spiralApproximation, false, false, curve, this._extendA0, this._extendA1, true); this._results.push(...CurveCurve.closeApproachProjectedXYPairs(curve, spiralApproximation, Geometry.smallMetricDistance)); + if (this._extendA0 || this._extendA1) + this.dispatchLineStringBSplineCurve(spiralApproximation, false, false, curve, this._extendA0, this._extendA1, true); const numberOfNewResults = this._results.length - numPreviousResults; this.refineSpiralResultsByNewton(curve, this._geometryB, numberOfNewResults); } else if (this._geometryB instanceof CurveCollection) { @@ -1182,6 +1186,8 @@ export class CurveCurveIntersectXY extends RecurseToCurvesGeometryHandler { for (const detail of resultsToBeRefined) { let spiralFraction = reversed ? detail.detailA.fraction : detail.detailB.fraction; let otherFraction = reversed ? detail.detailB.fraction : detail.detailA.fraction; + const otherExtend0 = reversed ? this._extendB0 : this._extendA0; + const otherExtend1 = reversed ? this._extendB1 : this._extendA1; if (spiral.fractionToPoint(spiralFraction).isAlmostEqualXY(otherCurve.fractionToPoint(otherFraction))) { // already an accurate intersection this.recordPointWithLocalFractions(otherFraction, otherCurve, 0, 1, spiralFraction, spiral, 0, 1, reversed); continue; @@ -1192,7 +1198,7 @@ export class CurveCurveIntersectXY extends RecurseToCurvesGeometryHandler { if (newtonSearcher.runIterations()) { spiralFraction = newtonSearcher.getU(); otherFraction = newtonSearcher.getV(); - if (this.acceptFraction(false, spiralFraction, false) && this.acceptFraction(false, otherFraction, false)) + if (this.acceptFraction(false, spiralFraction, false) && this.acceptFraction(otherExtend0, otherFraction, otherExtend1)) this.recordPointWithLocalFractions(otherFraction, otherCurve, 0, 1, spiralFraction, spiral, 0, 1, reversed); } } @@ -1208,11 +1214,12 @@ export class CurveCurveIntersectXY extends RecurseToCurvesGeometryHandler { if (this._geometryB instanceof TransitionSpiral3d) this._geometryB.emitStrokes(geomBApproximation); const numPreviousResults = this._results.length; - // this.handleLineString3d(spiralApproximation); if (this._geometryB instanceof TransitionSpiral3d) this._results.push(...CurveCurve.closeApproachProjectedXYPairs(spiralApproximation, geomBApproximation, Geometry.smallMetricDistance)); else this._results.push(...CurveCurve.closeApproachProjectedXYPairs(spiralApproximation, this._geometryB, Geometry.smallMetricDistance)); + if (this._extendB0 || this._extendB1) + this.handleLineString3d(spiralApproximation); const numberOfNewResults = this._results.length - numPreviousResults; this.refineSpiralResultsByNewton(this._geometryB, spiral, numberOfNewResults, true); } else if (this._geometryB instanceof CurveCollection) { diff --git a/core/geometry/src/test/curve/CurveCurveIntersectXY.test.ts b/core/geometry/src/test/curve/CurveCurveIntersectXY.test.ts index c346b902a1f..d8e7bf11beb 100644 --- a/core/geometry/src/test/curve/CurveCurveIntersectXY.test.ts +++ b/core/geometry/src/test/curve/CurveCurveIntersectXY.test.ts @@ -2373,6 +2373,7 @@ describe("CurveCurveIntersectXY", () => { ck: Checker, allGeometry: GeometryQuery[], curve0: AnyCurve, curve1: AnyCurve, numExpected: number, dx: number, dy: number, + extend0 = false, extend1 = false, ) { GeometryCoreTestIO.captureCloneGeometry(allGeometry, curve0, dx, dy); GeometryCoreTestIO.captureCloneGeometry(allGeometry, curve1, dx, dy); @@ -2427,9 +2428,9 @@ describe("CurveCurveIntersectXY", () => { } } // test both paths - const intersectionsAB = CurveCurve.intersectionXYPairs(curve0, false, curve1, false); + const intersectionsAB = CurveCurve.intersectionXYPairs(curve0, extend0, curve1, extend1); testSpiralIntersection(intersectionsAB, false); - const intersectionsBA = CurveCurve.intersectionXYPairs(curve1, false, curve0, false); + const intersectionsBA = CurveCurve.intersectionXYPairs(curve1, extend1, curve0, extend0); testSpiralIntersection(intersectionsBA, true); }; @@ -2540,9 +2541,9 @@ describe("CurveCurveIntersectXY", () => { const curveChain1 = CurveChainWithDistanceIndex.createCapture(path1); const bagOfCurves = BagOfCurves.create(path0, arc0, lineString0); - const curves: AnyCurve[] = [ - lineSegment1, + let curves: AnyCurve[] = [ lineSegment0, + lineSegment1, lineString0, arc0, arc1, @@ -2584,7 +2585,7 @@ describe("CurveCurveIntersectXY", () => { curveChain0, bagOfCurves, ]; - const numExpectedIntersections = [ + let numExpectedIntersections = [ 1, 1, 3, 2, 1, 1, // curve primitives other than spirals 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, // rotated and non-planar integrated spirals 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, // 1, 1, 1, 1, 1, 1, 1, 1 // rotated and non-planar direct spirals @@ -2592,14 +2593,14 @@ describe("CurveCurveIntersectXY", () => { ]; ck.testCoordinate(curves.length, numExpectedIntersections.length, "matching arrays"); // spiral vs all curves - const test0 = (spiral: TransitionSpiral3d) => { + const test0 = (spiral: TransitionSpiral3d, ddy = 0, extend = false) => { for (let j = 0; j < curves.length; j++) { const curve = curves[j]; const numExpectedIntersection = numExpectedIntersections[j]; - visualizeAndTestSpiralIntersection(ck, allGeometry, spiral, curve, numExpectedIntersection, dx, dy); + visualizeAndTestSpiralIntersection(ck, allGeometry, spiral, curve, numExpectedIntersection, dx, dy, false, extend); dy += 200; } - dy = 0; + dy = ddy; dx += 200; } for (let i = 0; i < integratedSpirals.length; i++) // skip rotated and non-planar integrated spirals @@ -2638,6 +2639,30 @@ describe("CurveCurveIntersectXY", () => { for (let i = 0; i < directSpirals.length; i++) // skip rotated and non-planar direct spirals if (i % 3 === 0) test1(directSpirals[i]); + // extend curve primitive + dx = 0; + dy += 200; + const lineSegment5 = LineSegment3d.create(Point3d.create(50, 100), Point3d.create(60, 50)); + const arc4 = Arc3d.createXY(Point3d.create(50, 10), 30, AngleSweep.createStartEndDegrees(90, 180)); + const lineString4 = LineString3d.create( + Point3d.create(20, 20), Point3d.create(20, 40), Point3d.create(40, 40), Point3d.create(40, 20), Point3d.create(70, 20), + ); + const lineString5 = LineString3d.create( + Point3d.create(10, 40), Point3d.create(10, 80), Point3d.create(30, 80), + ); + const arc5 = Arc3d.createXY(Point3d.create(30, 60), 20, AngleSweep.createStartEndDegrees(90, -90)); + const lineSegment6 = LineSegment3d.create(Point3d.create(30, 40), Point3d.create(50, 30)); + const path2 = Path.create(lineString5, arc5, lineSegment6); + curves = [lineSegment5, arc4, lineString4, path2]; + numExpectedIntersections = [1, 2, 2, 2]; + ck.testCoordinate(curves.length, numExpectedIntersections.length, "matching arrays"); + for (let i = 0; i < integratedSpirals.length; i++) // skip rotated and non-planar integrated spirals + if (i % 3 === 0) + test0(integratedSpirals[i], 7400, true); + dx += 250; + for (let i = 0; i < directSpirals.length; i++) // skip rotated and non-planar direct spirals + if (i % 3 === 0) + test0(directSpirals[i], 7400, true); GeometryCoreTestIO.saveGeometry(allGeometry, "CurveCurveIntersectXY", "SpiralIntersection"); expect(ck.getNumErrors()).toBe(0); From dfd4385ff179728f9b5b62620b2cdd61601150fb Mon Sep 17 00:00:00 2001 From: dassaf4 <68340676+dassaf4@users.noreply.github.com> Date: Fri, 21 Nov 2025 16:05:43 -0600 Subject: [PATCH 03/15] - make geomA/B arg order consistent - carefully handle `reverse` flag when computing additional seeds for extended curveA - `refineSpiralResultsByNewton`: instead of passing size of array tail, it is simpler to pass first index of tail, like in `CurveChainWithDistanceIndex.convertChildDetailToChainDetail` - consolidate duplicated code - add helpers for resetting geomA/B - increase maxIters for Newton to handle tangent intersections, which have sub-quadratic convergence - use appropriate tolerances for point/fraction comparison --- .../internalContexts/CurveCurveIntersectXY.ts | 279 ++++++++++-------- 1 file changed, 154 insertions(+), 125 deletions(-) diff --git a/core/geometry/src/curve/internalContexts/CurveCurveIntersectXY.ts b/core/geometry/src/curve/internalContexts/CurveCurveIntersectXY.ts index c2b2c2e03e2..c01cb74c82b 100644 --- a/core/geometry/src/curve/internalContexts/CurveCurveIntersectXY.ts +++ b/core/geometry/src/curve/internalContexts/CurveCurveIntersectXY.ts @@ -106,9 +106,29 @@ export class CurveCurveIntersectXY extends RecurseToCurvesGeometryHandler { this._coincidentGeometryContext = CoincidentGeometryQuery.create(tolerance); this._results = []; } - /** Reset the geometry, leaving all other parts unchanged (and preserving accumulated intersections). */ - public resetGeometry(geometryB: AnyCurve): void { - this._geometryB = geometryB; + /** + * Reset geometryA extension flags. + * * Undefined inputs are ignored. + * * All other instance data is unchanged, including accumulated intersections. + */ + public resetGeometryA(extendA0?: boolean, extendA1?: boolean): void { + if (extendA0 !== undefined) + this._extendA0 = extendA0; + if (extendA1 !== undefined) + this._extendA1 = extendA1; + } + /** + * Reset geometryB and its extension flags. + * * Undefined inputs are ignored. + * * All other instance data is unchanged, including accumulated intersections. + */ + public resetGeometryB(geometryB?: AnyCurve, extendB0?: boolean, extendB1?: boolean): void { + if (geometryB) + this._geometryB = geometryB; + if (extendB0 !== undefined) + this._extendB0 = extendB0; + if (extendB1 !== undefined) + this._extendB1 = extendB1; } /** Accept the fraction if it falls inside (possibly extended) fraction range. */ private acceptFraction(extend0: boolean, fraction: number, extend1: boolean, fractionTol: number = 1.0e-12): boolean { @@ -156,6 +176,7 @@ export class CurveCurveIntersectXY extends RecurseToCurvesGeometryHandler { * @param fractionB1 end of the subcurve of cpB * @param reversed whether to reverse the details in the recorded intersection pair * @param intervalDetails optional data for a coincident segment intersection + * @param fractionTolerance relative tolerance for comparing fractions to avoid duplicating the last intersection. Defaults to [[Geometry.smallAngleRadians]]. */ private recordPointWithLocalFractions( localFractionA: number, @@ -167,7 +188,8 @@ export class CurveCurveIntersectXY extends RecurseToCurvesGeometryHandler { fractionB0: number, fractionB1: number, reversed: boolean, - intervalDetails?: undefined | CurveLocationDetailPair, + intervalDetails?: CurveLocationDetailPair, + fractionTolerance?: number, ): void { let globalFractionA, globalFractionB; let globalFractionA1, globalFractionB1; @@ -189,12 +211,12 @@ export class CurveCurveIntersectXY extends RecurseToCurvesGeometryHandler { const oldDetailA = this._results[numPrevious - 1].detailA; const oldDetailB = this._results[numPrevious - 1].detailB; if (reversed) { - if (oldDetailB.isSameCurveAndFraction({ curve: cpA, fraction: globalFractionA }) && - oldDetailA.isSameCurveAndFraction({ curve: cpB, fraction: globalFractionB })) + if (oldDetailB.isSameCurveAndFraction({ curve: cpA, fraction: globalFractionA }, fractionTolerance) && + oldDetailA.isSameCurveAndFraction({ curve: cpB, fraction: globalFractionB }, fractionTolerance)) return; } else { - if (oldDetailA.isSameCurveAndFraction({ curve: cpA, fraction: globalFractionA }) && - oldDetailB.isSameCurveAndFraction({ curve: cpB, fraction: globalFractionB })) + if (oldDetailA.isSameCurveAndFraction({ curve: cpA, fraction: globalFractionA }, fractionTolerance) && + oldDetailB.isSameCurveAndFraction({ curve: cpB, fraction: globalFractionB }, fractionTolerance)) return; } } @@ -963,29 +985,29 @@ export class CurveCurveIntersectXY extends RecurseToCurvesGeometryHandler { * We take care of extend variables of geometry's children here if geometry is Path or Loop. */ private dispatchCurveCollection(geomA: AnyCurve, geomAHandler: (geomA: any) => any): void { - const geomB = this._geometryB; // save + const geomB = this._geometryB; if (!geomB || !geomB.children || !(geomB instanceof CurveCollection)) return; const children = geomB.children; - const extendB0 = this._extendB0; // save - const extendB1 = this._extendB1; // save + const saveExtendB0 = this._extendB0; + const saveExtendB1 = this._extendB1; for (let i = 0; i < children.length; i++) { - this.resetGeometry(children[i]); + let extendB0 = saveExtendB0; + let extendB1 = saveExtendB1; if (geomB instanceof Path && children.length > 1) { if (i === 0) - this._extendB1 = false; // first child can only extend from start + extendB1 = false; // first child can only extend from start else if (i === children.length - 1) - this._extendB0 = false; // last child can only extend from end + extendB0 = false; // last child can only extend from end else - this._extendB0 = this._extendB1 = false; // middle children cannot extend + extendB0 = extendB1 = false; // middle children cannot extend } else if (geomB instanceof Loop) { - this._extendB0 = this._extendB1 = false; // Loops cannot extend + extendB0 = extendB1 = false; // Loops cannot extend } + this.resetGeometryB(children[i], extendB0, extendB1); geomAHandler(geomA); - this._extendB0 = extendB0; // restore - this._extendB1 = extendB1; // restore } - this.resetGeometry(geomB); // restore + this.resetGeometryB(geomB, saveExtendB0, saveExtendB1); } /** * Low level dispatch of CurveChainWithDistanceIndex. @@ -999,37 +1021,40 @@ export class CurveCurveIntersectXY extends RecurseToCurvesGeometryHandler { } const index0 = this._results.length; const geomB = this._geometryB; // save - this.resetGeometry(geomB.path); + this.resetGeometryB(geomB.path); this.dispatchCurveCollection(geomA, geomAHandler); - this.resetGeometry(geomB); // restore + this.resetGeometryB(geomB); // restore this._results = CurveChainWithDistanceIndex.convertChildDetailToChainDetail( this._results, index0, undefined, geomB, true, ); } /** - * Invoke `child.dispatchToGeometryHandler(this)` for each child in the array returned by the query `g.children`. - * We take care of extend variables of geometry's children here if geometry is Path or Loop. + * Invoke dispatch on each child of `g` as "geometryA". + * * If `g` is a `Path` or `Loop`, adjust extension flags for geometryA accordingly. */ public override handleChildren(g: GeometryQuery): any { const children = g.children; - const extendA0 = this._extendA0; // save - const extendA1 = this._extendA1; // save - if (children) - for (let i = 0; i < children.length; i++) { - if (g instanceof Path && children.length > 1) { - if (i === 0) - this._extendA1 = false; // first child can only extend from start - else if (i === children.length - 1) - this._extendA0 = false; // last child can only extend from end - else - this._extendA0 = this._extendA1 = false; // middle children cannot extend - } else if (g instanceof Loop) { - this._extendA0 = this._extendA1 = false; // Loops cannot extend - } - children[i].dispatchToGeometryHandler(this); - this._extendA0 = extendA0; // restore - this._extendA1 = extendA1; // restore + if (!children) + return; + const saveExtendA0 = this._extendA0; + const saveExtendA1 = this._extendA1; + for (let i = 0; i < children.length; i++) { + let extendA0 = saveExtendA0; + let extendA1 = saveExtendA1; + if (g instanceof Path && children.length > 1) { + if (i === 0) + extendA1 = false; // first child can only extend from start + else if (i === children.length - 1) + extendA0 = false; // last child can only extend from end + else + extendA0 = extendA1 = false; // middle children cannot extend + } else if (g instanceof Loop) { + extendA0 = extendA1 = false; // Loops cannot extend } + this.resetGeometryA(extendA0, extendA1); + children[i].dispatchToGeometryHandler(this); + } + this.resetGeometryA(saveExtendA0, saveExtendA1); } /** Double dispatch handler for strongly typed segment. */ public override handleLineSegment3d(segmentA: LineSegment3d): any { @@ -1055,14 +1080,7 @@ export class CurveCurveIntersectXY extends RecurseToCurvesGeometryHandler { this._geometryB, this._extendB0, this._extendB1, false, ); } else if (this._geometryB instanceof TransitionSpiral3d) { - const spiralApproximation = LineString3d.create(); - this._geometryB.emitStrokes(spiralApproximation); - const numPreviousResults = this._results.length; - this._results.push(...CurveCurve.closeApproachProjectedXYPairs(segmentA, spiralApproximation, Geometry.smallMetricDistance)); - if (this._extendA0 || this._extendA1) - this.computeSegmentLineString(segmentA, this._extendA0, this._extendA1, spiralApproximation, false, false, false); - const numberOfNewResults = this._results.length - numPreviousResults; - this.refineSpiralResultsByNewton(segmentA, this._geometryB, numberOfNewResults); + this.dispatchCurveSpiral(segmentA, this._extendA0, this._extendA1, this._geometryB, false); } else if (this._geometryB instanceof CurveCollection) { this.dispatchCurveCollection(segmentA, this.handleLineSegment3d.bind(this)); } else if (this._geometryB instanceof CurveChainWithDistanceIndex) { @@ -1088,14 +1106,7 @@ export class CurveCurveIntersectXY extends RecurseToCurvesGeometryHandler { lsA, this._extendA0, this._extendA1, this._geometryB, this._extendB0, this._extendB1, false, ); } else if (this._geometryB instanceof TransitionSpiral3d) { - const spiralApproximation = LineString3d.create(); - this._geometryB.emitStrokes(spiralApproximation); - const numPreviousResults = this._results.length; - this._results.push(...CurveCurve.closeApproachProjectedXYPairs(lsA, spiralApproximation, Geometry.smallMetricDistance)); - if (this._extendA0 || this._extendA1) - this.computeLineStringLineString(lsA, this._extendA0, this._extendA1, spiralApproximation, false, false, false); - const numberOfNewResults = this._results.length - numPreviousResults; - this.refineSpiralResultsByNewton(lsA, this._geometryB, numberOfNewResults); + this.dispatchCurveSpiral(lsA, this._extendA0, this._extendA1, this._geometryB, false); } else if (this._geometryB instanceof CurveCollection) { this.dispatchCurveCollection(lsA, this.handleLineString3d.bind(this)); } else if (this._geometryB instanceof CurveChainWithDistanceIndex) { @@ -1104,124 +1115,142 @@ export class CurveCurveIntersectXY extends RecurseToCurvesGeometryHandler { return undefined; } /** Double dispatch handler for strongly typed arc. */ - public override handleArc3d(arc0: Arc3d): any { + public override handleArc3d(arcA: Arc3d): any { if (this._geometryB instanceof LineSegment3d) { this.dispatchSegmentArc( this._geometryB, this._extendB0, this._geometryB.point0Ref, 0.0, this._geometryB.point1Ref, 1.0, this._extendB1, - arc0, this._extendA0, this._extendA1, true, + arcA, this._extendA0, this._extendA1, true, ); } else if (this._geometryB instanceof LineString3d) { this.computeArcLineString( - arc0, this._extendA0, this._extendA1, this._geometryB, this._extendB0, this._extendB1, false, + arcA, this._extendA0, this._extendA1, this._geometryB, this._extendB0, this._extendB1, false, ); } else if (this._geometryB instanceof Arc3d) { - this.dispatchArcArc(arc0, this._extendA0, this._extendA1, this._geometryB, this._extendB0, this._extendB1, false); + this.dispatchArcArc(arcA, this._extendA0, this._extendA1, this._geometryB, this._extendB0, this._extendB1, false); } else if (this._geometryB instanceof BSplineCurve3d) { this.dispatchArcBsplineCurve3d( - arc0, this._extendA0, this._extendA1, this._geometryB, this._extendB0, this._extendB1, false, + arcA, this._extendA0, this._extendA1, this._geometryB, this._extendB0, this._extendB1, false, ); } else if (this._geometryB instanceof TransitionSpiral3d) { - const spiralApproximation = LineString3d.create(); - this._geometryB.emitStrokes(spiralApproximation); - const numPreviousResults = this._results.length; - this._results.push(...CurveCurve.closeApproachProjectedXYPairs(arc0, spiralApproximation, Geometry.smallMetricDistance)); - if (this._extendA0 || this._extendA1) - this.computeArcLineString(arc0, this._extendA0, this._extendA1, spiralApproximation, false, false, false); - const numberOfNewResults = this._results.length - numPreviousResults; - this.refineSpiralResultsByNewton(arc0, this._geometryB, numberOfNewResults); + this.dispatchCurveSpiral(arcA, this._extendA0, this._extendA1, this._geometryB, false); } else if (this._geometryB instanceof CurveCollection) { - this.dispatchCurveCollection(arc0, this.handleArc3d.bind(this)); + this.dispatchCurveCollection(arcA, this.handleArc3d.bind(this)); } else if (this._geometryB instanceof CurveChainWithDistanceIndex) { - this.dispatchCurveChainWithDistanceIndex(arc0, this.handleArc3d.bind(this)); + this.dispatchCurveChainWithDistanceIndex(arcA, this.handleArc3d.bind(this)); } return undefined; } /** Double dispatch handler for strongly typed bspline curve. */ - public override handleBSplineCurve3d(curve: BSplineCurve3d): any { + public override handleBSplineCurve3d(curveA: BSplineCurve3d): any { if (this._geometryB instanceof LineSegment3d) { this.dispatchSegmentBsplineCurve( this._geometryB, this._extendB0, this._geometryB.point0Ref, 0.0, this._geometryB.point1Ref, 1.0, this._extendB1, - curve, this._extendA0, this._extendA1, true, + curveA, this._extendA0, this._extendA1, true, ); } else if (this._geometryB instanceof LineString3d) { this.dispatchLineStringBSplineCurve( - this._geometryB, this._extendB0, this._extendB1, curve, this._extendA0, this._extendA1, true, + this._geometryB, this._extendB0, this._extendB1, curveA, this._extendA0, this._extendA1, true, ); } else if (this._geometryB instanceof Arc3d) { this.dispatchArcBsplineCurve3d( - this._geometryB, this._extendB0, this._extendB1, curve, this._extendA0, this._extendA1, true, + this._geometryB, this._extendB0, this._extendB1, curveA, this._extendA0, this._extendA1, true, ); } else if (this._geometryB instanceof BSplineCurve3dBase) { - this.dispatchBSplineCurve3dBSplineCurve3d(curve, this._geometryB, false); + this.dispatchBSplineCurve3dBSplineCurve3d(curveA, this._geometryB, false); } else if (this._geometryB instanceof TransitionSpiral3d) { - const spiralApproximation = LineString3d.create(); - this._geometryB.emitStrokes(spiralApproximation); - const numPreviousResults = this._results.length; - this._results.push(...CurveCurve.closeApproachProjectedXYPairs(curve, spiralApproximation, Geometry.smallMetricDistance)); - if (this._extendA0 || this._extendA1) - this.dispatchLineStringBSplineCurve(spiralApproximation, false, false, curve, this._extendA0, this._extendA1, true); - const numberOfNewResults = this._results.length - numPreviousResults; - this.refineSpiralResultsByNewton(curve, this._geometryB, numberOfNewResults); + this.dispatchCurveSpiral(curveA, this._extendA0, this._extendA1, this._geometryB, false); } else if (this._geometryB instanceof CurveCollection) { - this.dispatchCurveCollection(curve, this.handleBSplineCurve3d.bind(this)); + this.dispatchCurveCollection(curveA, this.handleBSplineCurve3d.bind(this)); } else if (this._geometryB instanceof CurveChainWithDistanceIndex) { - this.dispatchCurveChainWithDistanceIndex(curve, this.handleBSplineCurve3d.bind(this)); + this.dispatchCurveChainWithDistanceIndex(curveA, this.handleBSplineCurve3d.bind(this)); } return undefined; } /** - * Assuming the tail of `this._results` holds the given number of approximate solutions to the curve-spiral - * xy-intersection problem, replace each with its Newton refinement, unless it doesn't converge, in which case - * it is removed. - * @param spiral The transition spiral. - * @param otherCurve The other curve primitive. - * @param numberOfNewResults The number of results in the tail of `this._results` to be refined. - * @param reversed Whether the spiral is geometryA (true) or geometryB (false). + * Process tail of `this._results` for xy-intersections between the curve and spiral. + * * If a result is not already an intersection, refine it via Newton iteration unless it doesn't converge, in which case remove it. + * @param curveA The other curve primitive. May also be a transition spiral. + * @param spiralB The transition spiral. + * @param index0 index of first entry in tail of `this._results` to refine. + * @param reversed Whether `spiralB` is geometryA (true) or geometryB (false). */ - private refineSpiralResultsByNewton( - otherCurve: CurvePrimitive, spiral: TransitionSpiral3d, numberOfNewResults: number, reversed = false, - ): void { - const resultsToBeRefined = this._results.slice(this._results.length - numberOfNewResults); - this._results.length -= numberOfNewResults; // keep already refined results + private refineSpiralResultsByNewton(curveA: CurvePrimitive, spiralB: TransitionSpiral3d, index0: number, reversed = false): void { + const resultsToBeRefined = this._results.slice(index0); + this._results.length -= resultsToBeRefined.length; // keep already refined results for (const detail of resultsToBeRefined) { - let spiralFraction = reversed ? detail.detailA.fraction : detail.detailB.fraction; - let otherFraction = reversed ? detail.detailB.fraction : detail.detailA.fraction; - const otherExtend0 = reversed ? this._extendB0 : this._extendA0; - const otherExtend1 = reversed ? this._extendB1 : this._extendA1; - if (spiral.fractionToPoint(spiralFraction).isAlmostEqualXY(otherCurve.fractionToPoint(otherFraction))) { // already an accurate intersection - this.recordPointWithLocalFractions(otherFraction, otherCurve, 0, 1, spiralFraction, spiral, 0, 1, reversed); + const detailA = reversed ? detail.detailB : detail.detailA; + const detailB = reversed ? detail.detailA : detail.detailB; + assert(detailB.curve instanceof LineString3d, "Caller has discretized the spiral"); + let fractionA = detailA.fraction; + // convert linestring parameter to fractional length to more closely agree with spiral parameter + let fractionB = Geometry.safeDivideFraction(detailB.curve.curveLengthBetweenFractions(0, detailB.fraction), detailB.curve.curveLength(), detailB.fraction); + // Newton iteration fails when the (partial) derivative is small. Check for a root there to avoid missing it. + if (spiralB.fractionToPoint(fractionB).isAlmostEqualXY(detailA.point, this._coincidentGeometryContext.tolerance)) { + this.recordPointWithLocalFractions(fractionA, curveA, 0, 1, fractionB, spiralB, 0, 1, reversed); continue; } - const xyMatchingFunction = new CurveCurveIntersectionXYRRToRRD(spiral, otherCurve); - const newtonSearcher = new Newton2dUnboundedWithDerivative(xyMatchingFunction); - newtonSearcher.setUV(spiralFraction, otherFraction); + const extendA0 = reversed ? this._extendB0 : this._extendA0; + const extendA1 = reversed ? this._extendB1 : this._extendA1; + const xyMatchingFunction = new CurveCurveIntersectionXYRRToRRD(curveA, spiralB); + const maxIterations = 64; // observed: 56 iterations for tangent intersection (linear convergence) + const newtonSearcher = new Newton2dUnboundedWithDerivative(xyMatchingFunction, maxIterations); + const fractionTol = 2 * newtonSearcher.stepSizeTolerance; // cluster diameter for Newton solutions + newtonSearcher.setUV(fractionA, fractionB); if (newtonSearcher.runIterations()) { - spiralFraction = newtonSearcher.getU(); - otherFraction = newtonSearcher.getV(); - if (this.acceptFraction(false, spiralFraction, false) && this.acceptFraction(otherExtend0, otherFraction, otherExtend1)) - this.recordPointWithLocalFractions(otherFraction, otherCurve, 0, 1, spiralFraction, spiral, 0, 1, reversed); + fractionA = newtonSearcher.getU(); + fractionB = newtonSearcher.getV(); + if (this.acceptFraction(extendA0, fractionA, extendA1) && this.acceptFraction(false, fractionB, false)) + this.recordPointWithLocalFractions(fractionA, curveA, 0, 1, fractionB, spiralB, 0, 1, reversed, undefined, fractionTol); + } + } + } + /** + * Compute the xy-intersection of a curve and a spiral. + * * When `curveA` is extended, duplicate solutions may be returned. + * @param curveA curve to intersect with spiralB. May also be a transition spiral. + * @param extendA0 whether to compute xy-intersections with curveA extended beyond its start + * @param extendA1 whether to compute xy-intersections with curveA extended beyond its end + * @param spiralB transition spiral to intersect with curveA + * @param reversed whether spiralB data is returned in detailA of each pair + */ + private dispatchCurveSpiral(curveA: CurvePrimitive, extendA0: boolean, extendA1: boolean, spiralB: TransitionSpiral3d, reversed: boolean): void { + const index0 = this._results.length; + // approximate spiral(s) + const cpA = curveA instanceof TransitionSpiral3d ? LineString3d.create() : curveA; + const cpB = LineString3d.create(); + spiralB.emitStrokes(cpB); + let spiralLength = spiralB.quickLength(); + if (curveA instanceof TransitionSpiral3d && cpA instanceof LineString3d) { + curveA.emitStrokes(cpA); + spiralLength += curveA.quickLength(); + } + // append seeds using unextended curveA and stroked spiral + const maxDistance = Geometry.smallMetricDistance * (1 + 100 * spiralLength); // HEURISTIC: not too tight + this._results.push(...CurveCurve.closeApproachProjectedXYPairs(reversed ? cpB : cpA, reversed ? cpA : cpB, maxDistance)); + // append seeds using extended curveA and stroked spiral + if (curveA.isExtensibleFractionSpace && (extendA0 || extendA1) && this._geometryB) { + const geomB = this._geometryB; + const extendB0 = this._extendB0; + const extendB1 = this._extendB1; + const i0 = this._results.length; + this.resetGeometryA(false, false); + this.resetGeometryB(curveA, extendA0, extendA1); + this.handleLineString3d(cpB); // puts curveA results in detailB + if (!reversed) { + for (let i = i0; i < this._results.length; i++) + this._results[i].swapDetails(); // now curveA results are in detailA } + this.resetGeometryA(extendA0, extendA1); + this.resetGeometryB(geomB, extendB0, extendB1); } + this.refineSpiralResultsByNewton(curveA, spiralB, index0, reversed); } /** Double dispatch handler for strongly typed spiral curve. */ public override handleTransitionSpiral(spiral: TransitionSpiral3d): any { if (this._geometryB instanceof CurveChainWithDistanceIndex) { this.dispatchCurveChainWithDistanceIndex(spiral, this.handleTransitionSpiral.bind(this)); } else if (this._geometryB instanceof CurvePrimitive) { - const spiralApproximation = LineString3d.create(); - spiral.emitStrokes(spiralApproximation); - const geomBApproximation = LineString3d.create(); - if (this._geometryB instanceof TransitionSpiral3d) - this._geometryB.emitStrokes(geomBApproximation); - const numPreviousResults = this._results.length; - if (this._geometryB instanceof TransitionSpiral3d) - this._results.push(...CurveCurve.closeApproachProjectedXYPairs(spiralApproximation, geomBApproximation, Geometry.smallMetricDistance)); - else - this._results.push(...CurveCurve.closeApproachProjectedXYPairs(spiralApproximation, this._geometryB, Geometry.smallMetricDistance)); - if (this._extendB0 || this._extendB1) - this.handleLineString3d(spiralApproximation); - const numberOfNewResults = this._results.length - numPreviousResults; - this.refineSpiralResultsByNewton(this._geometryB, spiral, numberOfNewResults, true); + this.dispatchCurveSpiral(this._geometryB, this._extendB0, this._extendB1, spiral, true); } else if (this._geometryB instanceof CurveCollection) { this.dispatchCurveCollection(spiral, this.handleTransitionSpiral.bind(this)); } From ea3566188c271fd4d516fa6765feca5746c9d589 Mon Sep 17 00:00:00 2001 From: dassaf4 <68340676+dassaf4@users.noreply.github.com> Date: Fri, 21 Nov 2025 16:06:27 -0600 Subject: [PATCH 04/15] use `Checker.testExactNumber` for 0-tolerance comparisons (e.g., for counts) --- .../geometry/src/test/curve/CurveCurveIntersectXY.test.ts | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/core/geometry/src/test/curve/CurveCurveIntersectXY.test.ts b/core/geometry/src/test/curve/CurveCurveIntersectXY.test.ts index d8e7bf11beb..071b981b519 100644 --- a/core/geometry/src/test/curve/CurveCurveIntersectXY.test.ts +++ b/core/geometry/src/test/curve/CurveCurveIntersectXY.test.ts @@ -2357,7 +2357,7 @@ describe("CurveCurveIntersectXY", () => { GeometryCoreTestIO.captureCloneGeometry(allGeometry, [circle, curve]); // this requires a call to UnivariateBezier.runNewton that needs 20 iterations to converge! const intersections = CurveCurve.intersectionProjectedXYPairs(undefined, circle, false, curve, false); - if (ck.testCoordinate(1, intersections.length, "Expected one intersection between the circle and the bspline")) { + if (ck.testExactNumber(1, intersections.length, "Expected one intersection between the circle and the bspline")) { const intersection = intersections[0].detailA.point; if (ck.testPoint3d(intersection, intersections[0].detailB.point, "report same intersection point on both curves")) { GeometryCoreTestIO.createAndCaptureXYCircle(allGeometry, intersection, 0.0003); @@ -2382,7 +2382,7 @@ describe("CurveCurveIntersectXY", () => { GeometryCoreTestIO.captureCurveLocationDetails(allGeometry, intersections, 5, dx, dy); const curveName0 = curve0.constructor.name; const curveName1 = curve1.constructor.name; - ck.testCoordinate( + ck.testExactNumber( numExpected, intersections.length, `expect ${numExpected} intersection(s) between ${curveName1} and ${curveName0}`, @@ -2591,7 +2591,7 @@ describe("CurveCurveIntersectXY", () => { 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, // 1, 1, 1, 1, 1, 1, 1, 1 // rotated and non-planar direct spirals 5, 6, 5, 10, // path, loop, curve chain, and bag of curves ]; - ck.testCoordinate(curves.length, numExpectedIntersections.length, "matching arrays"); + ck.testExactNumber(curves.length, numExpectedIntersections.length, "matching arrays"); // spiral vs all curves const test0 = (spiral: TransitionSpiral3d, ddy = 0, extend = false) => { for (let j = 0; j < curves.length; j++) { @@ -2655,7 +2655,7 @@ describe("CurveCurveIntersectXY", () => { const path2 = Path.create(lineString5, arc5, lineSegment6); curves = [lineSegment5, arc4, lineString4, path2]; numExpectedIntersections = [1, 2, 2, 2]; - ck.testCoordinate(curves.length, numExpectedIntersections.length, "matching arrays"); + ck.testExactNumber(curves.length, numExpectedIntersections.length, "matching arrays"); for (let i = 0; i < integratedSpirals.length; i++) // skip rotated and non-planar integrated spirals if (i % 3 === 0) test0(integratedSpirals[i], 7400, true); From a05819cc676fcf564563bf930ac6040a6d31c08e Mon Sep 17 00:00:00 2001 From: dassaf4 <68340676+dassaf4@users.noreply.github.com> Date: Fri, 21 Nov 2025 16:06:49 -0600 Subject: [PATCH 05/15] helper methods --- core/geometry/src/Geometry.ts | 9 +++++++++ core/geometry/src/curve/CurveCurve.ts | 2 +- core/geometry/src/curve/CurveLocationDetail.ts | 9 ++++++--- core/geometry/src/numerics/Newton.ts | 4 ++++ core/geometry/src/test/Checker.ts | 4 ++-- 5 files changed, 22 insertions(+), 6 deletions(-) diff --git a/core/geometry/src/Geometry.ts b/core/geometry/src/Geometry.ts index 09b57803f67..05ac9d1c2cd 100644 --- a/core/geometry/src/Geometry.ts +++ b/core/geometry/src/Geometry.ts @@ -353,6 +353,15 @@ export class Geometry { d = -d; return d <= tolerance; } + /** + * Toleranced test for equivalent fractions. + * @param x first fraction to compare + * @param y second fraction to compare + * @param tolerance maximum difference between fractions considered equivalent, defaulting to [[smallFraction]]. + */ + public static isSameFraction(x: number, y: number, tolerance: number = Geometry.smallFraction): boolean { + return this.isSameCoordinate(x, y, tolerance); + } /** * Boolean test for metric coordinate near-equality (i.e., if `x` and `y` are almost equal) using * `tolerance = toleranceFactor * smallMetricDistance` diff --git a/core/geometry/src/curve/CurveCurve.ts b/core/geometry/src/curve/CurveCurve.ts index 4c73a173baf..7293e1e06eb 100644 --- a/core/geometry/src/curve/CurveCurve.ts +++ b/core/geometry/src/curve/CurveCurve.ts @@ -89,7 +89,7 @@ export class CurveCurve { for (let i = 0; i < primitives.length; i++) { const curveA = primitives[i]; for (let j = i + 1; j < primitives.length; j++) { - handler.resetGeometry(primitives[j]); + handler.resetGeometryB(primitives[j]); curveA.dispatchToGeometryHandler(handler); } } diff --git a/core/geometry/src/curve/CurveLocationDetail.ts b/core/geometry/src/curve/CurveLocationDetail.ts index 3810735b902..3223d9d5728 100644 --- a/core/geometry/src/curve/CurveLocationDetail.ts +++ b/core/geometry/src/curve/CurveLocationDetail.ts @@ -398,9 +398,12 @@ export class CurveLocationDetail { } return detailB; } - /** Compare only the curve and fraction of this detail with `other`. */ - public isSameCurveAndFraction(other: CurveLocationDetail | { curve: CurvePrimitive, fraction: number }): boolean { - return this.curve === other.curve && Geometry.isAlmostEqualNumber(this.fraction, other.fraction); + /** + * Compare only the curve and fraction of this detail with `other`. + * @param tolerance optional relative tolerance for comparing fractions with [[Geometry.isAlmostEqualNumber]]. + */ + public isSameCurveAndFraction(other: CurveLocationDetail | { curve: CurvePrimitive, fraction: number }, tolerance?: number): boolean { + return this.curve === other.curve && Geometry.isAlmostEqualNumber(this.fraction, other.fraction, tolerance); } /** * Transform the detail in place. diff --git a/core/geometry/src/numerics/Newton.ts b/core/geometry/src/numerics/Newton.ts index 8ce480e3b43..c01b63b965b 100644 --- a/core/geometry/src/numerics/Newton.ts +++ b/core/geometry/src/numerics/Newton.ts @@ -316,6 +316,10 @@ export class Newton2dUnboundedWithDerivative extends AbstractNewtonIterator { public getV(): number { return this._currentUV.y; } + /** Get the relative tolerance for comparing iterations in [[testConvergence]]. */ + public get stepSizeTolerance(): number { + return this._stepSizeTolerance; + } /** Update the current uv parameter by currentStep, i.e., compute `X_{n+1} := X_n - dX = (u_n - du, v_n - dv)`. */ public applyCurrentStep(): boolean { // console.log("(" + (this._currentUV.x - this._currentStep.x) + "," + (this._currentUV.y - this._currentStep.y) + ")"); diff --git a/core/geometry/src/test/Checker.ts b/core/geometry/src/test/Checker.ts index c8126878d50..3c84489112c 100644 --- a/core/geometry/src/test/Checker.ts +++ b/core/geometry/src/test/Checker.ts @@ -370,9 +370,9 @@ export class Checker { * * See also [[testExactNumber]], [[testNearNumber]], [[testSmallRelative]], [[testCoordinate]], [[testCoordinateWithToleranceFactor]] */ public testFraction(dataA: number, dataB: number, ...params: any[]): boolean { - if (Geometry.isSameCoordinate(dataA, dataB, Geometry.smallFraction)) + if (Geometry.isSameFraction(dataA, dataB)) return this.announceOK(); - return this.announceError("Expect same coordinate", dataA, dataB, params); + return this.announceError("Expect same fraction", dataA, dataB, params); } /** * Return true if 2 numbers are almost equal within default metric tolerance. From 7d88e0c262f3d66de431d06993b38769066f2a80 Mon Sep 17 00:00:00 2001 From: dassaf4 <68340676+dassaf4@users.noreply.github.com> Date: Fri, 21 Nov 2025 16:07:07 -0600 Subject: [PATCH 06/15] use linestring fraction as spiral fraction seed; it evaluates to a closer spiral point compared to using the linestring fractional arclength as spiral fraction seed. --- .../geometry/src/curve/CurveLocationDetail.ts | 6 +-- .../internalContexts/CurveCurveIntersectXY.ts | 46 +++++++++---------- 2 files changed, 26 insertions(+), 26 deletions(-) diff --git a/core/geometry/src/curve/CurveLocationDetail.ts b/core/geometry/src/curve/CurveLocationDetail.ts index 3223d9d5728..f81cde79000 100644 --- a/core/geometry/src/curve/CurveLocationDetail.ts +++ b/core/geometry/src/curve/CurveLocationDetail.ts @@ -400,10 +400,10 @@ export class CurveLocationDetail { } /** * Compare only the curve and fraction of this detail with `other`. - * @param tolerance optional relative tolerance for comparing fractions with [[Geometry.isAlmostEqualNumber]]. + * @param fractionTol optional relative tolerance for comparing fractions with [[Geometry.isAlmostEqualNumber]]. */ - public isSameCurveAndFraction(other: CurveLocationDetail | { curve: CurvePrimitive, fraction: number }, tolerance?: number): boolean { - return this.curve === other.curve && Geometry.isAlmostEqualNumber(this.fraction, other.fraction, tolerance); + public isSameCurveAndFraction(other: CurveLocationDetail | { curve: CurvePrimitive, fraction: number }, fractionTol?: number): boolean { + return this.curve === other.curve && Geometry.isAlmostEqualNumber(this.fraction, other.fraction, fractionTol); } /** * Transform the detail in place. diff --git a/core/geometry/src/curve/internalContexts/CurveCurveIntersectXY.ts b/core/geometry/src/curve/internalContexts/CurveCurveIntersectXY.ts index c01cb74c82b..31b98e426d9 100644 --- a/core/geometry/src/curve/internalContexts/CurveCurveIntersectXY.ts +++ b/core/geometry/src/curve/internalContexts/CurveCurveIntersectXY.ts @@ -176,7 +176,7 @@ export class CurveCurveIntersectXY extends RecurseToCurvesGeometryHandler { * @param fractionB1 end of the subcurve of cpB * @param reversed whether to reverse the details in the recorded intersection pair * @param intervalDetails optional data for a coincident segment intersection - * @param fractionTolerance relative tolerance for comparing fractions to avoid duplicating the last intersection. Defaults to [[Geometry.smallAngleRadians]]. + * @param fractionTol relative tolerance for comparing fractions to avoid duplicating the last intersection. Defaults to [[Geometry.smallAngleRadians]]. */ private recordPointWithLocalFractions( localFractionA: number, @@ -189,7 +189,7 @@ export class CurveCurveIntersectXY extends RecurseToCurvesGeometryHandler { fractionB1: number, reversed: boolean, intervalDetails?: CurveLocationDetailPair, - fractionTolerance?: number, + fractionTol?: number, ): void { let globalFractionA, globalFractionB; let globalFractionA1, globalFractionB1; @@ -211,12 +211,12 @@ export class CurveCurveIntersectXY extends RecurseToCurvesGeometryHandler { const oldDetailA = this._results[numPrevious - 1].detailA; const oldDetailB = this._results[numPrevious - 1].detailB; if (reversed) { - if (oldDetailB.isSameCurveAndFraction({ curve: cpA, fraction: globalFractionA }, fractionTolerance) && - oldDetailA.isSameCurveAndFraction({ curve: cpB, fraction: globalFractionB }, fractionTolerance)) + if (oldDetailB.isSameCurveAndFraction({ curve: cpA, fraction: globalFractionA }, fractionTol) && + oldDetailA.isSameCurveAndFraction({ curve: cpB, fraction: globalFractionB }, fractionTol)) return; } else { - if (oldDetailA.isSameCurveAndFraction({ curve: cpA, fraction: globalFractionA }, fractionTolerance) && - oldDetailB.isSameCurveAndFraction({ curve: cpB, fraction: globalFractionB }, fractionTolerance)) + if (oldDetailA.isSameCurveAndFraction({ curve: cpA, fraction: globalFractionA }, fractionTol) && + oldDetailB.isSameCurveAndFraction({ curve: cpB, fraction: globalFractionB }, fractionTol)) return; } } @@ -1172,18 +1172,17 @@ export class CurveCurveIntersectXY extends RecurseToCurvesGeometryHandler { * @param curveA The other curve primitive. May also be a transition spiral. * @param spiralB The transition spiral. * @param index0 index of first entry in tail of `this._results` to refine. - * @param reversed Whether `spiralB` is geometryA (true) or geometryB (false). + * @param reversed Whether `spiralB` data is in `detailA` of each recorded pair, and `curveA` data in `detailB`. */ private refineSpiralResultsByNewton(curveA: CurvePrimitive, spiralB: TransitionSpiral3d, index0: number, reversed = false): void { const resultsToBeRefined = this._results.slice(index0); this._results.length -= resultsToBeRefined.length; // keep already refined results - for (const detail of resultsToBeRefined) { - const detailA = reversed ? detail.detailB : detail.detailA; - const detailB = reversed ? detail.detailA : detail.detailB; + for (const pair of resultsToBeRefined) { + const detailA = reversed ? pair.detailB : pair.detailA; + const detailB = reversed ? pair.detailA : pair.detailB; assert(detailB.curve instanceof LineString3d, "Caller has discretized the spiral"); let fractionA = detailA.fraction; - // convert linestring parameter to fractional length to more closely agree with spiral parameter - let fractionB = Geometry.safeDivideFraction(detailB.curve.curveLengthBetweenFractions(0, detailB.fraction), detailB.curve.curveLength(), detailB.fraction); + let fractionB = detailB.fraction; // use linestring fraction as the spiral fraction; it is closer than fractional length! // Newton iteration fails when the (partial) derivative is small. Check for a root there to avoid missing it. if (spiralB.fractionToPoint(fractionB).isAlmostEqualXY(detailA.point, this._coincidentGeometryContext.tolerance)) { this.recordPointWithLocalFractions(fractionA, curveA, 0, 1, fractionB, spiralB, 0, 1, reversed); @@ -1192,9 +1191,9 @@ export class CurveCurveIntersectXY extends RecurseToCurvesGeometryHandler { const extendA0 = reversed ? this._extendB0 : this._extendA0; const extendA1 = reversed ? this._extendB1 : this._extendA1; const xyMatchingFunction = new CurveCurveIntersectionXYRRToRRD(curveA, spiralB); - const maxIterations = 64; // observed: 56 iterations for tangent intersection (linear convergence) + const maxIterations = 64; // observed: 55 iterations for tangent intersection (linear convergence) const newtonSearcher = new Newton2dUnboundedWithDerivative(xyMatchingFunction, maxIterations); - const fractionTol = 2 * newtonSearcher.stepSizeTolerance; // cluster diameter for Newton solutions + const fractionTol = 2 * newtonSearcher.stepSizeTolerance; // relative cluster diameter for Newton convergence newtonSearcher.setUV(fractionA, fractionB); if (newtonSearcher.runIterations()) { fractionA = newtonSearcher.getU(); @@ -1208,11 +1207,11 @@ export class CurveCurveIntersectXY extends RecurseToCurvesGeometryHandler { * Compute the xy-intersection of a curve and a spiral. * * When `curveA` is extended, duplicate solutions may be returned. * @param curveA curve to intersect with spiralB. May also be a transition spiral. - * @param extendA0 whether to compute xy-intersections with curveA extended beyond its start - * @param extendA1 whether to compute xy-intersections with curveA extended beyond its end - * @param spiralB transition spiral to intersect with curveA - * @param reversed whether spiralB data is returned in detailA of each pair - */ + * @param extendA0 whether to compute xy-intersections with curveA extended beyond its start. + * @param extendA1 whether to compute xy-intersections with curveA extended beyond its end. + * @param spiralB transition spiral to intersect with curveA. + * @param reversed Whether `spiralB` data is in `detailA` of each recorded pair, and `curveA` data in `detailB`. + */ private dispatchCurveSpiral(curveA: CurvePrimitive, extendA0: boolean, extendA1: boolean, spiralB: TransitionSpiral3d, reversed: boolean): void { const index0 = this._results.length; // approximate spiral(s) @@ -1233,12 +1232,13 @@ export class CurveCurveIntersectXY extends RecurseToCurvesGeometryHandler { const extendB0 = this._extendB0; const extendB1 = this._extendB1; const i0 = this._results.length; - this.resetGeometryA(false, false); - this.resetGeometryB(curveA, extendA0, extendA1); - this.handleLineString3d(cpB); // puts curveA results in detailB + // handleLineString3d requires us to swap geometries: + this.resetGeometryA(false, false); // spiral approximation cpB is sent in as geometryA + this.resetGeometryB(curveA, extendA0, extendA1); // curveA is installed as geometryB + this.handleLineString3d(cpB); // this puts spiral data in detailA (expected when reversed is true) if (!reversed) { for (let i = i0; i < this._results.length; i++) - this._results[i].swapDetails(); // now curveA results are in detailA + this._results[i].swapDetails(); // now spiral data is in detailB } this.resetGeometryA(extendA0, extendA1); this.resetGeometryB(geomB, extendB0, extendB1); From acce963a79d75c1d0a45981ebd4a3f002de97a78 Mon Sep 17 00:00:00 2001 From: dassaf4 <68340676+dassaf4@users.noreply.github.com> Date: Fri, 21 Nov 2025 16:54:44 -0600 Subject: [PATCH 07/15] ClosestPointStrokeHandler ctor bug when result detail is provided --- .../src/curve/internalContexts/ClosestPointStrokeHandler.ts | 2 ++ 1 file changed, 2 insertions(+) diff --git a/core/geometry/src/curve/internalContexts/ClosestPointStrokeHandler.ts b/core/geometry/src/curve/internalContexts/ClosestPointStrokeHandler.ts index b4893c3fbb1..cc9127303ab 100644 --- a/core/geometry/src/curve/internalContexts/ClosestPointStrokeHandler.ts +++ b/core/geometry/src/curve/internalContexts/ClosestPointStrokeHandler.ts @@ -43,6 +43,8 @@ export class ClosestPointStrokeHandler extends NewtonRtoRStrokeHandler implement this._workPoint = Point3d.create(); this._workRay = Ray3d.createZero(); this._closestPoint = result; + if (this._closestPoint) + this._closestPoint.a = Geometry.largeCoordinateResult this._extend = extend ?? false; this.startCurvePrimitive(undefined); this._newtonSolver = new Newton1dUnboundedApproximateDerivative(this); From 33485b95ed6b163d6e42e91b0ae81c13d8ec4223 Mon Sep 17 00:00:00 2001 From: dassaf4 <68340676+dassaf4@users.noreply.github.com> Date: Tue, 25 Nov 2025 14:26:45 -0600 Subject: [PATCH 08/15] remove unused methods and interval logic --- .../CurveCurveCloseApproachXY.ts | 92 ++++--------------- 1 file changed, 19 insertions(+), 73 deletions(-) diff --git a/core/geometry/src/curve/internalContexts/CurveCurveCloseApproachXY.ts b/core/geometry/src/curve/internalContexts/CurveCurveCloseApproachXY.ts index faec86c38ec..eb62bb4f314 100644 --- a/core/geometry/src/curve/internalContexts/CurveCurveCloseApproachXY.ts +++ b/core/geometry/src/curve/internalContexts/CurveCurveCloseApproachXY.ts @@ -124,8 +124,9 @@ export class CurveCurveCloseApproachXY extends RecurseToCurvesGeometryHandler { ): void { const d2 = pointA.distanceSquaredXY(pointB); if (d2 < this._maxDistanceSquared) { - const detailA = CurveLocationDetail.createCurveFractionPoint(cpA, fA, pointA); - const detailB = CurveLocationDetail.createCurveFractionPoint(cpB, fB, pointB); + const d = Math.sqrt(d2); + const detailA = CurveLocationDetail.createCurveFractionPointDistance(cpA, fA, pointA, d); + const detailB = CurveLocationDetail.createCurveFractionPointDistance(cpB, fB, pointB, d); const pair = CurveLocationDetailPair.createCapture(detailA, detailB); if (reversed) pair.swapDetails(); @@ -143,7 +144,6 @@ export class CurveCurveCloseApproachXY extends RecurseToCurvesGeometryHandler { * @param fractionB0 start of the second curve * @param fractionB1 end of the second curve * @param reversed whether to reverse the details in the pair (e.g., so that detailB refers to geometryB). - * @param intervalDetails optional CurveLocationDetailPair */ private recordPointWithLocalFractions( localFractionA: number, @@ -155,25 +155,12 @@ export class CurveCurveCloseApproachXY extends RecurseToCurvesGeometryHandler { fractionB0: number, fractionB1: number, reversed: boolean, - intervalDetails?: undefined | CurveLocationDetailPair, ): void { - let globalFractionA, globalFractionB; - let globalFractionA1, globalFractionB1; - const isInterval = intervalDetails !== undefined && - intervalDetails.detailA.hasFraction1 && - intervalDetails.detailB.hasFraction1; - if (isInterval) { - globalFractionA = Geometry.interpolate(fractionA0, intervalDetails.detailA.fraction, fractionA1); - globalFractionB = Geometry.interpolate(fractionB0, intervalDetails.detailB.fraction, fractionB1); - globalFractionA1 = Geometry.interpolate(fractionA0, intervalDetails.detailA.fraction1!, fractionA1); - globalFractionB1 = Geometry.interpolate(fractionB0, intervalDetails.detailB.fraction1!, fractionB1); - } else { - globalFractionA = globalFractionA1 = Geometry.interpolate(fractionA0, localFractionA, fractionA1); - globalFractionB = globalFractionB1 = Geometry.interpolate(fractionB0, localFractionB, fractionB1); - } + const globalFractionA = Geometry.interpolate(fractionA0, localFractionA, fractionA1); + const globalFractionB = Geometry.interpolate(fractionB0, localFractionB, fractionB1); // ignore duplicate of most recent approach const numPrevious = this._results.length; - if (numPrevious > 0 && !isInterval) { + if (numPrevious > 0) { const oldDetailA = this._results[numPrevious - 1].detailA; const oldDetailB = this._results[numPrevious - 1].detailB; if (reversed) { @@ -186,27 +173,20 @@ export class CurveCurveCloseApproachXY extends RecurseToCurvesGeometryHandler { return; } } - const detailA = CurveLocationDetail.createCurveFractionPoint( - cpA, globalFractionA, cpA.fractionToPoint(globalFractionA), - ); - const detailB = CurveLocationDetail.createCurveFractionPoint( - cpB, globalFractionB, cpB.fractionToPoint(globalFractionB), - ); - if (isInterval) { - detailA.captureFraction1Point1(globalFractionA1, cpA.fractionToPoint(globalFractionA1)); - detailB.captureFraction1Point1(globalFractionB1, cpB.fractionToPoint(globalFractionB1)); - } else { - const d2 = detailA.point.distanceSquaredXY(detailB.point); - if (d2 > this._maxDistanceSquared) - return; - detailA.setIntervalRole(CurveIntervalRole.isolated); - detailB.setIntervalRole(CurveIntervalRole.isolated); - } - if (reversed) { + const pointA = cpA.fractionToPoint(globalFractionA); + const pointB = cpB.fractionToPoint(globalFractionB); + const d2 = pointA.distanceSquaredXY(pointB); + if (d2 > this._maxDistanceSquared) + return; + const d = Math.sqrt(d2); + const detailA = CurveLocationDetail.createCurveFractionPointDistance(cpA, globalFractionA, pointA, d); + const detailB = CurveLocationDetail.createCurveFractionPointDistance(cpB, globalFractionB, pointB, d); + detailA.setIntervalRole(CurveIntervalRole.isolated); + detailB.setIntervalRole(CurveIntervalRole.isolated); + if (reversed) this._results.push(new CurveLocationDetailPair(detailB, detailA)); - } else { + else this._results.push(new CurveLocationDetailPair(detailA, detailB)); - } } /** * Capture a close approach pair that has point and local fraction but not curve. @@ -251,47 +231,13 @@ export class CurveCurveCloseApproachXY extends RecurseToCurvesGeometryHandler { // recompute the points just in case CurveLocationDetail.createCurveEvaluatedFraction(cpA, globalFractionA, pair.detailA); CurveLocationDetail.createCurveEvaluatedFraction(cpB, globalFractionB, pair.detailB); + pair.detailA.a = pair.detailB.a = pair.detailA.point.distanceXY(pair.detailB.point); pair.detailA.setIntervalRole(CurveIntervalRole.isolated); pair.detailB.setIntervalRole(CurveIntervalRole.isolated); if (reversed) pair.swapDetails(); this._results.push(pair); } - /** - * Emit recordPoint for multiple pairs (on full curve) if within maxDistance. - * @param cpA first curve primitive (possibly different from curve in detailA, but fraction compatible) - * @param cpB second curve primitive (possibly different from curve in detailA, but fraction compatible) - * @param pairs array of pairs - * @param reversed whether to reverse the details in the pair (e.g., so that detailB refers to geometryB). - */ - public recordPairs( - cpA: CurvePrimitive, cpB: CurvePrimitive, pairs: CurveLocationDetailPair[] | undefined, reversed: boolean, - ): void { - if (pairs !== undefined) { - for (const p of pairs) { - this.recordPointWithLocalFractions( - p.detailA.fraction, cpA, 0, 1, p.detailB.fraction, cpB, 0, 1, reversed, p, - ); - } - } - } - /** - * Record fully assembled (but possibly reversed) detail pair. - * @param detailA first detail - * @param detailB second detail - * @param reversed whether to reverse the details in the pair (e.g., so that detailB refers to geometryB). - */ - public captureDetailPair( - detailA: CurveLocationDetail | undefined, detailB: CurveLocationDetail | undefined, reversed: boolean, - ): void { - if (detailA && detailB) { - if (reversed) { - this._results.push(CurveLocationDetailPair.createCapture(detailB, detailA)); - } else { - this._results.push(CurveLocationDetailPair.createCapture(detailA, detailB)); - } - } - } private static updatePointToSegmentDistance( fractionA: number, pointA: Point3d, From 6c4e6dc2b8a24210c9a3baeaf0872279502385fc Mon Sep 17 00:00:00 2001 From: dassaf4 <68340676+dassaf4@users.noreply.github.com> Date: Tue, 25 Nov 2025 14:58:46 -0600 Subject: [PATCH 09/15] doc --- core/geometry/src/curve/CurvePrimitive.ts | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/core/geometry/src/curve/CurvePrimitive.ts b/core/geometry/src/curve/CurvePrimitive.ts index ee062c62c00..b06b066fd2a 100644 --- a/core/geometry/src/curve/CurvePrimitive.ts +++ b/core/geometry/src/curve/CurvePrimitive.ts @@ -174,8 +174,9 @@ export abstract class CurvePrimitive extends GeometryQuery { public abstract fractionToPointAndDerivative(fraction: number, result?: Ray3d): Ray3d; /** * Returns a ray whose origin is the curve point and direction is the unit tangent. - * @param fraction fractional position on the curve + * @param fraction fractional position on the curve. * @param result optional preallocated ray. + * @returns tangent ray with normalized direction or zero vector if the derivative vanishes. */ public fractionToPointAndUnitTangent(fraction: number, result?: Ray3d): Ray3d { const ray = this.fractionToPointAndDerivative(fraction, result); From 23e964a14d370d07d7048d781871d0d317f41b1b Mon Sep 17 00:00:00 2001 From: dassaf4 <68340676+dassaf4@users.noreply.github.com> Date: Wed, 26 Nov 2025 09:38:04 -0600 Subject: [PATCH 10/15] doc --- core/geometry/src/numerics/Newton.ts | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/core/geometry/src/numerics/Newton.ts b/core/geometry/src/numerics/Newton.ts index c01b63b965b..145278e0947 100644 --- a/core/geometry/src/numerics/Newton.ts +++ b/core/geometry/src/numerics/Newton.ts @@ -37,7 +37,7 @@ export abstract class AbstractNewtonIterator { public abstract applyCurrentStep(isFinalStep: boolean): boolean; /** * The constructor. - * @param stepSizeTarget tolerance to consider a single step converged. + * @param stepSizeTolerance tolerance to consider a single step converged. * This number should be "moderately" strict. Because 2 successive convergences are required, * it is expected that a first "accept" for (say) 10 to 14 digit step will be followed by another * iteration. A well behaved newton would then hypothetically double the number of digits to From 24ffc64210dfa89e11dbdb46bccad404f3b88aec Mon Sep 17 00:00:00 2001 From: dassaf4 <68340676+dassaf4@users.noreply.github.com> Date: Wed, 26 Nov 2025 10:20:56 -0600 Subject: [PATCH 11/15] finally added a tol arg for comparing points --- core/geometry/src/test/Checker.ts | 24 ++++++++++++++++++------ 1 file changed, 18 insertions(+), 6 deletions(-) diff --git a/core/geometry/src/test/Checker.ts b/core/geometry/src/test/Checker.ts index 3c84489112c..f1ff04c26fa 100644 --- a/core/geometry/src/test/Checker.ts +++ b/core/geometry/src/test/Checker.ts @@ -126,9 +126,13 @@ export class Checker { this._numOK++; return true; } - /** Test if 2 Point3ds are almost equal. */ + /** + * Test if 2 Point3ds are almost equal. + * @param params args to print on error; if first arg is a number, it is used as the comparison tolerance. + */ public testPoint3d(dataA: Point3d, dataB: Point3d, ...params: any[]): boolean { - if (Geometry.isSamePoint3d(dataA, dataB)) + const tol = (params.length > 0 && Number.isFinite(params[0])) ? params.splice(0, 1)[0] : undefined; + if (Geometry.isSamePoint3d(dataA, dataB, tol)) return this.announceOK(); this.announceError("expect same Point3d", dataA, dataB, params); return false; @@ -207,16 +211,24 @@ export class Checker { this.announceError("expect same Range2d", dataA, dataB, params); return false; } - /** Test if 2 Point3ds have almost equal X and Y parts. */ + /** + * Test if 2 Point3ds have almost equal X and Y parts. + * @param params args to print on error; if first arg is a number, it is used as the comparison tolerance. + */ public testPoint3dXY(dataA: Point3d, dataB: Point3d, ...params: any[]): boolean { - if (Geometry.isSamePoint3dXY(dataA, dataB)) + const tol = (params.length > 0 && Number.isFinite(params[0])) ? params.splice(0, 1)[0] : undefined; + if (Geometry.isSamePoint3dXY(dataA, dataB, tol)) return this.announceOK(); this.announceError("expect same Point3d XY", dataA, dataB, params); return false; } - /** Test if 2 Point2ds are almost equal. */ + /** + * Test if 2 Point2ds are almost equal. + * @param params args to print on error; if first arg is a number, it is used as the comparison tolerance. + */ public testPoint2d(dataA: Point2d, dataB: Point2d, ...params: any[]): boolean { - if (Geometry.isSamePoint2d(dataA, dataB)) + const tol = (params.length > 0 && Number.isFinite(params[0])) ? params.splice(0, 1)[0] : undefined; + if (Geometry.isSamePoint2d(dataA, dataB, tol)) return this.announceOK(); this.announceError("expect same Point2d", dataA, dataB, params); return false; From 4d02a6c25db4c7ff5e8aabd140fd098b79e2a1b7 Mon Sep 17 00:00:00 2001 From: dassaf4 <68340676+dassaf4@users.noreply.github.com> Date: Wed, 26 Nov 2025 14:49:41 -0600 Subject: [PATCH 12/15] fix anomalous testPoint3d call --- core/geometry/src/test/clipping/ClipPlanes.test.ts | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/core/geometry/src/test/clipping/ClipPlanes.test.ts b/core/geometry/src/test/clipping/ClipPlanes.test.ts index c02a609007d..298a5066027 100644 --- a/core/geometry/src/test/clipping/ClipPlanes.test.ts +++ b/core/geometry/src/test/clipping/ClipPlanes.test.ts @@ -438,8 +438,8 @@ describe("ClipPlaneSet", () => { convexSet0.reloadSweptPolygon(triangle, sweepDirection, 0); convexSetA.clipUnboundedSegment(linePointA, linePointB, (fA: number, fB: number) => { - ck.testPoint3d(edgePoint01, linePointA.interpolate(fA, linePointB), fA, "unboundedLine clip pointA"); - ck.testPoint3d(edgePoint12, linePointA.interpolate(fB, linePointB), fB, "unboundedLine clip pointB"); + ck.testPoint3d(edgePoint01, linePointA.interpolate(fA, linePointB), `unboundedLine clip pointA at fraction ${fA}`); + ck.testPoint3d(edgePoint12, linePointA.interpolate(fB, linePointB), `unboundedLine clip pointB at fraction ${fB}`); }); const tolerance = 1.0e-10; const values: number[] = [-0.5, 0.3, 0.5, 0.8, 1.1]; From adbd02c63a27efd34b97e5238cc44cfe56fd3163 Mon Sep 17 00:00:00 2001 From: dassaf4 <68340676+dassaf4@users.noreply.github.com> Date: Wed, 26 Nov 2025 16:06:19 -0600 Subject: [PATCH 13/15] spiral intersectXY: - compute two types of Newton seeds from solving the discrete spiral intersectionXY and closeApproachXY problems. - order seeds carefully so that most accurate is processed first. - cull convergent solutions with a dual tolerance filter. Subsequent equivalent solutions are not recorded. - remove check for a solution seed; this can install a less-accurate result. Instead, check the last iterate if Newton fails early on a zero derivative. - compute the max stroke error for a spiral ~ close to the midpoint of first or last stroke segment, whichever spans higher spiral curvature. This is used to limit the results of the discrete spiral closeApproachXY problem. Previous value of 1 micron was way too small, and will fail the new test. --- .../internalContexts/CurveCurveIntersectXY.ts | 170 ++++++++++++------ .../test/curve/CurveCurveIntersectXY.test.ts | 104 ++++++++--- 2 files changed, 201 insertions(+), 73 deletions(-) diff --git a/core/geometry/src/curve/internalContexts/CurveCurveIntersectXY.ts b/core/geometry/src/curve/internalContexts/CurveCurveIntersectXY.ts index 31b98e426d9..c617967d529 100644 --- a/core/geometry/src/curve/internalContexts/CurveCurveIntersectXY.ts +++ b/core/geometry/src/curve/internalContexts/CurveCurveIntersectXY.ts @@ -7,7 +7,7 @@ * @module Curve */ -import { assert } from "@itwin/core-bentley"; +import { assert, DuplicatePolicy, OrderedComparator, SortedArray } from "@itwin/core-bentley"; import { BezierCurve3dH } from "../../bspline/BezierCurve3dH"; import { BezierCurveBase } from "../../bspline/BezierCurveBase"; import { BSplineCurve3d, BSplineCurve3dBase } from "../../bspline/BSplineCurve"; @@ -30,6 +30,7 @@ import { SmallSystem } from "../../numerics/SmallSystem"; import { Arc3d } from "../Arc3d"; import { CurveChainWithDistanceIndex } from "../CurveChainWithDistanceIndex"; import { CurveCollection } from "../CurveCollection"; +import { CurveCurve } from "../CurveCurve"; import { CurveIntervalRole, CurveLocationDetail, CurveLocationDetailPair } from "../CurveLocationDetail"; import { CurvePrimitive } from "../CurvePrimitive"; import { AnyCurve } from "../CurveTypes"; @@ -40,7 +41,7 @@ import { Loop } from "../Loop"; import { Path } from "../Path"; import { ProxyCurve } from "../ProxyCurve"; import { TransitionSpiral3d } from "../spiral/TransitionSpiral3d"; -import { CurveCurve } from "../CurveCurve"; +import { StrokeOptions } from "../StrokeOptions"; // cspell:word XYRR /** @@ -1175,33 +1176,121 @@ export class CurveCurveIntersectXY extends RecurseToCurvesGeometryHandler { * @param reversed Whether `spiralB` data is in `detailA` of each recorded pair, and `curveA` data in `detailB`. */ private refineSpiralResultsByNewton(curveA: CurvePrimitive, spiralB: TransitionSpiral3d, index0: number, reversed = false): void { - const resultsToBeRefined = this._results.slice(index0); - this._results.length -= resultsToBeRefined.length; // keep already refined results - for (const pair of resultsToBeRefined) { + if (index0 >= this._results.length) + return; + // ASSUME: seeds in results tail are ordered by most accurate first, as only the first convergence within tolerance is recorded. + const xyMatchingFunction = new CurveCurveIntersectionXYRRToRRD(curveA, spiralB); + const maxIterations = 100; // observed 73 iterations to convergence in tangent case + const newtonSearcher = new Newton2dUnboundedWithDerivative(xyMatchingFunction, maxIterations); + const fractionTol = 2 * newtonSearcher.stepSizeTolerance; // relative cluster diameter for Newton convergence + const comparePairs: OrderedComparator = (a: CurveLocationDetailPair, b: CurveLocationDetailPair): number => { + assert(() => a.detailA.curve === b.detailA.curve && a.detailB.curve === b.detailB.curve, "pairs are compatible"); + // sort on either fraction, then on the point, using appropriate tolerances for each + if (Geometry.isAlmostEqualNumber(a.detailA.fraction, b.detailA.fraction, fractionTol)) + return 0; + if (a.detailA.point.isAlmostEqualXY(b.detailA.point, this._coincidentGeometryContext.tolerance)) + return 0; + return a.detailA.fraction - b.detailA.fraction; + }; + const myResults = new SortedArray(comparePairs, DuplicatePolicy.Retain); + const pushToMyResults = (cpA: CurvePrimitive, fA: number, cpB: CurvePrimitive, fB: number): boolean => { + const detailA = CurveLocationDetail.createCurveFractionPoint(cpA, fA, cpA.fractionToPoint(fA)); + const detailB = CurveLocationDetail.createCurveFractionPoint(cpB, fB, cpB.fractionToPoint(fB)); + detailA.setIntervalRole(CurveIntervalRole.isolated); + detailB.setIntervalRole(CurveIntervalRole.isolated); + let pushed = false; + myResults.insert(new CurveLocationDetailPair(reversed ? detailB : detailA, reversed ? detailA : detailB), () => pushed = true); + return pushed; + }; + for (let i = index0; i < this._results.length; i++) { + const pair = this._results[i]; const detailA = reversed ? pair.detailB : pair.detailA; const detailB = reversed ? pair.detailA : pair.detailB; assert(detailB.curve instanceof LineString3d, "Caller has discretized the spiral"); - let fractionA = detailA.fraction; - let fractionB = detailB.fraction; // use linestring fraction as the spiral fraction; it is closer than fractional length! - // Newton iteration fails when the (partial) derivative is small. Check for a root there to avoid missing it. - if (spiralB.fractionToPoint(fractionB).isAlmostEqualXY(detailA.point, this._coincidentGeometryContext.tolerance)) { - this.recordPointWithLocalFractions(fractionA, curveA, 0, 1, fractionB, spiralB, 0, 1, reversed); - continue; - } const extendA0 = reversed ? this._extendB0 : this._extendA0; const extendA1 = reversed ? this._extendB1 : this._extendA1; - const xyMatchingFunction = new CurveCurveIntersectionXYRRToRRD(curveA, spiralB); - const maxIterations = 64; // observed: 55 iterations for tangent intersection (linear convergence) - const newtonSearcher = new Newton2dUnboundedWithDerivative(xyMatchingFunction, maxIterations); - const fractionTol = 2 * newtonSearcher.stepSizeTolerance; // relative cluster diameter for Newton convergence - newtonSearcher.setUV(fractionA, fractionB); + newtonSearcher.setUV(detailA.fraction, detailB.fraction); // use linestring fraction as spiral param; it generally yields a closer point than fractional length! if (newtonSearcher.runIterations()) { - fractionA = newtonSearcher.getU(); - fractionB = newtonSearcher.getV(); + const fractionA = newtonSearcher.getU(); + const fractionB = newtonSearcher.getV(); if (this.acceptFraction(extendA0, fractionA, extendA1) && this.acceptFraction(false, fractionB, false)) - this.recordPointWithLocalFractions(fractionA, curveA, 0, 1, fractionB, spiralB, 0, 1, reversed, undefined, fractionTol); + pushToMyResults(curveA, fractionA, spiralB, fractionB); + } else if (newtonSearcher.numIterations < 10) { + // if Newton failed early due to vanishing (partial) derivative, check for a root there + const fractionA = newtonSearcher.getU(); + const fractionB = newtonSearcher.getV(); + if (curveA.fractionToPoint(fractionA).isAlmostEqualXY(spiralB.fractionToPoint(fractionB), this._coincidentGeometryContext.tolerance)) + pushToMyResults(curveA, fractionA, spiralB, fractionB); + } + } + this._results.splice(index0, this._results.length - index0, ...myResults.extractArray()); + } + /** + * Append stroke points and return the line string. + * * This is a convenient wrapper for [[CurvePrimitive.emitStrokes]] but the analogous instance method cannot be added + * to that class due to the ensuing recursion with subclass [[LineString3d]]. + * @param options options for stroking the instance curve. + * @param result object to receive appended stroke points; if omitted, a new object is created, populated, and returned. + */ + private strokeCurve(curve: CurvePrimitive, options?: StrokeOptions, result?: LineString3d): LineString3d { + const ls = result ? result : LineString3d.create(); + curve.emitStrokes(ls, options); + return ls; + } + /** Compute an approximation to the max chord height error of the stroked spiral, and another curve if a spiral. */ + private computeMaxSpiralStrokeError(spiral0: TransitionSpiral3d, ls0: LineString3d, spiral1?: CurvePrimitive, ls1?: CurvePrimitive): number { + let maxError = 0; + if (ls0.numEdges() > 0) { + // the max error occurs at the spiral end with higher curvature + const k0 = spiral0.fractionToCurvature(0); + const k1 = spiral0.fractionToCurvature(1); + const iChord = (k0 !== undefined && k1 !== undefined && Math.abs(k0) > Math.abs(k1)) ? 0 : ls0.numEdges() - 1; + const midPoint = Point3d.create(); + const detail = CurveLocationDetail.create(); + if (ls0.packedPoints.interpolateIndexIndex(iChord, 0.5, iChord + 1, midPoint)) { + if (spiral0.closestPoint(midPoint, false, detail)) + maxError = detail.a; + } + if (spiral1 && spiral1 instanceof TransitionSpiral3d && ls1 && ls1 instanceof LineString3d) { + const maxError1 = this.computeMaxSpiralStrokeError(spiral1, ls1); + if (maxError1 > maxError) + maxError = maxError1; } } + return maxError; + } + /** + * Solve the intersection problem for stroked, unextended curveB. + * * @return the number of results appended. + */ + private appendDiscreteIntersectionResults(curveA: CurvePrimitive, extendA0: boolean, extendA1: boolean, lsB: LineString3d, reversed: boolean): number { + const i0 = this._results.length; + // handleLineString3d requires us to swap geometries: + const geomB = this._geometryB; + const extendB0 = this._extendB0; + const extendB1 = this._extendB1; + this.resetGeometryA(false, false); // lsB is never extended + this.resetGeometryB(curveA, extendA0, extendA1); + this.handleLineString3d(lsB); // this puts lsB data in detailA, as expected when reversed is true + if (!reversed) { // swap lsB data to detailB + for (let i = i0; i < this._results.length; i++) + this._results[i].swapDetails(); + } + this.resetGeometryA(extendA0, extendA1); + this.resetGeometryB(geomB, extendB0, extendB1); + return this._results.length - i0; + } + /** + * Solve the close approach problem for stroked, unextended curveB. + * * Sort the results shortest projection distance first. + * @return the number of results appended. + */ + private appendDiscreteCloseApproachResults(curveA: CurvePrimitive, lsB: LineString3d, maxDistance: number, reversed: boolean): number { + const i0 = this._results.length; + const closeApproachPairs = CurveCurve.closeApproachProjectedXYPairs(reversed ? lsB : curveA, reversed ? curveA : lsB, maxDistance); + closeApproachPairs.sort((p0: CurveLocationDetailPair, p1: CurveLocationDetailPair) => p0.detailA.a - p1.detailA.a); + this._results.push(...closeApproachPairs); + return this._results.length - i0; } /** * Compute the xy-intersection of a curve and a spiral. @@ -1210,39 +1299,20 @@ export class CurveCurveIntersectXY extends RecurseToCurvesGeometryHandler { * @param extendA0 whether to compute xy-intersections with curveA extended beyond its start. * @param extendA1 whether to compute xy-intersections with curveA extended beyond its end. * @param spiralB transition spiral to intersect with curveA. - * @param reversed Whether `spiralB` data is in `detailA` of each recorded pair, and `curveA` data in `detailB`. + * @param reversed whether `spiralB` data will be recorded in `detailA` of each result, and `curveA` data in `detailB`. */ private dispatchCurveSpiral(curveA: CurvePrimitive, extendA0: boolean, extendA1: boolean, spiralB: TransitionSpiral3d, reversed: boolean): void { - const index0 = this._results.length; - // approximate spiral(s) - const cpA = curveA instanceof TransitionSpiral3d ? LineString3d.create() : curveA; - const cpB = LineString3d.create(); - spiralB.emitStrokes(cpB); - let spiralLength = spiralB.quickLength(); - if (curveA instanceof TransitionSpiral3d && cpA instanceof LineString3d) { - curveA.emitStrokes(cpA); - spiralLength += curveA.quickLength(); - } - // append seeds using unextended curveA and stroked spiral - const maxDistance = Geometry.smallMetricDistance * (1 + 100 * spiralLength); // HEURISTIC: not too tight - this._results.push(...CurveCurve.closeApproachProjectedXYPairs(reversed ? cpB : cpA, reversed ? cpA : cpB, maxDistance)); - // append seeds using extended curveA and stroked spiral - if (curveA.isExtensibleFractionSpace && (extendA0 || extendA1) && this._geometryB) { - const geomB = this._geometryB; - const extendB0 = this._extendB0; - const extendB1 = this._extendB1; - const i0 = this._results.length; - // handleLineString3d requires us to swap geometries: - this.resetGeometryA(false, false); // spiral approximation cpB is sent in as geometryA - this.resetGeometryB(curveA, extendA0, extendA1); // curveA is installed as geometryB - this.handleLineString3d(cpB); // this puts spiral data in detailA (expected when reversed is true) - if (!reversed) { - for (let i = i0; i < this._results.length; i++) - this._results[i].swapDetails(); // now spiral data is in detailB - } - this.resetGeometryA(extendA0, extendA1); - this.resetGeometryB(geomB, extendB0, extendB1); + let cpA = curveA; + if (curveA instanceof TransitionSpiral3d) { + cpA = this.strokeCurve(curveA); + extendA0 = extendA1 = false; } + const cpB = this.strokeCurve(spiralB); + const maxError = 1.01 * this.computeMaxSpiralStrokeError(spiralB, cpB, curveA, cpA); + const index0 = this._results.length; + // append seeds computed by solving discretized spiral problems, then refine the seeds via Newton + this.appendDiscreteIntersectionResults(cpA, extendA0, extendA1, cpB, reversed); // recorded first because more accurate + this.appendDiscreteCloseApproachResults(cpA, cpB, maxError, reversed); // seeds for finding tangent intersections this.refineSpiralResultsByNewton(curveA, spiralB, index0, reversed); } /** Double dispatch handler for strongly typed spiral curve. */ diff --git a/core/geometry/src/test/curve/CurveCurveIntersectXY.test.ts b/core/geometry/src/test/curve/CurveCurveIntersectXY.test.ts index 071b981b519..3c19a950966 100644 --- a/core/geometry/src/test/curve/CurveCurveIntersectXY.test.ts +++ b/core/geometry/src/test/curve/CurveCurveIntersectXY.test.ts @@ -8,7 +8,7 @@ import { Arc3d } from "../../curve/Arc3d"; import { CurveChainWithDistanceIndex } from "../../curve/CurveChainWithDistanceIndex"; import { BagOfCurves, CurveCollection } from "../../curve/CurveCollection"; import { CurveCurve } from "../../curve/CurveCurve"; -import { CurveLocationDetailPair } from "../../curve/CurveLocationDetail"; +import { CurveLocationDetail, CurveLocationDetailPair } from "../../curve/CurveLocationDetail"; import { CurvePrimitive } from "../../curve/CurvePrimitive"; import { GeometryQuery } from "../../curve/GeometryQuery"; import { LineSegment3d } from "../../curve/LineSegment3d"; @@ -2373,20 +2373,17 @@ describe("CurveCurveIntersectXY", () => { ck: Checker, allGeometry: GeometryQuery[], curve0: AnyCurve, curve1: AnyCurve, numExpected: number, dx: number, dy: number, - extend0 = false, extend1 = false, + extend0: boolean = false, extend1: boolean = false, + pointTol?: number, isTangent: boolean = false, ) { GeometryCoreTestIO.captureCloneGeometry(allGeometry, curve0, dx, dy); GeometryCoreTestIO.captureCloneGeometry(allGeometry, curve1, dx, dy); const testSpiralIntersection = (intersections: CurveLocationDetailPair[], reversed: boolean) => { GeometryCoreTestIO.captureCurveLocationDetails(allGeometry, intersections, 5, dx, dy); - const curveName0 = curve0.constructor.name; - const curveName1 = curve1.constructor.name; - ck.testExactNumber( - numExpected, - intersections.length, - `expect ${numExpected} intersection(s) between ${curveName1} and ${curveName0}`, - ); + const curveName0 = curve0 instanceof TransitionSpiral3d ? curve0.spiralType : curve0.constructor.name; + const curveName1 = curve1 instanceof TransitionSpiral3d ? curve1.spiralType : curve1.constructor.name; + ck.testExactNumber(numExpected, intersections.length, `expect ${numExpected} ${isTangent ? "tangent" : ""} intersection${numExpected === 1 ? "" : "(s)"} between ${curveName0} and ${curveName1}`) for (const intersection of intersections) { let curve0Detail = intersection.detailA; let curve1Detail = intersection.detailB; @@ -2428,12 +2425,33 @@ describe("CurveCurveIntersectXY", () => { } } // test both paths - const intersectionsAB = CurveCurve.intersectionXYPairs(curve0, extend0, curve1, extend1); + const intersectionsAB = CurveCurve.intersectionXYPairs(curve0, extend0, curve1, extend1, pointTol); testSpiralIntersection(intersectionsAB, false); - const intersectionsBA = CurveCurve.intersectionXYPairs(curve1, extend1, curve0, extend0); + const intersectionsBA = CurveCurve.intersectionXYPairs(curve1, extend1, curve0, extend0, pointTol); testSpiralIntersection(intersectionsBA, true); }; + /** Stroke the spiral, compute midpoint of highest curvature edge, project it to the spiral, return the fraction. */ + function computeSpiralFractionAtMaxStrokeError(spiral: TransitionSpiral3d, strokes?: LineString3d, options?: StrokeOptions): number | undefined { + const ls = strokes ? strokes : LineString3d.create(); + if (ls.numEdges() < 1) + spiral.emitStrokes(ls, options); + if (ls.numEdges() < 1) + return undefined; + const k0 = spiral.fractionToCurvature(0); + const k1 = spiral.fractionToCurvature(1); + if (k0 === undefined || k1 === undefined) + return undefined; + const iChord = (Math.abs(k0) > Math.abs(k1)) ? 0 : ls.numEdges() - 1; + const midPoint = Point3d.create(); + const detail = CurveLocationDetail.create(); + if (!ls.packedPoints.interpolateIndexIndex(iChord, 0.5, iChord + 1, midPoint)) + return undefined; + if (!spiral.closestPoint(midPoint, false, detail)) + return undefined; + return detail.fraction; + } + it("SpiralIntersection", () => { const ck = new Checker(); const allGeometry: GeometryQuery[] = []; @@ -2592,8 +2610,9 @@ describe("CurveCurveIntersectXY", () => { 5, 6, 5, 10, // path, loop, curve chain, and bag of curves ]; ck.testExactNumber(curves.length, numExpectedIntersections.length, "matching arrays"); - // spiral vs all curves - const test0 = (spiral: TransitionSpiral3d, ddy = 0, extend = false) => { + + // spiral vs curve + const spiralIntersectCurveTest = (spiral: TransitionSpiral3d, ddy = 0, extend = false) => { for (let j = 0; j < curves.length; j++) { const curve = curves[j]; const numExpectedIntersection = numExpectedIntersections[j]; @@ -2605,15 +2624,16 @@ describe("CurveCurveIntersectXY", () => { } for (let i = 0; i < integratedSpirals.length; i++) // skip rotated and non-planar integrated spirals if (i % 3 === 0) - test0(integratedSpirals[i]); + spiralIntersectCurveTest(integratedSpirals[i]); dx += 250; for (let i = 0; i < directSpirals.length; i++) // skip rotated and non-planar direct spirals if (i % 3 === 0) - test0(directSpirals[i]); + spiralIntersectCurveTest(directSpirals[i]); + + // curve chain/collection vs curve chain/collection dx = 0; dy = 6400; let numExpected = 12; - // curve chain/collection vs curve chain/collection visualizeAndTestSpiralIntersection(ck, allGeometry, curveChain0, curveChain1, numExpected, dx, dy); dy += 200; visualizeAndTestSpiralIntersection(ck, allGeometry, path0, path1, numExpected, dx, dy); @@ -2621,24 +2641,34 @@ describe("CurveCurveIntersectXY", () => { visualizeAndTestSpiralIntersection(ck, allGeometry, curveChain0, path1, numExpected, dx, dy); dy += 200; visualizeAndTestSpiralIntersection(ck, allGeometry, curveChain1, path0, numExpected, dx, dy); + // tangency at the interior of the spiral dy += 200; numExpected = 1; - const test1 = (spiral: TransitionSpiral3d) => { + const spiralTangencyTest = (spiral: TransitionSpiral3d) => { const ray = spiral.fractionToPointAndDerivative(0.5); - const ls = LineString3d.create( + const tangentLine = LineSegment3d.create( ray.origin.plusScaled(ray.direction.normalize()!, 50), ray.origin.plusScaled(ray.direction.normalize()!, -50) ); - visualizeAndTestSpiralIntersection(ck, allGeometry, spiral, ls, numExpected, dx, dy); + let tol: number | undefined; + // Note 1: This spiral has the smallest curvature (0.007) at tangency of the spirals in the tangency tests. + // The Newton solver really has two tolerances, a fraction tol (for convergence) and a point tol (for solution + // de-duping). These tolerances are uncoupled, so you can end up with situations like this one, where the spiral + // is so flat that Newton fractions converge faster than their corresponding points. Here we loosen the point + // tolerance to cull more duplicate solutions. See also Note 2 below. + if (spiral.spiralType === "AustralianRailCorp") + tol = 10 * Geometry.smallMetricDistance; + visualizeAndTestSpiralIntersection(ck, allGeometry, spiral, tangentLine, numExpected, dx, dy, undefined, undefined, tol, true); dx += 200; } for (let i = 0; i < integratedSpirals.length; i++) // skip rotated and non-planar integrated spirals if (i % 3 === 0) - test1(integratedSpirals[i]); + spiralTangencyTest(integratedSpirals[i]); dx += 250; for (let i = 0; i < directSpirals.length; i++) // skip rotated and non-planar direct spirals if (i % 3 === 0) - test1(directSpirals[i]); + spiralTangencyTest(directSpirals[i]); + // extend curve primitive dx = 0; dy += 200; @@ -2658,11 +2688,39 @@ describe("CurveCurveIntersectXY", () => { ck.testExactNumber(curves.length, numExpectedIntersections.length, "matching arrays"); for (let i = 0; i < integratedSpirals.length; i++) // skip rotated and non-planar integrated spirals if (i % 3 === 0) - test0(integratedSpirals[i], 7400, true); + spiralIntersectCurveTest(integratedSpirals[i], 7400, true); dx += 250; for (let i = 0; i < directSpirals.length; i++) // skip rotated and non-planar direct spirals if (i % 3 === 0) - test0(directSpirals[i], 7400, true); + spiralIntersectCurveTest(directSpirals[i], 7400, true); + + // Worst-case tangent test: justifies maxDistance used in computing closeApproachXY seeds for Newton. + dx = 0; + dy = -300; + const allSpirals = integratedSpirals.concat(directSpirals); + for (let i = 0; i < allSpirals.length; i++) { + const spiral = allSpirals[i]; + const fractionAtTangent = computeSpiralFractionAtMaxStrokeError(spiral); + if (ck.testDefined(fractionAtTangent, "found spiral fraction at max stroke error")) { + const ray = spiral.fractionToPointAndUnitTangent(fractionAtTangent); + ck.testTrue(!ray.direction.isZero, "max error point found"); + const line = LineSegment3d.create(ray.origin.plusScaled(ray.direction, 50), ray.origin.plusScaled(ray.direction, -50)); + // Note 2: On a tangency intersection, Newton's method suffers from linear convergence, and hence iterates + // can converge faster to each other than to a known intersection. To account for this we goose the tolerances. + // See also Note 1 above. + const pointTol = 10 * Geometry.smallMetricDistance; + const fractionTol = 1000 * Geometry.smallFraction; + const intersections = CurveCurve.intersectionXYPairs(spiral, false, line, false, pointTol); + GeometryCoreTestIO.captureCloneGeometry(allGeometry, [spiral, line], dx, dy); + GeometryCoreTestIO.captureCurveLocationDetails(allGeometry, intersections, 5, dx, dy); + if (ck.testExactNumber(1, intersections.length, `${spiral.spiralType}[${i}]: expect 1 tangent intersection`)) { + ck.testPoint3dXY(intersections[0].detailA.point, intersections[0].detailB.point, `${spiral.spiralType}[${i}]: spiral and line detail points match`); + ck.testPoint3dXY(ray.origin, intersections[0].detailA.point, pointTol, `${spiral.spiralType}[${i}]: expected spiral point at tangent intersection`); + ck.testNearNumber(fractionAtTangent, intersections[0].detailA.fraction, fractionTol, `${spiral.spiralType}[${i}]: expected spiral fraction at tangent intersection`); + } + dx += 300; + } + } GeometryCoreTestIO.saveGeometry(allGeometry, "CurveCurveIntersectXY", "SpiralIntersection"); expect(ck.getNumErrors()).toBe(0); From 05de9d5f60cb9783b5e7de331c4429bd113fba83 Mon Sep 17 00:00:00 2001 From: dassaf4 <68340676+dassaf4@users.noreply.github.com> Date: Wed, 26 Nov 2025 16:21:21 -0600 Subject: [PATCH 14/15] extract-api --- common/api/core-geometry.api.md | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/common/api/core-geometry.api.md b/common/api/core-geometry.api.md index 9b0ea5a3e16..7fa1ff884ab 100644 --- a/common/api/core-geometry.api.md +++ b/common/api/core-geometry.api.md @@ -1711,7 +1711,7 @@ export class CurveLocationDetail { isSameCurveAndFraction(other: CurveLocationDetail | { curve: CurvePrimitive; fraction: number; - }): boolean; + }, fractionTol?: number): boolean; point: Point3d; point1?: Point3d; pointQ: Point3d; @@ -2224,6 +2224,7 @@ export class Geometry { static isSameCoordinateSquared(x: number, y: number, tolerance?: number): boolean; static isSameCoordinateWithToleranceFactor(x: number, y: number, toleranceFactor: number): boolean; static isSameCoordinateXY(x0: number, y0: number, x1: number, y1: number, tolerance?: number): boolean; + static isSameFraction(x: number, y: number, tolerance?: number): boolean; static isSamePoint2d(dataA: Point2d, dataB: Point2d, tolerance?: number): boolean; static isSamePoint3d(dataA: Point3d, dataB: Point3d, tolerance?: number): boolean; static isSamePoint3dXY(dataA: Point3d, dataB: Point3d, tolerance?: number): boolean; From 046dd223796fde0822680f3bfe89f29581a22f78 Mon Sep 17 00:00:00 2001 From: "Saeed.Torabi" <114527638+saeeedtorabi@users.noreply.github.com> Date: Thu, 27 Nov 2025 14:15:10 -0500 Subject: [PATCH 15/15] Broke long lines to shorter lines to increase readability --- .../internalContexts/CurveCurveIntersectXY.ts | 35 ++++++++++++++----- .../test/curve/CurveCurveIntersectXY.test.ts | 31 ++++++++++++---- 2 files changed, 51 insertions(+), 15 deletions(-) diff --git a/core/geometry/src/curve/internalContexts/CurveCurveIntersectXY.ts b/core/geometry/src/curve/internalContexts/CurveCurveIntersectXY.ts index c617967d529..38e8400cb26 100644 --- a/core/geometry/src/curve/internalContexts/CurveCurveIntersectXY.ts +++ b/core/geometry/src/curve/internalContexts/CurveCurveIntersectXY.ts @@ -1169,13 +1169,16 @@ export class CurveCurveIntersectXY extends RecurseToCurvesGeometryHandler { } /** * Process tail of `this._results` for xy-intersections between the curve and spiral. - * * If a result is not already an intersection, refine it via Newton iteration unless it doesn't converge, in which case remove it. + * * If a result is not already an intersection, refine it via Newton iteration unless it doesn't converge, in which + * case remove it. * @param curveA The other curve primitive. May also be a transition spiral. * @param spiralB The transition spiral. * @param index0 index of first entry in tail of `this._results` to refine. * @param reversed Whether `spiralB` data is in `detailA` of each recorded pair, and `curveA` data in `detailB`. */ - private refineSpiralResultsByNewton(curveA: CurvePrimitive, spiralB: TransitionSpiral3d, index0: number, reversed = false): void { + private refineSpiralResultsByNewton( + curveA: CurvePrimitive, spiralB: TransitionSpiral3d, index0: number, reversed = false, + ): void { if (index0 >= this._results.length) return; // ASSUME: seeds in results tail are ordered by most accurate first, as only the first convergence within tolerance is recorded. @@ -1183,7 +1186,9 @@ export class CurveCurveIntersectXY extends RecurseToCurvesGeometryHandler { const maxIterations = 100; // observed 73 iterations to convergence in tangent case const newtonSearcher = new Newton2dUnboundedWithDerivative(xyMatchingFunction, maxIterations); const fractionTol = 2 * newtonSearcher.stepSizeTolerance; // relative cluster diameter for Newton convergence - const comparePairs: OrderedComparator = (a: CurveLocationDetailPair, b: CurveLocationDetailPair): number => { + const comparePairs: OrderedComparator = ( + a: CurveLocationDetailPair, b: CurveLocationDetailPair, + ): number => { assert(() => a.detailA.curve === b.detailA.curve && a.detailB.curve === b.detailB.curve, "pairs are compatible"); // sort on either fraction, then on the point, using appropriate tolerances for each if (Geometry.isAlmostEqualNumber(a.detailA.fraction, b.detailA.fraction, fractionTol)) @@ -1199,7 +1204,9 @@ export class CurveCurveIntersectXY extends RecurseToCurvesGeometryHandler { detailA.setIntervalRole(CurveIntervalRole.isolated); detailB.setIntervalRole(CurveIntervalRole.isolated); let pushed = false; - myResults.insert(new CurveLocationDetailPair(reversed ? detailB : detailA, reversed ? detailA : detailB), () => pushed = true); + myResults.insert( + new CurveLocationDetailPair(reversed ? detailB : detailA, reversed ? detailA : detailB), () => pushed = true, + ); return pushed; }; for (let i = index0; i < this._results.length; i++) { @@ -1238,7 +1245,9 @@ export class CurveCurveIntersectXY extends RecurseToCurvesGeometryHandler { return ls; } /** Compute an approximation to the max chord height error of the stroked spiral, and another curve if a spiral. */ - private computeMaxSpiralStrokeError(spiral0: TransitionSpiral3d, ls0: LineString3d, spiral1?: CurvePrimitive, ls1?: CurvePrimitive): number { + private computeMaxSpiralStrokeError( + spiral0: TransitionSpiral3d, ls0: LineString3d, spiral1?: CurvePrimitive, ls1?: CurvePrimitive, + ): number { let maxError = 0; if (ls0.numEdges() > 0) { // the max error occurs at the spiral end with higher curvature @@ -1263,7 +1272,9 @@ export class CurveCurveIntersectXY extends RecurseToCurvesGeometryHandler { * Solve the intersection problem for stroked, unextended curveB. * * @return the number of results appended. */ - private appendDiscreteIntersectionResults(curveA: CurvePrimitive, extendA0: boolean, extendA1: boolean, lsB: LineString3d, reversed: boolean): number { + private appendDiscreteIntersectionResults( + curveA: CurvePrimitive, extendA0: boolean, extendA1: boolean, lsB: LineString3d, reversed: boolean, + ): number { const i0 = this._results.length; // handleLineString3d requires us to swap geometries: const geomB = this._geometryB; @@ -1285,9 +1296,13 @@ export class CurveCurveIntersectXY extends RecurseToCurvesGeometryHandler { * * Sort the results shortest projection distance first. * @return the number of results appended. */ - private appendDiscreteCloseApproachResults(curveA: CurvePrimitive, lsB: LineString3d, maxDistance: number, reversed: boolean): number { + private appendDiscreteCloseApproachResults( + curveA: CurvePrimitive, lsB: LineString3d, maxDistance: number, reversed: boolean, + ): number { const i0 = this._results.length; - const closeApproachPairs = CurveCurve.closeApproachProjectedXYPairs(reversed ? lsB : curveA, reversed ? curveA : lsB, maxDistance); + const closeApproachPairs = CurveCurve.closeApproachProjectedXYPairs( + reversed ? lsB : curveA, reversed ? curveA : lsB, maxDistance, + ); closeApproachPairs.sort((p0: CurveLocationDetailPair, p1: CurveLocationDetailPair) => p0.detailA.a - p1.detailA.a); this._results.push(...closeApproachPairs); return this._results.length - i0; @@ -1301,7 +1316,9 @@ export class CurveCurveIntersectXY extends RecurseToCurvesGeometryHandler { * @param spiralB transition spiral to intersect with curveA. * @param reversed whether `spiralB` data will be recorded in `detailA` of each result, and `curveA` data in `detailB`. */ - private dispatchCurveSpiral(curveA: CurvePrimitive, extendA0: boolean, extendA1: boolean, spiralB: TransitionSpiral3d, reversed: boolean): void { + private dispatchCurveSpiral( + curveA: CurvePrimitive, extendA0: boolean, extendA1: boolean, spiralB: TransitionSpiral3d, reversed: boolean, + ): void { let cpA = curveA; if (curveA instanceof TransitionSpiral3d) { cpA = this.strokeCurve(curveA); diff --git a/core/geometry/src/test/curve/CurveCurveIntersectXY.test.ts b/core/geometry/src/test/curve/CurveCurveIntersectXY.test.ts index 3c19a950966..155bc077725 100644 --- a/core/geometry/src/test/curve/CurveCurveIntersectXY.test.ts +++ b/core/geometry/src/test/curve/CurveCurveIntersectXY.test.ts @@ -2658,7 +2658,9 @@ describe("CurveCurveIntersectXY", () => { // tolerance to cull more duplicate solutions. See also Note 2 below. if (spiral.spiralType === "AustralianRailCorp") tol = 10 * Geometry.smallMetricDistance; - visualizeAndTestSpiralIntersection(ck, allGeometry, spiral, tangentLine, numExpected, dx, dy, undefined, undefined, tol, true); + visualizeAndTestSpiralIntersection( + ck, allGeometry, spiral, tangentLine, numExpected, dx, dy, undefined, undefined, tol, true, + ); dx += 200; } for (let i = 0; i < integratedSpirals.length; i++) // skip rotated and non-planar integrated spirals @@ -2675,7 +2677,8 @@ describe("CurveCurveIntersectXY", () => { const lineSegment5 = LineSegment3d.create(Point3d.create(50, 100), Point3d.create(60, 50)); const arc4 = Arc3d.createXY(Point3d.create(50, 10), 30, AngleSweep.createStartEndDegrees(90, 180)); const lineString4 = LineString3d.create( - Point3d.create(20, 20), Point3d.create(20, 40), Point3d.create(40, 40), Point3d.create(40, 20), Point3d.create(70, 20), + Point3d.create(20, 20), Point3d.create(20, 40), Point3d.create(40, 40), + Point3d.create(40, 20), Point3d.create(70, 20), ); const lineString5 = LineString3d.create( Point3d.create(10, 40), Point3d.create(10, 80), Point3d.create(30, 80), @@ -2704,7 +2707,9 @@ describe("CurveCurveIntersectXY", () => { if (ck.testDefined(fractionAtTangent, "found spiral fraction at max stroke error")) { const ray = spiral.fractionToPointAndUnitTangent(fractionAtTangent); ck.testTrue(!ray.direction.isZero, "max error point found"); - const line = LineSegment3d.create(ray.origin.plusScaled(ray.direction, 50), ray.origin.plusScaled(ray.direction, -50)); + const line = LineSegment3d.create( + ray.origin.plusScaled(ray.direction, 50), ray.origin.plusScaled(ray.direction, -50), + ); // Note 2: On a tangency intersection, Newton's method suffers from linear convergence, and hence iterates // can converge faster to each other than to a known intersection. To account for this we goose the tolerances. // See also Note 1 above. @@ -2714,9 +2719,23 @@ describe("CurveCurveIntersectXY", () => { GeometryCoreTestIO.captureCloneGeometry(allGeometry, [spiral, line], dx, dy); GeometryCoreTestIO.captureCurveLocationDetails(allGeometry, intersections, 5, dx, dy); if (ck.testExactNumber(1, intersections.length, `${spiral.spiralType}[${i}]: expect 1 tangent intersection`)) { - ck.testPoint3dXY(intersections[0].detailA.point, intersections[0].detailB.point, `${spiral.spiralType}[${i}]: spiral and line detail points match`); - ck.testPoint3dXY(ray.origin, intersections[0].detailA.point, pointTol, `${spiral.spiralType}[${i}]: expected spiral point at tangent intersection`); - ck.testNearNumber(fractionAtTangent, intersections[0].detailA.fraction, fractionTol, `${spiral.spiralType}[${i}]: expected spiral fraction at tangent intersection`); + ck.testPoint3dXY( + intersections[0].detailA.point, + intersections[0].detailB.point, + `${spiral.spiralType}[${i}]: spiral and line detail points match`, + ); + ck.testPoint3dXY( + ray.origin, + intersections[0].detailA.point, + pointTol, + `${spiral.spiralType}[${i}]: expected spiral point at tangent intersection`, + ); + ck.testNearNumber( + fractionAtTangent, + intersections[0].detailA.fraction, + fractionTol, + `${spiral.spiralType}[${i}]: expected spiral fraction at tangent intersection`, + ); } dx += 300; }