diff --git a/common/api/core-geometry.api.md b/common/api/core-geometry.api.md index f2d8b7804cbf..25dbcc1f1e8c 100644 --- a/common/api/core-geometry.api.md +++ b/common/api/core-geometry.api.md @@ -156,8 +156,8 @@ export class AngleSweep implements BeJSONFunctions { static fromJSON(json?: AngleSweepProps): AngleSweep; interpolate(fraction: number, other: AngleSweep): AngleSweep; isAlmostEqual(other: AngleSweep): boolean; - isAlmostEqualAllowPeriodShift(other: AngleSweep): boolean; - isAlmostEqualNoPeriodShift(other: AngleSweep): boolean; + isAlmostEqualAllowPeriodShift(other: AngleSweep, radianTol?: number): boolean; + isAlmostEqualNoPeriodShift(other: AngleSweep, radianTol?: number): boolean; isAngleInSweep(angle: Angle): boolean; get isCCW(): boolean; get isEmpty(): boolean; @@ -257,10 +257,11 @@ export class Arc3d extends CurvePrimitive implements BeJSONFunctions { static createCircularStartEndRadius(start: Point3d, end: Point3d, radius: number, helper: Point3d | Vector3d): Arc3d | undefined; static createCircularStartMiddleEnd(pointA: XYAndZ, pointB: XYAndZ, pointC: XYAndZ, result?: Arc3d): Arc3d | LineString3d; static createCircularStartTangentEnd(start: Point3d, tangentAtStart: Vector3d, end: Point3d, result?: Arc3d): Arc3d | LineSegment3d; + static createCircularStartTangentRadius(start: Point3d, tangentAtStart: Vector3d, radius: number, upVector?: Vector3d, sweep?: Angle | AngleSweep): Arc3d | undefined; static createFilletArc(point0: Point3d, point1: Point3d, point2: Point3d, radius: number): ArcBlendData; static createRefs(center: Point3d, matrix: Matrix3d, sweep: AngleSweep, result?: Arc3d): Arc3d; static createScaledXYColumns(center: Point3d | undefined, matrix: Matrix3d, radius0: number, radius90: number, sweep?: AngleSweep, result?: Arc3d): Arc3d; - static createStartMiddleEnd(point0: XYAndZ, point1: XYAndZ, point2: XYAndZ, sweep?: AngleSweep, result?: Arc3d): Arc3d | undefined; + static createStartMiddleEnd(start: XYAndZ, middle: XYAndZ, end: XYAndZ, sweep?: AngleSweep, result?: Arc3d): Arc3d | undefined; static createUnitCircle(): Arc3d; static createXY(center: Point3d, radius: number, sweep?: AngleSweep): Arc3d; static createXYEllipse(center: Point3d, radiusA: number, radiusB: number, sweep?: AngleSweep): Arc3d; @@ -281,7 +282,7 @@ export class Arc3d extends CurvePrimitive implements BeJSONFunctions { getFractionToDistanceScale(): number | undefined; // @internal getPlaneAltitudeSineCosinePolynomial(plane: PlaneAltitudeEvaluator, result?: SineCosinePolynomial): SineCosinePolynomial; - isAlmostEqual(otherGeometry: GeometryQuery): boolean; + isAlmostEqual(otherGeometry: GeometryQuery, distanceTol?: number, radianTol?: number): boolean; get isCircular(): boolean; get isExtensibleFractionSpace(): boolean; isInPlane(plane: Plane3dByOriginAndUnitNormal): boolean; @@ -1640,8 +1641,8 @@ export class CurveExtendOptions { export class CurveFactory { static appendToArcInPlace(arcA: Arc3d, arcB: Arc3d, allowReverse?: boolean): boolean; static assembleArcChainOnEllipsoid(ellipsoid: Ellipsoid, pathPoints: GeodesicPathPoint[], fractionForIntermediateNormal?: number): Path; - static createArcPointTangentPoint(pointA: Point3d, tangentA: Vector3d, pointB: Point3d): Arc3d | undefined; - static createArcPointTangentRadius(pointA: Point3d, tangentA: Vector3d, radius: number, upVector?: Vector3d, sweep?: Angle | AngleSweep): Arc3d | undefined; + static createArcPointTangentPoint(start: Point3d, tangentAtStart: Vector3d, end: Point3d): Arc3d | undefined; + static createArcPointTangentRadius(start: Point3d, tangentAtStart: Vector3d, radius: number, upVector?: Vector3d, sweep?: Angle | AngleSweep): Arc3d | undefined; static createFilletsInLineString(points: LineString3d | IndexedXYZCollection | Point3d[], radius: number | number[], allowBackupAlongEdge?: boolean): Path | undefined; static createLineSpiralArcSpiralLine(spiralType: IntegratedSpiralTypeName, pointA: Point3d, pointB: Point3d, pointC: Point3d, lengthA: number, lengthB: number, arcRadius: number): GeometryQuery[] | undefined; static createLineSpiralSpiralLine(spiralType: IntegratedSpiralTypeName, startPoint: Point3d, shoulderPoint: Point3d, targetPoint: Point3d): GeometryQuery[] | undefined; @@ -4611,11 +4612,11 @@ export class PolyfaceBuilder extends NullGeometryHandler { // @public export class PolyfaceClip { static clipPolyface(polyface: Polyface, clipper: ClipPlane | ConvexClipPlaneSet): Polyface | undefined; - static clipPolyfaceClipPlane(polyface: Polyface, clipper: ClipPlane, insideClip?: boolean, buildClosureFaces?: boolean): Polyface; + static clipPolyfaceClipPlane(polyface: Polyface, clipper: ClipPlane, insideClip?: boolean, buildClosureFaces?: boolean): IndexedPolyface; // @internal static clipPolyfaceClipPlaneToBuilders(polyface: Polyface, clipper: PlaneAltitudeEvaluator, destination: ClippedPolyfaceBuilders): void; - static clipPolyfaceClipPlaneWithClosureFace(polyface: Polyface, clipper: ClipPlane, insideClip?: boolean, buildClosureFaces?: boolean): Polyface; - static clipPolyfaceConvexClipPlaneSet(polyface: Polyface, clipper: ConvexClipPlaneSet): Polyface; + static clipPolyfaceClipPlaneWithClosureFace(polyface: Polyface, clipper: ClipPlane, insideClip?: boolean, buildClosureFaces?: boolean): IndexedPolyface; + static clipPolyfaceConvexClipPlaneSet(polyface: Polyface, clipper: ConvexClipPlaneSet): IndexedPolyface; // @internal static clipPolyfaceConvexClipPlaneSetToBuilders(polyface: Polyface, clipper: ConvexClipPlaneSet, destination: ClippedPolyfaceBuilders): void; static clipPolyfaceInsideOutside(polyface: Polyface, clipper: ClipPlane | ConvexClipPlaneSet | UnionOfConvexClipPlaneSets, destination: ClippedPolyfaceBuilders, outputSelect?: number): void; diff --git a/common/changes/@itwin/core-geometry/saeed-torabi-circularArcMethods_2024-10-04-20-42.json b/common/changes/@itwin/core-geometry/saeed-torabi-circularArcMethods_2024-10-04-20-42.json new file mode 100644 index 000000000000..9d7216f28509 --- /dev/null +++ b/common/changes/@itwin/core-geometry/saeed-torabi-circularArcMethods_2024-10-04-20-42.json @@ -0,0 +1,10 @@ +{ + "changes": [ + { + "packageName": "@itwin/core-geometry", + "comment": "", + "type": "none" + } + ], + "packageName": "@itwin/core-geometry" +} \ No newline at end of file diff --git a/core/geometry/internaldocs/Arc3d.md b/core/geometry/internaldocs/Arc3d.md new file mode 100644 index 000000000000..fa96deb0fda7 --- /dev/null +++ b/core/geometry/internaldocs/Arc3d.md @@ -0,0 +1,22 @@ +# Create Circular Arc using Start, TangentAtStart, and End + +Below we explain the algorithm used in `Arc3d.createCircularStartTangentEnd` to create a circular arc using start point, tangent at start, and end point. + +We first set up a frame of three perpendicular unit vectors using `Matrix3d.createRigidFromColumns`: +* `frameColX` lies along `tangentAtStart`, +* the circle normal lies along the cross product of `tangentAtStart` and the vector `startToEnd`, and +* `frameColY` lies along the cross product of the normal and `tangentAtStart`. + +We seek a formula for the radius of the circle. From the radius, we can find the circle center by moving a distance of radius from `start` along `frameColY`. From the center and the two input points, we can compute the arc sweep. Then we are done. + +The key to finding the radius is the inscribed right triangle pictured below, with one vertex at `start`, with hypotenuse along a diameter, and with leg along `startToEnd`. We know this is a right triangle from the inscribed angle theorem of classical geometry. + +Suppose `v = startToEnd`, `w = frameColY`, `r` is the radius we seek, and $\theta$ is the right triangle's angle at `start`. Then with + +$$v\cdot w = ||v|| ||w|| \cos\theta = ||v|| \frac{||v||}{2r} = \frac{||v||^2}{2r},$$ + +we have: + +$$\frac{v\cdot v}{2v\cdot w} = \frac{||v||^2}{2\frac{||v||^2}{2r}} = r.$$ + +![>](./figs/Arc3d/createCircularStartTangentEnd.png) \ No newline at end of file diff --git a/core/geometry/internaldocs/figs/Arc3d/createCircularStartTangentEnd.png b/core/geometry/internaldocs/figs/Arc3d/createCircularStartTangentEnd.png new file mode 100644 index 000000000000..ea49490926a2 Binary files /dev/null and b/core/geometry/internaldocs/figs/Arc3d/createCircularStartTangentEnd.png differ diff --git a/core/geometry/internaldocs/quarticRoots.md b/core/geometry/internaldocs/quarticRoots.md new file mode 100644 index 000000000000..425c828c4dcd --- /dev/null +++ b/core/geometry/internaldocs/quarticRoots.md @@ -0,0 +1,36 @@ +# Finding Roots of a Quartic Polynomial + +In `AnalyticRoots.appendQuarticRoots` we compute the solutions of the general quartic equation: +$$a_4x^4 + a_3x^3 + a_2x^2 + a_1x + a_0 = 0.\tag{1}$$ + +Without loss of generality, assume $a_4 = 1$. We further simplify by performing the substitution $x = y - a_3/4$ to obtain the depressed quartic: +$$P(x) = x^4 + px^2 + qx + r.$$ + +Note that the solutions of $(1)$ are obtained by subtracting $a_3/4$ from each root of $P$. + +Most classical solutions to $(1)$ derive a so-called _resolvent_ cubic polynomial $R$ from the depressed quartic, and use one of the roots of $R$ to construct a pair of quadratic equations whose roots are the roots of $P$. We will show this, using a resolvent cubic of this form: +$$R(y) = y^3 - \frac{1}{2}py^2 - ry + \frac{1}{2}rp - \frac{1}{8}q^2.$$ + +Let $x$ be a root of $P$. Then, introducing a quantity $y$ to complete the square, we have the following equivalences: + +$$\begin{align} +\notag{}P(x) = 0 &\Leftrightarrow x^4 = -px^2 - qx - r \\ +\notag{}&\Leftrightarrow (x^2 + y)^2 = -px^2 - qx - r + 2yx^2 + y^2 \\ +&\Leftrightarrow (x^2 + y)^2 = (2y - p)x^2 - qx + (y^2 - r)\tag{2} \\ +&\Leftrightarrow (x^2 + y)^2 = \left(\sqrt{2y - p}x - \frac{q}{2\sqrt{2y - p}}\right)^2\tag{3} \\ +&\Leftrightarrow \frac{q^2}{4(2y - p)} = y^2 - r\tag{4} \\ +\notag{}&\Leftrightarrow q^2 = 4(y^2 - r)(2y - p) \\ +\notag{}&\Leftrightarrow 0 = 8 R(y), +\end{align}$$ + +where in $(3)$ we have rewritten the right hand side of $(2)$ as a square (since the left hand side of $(2)$ is a square), and in $(4)$ we have equated the constant terms on the right hand sides of $(2)$ and $(3)$. Thus it is seen that the quantity $y$ is actually a root of the resolvent. We can compute the roots of $R$ with `AnalyticRoots.appendCubicRoots`. + +Lastly, choose a root $z$ of $R$. We need only one, so `AnalyticRoots.mostDistantFromMean` picks one employing a numerical stability criterion. Then from the above equivalences, and taking the square root of both sides of $(3)$, we have two quadratics in $x$: + +$$\begin{align} +\notag{}R(z) = 0 &\Leftrightarrow x^2 + z = \pm\left(\sqrt{2z - p}x - \frac{q}{2\sqrt{2z - p}}\right) \\ +&\Leftrightarrow x^2 \pm\sqrt{2z - p}x + z \mp\sqrt{z^2 - r} = 0,\tag{5} +\end{align}$$ + +where we have used the equality $(4)$ to simplify the constant term. We solve these two quadratics with `AnalyticRoots.appendQuadraticRoots`, yielding 0-4 values of $x$. By the above equivalences, each of these is a root of $P$. + diff --git a/core/geometry/src/curve/Arc3d.ts b/core/geometry/src/curve/Arc3d.ts index 671a296a4175..d07ca6322568 100644 --- a/core/geometry/src/curve/Arc3d.ts +++ b/core/geometry/src/curve/Arc3d.ts @@ -384,21 +384,21 @@ export class Arc3d extends CurvePrimitive implements BeJSONFunctions { } /** * Create an elliptical arc from three points on the ellipse: two points on an axis and one in between. - * @param point0 start of arc, on an axis - * @param point1 point on arc somewhere between `point0` and `point2` - * @param point2 point on arc directly opposite `point0` - * @param sweep angular sweep, measured from `point0` in the direction of `point1`. - * For a half-ellipse from `point0` to `point2` passing through `point1`, pass `AngleSweep.createStartEndDegrees(0,180)`. + * @param start start of arc, on an axis + * @param middle point on arc somewhere between `start` and `end` + * @param end point on arc directly opposite `start` + * @param sweep angular sweep, measured from `start` in the direction of `middle`. + * For a half-ellipse from `start` to `end` passing through `middle`, pass `AngleSweep.createStartEndDegrees(0,180)`. * Default value is full sweep to create the entire ellipse. * @param result optional preallocated result * @returns elliptical arc, or undefined if construction impossible. */ public static createStartMiddleEnd( - point0: XYAndZ, point1: XYAndZ, point2: XYAndZ, sweep?: AngleSweep, result?: Arc3d, + start: XYAndZ, middle: XYAndZ, end: XYAndZ, sweep?: AngleSweep, result?: Arc3d, ): Arc3d | undefined { - const center = Point3d.createAdd2Scaled(point0, 0.5, point2, 0.5); - const vector0 = Vector3d.createStartEnd(center, point0); - const vector1 = Vector3d.createStartEnd(center, point1); + const center = Point3d.createAdd2Scaled(start, 0.5, end, 0.5); + const vector0 = Vector3d.createStartEnd(center, start); + const vector1 = Vector3d.createStartEnd(center, middle); const v0DotV1 = vector0.dotProduct(vector1); const v0Len2 = vector0.magnitudeSquared(); if (Math.abs(v0DotV1) >= v0Len2) @@ -406,7 +406,7 @@ export class Arc3d extends CurvePrimitive implements BeJSONFunctions { const normal = vector0.crossProduct(vector1); const vector90 = normal.unitCrossProductWithDefault(vector0, 0, 0, 0); const v1DotV90 = vector1.dotProduct(vector90); - // Solve the standard ellipse equation for the unknown axis length, given local coords of point1 (v0.v1/||v0||, v90.v1) + // solve the standard ellipse equation for the unknown axis length, given local coords of middle (v0.v1/||v0||, v90.v1) const v90Len = Geometry.safeDivideFraction(v0Len2 * v1DotV90, Math.sqrt(v0Len2 * v0Len2 - v0DotV1 * v0DotV1), 0); if (Geometry.isSmallMetricDistanceSquared(v90Len)) return undefined; @@ -415,37 +415,59 @@ export class Arc3d extends CurvePrimitive implements BeJSONFunctions { } /** * Create a circular arc defined by start point, tangent at start point, and end point. - * If tangent is parallel to line segment from start to end, return the line segment. + * * The circular arc is swept from `start` to `end` in the direction of `tangentAtStart`. + * * If `tangentAtStart` is parallel to the line segment from `start` to `end`, return the line segment. */ public static createCircularStartTangentEnd( start: Point3d, tangentAtStart: Vector3d, end: Point3d, result?: Arc3d, ): Arc3d | LineSegment3d { - // To find the circle passing through start and end with tangentAtStart at start: - // - find line 1: the perpendicular bisector of the line from start to end. - // - find line 2: the perpendicular to the tangentAtStart. - // - intersection of the two lines would be the circle center. - const vector = Vector3d.createStartEnd(start, end); - const normal = tangentAtStart.crossProduct(vector).normalize(); - if (normal) { - const vectorPerp = normal.crossProduct(vector); - const tangentPerp = normal.crossProduct(tangentAtStart); - const midPoint = start.plusScaled(vector, 0.5); - - const lineSeg1 = LineSegment3d.create(start, start.plusScaled(tangentPerp, 1)); - const lineSeg2 = LineSegment3d.create(midPoint, midPoint.plusScaled(vectorPerp, 1)); - const intersection = LineSegment3d.closestApproach(lineSeg1, true, lineSeg2, true); - - if (intersection) { - const center = intersection.detailA.point; - const vector0 = Vector3d.createStartEnd(center, start); - const vector90 = normal.crossProduct(vector0); - const endVector = Vector3d.createStartEnd(center, end); - const sweep = AngleSweep.create(vector0.signedAngleTo(endVector, normal)); + // see itwinjs-core\core\geometry\internaldocs\Arc3d.md to clarify below algorithm + const startToEnd = Vector3d.createStartEnd(start, end); + const frame = Matrix3d.createRigidFromColumns(tangentAtStart, startToEnd, AxisOrder.XYZ); + if (frame !== undefined) { + const vv = startToEnd.dotProduct(startToEnd); + const vw = frame.dotColumnY(startToEnd); + const radius = Geometry.conditionalDivideCoordinate(vv, 2 * vw); + if (radius !== undefined) { + const vector0 = frame.columnY(); + vector0.scaleInPlace(-radius); // center to start + const vector90 = frame.columnX(); + vector90.scaleInPlace(radius); + const centerToEnd = vector0.plus(startToEnd); + const sweepAngle = vector0.angleTo(centerToEnd); + let sweepRadians = sweepAngle.radians; // always positive and less than PI + if (tangentAtStart.dotProduct(centerToEnd) < 0.0) // sweepRadians is the wrong way + sweepRadians = 2.0 * Math.PI - sweepRadians; + const center = start.plusScaled(vector0, -1.0); + const sweep = AngleSweep.createStartEndRadians(0.0, sweepRadians); return Arc3d.create(center, vector0, vector90, sweep, result); } } return LineSegment3d.create(start, end); } + /** + * Create a circular arc from start point, tangent at start, radius, optional plane normal, arc sweep. + * * The vector from start point to center is in the direction of upVector crossed with tangentA. + * @param start start point. + * @param tangentAtStart vector in tangent direction at the start. + * @param radius signed radius. + * @param upVector optional out-of-plane vector. Defaults to positive Z. + * @param sweep angular range. If single `Angle` is given, start angle is at 0 degrees (the start point). + */ + public static createCircularStartTangentRadius( + start: Point3d, tangentAtStart: Vector3d, radius: number, upVector?: Vector3d, sweep?: Angle | AngleSweep, + ): Arc3d | undefined { + if (upVector === undefined) + upVector = Vector3d.unitZ(); + const vector0 = upVector.unitCrossProduct(tangentAtStart); + if (vector0 === undefined) + return undefined; + const center = start.plusScaled(vector0, radius); + // reverse the A-to-center vector and bring it up to scale + vector0.scaleInPlace(-radius); + const vector90 = tangentAtStart.scaleToLength(Math.abs(radius))!; // cannot fail; prior unitCrossProduct would have failed first + return Arc3d.create(center, vector0, vector90, AngleSweep.create(sweep)); + } /** * Create a circular arc defined by start and end points and radius. * @param start start point of the arc @@ -910,7 +932,7 @@ export class Arc3d extends CurvePrimitive implements BeJSONFunctions { return this.isCircular ? this._matrix.columnXMagnitude() : undefined; } - /** Return the larger of the two defining vectors. */ + /** Return the larger length of the two defining vectors. */ public maxVectorLength(): number { return Math.max(this._matrix.columnXMagnitude(), this._matrix.columnYMagnitude()); } @@ -1138,12 +1160,12 @@ export class Arc3d extends CurvePrimitive implements BeJSONFunctions { }; } /** Test if this arc is almost equal to another GeometryQuery object */ - public override isAlmostEqual(otherGeometry: GeometryQuery): boolean { + public override isAlmostEqual(otherGeometry: GeometryQuery, distanceTol: number = Geometry.smallMetricDistance, radianTol: number = Geometry.smallAngleRadians): boolean { if (otherGeometry instanceof Arc3d) { const other = otherGeometry; - return this._center.isAlmostEqual(other._center) - && this._matrix.isAlmostEqual(other._matrix) - && this._sweep.isAlmostEqualAllowPeriodShift(other._sweep); + return this._center.isAlmostEqual(other._center, distanceTol) + && this._matrix.isAlmostEqual(other._matrix, distanceTol) + && this._sweep.isAlmostEqualAllowPeriodShift(other._sweep, radianTol); } return false; } @@ -1257,7 +1279,7 @@ export class Arc3d extends CurvePrimitive implements BeJSONFunctions { * * `point` is the `point1` input. * * both fractions are zero * * `arc` is undefined. - * @param point0 first point of path. (the point before the point of inflection) + * @param point0 first point of path (the point before the point of inflection) * @param point1 second point of path (the point of inflection) * @param point2 third point of path (the point after the point of inflection) * @param radius arc radius diff --git a/core/geometry/src/curve/CurveFactory.ts b/core/geometry/src/curve/CurveFactory.ts index 3a3fe29fc942..e789208ced16 100644 --- a/core/geometry/src/curve/CurveFactory.ts +++ b/core/geometry/src/curve/CurveFactory.ts @@ -103,37 +103,18 @@ export class CurveFactory { path.tryAddChild(LineSegment3d.create(pointA.interpolate(fraction0, pointB), pointA.interpolate(fraction1, pointB))); } } - /** - * Create a circular arc from start point, tangent at start, and another point (endpoint) on the arc. - * @param pointA - * @param tangentA - * @param pointB + * Create a circular arc defined by start point, tangent at start point, and end point. + * * The circular arc is swept from start to end toward direction of the `tangentAtStart`. + * * If tangent is parallel to line segment from start to end, return `undefined`. */ - public static createArcPointTangentPoint(pointA: Point3d, tangentA: Vector3d, pointB: Point3d): Arc3d | undefined { - const vectorV = Vector3d.createStartEnd(pointA, pointB); - const frame = Matrix3d.createRigidFromColumns(tangentA, vectorV, AxisOrder.XYZ); - if (frame !== undefined) { - const vv = vectorV.dotProduct(vectorV); - const vw = frame.dotColumnY(vectorV); - const alpha = Geometry.conditionalDivideCoordinate(vv, 2 * vw); - if (alpha !== undefined) { - const vector0 = frame.columnY(); - vector0.scaleInPlace(-alpha); - const vector90 = frame.columnX(); - vector90.scaleInPlace(alpha); - const centerToEnd = vector0.plus(vectorV); - const sweepAngle = vector0.angleTo(centerToEnd); - let sweepRadians = sweepAngle.radians; // That's always positive and less than PI. - if (tangentA.dotProduct(centerToEnd) < 0.0) // ah, sweepRadians is the wrong way - sweepRadians = 2.0 * Math.PI - sweepRadians; - const center = pointA.plusScaled(vector0, -1.0); - return Arc3d.create(center, vector0, vector90, AngleSweep.createStartEndRadians(0.0, sweepRadians)); - } - } - return undefined; + public static createArcPointTangentPoint(start: Point3d, tangentAtStart: Vector3d, end: Point3d): Arc3d | undefined { + const ret = Arc3d.createCircularStartTangentEnd(start, tangentAtStart, end); + if (ret instanceof Arc3d) + return ret; + else + return undefined; } - /** * Construct a sequence of alternating lines and arcs with the arcs creating tangent transition between consecutive edges. * * If the radius parameter is a number, that radius is used throughout. @@ -466,26 +447,18 @@ export class CurveFactory { } /** - * Create a circular arc from start point, tangent at start, radius, optional plane normal, arc sweep + * Create a circular arc from start point, tangent at start, radius, optional plane normal, arc sweep. * * The vector from start point to center is in the direction of upVector crossed with tangentA. - * @param pointA start point - * @param tangentA vector in tangent direction at the start + * @param start start point. + * @param tangentAtStart vector in tangent direction at the start. * @param radius signed radius. - * @param upVector optional out-of-plane vector. Defaults to positive Z - * @param sweep angular range. If single `Angle` is given, start angle is at 0 degrees (the start point). - * + * @param upVector optional out-of-plane vector. Defaults to positive Z. + * @param sweep angular range. If single `Angle` is given, start angle is at 0 degrees (the start point). */ - public static createArcPointTangentRadius(pointA: Point3d, tangentA: Vector3d, radius: number, upVector?: Vector3d, sweep?: Angle | AngleSweep): Arc3d | undefined { - if (upVector === undefined) - upVector = Vector3d.unitZ(); - const vector0 = upVector.unitCrossProduct(tangentA); - if (vector0 === undefined) - return undefined; - const center = pointA.plusScaled(vector0, radius); - // reverse the A-to-center vector and bring it up to scale ... - vector0.scaleInPlace(-radius); - const vector90 = tangentA.scaleToLength(Math.abs(radius))!; // (Cannot fail -- prior unitCrossProduct would have failed first) - return Arc3d.create(center, vector0, vector90, AngleSweep.create(sweep)); + public static createArcPointTangentRadius( + start: Point3d, tangentAtStart: Vector3d, radius: number, upVector?: Vector3d, sweep?: Angle | AngleSweep, + ): Arc3d | undefined { + return Arc3d.createCircularStartTangentRadius(start, tangentAtStart, radius, upVector, sweep); } /** diff --git a/core/geometry/src/curve/internalContexts/CurveCurveCloseApproachXY.ts b/core/geometry/src/curve/internalContexts/CurveCurveCloseApproachXY.ts index 172d0234fca9..fb65db5426dc 100644 --- a/core/geometry/src/curve/internalContexts/CurveCurveCloseApproachXY.ts +++ b/core/geometry/src/curve/internalContexts/CurveCurveCloseApproachXY.ts @@ -436,7 +436,7 @@ export class CurveCurveCloseApproachXY extends RecurseToCurvesGeometryHandler { */ private getPointCurveClosestApproachXYNewton(curveP: CurvePrimitive, pointQ: Point3d): CurveLocationDetail | undefined { if (!(curveP instanceof Arc3d) && !(curveP instanceof LineSegment3d)) { - assert(!"getPointCurveClosestApproachXYNewton only supports Arc3d and LineSegment"); + assert(false, "getPointCurveClosestApproachXYNewton only supports Arc3d and LineSegment"); return undefined; } const seeds = [0.2, 0.4, 0.6, 0.8]; // HEURISTIC: arcs have up to 4 perpendiculars; lines have only 1 @@ -757,7 +757,7 @@ export class CurveCurveCloseApproachXY extends RecurseToCurvesGeometryHandler { if (!this._geometryB || !(this._geometryB instanceof CurveChainWithDistanceIndex)) return; if (geomA instanceof CurveChainWithDistanceIndex) { - assert(!"call handleCurveChainWithDistanceIndex(geomA) instead"); + assert(false, "call handleCurveChainWithDistanceIndex(geomA) instead"); return; } const index0 = this._results.length; diff --git a/core/geometry/src/curve/internalContexts/CurveCurveIntersectXY.ts b/core/geometry/src/curve/internalContexts/CurveCurveIntersectXY.ts index 966da18972ea..8cfd4f65d3c9 100644 --- a/core/geometry/src/curve/internalContexts/CurveCurveIntersectXY.ts +++ b/core/geometry/src/curve/internalContexts/CurveCurveIntersectXY.ts @@ -688,7 +688,7 @@ export class CurveCurveIntersectXY extends RecurseToCurvesGeometryHandler { } const bcurveAFraction = bezierA.fractionToParentFraction(bezierAFraction); const bcurveBFraction = bezierB.fractionToParentFraction(bezierBFraction); - if (!"verify results") { + if (false) { // verify results const xyzA0 = bezierA.fractionToPoint(bezierAFraction); const xyzA1 = bcurveA.fractionToPoint(bcurveAFraction); const xyzB0 = bezierB.fractionToPoint(bezierBFraction); @@ -966,7 +966,7 @@ export class CurveCurveIntersectXY extends RecurseToCurvesGeometryHandler { if (!this._geometryB || !(this._geometryB instanceof CurveChainWithDistanceIndex)) return; if (geomA instanceof CurveChainWithDistanceIndex) { - assert(!"call handleCurveChainWithDistanceIndex(geomA) instead"); + assert(false, "call handleCurveChainWithDistanceIndex(geomA) instead"); return; } const index0 = this._results.length; diff --git a/core/geometry/src/curve/internalContexts/CurveCurveIntersectXYZ.ts b/core/geometry/src/curve/internalContexts/CurveCurveIntersectXYZ.ts index 4a74191aaa19..33a0881c53fe 100644 --- a/core/geometry/src/curve/internalContexts/CurveCurveIntersectXYZ.ts +++ b/core/geometry/src/curve/internalContexts/CurveCurveIntersectXYZ.ts @@ -796,7 +796,7 @@ export class CurveCurveIntersectXYZ extends RecurseToCurvesGeometryHandler { if (!this._geometryB || !(this._geometryB instanceof CurveChainWithDistanceIndex)) return; if (geomA instanceof CurveChainWithDistanceIndex) { - assert(!"call handleCurveChainWithDistanceIndex(geomA) instead"); + assert(false, "call handleCurveChainWithDistanceIndex(geomA) instead"); return; } const index0 = this._results.length; diff --git a/core/geometry/src/curve/internalContexts/EllipticalArcApproximationContext.ts b/core/geometry/src/curve/internalContexts/EllipticalArcApproximationContext.ts index 4e9ec9488e28..5f618b319b03 100644 --- a/core/geometry/src/curve/internalContexts/EllipticalArcApproximationContext.ts +++ b/core/geometry/src/curve/internalContexts/EllipticalArcApproximationContext.ts @@ -342,7 +342,7 @@ class AdaptiveSubdivisionQ1IntervalErrorProcessor extends QuadrantFractionsProce /** Remember the initial value of the fraction f to be perturbed. */ public override announceQuadrantBegin(q: QuadrantFractions, reversed: boolean): boolean { assert(q.quadrant === 1); - assert(!reversed); // ASSUME bracket and q.fractions have same ordering + assert(!reversed); // ASSUME [bracket0, bracket1] and q.fractions have the same ordering // the first fraction might be an extra point for computing the first 3-pt arc. assert(q.fractions.length === 4 || (q.fractions.length === 3 && q.interpolateStartTangent)); this._error0 = this._error1 = Geometry.largeCoordinateResult; @@ -454,7 +454,7 @@ class AdaptiveSubdivisionQ1ErrorProcessor extends QuadrantFractionsProcessor { const interpolateStartTangent = Geometry.isAlmostEqualEitherNumber(f0, this._fractionRangeQ1.low, this._fractionRangeQ1.high, 0); const interpolateEndTangent = Geometry.isAlmostEqualEitherNumber(f1, this._fractionRangeQ1.low, this._fractionRangeQ1.high, 0); if (!interpolateStartTangent && undefined === fPrev) - fPrev = this.getPreviousFraction(f0); // createLastArc caller doesn't supply fPrev + fPrev = this.getPreviousFraction(f0); // createLastArc doesn't supply fPrev to announceArc const fractions = (undefined === fPrev) ? [f0, f, f1] : [fPrev, f0, f, f1]; const q1 = [QuadrantFractions.create(1, fractions, interpolateStartTangent, interpolateEndTangent)]; const processor = AdaptiveSubdivisionQ1IntervalErrorProcessor.create(this.fullEllipseXY, f0, f, f1); @@ -817,12 +817,13 @@ export class EllipticalArcApproximationContext { arc.sweep.setStartEndRadians(startAngle.radians, arc.sweep.endRadians); return arc; // returned arc starts at arcStart, ends at arcEnd }; - const createFirstArc = (f0: number, f1: number, reverse: boolean): void => { + const createFirstArc = (f0: number, f1: number): void => { // This arc starts at the first sample f0 and ends at f1. + // This arc interpolates point and tangent at f0, but only point at f1. ellipticalArc.fractionToPointAndDerivative(f0, ray); ellipticalArc.fractionToPoint(f1, pt1); - if (reverse) - ray.direction.scaleInPlace(-1); + if (f0 > f1) + ray.direction.scaleInPlace(-1); // computed arc is retrograde const arc = arcBetween2Samples(ray, pt1, false); if (arc) processor.announceArc(arc, undefined, f0, f1); @@ -838,12 +839,14 @@ export class EllipticalArcApproximationContext { if (arc) processor.announceArc(arc, fPrev, f1, f2); }; - const createLastArc = (f0: number, f1: number, reverse: boolean): void => { - // This arc starts at f0 and ends at the last sample f1. It is the only arc to use f1. + const createLastArc = (f0: number, f1: number): void => { + // This arc starts at f0 and ends at the last sample f1. + // This arc interpolates point and tangent at f1, but only point at f0. ellipticalArc.fractionToPoint(f0, pt0); ellipticalArc.fractionToPointAndDerivative(f1, ray); - if (!reverse) - ray.direction.scaleInPlace(-1); + if (f1 > f0) + ray.direction.scaleInPlace(-1); // computed arc is retrograde + // compute last arc from f1 to f0, then reverse const arc = arcBetween2Samples(ray, pt0, true); if (arc) processor.announceArc(arc, undefined, f0, f1); @@ -874,13 +877,13 @@ export class EllipticalArcApproximationContext { if (!processor.announceQuadrantBegin(q, reversed)) continue; if (q.interpolateStartTangent) - createFirstArc(q.fractions[0], q.fractions[1], reversed); + createFirstArc(q.fractions[0], q.fractions[1]); // the first inner arc approximates the ellipse over [f[1],f[2]]; the last inner arc, over [f[n-3],f[n-2]] for (let i = 0; i + 2 < n - 1; ++i) createInnerArc(q.fractions[i], q.fractions[i + 1], q.fractions[i + 2]); if (n > 2) { // the final arc approximates [f[n-2],f[n-1]] if (q.interpolateEndTangent) - createLastArc(q.fractions[n - 2], q.fractions[n - 1], reversed); + createLastArc(q.fractions[n - 2], q.fractions[n - 1]); else createInnerArc(q.fractions[n - 3], q.fractions[n - 2], q.fractions[n - 1]); } diff --git a/core/geometry/src/geometry3d/AngleSweep.ts b/core/geometry/src/geometry3d/AngleSweep.ts index 44e689c8651f..8cda7035aa02 100644 --- a/core/geometry/src/geometry3d/AngleSweep.ts +++ b/core/geometry/src/geometry3d/AngleSweep.ts @@ -523,28 +523,31 @@ export class AngleSweep implements BeJSONFunctions { } /** * Convert an AngleSweep to a JSON object. - * @return {*} {degrees: [startAngleInDegrees, endAngleInDegrees} + * @return {*} [startAngleInDegrees, endAngleInDegrees] */ public toJSON(): any { return [this.startDegrees, this.endDegrees]; } /** - * Test if this angle sweep and other angle sweep match with radians tolerance. - * * Period shifts are allowed. + * Test if two angle sweeps match within the given tolerance. + * * Period shifts are allowed, but orientations must be the same. + * @param other sweep to compare to this instance + * @param radianTol optional radian tolerance, default value `Geometry.smallAngleRadians` */ - public isAlmostEqualAllowPeriodShift(other: AngleSweep): boolean { - // We compare angle sweeps by checking if start angle and sweep match. We cannot compare start and end because for - // example (0, 90) and (360, 90) have the same start (we allow period shift) and end but are not same angle sweeps. - return Angle.isAlmostEqualRadiansAllowPeriodShift(this._radians0, other._radians0) - && Angle.isAlmostEqualRadiansAllowPeriodShift(this._radians1 - this._radians0, other._radians1 - other._radians0); + public isAlmostEqualAllowPeriodShift(other: AngleSweep, radianTol: number = Geometry.smallAngleRadians): boolean { + return this.isCCW === other.isCCW // this rules out equating opposite sweeps like [0,-100] and [0,260] + && Angle.isAlmostEqualRadiansAllowPeriodShift(this._radians0, other._radians0, radianTol) + && Angle.isAlmostEqualRadiansAllowPeriodShift(this._radians1 - this._radians0, other._radians1 - other._radians0, radianTol); } /** - * Test if this angle sweep and other angle sweep match with radians tolerance. + * Test if two angle sweeps match within the given tolerance. * * Period shifts are not allowed. + * @param other sweep to compare to this instance + * @param radianTol optional radian tolerance, default value `Geometry.smallAngleRadians` */ - public isAlmostEqualNoPeriodShift(other: AngleSweep): boolean { - return Angle.isAlmostEqualRadiansNoPeriodShift(this._radians0, other._radians0) - && Angle.isAlmostEqualRadiansNoPeriodShift(this._radians1 - this._radians0, other._radians1 - other._radians0); + public isAlmostEqualNoPeriodShift(other: AngleSweep, radianTol: number = Geometry.smallAngleRadians): boolean { + return Angle.isAlmostEqualRadiansNoPeriodShift(this._radians0, other._radians0, radianTol) + && Angle.isAlmostEqualRadiansNoPeriodShift(this._radians1 - this._radians0, other._radians1 - other._radians0, radianTol); } /** * Test if start and end angles match with radians tolerance. diff --git a/core/geometry/src/geometry3d/PolygonOps.ts b/core/geometry/src/geometry3d/PolygonOps.ts index a13988b44fe6..c8d8bacb3ed5 100644 --- a/core/geometry/src/geometry3d/PolygonOps.ts +++ b/core/geometry/src/geometry3d/PolygonOps.ts @@ -1262,7 +1262,7 @@ export class PolygonOps { const areaOfNormalParallelogram = Math.abs(outwardUnitNormalOfPrevEdge.crossProductXY(outwardUnitNormalOfEdge)); const coord = Geometry.conditionalDivideCoordinate(areaOfNormalParallelogram, projToPrevEdge.x * projToEdge.x, largestResult); if (undefined === coord) { - assert(!"unexpectedly small projection distance to an edge"); + assert(false, "unexpectedly small projection distance to an edge"); return undefined; // shouldn't happen due to chopping in computeEdgeDataXY: area/(dist*dist) <= 1/tol^2 = largestResult } coords[i] = coord; @@ -1273,7 +1273,7 @@ export class PolygonOps { } const scale = Geometry.conditionalDivideCoordinate(1.0, coordSum); if (undefined === scale) { - assert(!"unexpected zero barycentric coordinate sum"); + assert(false, "unexpected zero barycentric coordinate sum"); return undefined; } for (let i = 0; i < n; ++i) diff --git a/core/geometry/src/geometry4d/Map4d.ts b/core/geometry/src/geometry4d/Map4d.ts index ee151c8c74f0..60817c81d0ce 100644 --- a/core/geometry/src/geometry4d/Map4d.ts +++ b/core/geometry/src/geometry4d/Map4d.ts @@ -111,11 +111,11 @@ export class Map4d implements BeJSONFunctions { return this._matrix0.isAlmostEqual(other._matrix0) && this._matrix1.isAlmostEqual(other._matrix1); } /** - * Create a map between a frustum and world coordinates. - * @param origin lower left of frustum - * @param uVector Vector from lower left rear to lower right rear. - * @param vVector Vector from lower left rear to upper left rear. - * @param wVector Vector from lower left rear to lower left front, i.e. lower left rear towards eye. + * Create a world to NPC map that maps between world coordinates and the given frustum. + * @param origin lower left rear of frustum + * @param uVector Vector from origin to lower right rear. + * @param vVector Vector from origin to upper left rear. + * @param wVector Vector from origin to lower left front, i.e. origin towards eye. * @param fraction front size divided by rear size. */ public static createVectorFrustum( diff --git a/core/geometry/src/geometry4d/Matrix4d.ts b/core/geometry/src/geometry4d/Matrix4d.ts index 4a39098193ed..f5e7d3b19736 100644 --- a/core/geometry/src/geometry4d/Matrix4d.ts +++ b/core/geometry/src/geometry4d/Matrix4d.ts @@ -30,7 +30,7 @@ export type Matrix4dProps = Point4dProps[]; * * indices 8,9,10,11 are the "z row" They may be called the zx,zy,zz,zw entries * * indices 12,13,14,15 are the "w row". They may be called the wx,wy,wz,ww entries * * If "w row" contains numeric values 0,0,0,1, the Matrix4d is equivalent to a Transform with - * * The upper right 3x3 matrix (entries 0,1,2,4,5,6,8,9,10) are the 3x3 matrix part of the transform + * * The upper left 3x3 matrix (entries 0,1,2,4,5,6,8,9,10) are the 3x3 matrix part of the transform * * The far right column entries xw,yw,zw are the "origin" (sometimes called "translation") part of the transform. * @public */ @@ -190,7 +190,7 @@ export class Matrix4d implements BeJSONFunctions { return Matrix4d.createRowValues(scaleX, 0, 0, tx, 0, scaleY, 0, ty, 0, 0, scaleZ, tz, 0, 0, 0, 1, result); } /** - * Create a mapping the scales and translates (no rotation) from box A to boxB + * Create a mapping that scales and translates (no rotation) from box A to box B * @param lowA low point of box A * @param highA high point of box A * @param lowB low point of box B @@ -502,7 +502,7 @@ export class Matrix4d implements BeJSONFunctions { } /** multiply each matrix * points[i]. This produces a weighted xyzw. * Immediately renormalize back to xyz and replace the original point. - * If zero weight appears in the result (i.e. input is on eyeplane)leave the mapped xyz untouched. + * If zero weight appears in the result (i.e. input is on eyeplane) leave the mapped xyz untouched. */ public multiplyPoint3dArrayQuietNormalize(points: Point3d[]) { points.forEach((point) => this.multiplyXYZWQuietRenormalize(point.x, point.y, point.z, 1.0, point)); @@ -686,7 +686,7 @@ export class Matrix4d implements BeJSONFunctions { this._coffs[15] += a * vectorV.w; } /** - * ADD (n place) scale*A*B*AT where + * Add (in place) scale*A*B*AT where * * A is a pure translation with final column [x,y,z,1] * * B is the given `matrixB` * * AT is the transpose of A. @@ -735,12 +735,9 @@ export class Matrix4d implements BeJSONFunctions { * * A is a pure translation with final column [x,y,z,1] * * this is this matrix. * * AT is the transpose of A. - * * scale is a multiplier. - * @param matrixB the middle matrix. * @param ax x part of translation * @param ay y part of translation * @param az z part of translation - * @param scale scale factor for entire product */ public multiplyTranslationSandwichInPlace(ax: number, ay: number, az: number) { const bx = this._coffs[3]; diff --git a/core/geometry/src/numerics/Polynomials.ts b/core/geometry/src/numerics/Polynomials.ts index b8629ede19c4..aec17e6c02a2 100644 --- a/core/geometry/src/numerics/Polynomials.ts +++ b/core/geometry/src/numerics/Polynomials.ts @@ -28,7 +28,7 @@ import { Point4d } from "../geometry4d/Point4d"; * @internal */ export class Degree2PowerPolynomial { - /** The three coefficients for the quartic */ + /** The three coefficients for the quadratic */ public coffs: number[]; constructor(c0: number = 0, c1: number = 0, c2: number = 0) { @@ -913,7 +913,7 @@ export class AnalyticRoots { // EDL April 5, 2020 replace classic GraphicsGems solver by RWDNickalls. // Don't know if improveRoots is needed. // Breaks in AnalyticRoots.test.ts checkQuartic suggest it indeed converts many e-16 errors to zero. - // e-13 cases are unaffected + // e-13 cases are unaffected this.improveRoots(c, 3, results, false); } else { this.appendQuadraticRoots(c, results); @@ -923,6 +923,7 @@ export class AnalyticRoots { } /** Compute roots of quartic `c[0] + c[1] * x + c[2] * x^2 + c[3] * x^3 + c[4] * x^4` */ public static appendQuarticRoots(c: Float64Array | number[], results: GrowableFloat64Array) { + // for details, see core\geometry\internaldocs\quarticRoots.md const coffs = new Float64Array(4); let u: number; let v: number; @@ -952,7 +953,7 @@ export class AnalyticRoots { results.push(0); this.addConstant(origin, results); // apply origin return; - } else { // solve the resolvent cubic; more info: https://en.wikipedia.org/wiki/Resolvent_cubic#Second_definition + } else { // solve the resolvent cubic coffs[0] = 0.5 * r * p - 0.125 * q * q; coffs[1] = -r; coffs[2] = -0.5 * p; @@ -985,7 +986,6 @@ export class AnalyticRoots { coffs[2] = 1; this.appendQuadraticRoots(coffs, results); } - // substitute this.addConstant(origin, results); // apply origin results.sort(); this.improveRoots(c, 4, results, true); @@ -1120,8 +1120,7 @@ export class TrigPolynomial { // tolerance for small angle decision. private static readonly _smallAngle: number = 1.0e-11; - // see itwinjs-core\core\geometry\internaldocs\unitCircleEllipseIntersection.md - // on how below variables are derived. + // see core\geometry\internaldocs\unitCircleEllipseIntersection.md for derivation of these coefficients. /** Standard Basis coefficients for the numerator of the y-coordinate y(t) = S(t)/W(t) in the rational semicircle parameterization. */ public static readonly S = Float64Array.from([0.0, 2.0, -2.0]); /** Standard Basis coefficients for the numerator of the x-coordinate x(t) = C(t)/W(t) in the rational semicircle parameterization. */ @@ -1146,7 +1145,7 @@ export class TrigPolynomial { /** * Solve a polynomial created from trigonometric condition using Trig.S, Trig.C, Trig.W. * * Polynomial is of degree 4: - * `coff[0] + coff[1] * t + coff[2] * t^2 + coff[3] * t^3 + coff[4] * t^4` + * `p(t) = coff[0] + coff[1] * t + coff[2] * t^2 + coff[3] * t^3 + coff[4] * t^4` * * Solution logic includes inferring angular roots corresponding zero leading coefficients * (roots at infinity). * @param coff coefficients. @@ -1176,14 +1175,12 @@ export class TrigPolynomial { degree--; const roots = new GrowableFloat64Array(); if (degree === -1) { - // Umm. Dunno. Nothing there. + // do nothing } else { if (degree === 0) { - // p(t) is a nonzero constant - // No roots, but not degenerate. + // p(t) is a nonzero constant; no roots but not degenerate. } else if (degree === 1) { - // p(t) = coff[1] * t + coff[0] - roots.push(- coff[0] / coff[1]); + roots.push(-coff[0] / coff[1]); // p(t) = coff[0] + coff[1] * t } else if (degree === 2) { AnalyticRoots.appendQuadraticRoots(coff, roots); } else if (degree === 3) { @@ -1194,17 +1191,15 @@ export class TrigPolynomial { // TODO: WORK WITH BEZIER SOLVER } if (roots.length > 0) { - // Each solution t represents an angle with - // Math.Cos(theta) = C(t)/W(t) and sin(theta) = S(t)/W(t) + // Each solution t represents an angle with Math.Cos(theta) = C(t)/W(t) and sin(theta) = S(t)/W(t) // Division by W has no effect on atan2 calculations, so we just compute S(t),C(t) for (let i = 0; i < roots.length; i++) { const ss = PowerPolynomial.evaluate(this.S, roots.atUncheckedIndex(i)); const cc = PowerPolynomial.evaluate(this.C, roots.atUncheckedIndex(i)); radians.push(Math.atan2(ss, cc)); } - // Each leading zero at the front of the coefficients corresponds to a root at -PI/2. - // Only make one entry.... - // for (int i = degree; i < nominalDegree; i++) + // each leading zero at the front of the coefficient array corresponds to a root at -PI/2. + // only make one entry because we don't report multiplicity. if (degree < nominalDegree) radians.push(-0.5 * Math.PI); } @@ -1230,8 +1225,7 @@ export class TrigPolynomial { const coffs = new Float64Array(5); PowerPolynomial.zero(coffs); let degree; - // see itwinjs-core\core\geometry\internaldocs\unitCircleEllipseIntersection.md - // on how coffs (coefficient array) is built. + // see core\geometry\internaldocs\unitCircleEllipseIntersection.md for derivation of these coefficients if (Geometry.hypotenuseXYZ(axx, axy, ayy) > TrigPolynomial._coefficientRelTol * Geometry.hypotenuseXYZ(ax, ay, a)) { PowerPolynomial.accumulate(coffs, this.CW, ax); PowerPolynomial.accumulate(coffs, this.SW, ay); @@ -1274,9 +1268,14 @@ export class TrigPolynomial { * @param ellipseRadians solution angles in ellipse parameter space * @param circleRadians solution angles in circle parameter space */ - public static solveUnitCircleEllipseIntersection(cx: number, cy: number, ux: number, uy: number, - vx: number, vy: number, ellipseRadians: number[], circleRadians: number[]): boolean { + public static solveUnitCircleEllipseIntersection( + cx: number, cy: number, + ux: number, uy: number, + vx: number, vy: number, + ellipseRadians: number[], circleRadians: number[], + ): boolean { circleRadians.length = 0; + // see core\geometry\internaldocs\unitCircleEllipseIntersection.md for derivation of these coefficients: const acc = ux * ux + uy * uy; const acs = 2.0 * (ux * vx + uy * vy); const ass = vx * vx + vy * vy; @@ -1315,8 +1314,7 @@ export class TrigPolynomial { ellipseRadians: number[], circleRadians: number[], ): boolean { circleRadians.length = 0; - // see itwinjs-core\core\geometry\internaldocs\unitCircleEllipseIntersection.md - // on how below variables are derived. + // see core\geometry\internaldocs\unitCircleEllipseIntersection.md for derivation of these coefficients: const acc = ux * ux + uy * uy - uw * uw; const acs = 2.0 * (ux * vx + uy * vy - uw * vw); const ass = vx * vx + vy * vy - vw * vw; diff --git a/core/geometry/src/polyface/PolyfaceClip.ts b/core/geometry/src/polyface/PolyfaceClip.ts index 63ce8aa170f3..62aa49fa296d 100644 --- a/core/geometry/src/polyface/PolyfaceClip.ts +++ b/core/geometry/src/polyface/PolyfaceClip.ts @@ -99,14 +99,14 @@ export class PolyfaceClip { * * Return all surviving clip as a new mesh. * * WARNING: The new mesh is "points only" -- parameters, normals, etc are not interpolated */ - public static clipPolyfaceClipPlaneWithClosureFace(polyface: Polyface, clipper: ClipPlane, insideClip: boolean = true, buildClosureFaces: boolean = true) { + public static clipPolyfaceClipPlaneWithClosureFace(polyface: Polyface, clipper: ClipPlane, insideClip: boolean = true, buildClosureFaces: boolean = true): IndexedPolyface { return this.clipPolyfaceClipPlane(polyface, clipper, insideClip, buildClosureFaces); } /** Clip each facet of polyface to the ClipPlane. * * Return all surviving clip as a new mesh. * * WARNING: The new mesh is "points only" -- parameters, normals, etc are not interpolated */ - public static clipPolyfaceClipPlane(polyface: Polyface, clipper: ClipPlane, insideClip: boolean = true, buildClosureFaces: boolean = false): Polyface { + public static clipPolyfaceClipPlane(polyface: Polyface, clipper: ClipPlane, insideClip: boolean = true, buildClosureFaces: boolean = false): IndexedPolyface { const builders = ClippedPolyfaceBuilders.create(insideClip, !insideClip, buildClosureFaces); this.clipPolyfaceInsideOutside(polyface, clipper, builders); return builders.claimPolyface(insideClip ? 0 : 1, true)!; @@ -116,7 +116,7 @@ export class PolyfaceClip { * * Return surviving clip as a new mesh. * * WARNING: The new mesh is "points only". */ - public static clipPolyfaceConvexClipPlaneSet(polyface: Polyface, clipper: ConvexClipPlaneSet): Polyface { + public static clipPolyfaceConvexClipPlaneSet(polyface: Polyface, clipper: ConvexClipPlaneSet): IndexedPolyface { const visitor = polyface.createVisitor(0); const builder = PolyfaceBuilder.create(); const work = new GrowableXYZArray(10); @@ -215,7 +215,7 @@ export class PolyfaceClip { for (const child of region.children) this.addRegion(builder, child); } else { - assert(!"unexpected region encountered"); + assert(false, "unexpected region encountered"); } } } diff --git a/core/geometry/src/test/curve/Arc3d.test.ts b/core/geometry/src/test/curve/Arc3d.test.ts index baaf4a98d8e3..ba53972fdc38 100644 --- a/core/geometry/src/test/curve/Arc3d.test.ts +++ b/core/geometry/src/test/curve/Arc3d.test.ts @@ -6,6 +6,7 @@ import { assert, describe, expect, it } from "vitest"; import { compareNumbers, OrderedSet } from "@itwin/core-bentley"; import { Constant } from "../../Constant"; +import { CurveFactory } from "../../core-geometry"; import { Arc3d, EllipticalArcApproximationOptions, EllipticalArcSampleMethod, FractionMapper } from "../../curve/Arc3d"; import { CoordinateXYZ } from "../../curve/CoordinateXYZ"; import { CurveChainWithDistanceIndex } from "../../curve/CurveChainWithDistanceIndex"; @@ -716,6 +717,15 @@ describe("Arc3d", () => { ck.testPoint3d(circularArc4.endPoint(), end4); ck.testPoint3d(circularArc4.center, Point3d.create(0.75, 0, 0.75)); + dx += 10; + const start5 = Point3d.create(10, 0, 0); + const end5 = Point3d.create(7.0710678118654755, 2.1213203435596424, 0); + const tangent5 = Vector3d.create(-0, -18.84955592153876, -0); + const circularArc5A = Arc3d.createCircularStartTangentEnd(start5, tangent5, end5) as Arc3d; + const circularArc5B = CurveFactory.createArcPointTangentPoint(start5, tangent5, end5) as Arc3d; + ck.testTrue(circularArc5A.isAlmostEqual(circularArc5B, 1.0e-13, 1.0e-13), "methods are equivalent"); + GeometryCoreTestIO.captureCloneGeometry(allGeometry, [circularArc5A, circularArc5B], dx); + GeometryCoreTestIO.saveGeometry(allGeometry, "Arc3d", "createCircularStartTangentEnd"); expect(ck.getNumErrors()).toBe(0); }); @@ -1054,7 +1064,6 @@ describe("ApproximateArc3d", () => { } } } - y0 += yDelta(yWidth); } if ( @@ -1229,13 +1238,13 @@ describe("ApproximateArc3d", () => { x += delta; } - // Observed: subdivision wins 90.9% of comparisons to n-sample methods (95.67% with enableLongTests) + // Observed: subdivision wins 90.9% of comparisons to n-sample methods (95.7% with enableLongTests) const winPct = 100 * Geometry.safeDivideFraction(nSubdivisionComparisonWins, nComparisons, 0); GeometryCoreTestIO.consoleLog(`Subdivision wins ${nSubdivisionComparisonWins} of ${nComparisons} comparisons (${winPct}%).`); const targetWinPct = GeometryCoreTestIO.enableLongTests ? 90 : 85; ck.testLE(targetWinPct, winPct, `Subdivision is more accurate than another n-sample method over ${targetWinPct}% of the time.`); - // Observed: subdivision is most accurate method in 64% of ellipses tested (82.76% with enableLongTests) + // Observed: subdivision is most accurate method in 63.6% of ellipses tested (82.8% with enableLongTests) const winOverallPct = 100 * Geometry.safeDivideFraction(nEllipses - nSubdivisionLosses, nEllipses, 0); GeometryCoreTestIO.consoleLog(`Subdivision wins overall for ${nEllipses - nSubdivisionLosses} of ${nEllipses} ellipses (${winOverallPct}%).`); const targetNSampleWinPct = GeometryCoreTestIO.enableLongTests ? 80 : 60; diff --git a/core/geometry/src/test/numerics/Polynomial.test.ts b/core/geometry/src/test/numerics/Polynomial.test.ts index 6a2d62ca85ca..117d18f3eda0 100644 --- a/core/geometry/src/test/numerics/Polynomial.test.ts +++ b/core/geometry/src/test/numerics/Polynomial.test.ts @@ -708,3 +708,67 @@ function sphereSnakeGridPoints(center: Point3d, radius: number, thetaArray: numb } return points; } + +function captureUnitCircleEllipseIntersections( + allGeometry: any[], ck: any, + cx: number, cy: number, + ux: number, uy: number, + vx: number, vy: number, + dx: number, expectedIntersections: number, +) { + const unitCircle = Arc3d.create( + Point3d.create(0, 0), Vector3d.create(1, 0), Vector3d.create(0, 1), AngleSweep.createStartEndDegrees(), + ); + const arc = Arc3d.create( + Point3d.create(cx, cy), Vector3d.create(ux, uy), Vector3d.create(vx, vy), + AngleSweep.createStartEndDegrees(), + ); + GeometryCoreTestIO.captureCloneGeometry(allGeometry, unitCircle, dx); + GeometryCoreTestIO.captureCloneGeometry(allGeometry, arc, dx); + + const ellipseRadians: number[] = []; + const circleRadians: number[] = []; + TrigPolynomial.solveUnitCircleEllipseIntersection(cx, cy, ux, uy, vx, vy, ellipseRadians, circleRadians); + const len = ellipseRadians.length; + for (let i = 0; i < len; i++) { + const fraction = arc.sweep.radiansToSignedFraction(ellipseRadians[i]); + GeometryCoreTestIO.createAndCaptureXYCircle(allGeometry, arc.fractionToPoint(fraction), 0.2, dx); + } + ck.testExactNumber(len, expectedIntersections, `${expectedIntersections} intersection(s) expected`); +} + +it("unitCircleEllipseIntersection", () => { + const ck = new Checker(); + const allGeometry: GeometryQuery[] = []; + let dx = 0; + // intersection at lower half of the ellipse + let cx = 0, cy = 1.5; // ellipse center + let ux = 2, uy = 0; // ellipse vector0 + let vx = 0, vy = 1; // ellipse vector90 + let expectedIntersections = 2; + captureUnitCircleEllipseIntersections(allGeometry, ck, cx, cy, ux, uy, vx, vy, dx, expectedIntersections); + // intersection at t = 0 + dx += 12; + cx = 3, cy = 0; + ux = 2, uy = 0; + vx = 0, vy = 1; + expectedIntersections = 1; + captureUnitCircleEllipseIntersections(allGeometry, ck, cx, cy, ux, uy, vx, vy, dx, expectedIntersections); + // intersection at t = 1 + dx += 12; + cx = -3, cy = 0; + ux = 2, uy = 0; + vx = 0, vy = 1; + expectedIntersections = 1; + captureUnitCircleEllipseIntersections(allGeometry, ck, cx, cy, ux, uy, vx, vy, dx, expectedIntersections); + // intersection at t = infinity + dx += 12; + cx = 0, cy = 2; + ux = 2, uy = 0; + vx = 0, vy = 1; + expectedIntersections = 1; + // captureUnitCircleEllipseIntersections(allGeometry, ck, cx, cy, ux, uy, vx, vy, dx, expectedIntersections); + + GeometryCoreTestIO.saveGeometry(allGeometry, "CurveCurveIntersectXY", "unitCircleEllipseIntersection"); + expect(ck.getNumErrors()).toBe(0); +});