Skip to content

Commit

Permalink
Consolidated circular arc construction methods into Arc3d (#7239)
Browse files Browse the repository at this point in the history
Co-authored-by: dassaf4 <68340676+dassaf4@users.noreply.github.com>
  • Loading branch information
saeeedtorabi and dassaf4 authored Oct 9, 2024
1 parent 0b08dba commit 3269edd
Show file tree
Hide file tree
Showing 19 changed files with 302 additions and 164 deletions.
19 changes: 10 additions & 9 deletions common/api/core-geometry.api.md
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand Down Expand Up @@ -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;
Expand All @@ -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;
Expand Down Expand Up @@ -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;
Expand Down Expand Up @@ -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;
Expand Down
Original file line number Diff line number Diff line change
@@ -0,0 +1,10 @@
{
"changes": [
{
"packageName": "@itwin/core-geometry",
"comment": "",
"type": "none"
}
],
"packageName": "@itwin/core-geometry"
}
22 changes: 22 additions & 0 deletions core/geometry/internaldocs/Arc3d.md
Original file line number Diff line number Diff line change
@@ -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)
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
36 changes: 36 additions & 0 deletions core/geometry/internaldocs/quarticRoots.md
Original file line number Diff line number Diff line change
@@ -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$.

98 changes: 60 additions & 38 deletions core/geometry/src/curve/Arc3d.ts
Original file line number Diff line number Diff line change
Expand Up @@ -384,29 +384,29 @@ 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)
return undefined; // middle point projects to end of axis or beyond (rules out negative under the radical below)
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;
Expand All @@ -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
Expand Down Expand Up @@ -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());
}
Expand Down Expand Up @@ -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;
}
Expand Down Expand Up @@ -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
Expand Down
Loading

0 comments on commit 3269edd

Please sign in to comment.