From 695233345d0a7d05859e6a4cc20292cf45d57b80 Mon Sep 17 00:00:00 2001 From: James Bradley Date: Tue, 18 Jul 2023 14:49:54 -0400 Subject: [PATCH 01/22] init Polyline implementation of IHasArcLength --- CHANGELOG.md | 2 + .../src/Geometry/Interfaces/IHasArcLength.cs | 38 ++++++ Elements/src/Geometry/Polyline.cs | 109 +++++++++++++++++- Elements/test/PolylineTests.cs | 71 ++++++++++++ 4 files changed, 218 insertions(+), 2 deletions(-) create mode 100644 Elements/src/Geometry/Interfaces/IHasArcLength.cs diff --git a/CHANGELOG.md b/CHANGELOG.md index 80b9c034b..43b2f9105 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -15,6 +15,7 @@ - `Ellipse` - `EllipticalArc` - `IndexedPolycurve` +- `IHasArcLength` - `Grid1d.GetCellDomains` - `Message.Info` - `Message.Error` @@ -27,6 +28,7 @@ - `Polyline` now inherits from `BoundedCurve`. - `Polyline` is now parameterized 0->length. +- `Polyline` now implements the `IHasArcLength` interface - `Arc` now inherits from `TrimmedCurve`. - `Arc` is now parameterized 0->2Pi - `Line` now inherits from `TrimmedCurve`. diff --git a/Elements/src/Geometry/Interfaces/IHasArcLength.cs b/Elements/src/Geometry/Interfaces/IHasArcLength.cs new file mode 100644 index 000000000..fe8a3cc90 --- /dev/null +++ b/Elements/src/Geometry/Interfaces/IHasArcLength.cs @@ -0,0 +1,38 @@ +namespace Elements.Geometry.Interfaces +{ + /// + /// Represents a curve with arc length-based operations. + /// Implementing classes define methods for computing points and performing operations based on arc length. + /// + public interface IHasArcLength + { + /// + /// Returns the point on the curve at the specified arc length. + /// + /// The arc length along the curve. + /// The point on the curve at the specified arc length. + Vector3 PointAtLength(double length); + + /// + /// Returns the point on the curve at the specified normalized length. + /// The normalized length is a value between 0 and 1 representing the relative position along the curve. + /// + /// The normalized length along the curve. + /// The point on the curve at the specified normalized length. + Vector3 PointAtNormalizedLength(double normalizedLength); + + /// + /// Returns the midpoint of the curve. + /// + /// The midpoint of the curve. + Vector3 MidPoint(); + + /// + /// Divides the curve into segments of the specified length and returns the points along the curve at those intervals. + /// + /// The desired length for dividing the curve. + /// A list of points representing the divisions along the curve. + Vector3[] DivideByLength(double length); + } + +} \ No newline at end of file diff --git a/Elements/src/Geometry/Polyline.cs b/Elements/src/Geometry/Polyline.cs index b33e77545..6544fb497 100644 --- a/Elements/src/Geometry/Polyline.cs +++ b/Elements/src/Geometry/Polyline.cs @@ -4,6 +4,7 @@ using System.Linq; using ClipperLib; using Elements.Search; +using Elements.Geometry.Interfaces; namespace Elements.Geometry { @@ -14,7 +15,7 @@ namespace Elements.Geometry /// /// [!code-csharp[Main](../../Elements/test/PolylineTests.cs?name=example)] /// - public class Polyline : IndexedPolycurve + public class Polyline : IndexedPolycurve, IHasArcLength { /// /// The domain of the curve. @@ -84,6 +85,66 @@ public virtual Line[] Segments() return SegmentsInternal(this.Vertices); } + + + /// + /// The mid point of the curve. + /// + /// The length based midpoint. + public virtual Vector3 MidPoint() + { + return PointAtNormalizedLength(0.5); + } + + /// + /// Returns the point on the polyline corresponding to the specified length value. + /// + /// The length value along the polyline. + /// The point on the polyline corresponding to the specified length value. + /// Thrown when the specified length is out of range. + public virtual Vector3 PointAtLength(double length) + { + double totalLength = ArcLength(this.Domain.Min, this.Domain.Max); // Calculate the total length of the Polyline + if (length < 0 || length > totalLength) + { + throw new ArgumentException("The specified length is out of range."); + } + + double accumulatedLength = 0.0; + foreach (Line segment in Segments()) + { + double segmentLength = segment.ArcLength(segment.Domain.Min, segment.Domain.Max); + + if (accumulatedLength + segmentLength >= length) + { + double remainingDistance = length - accumulatedLength; + double parameter = remainingDistance / segmentLength; + return segment.PointAtNormalized(parameter); + } + + accumulatedLength += segmentLength; + } + + // If we reach here, the desired length is equal to the total length, + // so return the end point of the Polyline. + return End; + } + + /// + /// Returns the point on the polyline corresponding to the specified normalized length-based parameter value. + /// + /// The normalized length-based parameter value, ranging from 0 to 1. + /// The point on the polyline corresponding to the specified normalized length-based parameter value. + /// Thrown when the specified parameter is out of range. + public virtual Vector3 PointAtNormalizedLength(double parameter) + { + if (parameter < 0 || parameter > 1) + { + throw new ArgumentException("The specified parameter is out of range."); + } + return PointAtLength(parameter * this.ArcLength(this.Domain.Min, this.Domain.Max)); + } + /// /// Get the transform at the specified parameter along the polyline. /// @@ -386,6 +447,50 @@ private Transform CreateOrthogonalTransform(int i, Vector3 origin, Vector3 up) return new Transform(origin, up.Cross(tangent), tangent); } + /// + /// Divides the polyline into segments of the specified length. + /// + /// The desired length of each segment. + /// A list of points representing the segments. + public Vector3[] DivideByLength(double divisionLength) + { + var segments = new List(); + + if (this.Vertices.Count < 2) + { + // Handle invalid polyline with insufficient vertices + return new Vector3[0]; + } + + var currentProgression = 0.0; + segments = new List { this.Vertices.FirstOrDefault() }; + + foreach (var currentSegment in this.Segments()) + { + // currentProgression from last segment before hitting end + if (currentProgression != 0.0) + { + currentProgression -= divisionLength; + } + while (currentSegment.ArcLength(currentSegment.Domain.Min, currentSegment.Domain.Max) >= currentProgression + divisionLength) + { + segments.Add(currentSegment.PointAt(currentProgression + divisionLength)); + currentProgression += divisionLength; + } + // Set currentProgression from divisionLength less distance from last segment point + currentProgression = divisionLength - segments.LastOrDefault().DistanceTo(currentSegment.End); + } + + // Add the last vertex of the polyline as the endpoint of the last segment if it + // is not already part of the list + if (!segments.LastOrDefault().IsAlmostEqualTo(this.Vertices.LastOrDefault())) + { + segments.Add(this.Vertices.LastOrDefault()); + } + + return segments.ToArray(); + } + /// /// Offset this polyline by the specified amount. /// @@ -795,7 +900,7 @@ protected void Split(IList points, bool closed = false) var b = closed && i == this.Vertices.Count - 1 ? this.Vertices[0] : this.Vertices[i + 1]; var edge = (a, b); - // An edge may have multiple split points. + // An edge may have multiple split points. // We store these in a list and sort it along the // direction of the edge, before inserting the points // into the vertex list and incrementing i by the correct diff --git a/Elements/test/PolylineTests.cs b/Elements/test/PolylineTests.cs index 2ce511e3b..db6c4fe16 100644 --- a/Elements/test/PolylineTests.cs +++ b/Elements/test/PolylineTests.cs @@ -257,6 +257,77 @@ public void GetParameterAt() Assert.Equal(1.5, resultParameter); var testPoint = polyline.PointAt(resultParameter); Assert.True(point.IsAlmostEqualTo(testPoint)); + + polyline = new Polyline( + new List() + { + new Vector3(1, 5, 0), + new Vector3(5, 20, 0), + new Vector3(5, 0, 0), + new Vector3(9, 5, 0), + } + ); + var parameter = 0.5; + var testMidpoint = new Vector3(5.0, 14.560525, 0); // Midpoint + var testParameterMidpoint = polyline.PointAtNormalizedLength(parameter); + var t = testParameterMidpoint.IsAlmostEqualTo(testMidpoint); + Assert.True(testParameterMidpoint.IsAlmostEqualTo(testMidpoint)); + + var midlength = polyline.Length() * parameter; + var testLengthMidpoint = polyline.PointAtLength(midlength); + Assert.True(testLengthMidpoint.IsAlmostEqualTo(testParameterMidpoint)); + + var midpoint = polyline.MidPoint(); + Assert.True(midpoint.IsAlmostEqualTo(testLengthMidpoint)); + } + + [Fact] + public void DivideByLength_SingleSegment() + { + // Arrange + var start = new Vector3(1, 2); + var end = new Vector3(1, 4); + var polyline = new Polyline(new[] { start, end }); + var divisionLength = 1.0; + + // Act + var segments = polyline.DivideByLength(divisionLength); + + // Assert + Assert.Equal(3, segments.Count()); + Assert.Equal(start, segments[0]); + Assert.Equal(new Vector3(1, 3), segments[1]); + Assert.Equal(end, segments[2]); + } + + [Fact] + public void DivideByLength_MultipleSegments() + { + // Arrange + var polyline = new Polyline(new List() + { + new Vector3(1, 5, 0), + new Vector3(5, 20, 0), + new Vector3(5, 0, 0), + new Vector3(9, 5, 0), + }); + var divisionLength = 5.0; + + // Act + var segments = polyline.DivideByLength(divisionLength); + + // Assert + Assert.Equal(10, segments.Count()); + Assert.Equal(new Vector3(1.0, 5.0, 0), segments[0]); + Assert.Equal(new Vector3(2.288313, 9.831174, 0), segments[1]); + Assert.Equal(new Vector3(3.576626, 14.662349, 0), segments[2]); + Assert.Equal(new Vector3(4.864939, 19.493524, 0), segments[3]); + Assert.Equal(new Vector3(5.0, 15.524174, 0), segments[4]); + Assert.Equal(new Vector3(5.0, 10.524174, 0), segments[5]); + Assert.Equal(new Vector3(5.0, 5.524174, 0), segments[6]); + Assert.Equal(new Vector3(5.0, 0.524174, 0), segments[7]); + Assert.Equal(new Vector3(7.796025, 3.495032, 0), segments[8]); + Assert.Equal(new Vector3(9, 5, 0), segments[9]); } [Fact] From 7309fc7adc116489c8c4edf4df6558e322e1fe36 Mon Sep 17 00:00:00 2001 From: James Bradley Date: Tue, 18 Jul 2023 17:36:40 -0400 Subject: [PATCH 02/22] Line implements Interface IHasArcLength --- CHANGELOG.md | 1 + Elements/src/Elements.csproj | 2 +- Elements/src/Geometry/Line.cs | 136 ++++++++++++++++++++++++------ Elements/src/Geometry/Polyline.cs | 2 - Elements/test/LineTests.cs | 21 +++-- Elements/test/PolylineTests.cs | 1 - 6 files changed, 128 insertions(+), 35 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 43b2f9105..807ba4f6a 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -33,6 +33,7 @@ - `Arc` is now parameterized 0->2Pi - `Line` now inherits from `TrimmedCurve`. - `Line` is now parameterized 0->length. +- `Line` now implements the `IHasArcLength` interface - `Bezier` now inherits from `BoundedCurve`. - `Polyline` is now parameterized 0->length. - `Circle` is now parameterized 0->2Pi. diff --git a/Elements/src/Elements.csproj b/Elements/src/Elements.csproj index cdd21865a..8bc4905ea 100644 --- a/Elements/src/Elements.csproj +++ b/Elements/src/Elements.csproj @@ -10,7 +10,7 @@ The Elements library provides object types for generating the built environment. https://github.com/hypar-io/elements https://github.com/hypar-io/elements - $(Version) + 21.21.21 en diff --git a/Elements/src/Geometry/Line.cs b/Elements/src/Geometry/Line.cs index cea31dd95..2c33a1438 100644 --- a/Elements/src/Geometry/Line.cs +++ b/Elements/src/Geometry/Line.cs @@ -4,6 +4,7 @@ using System.Linq; using Newtonsoft.Json; using Elements.Spatial; +using Elements.Geometry.Interfaces; namespace Elements.Geometry { @@ -15,7 +16,7 @@ namespace Elements.Geometry /// [!code-csharp[Main](../../Elements/test/LineTests.cs?name=example)] /// /// TODO: Rename this class to LineSegment - public class Line : TrimmedCurve, IEquatable + public class Line : TrimmedCurve, IEquatable, IHasArcLength { /// /// The domain of the curve. @@ -123,6 +124,48 @@ public override Vector3 PointAt(double u) return this.BasisCurve.PointAt(u); } + /// + /// The mid point of the curve. + /// + /// The length based midpoint. + public virtual Vector3 MidPoint() + { + return PointAtNormalizedLength(0.5); + } + + /// + /// Returns the point on the line corresponding to the specified length value. + /// + /// The length value along the line. + /// The point on the line corresponding to the specified length value. + /// Thrown when the specified length is out of range. + public virtual Vector3 PointAtLength(double length) + { + double totalLength = ArcLength(this.Domain.Min, this.Domain.Max); // Calculate the total length of the Line + + if (length < 0 || length > totalLength) + { + throw new ArgumentException("The specified length is out of range."); + } + var lengthParameter = length / totalLength; + return this.PointAtNormalized(lengthParameter); + } + + /// + /// Returns the point on the line corresponding to the specified normalized length-based parameter value. + /// + /// The normalized length-based parameter value, ranging from 0 to 1. + /// The point on the line corresponding to the specified normalized length-based parameter value. + /// Thrown when the specified parameter is out of range. + public virtual Vector3 PointAtNormalizedLength(double parameter) + { + if (parameter < 0 || parameter > 1) + { + throw new ArgumentException("The specified parameter is out of range."); + } + return PointAtLength(parameter * this.ArcLength(this.Domain.Min, this.Domain.Max)); + } + /// public override Curve Transformed(Transform transform) { @@ -541,38 +584,79 @@ public static bool PointOnLine(Vector3 point, Vector3 start, Vector3 end, bool i return false; } + // /// + // /// Divide the line into as many segments of the provided length as possible. + // /// + // /// The length. + // /// A flag indicating whether segments shorter than l should be removed. + // public List DivideByLength(double l, bool removeShortSegments = false) + // { + // var len = this.Length(); + // if (l > len) + // { + // return new List() { new Line(this.Start, this.End) }; + // } + + // var total = 0.0; + // var d = this.Direction(); + // var lines = new List(); + // while (total + l <= len) + // { + // var a = this.Start + d * total; + // var b = a + d * l; + // lines.Add(new Line(a, b)); + // total += l; + // } + // if (total < len && !removeShortSegments) + // { + // var a = this.Start + d * total; + // if (!a.IsAlmostEqualTo(End)) + // { + // lines.Add(new Line(a, End)); + // } + // } + // return lines; + // } + /// - /// Divide the line into as many segments of the provided length as possible. + /// Divides the line into segments of the specified length. /// - /// The length. - /// A flag indicating whether segments shorter than l should be removed. - public List DivideByLength(double l, bool removeShortSegments = false) + /// The desired length of each segment. + /// A list of points representing the segments. + public Vector3[] DivideByLength(double divisionLength) { - var len = this.Length(); - if (l > len) + var segments = new List(); + + if (this.ArcLength(this.Domain.Min, this.Domain.Max) < double.Epsilon) { - return new List() { new Line(this.Start, this.End) }; + // Handle invalid line with insufficient length + return new Vector3[0]; } - var total = 0.0; - var d = this.Direction(); - var lines = new List(); - while (total + l <= len) + var currentProgression = 0.0; + segments = new List { this.Start }; + + // currentProgression from last segment before hitting end + if (currentProgression != 0.0) { - var a = this.Start + d * total; - var b = a + d * l; - lines.Add(new Line(a, b)); - total += l; + currentProgression -= divisionLength; } - if (total < len && !removeShortSegments) + while (this.ArcLength(this.Domain.Min, this.Domain.Max) >= currentProgression + divisionLength) { - var a = this.Start + d * total; - if (!a.IsAlmostEqualTo(End)) - { - lines.Add(new Line(a, End)); - } + segments.Add(this.PointAt(currentProgression + divisionLength)); + currentProgression += divisionLength; } - return lines; + // Set currentProgression from divisionLength less distance from last segment point + currentProgression = divisionLength - segments.LastOrDefault().DistanceTo(this.End); + + // Add the last vertex of the polyline as the endpoint of the last segment if it + // is not already part of the list + if (!segments.LastOrDefault().IsAlmostEqualTo(this.End)) + { + segments.Add(this.End); + } + + return segments.ToArray(); } /// @@ -925,7 +1009,7 @@ public double DistanceTo(Line other) // line vectors are not collinear, their directions share the common plane. else { - // dStartStart length is distance to the common plane. + // dStartStart length is distance to the common plane. dStartStart = dStartStart.ProjectOnto(cross); Vector3 vStartStart = other.Start + dStartStart - this.Start; Vector3 vStartEnd = other.Start + dStartStart - this.End; @@ -1028,8 +1112,8 @@ public List Trim(Polygon polygon, out List outsideSegments, bool inc var B = intersectionsOrdered[i + 1]; if (A.IsAlmostEqualTo(B)) // skip duplicate points { - // it's possible that A is outside, but B is at an edge, even - // if they are within tolerance of each other. + // it's possible that A is outside, but B is at an edge, even + // if they are within tolerance of each other. // This can happen due to floating point error when the point is almost exactly // epsilon distance from the edge. // so if we have duplicate points, we have to update the containment value. diff --git a/Elements/src/Geometry/Polyline.cs b/Elements/src/Geometry/Polyline.cs index 6544fb497..81628dbdd 100644 --- a/Elements/src/Geometry/Polyline.cs +++ b/Elements/src/Geometry/Polyline.cs @@ -85,8 +85,6 @@ public virtual Line[] Segments() return SegmentsInternal(this.Vertices); } - - /// /// The mid point of the curve. /// diff --git a/Elements/test/LineTests.cs b/Elements/test/LineTests.cs index 1703cdedd..18eb03589 100644 --- a/Elements/test/LineTests.cs +++ b/Elements/test/LineTests.cs @@ -163,11 +163,11 @@ public void DivideIntoEqualSegmentsSingle() public void DivideByLength() { var l = new Line(Vector3.Origin, new Vector3(5, 0)); - var segments = l.DivideByLength(1.1, true); - Assert.Equal(4, segments.Count); + var segments = l.DivideByLength(1.1); + Assert.Equal(6, segments.Count()); - var segments1 = l.DivideByLength(1.1); - Assert.Equal(5, segments1.Count); + var segments1 = l.DivideByLength(2); + Assert.Equal(4, segments1.Count()); } [Fact] @@ -650,6 +650,17 @@ public void GetParameterAt() var expectedVector = line.PointAt(uValue); Assert.InRange(uValue, 0, line.Length()); Assert.True(vector.IsAlmostEqualTo(expectedVector)); + + var parameter = 0.5; + var testParameterMidpoint = line.PointAtNormalizedLength(parameter); + Assert.True(testParameterMidpoint.IsAlmostEqualTo(middle)); + + var midlength = line.Length() * parameter; + var testLengthMidpoint = line.PointAtLength(midlength); + Assert.True(testLengthMidpoint.IsAlmostEqualTo(testParameterMidpoint)); + + var midpoint = line.MidPoint(); + Assert.True(midpoint.IsAlmostEqualTo(testLengthMidpoint)); } [Theory] @@ -1110,7 +1121,7 @@ public void LineDistancePointsOnSkewLines() Assert.Equal(delta.Length(), (new Line(pt12, pt11)).DistanceTo(new Line(pt21, pt22)), 12); Assert.Equal(delta.Length(), (new Line(pt12, pt11)).DistanceTo(new Line(pt22, pt21)), 12); //The segments (pt12, pt13) and (pt21, pt22) does not intersect. - //The shortest distance is from an endpoint to another segment - difference between lines plus between endpoints. + //The shortest distance is from an endpoint to another segment - difference between lines plus between endpoints. var expected = (q12 * v1).DistanceTo(new Line(delta + q21 * v2, delta + q22 * v2)); Assert.Equal(expected, (new Line(pt12, pt13)).DistanceTo(new Line(pt21, pt22)), 12); Assert.Equal(expected, (new Line(pt12, pt13)).DistanceTo(new Line(pt22, pt21)), 12); diff --git a/Elements/test/PolylineTests.cs b/Elements/test/PolylineTests.cs index db6c4fe16..e60245f4a 100644 --- a/Elements/test/PolylineTests.cs +++ b/Elements/test/PolylineTests.cs @@ -270,7 +270,6 @@ public void GetParameterAt() var parameter = 0.5; var testMidpoint = new Vector3(5.0, 14.560525, 0); // Midpoint var testParameterMidpoint = polyline.PointAtNormalizedLength(parameter); - var t = testParameterMidpoint.IsAlmostEqualTo(testMidpoint); Assert.True(testParameterMidpoint.IsAlmostEqualTo(testMidpoint)); var midlength = polyline.Length() * parameter; From 2afac10dcf65f0ec9b497da5bd38b90764144755 Mon Sep 17 00:00:00 2001 From: James Bradley Date: Tue, 18 Jul 2023 21:41:01 -0400 Subject: [PATCH 03/22] Circle implements Interface IHasArcLength --- CHANGELOG.md | 2 +- Elements/src/Geometry/Circle.cs | 183 +++++++++++++++++++++++++++++++- Elements/test/CircleTests.cs | 141 ++++++++++++++++++++++++ Elements/test/LineTests.cs | 1 - 4 files changed, 323 insertions(+), 4 deletions(-) create mode 100644 Elements/test/CircleTests.cs diff --git a/CHANGELOG.md b/CHANGELOG.md index 807ba4f6a..99251a7fd 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -37,7 +37,7 @@ - `Bezier` now inherits from `BoundedCurve`. - `Polyline` is now parameterized 0->length. - `Circle` is now parameterized 0->2Pi. -- `Line` is now parameterized 0->length. +- `Circle` now implements the `IHasArcLength` interface - `Vector3.DistanceTo(Ray ray)` now returns positive infinity instead of throwing. - `Message`: removed obsolete `FromLine` method. - `AdaptiveGrid`: removed obsolete `TryGetVertexIndex` with `tolerance` parameter. diff --git a/Elements/src/Geometry/Circle.cs b/Elements/src/Geometry/Circle.cs index 4e075593a..70fce6816 100644 --- a/Elements/src/Geometry/Circle.cs +++ b/Elements/src/Geometry/Circle.cs @@ -1,3 +1,4 @@ +using Elements.Validators; using System; using System.Collections.Generic; using Elements.Geometry.Interfaces; @@ -6,10 +7,10 @@ namespace Elements.Geometry { /// - /// A circle. + /// A circle. /// Parameterization of the circle is 0 -> 2PI. /// - public class Circle : Curve, IConic + public class Circle : Curve, IConic, IHasArcLength { /// The center of the circle. [JsonProperty("Center", Required = Required.AllowNull)] @@ -26,6 +27,17 @@ public Vector3 Center [System.ComponentModel.DataAnnotations.Range(0.0D, double.MaxValue)] public double Radius { get; protected set; } + /// The circumference of the circle. + [JsonIgnore] + [System.ComponentModel.DataAnnotations.Range(0.0D, double.MaxValue)] + public double Circumference { get; protected set; } + + /// + /// The domain of the curve. + /// + [JsonIgnore] + public Domain1d Domain => new Domain1d(0, 2 * Math.PI); + /// /// The coordinate system of the plane containing the circle. /// @@ -39,7 +51,15 @@ public Vector3 Center [JsonConstructor] public Circle(Vector3 center, double radius = 1.0) { + if (!Validator.DisableValidationOnConstruction) + { + if (Math.Abs(radius - 0.0) < double.Epsilon ? true : false) + { + throw new ArgumentException($"The circle could not be created. The radius of the circle cannot be the zero: radius {radius}"); + } + } this.Radius = radius; + this.Circumference = 2 * Math.PI * this.Radius; this.Transform = new Transform(center); } @@ -49,7 +69,15 @@ public Circle(Vector3 center, double radius = 1.0) /// The radius of the circle. public Circle(double radius = 1.0) { + if (!Validator.DisableValidationOnConstruction) + { + if (Math.Abs(radius - 0.0) < double.Epsilon ? true : false) + { + throw new ArgumentException($"The circle could not be created. The radius of the circle cannot be the zero: radius {radius}"); + } + } this.Radius = radius; + this.Circumference = 2 * Math.PI * this.Radius; this.Transform = new Transform(); } @@ -58,8 +86,65 @@ public Circle(double radius = 1.0) /// public Circle(Transform transform, double radius = 1.0) { + if (!Validator.DisableValidationOnConstruction) + { + if (Math.Abs(radius - 0.0) < double.Epsilon ? true : false) + { + throw new ArgumentException($"The circle could not be created. The radius of the circle cannot be the zero: radius {radius}"); + } + } this.Transform = transform; this.Radius = radius; + this.Circumference = 2 * Math.PI * this.Radius; + } + + /// + /// Calculate the length of the circle between two parameters. + /// + public double ArcLength(double start, double end) + { + // Convert start and end parameters from radians to degrees + double _startAngle = start * 180.0 / Math.PI; + double _endAngle = end * 180.0 / Math.PI; + + // Ensure the start angle is within the valid domain range of 0 to 360 degrees + double startAngle = _startAngle % 360; + if (startAngle < 0) + { + startAngle += 360; + } + + // Ensure the end angle is within the valid domain range of 0 to 360 degrees + double endAngle = _endAngle % 360; + if (endAngle < 0) + { + endAngle += 360; + } + else if (endAngle == 0 && Math.Abs(_endAngle) >= 2 * Math.PI) + { + endAngle = 360; + } + + // Calculate the difference in angles + double angleDifference = endAngle - startAngle; + + // Adjust the angle difference if it crosses the 360-degree boundary + if (angleDifference < 0) + { + angleDifference += 360; + } + else if (angleDifference >= 2 * Math.PI) + { + return Circumference; // Full circle, return circumference + } + + // Convert the angle difference back to radians + double angleDifferenceRadians = angleDifference * Math.PI / 180.0; + + // Calculate the arc length using the formula: arc length = radius * angle + double arcLength = Radius * angleDifferenceRadians; + + return arcLength; } /// @@ -83,6 +168,14 @@ public Polygon ToPolygon(int divisions = 10) return new Polygon(pts, true); } + /// + /// Are the two circles almost equal? + /// + public bool IsAlmostEqualTo(Circle other, double tolerance = Vector3.EPSILON) + { + return (Center.IsAlmostEqualTo(other.Center, tolerance) && Math.Abs(Radius - other.Radius) < tolerance ? true : false); + } + /// /// Convert a circle to a circular arc. /// @@ -94,6 +187,15 @@ public Polygon ToPolygon(int divisions = 10) /// The bounded curve to convert. public static implicit operator ModelCurve(Circle c) => new ModelCurve(c); + /// + /// Calculates and returns the midpoint of the circle. + /// + /// The midpoint of the circle. + public Vector3 MidPoint() + { + return PointAt(Math.PI); + } + /// /// Return the point at parameter u on the arc. /// @@ -111,6 +213,46 @@ private Vector3 PointAtUntransformed(double u) return new Vector3(x, y); } + /// + /// Calculates and returns the point on the circle at a specific arc length. + /// + /// The arc length along the circumference of the circle. + /// The point on the circle at the specified arc length. + public Vector3 PointAtLength(double length) + { + double parameter = (length / Circumference) * 2 * Math.PI; + return PointAt(parameter); + } + + /// + /// Calculates and returns the point on the circle at a normalized arc length. + /// + /// The normalized arc length between 0 and 1. + /// The point on the circle at the specified normalized arc length. + public Vector3 PointAtNormalizedLength(double normalizedLength) + { + double parameter = normalizedLength * 2 * Math.PI; + return PointAt(parameter); + } + + /// + /// Calculates the parameter within the range of 0 to 2π at a given point on the circle. + /// + /// The point on the circle. + /// The parameter within the range of 0 to 2π at the given point on the circle. + public double GetParameterAt(Vector3 point) + { + Vector3 relativePoint = point - Center; + + double theta = Math.Atan2(relativePoint.Y, relativePoint.X); + + if (theta < 0) + { + theta += 2 * Math.PI; + } + return theta; + } + /// /// Return transform on the arc at parameter u. /// @@ -148,5 +290,42 @@ public override double ParameterAtDistanceFromParameter(double distance, double return start + theta; } + + /// + /// Divides the circle into segments of the specified length and returns a list of points representing the division. + /// + /// The length of each segment. + /// A list of points representing the division of the circle. + public Vector3[] DivideByLength(double length) + { + List points = new List(); + double circumference = 2 * Math.PI * Radius; + int segmentCount = (int)Math.Ceiling(circumference / length); + double segmentLength = circumference / segmentCount; + + for (int i = 0; i < segmentCount; i++) + { + double parameter = i * segmentLength / circumference; + points.Add(PointAtNormalizedLength(parameter)); + } + + return points.ToArray(); + } + + /// + /// Checks if a given point lies on a circle within a specified tolerance. + /// + /// The point to be checked. + /// The circle to check against. + /// The tolerance value (optional). Default is 1E-05. + /// True if the point lies on the circle within the tolerance, otherwise false. + public static bool PointOnCircle(Vector3 point, Circle circle, double tolerance = 1E-05) + { + Vector3 centerToPoint = point - circle.Center; + double distanceToCenter = centerToPoint.Length(); + + // Check if the distance from the point to the center is within the tolerance of the circle's radius + return Math.Abs(distanceToCenter - circle.Radius) < tolerance; + } } } \ No newline at end of file diff --git a/Elements/test/CircleTests.cs b/Elements/test/CircleTests.cs new file mode 100644 index 000000000..4e423ec76 --- /dev/null +++ b/Elements/test/CircleTests.cs @@ -0,0 +1,141 @@ +using System; +using System.Collections.Generic; +using System.IO; +using System.Linq; +using Elements.Tests; +using Newtonsoft.Json; +using Xunit; + +namespace Elements.Geometry.Tests +{ + public class CircleTests : ModelTest + { + public CircleTests() + { + this.GenerateIfc = false; + } + + [Fact, Trait("Category", "Examples")] + public void CircleExample() + { + this.Name = "Elements_Geometry_Circle"; + + // + var a = new Vector3(); + var b = 1.0; + var c = new Circle(a, b); + // + + this.Model.AddElement(c); + } + + [Fact] + public void Equality() + { + var p = 1.0; + var circleA = new Circle(Vector3.Origin, p); + var circleB = new Circle(Vector3.Origin, p + 1E-4); + var circleC = new Circle(Vector3.Origin, p + 1E-6); + + Assert.False(circleA.IsAlmostEqualTo(circleB)); + Assert.True(circleA.IsAlmostEqualTo(circleB, 1E-3)); + Assert.True(circleA.IsAlmostEqualTo(circleC)); + } + + [Fact] + public void Construct() + { + var a = new Vector3(); + var b = 1.0; + var c = new Circle(a, b); + Assert.Equal(1.0, c.Radius); + Assert.Equal(new Vector3(0, 0), c.Center); + Assert.Equal(a + new Vector3(b, 0, 0), c.PointAt(-1e-10)); + } + + [Fact] + public void ZeroRadius_ThrowsException() + { + var a = new Vector3(); + Assert.Throws(() => new Circle(a, 0)); + } + + [Fact] + public void DivideIntoEqualSegments() + { + var l = new Line(Vector3.Origin, new Vector3(100, 0)); + var segments = l.DivideIntoEqualSegments(41); + var len = l.Length(); + Assert.Equal(41, segments.Count); + foreach (var s in segments) + { + Assert.Equal(s.Length(), len / 41, 5); + } + } + + [Fact] + public void DivideIntoEqualSegmentsSingle() + { + var l = new Line(Vector3.Origin, new Vector3(100, 0)); + var segments = l.DivideIntoEqualSegments(1); + Assert.Single(segments); + Assert.True(segments.First().Start.IsAlmostEqualTo(l.Start, 1e-10)); + Assert.True(segments.First().End.IsAlmostEqualTo(l.End, 1e-10)); + } + + [Fact] + public void DivideByLength() + { + var l = new Line(Vector3.Origin, new Vector3(5, 0)); + var segments = l.DivideByLength(1.1); + Assert.Equal(6, segments.Count()); + + var segments1 = l.DivideByLength(2); + Assert.Equal(4, segments1.Count()); + } + + [Fact] + public void GetParameterAt() + { + var center = Vector3.Origin; + var radius = 5.0; + var circle = new Circle(center, radius); + var start = new Vector3(5.0, 0.0, 0.0); + var mid = new Vector3(-5.0, 0.0, 0.0); + + Assert.Equal(0, circle.GetParameterAt(start)); + + var almostEqualStart = new Vector3(5.000001, 0.000005, 0); + Assert.True(start.IsAlmostEqualTo(almostEqualStart)); + + Assert.True(Math.Abs(circle.GetParameterAt(mid) - Math.PI) < double.Epsilon ? true : false); + + var vector = new Vector3(3.535533, 3.535533, 0.0); + var uValue = circle.GetParameterAt(vector); + var expectedVector = circle.PointAt(uValue); + Assert.InRange(uValue, 0, circle.Circumference); + Assert.True(vector.IsAlmostEqualTo(expectedVector)); + + var parameter = 0.5; + var testParameterMidpoint = circle.PointAtNormalizedLength(parameter); + Assert.True(testParameterMidpoint.IsAlmostEqualTo(mid)); + + var midlength = circle.Circumference * parameter; + var testLengthMidpoint = circle.PointAtLength(midlength); + Assert.True(testLengthMidpoint.IsAlmostEqualTo(testParameterMidpoint)); + + var midpoint = circle.MidPoint(); + Assert.True(midpoint.IsAlmostEqualTo(testLengthMidpoint)); + } + + [Fact] + public void PointOnLine() + { + Circle circle = new Circle(Vector3.Origin, 5.0); + + Assert.False(Circle.PointOnCircle(Vector3.Origin, circle)); + Assert.False(Circle.PointOnCircle(new Vector3(4, 0, 0), circle)); + Assert.True(Circle.PointOnCircle(circle.PointAtNormalizedLength(0.5), circle)); + } + } +} diff --git a/Elements/test/LineTests.cs b/Elements/test/LineTests.cs index 18eb03589..19e9ddc3b 100644 --- a/Elements/test/LineTests.cs +++ b/Elements/test/LineTests.cs @@ -618,7 +618,6 @@ public void TryGetOverlap() Assert.Equal(new Line((5, 0, 0), (10, 0, 0)), overlap); } - [Fact] public void GetParameterAt() { From 374fe3466c01d050a83b44574847491d39978cde Mon Sep 17 00:00:00 2001 From: James Bradley Date: Fri, 21 Jul 2023 15:46:54 -0400 Subject: [PATCH 04/22] Polish Bezier, IHasCurveLength --- CHANGELOG.md | 18 +- Elements/src/Geometry/Bezier.cs | 586 +++++++++++++++++++++++++++++- Elements/src/Geometry/Polyline.cs | 4 +- Elements/test/BezierTests.cs | 161 ++++++++ 4 files changed, 751 insertions(+), 18 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 99251a7fd..5a9edd9ff 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -23,21 +23,31 @@ - `Topography.Trimmed` - `new Topography(Topography other)` - `Topography.TopMesh()` +- `Bezier.PointAtLength()` +- `Bezier.PointAtNormalizedLength()` +- `Bezier.ParameterAt()` +- `Bezier.DivideByLength()` +- `Bezier.Split()` +- `Bezier.SplitAt()` +- `Bezier.SplitByLength()` +- `Bezier.ConstructPiecewiseCubicBezier()` ### Changed - `Polyline` now inherits from `BoundedCurve`. - `Polyline` is now parameterized 0->length. -- `Polyline` now implements the `IHasArcLength` interface +- `Polyline` now implements the `IHasArcLength` interface. - `Arc` now inherits from `TrimmedCurve`. -- `Arc` is now parameterized 0->2Pi +- `Arc` is now parameterized 0->2Pi. - `Line` now inherits from `TrimmedCurve`. - `Line` is now parameterized 0->length. -- `Line` now implements the `IHasArcLength` interface +- `Line` now implements the `IHasArcLength` interface. - `Bezier` now inherits from `BoundedCurve`. +- `Bezier` now implements the `IHasArcLength` interface. +- `Bezier.ArcLength()` now uses Gauss quadrature approximation vs linear sampling. - `Polyline` is now parameterized 0->length. - `Circle` is now parameterized 0->2Pi. -- `Circle` now implements the `IHasArcLength` interface +- `Circle` now implements the `IHasArcLength` interface. - `Vector3.DistanceTo(Ray ray)` now returns positive infinity instead of throwing. - `Message`: removed obsolete `FromLine` method. - `AdaptiveGrid`: removed obsolete `TryGetVertexIndex` with `tolerance` parameter. diff --git a/Elements/src/Geometry/Bezier.cs b/Elements/src/Geometry/Bezier.cs index aa456041a..8d9934ce2 100644 --- a/Elements/src/Geometry/Bezier.cs +++ b/Elements/src/Geometry/Bezier.cs @@ -1,6 +1,8 @@ using System; using System.Collections.Generic; +using Elements.Geometry.Interfaces; using Newtonsoft.Json; +using System.Linq; namespace Elements.Geometry { @@ -9,7 +11,7 @@ namespace Elements.Geometry // http://webhome.cs.uvic.ca/~blob/courses/305/notes/pdf/ref-frames.pdf /// - /// The frame type to be used for operations requiring + /// The frame type to be used for operations requiring /// a moving frame around the curve. /// public enum FrameType @@ -31,7 +33,7 @@ public enum FrameType /// /// [!code-csharp[Main](../../Elements/test/BezierTests.cs?name=example)] /// - public class Bezier : BoundedCurve + public class Bezier : BoundedCurve, IHasArcLength { private readonly int _lengthSamples = 500; @@ -103,23 +105,156 @@ public override double Length() } /// - /// Calculate the length of the bezier between start and end parameters. + /// Constants for Gauss quadrature points and weights (n = 24) + /// https://pomax.github.io/bezierinfo/legendre-gauss.html /// - /// The length of the bezier between start and end. + private static readonly double[] T = new double[] + { + -0.0640568928626056260850430826247450385909, + 0.0640568928626056260850430826247450385909, + -0.1911188674736163091586398207570696318404, + 0.1911188674736163091586398207570696318404, + -0.3150426796961633743867932913198102407864, + 0.3150426796961633743867932913198102407864, + -0.4337935076260451384870842319133497124524, + 0.4337935076260451384870842319133497124524, + -0.5454214713888395356583756172183723700107, + 0.5454214713888395356583756172183723700107, + -0.6480936519369755692524957869107476266696, + 0.6480936519369755692524957869107476266696, + -0.7401241915785543642438281030999784255232, + 0.7401241915785543642438281030999784255232, + -0.8200019859739029219539498726697452080761, + 0.8200019859739029219539498726697452080761, + -0.8864155270044010342131543419821967550873, + 0.8864155270044010342131543419821967550873, + -0.9382745520027327585236490017087214496548, + 0.9382745520027327585236490017087214496548, + -0.9747285559713094981983919930081690617411, + 0.9747285559713094981983919930081690617411, + -0.9951872199970213601799974097007368118745, + 0.9951872199970213601799974097007368118745 + }; + + /// + /// Constants for Gauss quadrature weights corresponding to the points (n = 24) + /// https://pomax.github.io/bezierinfo/legendre-gauss.html + /// + private static readonly double[] C = new double[] + { + 0.1279381953467521569740561652246953718517, + 0.1279381953467521569740561652246953718517, + 0.1258374563468282961213753825111836887264, + 0.1258374563468282961213753825111836887264, + 0.121670472927803391204463153476262425607, + 0.121670472927803391204463153476262425607, + 0.1155056680537256013533444839067835598622, + 0.1155056680537256013533444839067835598622, + 0.1074442701159656347825773424466062227946, + 0.1074442701159656347825773424466062227946, + 0.0976186521041138882698806644642471544279, + 0.0976186521041138882698806644642471544279, + 0.086190161531953275917185202983742667185, + 0.086190161531953275917185202983742667185, + 0.0733464814110803057340336152531165181193, + 0.0733464814110803057340336152531165181193, + 0.0592985849154367807463677585001085845412, + 0.0592985849154367807463677585001085845412, + 0.0442774388174198061686027482113382288593, + 0.0442774388174198061686027482113382288593, + 0.0285313886289336631813078159518782864491, + 0.0285313886289336631813078159518782864491, + 0.0123412297999871995468056670700372915759, + 0.0123412297999871995468056670700372915759 + }; + + /// + /// Computes the arc length of the Bézier curve between the given parameter values start and end. + /// https://pomax.github.io/bezierinfo/#arclength + /// + /// The starting parameter value of the Bézier curve. + /// The ending parameter value of the Bézier curve. + /// The arc length between the specified parameter values. public override double ArcLength(double start, double end) { - if (!Domain.Includes(start, true)) + double z = 0.5; // Scaling factor for the Legendre-Gauss quadrature + int len = T.Length; // Number of points in the Legendre-Gauss quadrature + + double sum = 0; // Accumulated sum for the arc length calculation + + // Iterating through the Legendre-Gauss quadrature points and weights + for (int i = 0; i < len; i++) { - throw new ArgumentOutOfRangeException("start", $"The start parameter {start} must be between {Domain.Min} and {Domain.Max}."); + double t = z * T[i] + z; // Mapping the quadrature point to the Bézier parameter range [0, 1] + Vector3 derivative = Derivative(t); // Calculating the derivative of the Bézier curve at parameter t + sum += C[i] * ArcFn(t, derivative); // Adding the weighted arc length contribution to the sum } - if (!Domain.Includes(end, true)) + + // Scaling the sum by the scaling factor and the parameter interval (end - start) to get the arc length between start and end. + return z * sum * (end - start); + } + + /// + /// Calculates the arc length contribution at parameter t based on the derivative of the Bézier curve. + /// + /// The parameter value of the Bézier curve. + /// The derivative of the Bézier curve at parameter t as a Vector3. + /// The arc length contribution at parameter t. + private double ArcFn(double t, Vector3 d) + { + // Compute the Euclidean distance of the derivative vector (d) at parameter t + return Math.Sqrt(d.X * d.X + d.Y * d.Y); + } + + /// + /// Computes the derivative of the Bézier curve at parameter t. + /// https://pomax.github.io/bezierinfo/#derivatives + /// + /// The parameter value of the Bézier curve. + /// The derivative of the Bézier curve as a Vector3. + private Vector3 Derivative(double t) + { + int n = ControlPoints.Count - 1; // Degree of the Bézier curve + Vector3[] derivatives = new Vector3[n]; // Array to store the derivative control points + + // Calculating the derivative control points using the given formula + for (int i = 0; i < n; i++) { - throw new ArgumentOutOfRangeException("end", $"The end parameter {end} must be between {Domain.Min} and {Domain.Max}."); + derivatives[i] = n * (ControlPoints[i + 1] - ControlPoints[i]); } - // TODO: We use max value here so that the calculation will continue - // until at least the end of the curve. This is not a nice solution. - return ArcLengthUntil(start, double.MaxValue, out end); + // Using the derivative control points to construct an (n-1)th degree Bézier curve at parameter t. + return BezierCurveValue(t, derivatives); + } + + /// + /// Evaluates the value of an (n-1)th degree Bézier curve at parameter t using the given control points. + /// + /// The parameter value of the Bézier curve. + /// The control points for the Bézier curve. + /// The value of the Bézier curve at parameter t as a Vector3. + private Vector3 BezierCurveValue(double t, Vector3[] controlPoints) + { + int n = controlPoints.Length - 1; // Degree of the Bézier curve + Vector3[] points = new Vector3[n + 1]; + + // Initialize the points array with the provided control points + for (int i = 0; i <= n; i++) + { + points[i] = controlPoints[i]; + } + + // De Casteljau's algorithm to evaluate the value of the Bézier curve at parameter t + for (int r = 1; r <= n; r++) + { + for (int i = 0; i <= n - r; i++) + { + points[i] = (1 - t) * points[i] + t * points[i + 1]; + } + } + + // The first element of the points array contains the value of the Bézier curve at parameter t. + return points[0]; } private double ArcLengthUntil(double start, double distance, out double end) @@ -149,6 +284,15 @@ private double ArcLengthUntil(double start, double distance, out double end) return length; } + /// + /// The mid point of the curve. + /// + /// The length based midpoint. + public virtual Vector3 MidPoint() + { + return PointAtNormalizedLength(0.5); + } + /// /// Get the point on the curve at parameter u. /// @@ -168,6 +312,157 @@ public override Vector3 PointAt(double u) return p; } + /// + /// Returns the point on the bezier corresponding to the specified length value. + /// + /// The length value along the bezier. + /// The point on the bezier corresponding to the specified length value. + /// Thrown when the specified length is out of range. + public virtual Vector3 PointAtLength(double length) + { + double totalLength = ArcLength(this.Domain.Min, this.Domain.Max); // Calculate the total length of the Bezier + if (length < 0 || length > totalLength) + { + throw new ArgumentException("The specified length is out of range."); + } + return PointAt(ParameterAtDistanceFromParameter(length, Domain.Min)); + } + + /// + /// Returns the point on the bezier corresponding to the specified normalized length-based parameter value. + /// + /// The normalized length-based parameter value, ranging from 0 to 1. + /// The point on the bezier corresponding to the specified normalized length-based parameter value. + /// Thrown when the specified parameter is out of range. + public virtual Vector3 PointAtNormalizedLength(double parameter) + { + if (parameter < 0 || parameter > 1) + { + throw new ArgumentException("The specified parameter is out of range."); + } + return PointAtLength(parameter * this.ArcLength(this.Domain.Min, this.Domain.Max)); + } + + /// + /// Finds the parameter value on the Bezier curve that corresponds to the given 3D point within a specified threshold. + /// + /// The 3D point to find the corresponding parameter for. + /// The maximum distance threshold to consider a match between the projected point and the original point. + /// The parameter value on the Bezier curve if the distance between the projected point and the original point is within the threshold, otherwise returns null. + public double? ParameterAt(Vector3 point, double threshold = 0.0001) + { + // Find the parameter corresponding to the projected point on the Bezier curve + var parameter = ProjectedPoint(point, threshold); + + if (parameter == null) + { + // If the projected point does not return a relevant parameter return null + return null; + } + // Find the 3D point on the Bezier curve at the obtained parameter value + var projection = PointAt((double)parameter); + + // Check if the distance between the projected point and the original point is within + // a tolerence of the threshold + if (projection.DistanceTo(point) < (threshold * 10) - threshold) + { + // If the distance is within the threshold, return the parameter value + return parameter; + } + else + { + // If the distance exceeds the threshold, consider the point as not on the Bezier curve and return null + return null; + } + } + + /// + /// Projects a 3D point onto the Bezier curve to find the parameter value of the closest point on the curve. + /// + /// The 3D point to project onto the Bezier curve. + /// The maximum threshold to refine the projection and find the closest point. + /// The parameter value on the Bezier curve corresponding to the projected point, or null if the projection is not within the specified threshold. + public double? ProjectedPoint(Vector3 point, double threshold = 0.001) + { + // https://pomax.github.io/bezierinfo/#projections + // Generate a lookup table (LUT) of points and their corresponding parameter values on the Bezier curve + List<(Vector3 point, double t)> lut = GenerateLookupTable(); + + // Initialize variables to store the closest distance (d) and the index (index) of the closest point in the lookup table + double d = double.MaxValue; + int index = 0; + + // Find the closest point to the input point in the lookup table (LUT) using Euclidean distance + for (int i = 0; i < lut.Count; i++) + { + double q = Math.Sqrt(Math.Pow((point - lut[i].point).X, 2) + Math.Pow((point - lut[i].point).Y, 2) + Math.Pow((point - lut[i].point).Z, 2)); + if (q < d) + { + d = q; + index = i; + } + } + + // Obtain the parameter values of the neighboring points in the LUT for further refinement + double t1 = lut[Math.Max(index - 1, 0)].t; + double t2 = lut[Math.Min(index + 1, lut.Count - 1)].t; + double v = t2 - t1; + + // Refine the projection by iteratively narrowing down the parameter range to find the closest point + while (v > threshold) + { + // Calculate intermediate parameter values + double t0 = t1 + v / 4; + double t3 = t2 - v / 4; + + // Calculate corresponding points on the Bezier curve using the intermediate parameter values + Vector3 p0 = PointAt(t0); + Vector3 p3 = PointAt(t3); + + // Calculate the distances between the input point and the points on the Bezier curve + double d0 = Math.Sqrt(Math.Pow((point - p0).X, 2) + Math.Pow((point - p0).Y, 2) + Math.Pow((point - p0).Z, 2)); + double d3 = Math.Sqrt(Math.Pow((point - p3).X, 2) + Math.Pow((point - p3).Y, 2) + Math.Pow((point - p3).Z, 2)); + + // Choose the sub-range that is closer to the input point and update the range + if (d0 < d3) + { + t2 = t3; + } + else + { + t1 = t0; + } + + // Update the range difference for the next iteration + v = t2 - t1; + } + + // Return the average of the refined parameter values as the projection of the input point on the Bezier curve + return (t1 + t2) / 2; + } + + /// + /// Generates a lookup table of points and their corresponding parameter values on the Bezier curve. + /// + /// Number of samples to take along the curve. + /// A list of tuples containing the sampled points and their corresponding parameter values on the Bezier curve. + private List<(Vector3 point, double t)> GenerateLookupTable(int numSamples = 100) + { + // Initialize an empty list to store the lookup table (LUT) + List<(Vector3 point, double t)> lut = new List<(Vector3 point, double t)>(); + + // Generate lookup table by sampling points on the Bezier curve + for (int i = 0; i <= numSamples; i++) + { + double t = (double)i / numSamples; // Calculate the parameter value based on the current sample index + Vector3 point = PointAt(t); // Get the 3D point on the Bezier curve corresponding to the current parameter value + lut.Add((point, t)); // Add the sampled point and its corresponding parameter value to the lookup table (LUT) + } + + // Return the completed lookup table (LUT) + return lut; + } + private double BinomialCoefficient(int n, int i) { return Factorial(n) / (Factorial(i) * Factorial(n - i)); @@ -341,5 +636,274 @@ public override double ParameterAtDistanceFromParameter(double distance, double ArcLengthUntil(start, distance, out var end); return end; } + + /// + /// Divides the bezier into segments of the specified length. + /// + /// The desired length of each segment. + /// A list of points representing the segment divisions. + public Vector3[] DivideByLength(double divisionLength) + { + var totalLength = this.ArcLength(Domain.Min, Domain.Max); + if (totalLength <= 0) + { + // Handle invalid bezier with insufficient length + return new Vector3[0]; + } + var parameter = ParameterAtDistanceFromParameter(divisionLength, Domain.Min); + var segments = new List { this.Start }; + + while (parameter < Domain.Max) + { + segments.Add(PointAt(parameter)); + var newParameter = ParameterAtDistanceFromParameter(divisionLength, parameter); + parameter = newParameter != parameter ? newParameter : Domain.Max; + } + + // Add the last vertex of the bezier as the endpoint of the last segment if it + // is not already part of the list + if (!segments[segments.Count - 1].IsAlmostEqualTo(this.End)) + { + segments.Add(this.End); + } + + return segments.ToArray(); + } + + /// + /// Divides the bezier into segments of the specified length. + /// + /// The desired length of each segment. + /// A list of beziers representing the segments. + public List SplitByLength(double divisionLength) + { + var totalLength = this.ArcLength(Domain.Min, Domain.Max); + if (totalLength <= 0) + { + // Handle invalid bezier with insufficient length + return null; + } + var currentParameter = ParameterAtDistanceFromParameter(divisionLength, Domain.Min); + var parameters = new List { this.Domain.Min }; + + while (currentParameter < Domain.Max) + { + parameters.Add(currentParameter); + var newParameter = ParameterAtDistanceFromParameter(divisionLength, currentParameter); + currentParameter = newParameter != currentParameter ? newParameter : Domain.Max; + } + + // Add the last vertex of the bezier as the endpoint of the last segment if it + // is not already part of the list + if (!parameters[parameters.Count - 1].ApproximatelyEquals(this.Domain.Max)) + { + parameters.Add(this.Domain.Max); + } + + return Split(parameters); + } + + /// + /// Splits the Bezier curve into segments at specified parameter values. + /// + /// The list of parameter values to split the curve at. + /// If true the parameters will be length normalized. + /// A list of Bezier segments obtained after splitting. + public List Split(List parameters, bool normalize = false) + { + // Calculate the total length of the Bezier curve + var totalLength = this.ArcLength(Domain.Min, Domain.Max); + + // Check for invalid curve with insufficient length + if (totalLength <= 0) + { + throw new InvalidOperationException($"Invalid bezier with insufficient length. Total Length = {totalLength}"); + } + + // Check if the list of parameters is empty or null + if (parameters == null || parameters.Count == 0) + { + throw new ArgumentException("No split points provided."); + } + + // Initialize a list to store the resulting Bezier segments + var segments = new List(); + var bezier = this; // Create a reference to the original Bezier curve + + if (normalize) + { + parameters = parameters.Select(parameter => ParameterAtDistanceFromParameter(parameter * this.ArcLength(this.Domain.Min, this.Domain.Max), Domain.Min)).ToList(); + } + parameters.Sort(); // Sort the parameters in ascending order + + // Iterate through each parameter to split the curve + for (int i = 0; i < parameters.Count; i++) + { + double t = (Domain.Min <= parameters[i] && parameters[i] <= Domain.Max) + ? parameters[i] // Ensure the parameter is within the domain + : throw new ArgumentException($"Parameter {parameters[i]} is not within the domain ({Domain.Min}->{Domain.Max}) of the Bezier curve."); + + // Check if the parameter is within the valid range [0, 1] + if (t >= 0 && t <= 1) + { + // Split the curve at the given parameter and obtain the two resulting Bezier segments + var tuple = bezier.SplitAt(t); + + // Store the first split Bezier in the list + segments.Add(tuple.Item1); + + // Update bezier to the second split Bezier to continue splitting + bezier = tuple.Item2; + + // Remap subsequent parameters to the new Bezier curve's parameter space + for (int j = i + 1; j < parameters.Count; j++) + { + parameters[j] = (parameters[j] - t) / (1 - t); + } + } + } + + segments.Add(bezier); + + // Return the list of Bezier segments obtained after splitting + return segments; + } + + /// + /// Splits the bezier curve at the given parameter value. + /// + /// The parameter value at which to split the curve. + /// A tuple containing two split bezier curves. + public Tuple SplitAt(double t) + { + // Extract the control points from the input bezier + var startPoint = ControlPoints[0]; + var controlPoint1 = ControlPoints[1]; + var controlPoint2 = ControlPoints[2]; + var endPoint = ControlPoints[3]; + + // Compute the intermediate points using de Casteljau's algorithm + var q0 = (1 - t) * startPoint + t * controlPoint1; + var q1 = (1 - t) * controlPoint1 + t * controlPoint2; + var q2 = (1 - t) * controlPoint2 + t * endPoint; + + var r0 = (1 - t) * q0 + t * q1; + var r1 = (1 - t) * q1 + t * q2; + + // Compute the split point on the bezier curve + var splitPoint = (1 - t) * r0 + t * r1; + + // Construct the first split bezier curve + var subBezier1 = new Bezier(new List() { startPoint, q0, r0, splitPoint }); + + // Construct the second split bezier curve + var subBezier2 = new Bezier(new List() { splitPoint, r1, q2, endPoint }); + + // Return a tuple containing the split bezier curves + return new Tuple(subBezier1, subBezier2); + } + + /// + /// Constructs piecewise cubic Bézier curves from a list of points using control points calculated with the specified looseness. + /// + /// The list of points defining the path. + /// The looseness factor used to calculate control points. A higher value results in smoother curves. + /// If true, the path will be closed, connecting the last point with the first one. + /// A list of piecewise cubic Bézier curves approximating the path defined by the input points. + public static List ConstructPiecewiseCubicBezier(List points, double looseness = 6.0, bool close = false) + { + List beziers = new List(); + + // Calculate the control points. + List controlPoints = CalculateControlPoints(points, looseness); + + // Create the start Bezier curve. + Bezier startBezier = new Bezier( + new List + { + points[0], + controlPoints[0][1], + points[1] + } + ); + + // Add the start Bezier curve to the list. + beziers.Add(startBezier); + + // Iterate through pairs of points. + for (int i = 1; i < points.Count - 2; i++) + { + // Create the control points. + List bezierControlPoints = new List + { + points[i], + controlPoints[i - 1][0], + controlPoints[i][1], + points[i + 1] + }; + + // Create the Bezier curve. + Bezier bezier = new Bezier(bezierControlPoints); + + // Add the Bezier curve to the list. + beziers.Add(bezier); + } + + // Create the end Bezier curve. + Bezier endBezier = new Bezier( + new List + { + points[points.Count() - 1], + controlPoints[controlPoints.Count() - 1][0], + points[points.Count() - 2] + } + ); + + // Add the end Bezier curve to the list. + beziers.Add(endBezier); + + // Return the list of Bezier curves. + return beziers; + } + + /// + /// Calculates the control points for constructing piecewise cubic Bézier curves from the given list of points and looseness factor. + /// + /// The list of points defining the path. + /// The looseness factor used to calculate control points. A higher value results in smoother curves. + /// A list of control points (pairs of Vector3) for the piecewise cubic Bézier curves. + private static List CalculateControlPoints(List points, double looseness) + { + List controlPoints = new List(); + + for (int i = 1; i < points.Count - 1; i++) + { + // Calculate the differences in x and y coordinates. + var dx = points[i - 1].X - points[i + 1].X; + var dy = points[i - 1].Y - points[i + 1].Y; + + // Calculate the control point coordinates. + var controlPointX1 = points[i].X - dx * (1 / looseness); + var controlPointY1 = points[i].Y - dy * (1 / looseness); + var controlPoint1 = new Vector3(controlPointX1, controlPointY1, 0); + + var controlPointX2 = points[i].X + dx * (1 / looseness); + var controlPointY2 = points[i].Y + dy * (1 / looseness); + var controlPoint2 = new Vector3(controlPointX2, controlPointY2, 0); + + // Create an array to store the control points. + Vector3[] controlPointArray = new Vector3[] + { + controlPoint1, + controlPoint2 + }; + + // Add the control points to the list. + controlPoints.Add(controlPointArray); + } + + // Return the list of control points. + return controlPoints; + } } } \ No newline at end of file diff --git a/Elements/src/Geometry/Polyline.cs b/Elements/src/Geometry/Polyline.cs index 81628dbdd..89827831d 100644 --- a/Elements/src/Geometry/Polyline.cs +++ b/Elements/src/Geometry/Polyline.cs @@ -452,8 +452,6 @@ private Transform CreateOrthogonalTransform(int i, Vector3 origin, Vector3 up) /// A list of points representing the segments. public Vector3[] DivideByLength(double divisionLength) { - var segments = new List(); - if (this.Vertices.Count < 2) { // Handle invalid polyline with insufficient vertices @@ -461,7 +459,7 @@ public Vector3[] DivideByLength(double divisionLength) } var currentProgression = 0.0; - segments = new List { this.Vertices.FirstOrDefault() }; + var segments = new List { this.Vertices.FirstOrDefault() }; foreach (var currentSegment in this.Segments()) { diff --git a/Elements/test/BezierTests.cs b/Elements/test/BezierTests.cs index 359212336..feb6019d8 100644 --- a/Elements/test/BezierTests.cs +++ b/Elements/test/BezierTests.cs @@ -67,5 +67,166 @@ public void Bezier_Length_OffsetFromOrigin() var polylineLength = bezier.ToPolyline(divisions).Length(); Assert.Equal(polylineLength, bezier.Length()); } + + [Fact] + public void Bezier_ArcLength() + { + var a = new Vector3(50, 150, 0); + var b = new Vector3(105, 66, 0); + var c = new Vector3(170, 230, 0); + var d = new Vector3(200, 150, 0); + var ctrlPts = new List { a, b, c, d }; + var bezier = new Bezier(ctrlPts); + + var expectedLength = 184.38886379602502; // approximation as the integral function used for calculating length is not 100% accurate + Assert.Equal(expectedLength, bezier.ArcLength(0.0, 1.0), 2); + } + + [Fact] + public void GetParameterAt() + { + var tolerance = 0.00001; + + var a = new Vector3(1, 5, 0); + var b = new Vector3(5, 20, 0); + var c = new Vector3(5, -10, 0); + var d = new Vector3(9, 5, 0); + var ctrlPts = new List { a, b, c, d }; + var bezier = new Bezier(ctrlPts); + + var samplePt = new Vector3(0, 0, 0); + Assert.Null(bezier.ParameterAt(samplePt, tolerance)); + + samplePt = new Vector3(6.625, 0.7812, 0.0); + Assert.True((double)bezier.ParameterAt(samplePt, tolerance) - 0.75 <= tolerance * 10); + } + + [Fact] + public void GetPointAt() + { + var a = new Vector3(1, 5, 0); + var b = new Vector3(5, 20, 0); + var c = new Vector3(5, -10, 0); + var d = new Vector3(9, 5, 0); + var ctrlPts = new List { a, b, c, d }; + var bezier = new Bezier(ctrlPts); + + var testPt = new Vector3(3.3750, 9.21875, 0.0000); + Assert.True(testPt.Equals(bezier.PointAt(0.25))); + + testPt = new Vector3(4.699, 6.11375, 0.0000); + Assert.True(testPt.Equals(bezier.PointAtNormalized(0.45))); + + testPt = new Vector3(4.515904, 6.75392, 0.0000); + Assert.True(testPt.Equals(bezier.PointAtLength(8.0))); + + testPt = new Vector3(3.048823, 9.329262, 0.0000); + Assert.True(testPt.Equals(bezier.PointAtNormalizedLength(0.25))); + } + + [Fact] + public void DivideByLength() + { + var a = new Vector3(50, 150, 0); + var b = new Vector3(105, 66, 0); + var c = new Vector3(170, 230, 0); + var d = new Vector3(200, 150, 0); + var ctrlPts = new List { a, b, c, d }; + var bezier = new Bezier(ctrlPts); + + var testPts = new List(){ + new Vector3(50.00, 150.00, 0.00), + new Vector3(90.705919, 125.572992, 0.00), + new Vector3(134.607122, 148.677130, 0.00), + new Vector3(177.600231, 172.675064, 0.00), + new Vector3(200.00, 150.00, 0.00), + }; + + var ptsFromBezier = bezier.DivideByLength(50.0); + + for (int i = 0; i < ptsFromBezier.Length; i++) + { + Assert.True(ptsFromBezier[i].Equals(testPts[i])); + } + } + + [Fact] + public void Split() + { + var a = new Vector3(50, 150, 0); + var b = new Vector3(105, 66, 0); + var c = new Vector3(170, 230, 0); + var d = new Vector3(200, 150, 0); + var ctrlPts = new List { a, b, c, d }; + var bezier = new Bezier(ctrlPts); + + var testBeziers = new List(){ + new Bezier(new List() { + new Vector3(50.00, 150.00, 0.00), + new Vector3(63.75, 129.00, 0.00), + new Vector3(78.1250, 123.50, 0.00), + new Vector3(92.421875, 125.8125, 0.00) + } + ), + new Bezier(new List() { + new Vector3(92.421875, 125.8125, 0.00), + new Vector3(121.015625, 130.4375, 0.00), + new Vector3(149.296875, 166.3125, 0.00), + new Vector3(171.640625, 171.9375, 0.00) + } + ), + new Bezier(new List() { + new Vector3(171.640625, 171.9375, 0.00), + new Vector3(182.8125, 174.7500, 0.00), + new Vector3(192.50, 170.00, 0.00), + new Vector3(200.00, 150.00, 0.00) + } + ), + }; + + var beziers = bezier.Split(new List() { 0.25, 0.75 }); + + for (int i = 0; i < beziers.Count; i++) + { + Assert.True(beziers[i].ControlPoints[0].Equals(testBeziers[i].ControlPoints[0])); + Assert.True(beziers[i].ControlPoints[1].Equals(testBeziers[i].ControlPoints[1])); + Assert.True(beziers[i].ControlPoints[2].Equals(testBeziers[i].ControlPoints[2])); + Assert.True(beziers[i].ControlPoints[3].Equals(testBeziers[i].ControlPoints[3])); + } + + testBeziers = new List(){ + new Bezier(new List() { + new Vector3(50.00, 150.00, 0.00), + new Vector3(61.88, 131.856, 0.00), + new Vector3(74.22656, 125.282688, 0.00), + new Vector3(86.586183, 125.321837, 0.00) + } + ), + new Bezier(new List() { + new Vector3(86.586183, 125.321837, 0.00), + new Vector3(114.738659, 125.411011, 0.00), + new Vector3(142.958913, 159.807432, 0.00), + new Vector3(165.887648, 169.916119, 0.00) + } + ), + new Bezier(new List() { + new Vector3(165.887648, 169.916119, 0.00), + new Vector3(179.49576, 175.915584, 0.00), + new Vector3(191.24, 173.36, 0.00), + new Vector3(200.00, 150.00, 0.00) + } + ), + }; + + var normalizedBeziers = bezier.Split(new List() { 0.25, 0.75 }, true); + + for (int i = 0; i < normalizedBeziers.Count; i++) + { + Assert.True(normalizedBeziers[i].ControlPoints[0].Equals(testBeziers[i].ControlPoints[0])); + Assert.True(normalizedBeziers[i].ControlPoints[1].Equals(testBeziers[i].ControlPoints[1])); + Assert.True(normalizedBeziers[i].ControlPoints[2].Equals(testBeziers[i].ControlPoints[2])); + Assert.True(normalizedBeziers[i].ControlPoints[3].Equals(testBeziers[i].ControlPoints[3])); + } + } } } \ No newline at end of file From 8a1675ac64f395cf7119cbedfda2d1eb765f3d36 Mon Sep 17 00:00:00 2001 From: James Bradley Date: Fri, 21 Jul 2023 18:05:33 -0400 Subject: [PATCH 05/22] forgot the Z --- Elements/src/Geometry/Bezier.cs | 11 +++++++---- 1 file changed, 7 insertions(+), 4 deletions(-) diff --git a/Elements/src/Geometry/Bezier.cs b/Elements/src/Geometry/Bezier.cs index 8d9934ce2..7b0716506 100644 --- a/Elements/src/Geometry/Bezier.cs +++ b/Elements/src/Geometry/Bezier.cs @@ -881,21 +881,24 @@ private static List CalculateControlPoints(List points, doub // Calculate the differences in x and y coordinates. var dx = points[i - 1].X - points[i + 1].X; var dy = points[i - 1].Y - points[i + 1].Y; + var dz = points[i - 1].Z - points[i + 1].Z; // Calculate the control point coordinates. var controlPointX1 = points[i].X - dx * (1 / looseness); var controlPointY1 = points[i].Y - dy * (1 / looseness); - var controlPoint1 = new Vector3(controlPointX1, controlPointY1, 0); + var controlPointZ1 = points[i].Z - dz * (1 / looseness); + var controlPoint1 = new Vector3(controlPointX1, controlPointY1, controlPointZ1); var controlPointX2 = points[i].X + dx * (1 / looseness); var controlPointY2 = points[i].Y + dy * (1 / looseness); - var controlPoint2 = new Vector3(controlPointX2, controlPointY2, 0); + var controlPointZ2 = points[i].Z + dz * (1 / looseness); + var controlPoint2 = new Vector3(controlPointX2, controlPointY2, controlPointZ2); // Create an array to store the control points. Vector3[] controlPointArray = new Vector3[] { - controlPoint1, - controlPoint2 + controlPoint1, + controlPoint2 }; // Add the control points to the list. From 5324e71f91d738da7db2dd4beae206900aaed2fa Mon Sep 17 00:00:00 2001 From: James Bradley Date: Tue, 18 Jul 2023 14:49:54 -0400 Subject: [PATCH 06/22] init Polyline implementation of IHasArcLength --- CHANGELOG.md | 2 + .../src/Geometry/Interfaces/IHasArcLength.cs | 38 ++++++ Elements/src/Geometry/Polyline.cs | 109 +++++++++++++++++- Elements/test/PolylineTests.cs | 71 ++++++++++++ 4 files changed, 218 insertions(+), 2 deletions(-) create mode 100644 Elements/src/Geometry/Interfaces/IHasArcLength.cs diff --git a/CHANGELOG.md b/CHANGELOG.md index 22541dfa7..d85bc0136 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -81,6 +81,7 @@ - `Ellipse` - `EllipticalArc` - `IndexedPolycurve` +- `IHasArcLength` - `Grid1d.GetCellDomains` - `Message.Info` - `Message.Error` @@ -94,6 +95,7 @@ - `Polyline` now inherits from `BoundedCurve`. - `Polyline` is now parameterized 0->length. +- `Polyline` now implements the `IHasArcLength` interface - `Arc` now inherits from `TrimmedCurve`. - `Arc` is now parameterized 0->2Pi - `Arc` now automatically corrects decreasing angle domains to be increasing, while preserving direction. diff --git a/Elements/src/Geometry/Interfaces/IHasArcLength.cs b/Elements/src/Geometry/Interfaces/IHasArcLength.cs new file mode 100644 index 000000000..fe8a3cc90 --- /dev/null +++ b/Elements/src/Geometry/Interfaces/IHasArcLength.cs @@ -0,0 +1,38 @@ +namespace Elements.Geometry.Interfaces +{ + /// + /// Represents a curve with arc length-based operations. + /// Implementing classes define methods for computing points and performing operations based on arc length. + /// + public interface IHasArcLength + { + /// + /// Returns the point on the curve at the specified arc length. + /// + /// The arc length along the curve. + /// The point on the curve at the specified arc length. + Vector3 PointAtLength(double length); + + /// + /// Returns the point on the curve at the specified normalized length. + /// The normalized length is a value between 0 and 1 representing the relative position along the curve. + /// + /// The normalized length along the curve. + /// The point on the curve at the specified normalized length. + Vector3 PointAtNormalizedLength(double normalizedLength); + + /// + /// Returns the midpoint of the curve. + /// + /// The midpoint of the curve. + Vector3 MidPoint(); + + /// + /// Divides the curve into segments of the specified length and returns the points along the curve at those intervals. + /// + /// The desired length for dividing the curve. + /// A list of points representing the divisions along the curve. + Vector3[] DivideByLength(double length); + } + +} \ No newline at end of file diff --git a/Elements/src/Geometry/Polyline.cs b/Elements/src/Geometry/Polyline.cs index 0af338570..0357fd4ad 100644 --- a/Elements/src/Geometry/Polyline.cs +++ b/Elements/src/Geometry/Polyline.cs @@ -4,6 +4,7 @@ using System.Linq; using ClipperLib; using Elements.Search; +using Elements.Geometry.Interfaces; namespace Elements.Geometry { @@ -14,7 +15,7 @@ namespace Elements.Geometry /// /// [!code-csharp[Main](../../Elements/test/PolylineTests.cs?name=example)] /// - public class Polyline : IndexedPolycurve + public class Polyline : IndexedPolycurve, IHasArcLength { /// /// The domain of the curve. @@ -84,6 +85,66 @@ public virtual Line[] Segments() return SegmentsInternal(this.Vertices); } + + + /// + /// The mid point of the curve. + /// + /// The length based midpoint. + public virtual Vector3 MidPoint() + { + return PointAtNormalizedLength(0.5); + } + + /// + /// Returns the point on the polyline corresponding to the specified length value. + /// + /// The length value along the polyline. + /// The point on the polyline corresponding to the specified length value. + /// Thrown when the specified length is out of range. + public virtual Vector3 PointAtLength(double length) + { + double totalLength = ArcLength(this.Domain.Min, this.Domain.Max); // Calculate the total length of the Polyline + if (length < 0 || length > totalLength) + { + throw new ArgumentException("The specified length is out of range."); + } + + double accumulatedLength = 0.0; + foreach (Line segment in Segments()) + { + double segmentLength = segment.ArcLength(segment.Domain.Min, segment.Domain.Max); + + if (accumulatedLength + segmentLength >= length) + { + double remainingDistance = length - accumulatedLength; + double parameter = remainingDistance / segmentLength; + return segment.PointAtNormalized(parameter); + } + + accumulatedLength += segmentLength; + } + + // If we reach here, the desired length is equal to the total length, + // so return the end point of the Polyline. + return End; + } + + /// + /// Returns the point on the polyline corresponding to the specified normalized length-based parameter value. + /// + /// The normalized length-based parameter value, ranging from 0 to 1. + /// The point on the polyline corresponding to the specified normalized length-based parameter value. + /// Thrown when the specified parameter is out of range. + public virtual Vector3 PointAtNormalizedLength(double parameter) + { + if (parameter < 0 || parameter > 1) + { + throw new ArgumentException("The specified parameter is out of range."); + } + return PointAtLength(parameter * this.ArcLength(this.Domain.Min, this.Domain.Max)); + } + /// /// Get the transform at the specified parameter along the polyline. /// @@ -422,6 +483,50 @@ private Transform CreateOrthogonalTransform(int i, Vector3 origin, Vector3 up) return new Transform(origin, up.Cross(tangent), tangent); } + /// + /// Divides the polyline into segments of the specified length. + /// + /// The desired length of each segment. + /// A list of points representing the segments. + public Vector3[] DivideByLength(double divisionLength) + { + var segments = new List(); + + if (this.Vertices.Count < 2) + { + // Handle invalid polyline with insufficient vertices + return new Vector3[0]; + } + + var currentProgression = 0.0; + segments = new List { this.Vertices.FirstOrDefault() }; + + foreach (var currentSegment in this.Segments()) + { + // currentProgression from last segment before hitting end + if (currentProgression != 0.0) + { + currentProgression -= divisionLength; + } + while (currentSegment.ArcLength(currentSegment.Domain.Min, currentSegment.Domain.Max) >= currentProgression + divisionLength) + { + segments.Add(currentSegment.PointAt(currentProgression + divisionLength)); + currentProgression += divisionLength; + } + // Set currentProgression from divisionLength less distance from last segment point + currentProgression = divisionLength - segments.LastOrDefault().DistanceTo(currentSegment.End); + } + + // Add the last vertex of the polyline as the endpoint of the last segment if it + // is not already part of the list + if (!segments.LastOrDefault().IsAlmostEqualTo(this.Vertices.LastOrDefault())) + { + segments.Add(this.Vertices.LastOrDefault()); + } + + return segments.ToArray(); + } + /// /// Offset this polyline by the specified amount. /// @@ -831,7 +936,7 @@ protected void Split(IList points, bool closed = false) var b = closed && i == this.Vertices.Count - 1 ? this.Vertices[0] : this.Vertices[i + 1]; var edge = (a, b); - // An edge may have multiple split points. + // An edge may have multiple split points. // We store these in a list and sort it along the // direction of the edge, before inserting the points // into the vertex list and incrementing i by the correct diff --git a/Elements/test/PolylineTests.cs b/Elements/test/PolylineTests.cs index 645c406f1..5bf12857b 100644 --- a/Elements/test/PolylineTests.cs +++ b/Elements/test/PolylineTests.cs @@ -257,6 +257,77 @@ public void GetParameterAt() Assert.Equal(1.5, resultParameter); var testPoint = polyline.PointAt(resultParameter); Assert.True(point.IsAlmostEqualTo(testPoint)); + + polyline = new Polyline( + new List() + { + new Vector3(1, 5, 0), + new Vector3(5, 20, 0), + new Vector3(5, 0, 0), + new Vector3(9, 5, 0), + } + ); + var parameter = 0.5; + var testMidpoint = new Vector3(5.0, 14.560525, 0); // Midpoint + var testParameterMidpoint = polyline.PointAtNormalizedLength(parameter); + var t = testParameterMidpoint.IsAlmostEqualTo(testMidpoint); + Assert.True(testParameterMidpoint.IsAlmostEqualTo(testMidpoint)); + + var midlength = polyline.Length() * parameter; + var testLengthMidpoint = polyline.PointAtLength(midlength); + Assert.True(testLengthMidpoint.IsAlmostEqualTo(testParameterMidpoint)); + + var midpoint = polyline.MidPoint(); + Assert.True(midpoint.IsAlmostEqualTo(testLengthMidpoint)); + } + + [Fact] + public void DivideByLength_SingleSegment() + { + // Arrange + var start = new Vector3(1, 2); + var end = new Vector3(1, 4); + var polyline = new Polyline(new[] { start, end }); + var divisionLength = 1.0; + + // Act + var segments = polyline.DivideByLength(divisionLength); + + // Assert + Assert.Equal(3, segments.Count()); + Assert.Equal(start, segments[0]); + Assert.Equal(new Vector3(1, 3), segments[1]); + Assert.Equal(end, segments[2]); + } + + [Fact] + public void DivideByLength_MultipleSegments() + { + // Arrange + var polyline = new Polyline(new List() + { + new Vector3(1, 5, 0), + new Vector3(5, 20, 0), + new Vector3(5, 0, 0), + new Vector3(9, 5, 0), + }); + var divisionLength = 5.0; + + // Act + var segments = polyline.DivideByLength(divisionLength); + + // Assert + Assert.Equal(10, segments.Count()); + Assert.Equal(new Vector3(1.0, 5.0, 0), segments[0]); + Assert.Equal(new Vector3(2.288313, 9.831174, 0), segments[1]); + Assert.Equal(new Vector3(3.576626, 14.662349, 0), segments[2]); + Assert.Equal(new Vector3(4.864939, 19.493524, 0), segments[3]); + Assert.Equal(new Vector3(5.0, 15.524174, 0), segments[4]); + Assert.Equal(new Vector3(5.0, 10.524174, 0), segments[5]); + Assert.Equal(new Vector3(5.0, 5.524174, 0), segments[6]); + Assert.Equal(new Vector3(5.0, 0.524174, 0), segments[7]); + Assert.Equal(new Vector3(7.796025, 3.495032, 0), segments[8]); + Assert.Equal(new Vector3(9, 5, 0), segments[9]); } [Fact] From af5f18d7f79eb6508425082305714d38d9dac59c Mon Sep 17 00:00:00 2001 From: James Bradley Date: Tue, 18 Jul 2023 17:36:40 -0400 Subject: [PATCH 07/22] Line implements Interface IHasArcLength --- CHANGELOG.md | 1 + Elements/src/Elements.csproj | 2 +- Elements/src/Geometry/Line.cs | 130 ++++++++++++++++++++++++------ Elements/src/Geometry/Polyline.cs | 2 - Elements/test/LineTests.cs | 21 +++-- Elements/test/PolylineTests.cs | 1 - 6 files changed, 125 insertions(+), 32 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index d85bc0136..c3ffa80bb 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -101,6 +101,7 @@ - `Arc` now automatically corrects decreasing angle domains to be increasing, while preserving direction. - `Line` now inherits from `TrimmedCurve`. - `Line` is now parameterized 0->length. +- `Line` now implements the `IHasArcLength` interface - `Bezier` now inherits from `BoundedCurve`. - `Polyline` is now parameterized 0->length. - `Circle` is now parameterized 0->2Pi. diff --git a/Elements/src/Elements.csproj b/Elements/src/Elements.csproj index cdd21865a..8bc4905ea 100644 --- a/Elements/src/Elements.csproj +++ b/Elements/src/Elements.csproj @@ -10,7 +10,7 @@ The Elements library provides object types for generating the built environment. https://github.com/hypar-io/elements https://github.com/hypar-io/elements - $(Version) + 21.21.21 en diff --git a/Elements/src/Geometry/Line.cs b/Elements/src/Geometry/Line.cs index 068292ef5..d23279878 100644 --- a/Elements/src/Geometry/Line.cs +++ b/Elements/src/Geometry/Line.cs @@ -4,6 +4,7 @@ using System.Linq; using Newtonsoft.Json; using Elements.Spatial; +using Elements.Geometry.Interfaces; namespace Elements.Geometry { @@ -15,7 +16,7 @@ namespace Elements.Geometry /// [!code-csharp[Main](../../Elements/test/LineTests.cs?name=example)] /// /// TODO: Rename this class to LineSegment - public class Line : TrimmedCurve, IEquatable + public class Line : TrimmedCurve, IEquatable, IHasArcLength { /// /// The domain of the curve. @@ -123,6 +124,48 @@ public override Vector3 PointAt(double u) return this.BasisCurve.PointAt(u); } + /// + /// The mid point of the curve. + /// + /// The length based midpoint. + public virtual Vector3 MidPoint() + { + return PointAtNormalizedLength(0.5); + } + + /// + /// Returns the point on the line corresponding to the specified length value. + /// + /// The length value along the line. + /// The point on the line corresponding to the specified length value. + /// Thrown when the specified length is out of range. + public virtual Vector3 PointAtLength(double length) + { + double totalLength = ArcLength(this.Domain.Min, this.Domain.Max); // Calculate the total length of the Line + + if (length < 0 || length > totalLength) + { + throw new ArgumentException("The specified length is out of range."); + } + var lengthParameter = length / totalLength; + return this.PointAtNormalized(lengthParameter); + } + + /// + /// Returns the point on the line corresponding to the specified normalized length-based parameter value. + /// + /// The normalized length-based parameter value, ranging from 0 to 1. + /// The point on the line corresponding to the specified normalized length-based parameter value. + /// Thrown when the specified parameter is out of range. + public virtual Vector3 PointAtNormalizedLength(double parameter) + { + if (parameter < 0 || parameter > 1) + { + throw new ArgumentException("The specified parameter is out of range."); + } + return PointAtLength(parameter * this.ArcLength(this.Domain.Min, this.Domain.Max)); + } + /// public override Curve Transformed(Transform transform) { @@ -553,38 +596,79 @@ public static bool PointOnLine(Vector3 point, Vector3 start, Vector3 end, bool i return false; } + // /// + // /// Divide the line into as many segments of the provided length as possible. + // /// + // /// The length. + // /// A flag indicating whether segments shorter than l should be removed. + // public List DivideByLength(double l, bool removeShortSegments = false) + // { + // var len = this.Length(); + // if (l > len) + // { + // return new List() { new Line(this.Start, this.End) }; + // } + + // var total = 0.0; + // var d = this.Direction(); + // var lines = new List(); + // while (total + l <= len) + // { + // var a = this.Start + d * total; + // var b = a + d * l; + // lines.Add(new Line(a, b)); + // total += l; + // } + // if (total < len && !removeShortSegments) + // { + // var a = this.Start + d * total; + // if (!a.IsAlmostEqualTo(End)) + // { + // lines.Add(new Line(a, End)); + // } + // } + // return lines; + // } + /// - /// Divide the line into as many segments of the provided length as possible. + /// Divides the line into segments of the specified length. /// - /// The length. - /// A flag indicating whether segments shorter than l should be removed. - public List DivideByLength(double l, bool removeShortSegments = false) + /// The desired length of each segment. + /// A list of points representing the segments. + public Vector3[] DivideByLength(double divisionLength) { - var len = this.Length(); - if (l > len) + var segments = new List(); + + if (this.ArcLength(this.Domain.Min, this.Domain.Max) < double.Epsilon) { - return new List() { new Line(this.Start, this.End) }; + // Handle invalid line with insufficient length + return new Vector3[0]; } - var total = 0.0; - var d = this.Direction(); - var lines = new List(); - while (total + l <= len) + var currentProgression = 0.0; + segments = new List { this.Start }; + + // currentProgression from last segment before hitting end + if (currentProgression != 0.0) { - var a = this.Start + d * total; - var b = a + d * l; - lines.Add(new Line(a, b)); - total += l; + currentProgression -= divisionLength; } - if (total < len && !removeShortSegments) + while (this.ArcLength(this.Domain.Min, this.Domain.Max) >= currentProgression + divisionLength) { - var a = this.Start + d * total; - if (!a.IsAlmostEqualTo(End)) - { - lines.Add(new Line(a, End)); - } + segments.Add(this.PointAt(currentProgression + divisionLength)); + currentProgression += divisionLength; } - return lines; + // Set currentProgression from divisionLength less distance from last segment point + currentProgression = divisionLength - segments.LastOrDefault().DistanceTo(this.End); + + // Add the last vertex of the polyline as the endpoint of the last segment if it + // is not already part of the list + if (!segments.LastOrDefault().IsAlmostEqualTo(this.End)) + { + segments.Add(this.End); + } + + return segments.ToArray(); } /// diff --git a/Elements/src/Geometry/Polyline.cs b/Elements/src/Geometry/Polyline.cs index 0357fd4ad..cdfedc343 100644 --- a/Elements/src/Geometry/Polyline.cs +++ b/Elements/src/Geometry/Polyline.cs @@ -85,8 +85,6 @@ public virtual Line[] Segments() return SegmentsInternal(this.Vertices); } - - /// /// The mid point of the curve. /// diff --git a/Elements/test/LineTests.cs b/Elements/test/LineTests.cs index 507c7b61c..4d1888975 100644 --- a/Elements/test/LineTests.cs +++ b/Elements/test/LineTests.cs @@ -262,11 +262,11 @@ public void DivideIntoEqualSegmentsSingle() public void DivideByLength() { var l = new Line(Vector3.Origin, new Vector3(5, 0)); - var segments = l.DivideByLength(1.1, true); - Assert.Equal(4, segments.Count); + var segments = l.DivideByLength(1.1); + Assert.Equal(6, segments.Count()); - var segments1 = l.DivideByLength(1.1); - Assert.Equal(5, segments1.Count); + var segments1 = l.DivideByLength(2); + Assert.Equal(4, segments1.Count()); } [Fact] @@ -749,6 +749,17 @@ public void GetParameterAt() var expectedVector = line.PointAt(uValue); Assert.InRange(uValue, 0, line.Length()); Assert.True(vector.IsAlmostEqualTo(expectedVector)); + + var parameter = 0.5; + var testParameterMidpoint = line.PointAtNormalizedLength(parameter); + Assert.True(testParameterMidpoint.IsAlmostEqualTo(middle)); + + var midlength = line.Length() * parameter; + var testLengthMidpoint = line.PointAtLength(midlength); + Assert.True(testLengthMidpoint.IsAlmostEqualTo(testParameterMidpoint)); + + var midpoint = line.MidPoint(); + Assert.True(midpoint.IsAlmostEqualTo(testLengthMidpoint)); } [Theory] @@ -1209,7 +1220,7 @@ public void LineDistancePointsOnSkewLines() Assert.Equal(delta.Length(), (new Line(pt12, pt11)).DistanceTo(new Line(pt21, pt22)), 12); Assert.Equal(delta.Length(), (new Line(pt12, pt11)).DistanceTo(new Line(pt22, pt21)), 12); //The segments (pt12, pt13) and (pt21, pt22) does not intersect. - //The shortest distance is from an endpoint to another segment - difference between lines plus between endpoints. + //The shortest distance is from an endpoint to another segment - difference between lines plus between endpoints. var expected = (q12 * v1).DistanceTo(new Line(delta + q21 * v2, delta + q22 * v2)); Assert.Equal(expected, (new Line(pt12, pt13)).DistanceTo(new Line(pt21, pt22)), 12); Assert.Equal(expected, (new Line(pt12, pt13)).DistanceTo(new Line(pt22, pt21)), 12); diff --git a/Elements/test/PolylineTests.cs b/Elements/test/PolylineTests.cs index 5bf12857b..ebf69e91a 100644 --- a/Elements/test/PolylineTests.cs +++ b/Elements/test/PolylineTests.cs @@ -270,7 +270,6 @@ public void GetParameterAt() var parameter = 0.5; var testMidpoint = new Vector3(5.0, 14.560525, 0); // Midpoint var testParameterMidpoint = polyline.PointAtNormalizedLength(parameter); - var t = testParameterMidpoint.IsAlmostEqualTo(testMidpoint); Assert.True(testParameterMidpoint.IsAlmostEqualTo(testMidpoint)); var midlength = polyline.Length() * parameter; From 4edac7aa8e2866c0961092a71b36c1c8bce93edc Mon Sep 17 00:00:00 2001 From: James Bradley Date: Tue, 18 Jul 2023 21:41:01 -0400 Subject: [PATCH 08/22] Circle implements Interface IHasArcLength --- CHANGELOG.md | 2 +- Elements/src/Geometry/Circle.cs | 285 ++++++++++++-------------------- Elements/test/CircleTests.cs | 127 +++++++++++++- Elements/test/LineTests.cs | 1 - 4 files changed, 230 insertions(+), 185 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index c3ffa80bb..a32f60829 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -105,7 +105,7 @@ - `Bezier` now inherits from `BoundedCurve`. - `Polyline` is now parameterized 0->length. - `Circle` is now parameterized 0->2Pi. -- `Line` is now parameterized 0->length. +- `Circle` now implements the `IHasArcLength` interface - `Vector3.DistanceTo(Ray ray)` now returns positive infinity instead of throwing. - `Message`: removed obsolete `FromLine` method. - `AdaptiveGrid`: removed obsolete `TryGetVertexIndex` with `tolerance` parameter. diff --git a/Elements/src/Geometry/Circle.cs b/Elements/src/Geometry/Circle.cs index 69296293e..def74117b 100644 --- a/Elements/src/Geometry/Circle.cs +++ b/Elements/src/Geometry/Circle.cs @@ -1,3 +1,4 @@ +using Elements.Validators; using System; using System.Collections.Generic; using System.Linq; @@ -7,10 +8,10 @@ namespace Elements.Geometry { /// - /// A circle. + /// A circle. /// Parameterization of the circle is 0 -> 2PI. /// - public class Circle : Curve, IConic + public class Circle : Curve, IConic, IHasArcLength { /// The center of the circle. [JsonProperty("Center", Required = Required.AllowNull)] @@ -27,6 +28,17 @@ public Vector3 Center [System.ComponentModel.DataAnnotations.Range(0.0D, double.MaxValue)] public double Radius { get; protected set; } + /// The circumference of the circle. + [JsonIgnore] + [System.ComponentModel.DataAnnotations.Range(0.0D, double.MaxValue)] + public double Circumference { get; protected set; } + + /// + /// The domain of the curve. + /// + [JsonIgnore] + public Domain1d Domain => new Domain1d(0, 2 * Math.PI); + /// /// The coordinate system of the plane containing the circle. /// @@ -50,7 +62,15 @@ public Vector3 Normal [JsonConstructor] public Circle(Vector3 center, double radius = 1.0) { + if (!Validator.DisableValidationOnConstruction) + { + if (Math.Abs(radius - 0.0) < double.Epsilon ? true : false) + { + throw new ArgumentException($"The circle could not be created. The radius of the circle cannot be the zero: radius {radius}"); + } + } this.Radius = radius; + this.Circumference = 2 * Math.PI * this.Radius; this.Transform = new Transform(center); } @@ -60,7 +80,15 @@ public Circle(Vector3 center, double radius = 1.0) /// The radius of the circle. public Circle(double radius = 1.0) { + if (!Validator.DisableValidationOnConstruction) + { + if (Math.Abs(radius - 0.0) < double.Epsilon ? true : false) + { + throw new ArgumentException($"The circle could not be created. The radius of the circle cannot be the zero: radius {radius}"); + } + } this.Radius = radius; + this.Circumference = 2 * Math.PI * this.Radius; this.Transform = new Transform(); } @@ -69,8 +97,65 @@ public Circle(double radius = 1.0) /// public Circle(Transform transform, double radius = 1.0) { + if (!Validator.DisableValidationOnConstruction) + { + if (Math.Abs(radius - 0.0) < double.Epsilon ? true : false) + { + throw new ArgumentException($"The circle could not be created. The radius of the circle cannot be the zero: radius {radius}"); + } + } this.Transform = transform; this.Radius = radius; + this.Circumference = 2 * Math.PI * this.Radius; + } + + /// + /// Calculate the length of the circle between two parameters. + /// + public double ArcLength(double start, double end) + { + // Convert start and end parameters from radians to degrees + double _startAngle = start * 180.0 / Math.PI; + double _endAngle = end * 180.0 / Math.PI; + + // Ensure the start angle is within the valid domain range of 0 to 360 degrees + double startAngle = _startAngle % 360; + if (startAngle < 0) + { + startAngle += 360; + } + + // Ensure the end angle is within the valid domain range of 0 to 360 degrees + double endAngle = _endAngle % 360; + if (endAngle < 0) + { + endAngle += 360; + } + else if (endAngle == 0 && Math.Abs(_endAngle) >= 2 * Math.PI) + { + endAngle = 360; + } + + // Calculate the difference in angles + double angleDifference = endAngle - startAngle; + + // Adjust the angle difference if it crosses the 360-degree boundary + if (angleDifference < 0) + { + angleDifference += 360; + } + else if (angleDifference >= 2 * Math.PI) + { + return Circumference; // Full circle, return circumference + } + + // Convert the angle difference back to radians + double angleDifferenceRadians = angleDifference * Math.PI / 180.0; + + // Calculate the arc length using the formula: arc length = radius * angle + double arcLength = Radius * angleDifferenceRadians; + + return arcLength; } /// @@ -94,6 +179,14 @@ public Polygon ToPolygon(int divisions = 10) return new Polygon(pts, true); } + /// + /// Are the two circles almost equal? + /// + public bool IsAlmostEqualTo(Circle other, double tolerance = Vector3.EPSILON) + { + return (Center.IsAlmostEqualTo(other.Center, tolerance) && Math.Abs(Radius - other.Radius) < tolerance ? true : false); + } + /// /// Convert a circle to a circular arc. /// @@ -105,6 +198,15 @@ public Polygon ToPolygon(int divisions = 10) /// The bounded curve to convert. public static implicit operator ModelCurve(Circle c) => new ModelCurve(c); + /// + /// Calculates and returns the midpoint of the circle. + /// + /// The midpoint of the circle. + public Vector3 MidPoint() + { + return PointAt(Math.PI); + } + /// /// Return the point at parameter u on the arc. /// @@ -122,33 +224,6 @@ private Vector3 PointAtUntransformed(double u) return new Vector3(x, y); } - /// - /// Check if certain point is on the circle. - /// - /// Point to check. - /// Calculated parameter of point on circle. - /// True if point lays on the circle. - public bool ParameterAt(Vector3 pt, out double t) - { - var local = Transform.Inverted().OfPoint(pt); - if (local.Z.ApproximatelyEquals(0) && - local.LengthSquared().ApproximatelyEquals( - Radius * Radius, Vector3.EPSILON * Vector3.EPSILON)) - { - t = ParameterAtUntransformed(local); - return true; - } - - t = 0; - return false; - } - - private double ParameterAtUntransformed(Vector3 pt) - { - var v = pt / Radius; - return Math.Atan2(v.Y, v.X); - } - /// /// Return transform on the arc at parameter u. /// @@ -186,157 +261,5 @@ public override double ParameterAtDistanceFromParameter(double distance, double return start + theta; } - - /// - public override bool Intersects(ICurve curve, out List results) - { - switch (curve) - { - case BoundedCurve boundedCurve: - return boundedCurve.Intersects(this, out results); - case InfiniteLine line : - return Intersects(line, out results); - case Circle circle: - return Intersects(circle, out results); - case Ellipse elliplse: - return Intersects(elliplse, out results); - default: - throw new NotImplementedException(); - } - } - - /// - /// Does this circle intersects with other circle? - /// Circles with the same positions and radii are not considered as intersecting. - /// - /// Other circle to intersect. - /// List containing up to two intersection points. - /// True if any intersections exist, otherwise false. - public bool Intersects(Circle other, out List results) - { - results = new List(); - - Plane planeA = new Plane(Center, Normal); - Plane planeB = new Plane(other.Center, other.Normal); - - // Check if two circles are on the same plane. - if (Normal.IsParallelTo(other.Normal, Vector3.EPSILON * Vector3.EPSILON) && - other.Center.DistanceTo(planeA).ApproximatelyEquals(0)) - { - var delta = other.Center - Center; - var dist = delta.Length(); - // Check if circles are on correct distance for intersection to happen. - if (dist.ApproximatelyEquals(0) || - dist > Radius + other.Radius || dist < Math.Abs(Radius - other.Radius)) - { - return false; - } - - // Build triangle with center of one circle and two intersection points. - var r1squre = Radius * Radius; - var r2squre = other.Radius * other.Radius; - var lineDist = (r1squre - r2squre + dist * dist) / (2 * dist); - var linePoint = Center + lineDist * delta.Unitized(); - double perpDistance = Math.Sqrt(r1squre - lineDist * lineDist); - // If triangle side is 0 - circles touches. Only one intersection recorded. - if (perpDistance.ApproximatelyEquals(0)) - { - results.Add(linePoint); - } - else - { - Vector3 perpDirection = delta.Cross(Normal).Unitized(); - results.Add(linePoint + perpDirection * perpDistance); - results.Add(linePoint - perpDirection * perpDistance); - } - } - // Ignore circles on parallel planes. - // Find intersection line between two planes. - else if (planeA.Intersects(planeB, out var line) && - Intersects(line, out var candidates)) - { - foreach (var item in candidates) - { - // Check each point that lays on intersection line and one of the circles. - // They are on both if they have correct distance to circle centers. - if (item.DistanceTo(other.Center).ApproximatelyEquals(other.Radius)) - { - results.Add(item); - } - } - } - - return results.Any(); - } - - /// - /// Does this circle intersects with an infinite line? - /// - /// Infinite line to intersect. - /// List containing up to two intersection points. - /// True if any intersections exist, otherwise false. - public bool Intersects(InfiniteLine line, out List results) - { - results = new List(); - - Plane circlePlane = new Plane(Center, Normal); - Vector3 closestPoint; - bool lineOnPlane = line.Origin.DistanceTo(circlePlane).ApproximatelyEquals(0) && - line.Direction.Dot(Normal).ApproximatelyEquals(0); - - // If line share a plane with circle - find closest point on it to circle center. - // If not - check if there an intersection between line and circle plane. - if (lineOnPlane) - { - closestPoint = Center.ClosestPointOn(line); - } - else if (!line.Intersects(circlePlane, out closestPoint)) - { - return false; - } - - var delta = closestPoint - Center; - var lengthSquared = delta.LengthSquared(); - var radiusSquared = Radius * Radius; - var toleranceSquared = Vector3.EPSILON * Vector3.EPSILON; - // if line not on circle plane - only one intersection is possible if it's radius away. - // this will also happen if line is on plane but only touches the circle. - if (lengthSquared.ApproximatelyEquals(radiusSquared, toleranceSquared)) - { - results.Add(closestPoint); - } - else if (lineOnPlane && lengthSquared < radiusSquared) - { - var distance = Math.Sqrt(radiusSquared - lengthSquared); - results.Add(closestPoint + line.Direction * distance); - results.Add(closestPoint - line.Direction * distance); - } - - return results.Any(); - } - - /// - /// Does this circle intersects with an ellipse? - /// Circle and ellipse that are coincides are not considered as intersecting. - /// - /// - /// Ellipse to intersect. - /// List containing up to four intersection points. - /// True if any intersections exist, otherwise false. - public bool Intersects(Ellipse ellipse, out List results) - { - return ellipse.Intersects(this, out results); - } - - /// - /// Does this circle intersects with a bounded curve? - /// - /// Curve to intersect. - /// List containing intersection points. - /// True if any intersections exist, otherwise false. - public bool Intersects(BoundedCurve curve, out List results) - { - return curve.Intersects(this, out results); - } } } \ No newline at end of file diff --git a/Elements/test/CircleTests.cs b/Elements/test/CircleTests.cs index 5be213dde..4eb62dfda 100644 --- a/Elements/test/CircleTests.cs +++ b/Elements/test/CircleTests.cs @@ -1,4 +1,4 @@ -using Elements; +using Elements; using Elements.Geometry; using Elements.Tests; using System; @@ -14,7 +14,7 @@ public class CircleTests [Fact] public void CirceIntersectsCircle() { - // Planar intersecting circles + // Planar intersecting circles Circle c0 = new Circle(5); Circle c1 = new Circle(new Vector3(8, 0, 0), 5); Assert.True(c0.Intersects(c1, out var results)); @@ -135,5 +135,128 @@ public void AdjustRadian() Assert.Equal(Units.DegreesToRadians(60), Units.AdjustRadian(Units.DegreesToRadians(420), reference), 6); Assert.Equal(Units.DegreesToRadians(-100), Units.AdjustRadian(Units.DegreesToRadians(-100), reference), 6); } + + [Fact, Trait("Category", "Examples")] + public void CircleExample() + { + this.Name = "Elements_Geometry_Circle"; + + // + var a = new Vector3(); + var b = 1.0; + var c = new Circle(a, b); + // + + this.Model.AddElement(c); + } + + [Fact] + public void Equality() + { + var p = 1.0; + var circleA = new Circle(Vector3.Origin, p); + var circleB = new Circle(Vector3.Origin, p + 1E-4); + var circleC = new Circle(Vector3.Origin, p + 1E-6); + + Assert.False(circleA.IsAlmostEqualTo(circleB)); + Assert.True(circleA.IsAlmostEqualTo(circleB, 1E-3)); + Assert.True(circleA.IsAlmostEqualTo(circleC)); + } + + [Fact] + public void Construct() + { + var a = new Vector3(); + var b = 1.0; + var c = new Circle(a, b); + Assert.Equal(1.0, c.Radius); + Assert.Equal(new Vector3(0, 0), c.Center); + Assert.Equal(a + new Vector3(b, 0, 0), c.PointAt(-1e-10)); + } + + [Fact] + public void ZeroRadius_ThrowsException() + { + var a = new Vector3(); + Assert.Throws(() => new Circle(a, 0)); + } + + [Fact] + public void DivideIntoEqualSegments() + { + var l = new Line(Vector3.Origin, new Vector3(100, 0)); + var segments = l.DivideIntoEqualSegments(41); + var len = l.Length(); + Assert.Equal(41, segments.Count); + foreach (var s in segments) + { + Assert.Equal(s.Length(), len / 41, 5); + } + } + + [Fact] + public void DivideIntoEqualSegmentsSingle() + { + var l = new Line(Vector3.Origin, new Vector3(100, 0)); + var segments = l.DivideIntoEqualSegments(1); + Assert.Single(segments); + Assert.True(segments.First().Start.IsAlmostEqualTo(l.Start, 1e-10)); + Assert.True(segments.First().End.IsAlmostEqualTo(l.End, 1e-10)); + } + + [Fact] + public void DivideByLength() + { + var l = new Line(Vector3.Origin, new Vector3(5, 0)); + var segments = l.DivideByLength(1.1); + Assert.Equal(6, segments.Count()); + + var segments1 = l.DivideByLength(2); + Assert.Equal(4, segments1.Count()); + } + + [Fact] + public void GetParameterAt() + { + var center = Vector3.Origin; + var radius = 5.0; + var circle = new Circle(center, radius); + var start = new Vector3(5.0, 0.0, 0.0); + var mid = new Vector3(-5.0, 0.0, 0.0); + + Assert.Equal(0, circle.GetParameterAt(start)); + + var almostEqualStart = new Vector3(5.000001, 0.000005, 0); + Assert.True(start.IsAlmostEqualTo(almostEqualStart)); + + Assert.True(Math.Abs(circle.GetParameterAt(mid) - Math.PI) < double.Epsilon ? true : false); + + var vector = new Vector3(3.535533, 3.535533, 0.0); + var uValue = circle.GetParameterAt(vector); + var expectedVector = circle.PointAt(uValue); + Assert.InRange(uValue, 0, circle.Circumference); + Assert.True(vector.IsAlmostEqualTo(expectedVector)); + + var parameter = 0.5; + var testParameterMidpoint = circle.PointAtNormalizedLength(parameter); + Assert.True(testParameterMidpoint.IsAlmostEqualTo(mid)); + + var midlength = circle.Circumference * parameter; + var testLengthMidpoint = circle.PointAtLength(midlength); + Assert.True(testLengthMidpoint.IsAlmostEqualTo(testParameterMidpoint)); + + var midpoint = circle.MidPoint(); + Assert.True(midpoint.IsAlmostEqualTo(testLengthMidpoint)); + } + + [Fact] + public void PointOnLine() + { + Circle circle = new Circle(Vector3.Origin, 5.0); + + Assert.False(Circle.PointOnCircle(Vector3.Origin, circle)); + Assert.False(Circle.PointOnCircle(new Vector3(4, 0, 0), circle)); + Assert.True(Circle.PointOnCircle(circle.PointAtNormalizedLength(0.5), circle)); + } } } diff --git a/Elements/test/LineTests.cs b/Elements/test/LineTests.cs index 4d1888975..1a9fad806 100644 --- a/Elements/test/LineTests.cs +++ b/Elements/test/LineTests.cs @@ -717,7 +717,6 @@ public void TryGetOverlap() Assert.Equal(new Line((5, 0, 0), (10, 0, 0)), overlap); } - [Fact] public void GetParameterAt() { From 67f1a8fb3a9476c25ef2e4b024e16d22eec7ce8e Mon Sep 17 00:00:00 2001 From: James Bradley Date: Fri, 21 Jul 2023 15:46:54 -0400 Subject: [PATCH 09/22] Polish Bezier, IHasCurveLength --- CHANGELOG.md | 17 +- Elements/src/Geometry/Bezier.cs | 625 +++++++++++++++--------------- Elements/src/Geometry/Polyline.cs | 4 +- Elements/test/BezierTests.cs | 165 +++++++- 4 files changed, 482 insertions(+), 329 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index a32f60829..067fe8328 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -90,22 +90,31 @@ - `new Topography(Topography other)` - `Topography.TopMesh()` - `UpdateElementRepresentations` flag to all serialization methods +- `Bezier.PointAtLength()` +- `Bezier.PointAtNormalizedLength()` +- `Bezier.ParameterAt()` +- `Bezier.DivideByLength()` +- `Bezier.Split()` +- `Bezier.SplitAt()` +- `Bezier.SplitByLength()` +- `Bezier.ConstructPiecewiseCubicBezier()` ### Changed - `Polyline` now inherits from `BoundedCurve`. - `Polyline` is now parameterized 0->length. -- `Polyline` now implements the `IHasArcLength` interface +- `Polyline` now implements the `IHasArcLength` interface. - `Arc` now inherits from `TrimmedCurve`. - `Arc` is now parameterized 0->2Pi -- `Arc` now automatically corrects decreasing angle domains to be increasing, while preserving direction. - `Line` now inherits from `TrimmedCurve`. - `Line` is now parameterized 0->length. -- `Line` now implements the `IHasArcLength` interface +- `Line` now implements the `IHasArcLength` interface. - `Bezier` now inherits from `BoundedCurve`. +- `Bezier` now implements the `IHasArcLength` interface. +- `Bezier.ArcLength()` now uses Gauss quadrature approximation vs linear sampling. - `Polyline` is now parameterized 0->length. - `Circle` is now parameterized 0->2Pi. -- `Circle` now implements the `IHasArcLength` interface +- `Circle` now implements the `IHasArcLength` interface. - `Vector3.DistanceTo(Ray ray)` now returns positive infinity instead of throwing. - `Message`: removed obsolete `FromLine` method. - `AdaptiveGrid`: removed obsolete `TryGetVertexIndex` with `tolerance` parameter. diff --git a/Elements/src/Geometry/Bezier.cs b/Elements/src/Geometry/Bezier.cs index c9fbad96a..af7248b72 100644 --- a/Elements/src/Geometry/Bezier.cs +++ b/Elements/src/Geometry/Bezier.cs @@ -3,6 +3,7 @@ using System.Linq; using Elements.Geometry.Interfaces; using Newtonsoft.Json; +using System.Linq; namespace Elements.Geometry { @@ -11,7 +12,7 @@ namespace Elements.Geometry // http://webhome.cs.uvic.ca/~blob/courses/305/notes/pdf/ref-frames.pdf /// - /// The frame type to be used for operations requiring + /// The frame type to be used for operations requiring /// a moving frame around the curve. /// public enum FrameType @@ -33,7 +34,7 @@ public enum FrameType /// /// [!code-csharp[Main](../../Elements/test/BezierTests.cs?name=example)] /// - public class Bezier : BoundedCurve + public class Bezier : BoundedCurve, IHasArcLength { private readonly int _lengthSamples = 500; @@ -105,23 +106,156 @@ public override double Length() } /// - /// Calculate the length of the bezier between start and end parameters. + /// Constants for Gauss quadrature points and weights (n = 24) + /// https://pomax.github.io/bezierinfo/legendre-gauss.html /// - /// The length of the bezier between start and end. + private static readonly double[] T = new double[] + { + -0.0640568928626056260850430826247450385909, + 0.0640568928626056260850430826247450385909, + -0.1911188674736163091586398207570696318404, + 0.1911188674736163091586398207570696318404, + -0.3150426796961633743867932913198102407864, + 0.3150426796961633743867932913198102407864, + -0.4337935076260451384870842319133497124524, + 0.4337935076260451384870842319133497124524, + -0.5454214713888395356583756172183723700107, + 0.5454214713888395356583756172183723700107, + -0.6480936519369755692524957869107476266696, + 0.6480936519369755692524957869107476266696, + -0.7401241915785543642438281030999784255232, + 0.7401241915785543642438281030999784255232, + -0.8200019859739029219539498726697452080761, + 0.8200019859739029219539498726697452080761, + -0.8864155270044010342131543419821967550873, + 0.8864155270044010342131543419821967550873, + -0.9382745520027327585236490017087214496548, + 0.9382745520027327585236490017087214496548, + -0.9747285559713094981983919930081690617411, + 0.9747285559713094981983919930081690617411, + -0.9951872199970213601799974097007368118745, + 0.9951872199970213601799974097007368118745 + }; + + /// + /// Constants for Gauss quadrature weights corresponding to the points (n = 24) + /// https://pomax.github.io/bezierinfo/legendre-gauss.html + /// + private static readonly double[] C = new double[] + { + 0.1279381953467521569740561652246953718517, + 0.1279381953467521569740561652246953718517, + 0.1258374563468282961213753825111836887264, + 0.1258374563468282961213753825111836887264, + 0.121670472927803391204463153476262425607, + 0.121670472927803391204463153476262425607, + 0.1155056680537256013533444839067835598622, + 0.1155056680537256013533444839067835598622, + 0.1074442701159656347825773424466062227946, + 0.1074442701159656347825773424466062227946, + 0.0976186521041138882698806644642471544279, + 0.0976186521041138882698806644642471544279, + 0.086190161531953275917185202983742667185, + 0.086190161531953275917185202983742667185, + 0.0733464814110803057340336152531165181193, + 0.0733464814110803057340336152531165181193, + 0.0592985849154367807463677585001085845412, + 0.0592985849154367807463677585001085845412, + 0.0442774388174198061686027482113382288593, + 0.0442774388174198061686027482113382288593, + 0.0285313886289336631813078159518782864491, + 0.0285313886289336631813078159518782864491, + 0.0123412297999871995468056670700372915759, + 0.0123412297999871995468056670700372915759 + }; + + /// + /// Computes the arc length of the Bézier curve between the given parameter values start and end. + /// https://pomax.github.io/bezierinfo/#arclength + /// + /// The starting parameter value of the Bézier curve. + /// The ending parameter value of the Bézier curve. + /// The arc length between the specified parameter values. public override double ArcLength(double start, double end) { - if (!Domain.Includes(start, true)) + double z = 0.5; // Scaling factor for the Legendre-Gauss quadrature + int len = T.Length; // Number of points in the Legendre-Gauss quadrature + + double sum = 0; // Accumulated sum for the arc length calculation + + // Iterating through the Legendre-Gauss quadrature points and weights + for (int i = 0; i < len; i++) { - throw new ArgumentOutOfRangeException("start", $"The start parameter {start} must be between {Domain.Min} and {Domain.Max}."); + double t = z * T[i] + z; // Mapping the quadrature point to the Bézier parameter range [0, 1] + Vector3 derivative = Derivative(t); // Calculating the derivative of the Bézier curve at parameter t + sum += C[i] * ArcFn(t, derivative); // Adding the weighted arc length contribution to the sum } - if (!Domain.Includes(end, true)) + + // Scaling the sum by the scaling factor and the parameter interval (end - start) to get the arc length between start and end. + return z * sum * (end - start); + } + + /// + /// Calculates the arc length contribution at parameter t based on the derivative of the Bézier curve. + /// + /// The parameter value of the Bézier curve. + /// The derivative of the Bézier curve at parameter t as a Vector3. + /// The arc length contribution at parameter t. + private double ArcFn(double t, Vector3 d) + { + // Compute the Euclidean distance of the derivative vector (d) at parameter t + return Math.Sqrt(d.X * d.X + d.Y * d.Y); + } + + /// + /// Computes the derivative of the Bézier curve at parameter t. + /// https://pomax.github.io/bezierinfo/#derivatives + /// + /// The parameter value of the Bézier curve. + /// The derivative of the Bézier curve as a Vector3. + private Vector3 Derivative(double t) + { + int n = ControlPoints.Count - 1; // Degree of the Bézier curve + Vector3[] derivatives = new Vector3[n]; // Array to store the derivative control points + + // Calculating the derivative control points using the given formula + for (int i = 0; i < n; i++) { - throw new ArgumentOutOfRangeException("end", $"The end parameter {end} must be between {Domain.Min} and {Domain.Max}."); + derivatives[i] = n * (ControlPoints[i + 1] - ControlPoints[i]); } - // TODO: We use max value here so that the calculation will continue - // until at least the end of the curve. This is not a nice solution. - return ArcLengthUntil(start, double.MaxValue, out end); + // Using the derivative control points to construct an (n-1)th degree Bézier curve at parameter t. + return BezierCurveValue(t, derivatives); + } + + /// + /// Evaluates the value of an (n-1)th degree Bézier curve at parameter t using the given control points. + /// + /// The parameter value of the Bézier curve. + /// The control points for the Bézier curve. + /// The value of the Bézier curve at parameter t as a Vector3. + private Vector3 BezierCurveValue(double t, Vector3[] controlPoints) + { + int n = controlPoints.Length - 1; // Degree of the Bézier curve + Vector3[] points = new Vector3[n + 1]; + + // Initialize the points array with the provided control points + for (int i = 0; i <= n; i++) + { + points[i] = controlPoints[i]; + } + + // De Casteljau's algorithm to evaluate the value of the Bézier curve at parameter t + for (int r = 1; r <= n; r++) + { + for (int i = 0; i <= n - r; i++) + { + points[i] = (1 - t) * points[i] + t * points[i + 1]; + } + } + + // The first element of the points array contains the value of the Bézier curve at parameter t. + return points[0]; } private double ArcLengthUntil(double start, double distance, out double end) @@ -151,6 +285,15 @@ private double ArcLengthUntil(double start, double distance, out double end) return length; } + /// + /// The mid point of the curve. + /// + /// The length based midpoint. + public virtual Vector3 MidPoint() + { + return PointAtNormalizedLength(0.5); + } + /// /// Get the point on the curve at parameter u. /// @@ -170,6 +313,157 @@ public override Vector3 PointAt(double u) return p; } + /// + /// Returns the point on the bezier corresponding to the specified length value. + /// + /// The length value along the bezier. + /// The point on the bezier corresponding to the specified length value. + /// Thrown when the specified length is out of range. + public virtual Vector3 PointAtLength(double length) + { + double totalLength = ArcLength(this.Domain.Min, this.Domain.Max); // Calculate the total length of the Bezier + if (length < 0 || length > totalLength) + { + throw new ArgumentException("The specified length is out of range."); + } + return PointAt(ParameterAtDistanceFromParameter(length, Domain.Min)); + } + + /// + /// Returns the point on the bezier corresponding to the specified normalized length-based parameter value. + /// + /// The normalized length-based parameter value, ranging from 0 to 1. + /// The point on the bezier corresponding to the specified normalized length-based parameter value. + /// Thrown when the specified parameter is out of range. + public virtual Vector3 PointAtNormalizedLength(double parameter) + { + if (parameter < 0 || parameter > 1) + { + throw new ArgumentException("The specified parameter is out of range."); + } + return PointAtLength(parameter * this.ArcLength(this.Domain.Min, this.Domain.Max)); + } + + /// + /// Finds the parameter value on the Bezier curve that corresponds to the given 3D point within a specified threshold. + /// + /// The 3D point to find the corresponding parameter for. + /// The maximum distance threshold to consider a match between the projected point and the original point. + /// The parameter value on the Bezier curve if the distance between the projected point and the original point is within the threshold, otherwise returns null. + public double? ParameterAt(Vector3 point, double threshold = 0.0001) + { + // Find the parameter corresponding to the projected point on the Bezier curve + var parameter = ProjectedPoint(point, threshold); + + if (parameter == null) + { + // If the projected point does not return a relevant parameter return null + return null; + } + // Find the 3D point on the Bezier curve at the obtained parameter value + var projection = PointAt((double)parameter); + + // Check if the distance between the projected point and the original point is within + // a tolerence of the threshold + if (projection.DistanceTo(point) < (threshold * 10) - threshold) + { + // If the distance is within the threshold, return the parameter value + return parameter; + } + else + { + // If the distance exceeds the threshold, consider the point as not on the Bezier curve and return null + return null; + } + } + + /// + /// Projects a 3D point onto the Bezier curve to find the parameter value of the closest point on the curve. + /// + /// The 3D point to project onto the Bezier curve. + /// The maximum threshold to refine the projection and find the closest point. + /// The parameter value on the Bezier curve corresponding to the projected point, or null if the projection is not within the specified threshold. + public double? ProjectedPoint(Vector3 point, double threshold = 0.001) + { + // https://pomax.github.io/bezierinfo/#projections + // Generate a lookup table (LUT) of points and their corresponding parameter values on the Bezier curve + List<(Vector3 point, double t)> lut = GenerateLookupTable(); + + // Initialize variables to store the closest distance (d) and the index (index) of the closest point in the lookup table + double d = double.MaxValue; + int index = 0; + + // Find the closest point to the input point in the lookup table (LUT) using Euclidean distance + for (int i = 0; i < lut.Count; i++) + { + double q = Math.Sqrt(Math.Pow((point - lut[i].point).X, 2) + Math.Pow((point - lut[i].point).Y, 2) + Math.Pow((point - lut[i].point).Z, 2)); + if (q < d) + { + d = q; + index = i; + } + } + + // Obtain the parameter values of the neighboring points in the LUT for further refinement + double t1 = lut[Math.Max(index - 1, 0)].t; + double t2 = lut[Math.Min(index + 1, lut.Count - 1)].t; + double v = t2 - t1; + + // Refine the projection by iteratively narrowing down the parameter range to find the closest point + while (v > threshold) + { + // Calculate intermediate parameter values + double t0 = t1 + v / 4; + double t3 = t2 - v / 4; + + // Calculate corresponding points on the Bezier curve using the intermediate parameter values + Vector3 p0 = PointAt(t0); + Vector3 p3 = PointAt(t3); + + // Calculate the distances between the input point and the points on the Bezier curve + double d0 = Math.Sqrt(Math.Pow((point - p0).X, 2) + Math.Pow((point - p0).Y, 2) + Math.Pow((point - p0).Z, 2)); + double d3 = Math.Sqrt(Math.Pow((point - p3).X, 2) + Math.Pow((point - p3).Y, 2) + Math.Pow((point - p3).Z, 2)); + + // Choose the sub-range that is closer to the input point and update the range + if (d0 < d3) + { + t2 = t3; + } + else + { + t1 = t0; + } + + // Update the range difference for the next iteration + v = t2 - t1; + } + + // Return the average of the refined parameter values as the projection of the input point on the Bezier curve + return (t1 + t2) / 2; + } + + /// + /// Generates a lookup table of points and their corresponding parameter values on the Bezier curve. + /// + /// Number of samples to take along the curve. + /// A list of tuples containing the sampled points and their corresponding parameter values on the Bezier curve. + private List<(Vector3 point, double t)> GenerateLookupTable(int numSamples = 100) + { + // Initialize an empty list to store the lookup table (LUT) + List<(Vector3 point, double t)> lut = new List<(Vector3 point, double t)>(); + + // Generate lookup table by sampling points on the Bezier curve + for (int i = 0; i <= numSamples; i++) + { + double t = (double)i / numSamples; // Calculate the parameter value based on the current sample index + Vector3 point = PointAt(t); // Get the 3D point on the Bezier curve corresponding to the current parameter value + lut.Add((point, t)); // Add the sampled point and its corresponding parameter value to the lookup table (LUT) + } + + // Return the completed lookup table (LUT) + return lut; + } + private double BinomialCoefficient(int n, int i) { return Factorial(n) / (Factorial(i) * Factorial(n - i)); @@ -343,314 +637,5 @@ public override double ParameterAtDistanceFromParameter(double distance, double ArcLengthUntil(start, distance, out var end); return end; } - - /// - public override bool Intersects(ICurve curve, out List results) - { - switch (curve) - { - case Line line: - return line.Intersects(this, out results); - case Arc arc: - return arc.Intersects(this, out results); - case EllipticalArc ellipticArc: - return ellipticArc.Intersects(this, out results); - case InfiniteLine line: - return Intersects(line, out results); - case Circle circle: - return Intersects(circle, out results); - case Ellipse ellipse: - return Intersects(ellipse, out results); - case IndexedPolycurve polycurve: - return polycurve.Intersects(this, out results); - case Bezier bezier: - return Intersects(bezier, out results); - default: - throw new NotImplementedException(); - } - } - - /// - /// Does this bezier curve intersects with an infinite line? - /// Iterative approximation is used to find intersections. - /// - /// Infinite line to intersect. - /// List containing intersection points. - /// True if intersection exists, otherwise false. - public bool Intersects(InfiniteLine line, out List results) - { - BBox3 box = new BBox3(ControlPoints); - // Bezier curve always inside it's bounding box. - // Rough check if line intersects curve box. - if (!new Line(line.Origin, line.Origin + line.Direction).Intersects( - box, out _, true)) - { - results = new List(); - return false; - } - - var l = this.Length(); - var div = (int)Math.Round(l / DefaultMinimumChordLength); - div = Math.Min(div, _lengthSamples); - - // Iteratively, find points on Bezier with 0 distance to the line. - // It Bezier was limited to 4 points - more effective approach could be used. - var roots = Equations.SolveIterative(Domain.Min, Domain.Max, div, - new Func((t) => - { - var p = PointAt(t); - return (p - p.ClosestPointOn(line)).LengthSquared(); - }), Vector3.EPSILON * Vector3.EPSILON); - - results = roots.Select(r => PointAt(r)).UniqueAverageWithinTolerance( - Vector3.EPSILON * 2).ToList(); - return results.Any(); - } - - /// - /// Does this Bezier curve intersects with a circle? - /// Iterative approximation is used to find intersections. - /// - /// Circle to intersect. - /// List containing intersection points. - /// True if any intersections exist, otherwise false. - public bool Intersects(Circle circle, out List results) - { - BBox3 box = new BBox3(ControlPoints); - // Bezier curve always inside it's bounding box. - // Rough check if curve is too far away. - var boxCenter = box.Center(); - if (circle.Center.DistanceTo(boxCenter) > - circle.Radius + (box.Max - boxCenter).Length()) - { - results = new List(); - return false; - } - - var l = this.Length(); - var div = (int)Math.Round(l / DefaultMinimumChordLength); - div = Math.Min(div, _lengthSamples); - - // Iteratively, find points on Bezier with radius distance to the circle. - var invertedT = circle.Transform.Inverted(); - var roots = Equations.SolveIterative(Domain.Min, Domain.Max, div, - new Func((t) => - { - var p = PointAt(t); - var local = invertedT.OfPoint(p); - return local.LengthSquared() - circle.Radius * circle.Radius; - }), Vector3.EPSILON * Vector3.EPSILON); - - results = roots.Select(r => PointAt(r)).UniqueAverageWithinTolerance( - Vector3.EPSILON * 2).ToList(); - return results.Any(); - } - - /// - /// Does this Bezier curve intersects with an ellipse? - /// Iterative approximation is used to find intersections. - /// - /// Ellipse to intersect. - /// List containing intersection points. - /// True if any intersections exist, otherwise false. - public bool Intersects(Ellipse ellipse, out List results) - { - BBox3 box = new BBox3(ControlPoints); - // Bezier curve always inside it's bounding box. - // Rough check if curve is too far away. - var boxCenter = box.Center(); - if (ellipse.Center.DistanceTo(boxCenter) > - Math.Max(ellipse.MajorAxis, ellipse.MinorAxis) + (box.Max - boxCenter).Length()) - { - results = new List(); - return false; - } - - var l = this.Length(); - var div = (int)Math.Round(l / DefaultMinimumChordLength); - div = Math.Min(div, _lengthSamples); - - // Iteratively, find points on ellipse with distance - // to other ellipse equal to its focal distance. - var invertedT = ellipse.Transform.Inverted(); - var roots = Equations.SolveIterative(Domain.Min, Domain.Max, div, - new Func((t) => - { - var p = PointAt(t); - var local = invertedT.OfPoint(p); - var dx = Math.Pow(local.X / ellipse.MajorAxis, 2); - var dy = Math.Pow(local.Y / ellipse.MinorAxis, 2); - return dx + dy + local.Z * local.Z - 1; - }), Vector3.EPSILON * Vector3.EPSILON); - - results = roots.Select(r => PointAt(r)).UniqueAverageWithinTolerance( - Vector3.EPSILON * 2).ToList(); - return results.Any(); - } - - /// - /// Does this Bezier curve intersects with other Bezier curve? - /// Iterative approximation is used to find intersections. - /// - /// Other Bezier curve to intersect. - /// List containing intersection points. - /// True if any intersections exist, otherwise false. - public bool Intersects(Bezier other, out List results) - { - results = new List(); - int leftSteps = BoxApproximationStepsCount(); - int rightSteps = other.BoxApproximationStepsCount(); - - var leftCache = new Dictionary(); - var rightCache = new Dictionary(); - - BBox3 box = CurveBox(leftSteps, leftCache); - BBox3 otherBox = other.CurveBox(rightSteps, rightCache); - - Intersects(other, - (box, Domain), - (otherBox, other.Domain), - leftCache, - rightCache, - ref results); - - // Subdivision algorithm produces duplicates, all tolerance away from real answer. - // Grouping and averaging them improves output as we as eliminates duplications. - results = results.UniqueAverageWithinTolerance().ToList(); - return results.Any(); - } - - private BBox3 CurveBox(int numSteps, Dictionary cache) - { - List points = new List(); - double step = Domain.Length / numSteps; - for (int i = 0; i < numSteps; i++) - { - var t = Domain.Min + i * step; - points.Add(PointAtCached(t, cache)); - } - points.Add(PointAtCached(Domain.Max, cache)); - BBox3 box = new BBox3(points); - return box; - } - - private void Intersects(Bezier other, - (BBox3 Box, Domain1d Domain) left, - (BBox3 Box, Domain1d Domain) right, - Dictionary leftCache, - Dictionary rightCache, - ref List results) - { - // If bounding boxes of two curves (not control points) are not intersect - // curves not intersect. - if (!left.Box.Intersects(right.Box)) - { - return; - } - - var leftSplit = SplitCurveBox(left, leftCache, out var loLeft, out var hiLeft); - var rightSplit = other.SplitCurveBox(right, rightCache, out var loRight, out var hiRight); - - // If boxes of two curves are tolerance sized - - // average point of their centers is treated as intersection point. - if (!leftSplit && !rightSplit) - { - results.Add(left.Box.Center().Average(right.Box.Center())); - return; - } - - // Recursively repeat process on box of subdivided curves until they are small enough. - // Each pair, which bounding boxes are not intersecting are discarded. - if (!leftSplit) - { - Intersects(other, left, loRight, leftCache, rightCache, ref results); - Intersects(other, left, hiRight, leftCache, rightCache, ref results); - } - else if (!rightSplit) - { - Intersects(other, loLeft, right, leftCache, rightCache, ref results); - Intersects(other, hiLeft, right, leftCache, rightCache, ref results); - } - else - { - Intersects(other, loLeft, loRight, leftCache, rightCache, ref results); - Intersects(other, hiLeft, loRight, leftCache, rightCache, ref results); - Intersects(other, loLeft, hiRight, leftCache, rightCache, ref results); - Intersects(other, hiLeft, hiRight, leftCache, rightCache, ref results); - } - } - - private bool SplitCurveBox((BBox3 Box, Domain1d Domain) def, - Dictionary cache, - out (BBox3 Box, Domain1d Domain) low, - out (BBox3 Box, Domain1d Domain) high) - { - low = (default, default); - high = (default, default); - - // If curve bounding box is tolerance size - it's considered as intersection. - // Otherwise calculate new boxes of two halves of the curve. - var epsilon2 = Vector3.EPSILON * Vector3.EPSILON; - var leftConvergent = (def.Box.Max - def.Box.Min).LengthSquared() < epsilon2 * 2; - if (leftConvergent) - { - return false; - } - - // If curve bounding box is tolerance size - it's considered as intersection. - // Otherwise calculate new boxes of two halves of the curve. - low = CurveBoxHalf(def, cache, true); - high = CurveBoxHalf(def, cache, false); - return true; - } - - /// - /// Get bounding box of curve segment (not control points) for half of given domain. - /// - /// Curve segment definition - box, domain and number of subdivisions. - /// Dictionary of precomputed points at parameter. - /// Take lower of higher part of curve. - /// Definition of curve segment that is half of given. - private (BBox3 Box, Domain1d Domain) CurveBoxHalf( - (BBox3 Box, Domain1d Domain) def, - Dictionary cache, - bool low) - { - var steps = BoxApproximationStepsCount(); - double step = (def.Domain.Length / 2) / steps; - var min = low ? def.Domain.Min : def.Domain.Min + def.Domain.Length / 2; - var max = low ? def.Domain.Min + def.Domain.Length / 2 : def.Domain.Max; - - List points = new List(); - for (int i = 0; i < steps; i++) - { - var t = min + i * step; - points.Add(PointAtCached(t, cache)); - } - points.Add(PointAtCached(max, cache)); - - var box = new BBox3(points); - var domain = new Domain1d(min, max); - return (box, domain); - } - - private Vector3 PointAtCached(double t, Dictionary cache) - { - if (cache.TryGetValue(t, out var p)) - { - return p; - } - else - { - p = PointAt(t); - cache.Add(t, p); - return p; - } - } - - private int BoxApproximationStepsCount() - { - return ControlPoints.Count * 8 - 1; - } } } \ No newline at end of file diff --git a/Elements/src/Geometry/Polyline.cs b/Elements/src/Geometry/Polyline.cs index cdfedc343..543a240b8 100644 --- a/Elements/src/Geometry/Polyline.cs +++ b/Elements/src/Geometry/Polyline.cs @@ -488,8 +488,6 @@ private Transform CreateOrthogonalTransform(int i, Vector3 origin, Vector3 up) /// A list of points representing the segments. public Vector3[] DivideByLength(double divisionLength) { - var segments = new List(); - if (this.Vertices.Count < 2) { // Handle invalid polyline with insufficient vertices @@ -497,7 +495,7 @@ public Vector3[] DivideByLength(double divisionLength) } var currentProgression = 0.0; - segments = new List { this.Vertices.FirstOrDefault() }; + var segments = new List { this.Vertices.FirstOrDefault() }; foreach (var currentSegment in this.Segments()) { diff --git a/Elements/test/BezierTests.cs b/Elements/test/BezierTests.cs index 2df434b0e..d639ed100 100644 --- a/Elements/test/BezierTests.cs +++ b/Elements/test/BezierTests.cs @@ -147,8 +147,8 @@ public void IntersectsPolycurve() var polygon = new Polygon(new Vector3[]{ (0, 3), (6, 3), (4, 1), (-2, 1) - }); - + }); + Assert.True(bezier.Intersects(polygon, out var results)); Assert.Equal(4, results.Count); Assert.Contains(new Vector3(0.93475, 3), results); @@ -173,5 +173,166 @@ public void IntersectsBezier() Assert.Contains(new Vector3(0.5755, 2.5), results); Assert.Contains(new Vector3(4.4245, 2.5), results); } + + public void Bezier_ArcLength() + { + var a = new Vector3(50, 150, 0); + var b = new Vector3(105, 66, 0); + var c = new Vector3(170, 230, 0); + var d = new Vector3(200, 150, 0); + var ctrlPts = new List { a, b, c, d }; + var bezier = new Bezier(ctrlPts); + + var expectedLength = 184.38886379602502; // approximation as the integral function used for calculating length is not 100% accurate + Assert.Equal(expectedLength, bezier.ArcLength(0.0, 1.0), 2); + } + + [Fact] + public void GetParameterAt() + { + var tolerance = 0.00001; + + var a = new Vector3(1, 5, 0); + var b = new Vector3(5, 20, 0); + var c = new Vector3(5, -10, 0); + var d = new Vector3(9, 5, 0); + var ctrlPts = new List { a, b, c, d }; + var bezier = new Bezier(ctrlPts); + + var samplePt = new Vector3(0, 0, 0); + Assert.Null(bezier.ParameterAt(samplePt, tolerance)); + + samplePt = new Vector3(6.625, 0.7812, 0.0); + Assert.True((double)bezier.ParameterAt(samplePt, tolerance) - 0.75 <= tolerance * 10); + } + + [Fact] + public void GetPointAt() + { + var a = new Vector3(1, 5, 0); + var b = new Vector3(5, 20, 0); + var c = new Vector3(5, -10, 0); + var d = new Vector3(9, 5, 0); + var ctrlPts = new List { a, b, c, d }; + var bezier = new Bezier(ctrlPts); + + var testPt = new Vector3(3.3750, 9.21875, 0.0000); + Assert.True(testPt.Equals(bezier.PointAt(0.25))); + + testPt = new Vector3(4.699, 6.11375, 0.0000); + Assert.True(testPt.Equals(bezier.PointAtNormalized(0.45))); + + testPt = new Vector3(4.515904, 6.75392, 0.0000); + Assert.True(testPt.Equals(bezier.PointAtLength(8.0))); + + testPt = new Vector3(3.048823, 9.329262, 0.0000); + Assert.True(testPt.Equals(bezier.PointAtNormalizedLength(0.25))); + } + + [Fact] + public void DivideByLength() + { + var a = new Vector3(50, 150, 0); + var b = new Vector3(105, 66, 0); + var c = new Vector3(170, 230, 0); + var d = new Vector3(200, 150, 0); + var ctrlPts = new List { a, b, c, d }; + var bezier = new Bezier(ctrlPts); + + var testPts = new List(){ + new Vector3(50.00, 150.00, 0.00), + new Vector3(90.705919, 125.572992, 0.00), + new Vector3(134.607122, 148.677130, 0.00), + new Vector3(177.600231, 172.675064, 0.00), + new Vector3(200.00, 150.00, 0.00), + }; + + var ptsFromBezier = bezier.DivideByLength(50.0); + + for (int i = 0; i < ptsFromBezier.Length; i++) + { + Assert.True(ptsFromBezier[i].Equals(testPts[i])); + } + } + + [Fact] + public void Split() + { + var a = new Vector3(50, 150, 0); + var b = new Vector3(105, 66, 0); + var c = new Vector3(170, 230, 0); + var d = new Vector3(200, 150, 0); + var ctrlPts = new List { a, b, c, d }; + var bezier = new Bezier(ctrlPts); + + var testBeziers = new List(){ + new Bezier(new List() { + new Vector3(50.00, 150.00, 0.00), + new Vector3(63.75, 129.00, 0.00), + new Vector3(78.1250, 123.50, 0.00), + new Vector3(92.421875, 125.8125, 0.00) + } + ), + new Bezier(new List() { + new Vector3(92.421875, 125.8125, 0.00), + new Vector3(121.015625, 130.4375, 0.00), + new Vector3(149.296875, 166.3125, 0.00), + new Vector3(171.640625, 171.9375, 0.00) + } + ), + new Bezier(new List() { + new Vector3(171.640625, 171.9375, 0.00), + new Vector3(182.8125, 174.7500, 0.00), + new Vector3(192.50, 170.00, 0.00), + new Vector3(200.00, 150.00, 0.00) + } + ), + }; + + var beziers = bezier.Split(new List() { 0.25, 0.75 }); + + for (int i = 0; i < beziers.Count; i++) + { + Assert.True(beziers[i].ControlPoints[0].Equals(testBeziers[i].ControlPoints[0])); + Assert.True(beziers[i].ControlPoints[1].Equals(testBeziers[i].ControlPoints[1])); + Assert.True(beziers[i].ControlPoints[2].Equals(testBeziers[i].ControlPoints[2])); + Assert.True(beziers[i].ControlPoints[3].Equals(testBeziers[i].ControlPoints[3])); + } + + testBeziers = new List(){ + new Bezier(new List() { + new Vector3(50.00, 150.00, 0.00), + new Vector3(61.88, 131.856, 0.00), + new Vector3(74.22656, 125.282688, 0.00), + new Vector3(86.586183, 125.321837, 0.00) + } + ), + new Bezier(new List() { + new Vector3(86.586183, 125.321837, 0.00), + new Vector3(114.738659, 125.411011, 0.00), + new Vector3(142.958913, 159.807432, 0.00), + new Vector3(165.887648, 169.916119, 0.00) + } + ), + new Bezier(new List() { + new Vector3(165.887648, 169.916119, 0.00), + new Vector3(179.49576, 175.915584, 0.00), + new Vector3(191.24, 173.36, 0.00), + new Vector3(200.00, 150.00, 0.00) + } + ), + }; + + var normalizedBeziers = bezier.Split(new List() { 0.25, 0.75 }, true); + + for (int i = 0; i < normalizedBeziers.Count; i++) + { + Assert.True(normalizedBeziers[i].ControlPoints[0].Equals(testBeziers[i].ControlPoints[0])); + Assert.True(normalizedBeziers[i].ControlPoints[1].Equals(testBeziers[i].ControlPoints[1])); + Assert.True(normalizedBeziers[i].ControlPoints[2].Equals(testBeziers[i].ControlPoints[2])); + Assert.True(normalizedBeziers[i].ControlPoints[3].Equals(testBeziers[i].ControlPoints[3])); + } +>>>>>>> 374fe346(Polish Bezier, IHasCurveLength) + } } } \ No newline at end of file From 9d6908c239b175ab603c8fc8da518509094bdcd6 Mon Sep 17 00:00:00 2001 From: James Bradley Date: Fri, 21 Jul 2023 18:05:33 -0400 Subject: [PATCH 10/22] forgot the Z --- Elements/src/Geometry/Bezier.cs | 272 ++++++++++++++++++++++++++++++++ 1 file changed, 272 insertions(+) diff --git a/Elements/src/Geometry/Bezier.cs b/Elements/src/Geometry/Bezier.cs index af7248b72..2ba19994e 100644 --- a/Elements/src/Geometry/Bezier.cs +++ b/Elements/src/Geometry/Bezier.cs @@ -637,5 +637,277 @@ public override double ParameterAtDistanceFromParameter(double distance, double ArcLengthUntil(start, distance, out var end); return end; } + + /// + /// Divides the bezier into segments of the specified length. + /// + /// The desired length of each segment. + /// A list of points representing the segment divisions. + public Vector3[] DivideByLength(double divisionLength) + { + var totalLength = this.ArcLength(Domain.Min, Domain.Max); + if (totalLength <= 0) + { + // Handle invalid bezier with insufficient length + return new Vector3[0]; + } + var parameter = ParameterAtDistanceFromParameter(divisionLength, Domain.Min); + var segments = new List { this.Start }; + + while (parameter < Domain.Max) + { + segments.Add(PointAt(parameter)); + var newParameter = ParameterAtDistanceFromParameter(divisionLength, parameter); + parameter = newParameter != parameter ? newParameter : Domain.Max; + } + + // Add the last vertex of the bezier as the endpoint of the last segment if it + // is not already part of the list + if (!segments[segments.Count - 1].IsAlmostEqualTo(this.End)) + { + segments.Add(this.End); + } + + return segments.ToArray(); + } + + /// + /// Divides the bezier into segments of the specified length. + /// + /// The desired length of each segment. + /// A list of beziers representing the segments. + public List SplitByLength(double divisionLength) + { + var totalLength = this.ArcLength(Domain.Min, Domain.Max); + if (totalLength <= 0) + { + // Handle invalid bezier with insufficient length + return null; + } + var currentParameter = ParameterAtDistanceFromParameter(divisionLength, Domain.Min); + var parameters = new List { this.Domain.Min }; + + while (currentParameter < Domain.Max) + { + parameters.Add(currentParameter); + var newParameter = ParameterAtDistanceFromParameter(divisionLength, currentParameter); + currentParameter = newParameter != currentParameter ? newParameter : Domain.Max; + } + + // Add the last vertex of the bezier as the endpoint of the last segment if it + // is not already part of the list + if (!parameters[parameters.Count - 1].ApproximatelyEquals(this.Domain.Max)) + { + parameters.Add(this.Domain.Max); + } + + return Split(parameters); + } + + /// + /// Splits the Bezier curve into segments at specified parameter values. + /// + /// The list of parameter values to split the curve at. + /// If true the parameters will be length normalized. + /// A list of Bezier segments obtained after splitting. + public List Split(List parameters, bool normalize = false) + { + // Calculate the total length of the Bezier curve + var totalLength = this.ArcLength(Domain.Min, Domain.Max); + + // Check for invalid curve with insufficient length + if (totalLength <= 0) + { + throw new InvalidOperationException($"Invalid bezier with insufficient length. Total Length = {totalLength}"); + } + + // Check if the list of parameters is empty or null + if (parameters == null || parameters.Count == 0) + { + throw new ArgumentException("No split points provided."); + } + + // Initialize a list to store the resulting Bezier segments + var segments = new List(); + var bezier = this; // Create a reference to the original Bezier curve + + if (normalize) + { + parameters = parameters.Select(parameter => ParameterAtDistanceFromParameter(parameter * this.ArcLength(this.Domain.Min, this.Domain.Max), Domain.Min)).ToList(); + } + parameters.Sort(); // Sort the parameters in ascending order + + // Iterate through each parameter to split the curve + for (int i = 0; i < parameters.Count; i++) + { + double t = (Domain.Min <= parameters[i] && parameters[i] <= Domain.Max) + ? parameters[i] // Ensure the parameter is within the domain + : throw new ArgumentException($"Parameter {parameters[i]} is not within the domain ({Domain.Min}->{Domain.Max}) of the Bezier curve."); + + // Check if the parameter is within the valid range [0, 1] + if (t >= 0 && t <= 1) + { + // Split the curve at the given parameter and obtain the two resulting Bezier segments + var tuple = bezier.SplitAt(t); + + // Store the first split Bezier in the list + segments.Add(tuple.Item1); + + // Update bezier to the second split Bezier to continue splitting + bezier = tuple.Item2; + + // Remap subsequent parameters to the new Bezier curve's parameter space + for (int j = i + 1; j < parameters.Count; j++) + { + parameters[j] = (parameters[j] - t) / (1 - t); + } + } + } + + segments.Add(bezier); + + // Return the list of Bezier segments obtained after splitting + return segments; + } + + /// + /// Splits the bezier curve at the given parameter value. + /// + /// The parameter value at which to split the curve. + /// A tuple containing two split bezier curves. + public Tuple SplitAt(double t) + { + // Extract the control points from the input bezier + var startPoint = ControlPoints[0]; + var controlPoint1 = ControlPoints[1]; + var controlPoint2 = ControlPoints[2]; + var endPoint = ControlPoints[3]; + + // Compute the intermediate points using de Casteljau's algorithm + var q0 = (1 - t) * startPoint + t * controlPoint1; + var q1 = (1 - t) * controlPoint1 + t * controlPoint2; + var q2 = (1 - t) * controlPoint2 + t * endPoint; + + var r0 = (1 - t) * q0 + t * q1; + var r1 = (1 - t) * q1 + t * q2; + + // Compute the split point on the bezier curve + var splitPoint = (1 - t) * r0 + t * r1; + + // Construct the first split bezier curve + var subBezier1 = new Bezier(new List() { startPoint, q0, r0, splitPoint }); + + // Construct the second split bezier curve + var subBezier2 = new Bezier(new List() { splitPoint, r1, q2, endPoint }); + + // Return a tuple containing the split bezier curves + return new Tuple(subBezier1, subBezier2); + } + + /// + /// Constructs piecewise cubic Bézier curves from a list of points using control points calculated with the specified looseness. + /// + /// The list of points defining the path. + /// The looseness factor used to calculate control points. A higher value results in smoother curves. + /// If true, the path will be closed, connecting the last point with the first one. + /// A list of piecewise cubic Bézier curves approximating the path defined by the input points. + public static List ConstructPiecewiseCubicBezier(List points, double looseness = 6.0, bool close = false) + { + List beziers = new List(); + + // Calculate the control points. + List controlPoints = CalculateControlPoints(points, looseness); + + // Create the start Bezier curve. + Bezier startBezier = new Bezier( + new List + { + points[0], + controlPoints[0][1], + points[1] + } + ); + + // Add the start Bezier curve to the list. + beziers.Add(startBezier); + + // Iterate through pairs of points. + for (int i = 1; i < points.Count - 2; i++) + { + // Create the control points. + List bezierControlPoints = new List + { + points[i], + controlPoints[i - 1][0], + controlPoints[i][1], + points[i + 1] + }; + + // Create the Bezier curve. + Bezier bezier = new Bezier(bezierControlPoints); + + // Add the Bezier curve to the list. + beziers.Add(bezier); + } + + // Create the end Bezier curve. + Bezier endBezier = new Bezier( + new List + { + points[points.Count() - 1], + controlPoints[controlPoints.Count() - 1][0], + points[points.Count() - 2] + } + ); + + // Add the end Bezier curve to the list. + beziers.Add(endBezier); + + // Return the list of Bezier curves. + return beziers; + } + + /// + /// Calculates the control points for constructing piecewise cubic Bézier curves from the given list of points and looseness factor. + /// + /// The list of points defining the path. + /// The looseness factor used to calculate control points. A higher value results in smoother curves. + /// A list of control points (pairs of Vector3) for the piecewise cubic Bézier curves. + private static List CalculateControlPoints(List points, double looseness) + { + List controlPoints = new List(); + + for (int i = 1; i < points.Count - 1; i++) + { + // Calculate the differences in x and y coordinates. + var dx = points[i - 1].X - points[i + 1].X; + var dy = points[i - 1].Y - points[i + 1].Y; + var dz = points[i - 1].Z - points[i + 1].Z; + + // Calculate the control point coordinates. + var controlPointX1 = points[i].X - dx * (1 / looseness); + var controlPointY1 = points[i].Y - dy * (1 / looseness); + var controlPointZ1 = points[i].Z - dz * (1 / looseness); + var controlPoint1 = new Vector3(controlPointX1, controlPointY1, controlPointZ1); + + var controlPointX2 = points[i].X + dx * (1 / looseness); + var controlPointY2 = points[i].Y + dy * (1 / looseness); + var controlPointZ2 = points[i].Z + dz * (1 / looseness); + var controlPoint2 = new Vector3(controlPointX2, controlPointY2, controlPointZ2); + + // Create an array to store the control points. + Vector3[] controlPointArray = new Vector3[] + { + controlPoint1, + controlPoint2 + }; + + // Add the control points to the list. + controlPoints.Add(controlPointArray); + } + + // Return the list of control points. + return controlPoints; + } } } \ No newline at end of file From 0a255a7b44a295456b7396f928b176285d4ef1c2 Mon Sep 17 00:00:00 2001 From: James Bradley Date: Mon, 13 Nov 2023 09:25:39 -0500 Subject: [PATCH 11/22] forgot the Z --- Elements/src/Geometry/Bezier.cs | 272 ++++++++++++++++++++++++++++++++ 1 file changed, 272 insertions(+) diff --git a/Elements/src/Geometry/Bezier.cs b/Elements/src/Geometry/Bezier.cs index af7248b72..2ba19994e 100644 --- a/Elements/src/Geometry/Bezier.cs +++ b/Elements/src/Geometry/Bezier.cs @@ -637,5 +637,277 @@ public override double ParameterAtDistanceFromParameter(double distance, double ArcLengthUntil(start, distance, out var end); return end; } + + /// + /// Divides the bezier into segments of the specified length. + /// + /// The desired length of each segment. + /// A list of points representing the segment divisions. + public Vector3[] DivideByLength(double divisionLength) + { + var totalLength = this.ArcLength(Domain.Min, Domain.Max); + if (totalLength <= 0) + { + // Handle invalid bezier with insufficient length + return new Vector3[0]; + } + var parameter = ParameterAtDistanceFromParameter(divisionLength, Domain.Min); + var segments = new List { this.Start }; + + while (parameter < Domain.Max) + { + segments.Add(PointAt(parameter)); + var newParameter = ParameterAtDistanceFromParameter(divisionLength, parameter); + parameter = newParameter != parameter ? newParameter : Domain.Max; + } + + // Add the last vertex of the bezier as the endpoint of the last segment if it + // is not already part of the list + if (!segments[segments.Count - 1].IsAlmostEqualTo(this.End)) + { + segments.Add(this.End); + } + + return segments.ToArray(); + } + + /// + /// Divides the bezier into segments of the specified length. + /// + /// The desired length of each segment. + /// A list of beziers representing the segments. + public List SplitByLength(double divisionLength) + { + var totalLength = this.ArcLength(Domain.Min, Domain.Max); + if (totalLength <= 0) + { + // Handle invalid bezier with insufficient length + return null; + } + var currentParameter = ParameterAtDistanceFromParameter(divisionLength, Domain.Min); + var parameters = new List { this.Domain.Min }; + + while (currentParameter < Domain.Max) + { + parameters.Add(currentParameter); + var newParameter = ParameterAtDistanceFromParameter(divisionLength, currentParameter); + currentParameter = newParameter != currentParameter ? newParameter : Domain.Max; + } + + // Add the last vertex of the bezier as the endpoint of the last segment if it + // is not already part of the list + if (!parameters[parameters.Count - 1].ApproximatelyEquals(this.Domain.Max)) + { + parameters.Add(this.Domain.Max); + } + + return Split(parameters); + } + + /// + /// Splits the Bezier curve into segments at specified parameter values. + /// + /// The list of parameter values to split the curve at. + /// If true the parameters will be length normalized. + /// A list of Bezier segments obtained after splitting. + public List Split(List parameters, bool normalize = false) + { + // Calculate the total length of the Bezier curve + var totalLength = this.ArcLength(Domain.Min, Domain.Max); + + // Check for invalid curve with insufficient length + if (totalLength <= 0) + { + throw new InvalidOperationException($"Invalid bezier with insufficient length. Total Length = {totalLength}"); + } + + // Check if the list of parameters is empty or null + if (parameters == null || parameters.Count == 0) + { + throw new ArgumentException("No split points provided."); + } + + // Initialize a list to store the resulting Bezier segments + var segments = new List(); + var bezier = this; // Create a reference to the original Bezier curve + + if (normalize) + { + parameters = parameters.Select(parameter => ParameterAtDistanceFromParameter(parameter * this.ArcLength(this.Domain.Min, this.Domain.Max), Domain.Min)).ToList(); + } + parameters.Sort(); // Sort the parameters in ascending order + + // Iterate through each parameter to split the curve + for (int i = 0; i < parameters.Count; i++) + { + double t = (Domain.Min <= parameters[i] && parameters[i] <= Domain.Max) + ? parameters[i] // Ensure the parameter is within the domain + : throw new ArgumentException($"Parameter {parameters[i]} is not within the domain ({Domain.Min}->{Domain.Max}) of the Bezier curve."); + + // Check if the parameter is within the valid range [0, 1] + if (t >= 0 && t <= 1) + { + // Split the curve at the given parameter and obtain the two resulting Bezier segments + var tuple = bezier.SplitAt(t); + + // Store the first split Bezier in the list + segments.Add(tuple.Item1); + + // Update bezier to the second split Bezier to continue splitting + bezier = tuple.Item2; + + // Remap subsequent parameters to the new Bezier curve's parameter space + for (int j = i + 1; j < parameters.Count; j++) + { + parameters[j] = (parameters[j] - t) / (1 - t); + } + } + } + + segments.Add(bezier); + + // Return the list of Bezier segments obtained after splitting + return segments; + } + + /// + /// Splits the bezier curve at the given parameter value. + /// + /// The parameter value at which to split the curve. + /// A tuple containing two split bezier curves. + public Tuple SplitAt(double t) + { + // Extract the control points from the input bezier + var startPoint = ControlPoints[0]; + var controlPoint1 = ControlPoints[1]; + var controlPoint2 = ControlPoints[2]; + var endPoint = ControlPoints[3]; + + // Compute the intermediate points using de Casteljau's algorithm + var q0 = (1 - t) * startPoint + t * controlPoint1; + var q1 = (1 - t) * controlPoint1 + t * controlPoint2; + var q2 = (1 - t) * controlPoint2 + t * endPoint; + + var r0 = (1 - t) * q0 + t * q1; + var r1 = (1 - t) * q1 + t * q2; + + // Compute the split point on the bezier curve + var splitPoint = (1 - t) * r0 + t * r1; + + // Construct the first split bezier curve + var subBezier1 = new Bezier(new List() { startPoint, q0, r0, splitPoint }); + + // Construct the second split bezier curve + var subBezier2 = new Bezier(new List() { splitPoint, r1, q2, endPoint }); + + // Return a tuple containing the split bezier curves + return new Tuple(subBezier1, subBezier2); + } + + /// + /// Constructs piecewise cubic Bézier curves from a list of points using control points calculated with the specified looseness. + /// + /// The list of points defining the path. + /// The looseness factor used to calculate control points. A higher value results in smoother curves. + /// If true, the path will be closed, connecting the last point with the first one. + /// A list of piecewise cubic Bézier curves approximating the path defined by the input points. + public static List ConstructPiecewiseCubicBezier(List points, double looseness = 6.0, bool close = false) + { + List beziers = new List(); + + // Calculate the control points. + List controlPoints = CalculateControlPoints(points, looseness); + + // Create the start Bezier curve. + Bezier startBezier = new Bezier( + new List + { + points[0], + controlPoints[0][1], + points[1] + } + ); + + // Add the start Bezier curve to the list. + beziers.Add(startBezier); + + // Iterate through pairs of points. + for (int i = 1; i < points.Count - 2; i++) + { + // Create the control points. + List bezierControlPoints = new List + { + points[i], + controlPoints[i - 1][0], + controlPoints[i][1], + points[i + 1] + }; + + // Create the Bezier curve. + Bezier bezier = new Bezier(bezierControlPoints); + + // Add the Bezier curve to the list. + beziers.Add(bezier); + } + + // Create the end Bezier curve. + Bezier endBezier = new Bezier( + new List + { + points[points.Count() - 1], + controlPoints[controlPoints.Count() - 1][0], + points[points.Count() - 2] + } + ); + + // Add the end Bezier curve to the list. + beziers.Add(endBezier); + + // Return the list of Bezier curves. + return beziers; + } + + /// + /// Calculates the control points for constructing piecewise cubic Bézier curves from the given list of points and looseness factor. + /// + /// The list of points defining the path. + /// The looseness factor used to calculate control points. A higher value results in smoother curves. + /// A list of control points (pairs of Vector3) for the piecewise cubic Bézier curves. + private static List CalculateControlPoints(List points, double looseness) + { + List controlPoints = new List(); + + for (int i = 1; i < points.Count - 1; i++) + { + // Calculate the differences in x and y coordinates. + var dx = points[i - 1].X - points[i + 1].X; + var dy = points[i - 1].Y - points[i + 1].Y; + var dz = points[i - 1].Z - points[i + 1].Z; + + // Calculate the control point coordinates. + var controlPointX1 = points[i].X - dx * (1 / looseness); + var controlPointY1 = points[i].Y - dy * (1 / looseness); + var controlPointZ1 = points[i].Z - dz * (1 / looseness); + var controlPoint1 = new Vector3(controlPointX1, controlPointY1, controlPointZ1); + + var controlPointX2 = points[i].X + dx * (1 / looseness); + var controlPointY2 = points[i].Y + dy * (1 / looseness); + var controlPointZ2 = points[i].Z + dz * (1 / looseness); + var controlPoint2 = new Vector3(controlPointX2, controlPointY2, controlPointZ2); + + // Create an array to store the control points. + Vector3[] controlPointArray = new Vector3[] + { + controlPoint1, + controlPoint2 + }; + + // Add the control points to the list. + controlPoints.Add(controlPointArray); + } + + // Return the list of control points. + return controlPoints; + } } } \ No newline at end of file From 823f47d486f62249a885519842b0d844ada41234 Mon Sep 17 00:00:00 2001 From: James Bradley Date: Tue, 18 Jul 2023 14:49:54 -0400 Subject: [PATCH 12/22] init Polyline implementation of IHasArcLength --- CHANGELOG.md | 2 +- Elements/src/Geometry/Polyline.cs | 70 +++++++++++++++++++++++++++++-- Elements/test/PolylineTests.cs | 1 + 3 files changed, 69 insertions(+), 4 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 067fe8328..e6e7b8bba 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -103,7 +103,7 @@ - `Polyline` now inherits from `BoundedCurve`. - `Polyline` is now parameterized 0->length. -- `Polyline` now implements the `IHasArcLength` interface. +- `Polyline` now implements the `IHasArcLength` interface - `Arc` now inherits from `TrimmedCurve`. - `Arc` is now parameterized 0->2Pi - `Line` now inherits from `TrimmedCurve`. diff --git a/Elements/src/Geometry/Polyline.cs b/Elements/src/Geometry/Polyline.cs index 543a240b8..db4d7f617 100644 --- a/Elements/src/Geometry/Polyline.cs +++ b/Elements/src/Geometry/Polyline.cs @@ -5,6 +5,7 @@ using ClipperLib; using Elements.Search; using Elements.Geometry.Interfaces; +using Elements.Geometry.Interfaces; namespace Elements.Geometry { @@ -85,6 +86,66 @@ public virtual Line[] Segments() return SegmentsInternal(this.Vertices); } + /// + /// The mid point of the curve. + /// + /// The length based midpoint. + public virtual Vector3 MidPoint() + { + return PointAtNormalizedLength(0.5); + } + + /// + /// Returns the point on the polyline corresponding to the specified length value. + /// + /// The length value along the polyline. + /// The point on the polyline corresponding to the specified length value. + /// Thrown when the specified length is out of range. + public virtual Vector3 PointAtLength(double length) + { + double totalLength = ArcLength(this.Domain.Min, this.Domain.Max); // Calculate the total length of the Polyline + if (length < 0 || length > totalLength) + { + throw new ArgumentException("The specified length is out of range."); + } + + double accumulatedLength = 0.0; + foreach (Line segment in Segments()) + { + double segmentLength = segment.ArcLength(segment.Domain.Min, segment.Domain.Max); + + if (accumulatedLength + segmentLength >= length) + { + double remainingDistance = length - accumulatedLength; + double parameter = remainingDistance / segmentLength; + return segment.PointAtNormalized(parameter); + } + + accumulatedLength += segmentLength; + } + + // If we reach here, the desired length is equal to the total length, + // so return the end point of the Polyline. + return End; + } + + /// + /// Returns the point on the polyline corresponding to the specified normalized length-based parameter value. + /// + /// The normalized length-based parameter value, ranging from 0 to 1. + /// The point on the polyline corresponding to the specified normalized length-based parameter value. + /// Thrown when the specified parameter is out of range. + public virtual Vector3 PointAtNormalizedLength(double parameter) + { + if (parameter < 0 || parameter > 1) + { + throw new ArgumentException("The specified parameter is out of range."); + } + return PointAtLength(parameter * this.ArcLength(this.Domain.Min, this.Domain.Max)); + } + + + /// /// The mid point of the curve. /// @@ -362,7 +423,7 @@ public override Transform[] Frames(double startSetbackDistance = 0.0, // Calculate number of frames. 2 frames corresponding to end parameters. // 1 if startIndex == endIndex. - var length = endIndex - startIndex + 3; + var length = endIndex - startIndex + 3; // startIndex is set to the first distinct vertex after startParam. if (startParam.ApproximatelyEquals(startIndex)) @@ -389,7 +450,7 @@ public override Transform[] Frames(double startSetbackDistance = 0.0, result[0] = new Transform(PointAt(startParam), normals[startIndex - 1].Cross(tangent), tangent); index++; } - + for (var i = startIndex; i <= endIndex; i++, index++) { result[index] = CreateOrthogonalTransform(i, Vertices[i], normals[i]); @@ -488,6 +549,8 @@ private Transform CreateOrthogonalTransform(int i, Vector3 origin, Vector3 up) /// A list of points representing the segments. public Vector3[] DivideByLength(double divisionLength) { + var segments = new List(); + if (this.Vertices.Count < 2) { // Handle invalid polyline with insufficient vertices @@ -495,7 +558,7 @@ public Vector3[] DivideByLength(double divisionLength) } var currentProgression = 0.0; - var segments = new List { this.Vertices.FirstOrDefault() }; + segments = new List { this.Vertices.FirstOrDefault() }; foreach (var currentSegment in this.Segments()) { @@ -932,6 +995,7 @@ protected void Split(IList points, bool closed = false) var b = closed && i == this.Vertices.Count - 1 ? this.Vertices[0] : this.Vertices[i + 1]; var edge = (a, b); + // An edge may have multiple split points. // An edge may have multiple split points. // We store these in a list and sort it along the // direction of the edge, before inserting the points diff --git a/Elements/test/PolylineTests.cs b/Elements/test/PolylineTests.cs index ebf69e91a..5bf12857b 100644 --- a/Elements/test/PolylineTests.cs +++ b/Elements/test/PolylineTests.cs @@ -270,6 +270,7 @@ public void GetParameterAt() var parameter = 0.5; var testMidpoint = new Vector3(5.0, 14.560525, 0); // Midpoint var testParameterMidpoint = polyline.PointAtNormalizedLength(parameter); + var t = testParameterMidpoint.IsAlmostEqualTo(testMidpoint); Assert.True(testParameterMidpoint.IsAlmostEqualTo(testMidpoint)); var midlength = polyline.Length() * parameter; From 6c6fe1df8a475eed3ec3b528677f244cf9d5e0c1 Mon Sep 17 00:00:00 2001 From: James Bradley Date: Tue, 18 Jul 2023 17:36:40 -0400 Subject: [PATCH 13/22] Line implements Interface IHasArcLength --- CHANGELOG.md | 2 +- Elements/src/Geometry/Polyline.cs | 60 ------------------------------- Elements/test/PolylineTests.cs | 1 - 3 files changed, 1 insertion(+), 62 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index e6e7b8bba..8cd9f5aa2 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -108,7 +108,7 @@ - `Arc` is now parameterized 0->2Pi - `Line` now inherits from `TrimmedCurve`. - `Line` is now parameterized 0->length. -- `Line` now implements the `IHasArcLength` interface. +- `Line` now implements the `IHasArcLength` interface - `Bezier` now inherits from `BoundedCurve`. - `Bezier` now implements the `IHasArcLength` interface. - `Bezier.ArcLength()` now uses Gauss quadrature approximation vs linear sampling. diff --git a/Elements/src/Geometry/Polyline.cs b/Elements/src/Geometry/Polyline.cs index db4d7f617..bc9267cd6 100644 --- a/Elements/src/Geometry/Polyline.cs +++ b/Elements/src/Geometry/Polyline.cs @@ -86,66 +86,6 @@ public virtual Line[] Segments() return SegmentsInternal(this.Vertices); } - /// - /// The mid point of the curve. - /// - /// The length based midpoint. - public virtual Vector3 MidPoint() - { - return PointAtNormalizedLength(0.5); - } - - /// - /// Returns the point on the polyline corresponding to the specified length value. - /// - /// The length value along the polyline. - /// The point on the polyline corresponding to the specified length value. - /// Thrown when the specified length is out of range. - public virtual Vector3 PointAtLength(double length) - { - double totalLength = ArcLength(this.Domain.Min, this.Domain.Max); // Calculate the total length of the Polyline - if (length < 0 || length > totalLength) - { - throw new ArgumentException("The specified length is out of range."); - } - - double accumulatedLength = 0.0; - foreach (Line segment in Segments()) - { - double segmentLength = segment.ArcLength(segment.Domain.Min, segment.Domain.Max); - - if (accumulatedLength + segmentLength >= length) - { - double remainingDistance = length - accumulatedLength; - double parameter = remainingDistance / segmentLength; - return segment.PointAtNormalized(parameter); - } - - accumulatedLength += segmentLength; - } - - // If we reach here, the desired length is equal to the total length, - // so return the end point of the Polyline. - return End; - } - - /// - /// Returns the point on the polyline corresponding to the specified normalized length-based parameter value. - /// - /// The normalized length-based parameter value, ranging from 0 to 1. - /// The point on the polyline corresponding to the specified normalized length-based parameter value. - /// Thrown when the specified parameter is out of range. - public virtual Vector3 PointAtNormalizedLength(double parameter) - { - if (parameter < 0 || parameter > 1) - { - throw new ArgumentException("The specified parameter is out of range."); - } - return PointAtLength(parameter * this.ArcLength(this.Domain.Min, this.Domain.Max)); - } - - - /// /// The mid point of the curve. /// diff --git a/Elements/test/PolylineTests.cs b/Elements/test/PolylineTests.cs index 5bf12857b..ebf69e91a 100644 --- a/Elements/test/PolylineTests.cs +++ b/Elements/test/PolylineTests.cs @@ -270,7 +270,6 @@ public void GetParameterAt() var parameter = 0.5; var testMidpoint = new Vector3(5.0, 14.560525, 0); // Midpoint var testParameterMidpoint = polyline.PointAtNormalizedLength(parameter); - var t = testParameterMidpoint.IsAlmostEqualTo(testMidpoint); Assert.True(testParameterMidpoint.IsAlmostEqualTo(testMidpoint)); var midlength = polyline.Length() * parameter; From c8da721f7c035cff3e1768cf1f6b5bfcf97ead19 Mon Sep 17 00:00:00 2001 From: James Bradley Date: Tue, 18 Jul 2023 21:41:01 -0400 Subject: [PATCH 14/22] Circle implements Interface IHasArcLength --- CHANGELOG.md | 2 +- Elements/src/Geometry/Circle.cs | 77 ++++++++++++++++++ Elements/test/CircleTests.cs | 135 ++------------------------------ 3 files changed, 85 insertions(+), 129 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 8cd9f5aa2..41ae264c2 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -114,7 +114,7 @@ - `Bezier.ArcLength()` now uses Gauss quadrature approximation vs linear sampling. - `Polyline` is now parameterized 0->length. - `Circle` is now parameterized 0->2Pi. -- `Circle` now implements the `IHasArcLength` interface. +- `Circle` now implements the `IHasArcLength` interface - `Vector3.DistanceTo(Ray ray)` now returns positive infinity instead of throwing. - `Message`: removed obsolete `FromLine` method. - `AdaptiveGrid`: removed obsolete `TryGetVertexIndex` with `tolerance` parameter. diff --git a/Elements/src/Geometry/Circle.cs b/Elements/src/Geometry/Circle.cs index def74117b..533982cc0 100644 --- a/Elements/src/Geometry/Circle.cs +++ b/Elements/src/Geometry/Circle.cs @@ -224,6 +224,46 @@ private Vector3 PointAtUntransformed(double u) return new Vector3(x, y); } + /// + /// Calculates and returns the point on the circle at a specific arc length. + /// + /// The arc length along the circumference of the circle. + /// The point on the circle at the specified arc length. + public Vector3 PointAtLength(double length) + { + double parameter = (length / Circumference) * 2 * Math.PI; + return PointAt(parameter); + } + + /// + /// Calculates and returns the point on the circle at a normalized arc length. + /// + /// The normalized arc length between 0 and 1. + /// The point on the circle at the specified normalized arc length. + public Vector3 PointAtNormalizedLength(double normalizedLength) + { + double parameter = normalizedLength * 2 * Math.PI; + return PointAt(parameter); + } + + /// + /// Calculates the parameter within the range of 0 to 2π at a given point on the circle. + /// + /// The point on the circle. + /// The parameter within the range of 0 to 2π at the given point on the circle. + public double GetParameterAt(Vector3 point) + { + Vector3 relativePoint = point - Center; + + double theta = Math.Atan2(relativePoint.Y, relativePoint.X); + + if (theta < 0) + { + theta += 2 * Math.PI; + } + return theta; + } + /// /// Return transform on the arc at parameter u. /// @@ -261,5 +301,42 @@ public override double ParameterAtDistanceFromParameter(double distance, double return start + theta; } + + /// + /// Divides the circle into segments of the specified length and returns a list of points representing the division. + /// + /// The length of each segment. + /// A list of points representing the division of the circle. + public Vector3[] DivideByLength(double length) + { + List points = new List(); + double circumference = 2 * Math.PI * Radius; + int segmentCount = (int)Math.Ceiling(circumference / length); + double segmentLength = circumference / segmentCount; + + for (int i = 0; i < segmentCount; i++) + { + double parameter = i * segmentLength / circumference; + points.Add(PointAtNormalizedLength(parameter)); + } + + return points.ToArray(); + } + + /// + /// Checks if a given point lies on a circle within a specified tolerance. + /// + /// The point to be checked. + /// The circle to check against. + /// The tolerance value (optional). Default is 1E-05. + /// True if the point lies on the circle within the tolerance, otherwise false. + public static bool PointOnCircle(Vector3 point, Circle circle, double tolerance = 1E-05) + { + Vector3 centerToPoint = point - circle.Center; + double distanceToCenter = centerToPoint.Length(); + + // Check if the distance from the point to the center is within the tolerance of the circle's radius + return Math.Abs(distanceToCenter - circle.Radius) < tolerance; + } } } \ No newline at end of file diff --git a/Elements/test/CircleTests.cs b/Elements/test/CircleTests.cs index 4eb62dfda..4e423ec76 100644 --- a/Elements/test/CircleTests.cs +++ b/Elements/test/CircleTests.cs @@ -1,139 +1,18 @@ -using Elements; -using Elements.Geometry; -using Elements.Tests; using System; +using System.Collections.Generic; +using System.IO; using System.Linq; -using System.Security.Cryptography; +using Elements.Tests; +using Newtonsoft.Json; using Xunit; -using Xunit.Abstractions; namespace Elements.Geometry.Tests { - public class CircleTests + public class CircleTests : ModelTest { - [Fact] - public void CirceIntersectsCircle() + public CircleTests() { - // Planar intersecting circles - Circle c0 = new Circle(5); - Circle c1 = new Circle(new Vector3(8, 0, 0), 5); - Assert.True(c0.Intersects(c1, out var results)); - Assert.Equal(2, results.Count()); - Assert.Contains(new Vector3(4, 3, 0), results); - Assert.Contains(new Vector3(4, -3, 0), results); - c1 = new Circle(new Vector3(3, 0, 0), 5); - Assert.True(c0.Intersects(c1, out results)); - Assert.Equal(2, results.Count()); - Assert.Contains(new Vector3(1.5, 4.769696, 0), results); - Assert.Contains(new Vector3(1.5, -4.769696, 0), results); - - // Planar intersecting circles with opposite normals - c1 = new Circle(new Transform(new Vector3(8, 0, 0), Vector3.ZAxis.Negate()), 5); - Assert.True(c0.Intersects(c1, out results)); - Assert.Equal(2, results.Count()); - Assert.Contains(new Vector3(4, 3, 0), results); - Assert.Contains(new Vector3(4, -3, 0), results); - - // Planar touching circles - c1 = new Circle(new Vector3(8, 0, 0), 3); - Assert.True(c0.Intersects(c1, out results)); - Assert.Single(results); - Assert.Contains(new Vector3(5, 0, 0), results); - c1 = new Circle(new Vector3(-3, 0, 0), 2); - Assert.True(c0.Intersects(c1, out results)); - Assert.Single(results); - Assert.Contains(new Vector3(-5, 0, 0), results); - - // Planar one circle inside other - c1 = new Circle(new Vector3(1, 0, 0), 3); - Assert.False(c0.Intersects(c1, out _)); - - // Planar non-intersecting circles - c1 = new Circle(new Vector3(10, 0, 0), 3); - Assert.False(c0.Intersects(c1, out _)); - - // Non-planar intersecting circles - var t = new Transform(new Vector3(4, 0, 4), Vector3.XAxis); - c1 = new Circle(t, 5); - Assert.True(c0.Intersects(c1, out results)); - Assert.Equal(2, results.Count()); - Assert.Contains(new Vector3(4, 3, 0), results); - Assert.Contains(new Vector3(4, -3, 0), results); - - // Non-planar touching circles - t = new Transform(new Vector3(0, 7, 2), new Vector3(0, 1, -1).Unitized()); - c1 = new Circle(t, Math.Sqrt(2) * 2); - Assert.True(c0.Intersects(c1, out results)); - Assert.Single(results); - Assert.Contains(new Vector3(0, 5, 0), results); - - // Non-planar non-intersecting circles - t = new Transform(new Vector3(0, 7, 2), new Vector3(0, -1, -1).Unitized()); - c1 = new Circle(t, 1); - Assert.False(c0.Intersects(c1, out _)); - - // Same normal different origins - c1 = new Circle(new Vector3(0, 0, 10), 5); - Assert.False(c0.Intersects(c1, out _)); - - // Same circle - c1 = new Circle(5); - Assert.False(c0.Intersects(c1, out _)); - } - - [Fact] - public void CircleIntersectsLine() - { - var circleDir = new Vector3(1, 0, 1); - var t = new Transform(new Vector3(1, 1, 1), circleDir); - var circle = new Circle(t, 2); - - // Planar intersecting - var line = new InfiniteLine(new Vector3(2, 0, 0), Vector3.YAxis); - Assert.True(circle.Intersects(line, out var results)); - Assert.Equal(2, results.Count()); - Assert.Contains(new Vector3(2, 1 + Math.Sqrt(2), 0), results); - Assert.Contains(new Vector3(2, 1 - Math.Sqrt(2), 0), results); - - // Planar touching - var dir = circleDir.Cross(Vector3.YAxis); - line = new InfiniteLine(new Vector3(2, -1, 0), dir); - Assert.True(circle.Intersects(line, out results)); - Assert.Single(results); - Assert.Contains(new Vector3(1, -1, 1), results); - - // Planar non-intersecting - line = new InfiniteLine(new Vector3(3, 0, 0), dir); - Assert.False(circle.Intersects(line, out _)); - - // Non-planar touching - line = new InfiniteLine(new Vector3(1 + Math.Sqrt(2), 1, 3), Vector3.ZAxis); - Assert.True(circle.Intersects(line, out results)); - Assert.Single(results); - Assert.Contains(new Vector3(1 + Math.Sqrt(2), 1, 1 - Math.Sqrt(2)), results); - - // Non-planar non-intersecting - line = new InfiniteLine(new Vector3(1, 1, 1), Vector3.ZAxis); - Assert.False(circle.Intersects(line, out _)); - } - - [Fact] - public void AdjustRadian() - { - double reference = Units.DegreesToRadians(45); - Assert.Equal(Units.DegreesToRadians(385), Units.AdjustRadian(Units.DegreesToRadians(25), reference), 6); - Assert.Equal(Units.DegreesToRadians(60), Units.AdjustRadian(Units.DegreesToRadians(420), reference), 6); - Assert.Equal(Units.DegreesToRadians(260), Units.AdjustRadian(Units.DegreesToRadians(-100), reference), 6); - - reference = Units.DegreesToRadians(400); - Assert.Equal(Units.DegreesToRadians(745), Units.AdjustRadian(Units.DegreesToRadians(25), reference), 6); - Assert.Equal(Units.DegreesToRadians(400), Units.AdjustRadian(Units.DegreesToRadians(400), reference), 6); - Assert.Equal(Units.DegreesToRadians(620), Units.AdjustRadian(Units.DegreesToRadians(-100), reference), 6); - - reference = Units.DegreesToRadians(-160); - Assert.Equal(Units.DegreesToRadians(25), Units.AdjustRadian(Units.DegreesToRadians(25), reference), 6); - Assert.Equal(Units.DegreesToRadians(60), Units.AdjustRadian(Units.DegreesToRadians(420), reference), 6); - Assert.Equal(Units.DegreesToRadians(-100), Units.AdjustRadian(Units.DegreesToRadians(-100), reference), 6); + this.GenerateIfc = false; } [Fact, Trait("Category", "Examples")] From b45c86dbd66f5e7b4fd63d26c052e4bdf1376b6a Mon Sep 17 00:00:00 2001 From: James Bradley Date: Fri, 21 Jul 2023 15:46:54 -0400 Subject: [PATCH 15/22] Polish Bezier, IHasCurveLength --- CHANGELOG.md | 11 +- Elements/src/Geometry/Bezier.cs | 250 +++++++++++++++++++++++++++++- Elements/src/Geometry/Polyline.cs | 4 +- Elements/test/BezierTests.cs | 106 ------------- 4 files changed, 249 insertions(+), 122 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 41ae264c2..3f54d3b93 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -89,7 +89,6 @@ - `Topography.Trimmed` - `new Topography(Topography other)` - `Topography.TopMesh()` -- `UpdateElementRepresentations` flag to all serialization methods - `Bezier.PointAtLength()` - `Bezier.PointAtNormalizedLength()` - `Bezier.ParameterAt()` @@ -103,18 +102,20 @@ - `Polyline` now inherits from `BoundedCurve`. - `Polyline` is now parameterized 0->length. -- `Polyline` now implements the `IHasArcLength` interface +- `Polyline` now implements the `IHasArcLength` interface. - `Arc` now inherits from `TrimmedCurve`. -- `Arc` is now parameterized 0->2Pi +- `Arc` is now parameterized 0->2Pi. - `Line` now inherits from `TrimmedCurve`. - `Line` is now parameterized 0->length. -- `Line` now implements the `IHasArcLength` interface +- `Line` now implements the `IHasArcLength` interface. - `Bezier` now inherits from `BoundedCurve`. - `Bezier` now implements the `IHasArcLength` interface. - `Bezier.ArcLength()` now uses Gauss quadrature approximation vs linear sampling. +- `Bezier` now implements the `IHasArcLength` interface. +- `Bezier.ArcLength()` now uses Gauss quadrature approximation vs linear sampling. - `Polyline` is now parameterized 0->length. - `Circle` is now parameterized 0->2Pi. -- `Circle` now implements the `IHasArcLength` interface +- `Circle` now implements the `IHasArcLength` interface. - `Vector3.DistanceTo(Ray ray)` now returns positive infinity instead of throwing. - `Message`: removed obsolete `FromLine` method. - `AdaptiveGrid`: removed obsolete `TryGetVertexIndex` with `tolerance` parameter. diff --git a/Elements/src/Geometry/Bezier.cs b/Elements/src/Geometry/Bezier.cs index 2ba19994e..105b676c1 100644 --- a/Elements/src/Geometry/Bezier.cs +++ b/Elements/src/Geometry/Bezier.cs @@ -1,9 +1,9 @@ using System; using System.Collections.Generic; -using System.Linq; using Elements.Geometry.Interfaces; using Newtonsoft.Json; using System.Linq; +using System.Linq; namespace Elements.Geometry { @@ -35,6 +35,7 @@ public enum FrameType /// [!code-csharp[Main](../../Elements/test/BezierTests.cs?name=example)] /// public class Bezier : BoundedCurve, IHasArcLength + public class Bezier : BoundedCurve, IHasArcLength { private readonly int _lengthSamples = 500; @@ -172,7 +173,76 @@ public override double Length() /// /// Computes the arc length of the Bézier curve between the given parameter values start and end. /// https://pomax.github.io/bezierinfo/#arclength + /// Constants for Gauss quadrature points and weights (n = 24) + /// https://pomax.github.io/bezierinfo/legendre-gauss.html + /// + private static readonly double[] T = new double[] + { + -0.0640568928626056260850430826247450385909, + 0.0640568928626056260850430826247450385909, + -0.1911188674736163091586398207570696318404, + 0.1911188674736163091586398207570696318404, + -0.3150426796961633743867932913198102407864, + 0.3150426796961633743867932913198102407864, + -0.4337935076260451384870842319133497124524, + 0.4337935076260451384870842319133497124524, + -0.5454214713888395356583756172183723700107, + 0.5454214713888395356583756172183723700107, + -0.6480936519369755692524957869107476266696, + 0.6480936519369755692524957869107476266696, + -0.7401241915785543642438281030999784255232, + 0.7401241915785543642438281030999784255232, + -0.8200019859739029219539498726697452080761, + 0.8200019859739029219539498726697452080761, + -0.8864155270044010342131543419821967550873, + 0.8864155270044010342131543419821967550873, + -0.9382745520027327585236490017087214496548, + 0.9382745520027327585236490017087214496548, + -0.9747285559713094981983919930081690617411, + 0.9747285559713094981983919930081690617411, + -0.9951872199970213601799974097007368118745, + 0.9951872199970213601799974097007368118745 + }; + + /// + /// Constants for Gauss quadrature weights corresponding to the points (n = 24) + /// https://pomax.github.io/bezierinfo/legendre-gauss.html /// + private static readonly double[] C = new double[] + { + 0.1279381953467521569740561652246953718517, + 0.1279381953467521569740561652246953718517, + 0.1258374563468282961213753825111836887264, + 0.1258374563468282961213753825111836887264, + 0.121670472927803391204463153476262425607, + 0.121670472927803391204463153476262425607, + 0.1155056680537256013533444839067835598622, + 0.1155056680537256013533444839067835598622, + 0.1074442701159656347825773424466062227946, + 0.1074442701159656347825773424466062227946, + 0.0976186521041138882698806644642471544279, + 0.0976186521041138882698806644642471544279, + 0.086190161531953275917185202983742667185, + 0.086190161531953275917185202983742667185, + 0.0733464814110803057340336152531165181193, + 0.0733464814110803057340336152531165181193, + 0.0592985849154367807463677585001085845412, + 0.0592985849154367807463677585001085845412, + 0.0442774388174198061686027482113382288593, + 0.0442774388174198061686027482113382288593, + 0.0285313886289336631813078159518782864491, + 0.0285313886289336631813078159518782864491, + 0.0123412297999871995468056670700372915759, + 0.0123412297999871995468056670700372915759 + }; + + /// + /// Computes the arc length of the Bézier curve between the given parameter values start and end. + /// https://pomax.github.io/bezierinfo/#arclength + /// + /// The starting parameter value of the Bézier curve. + /// The ending parameter value of the Bézier curve. + /// The arc length between the specified parameter values. /// The starting parameter value of the Bézier curve. /// The ending parameter value of the Bézier curve. /// The arc length between the specified parameter values. @@ -183,6 +253,13 @@ public override double ArcLength(double start, double end) double sum = 0; // Accumulated sum for the arc length calculation + // Iterating through the Legendre-Gauss quadrature points and weights + for (int i = 0; i < len; i++) + double z = 0.5; // Scaling factor for the Legendre-Gauss quadrature + int len = T.Length; // Number of points in the Legendre-Gauss quadrature + + double sum = 0; // Accumulated sum for the arc length calculation + // Iterating through the Legendre-Gauss quadrature points and weights for (int i = 0; i < len; i++) { @@ -294,6 +371,15 @@ public virtual Vector3 MidPoint() return PointAtNormalizedLength(0.5); } + /// + /// The mid point of the curve. + /// + /// The length based midpoint. + public virtual Vector3 MidPoint() + { + return PointAtNormalizedLength(0.5); + } + /// /// Get the point on the curve at parameter u. /// @@ -464,6 +550,157 @@ public virtual Vector3 PointAtNormalizedLength(double parameter) return lut; } + /// + /// Returns the point on the bezier corresponding to the specified length value. + /// + /// The length value along the bezier. + /// The point on the bezier corresponding to the specified length value. + /// Thrown when the specified length is out of range. + public virtual Vector3 PointAtLength(double length) + { + double totalLength = ArcLength(this.Domain.Min, this.Domain.Max); // Calculate the total length of the Bezier + if (length < 0 || length > totalLength) + { + throw new ArgumentException("The specified length is out of range."); + } + return PointAt(ParameterAtDistanceFromParameter(length, Domain.Min)); + } + + /// + /// Returns the point on the bezier corresponding to the specified normalized length-based parameter value. + /// + /// The normalized length-based parameter value, ranging from 0 to 1. + /// The point on the bezier corresponding to the specified normalized length-based parameter value. + /// Thrown when the specified parameter is out of range. + public virtual Vector3 PointAtNormalizedLength(double parameter) + { + if (parameter < 0 || parameter > 1) + { + throw new ArgumentException("The specified parameter is out of range."); + } + return PointAtLength(parameter * this.ArcLength(this.Domain.Min, this.Domain.Max)); + } + + /// + /// Finds the parameter value on the Bezier curve that corresponds to the given 3D point within a specified threshold. + /// + /// The 3D point to find the corresponding parameter for. + /// The maximum distance threshold to consider a match between the projected point and the original point. + /// The parameter value on the Bezier curve if the distance between the projected point and the original point is within the threshold, otherwise returns null. + public double? ParameterAt(Vector3 point, double threshold = 0.0001) + { + // Find the parameter corresponding to the projected point on the Bezier curve + var parameter = ProjectedPoint(point, threshold); + + if (parameter == null) + { + // If the projected point does not return a relevant parameter return null + return null; + } + // Find the 3D point on the Bezier curve at the obtained parameter value + var projection = PointAt((double)parameter); + + // Check if the distance between the projected point and the original point is within + // a tolerence of the threshold + if (projection.DistanceTo(point) < (threshold * 10) - threshold) + { + // If the distance is within the threshold, return the parameter value + return parameter; + } + else + { + // If the distance exceeds the threshold, consider the point as not on the Bezier curve and return null + return null; + } + } + + /// + /// Projects a 3D point onto the Bezier curve to find the parameter value of the closest point on the curve. + /// + /// The 3D point to project onto the Bezier curve. + /// The maximum threshold to refine the projection and find the closest point. + /// The parameter value on the Bezier curve corresponding to the projected point, or null if the projection is not within the specified threshold. + public double? ProjectedPoint(Vector3 point, double threshold = 0.001) + { + // https://pomax.github.io/bezierinfo/#projections + // Generate a lookup table (LUT) of points and their corresponding parameter values on the Bezier curve + List<(Vector3 point, double t)> lut = GenerateLookupTable(); + + // Initialize variables to store the closest distance (d) and the index (index) of the closest point in the lookup table + double d = double.MaxValue; + int index = 0; + + // Find the closest point to the input point in the lookup table (LUT) using Euclidean distance + for (int i = 0; i < lut.Count; i++) + { + double q = Math.Sqrt(Math.Pow((point - lut[i].point).X, 2) + Math.Pow((point - lut[i].point).Y, 2) + Math.Pow((point - lut[i].point).Z, 2)); + if (q < d) + { + d = q; + index = i; + } + } + + // Obtain the parameter values of the neighboring points in the LUT for further refinement + double t1 = lut[Math.Max(index - 1, 0)].t; + double t2 = lut[Math.Min(index + 1, lut.Count - 1)].t; + double v = t2 - t1; + + // Refine the projection by iteratively narrowing down the parameter range to find the closest point + while (v > threshold) + { + // Calculate intermediate parameter values + double t0 = t1 + v / 4; + double t3 = t2 - v / 4; + + // Calculate corresponding points on the Bezier curve using the intermediate parameter values + Vector3 p0 = PointAt(t0); + Vector3 p3 = PointAt(t3); + + // Calculate the distances between the input point and the points on the Bezier curve + double d0 = Math.Sqrt(Math.Pow((point - p0).X, 2) + Math.Pow((point - p0).Y, 2) + Math.Pow((point - p0).Z, 2)); + double d3 = Math.Sqrt(Math.Pow((point - p3).X, 2) + Math.Pow((point - p3).Y, 2) + Math.Pow((point - p3).Z, 2)); + + // Choose the sub-range that is closer to the input point and update the range + if (d0 < d3) + { + t2 = t3; + } + else + { + t1 = t0; + } + + // Update the range difference for the next iteration + v = t2 - t1; + } + + // Return the average of the refined parameter values as the projection of the input point on the Bezier curve + return (t1 + t2) / 2; + } + + /// + /// Generates a lookup table of points and their corresponding parameter values on the Bezier curve. + /// + /// Number of samples to take along the curve. + /// A list of tuples containing the sampled points and their corresponding parameter values on the Bezier curve. + private List<(Vector3 point, double t)> GenerateLookupTable(int numSamples = 100) + { + // Initialize an empty list to store the lookup table (LUT) + List<(Vector3 point, double t)> lut = new List<(Vector3 point, double t)>(); + + // Generate lookup table by sampling points on the Bezier curve + for (int i = 0; i <= numSamples; i++) + { + double t = (double)i / numSamples; // Calculate the parameter value based on the current sample index + Vector3 point = PointAt(t); // Get the 3D point on the Bezier curve corresponding to the current parameter value + lut.Add((point, t)); // Add the sampled point and its corresponding parameter value to the lookup table (LUT) + } + + // Return the completed lookup table (LUT) + return lut; + } + private double BinomialCoefficient(int n, int i) { return Factorial(n) / (Factorial(i) * Factorial(n - i)); @@ -882,24 +1119,21 @@ private static List CalculateControlPoints(List points, doub // Calculate the differences in x and y coordinates. var dx = points[i - 1].X - points[i + 1].X; var dy = points[i - 1].Y - points[i + 1].Y; - var dz = points[i - 1].Z - points[i + 1].Z; // Calculate the control point coordinates. var controlPointX1 = points[i].X - dx * (1 / looseness); var controlPointY1 = points[i].Y - dy * (1 / looseness); - var controlPointZ1 = points[i].Z - dz * (1 / looseness); - var controlPoint1 = new Vector3(controlPointX1, controlPointY1, controlPointZ1); + var controlPoint1 = new Vector3(controlPointX1, controlPointY1, 0); var controlPointX2 = points[i].X + dx * (1 / looseness); var controlPointY2 = points[i].Y + dy * (1 / looseness); - var controlPointZ2 = points[i].Z + dz * (1 / looseness); - var controlPoint2 = new Vector3(controlPointX2, controlPointY2, controlPointZ2); + var controlPoint2 = new Vector3(controlPointX2, controlPointY2, 0); // Create an array to store the control points. Vector3[] controlPointArray = new Vector3[] { - controlPoint1, - controlPoint2 + controlPoint1, + controlPoint2 }; // Add the control points to the list. diff --git a/Elements/src/Geometry/Polyline.cs b/Elements/src/Geometry/Polyline.cs index bc9267cd6..3246f5577 100644 --- a/Elements/src/Geometry/Polyline.cs +++ b/Elements/src/Geometry/Polyline.cs @@ -489,8 +489,6 @@ private Transform CreateOrthogonalTransform(int i, Vector3 origin, Vector3 up) /// A list of points representing the segments. public Vector3[] DivideByLength(double divisionLength) { - var segments = new List(); - if (this.Vertices.Count < 2) { // Handle invalid polyline with insufficient vertices @@ -498,7 +496,7 @@ public Vector3[] DivideByLength(double divisionLength) } var currentProgression = 0.0; - segments = new List { this.Vertices.FirstOrDefault() }; + var segments = new List { this.Vertices.FirstOrDefault() }; foreach (var currentSegment in this.Segments()) { diff --git a/Elements/test/BezierTests.cs b/Elements/test/BezierTests.cs index d639ed100..5fc2f243b 100644 --- a/Elements/test/BezierTests.cs +++ b/Elements/test/BezierTests.cs @@ -69,111 +69,6 @@ public void Bezier_Length_OffsetFromOrigin() } [Fact] - public void IntersectsLine() - { - var a = Vector3.Origin; - var b = new Vector3(0, 5); - var c = new Vector3(5, 5); - var d = new Vector3(5, 0); - var ctrlPts = new List { a, b, c, d }; - var bezier = new Bezier(ctrlPts); - - var line = new Line(new Vector3(0, 2), new Vector3(5, 2)); - Assert.True(bezier.Intersects(line, out var results)); - Assert.Equal(2, results.Count); - Assert.All(results, r => Assert.True(r.DistanceTo(line) < Vector3.EPSILON)); - - line = new Line(new Vector3(5, 0, 5), new Vector3(5, 0, 0)); - Assert.True(bezier.Intersects(line, out results)); - Assert.Single(results); - Assert.Contains(new Vector3(5, 0), results); - - line = new Line(new Vector3(5, 5), new Vector3(0, 5)); - Assert.False(bezier.Intersects(line, out results)); - } - - [Fact] - public void IntersectsCircle() - { - var a = Vector3.Origin; - var b = new Vector3(0, 5); - var c = new Vector3(5, 5); - var d = new Vector3(5, 0); - var ctrlPts = new List { a, b, c, d }; - var bezier = new Bezier(ctrlPts); - - var arc = new Arc(new Vector3(2.5, 2.5), 2.5, 180, 270); - Assert.True(bezier.Intersects(arc, out var results)); - Assert.Single(results); - - Transform t = new Transform(new Vector3(2.5, 0), Vector3.YAxis); - arc = new Arc(new Vector3(2.5, 0), 2.5, 0, -180); - Assert.True(bezier.Intersects(arc, out results)); - Assert.Equal(2, results.Count); - Assert.Contains(new Vector3(5, 0), results); - Assert.Contains(new Vector3(0, 0), results); - } - - [Fact] - public void IntersectsEllipse() - { - var a = Vector3.Origin; - var b = new Vector3(0, 5); - var c = new Vector3(5, 5); - var d = new Vector3(5, 0); - var ctrlPts = new List { a, b, c, d }; - var bezier = new Bezier(ctrlPts); - - var arc = new EllipticalArc(new Vector3(2.5, 2), 3.5, 1, 0, 270); - Assert.True(bezier.Intersects(arc, out var results)); - Assert.Equal(3, results.Count); - - arc = new EllipticalArc(new Vector3(2.5, 0), 2.5, 2, 0, -180); - Assert.True(bezier.Intersects(arc, out results)); - Assert.Equal(2, results.Count); - Assert.Contains(new Vector3(5, 0), results); - Assert.Contains(new Vector3(0, 0), results); - } - - [Fact] - public void IntersectsPolycurve() - { - var a = Vector3.Origin; - var b = new Vector3(0, 5); - var c = new Vector3(5, 5); - var d = new Vector3(5, 0); - var ctrlPts = new List { a, b, c, d }; - var bezier = new Bezier(ctrlPts); - - var polygon = new Polygon(new Vector3[]{ - (0, 3), (6, 3), (4, 1), (-2, 1) - }); - - Assert.True(bezier.Intersects(polygon, out var results)); - Assert.Equal(4, results.Count); - Assert.Contains(new Vector3(0.93475, 3), results); - Assert.Contains(new Vector3(4.06525, 3), results); - Assert.Contains(new Vector3(4.75132, 1.75133), results); - Assert.Contains(new Vector3(0.07367, 1), results); - } - - [Fact] - public void IntersectsBezier() - { - var a = Vector3.Origin; - var b = new Vector3(0, 5); - var c = new Vector3(5, 5); - var d = new Vector3(5, 0); - var ctrlPts = new List { a, b, c, d }; - var bezier = new Bezier(ctrlPts); - - var other = new Bezier(new List { b, a, d, c }); - Assert.True(bezier.Intersects(other, out var results)); - Assert.Equal(2, results.Count); - Assert.Contains(new Vector3(0.5755, 2.5), results); - Assert.Contains(new Vector3(4.4245, 2.5), results); - } - public void Bezier_ArcLength() { var a = new Vector3(50, 150, 0); @@ -332,7 +227,6 @@ public void Split() Assert.True(normalizedBeziers[i].ControlPoints[2].Equals(testBeziers[i].ControlPoints[2])); Assert.True(normalizedBeziers[i].ControlPoints[3].Equals(testBeziers[i].ControlPoints[3])); } ->>>>>>> 374fe346(Polish Bezier, IHasCurveLength) } } } \ No newline at end of file From cc9cff7eb79ba127439002e5cd6bf7fe0a412829 Mon Sep 17 00:00:00 2001 From: James Bradley Date: Fri, 21 Jul 2023 18:05:33 -0400 Subject: [PATCH 16/22] forgot the Z --- Elements/src/Geometry/Bezier.cs | 11 +++++++---- 1 file changed, 7 insertions(+), 4 deletions(-) diff --git a/Elements/src/Geometry/Bezier.cs b/Elements/src/Geometry/Bezier.cs index 105b676c1..5d29f3a8b 100644 --- a/Elements/src/Geometry/Bezier.cs +++ b/Elements/src/Geometry/Bezier.cs @@ -1119,21 +1119,24 @@ private static List CalculateControlPoints(List points, doub // Calculate the differences in x and y coordinates. var dx = points[i - 1].X - points[i + 1].X; var dy = points[i - 1].Y - points[i + 1].Y; + var dz = points[i - 1].Z - points[i + 1].Z; // Calculate the control point coordinates. var controlPointX1 = points[i].X - dx * (1 / looseness); var controlPointY1 = points[i].Y - dy * (1 / looseness); - var controlPoint1 = new Vector3(controlPointX1, controlPointY1, 0); + var controlPointZ1 = points[i].Z - dz * (1 / looseness); + var controlPoint1 = new Vector3(controlPointX1, controlPointY1, controlPointZ1); var controlPointX2 = points[i].X + dx * (1 / looseness); var controlPointY2 = points[i].Y + dy * (1 / looseness); - var controlPoint2 = new Vector3(controlPointX2, controlPointY2, 0); + var controlPointZ2 = points[i].Z + dz * (1 / looseness); + var controlPoint2 = new Vector3(controlPointX2, controlPointY2, controlPointZ2); // Create an array to store the control points. Vector3[] controlPointArray = new Vector3[] { - controlPoint1, - controlPoint2 + controlPoint1, + controlPoint2 }; // Add the control points to the list. From 6960afbcc1577080397e608575a1567e37630af8 Mon Sep 17 00:00:00 2001 From: James Bradley Date: Mon, 13 Nov 2023 10:13:04 -0500 Subject: [PATCH 17/22] re-implement --- Elements/src/Geometry/Bezier.cs | 551 +++++++++++++++++------------- Elements/src/Geometry/Circle.cs | 237 ++++++++++--- Elements/src/Geometry/Line.cs | 75 ++-- Elements/src/Geometry/Polyline.cs | 7 +- Elements/test/BezierTests.cs | 106 ++++++ 5 files changed, 648 insertions(+), 328 deletions(-) diff --git a/Elements/src/Geometry/Bezier.cs b/Elements/src/Geometry/Bezier.cs index 5d29f3a8b..a7f93c1c2 100644 --- a/Elements/src/Geometry/Bezier.cs +++ b/Elements/src/Geometry/Bezier.cs @@ -1,9 +1,8 @@ using System; using System.Collections.Generic; +using System.Linq; using Elements.Geometry.Interfaces; using Newtonsoft.Json; -using System.Linq; -using System.Linq; namespace Elements.Geometry { @@ -34,8 +33,7 @@ public enum FrameType /// /// [!code-csharp[Main](../../Elements/test/BezierTests.cs?name=example)] /// - public class Bezier : BoundedCurve, IHasArcLength - public class Bezier : BoundedCurve, IHasArcLength + public class Bezier : BoundedCurve { private readonly int _lengthSamples = 500; @@ -138,72 +136,6 @@ public override double Length() 0.9951872199970213601799974097007368118745 }; - /// - /// Constants for Gauss quadrature weights corresponding to the points (n = 24) - /// https://pomax.github.io/bezierinfo/legendre-gauss.html - /// - private static readonly double[] C = new double[] - { - 0.1279381953467521569740561652246953718517, - 0.1279381953467521569740561652246953718517, - 0.1258374563468282961213753825111836887264, - 0.1258374563468282961213753825111836887264, - 0.121670472927803391204463153476262425607, - 0.121670472927803391204463153476262425607, - 0.1155056680537256013533444839067835598622, - 0.1155056680537256013533444839067835598622, - 0.1074442701159656347825773424466062227946, - 0.1074442701159656347825773424466062227946, - 0.0976186521041138882698806644642471544279, - 0.0976186521041138882698806644642471544279, - 0.086190161531953275917185202983742667185, - 0.086190161531953275917185202983742667185, - 0.0733464814110803057340336152531165181193, - 0.0733464814110803057340336152531165181193, - 0.0592985849154367807463677585001085845412, - 0.0592985849154367807463677585001085845412, - 0.0442774388174198061686027482113382288593, - 0.0442774388174198061686027482113382288593, - 0.0285313886289336631813078159518782864491, - 0.0285313886289336631813078159518782864491, - 0.0123412297999871995468056670700372915759, - 0.0123412297999871995468056670700372915759 - }; - - /// - /// Computes the arc length of the Bézier curve between the given parameter values start and end. - /// https://pomax.github.io/bezierinfo/#arclength - /// Constants for Gauss quadrature points and weights (n = 24) - /// https://pomax.github.io/bezierinfo/legendre-gauss.html - /// - private static readonly double[] T = new double[] - { - -0.0640568928626056260850430826247450385909, - 0.0640568928626056260850430826247450385909, - -0.1911188674736163091586398207570696318404, - 0.1911188674736163091586398207570696318404, - -0.3150426796961633743867932913198102407864, - 0.3150426796961633743867932913198102407864, - -0.4337935076260451384870842319133497124524, - 0.4337935076260451384870842319133497124524, - -0.5454214713888395356583756172183723700107, - 0.5454214713888395356583756172183723700107, - -0.6480936519369755692524957869107476266696, - 0.6480936519369755692524957869107476266696, - -0.7401241915785543642438281030999784255232, - 0.7401241915785543642438281030999784255232, - -0.8200019859739029219539498726697452080761, - 0.8200019859739029219539498726697452080761, - -0.8864155270044010342131543419821967550873, - 0.8864155270044010342131543419821967550873, - -0.9382745520027327585236490017087214496548, - 0.9382745520027327585236490017087214496548, - -0.9747285559713094981983919930081690617411, - 0.9747285559713094981983919930081690617411, - -0.9951872199970213601799974097007368118745, - 0.9951872199970213601799974097007368118745 - }; - /// /// Constants for Gauss quadrature weights corresponding to the points (n = 24) /// https://pomax.github.io/bezierinfo/legendre-gauss.html @@ -243,9 +175,6 @@ public override double Length() /// The starting parameter value of the Bézier curve. /// The ending parameter value of the Bézier curve. /// The arc length between the specified parameter values. - /// The starting parameter value of the Bézier curve. - /// The ending parameter value of the Bézier curve. - /// The arc length between the specified parameter values. public override double ArcLength(double start, double end) { double z = 0.5; // Scaling factor for the Legendre-Gauss quadrature @@ -253,13 +182,6 @@ public override double ArcLength(double start, double end) double sum = 0; // Accumulated sum for the arc length calculation - // Iterating through the Legendre-Gauss quadrature points and weights - for (int i = 0; i < len; i++) - double z = 0.5; // Scaling factor for the Legendre-Gauss quadrature - int len = T.Length; // Number of points in the Legendre-Gauss quadrature - - double sum = 0; // Accumulated sum for the arc length calculation - // Iterating through the Legendre-Gauss quadrature points and weights for (int i = 0; i < len; i++) { @@ -371,15 +293,6 @@ public virtual Vector3 MidPoint() return PointAtNormalizedLength(0.5); } - /// - /// The mid point of the curve. - /// - /// The length based midpoint. - public virtual Vector3 MidPoint() - { - return PointAtNormalizedLength(0.5); - } - /// /// Get the point on the curve at parameter u. /// @@ -550,157 +463,6 @@ public virtual Vector3 PointAtNormalizedLength(double parameter) return lut; } - /// - /// Returns the point on the bezier corresponding to the specified length value. - /// - /// The length value along the bezier. - /// The point on the bezier corresponding to the specified length value. - /// Thrown when the specified length is out of range. - public virtual Vector3 PointAtLength(double length) - { - double totalLength = ArcLength(this.Domain.Min, this.Domain.Max); // Calculate the total length of the Bezier - if (length < 0 || length > totalLength) - { - throw new ArgumentException("The specified length is out of range."); - } - return PointAt(ParameterAtDistanceFromParameter(length, Domain.Min)); - } - - /// - /// Returns the point on the bezier corresponding to the specified normalized length-based parameter value. - /// - /// The normalized length-based parameter value, ranging from 0 to 1. - /// The point on the bezier corresponding to the specified normalized length-based parameter value. - /// Thrown when the specified parameter is out of range. - public virtual Vector3 PointAtNormalizedLength(double parameter) - { - if (parameter < 0 || parameter > 1) - { - throw new ArgumentException("The specified parameter is out of range."); - } - return PointAtLength(parameter * this.ArcLength(this.Domain.Min, this.Domain.Max)); - } - - /// - /// Finds the parameter value on the Bezier curve that corresponds to the given 3D point within a specified threshold. - /// - /// The 3D point to find the corresponding parameter for. - /// The maximum distance threshold to consider a match between the projected point and the original point. - /// The parameter value on the Bezier curve if the distance between the projected point and the original point is within the threshold, otherwise returns null. - public double? ParameterAt(Vector3 point, double threshold = 0.0001) - { - // Find the parameter corresponding to the projected point on the Bezier curve - var parameter = ProjectedPoint(point, threshold); - - if (parameter == null) - { - // If the projected point does not return a relevant parameter return null - return null; - } - // Find the 3D point on the Bezier curve at the obtained parameter value - var projection = PointAt((double)parameter); - - // Check if the distance between the projected point and the original point is within - // a tolerence of the threshold - if (projection.DistanceTo(point) < (threshold * 10) - threshold) - { - // If the distance is within the threshold, return the parameter value - return parameter; - } - else - { - // If the distance exceeds the threshold, consider the point as not on the Bezier curve and return null - return null; - } - } - - /// - /// Projects a 3D point onto the Bezier curve to find the parameter value of the closest point on the curve. - /// - /// The 3D point to project onto the Bezier curve. - /// The maximum threshold to refine the projection and find the closest point. - /// The parameter value on the Bezier curve corresponding to the projected point, or null if the projection is not within the specified threshold. - public double? ProjectedPoint(Vector3 point, double threshold = 0.001) - { - // https://pomax.github.io/bezierinfo/#projections - // Generate a lookup table (LUT) of points and their corresponding parameter values on the Bezier curve - List<(Vector3 point, double t)> lut = GenerateLookupTable(); - - // Initialize variables to store the closest distance (d) and the index (index) of the closest point in the lookup table - double d = double.MaxValue; - int index = 0; - - // Find the closest point to the input point in the lookup table (LUT) using Euclidean distance - for (int i = 0; i < lut.Count; i++) - { - double q = Math.Sqrt(Math.Pow((point - lut[i].point).X, 2) + Math.Pow((point - lut[i].point).Y, 2) + Math.Pow((point - lut[i].point).Z, 2)); - if (q < d) - { - d = q; - index = i; - } - } - - // Obtain the parameter values of the neighboring points in the LUT for further refinement - double t1 = lut[Math.Max(index - 1, 0)].t; - double t2 = lut[Math.Min(index + 1, lut.Count - 1)].t; - double v = t2 - t1; - - // Refine the projection by iteratively narrowing down the parameter range to find the closest point - while (v > threshold) - { - // Calculate intermediate parameter values - double t0 = t1 + v / 4; - double t3 = t2 - v / 4; - - // Calculate corresponding points on the Bezier curve using the intermediate parameter values - Vector3 p0 = PointAt(t0); - Vector3 p3 = PointAt(t3); - - // Calculate the distances between the input point and the points on the Bezier curve - double d0 = Math.Sqrt(Math.Pow((point - p0).X, 2) + Math.Pow((point - p0).Y, 2) + Math.Pow((point - p0).Z, 2)); - double d3 = Math.Sqrt(Math.Pow((point - p3).X, 2) + Math.Pow((point - p3).Y, 2) + Math.Pow((point - p3).Z, 2)); - - // Choose the sub-range that is closer to the input point and update the range - if (d0 < d3) - { - t2 = t3; - } - else - { - t1 = t0; - } - - // Update the range difference for the next iteration - v = t2 - t1; - } - - // Return the average of the refined parameter values as the projection of the input point on the Bezier curve - return (t1 + t2) / 2; - } - - /// - /// Generates a lookup table of points and their corresponding parameter values on the Bezier curve. - /// - /// Number of samples to take along the curve. - /// A list of tuples containing the sampled points and their corresponding parameter values on the Bezier curve. - private List<(Vector3 point, double t)> GenerateLookupTable(int numSamples = 100) - { - // Initialize an empty list to store the lookup table (LUT) - List<(Vector3 point, double t)> lut = new List<(Vector3 point, double t)>(); - - // Generate lookup table by sampling points on the Bezier curve - for (int i = 0; i <= numSamples; i++) - { - double t = (double)i / numSamples; // Calculate the parameter value based on the current sample index - Vector3 point = PointAt(t); // Get the 3D point on the Bezier curve corresponding to the current parameter value - lut.Add((point, t)); // Add the sampled point and its corresponding parameter value to the lookup table (LUT) - } - - // Return the completed lookup table (LUT) - return lut; - } - private double BinomialCoefficient(int n, int i) { return Factorial(n) / (Factorial(i) * Factorial(n - i)); @@ -1146,5 +908,314 @@ private static List CalculateControlPoints(List points, doub // Return the list of control points. return controlPoints; } + + /// + public override bool Intersects(ICurve curve, out List results) + { + switch (curve) + { + case Line line: + return line.Intersects(this, out results); + case Arc arc: + return arc.Intersects(this, out results); + case EllipticalArc ellipticArc: + return ellipticArc.Intersects(this, out results); + case InfiniteLine line: + return Intersects(line, out results); + case Circle circle: + return Intersects(circle, out results); + case Ellipse ellipse: + return Intersects(ellipse, out results); + case IndexedPolycurve polycurve: + return polycurve.Intersects(this, out results); + case Bezier bezier: + return Intersects(bezier, out results); + default: + throw new NotImplementedException(); + } + } + + /// + /// Does this bezier curve intersects with an infinite line? + /// Iterative approximation is used to find intersections. + /// + /// Infinite line to intersect. + /// List containing intersection points. + /// True if intersection exists, otherwise false. + public bool Intersects(InfiniteLine line, out List results) + { + BBox3 box = new BBox3(ControlPoints); + // Bezier curve always inside it's bounding box. + // Rough check if line intersects curve box. + if (!new Line(line.Origin, line.Origin + line.Direction).Intersects( + box, out _, true)) + { + results = new List(); + return false; + } + + var l = this.Length(); + var div = (int)Math.Round(l / DefaultMinimumChordLength); + div = Math.Min(div, _lengthSamples); + + // Iteratively, find points on Bezier with 0 distance to the line. + // It Bezier was limited to 4 points - more effective approach could be used. + var roots = Equations.SolveIterative(Domain.Min, Domain.Max, div, + new Func((t) => + { + var p = PointAt(t); + return (p - p.ClosestPointOn(line)).LengthSquared(); + }), Vector3.EPSILON * Vector3.EPSILON); + + results = roots.Select(r => PointAt(r)).UniqueAverageWithinTolerance( + Vector3.EPSILON * 2).ToList(); + return results.Any(); + } + + /// + /// Does this Bezier curve intersects with a circle? + /// Iterative approximation is used to find intersections. + /// + /// Circle to intersect. + /// List containing intersection points. + /// True if any intersections exist, otherwise false. + public bool Intersects(Circle circle, out List results) + { + BBox3 box = new BBox3(ControlPoints); + // Bezier curve always inside it's bounding box. + // Rough check if curve is too far away. + var boxCenter = box.Center(); + if (circle.Center.DistanceTo(boxCenter) > + circle.Radius + (box.Max - boxCenter).Length()) + { + results = new List(); + return false; + } + + var l = this.Length(); + var div = (int)Math.Round(l / DefaultMinimumChordLength); + div = Math.Min(div, _lengthSamples); + + // Iteratively, find points on Bezier with radius distance to the circle. + var invertedT = circle.Transform.Inverted(); + var roots = Equations.SolveIterative(Domain.Min, Domain.Max, div, + new Func((t) => + { + var p = PointAt(t); + var local = invertedT.OfPoint(p); + return local.LengthSquared() - circle.Radius * circle.Radius; + }), Vector3.EPSILON * Vector3.EPSILON); + + results = roots.Select(r => PointAt(r)).UniqueAverageWithinTolerance( + Vector3.EPSILON * 2).ToList(); + return results.Any(); + } + + /// + /// Does this Bezier curve intersects with an ellipse? + /// Iterative approximation is used to find intersections. + /// + /// Ellipse to intersect. + /// List containing intersection points. + /// True if any intersections exist, otherwise false. + public bool Intersects(Ellipse ellipse, out List results) + { + BBox3 box = new BBox3(ControlPoints); + // Bezier curve always inside it's bounding box. + // Rough check if curve is too far away. + var boxCenter = box.Center(); + if (ellipse.Center.DistanceTo(boxCenter) > + Math.Max(ellipse.MajorAxis, ellipse.MinorAxis) + (box.Max - boxCenter).Length()) + { + results = new List(); + return false; + } + + var l = this.Length(); + var div = (int)Math.Round(l / DefaultMinimumChordLength); + div = Math.Min(div, _lengthSamples); + + // Iteratively, find points on ellipse with distance + // to other ellipse equal to its focal distance. + var invertedT = ellipse.Transform.Inverted(); + var roots = Equations.SolveIterative(Domain.Min, Domain.Max, div, + new Func((t) => + { + var p = PointAt(t); + var local = invertedT.OfPoint(p); + var dx = Math.Pow(local.X / ellipse.MajorAxis, 2); + var dy = Math.Pow(local.Y / ellipse.MinorAxis, 2); + return dx + dy + local.Z * local.Z - 1; + }), Vector3.EPSILON * Vector3.EPSILON); + + results = roots.Select(r => PointAt(r)).UniqueAverageWithinTolerance( + Vector3.EPSILON * 2).ToList(); + return results.Any(); + } + + /// + /// Does this Bezier curve intersects with other Bezier curve? + /// Iterative approximation is used to find intersections. + /// + /// Other Bezier curve to intersect. + /// List containing intersection points. + /// True if any intersections exist, otherwise false. + public bool Intersects(Bezier other, out List results) + { + results = new List(); + int leftSteps = BoxApproximationStepsCount(); + int rightSteps = other.BoxApproximationStepsCount(); + + var leftCache = new Dictionary(); + var rightCache = new Dictionary(); + + BBox3 box = CurveBox(leftSteps, leftCache); + BBox3 otherBox = other.CurveBox(rightSteps, rightCache); + + Intersects(other, + (box, Domain), + (otherBox, other.Domain), + leftCache, + rightCache, + ref results); + + // Subdivision algorithm produces duplicates, all tolerance away from real answer. + // Grouping and averaging them improves output as we as eliminates duplications. + results = results.UniqueAverageWithinTolerance().ToList(); + return results.Any(); + } + + private BBox3 CurveBox(int numSteps, Dictionary cache) + { + List points = new List(); + double step = Domain.Length / numSteps; + for (int i = 0; i < numSteps; i++) + { + var t = Domain.Min + i * step; + points.Add(PointAtCached(t, cache)); + } + points.Add(PointAtCached(Domain.Max, cache)); + BBox3 box = new BBox3(points); + return box; + } + + private void Intersects(Bezier other, + (BBox3 Box, Domain1d Domain) left, + (BBox3 Box, Domain1d Domain) right, + Dictionary leftCache, + Dictionary rightCache, + ref List results) + { + // If bounding boxes of two curves (not control points) are not intersect + // curves not intersect. + if (!left.Box.Intersects(right.Box)) + { + return; + } + + var leftSplit = SplitCurveBox(left, leftCache, out var loLeft, out var hiLeft); + var rightSplit = other.SplitCurveBox(right, rightCache, out var loRight, out var hiRight); + + // If boxes of two curves are tolerance sized - + // average point of their centers is treated as intersection point. + if (!leftSplit && !rightSplit) + { + results.Add(left.Box.Center().Average(right.Box.Center())); + return; + } + + // Recursively repeat process on box of subdivided curves until they are small enough. + // Each pair, which bounding boxes are not intersecting are discarded. + if (!leftSplit) + { + Intersects(other, left, loRight, leftCache, rightCache, ref results); + Intersects(other, left, hiRight, leftCache, rightCache, ref results); + } + else if (!rightSplit) + { + Intersects(other, loLeft, right, leftCache, rightCache, ref results); + Intersects(other, hiLeft, right, leftCache, rightCache, ref results); + } + else + { + Intersects(other, loLeft, loRight, leftCache, rightCache, ref results); + Intersects(other, hiLeft, loRight, leftCache, rightCache, ref results); + Intersects(other, loLeft, hiRight, leftCache, rightCache, ref results); + Intersects(other, hiLeft, hiRight, leftCache, rightCache, ref results); + } + } + + private bool SplitCurveBox((BBox3 Box, Domain1d Domain) def, + Dictionary cache, + out (BBox3 Box, Domain1d Domain) low, + out (BBox3 Box, Domain1d Domain) high) + { + low = (default, default); + high = (default, default); + + // If curve bounding box is tolerance size - it's considered as intersection. + // Otherwise calculate new boxes of two halves of the curve. + var epsilon2 = Vector3.EPSILON * Vector3.EPSILON; + var leftConvergent = (def.Box.Max - def.Box.Min).LengthSquared() < epsilon2 * 2; + if (leftConvergent) + { + return false; + } + + // If curve bounding box is tolerance size - it's considered as intersection. + // Otherwise calculate new boxes of two halves of the curve. + low = CurveBoxHalf(def, cache, true); + high = CurveBoxHalf(def, cache, false); + return true; + } + + /// + /// Get bounding box of curve segment (not control points) for half of given domain. + /// + /// Curve segment definition - box, domain and number of subdivisions. + /// Dictionary of precomputed points at parameter. + /// Take lower of higher part of curve. + /// Definition of curve segment that is half of given. + private (BBox3 Box, Domain1d Domain) CurveBoxHalf( + (BBox3 Box, Domain1d Domain) def, + Dictionary cache, + bool low) + { + var steps = BoxApproximationStepsCount(); + double step = (def.Domain.Length / 2) / steps; + var min = low ? def.Domain.Min : def.Domain.Min + def.Domain.Length / 2; + var max = low ? def.Domain.Min + def.Domain.Length / 2 : def.Domain.Max; + + List points = new List(); + for (int i = 0; i < steps; i++) + { + var t = min + i * step; + points.Add(PointAtCached(t, cache)); + } + points.Add(PointAtCached(max, cache)); + + var box = new BBox3(points); + var domain = new Domain1d(min, max); + return (box, domain); + } + + private Vector3 PointAtCached(double t, Dictionary cache) + { + if (cache.TryGetValue(t, out var p)) + { + return p; + } + else + { + p = PointAt(t); + cache.Add(t, p); + return p; + } + } + + private int BoxApproximationStepsCount() + { + return ControlPoints.Count * 8 - 1; + } } } \ No newline at end of file diff --git a/Elements/src/Geometry/Circle.cs b/Elements/src/Geometry/Circle.cs index 533982cc0..f71aca339 100644 --- a/Elements/src/Geometry/Circle.cs +++ b/Elements/src/Geometry/Circle.cs @@ -1,4 +1,3 @@ -using Elements.Validators; using System; using System.Collections.Generic; using System.Linq; @@ -11,7 +10,7 @@ namespace Elements.Geometry /// A circle. /// Parameterization of the circle is 0 -> 2PI. /// - public class Circle : Curve, IConic, IHasArcLength + public class Circle : Curve, IConic { /// The center of the circle. [JsonProperty("Center", Required = Required.AllowNull)] @@ -28,17 +27,6 @@ public Vector3 Center [System.ComponentModel.DataAnnotations.Range(0.0D, double.MaxValue)] public double Radius { get; protected set; } - /// The circumference of the circle. - [JsonIgnore] - [System.ComponentModel.DataAnnotations.Range(0.0D, double.MaxValue)] - public double Circumference { get; protected set; } - - /// - /// The domain of the curve. - /// - [JsonIgnore] - public Domain1d Domain => new Domain1d(0, 2 * Math.PI); - /// /// The coordinate system of the plane containing the circle. /// @@ -62,15 +50,7 @@ public Vector3 Normal [JsonConstructor] public Circle(Vector3 center, double radius = 1.0) { - if (!Validator.DisableValidationOnConstruction) - { - if (Math.Abs(radius - 0.0) < double.Epsilon ? true : false) - { - throw new ArgumentException($"The circle could not be created. The radius of the circle cannot be the zero: radius {radius}"); - } - } this.Radius = radius; - this.Circumference = 2 * Math.PI * this.Radius; this.Transform = new Transform(center); } @@ -80,15 +60,7 @@ public Circle(Vector3 center, double radius = 1.0) /// The radius of the circle. public Circle(double radius = 1.0) { - if (!Validator.DisableValidationOnConstruction) - { - if (Math.Abs(radius - 0.0) < double.Epsilon ? true : false) - { - throw new ArgumentException($"The circle could not be created. The radius of the circle cannot be the zero: radius {radius}"); - } - } this.Radius = radius; - this.Circumference = 2 * Math.PI * this.Radius; this.Transform = new Transform(); } @@ -97,16 +69,8 @@ public Circle(double radius = 1.0) /// public Circle(Transform transform, double radius = 1.0) { - if (!Validator.DisableValidationOnConstruction) - { - if (Math.Abs(radius - 0.0) < double.Epsilon ? true : false) - { - throw new ArgumentException($"The circle could not be created. The radius of the circle cannot be the zero: radius {radius}"); - } - } this.Transform = transform; this.Radius = radius; - this.Circumference = 2 * Math.PI * this.Radius; } /// @@ -264,6 +228,49 @@ public double GetParameterAt(Vector3 point) return theta; } + /// + /// Check if certain point is on the circle. + /// + /// Point to check. + /// Calculated parameter of point on circle. + /// True if point lays on the circle. + public bool ParameterAt(Vector3 pt, out double t) + { + var local = Transform.Inverted().OfPoint(pt); + if (local.Z.ApproximatelyEquals(0) && + local.LengthSquared().ApproximatelyEquals( + Radius * Radius, Vector3.EPSILON * Vector3.EPSILON)) + { + t = ParameterAtUntransformed(local); + return true; + } + + t = 0; + return false; + } + + /// + /// Checks if a given point lies on a circle within a specified tolerance. + /// + /// The point to be checked. + /// The circle to check against. + /// The tolerance value (optional). Default is 1E-05. + /// True if the point lies on the circle within the tolerance, otherwise false. + public static bool PointOnCircle(Vector3 point, Circle circle, double tolerance = 1E-05) + { + Vector3 centerToPoint = point - circle.Center; + double distanceToCenter = centerToPoint.Length(); + + // Check if the distance from the point to the center is within the tolerance of the circle's radius + return Math.Abs(distanceToCenter - circle.Radius) < tolerance; + } + + private double ParameterAtUntransformed(Vector3 pt) + { + var v = pt / Radius; + return Math.Atan2(v.Y, v.X); + } + /// /// Return transform on the arc at parameter u. /// @@ -323,20 +330,156 @@ public Vector3[] DivideByLength(double length) return points.ToArray(); } + /// + public override bool Intersects(ICurve curve, out List results) + { + switch (curve) + { + case BoundedCurve boundedCurve: + return boundedCurve.Intersects(this, out results); + case InfiniteLine line: + return Intersects(line, out results); + case Circle circle: + return Intersects(circle, out results); + case Ellipse elliplse: + return Intersects(elliplse, out results); + default: + throw new NotImplementedException(); + } + } + /// - /// Checks if a given point lies on a circle within a specified tolerance. + /// Does this circle intersects with other circle? + /// Circles with the same positions and radii are not considered as intersecting. /// - /// The point to be checked. - /// The circle to check against. - /// The tolerance value (optional). Default is 1E-05. - /// True if the point lies on the circle within the tolerance, otherwise false. - public static bool PointOnCircle(Vector3 point, Circle circle, double tolerance = 1E-05) + /// Other circle to intersect. + /// List containing up to two intersection points. + /// True if any intersections exist, otherwise false. + public bool Intersects(Circle other, out List results) { - Vector3 centerToPoint = point - circle.Center; - double distanceToCenter = centerToPoint.Length(); + results = new List(); - // Check if the distance from the point to the center is within the tolerance of the circle's radius - return Math.Abs(distanceToCenter - circle.Radius) < tolerance; + Plane planeA = new Plane(Center, Normal); + Plane planeB = new Plane(other.Center, other.Normal); + + // Check if two circles are on the same plane. + if (Normal.IsParallelTo(other.Normal, Vector3.EPSILON * Vector3.EPSILON) && + other.Center.DistanceTo(planeA).ApproximatelyEquals(0)) + { + var delta = other.Center - Center; + var dist = delta.Length(); + // Check if circles are on correct distance for intersection to happen. + if (dist.ApproximatelyEquals(0) || + dist > Radius + other.Radius || dist < Math.Abs(Radius - other.Radius)) + { + return false; + } + + // Build triangle with center of one circle and two intersection points. + var r1squre = Radius * Radius; + var r2squre = other.Radius * other.Radius; + var lineDist = (r1squre - r2squre + dist * dist) / (2 * dist); + var linePoint = Center + lineDist * delta.Unitized(); + double perpDistance = Math.Sqrt(r1squre - lineDist * lineDist); + // If triangle side is 0 - circles touches. Only one intersection recorded. + if (perpDistance.ApproximatelyEquals(0)) + { + results.Add(linePoint); + } + else + { + Vector3 perpDirection = delta.Cross(Normal).Unitized(); + results.Add(linePoint + perpDirection * perpDistance); + results.Add(linePoint - perpDirection * perpDistance); + } + } + // Ignore circles on parallel planes. + // Find intersection line between two planes. + else if (planeA.Intersects(planeB, out var line) && + Intersects(line, out var candidates)) + { + foreach (var item in candidates) + { + // Check each point that lays on intersection line and one of the circles. + // They are on both if they have correct distance to circle centers. + if (item.DistanceTo(other.Center).ApproximatelyEquals(other.Radius)) + { + results.Add(item); + } + } + } + + return results.Any(); + } + + /// + /// Does this circle intersects with an infinite line? + /// + /// Infinite line to intersect. + /// List containing up to two intersection points. + /// True if any intersections exist, otherwise false. + public bool Intersects(InfiniteLine line, out List results) + { + results = new List(); + + Plane circlePlane = new Plane(Center, Normal); + Vector3 closestPoint; + bool lineOnPlane = line.Origin.DistanceTo(circlePlane).ApproximatelyEquals(0) && + line.Direction.Dot(Normal).ApproximatelyEquals(0); + + // If line share a plane with circle - find closest point on it to circle center. + // If not - check if there an intersection between line and circle plane. + if (lineOnPlane) + { + closestPoint = Center.ClosestPointOn(line); + } + else if (!line.Intersects(circlePlane, out closestPoint)) + { + return false; + } + + var delta = closestPoint - Center; + var lengthSquared = delta.LengthSquared(); + var radiusSquared = Radius * Radius; + var toleranceSquared = Vector3.EPSILON * Vector3.EPSILON; + // if line not on circle plane - only one intersection is possible if it's radius away. + // this will also happen if line is on plane but only touches the circle. + if (lengthSquared.ApproximatelyEquals(radiusSquared, toleranceSquared)) + { + results.Add(closestPoint); + } + else if (lineOnPlane && lengthSquared < radiusSquared) + { + var distance = Math.Sqrt(radiusSquared - lengthSquared); + results.Add(closestPoint + line.Direction * distance); + results.Add(closestPoint - line.Direction * distance); + } + + return results.Any(); + } + + /// + /// Does this circle intersects with an ellipse? + /// Circle and ellipse that are coincides are not considered as intersecting. + /// + /// + /// Ellipse to intersect. + /// List containing up to four intersection points. + /// True if any intersections exist, otherwise false. + public bool Intersects(Ellipse ellipse, out List results) + { + return ellipse.Intersects(this, out results); + } + + /// + /// Does this circle intersects with a bounded curve? + /// + /// Curve to intersect. + /// List containing intersection points. + /// True if any intersections exist, otherwise false. + public bool Intersects(BoundedCurve curve, out List results) + { + return curve.Intersects(this, out results); } } } \ No newline at end of file diff --git a/Elements/src/Geometry/Line.cs b/Elements/src/Geometry/Line.cs index d23279878..ae193c289 100644 --- a/Elements/src/Geometry/Line.cs +++ b/Elements/src/Geometry/Line.cs @@ -4,7 +4,6 @@ using System.Linq; using Newtonsoft.Json; using Elements.Spatial; -using Elements.Geometry.Interfaces; namespace Elements.Geometry { @@ -16,7 +15,7 @@ namespace Elements.Geometry /// [!code-csharp[Main](../../Elements/test/LineTests.cs?name=example)] /// /// TODO: Rename this class to LineSegment - public class Line : TrimmedCurve, IEquatable, IHasArcLength + public class Line : TrimmedCurve, IEquatable { /// /// The domain of the curve. @@ -596,39 +595,41 @@ public static bool PointOnLine(Vector3 point, Vector3 start, Vector3 end, bool i return false; } - // /// - // /// Divide the line into as many segments of the provided length as possible. - // /// - // /// The length. - // /// A flag indicating whether segments shorter than l should be removed. - // public List DivideByLength(double l, bool removeShortSegments = false) - // { - // var len = this.Length(); - // if (l > len) - // { - // return new List() { new Line(this.Start, this.End) }; - // } - - // var total = 0.0; - // var d = this.Direction(); - // var lines = new List(); - // while (total + l <= len) - // { - // var a = this.Start + d * total; - // var b = a + d * l; - // lines.Add(new Line(a, b)); - // total += l; - // } - // if (total < len && !removeShortSegments) - // { - // var a = this.Start + d * total; - // if (!a.IsAlmostEqualTo(End)) - // { - // lines.Add(new Line(a, End)); - // } - // } - // return lines; - // } + + + /// + /// Divide the line into as many segments of the provided length as possible. + /// + /// The length. + /// A flag indicating whether segments shorter than l should be removed. + public List DivideByLength(double l, bool removeShortSegments = false) + { + var len = this.Length(); + if (l > len) + { + return new List() { new Line(this.Start, this.End) }; + } + + var total = 0.0; + var d = this.Direction(); + var lines = new List(); + while (total + l <= len) + { + var a = this.Start + d * total; + var b = a + d * l; + lines.Add(new Line(a, b)); + total += l; + } + if (total < len && !removeShortSegments) + { + var a = this.Start + d * total; + if (!a.IsAlmostEqualTo(End)) + { + lines.Add(new Line(a, End)); + } + } + return lines; + } /// /// Divides the line into segments of the specified length. @@ -671,7 +672,9 @@ public Vector3[] DivideByLength(double divisionLength) return segments.ToArray(); } - /// + /// + /// The mid point of the line. + /// public override Vector3 Mid() { return Start.Average(End); diff --git a/Elements/src/Geometry/Polyline.cs b/Elements/src/Geometry/Polyline.cs index d140ef9e0..ece888b0d 100644 --- a/Elements/src/Geometry/Polyline.cs +++ b/Elements/src/Geometry/Polyline.cs @@ -4,8 +4,6 @@ using System.Linq; using ClipperLib; using Elements.Search; -using Elements.Geometry.Interfaces; -using Elements.Geometry.Interfaces; namespace Elements.Geometry { @@ -16,7 +14,7 @@ namespace Elements.Geometry /// /// [!code-csharp[Main](../../Elements/test/PolylineTests.cs?name=example)] /// - public class Polyline : IndexedPolycurve, IHasArcLength + public class Polyline : IndexedPolycurve { /// /// The domain of the curve. @@ -919,7 +917,6 @@ protected void Split(IList points, bool closed = false) var b = closed && i == this.Vertices.Count - 1 ? this.Vertices[0] : this.Vertices[i + 1]; var edge = (a, b); - // An edge may have multiple split points. // An edge may have multiple split points. // We store these in a list and sort it along the // direction of the edge, before inserting the points @@ -1353,4 +1350,4 @@ public IndexedPolycurve Fillet(double radius) return new IndexedPolycurve(curves); } } -} +} \ No newline at end of file diff --git a/Elements/test/BezierTests.cs b/Elements/test/BezierTests.cs index 5fc2f243b..aafdb9584 100644 --- a/Elements/test/BezierTests.cs +++ b/Elements/test/BezierTests.cs @@ -68,6 +68,30 @@ public void Bezier_Length_OffsetFromOrigin() Assert.Equal(polylineLength, bezier.Length()); } + [Fact] + public void IntersectsLine() + { + var a = Vector3.Origin; + var b = new Vector3(0, 5); + var c = new Vector3(5, 5); + var d = new Vector3(5, 0); + var ctrlPts = new List { a, b, c, d }; + var bezier = new Bezier(ctrlPts); + + var line = new Line(new Vector3(0, 2), new Vector3(5, 2)); + Assert.True(bezier.Intersects(line, out var results)); + Assert.Equal(2, results.Count); + Assert.All(results, r => Assert.True(r.DistanceTo(line) < Vector3.EPSILON)); + + line = new Line(new Vector3(5, 0, 5), new Vector3(5, 0, 0)); + Assert.True(bezier.Intersects(line, out results)); + Assert.Single(results); + Assert.Contains(new Vector3(5, 0), results); + + line = new Line(new Vector3(5, 5), new Vector3(0, 5)); + Assert.False(bezier.Intersects(line, out results)); + } + [Fact] public void Bezier_ArcLength() { @@ -228,5 +252,87 @@ public void Split() Assert.True(normalizedBeziers[i].ControlPoints[3].Equals(testBeziers[i].ControlPoints[3])); } } + + [Fact] + public void IntersectsCircle() + { + var a = Vector3.Origin; + var b = new Vector3(0, 5); + var c = new Vector3(5, 5); + var d = new Vector3(5, 0); + var ctrlPts = new List { a, b, c, d }; + var bezier = new Bezier(ctrlPts); + + var arc = new Arc(new Vector3(2.5, 2.5), 2.5, 180, 270); + Assert.True(bezier.Intersects(arc, out var results)); + Assert.Single(results); + + Transform t = new Transform(new Vector3(2.5, 0), Vector3.YAxis); + arc = new Arc(new Vector3(2.5, 0), 2.5, 0, -180); + Assert.True(bezier.Intersects(arc, out results)); + Assert.Equal(2, results.Count); + Assert.Contains(new Vector3(5, 0), results); + Assert.Contains(new Vector3(0, 0), results); + } + + [Fact] + public void IntersectsEllipse() + { + var a = Vector3.Origin; + var b = new Vector3(0, 5); + var c = new Vector3(5, 5); + var d = new Vector3(5, 0); + var ctrlPts = new List { a, b, c, d }; + var bezier = new Bezier(ctrlPts); + + var arc = new EllipticalArc(new Vector3(2.5, 2), 3.5, 1, 0, 270); + Assert.True(bezier.Intersects(arc, out var results)); + Assert.Equal(3, results.Count); + + arc = new EllipticalArc(new Vector3(2.5, 0), 2.5, 2, 0, -180); + Assert.True(bezier.Intersects(arc, out results)); + Assert.Equal(2, results.Count); + Assert.Contains(new Vector3(5, 0), results); + Assert.Contains(new Vector3(0, 0), results); + } + + [Fact] + public void IntersectsPolycurve() + { + var a = Vector3.Origin; + var b = new Vector3(0, 5); + var c = new Vector3(5, 5); + var d = new Vector3(5, 0); + var ctrlPts = new List { a, b, c, d }; + var bezier = new Bezier(ctrlPts); + + var polygon = new Polygon(new Vector3[]{ + (0, 3), (6, 3), (4, 1), (-2, 1) + }); + + Assert.True(bezier.Intersects(polygon, out var results)); + Assert.Equal(4, results.Count); + Assert.Contains(new Vector3(0.93475, 3), results); + Assert.Contains(new Vector3(4.06525, 3), results); + Assert.Contains(new Vector3(4.75132, 1.75133), results); + Assert.Contains(new Vector3(0.07367, 1), results); + } + + [Fact] + public void IntersectsBezier() + { + var a = Vector3.Origin; + var b = new Vector3(0, 5); + var c = new Vector3(5, 5); + var d = new Vector3(5, 0); + var ctrlPts = new List { a, b, c, d }; + var bezier = new Bezier(ctrlPts); + + var other = new Bezier(new List { b, a, d, c }); + Assert.True(bezier.Intersects(other, out var results)); + Assert.Equal(2, results.Count); + Assert.Contains(new Vector3(0.5755, 2.5), results); + Assert.Contains(new Vector3(4.4245, 2.5), results); + } } } \ No newline at end of file From ea33a75bd5bb65154ca0ca0e69d9c1e04a7fc7ee Mon Sep 17 00:00:00 2001 From: James Bradley Date: Mon, 13 Nov 2023 10:43:26 -0500 Subject: [PATCH 18/22] adress conflicts --- Elements.MEP/src/Fittings/Reducer.cs | 4 +- Elements/src/Geometry/Circle.cs | 5 + Elements/src/Geometry/Line.cs | 6 +- Elements/test/CircleTests.cs | 166 ++++++++++++++++++++------- Elements/test/LineTests.cs | 8 +- 5 files changed, 139 insertions(+), 50 deletions(-) diff --git a/Elements.MEP/src/Fittings/Reducer.cs b/Elements.MEP/src/Fittings/Reducer.cs index 2b48a7195..a80ce92d4 100644 --- a/Elements.MEP/src/Fittings/Reducer.cs +++ b/Elements.MEP/src/Fittings/Reducer.cs @@ -36,7 +36,7 @@ public static Reducer ReducerForPipe(StraightSegment pipe, double reducerLength, var path = reducerAtEnd ? pipe.Path.Segments()[0].Reversed() : pipe.Path.Segments()[0]; - var position = path.DivideByLength(distanceFromEnd)[0].End; + var position = path.DivideByLengthToSegments(distanceFromEnd)[0].End; var orientation = path.Direction(); // var fittingMaterial = new Material("green", new Color(0, 1, 0, 0.5); @@ -133,7 +133,7 @@ public void Move(Vector3 translation) } /// - /// Port with smaller diameter points to the +X axis. + /// Port with smaller diameter points to the +X axis. /// If there is eccentric transform, the smaller part will be shifted to the -Z axis. /// We point smaller diameter in the +X direction so that there is one reducer defined in the standard orientation, to which this transformation is then applied. /// This let's us just have one size 110/90 that is rotated into a 90/110 orientation when needed. diff --git a/Elements/src/Geometry/Circle.cs b/Elements/src/Geometry/Circle.cs index f71aca339..24178c883 100644 --- a/Elements/src/Geometry/Circle.cs +++ b/Elements/src/Geometry/Circle.cs @@ -27,6 +27,11 @@ public Vector3 Center [System.ComponentModel.DataAnnotations.Range(0.0D, double.MaxValue)] public double Radius { get; protected set; } + /// The circumference of the circle. + [JsonIgnore] + [System.ComponentModel.DataAnnotations.Range(0.0D, double.MaxValue)] + public double Circumference { get; protected set; } + /// /// The coordinate system of the plane containing the circle. /// diff --git a/Elements/src/Geometry/Line.cs b/Elements/src/Geometry/Line.cs index ae193c289..a96dca17e 100644 --- a/Elements/src/Geometry/Line.cs +++ b/Elements/src/Geometry/Line.cs @@ -595,14 +595,12 @@ public static bool PointOnLine(Vector3 point, Vector3 start, Vector3 end, bool i return false; } - - /// /// Divide the line into as many segments of the provided length as possible. /// /// The length. /// A flag indicating whether segments shorter than l should be removed. - public List DivideByLength(double l, bool removeShortSegments = false) + public List DivideByLengthToSegments(double l, bool removeShortSegments = false) { var len = this.Length(); if (l > len) @@ -636,7 +634,7 @@ public List DivideByLength(double l, bool removeShortSegments = false) /// /// The desired length of each segment. /// A list of points representing the segments. - public Vector3[] DivideByLength(double divisionLength) + public Vector3[] DivideByLengthToPoints(double divisionLength) { var segments = new List(); diff --git a/Elements/test/CircleTests.cs b/Elements/test/CircleTests.cs index 7a95f569a..1ee15a3d1 100644 --- a/Elements/test/CircleTests.cs +++ b/Elements/test/CircleTests.cs @@ -15,11 +15,6 @@ public CircleTests() this.GenerateIfc = false; } - public CircleTests() - { - this.GenerateIfc = false; - } - [Fact, Trait("Category", "Examples")] public void CircleExample() { @@ -65,40 +60,6 @@ public void ZeroRadius_ThrowsException() Assert.Throws(() => new Circle(a, 0)); } - [Fact] - public void DivideIntoEqualSegments() - { - var l = new Line(Vector3.Origin, new Vector3(100, 0)); - var segments = l.DivideIntoEqualSegments(41); - var len = l.Length(); - Assert.Equal(41, segments.Count); - foreach (var s in segments) - { - Assert.Equal(s.Length(), len / 41, 5); - } - } - - [Fact] - public void DivideIntoEqualSegmentsSingle() - { - var l = new Line(Vector3.Origin, new Vector3(100, 0)); - var segments = l.DivideIntoEqualSegments(1); - Assert.Single(segments); - Assert.True(segments.First().Start.IsAlmostEqualTo(l.Start, 1e-10)); - Assert.True(segments.First().End.IsAlmostEqualTo(l.End, 1e-10)); - } - - [Fact] - public void DivideByLength() - { - var l = new Line(Vector3.Origin, new Vector3(5, 0)); - var segments = l.DivideByLength(1.1); - Assert.Equal(6, segments.Count()); - - var segments1 = l.DivideByLength(2); - Assert.Equal(4, segments1.Count()); - } - [Fact] public void GetParameterAt() { @@ -134,7 +95,7 @@ public void GetParameterAt() } [Fact] - public void PointOnLine() + public void PointOnCircle() { Circle circle = new Circle(Vector3.Origin, 5.0); @@ -142,5 +103,130 @@ public void PointOnLine() Assert.False(Circle.PointOnCircle(new Vector3(4, 0, 0), circle)); Assert.True(Circle.PointOnCircle(circle.PointAtNormalizedLength(0.5), circle)); } + + [Fact] + public void CirceIntersectsCircle() + { + // Planar intersecting circles + Circle c0 = new Circle(5); + Circle c1 = new Circle(new Vector3(8, 0, 0), 5); + Assert.True(c0.Intersects(c1, out var results)); + Assert.Equal(2, results.Count()); + Assert.Contains(new Vector3(4, 3, 0), results); + Assert.Contains(new Vector3(4, -3, 0), results); + c1 = new Circle(new Vector3(3, 0, 0), 5); + Assert.True(c0.Intersects(c1, out results)); + Assert.Equal(2, results.Count()); + Assert.Contains(new Vector3(1.5, 4.769696, 0), results); + Assert.Contains(new Vector3(1.5, -4.769696, 0), results); + + // Planar intersecting circles with opposite normals + c1 = new Circle(new Transform(new Vector3(8, 0, 0), Vector3.ZAxis.Negate()), 5); + Assert.True(c0.Intersects(c1, out results)); + Assert.Equal(2, results.Count()); + Assert.Contains(new Vector3(4, 3, 0), results); + Assert.Contains(new Vector3(4, -3, 0), results); + + // Planar touching circles + c1 = new Circle(new Vector3(8, 0, 0), 3); + Assert.True(c0.Intersects(c1, out results)); + Assert.Single(results); + Assert.Contains(new Vector3(5, 0, 0), results); + c1 = new Circle(new Vector3(-3, 0, 0), 2); + Assert.True(c0.Intersects(c1, out results)); + Assert.Single(results); + Assert.Contains(new Vector3(-5, 0, 0), results); + + // Planar one circle inside other + c1 = new Circle(new Vector3(1, 0, 0), 3); + Assert.False(c0.Intersects(c1, out _)); + + // Planar non-intersecting circles + c1 = new Circle(new Vector3(10, 0, 0), 3); + Assert.False(c0.Intersects(c1, out _)); + + // Non-planar intersecting circles + var t = new Transform(new Vector3(4, 0, 4), Vector3.XAxis); + c1 = new Circle(t, 5); + Assert.True(c0.Intersects(c1, out results)); + Assert.Equal(2, results.Count()); + Assert.Contains(new Vector3(4, 3, 0), results); + Assert.Contains(new Vector3(4, -3, 0), results); + + // Non-planar touching circles + t = new Transform(new Vector3(0, 7, 2), new Vector3(0, 1, -1).Unitized()); + c1 = new Circle(t, Math.Sqrt(2) * 2); + Assert.True(c0.Intersects(c1, out results)); + Assert.Single(results); + Assert.Contains(new Vector3(0, 5, 0), results); + + // Non-planar non-intersecting circles + t = new Transform(new Vector3(0, 7, 2), new Vector3(0, -1, -1).Unitized()); + c1 = new Circle(t, 1); + Assert.False(c0.Intersects(c1, out _)); + + // Same normal different origins + c1 = new Circle(new Vector3(0, 0, 10), 5); + Assert.False(c0.Intersects(c1, out _)); + + // Same circle + c1 = new Circle(5); + Assert.False(c0.Intersects(c1, out _)); + } + + [Fact] + public void CircleIntersectsLine() + { + var circleDir = new Vector3(1, 0, 1); + var t = new Transform(new Vector3(1, 1, 1), circleDir); + var circle = new Circle(t, 2); + + // Planar intersecting + var line = new InfiniteLine(new Vector3(2, 0, 0), Vector3.YAxis); + Assert.True(circle.Intersects(line, out var results)); + Assert.Equal(2, results.Count()); + Assert.Contains(new Vector3(2, 1 + Math.Sqrt(2), 0), results); + Assert.Contains(new Vector3(2, 1 - Math.Sqrt(2), 0), results); + + // Planar touching + var dir = circleDir.Cross(Vector3.YAxis); + line = new InfiniteLine(new Vector3(2, -1, 0), dir); + Assert.True(circle.Intersects(line, out results)); + Assert.Single(results); + Assert.Contains(new Vector3(1, -1, 1), results); + + // Planar non-intersecting + line = new InfiniteLine(new Vector3(3, 0, 0), dir); + Assert.False(circle.Intersects(line, out _)); + + // Non-planar touching + line = new InfiniteLine(new Vector3(1 + Math.Sqrt(2), 1, 3), Vector3.ZAxis); + Assert.True(circle.Intersects(line, out results)); + Assert.Single(results); + Assert.Contains(new Vector3(1 + Math.Sqrt(2), 1, 1 - Math.Sqrt(2)), results); + + // Non-planar non-intersecting + line = new InfiniteLine(new Vector3(1, 1, 1), Vector3.ZAxis); + Assert.False(circle.Intersects(line, out _)); + } + + [Fact] + public void AdjustRadian() + { + double reference = Units.DegreesToRadians(45); + Assert.Equal(Units.DegreesToRadians(385), Units.AdjustRadian(Units.DegreesToRadians(25), reference), 6); + Assert.Equal(Units.DegreesToRadians(60), Units.AdjustRadian(Units.DegreesToRadians(420), reference), 6); + Assert.Equal(Units.DegreesToRadians(260), Units.AdjustRadian(Units.DegreesToRadians(-100), reference), 6); + + reference = Units.DegreesToRadians(400); + Assert.Equal(Units.DegreesToRadians(745), Units.AdjustRadian(Units.DegreesToRadians(25), reference), 6); + Assert.Equal(Units.DegreesToRadians(400), Units.AdjustRadian(Units.DegreesToRadians(400), reference), 6); + Assert.Equal(Units.DegreesToRadians(620), Units.AdjustRadian(Units.DegreesToRadians(-100), reference), 6); + + reference = Units.DegreesToRadians(-160); + Assert.Equal(Units.DegreesToRadians(25), Units.AdjustRadian(Units.DegreesToRadians(25), reference), 6); + Assert.Equal(Units.DegreesToRadians(60), Units.AdjustRadian(Units.DegreesToRadians(420), reference), 6); + Assert.Equal(Units.DegreesToRadians(-100), Units.AdjustRadian(Units.DegreesToRadians(-100), reference), 6); + } } } diff --git a/Elements/test/LineTests.cs b/Elements/test/LineTests.cs index 1a9fad806..abf0a06d1 100644 --- a/Elements/test/LineTests.cs +++ b/Elements/test/LineTests.cs @@ -141,7 +141,7 @@ public void IntersectsQuick() public void IntersectsCircle() { Circle c = new Circle(new Vector3(5, 5, 5), 5); - + // Intersects circle at one point and touches at other. Line l = new Line(new Vector3(0, 5, 5), new Vector3(15, 5, 5)); Assert.True(l.Intersects(c, out var results)); @@ -259,13 +259,13 @@ public void DivideIntoEqualSegmentsSingle() } [Fact] - public void DivideByLength() + public void DivideByLengthToSegments() { var l = new Line(Vector3.Origin, new Vector3(5, 0)); - var segments = l.DivideByLength(1.1); + var segments = l.DivideByLengthToSegments(1.1); Assert.Equal(6, segments.Count()); - var segments1 = l.DivideByLength(2); + var segments1 = l.DivideByLengthToSegments(2); Assert.Equal(4, segments1.Count()); } From 617d4a38064c94b60e93bdbc12903a4012373a39 Mon Sep 17 00:00:00 2001 From: James Bradley Date: Mon, 13 Nov 2023 11:08:16 -0500 Subject: [PATCH 19/22] fix line divide test --- Elements/test/LineTests.cs | 11 +++++++++++ 1 file changed, 11 insertions(+) diff --git a/Elements/test/LineTests.cs b/Elements/test/LineTests.cs index abf0a06d1..b3acd1836 100644 --- a/Elements/test/LineTests.cs +++ b/Elements/test/LineTests.cs @@ -269,6 +269,17 @@ public void DivideByLengthToSegments() Assert.Equal(4, segments1.Count()); } + [Fact] + public void DivideByLengthToPoints() + { + var l = new Line(Vector3.Origin, new Vector3(5, 0)); + var segments = l.DivideByLengthToSegments(1.1); + Assert.Equal(5, segments.Count()); + + var segments1 = l.DivideByLengthToSegments(2); + Assert.Equal(3, segments1.Count()); + } + [Fact] public void DivideByLengthFromCenter() { From cb8a1e2c4d32907517f06276a76673ddd69d9dbb Mon Sep 17 00:00:00 2001 From: James Bradley Date: Mon, 13 Nov 2023 11:20:51 -0500 Subject: [PATCH 20/22] Fix circumference and tests --- Elements/src/Geometry/Circle.cs | 10 ++++++---- Elements/test/CircleTests.cs | 1 + 2 files changed, 7 insertions(+), 4 deletions(-) diff --git a/Elements/src/Geometry/Circle.cs b/Elements/src/Geometry/Circle.cs index 24178c883..610de3717 100644 --- a/Elements/src/Geometry/Circle.cs +++ b/Elements/src/Geometry/Circle.cs @@ -56,6 +56,7 @@ public Vector3 Normal public Circle(Vector3 center, double radius = 1.0) { this.Radius = radius; + this.Circumference = 2 * Math.PI * this.Radius; this.Transform = new Transform(center); } @@ -66,6 +67,7 @@ public Circle(Vector3 center, double radius = 1.0) public Circle(double radius = 1.0) { this.Radius = radius; + this.Circumference = 2 * Math.PI * this.Radius; this.Transform = new Transform(); } @@ -76,6 +78,7 @@ public Circle(Transform transform, double radius = 1.0) { this.Transform = transform; this.Radius = radius; + this.Circumference = 2 * Math.PI * this.Radius; } /// @@ -322,13 +325,12 @@ public override double ParameterAtDistanceFromParameter(double distance, double public Vector3[] DivideByLength(double length) { List points = new List(); - double circumference = 2 * Math.PI * Radius; - int segmentCount = (int)Math.Ceiling(circumference / length); - double segmentLength = circumference / segmentCount; + int segmentCount = (int)Math.Ceiling(Circumference / length); + double segmentLength = Circumference / segmentCount; for (int i = 0; i < segmentCount; i++) { - double parameter = i * segmentLength / circumference; + double parameter = i * segmentLength / Circumference; points.Add(PointAtNormalizedLength(parameter)); } diff --git a/Elements/test/CircleTests.cs b/Elements/test/CircleTests.cs index 1ee15a3d1..75a9331cf 100644 --- a/Elements/test/CircleTests.cs +++ b/Elements/test/CircleTests.cs @@ -76,6 +76,7 @@ public void GetParameterAt() Assert.True(Math.Abs(circle.GetParameterAt(mid) - Math.PI) < double.Epsilon ? true : false); + Assert.Equal(circle.Circumference, 2 * Math.PI * radius); var vector = new Vector3(3.535533, 3.535533, 0.0); var uValue = circle.GetParameterAt(vector); var expectedVector = circle.PointAt(uValue); From ef80d463c8ddf1ec160d525fe2938c02d42a814e Mon Sep 17 00:00:00 2001 From: James Bradley Date: Mon, 13 Nov 2023 11:25:13 -0500 Subject: [PATCH 21/22] fix line tests --- Elements/test/LineTests.cs | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/Elements/test/LineTests.cs b/Elements/test/LineTests.cs index b3acd1836..1daf43de8 100644 --- a/Elements/test/LineTests.cs +++ b/Elements/test/LineTests.cs @@ -259,18 +259,18 @@ public void DivideIntoEqualSegmentsSingle() } [Fact] - public void DivideByLengthToSegments() + public void DivideByLengthToPoints() { var l = new Line(Vector3.Origin, new Vector3(5, 0)); - var segments = l.DivideByLengthToSegments(1.1); + var segments = l.DivideByLengthToPoints(1.1); Assert.Equal(6, segments.Count()); - var segments1 = l.DivideByLengthToSegments(2); + var segments1 = l.DivideByLengthToPoints(2); Assert.Equal(4, segments1.Count()); } [Fact] - public void DivideByLengthToPoints() + public void DivideByLengthToPSegments() { var l = new Line(Vector3.Origin, new Vector3(5, 0)); var segments = l.DivideByLengthToSegments(1.1); From 42e6a2411d6841210233256f6f3d2eec32e10c7d Mon Sep 17 00:00:00 2001 From: James Bradley Date: Mon, 13 Nov 2023 14:58:21 -0500 Subject: [PATCH 22/22] conditional bezier arclength for quad beziers only --- Elements/src/Geometry/Bezier.cs | 28 +++++++++++++++++++++++++++- Elements/src/Geometry/Circle.cs | 22 ++++++++++++++++++++++ 2 files changed, 49 insertions(+), 1 deletion(-) diff --git a/Elements/src/Geometry/Bezier.cs b/Elements/src/Geometry/Bezier.cs index a7f93c1c2..d99e639de 100644 --- a/Elements/src/Geometry/Bezier.cs +++ b/Elements/src/Geometry/Bezier.cs @@ -104,6 +104,32 @@ public override double Length() return length; } + /// + /// Calculate the length of the bezier between start and end parameters. + /// + /// The length of the bezier between start and end. + public override double ArcLength(double start, double end) + { + if (!Domain.Includes(start, true)) + { + throw new ArgumentOutOfRangeException("start", $"The start parameter {start} must be between {Domain.Min} and {Domain.Max}."); + } + if (!Domain.Includes(end, true)) + { + throw new ArgumentOutOfRangeException("end", $"The end parameter {end} must be between {Domain.Min} and {Domain.Max}."); + } + + // if we have a true Bezier, calculate a more accurate ArcLength using Gauss weights + if (ControlPoints.Count == 2) + { + return BezierArcLength(start, end); + } + + // TODO: We use max value here so that the calculation will continue + // until at least the end of the curve. This is not a nice solution. + return ArcLengthUntil(start, double.MaxValue, out end); + } + /// /// Constants for Gauss quadrature points and weights (n = 24) /// https://pomax.github.io/bezierinfo/legendre-gauss.html @@ -175,7 +201,7 @@ public override double Length() /// The starting parameter value of the Bézier curve. /// The ending parameter value of the Bézier curve. /// The arc length between the specified parameter values. - public override double ArcLength(double start, double end) + public double BezierArcLength(double start, double end) { double z = 0.5; // Scaling factor for the Legendre-Gauss quadrature int len = T.Length; // Number of points in the Legendre-Gauss quadrature diff --git a/Elements/src/Geometry/Circle.cs b/Elements/src/Geometry/Circle.cs index 610de3717..7a75056b8 100644 --- a/Elements/src/Geometry/Circle.cs +++ b/Elements/src/Geometry/Circle.cs @@ -1,6 +1,7 @@ using System; using System.Collections.Generic; using System.Linq; +using Elements.Validators; using Elements.Geometry.Interfaces; using Newtonsoft.Json; @@ -55,6 +56,13 @@ public Vector3 Normal [JsonConstructor] public Circle(Vector3 center, double radius = 1.0) { + if (!Validator.DisableValidationOnConstruction) + { + if (Math.Abs(radius - 0.0) < double.Epsilon ? true : false) + { + throw new ArgumentException($"The circle could not be created. The radius of the circle cannot be the zero: radius {radius}"); + } + } this.Radius = radius; this.Circumference = 2 * Math.PI * this.Radius; this.Transform = new Transform(center); @@ -66,6 +74,13 @@ public Circle(Vector3 center, double radius = 1.0) /// The radius of the circle. public Circle(double radius = 1.0) { + if (!Validator.DisableValidationOnConstruction) + { + if (Math.Abs(radius - 0.0) < double.Epsilon ? true : false) + { + throw new ArgumentException($"The circle could not be created. The radius of the circle cannot be the zero: radius {radius}"); + } + } this.Radius = radius; this.Circumference = 2 * Math.PI * this.Radius; this.Transform = new Transform(); @@ -76,6 +91,13 @@ public Circle(double radius = 1.0) /// public Circle(Transform transform, double radius = 1.0) { + if (!Validator.DisableValidationOnConstruction) + { + if (Math.Abs(radius - 0.0) < double.Epsilon ? true : false) + { + throw new ArgumentException($"The circle could not be created. The radius of the circle cannot be the zero: radius {radius}"); + } + } this.Transform = transform; this.Radius = radius; this.Circumference = 2 * Math.PI * this.Radius;