diff --git a/src/GShark.Test.XUnit/Analyze/CurveTests.cs b/src/GShark.Test.XUnit/Analyze/CurveTests.cs new file mode 100644 index 00000000..278bebb3 --- /dev/null +++ b/src/GShark.Test.XUnit/Analyze/CurveTests.cs @@ -0,0 +1,123 @@ +using FluentAssertions; +using GShark.Core; +using GShark.Geometry; +using GShark.Test.XUnit.Data; +using System.Collections.Generic; +using Xunit; +using Xunit.Abstractions; + +namespace GShark.Test.XUnit.Analyze +{ + public class CurveTests + { + private readonly ITestOutputHelper _testOutput; + + public CurveTests(ITestOutputHelper testOutput) + { + _testOutput = testOutput; + } + + [Fact] + public void It_Returns_The_Approximated_Length_Of_A_Bezier() + { + // Arrange + int degree = 3; + List pts = new List + { + new Point3(0, 0, 0), + new Point3(0.5, 0, 0), + new Point3(2.5, 0, 0), + new Point3(3, 0, 0) + }; + + NurbsCurve curve = new NurbsCurve(pts, degree); + double expectedLength = 3.0; + + // Act + double curveLength = curve.Length; + + // Assert + curveLength.Should().BeApproximately(expectedLength, GSharkMath.MaxTolerance); + } + + [Fact] + public void It_Returns_Parameter_At_The_Given_Length_Of_A_Bezier() + { + // Arrange + NurbsBase curve = NurbsBaseCollection.NurbsPlanarExample(); + double[] tValuesExpected = new[] { 0, 0.122941, 0.265156, 0.420293, 0.579707, 0.734844, 0.877059, 1 }; + + int steps = 7; + double length = curve.Length / steps; + double sumLengths = 0.0; + + for (int i = 0; i < steps + 1; i++) + { + // Act + double t = curve.ParameterAtLength(sumLengths); + double segmentLength = curve.LengthAt(t); + + // Assert + t.Should().BeApproximately(tValuesExpected[i], GSharkMath.MaxTolerance); + segmentLength.Should().BeApproximately(sumLengths, GSharkMath.MaxTolerance); + + sumLengths += length; + } + } + + [Fact] + public void It_Returns_The_Length_Of_The_Curve() + { + // Arrange + NurbsBase curve = NurbsBaseCollection.NurbsPlanarExample(); + double expectedLength = 50.334675; + + // Act + double crvLength = curve.Length; + + // Assert + crvLength.Should().BeApproximately(expectedLength, GSharkMath.MinTolerance); + } + + [Theory] + [InlineData(new double[] { 5, 7, 0 }, new double[] { 5.982099, 5.950299, 0 }, 0.021824)] + [InlineData(new double[] { 12, 10, 0 }, new double[] { 11.781824, 10.364244, 0 }, 0.150707)] + [InlineData(new double[] { 21, 17, 0 }, new double[] { 21.5726, 14.101932, 0 }, 0.36828)] + [InlineData(new double[] { 32, 15, 0 }, new double[] { 31.906562, 14.36387, 0 }, 0.597924)] + [InlineData(new double[] { 41, 8, 0 }, new double[] { 42.554645, 10.750437, 0 }, 0.834548)] + [InlineData(new double[] { 50, 5, 0 }, new double[] { 50, 5, 0 }, 1.0)] + public void It_Returns_The_Closest_Point_And_Parameter(double[] ptToCheck, double[] ptExpected, double tValExpected) + { + // Arrange + NurbsBase curve = NurbsBaseCollection.NurbsPlanarExample(); + Point3 testPt = new Point3(ptToCheck[0], ptToCheck[1], ptToCheck[2]); + Point3 expectedPt = new Point3(ptExpected[0], ptExpected[1], ptExpected[2]); + + // Act + Point3 pt = curve.ClosestPoint(testPt); + double parameter = curve.ClosestParameter(testPt); + + // Assert + parameter.Should().BeApproximately(tValExpected, GSharkMath.MaxTolerance); + pt.EpsilonEquals(expectedPt, GSharkMath.MaxTolerance).Should().BeTrue(); + } + + [Theory] + [InlineData(0, 0)] + [InlineData(15, 0.278127)] + [InlineData(33, 0.672164)] + [InlineData(46, 0.928308)] + [InlineData(50.334675, 1)] + public void It_Returns_Parameter_At_The_Given_Length(double segmentLength, double tValueExpected) + { + // Arrange + NurbsBase curve = NurbsBaseCollection.NurbsPlanarExample(); + + // Act + double parameter = curve.ParameterAtLength(segmentLength); + + // Assert + parameter.Should().BeApproximately(tValueExpected, GSharkMath.MinTolerance); + } + } +} diff --git a/src/GShark.Test.XUnit/Analyze/SurfaceTests.cs b/src/GShark.Test.XUnit/Analyze/SurfaceTests.cs new file mode 100644 index 00000000..5c7fb3a8 --- /dev/null +++ b/src/GShark.Test.XUnit/Analyze/SurfaceTests.cs @@ -0,0 +1,71 @@ +using FluentAssertions; +using GShark.Core; +using GShark.Enumerations; +using GShark.Geometry; +using GShark.Test.XUnit.Data; +using Xunit; +using Xunit.Abstractions; + +namespace GShark.Test.XUnit.Analyze +{ + public class SurfaceTests + { + private readonly ITestOutputHelper _testOutput; + + public SurfaceTests(ITestOutputHelper testOutput) + { + _testOutput = testOutput; + } + + [Theory] + [InlineData(0.204157623157292, 0.716170472509343, new double[] { 2.5, 7, 5 })] + [InlineData(0.237211551442712, 0.154628316784507, new double[] { 2.5, 1.5, 2 })] + [InlineData(0.910119163727208, 0.229417610613794, new double[] { 9, 2.5, 1 })] + [InlineData(0.50870054333679, 0.360138133269618, new double[] { 5, 5, 1 })] + public void It_Returns_Parameter_U_V_Of_A_Closest_Point(double u, double v, double[] testPt) + { + // Arrange + NurbsSurface surface = NurbsSurfaceCollection.SurfaceFromPoints(); + Point3 pt = new Point3(testPt[0], testPt[1], testPt[2]); + (double u, double v) expectedUV = (u, v); + + // Act + var closestParameter = surface.ClosestParameter(pt); + + // Assert + (closestParameter.U - expectedUV.u).Should().BeLessThan(GSharkMath.MaxTolerance); + (closestParameter.V - expectedUV.v).Should().BeLessThan(GSharkMath.MaxTolerance); + } + + [Fact] + public void Returns_The_Surface_Isocurve_At_U_Direction() + { + // Arrange + NurbsSurface surface = NurbsSurfaceCollection.SurfaceFromPoints(); + Point3 expectedPt = new Point3(3.591549, 10, 4.464789); + + // Act + NurbsBase Isocurve = surface.IsoCurve(0.3, SurfaceDirection.U); + + // Assert + Isocurve.ControlPointLocations[1].DistanceTo(expectedPt).Should().BeLessThan(GSharkMath.MinTolerance); + } + + [Fact] + public void Returns_The_Surface_Isocurve_At_V_Direction() + { + // Arrange + NurbsSurface surface = NurbsSurfaceCollection.SurfaceFromPoints(); + Point3 expectedPt = new Point3(5, 4.615385, 2.307692); + Point3 expectedPtAt = new Point3(5, 3.913043, 1.695652); + + // Act + NurbsBase Isocurve = surface.IsoCurve(0.3, SurfaceDirection.V); + Point3 ptAt = Isocurve.PointAt(0.5); + + // Assert + Isocurve.ControlPointLocations[1].DistanceTo(expectedPt).Should().BeLessThan(GSharkMath.MinTolerance); + ptAt.DistanceTo(expectedPtAt).Should().BeLessThan(GSharkMath.MinTolerance); + } + } +} diff --git a/src/GShark.Test.XUnit/Core/TrigonometryTests.cs b/src/GShark.Test.XUnit/Core/TrigonometryTests.cs index 2c8d63cd..18910c97 100644 --- a/src/GShark.Test.XUnit/Core/TrigonometryTests.cs +++ b/src/GShark.Test.XUnit/Core/TrigonometryTests.cs @@ -62,7 +62,7 @@ public void It_Returns_The_Closest_Point_On_A_Segment(double[] ptToCheck, double Point3 pt1 = new Point3(10, 10, 0); // Act - (double tValue, Point3 pt) closestPt = Analyze.ClosestPointToSegment(testPt, pt0, pt1, 0, 1); + (double tValue, Point3 pt) closestPt = Trigonometry.ClosestPointToSegment(testPt, pt0, pt1, 0, 1); // Assert closestPt.tValue.Should().BeApproximately(tValExpected, GSharkMath.MaxTolerance); diff --git a/src/GShark.Test.XUnit/Operation/IntersectionTests.cs b/src/GShark.Test.XUnit/Intersection/IntersectionTests.cs similarity index 99% rename from src/GShark.Test.XUnit/Operation/IntersectionTests.cs rename to src/GShark.Test.XUnit/Intersection/IntersectionTests.cs index 14720e26..de555325 100644 --- a/src/GShark.Test.XUnit/Operation/IntersectionTests.cs +++ b/src/GShark.Test.XUnit/Intersection/IntersectionTests.cs @@ -1,13 +1,12 @@ using FluentAssertions; using GShark.Core; using GShark.Geometry; -using GShark.Operation; -using System.Collections.Generic; using GShark.Intersection; +using System.Collections.Generic; using Xunit; using Xunit.Abstractions; -namespace GShark.Test.XUnit.Operation +namespace GShark.Test.XUnit.Intersection { public class IntersectionTests { diff --git a/src/GShark.Test.XUnit/Operation/AnalyzeTests.cs b/src/GShark.Test.XUnit/Operation/AnalyzeTests.cs deleted file mode 100644 index 47a7e883..00000000 --- a/src/GShark.Test.XUnit/Operation/AnalyzeTests.cs +++ /dev/null @@ -1,176 +0,0 @@ -using FluentAssertions; -using GShark.Core; -using GShark.Geometry; -using GShark.Operation; -using GShark.Test.XUnit.Data; -using System.Collections.Generic; -using GShark.Enumerations; -using Xunit; -using Xunit.Abstractions; - -namespace GShark.Test.XUnit.Operation -{ - public class AnalyzeTests - { - private readonly ITestOutputHelper _testOutput; - - public AnalyzeTests(ITestOutputHelper testOutput) - { - _testOutput = testOutput; - } - - [Fact] - public void RationalBezierCurveArcLength_Returns_The_Approximated_Length() - { - // Arrange - int degree = 3; - List pts = new List - { - new Point3(0, 0, 0), - new Point3(0.5, 0, 0), - new Point3(2.5, 0, 0), - new Point3(3, 0, 0) - }; - - NurbsCurve curve = new NurbsCurve(pts, degree); - double expectedLength = 3.0; - - // Act - double curveLength = Analyze.BezierCurveLength(curve, 1); - - // Assert - curveLength.Should().BeApproximately(expectedLength, GSharkMath.MaxTolerance); - } - - [Fact] - public void RationalBezierCurveParamAtLength_Returns_Parameters_At_Passed_Lengths() - { - // Arrange - NurbsBase curve = NurbsBaseCollection.NurbsPlanarExample(); - double[] tValuesExpected = new[] { 0, 0.122941, 0.265156, 0.420293, 0.579707, 0.734844, 0.877059, 1 }; - - int steps = 7; - double length = curve.Length / steps; - double sumLengths = 0.0; - - for (int i = 0; i < steps + 1; i++) - { - // Act - double t = Analyze.BezierCurveParamAtLength(curve, sumLengths, GSharkMath.MaxTolerance); - - double segmentLength = Analyze.BezierCurveLength(curve, t); - - // Assert - t.Should().BeApproximately(tValuesExpected[i], GSharkMath.MaxTolerance); - segmentLength.Should().BeApproximately(sumLengths, GSharkMath.MaxTolerance); - - sumLengths += length; - } - } - - [Fact] - public void It_Returns_The_Length_Of_The_Curve() - { - // Arrange - NurbsBase curve = NurbsBaseCollection.NurbsPlanarExample(); - double expectedLength = 50.334675; - - // Act - double crvLength = Analyze.CurveLength(curve); - - // Assert - crvLength.Should().BeApproximately(expectedLength, GSharkMath.MinTolerance); - } - - [Theory] - [InlineData(new double[] { 5, 7, 0 }, new double[] { 5.982099, 5.950299, 0 }, 0.021824)] - [InlineData(new double[] { 12, 10, 0 }, new double[] { 11.781824, 10.364244, 0 }, 0.150707)] - [InlineData(new double[] { 21, 17, 0 }, new double[] { 21.5726, 14.101932, 0 }, 0.36828)] - [InlineData(new double[] { 32, 15, 0 }, new double[] { 31.906562, 14.36387, 0 }, 0.597924)] - [InlineData(new double[] { 41, 8, 0 }, new double[] { 42.554645, 10.750437, 0 }, 0.834548)] - [InlineData(new double[] { 50, 5, 0 }, new double[] { 50, 5, 0 }, 1.0)] - public void It_Returns_The_Closest_Point_And_The_Parameter_t(double[] ptToCheck, double[] ptExpected, double tValExpected) - { - // Arrange - NurbsBase curve = NurbsBaseCollection.NurbsPlanarExample(); - Point3 testPt = new Point3(ptToCheck[0], ptToCheck[1], ptToCheck[2]); - Point3 expectedPt = new Point3(ptExpected[0], ptExpected[1], ptExpected[2]); - - // Act - var pt = Analyze.CurveClosestPoint(curve, testPt, out double t); - - // Assert - t.Should().BeApproximately(tValExpected, GSharkMath.MaxTolerance); - pt.EpsilonEquals(expectedPt, GSharkMath.MaxTolerance).Should().BeTrue(); - } - - [Theory] - [InlineData(0, 0)] - [InlineData(15, 0.278127)] - [InlineData(33, 0.672164)] - [InlineData(46, 0.928308)] - [InlineData(50.334675, 1)] - public void RationalCurveParameterAtLength_Returns_Parameter_t_At_The_Given_Length(double segmentLength, double tValueExpected) - { - // Arrange - NurbsBase curve = NurbsBaseCollection.NurbsPlanarExample(); - - // Act - double t = Analyze.CurveParameterAtLength(curve, segmentLength); - - // Assert - t.Should().BeApproximately(tValueExpected, 1e-5); - } - - [Theory] - [InlineData(0.204157623157292, 0.716170472509343, new double[] { 2.5, 7, 5 })] - [InlineData(0.237211551442712, 0.154628316784507, new double[] { 2.5, 1.5, 2 })] - [InlineData(0.910119163727208, 0.229417610613794, new double[] { 9, 2.5, 1 })] - [InlineData(0.50870054333679, 0.360138133269618, new double[] { 5, 5, 1 })] - public void RationalSurfaceClosestParam_Returns_Parameters_U_V_Of_A_Closest_Point(double u, double v, double[] testPt) - { - // Arrange - NurbsSurface surface = NurbsSurfaceCollection.SurfaceFromPoints(); - Point3 pt = new Point3(testPt[0], testPt[1], testPt[2]); - (double u, double v) expectedUV = (u, v); - - // Act - var closestParameter = Analyze.SurfaceClosestParameter(surface, pt); - - // Assert - (closestParameter.u - expectedUV.u).Should().BeLessThan(GSharkMath.MaxTolerance); - (closestParameter.v - expectedUV.v).Should().BeLessThan(GSharkMath.MaxTolerance); - } - - [Fact] - public void Returns_The_Surface_Isocurve_At_U_Direction() - { - // Arrange - NurbsSurface surface = NurbsSurfaceCollection.SurfaceFromPoints(); - Point3 expectedPt = new Point3(3.591549, 10, 4.464789); - - // Act - NurbsBase Isocurve = Analyze.Isocurve(surface, 0.3, SurfaceDirection.U); - - // Assert - Isocurve.ControlPointLocations[1].DistanceTo(expectedPt).Should().BeLessThan(GSharkMath.MinTolerance); - } - - [Fact] - public void Returns_The_Surface_Isocurve_At_V_Direction() - { - // Arrange - NurbsSurface surface = NurbsSurfaceCollection.SurfaceFromPoints(); - Point3 expectedPt = new Point3(5, 4.615385, 2.307692); - Point3 expectedPtAt = new Point3(5, 3.913043, 1.695652); - - // Act - NurbsBase Isocurve = Analyze.Isocurve(surface, 0.3, SurfaceDirection.V); - Point3 ptAt = Isocurve.PointAt(0.5); - - // Assert - Isocurve.ControlPointLocations[1].DistanceTo(expectedPt).Should().BeLessThan(GSharkMath.MinTolerance); - ptAt.DistanceTo(expectedPtAt).Should().BeLessThan(GSharkMath.MinTolerance); - } - } -} diff --git a/src/GShark.Test.XUnit/Operation/CheckTests.cs b/src/GShark.Test.XUnit/Operation/CheckTests.cs deleted file mode 100644 index 3a711177..00000000 --- a/src/GShark.Test.XUnit/Operation/CheckTests.cs +++ /dev/null @@ -1,6 +0,0 @@ -namespace GShark.Test.XUnit.Operation -{ - public class CheckTests - { - } -} diff --git a/src/GShark/Analyze/Curve.cs b/src/GShark/Analyze/Curve.cs new file mode 100644 index 00000000..6e7a6352 --- /dev/null +++ b/src/GShark/Analyze/Curve.cs @@ -0,0 +1,294 @@ +using GShark.Core; +using GShark.Geometry; +using GShark.Operation; +using System; +using System.Collections.Generic; +using System.Linq; + +namespace GShark.Analyze +{ + /// + /// Contains methods for analyzing curves. + /// + internal static class Curve + { + /// + /// Computes the approximate length of a rational curve by gaussian quadrature - assumes a smooth curve. + /// + /// The curve object. + /// The parameter at which to approximate the length. + /// + /// The degree of gaussian quadrature to perform. + /// A higher number yields a more exact result, default set to 17. + /// + /// The approximate length. + internal static double Length(NurbsBase curve, double u = -1.0, int gaussDegIncrease = 17) + { + double uSet = u < 0.0 ? curve.Knots.Last() : u; + + List crvs = Modify.Curve.DecomposeIntoBeziers(curve); + double sum = 0.0; + + foreach (NurbsBase bezier in crvs) + { + if (!(bezier.Knots[0] + GSharkMath.Epsilon < uSet)) + { + break; + } + double param = Math.Min(bezier.Knots.Last(), uSet); + sum += BezierLength(bezier, param, gaussDegIncrease); + } + + return sum; + } + + /// + /// Computes the curvature vector. + /// The vector has length equal to the radius of the curvature circle and with direction to the center of the circle. + /// https://github.com/mcneel/opennurbs/blob/484ba88836bbedff8fe0b9e574fcd6434b49c21c/opennurbs_math.cpp#L839 + /// + /// The first derivative. + /// The second derivative. + /// The curvature vector. + internal static Vector3 Curvature(Vector3 derivative1, Vector3 derivative2) + { + double d1 = derivative1.Length; + + // If the first derivative exists and is zero the curvature is a zero vector. + if (d1 == 0.0) + { + return Vector3.Zero; + } + + Vector3 tangent = derivative1 / d1; + double d2DotTang = (derivative2 * -1) * tangent; + d1 = 1.0 / (d1 * d1); + Vector3 curvature = d1 * (derivative2 + d2DotTang * tangent); // usually identified as k. + + double curvatureLength = curvature.Length; + if (curvatureLength < 1.490116119385E-08) // SqrtEpsilon value that is used when comparing square roots. + { + throw new Exception("Curvature is infinite."); + } + + double radius = (curvature / (curvatureLength * curvatureLength)).Length; + return curvature.Unitize().Amplify(radius); + } + + /// + /// Computes the approximate length of a rational bezier curve by gaussian quadrature - assumes a smooth curve. + /// + /// The curve object. + /// The parameter at which to approximate the length. + /// + /// the degree of gaussian quadrature to perform. + /// A higher number yields a more exact result, default set to 17. + /// + /// The approximate length of a bezier. + internal static double BezierLength(NurbsBase curve, double u = -1.0, int gaussDegIncrease = 17) + { + double uSet = u < 0.0 ? curve.Knots.Last() : u; + double z = (uSet - curve.Knots[0]) / 2; + double sum = 0.0; + int gaussDegree = curve.Degree + gaussDegIncrease; + + for (int i = 0; i < gaussDegree; i++) + { + double cu = z * LegendreGaussData.tValues[gaussDegree][i] + z + curve.Knots[0]; + List tan = Evaluation.RationalCurveDerivatives(curve, cu); + + sum += LegendreGaussData.cValues[gaussDegree][i] * tan[1].Length; + } + + return z * sum; + } + + /// + /// Computes the curve parameter at a given length. + /// + /// The curve object. + /// The length to find the parameter. + /// If set less or equal 0.0, the tolerance used is 1e-10. + /// The parameter at the given length. + internal static double BezierParameterAtLength(NurbsBase curve, double segmentLength, double tolerance) + { + if (segmentLength < 0) + { + return curve.Knots[0]; + } + + // We compute the whole length. + double curveLength = BezierLength(curve); + + if (segmentLength > curveLength) + { + return curve.Knots[curve.Knots.Count - 1]; + } + + // Divide and conquer. + double setTolerance = tolerance <= 0.0 ? GSharkMath.Epsilon : tolerance; + + double startT = curve.Knots[0]; + double startLength = 0.0; + + double endT = curve.Knots[curve.Knots.Count - 1]; + double endLength = curveLength; + + while (endLength - startLength > setTolerance) + { + double midT = (startT + endT) / 2; + double midLength = BezierLength(curve, midT); + + if (midLength > segmentLength) + { + endT = midT; + endLength = midLength; + } + else + { + startT = midT; + startLength = midLength; + } + } + + return (startT + endT) / 2; + } + + /// + /// Computes the closest parameters on a curve to a given point.
+ /// Corresponds to page 244 chapter six from The NURBS Book by Piegl and Tiller. + /// + /// The curve object. + /// Point to search from. + /// The closest parameter on the curve. + internal static double ClosestParameter(NurbsBase curve, Point3 point) + { + double minimumDistance = double.PositiveInfinity; + double tParameter = 0D; + List ctrlPts = curve.ControlPointLocations; + + (List tValues, List pts) = Sampling.Curve.RegularSample(curve, ctrlPts.Count * curve.Degree); + + for (int i = 0; i < pts.Count - 1; i++) + { + double t0 = tValues[i]; + double t1 = tValues[i + 1]; + + Point3 pt0 = pts[i]; + Point3 pt1 = pts[i + 1]; + + var (tValue, pt) = Trigonometry.ClosestPointToSegment(point, pt0, pt1, t0, t1); + double distance = pt.DistanceTo(point); + + if (!(distance < minimumDistance)) continue; + minimumDistance = distance; + tParameter = tValue; + } + + int maxIterations = 5; + int j = 0; + // Two zero tolerances can be used to indicate convergence: + double tol1 = GSharkMath.MaxTolerance; // a measure of Euclidean distance; + double tol2 = 0.0005; // a zero cosine measure. + double tVal0 = curve.Knots[0]; + double tVal1 = curve.Knots[curve.Knots.Count - 1]; + bool closedCurve = (ctrlPts[0] - ctrlPts[ctrlPts.Count - 1]).SquareLength < GSharkMath.Epsilon; + double Cu = tParameter; + + // To avoid infinite loop we limited the interaction. + while (j < maxIterations) + { + List e = Evaluation.RationalCurveDerivatives(curve, Cu, 2); + Vector3 diff = e[0] - new Vector3(point); // C(u) - P + + + // First condition, point coincidence: + // |C(u) - p| < e1 + double c1v = diff.Length; + bool c1 = c1v <= tol1; + + // Second condition, zero cosine: + // C'(u) * (C(u) - P) + // ------------------ < e2 + // |C'(u)| |C(u) - P| + double c2n = Vector3.DotProduct(e[1], diff); + double c2d = (e[1] * c1v).Length; + double c2v = c2n / c2d; + bool c2 = Math.Abs(c2v) <= tol2; + + // If at least one of these conditions is not satisfied, + // a new value, ui+l> is computed using the NewtonIteration. + // Then two more conditions are checked. + if (c1 && c2) return Cu; + double ct = NewtonIteration(Cu, e, diff); + + // Ensure that the parameter stays within the boundary of the curve. + if (ct < tVal0) ct = closedCurve ? tVal1 - (ct - tVal0) : tVal0; + if (ct > tVal1) ct = closedCurve ? tVal0 + (ct - tVal1) : tVal1; + + // the parameter does not change significantly, the point is off the end of the curve. + double c3v = (e[1] * (ct - Cu)).Length; + if (c3v < tol1) return Cu; + + Cu = ct; + j++; + } + + return Cu; + } + + /// + /// Newton iteration to minimize the distance between a point and a curve. + /// + /// The parameter obtained at the ith Newton iteration. + /// Point on curve identify as C'(u) + /// Representing the difference from C(u) - P. + /// The minimized parameter. + private static double NewtonIteration(double u, List derivativePts, Vector3 difference) + { + // The distance from P to C(u) is minimum when f(u) = 0, whether P is on the curve or not. + // C'(u) * ( C(u) - P ) = 0 = f(u) + // C(u) is the curve, p is the point, * is a dot product + double f = Vector3.DotProduct(derivativePts[1], difference); + + // f' = C"(u) * ( C(u) - p ) + C'(u) * C'(u) + double s0 = Vector3.DotProduct(derivativePts[2], difference); + double s1 = Vector3.DotProduct(derivativePts[1], derivativePts[1]); + double df = s0 + s1; + + return u - f / df; + } + + /// + /// Approximates the parameter at a given length on a curve. + /// + /// The curve object. + /// The arc length for which to do the procedure. + /// If set less or equal 0.0, the tolerance used is 1e-10. + /// The parameter on the curve. + internal static double ParameterAtLength(NurbsBase curve, double segmentLength, double tolerance = -1) + { + if (segmentLength < GSharkMath.Epsilon) return curve.Knots[0]; + if (Math.Abs(curve.Length - segmentLength) < GSharkMath.Epsilon) return curve.Knots[curve.Knots.Count - 1]; + + List curves = Modify.Curve.DecomposeIntoBeziers(curve); + int i = 0; + double curveLength = -GSharkMath.Epsilon; + double segmentLengthLeft = segmentLength; + + // Iterate through the curves consuming the bezier's, summing their length along the way. + while (curveLength < segmentLength && i < curves.Count) + { + double bezierLength = BezierLength(curves[i]); + curveLength += bezierLength; + + if (segmentLength < curveLength + GSharkMath.Epsilon) + return BezierParameterAtLength(curves[i], segmentLengthLeft, tolerance); + i++; + segmentLengthLeft -= bezierLength; + } + + return -1; + } + } +} \ No newline at end of file diff --git a/src/GShark/Analyze/Surface.cs b/src/GShark/Analyze/Surface.cs new file mode 100644 index 00000000..8a8ddc88 --- /dev/null +++ b/src/GShark/Analyze/Surface.cs @@ -0,0 +1,252 @@ +using GShark.Core; +using GShark.Enumerations; +using GShark.Geometry; +using GShark.Operation; +using System; +using System.Collections.Generic; +using System.Linq; + +namespace GShark.Analyze +{ + /// + /// Contains methods for analyzing surfaces. + /// + internal static class Surface + { + /// + /// Computes the closest parameters on a surface to a given point.
+ /// Corresponds to page 244 chapter six from The NURBS Book by Piegl and Tiller. + /// + /// The surface object. + /// Point to search from. + /// The closest parameter on the surface. + internal static (double u, double v) ClosestParameter(NurbsSurface surface, Point3 point) + { + double minimumDistance = double.PositiveInfinity; + (double u, double v) selectedUV = (0.5, 0.5); + NurbsSurface splitSrf = surface; + double param = 0.5; + int maxIterations = 5; + + for (int i = 0; i < maxIterations; i++) + { + NurbsSurface[] surfaces = splitSrf.SplitAt(0.5, SplitDirection.Both); + Point3[] pts = surfaces.Select(s => s.PointAt(0.5, 0.5)).ToArray(); + double[] distanceBetweenPts = pts.Select(point.DistanceTo).ToArray(); + if (distanceBetweenPts.All(d => d > minimumDistance)) break; + + (double, double)[] srfUV = DefiningUV(selectedUV, param); + + for (int j = 0; j < distanceBetweenPts.Length; j++) + { + if (!(distanceBetweenPts[j] < minimumDistance)) continue; + minimumDistance = distanceBetweenPts[j]; + selectedUV = srfUV[j]; + splitSrf = surfaces[j]; + } + + param *= 0.5; + } + + int t = 0; + // Two zero tolerances can be used to indicate convergence: + double tol1 = GSharkMath.MaxTolerance; // a measure of Euclidean distance; + double tol2 = 0.0005; // a zero cosine measure. + double minU = surface.KnotsU[0]; + double maxU = surface.KnotsU.Last(); + double minV = surface.KnotsV[0]; + double maxV = surface.KnotsV.Last(); + bool closedDirectionU = surface.IsClosed(SurfaceDirection.U); + bool closedDirectionV = surface.IsClosed(SurfaceDirection.V); + + // To avoid infinite loop we limited the interaction. + while (t < maxIterations) + { + // Get derivatives. + var eval = Evaluation.RationalSurfaceDerivatives(surface, selectedUV.u, selectedUV.v, 2); + + // Convergence criteria: + // First condition, point coincidence: + // |S(u,v) - p| < e1 + Vector3 diff = eval[0, 0] - new Vector3(point); + double c1v = diff.Length; + bool c1 = c1v <= tol1; + + // Second condition, zero cosine: + // |Su(u,v)*(S(u,v) - P)| + // ---------------------- < e2 + // |Su(u,v)| |S(u,v) - P| + // + // |Sv(u,v)*(S(u,v) - P)| + // ---------------------- < e2 + // |Sv(u,v)| |S(u,v) - P| + double c2an = eval[1, 0] * diff; + double c2ad = eval[1, 0].Length * c1v; + + double c2bn = eval[0, 1] * diff; + double c2bd = eval[0, 1].Length * c1v; + + double c2av = c2an / c2ad; + double c2bv = c2bn / c2bd; + + bool c2a = c2av <= tol2; + bool c2b = c2bv <= tol2; + + // If all the criteria are satisfied we are done. + if (c1 && c2a && c2b) + { + return selectedUV; + } + + // Otherwise a new value ( Ui + 1, Vi + 1) is computed using Eq. 6.7 + var ct = NewtonIteration(selectedUV, eval, diff); + + // Ensure the parameters stay in range (Ui+1 E [minU, maxU] and Vi+1 E [minV, maxV]). + if (ct.u < minU) + { + ct = (closedDirectionU) ? (maxU - (minU - ct.u), ct.v) : (minU, ct.v); + } + + if (ct.u > maxU) + { + ct = (closedDirectionU) ? (minU + (ct.u - maxU), ct.v) : (maxU, ct.v); + } + + if (ct.v < minV) + { + ct = (closedDirectionV) ? (ct.u, maxV - (minV - ct.v)) : (ct.u, minV); + } + + if (ct.v > maxV) + { + ct = (closedDirectionV) ? (ct.u, minV + (ct.v - maxV)) : (ct.u, maxV); + } + + // Parameters do not change significantly. + double c3u = (eval[1, 0] * (ct.u - selectedUV.u)).Length; + double c3v = (eval[0, 1] * (ct.v - selectedUV.v)).Length; + + if (c3u + c3v < tol1) + { + return selectedUV; + } + + selectedUV = ct; + t++; + } + + return selectedUV; + } + + /// + /// Newton iteration to minimize the distance between a point and a surface. + /// Corresponds to Eq. 6.5 at page 232 from The NURBS Book by Piegl and Tiller. + /// + /// The parameter uv obtained at the ith Newton iteration. + /// Derivatives of the surface identify as S(u,v). + /// Representing the difference from S(u,v) - P. + /// The minimized parameter. + private static (double u, double v) NewtonIteration((double u, double v) uv, Vector3[,] derivatives, Vector3 r) + { + Vector3 Su = derivatives[1, 0]; + Vector3 Sv = derivatives[0, 1]; + + Vector3 Suu = derivatives[2, 0]; + Vector3 Svv = derivatives[0, 2]; + + Vector3 Suv = derivatives[1, 1]; + Vector3 Svu = derivatives[1, 1]; + + double f = Su * r; + double g = Sv * r; + + // Eq. 6.5 + Vector k = new Vector { -f, -g }; + + Matrix J = new Matrix + { + new List {Su * Su + Suu * r, Su * Sv + Suv * r}, + new List {Su * Sv + Svu * r, Sv * Sv + Svv * r} + }; + + // Eq. 6.6 + Matrix matrixLu = Matrix.Decompose(J, out int[] permutation); + Vector d = Matrix.Solve(matrixLu, permutation, k); + + // Eq. 6.7 + return (d[0] + uv.u, d[1] + uv.v); + } + + /// + /// Defines the U and V parameters for a surface split in both direction, subtracting or adding half of the input parameter based on the quadrant. + /// + private static (double u, double v)[] DefiningUV((double u, double v) surfaceUV, double parameter) + { + double halfParameter = parameter * 0.5; + var UV = new (double u, double v)[4] + { + (surfaceUV.u + halfParameter, surfaceUV.v - halfParameter), + (surfaceUV.u + halfParameter, surfaceUV.v + halfParameter), + (surfaceUV.u - halfParameter, surfaceUV.v - halfParameter), + (surfaceUV.u - halfParameter, surfaceUV.v + halfParameter) + }; + + return UV; + } + + /// + /// Extracts the isoparametric curves (isocurves) at the given parameter and surface direction. + /// + /// The surface object to extract the isocurve. + /// The parameter between 0.0 to 1.0 whether the isocurve will be extracted. + /// The U or V direction whether the isocurve will be extracted. + /// The isocurve extracted. + internal static NurbsCurve Isocurve(NurbsSurface surface, double parameter, SurfaceDirection direction) + { + KnotVector knots = (direction == SurfaceDirection.V) ? surface.KnotsV : surface.KnotsU; + int degree = (direction == SurfaceDirection.V) ? surface.DegreeV : surface.DegreeU; + + Dictionary knotMultiplicity = knots.Multiplicities(); + // If the knotVector already exists in the array, don't make duplicates. + double knotKey = -1; + foreach (KeyValuePair keyValuePair in knotMultiplicity) + { + if (!(Math.Abs(parameter - keyValuePair.Key) < GSharkMath.Epsilon)) continue; + knotKey = keyValuePair.Key; + break; + } + + int knotToInsert = degree + 1; + if (knotKey >= 0) + { + knotToInsert = knotToInsert - knotMultiplicity[knotKey]; + } + + // Insert knots + NurbsSurface refinedSurface = surface; + if (knotToInsert > 0) + { + List knotsToInsert = CollectionHelpers.RepeatData(parameter, knotToInsert); + refinedSurface = KnotVector.Refine(surface, knotsToInsert, direction); + } + + // Obtain the correct index of control points to extract. + int span = knots.Span(degree, parameter); + + if (Math.Abs(parameter - knots[0]) < GSharkMath.Epsilon) + { + span = 0; + } + if (Math.Abs(parameter - knots.Last()) < GSharkMath.Epsilon) + { + span = (direction == SurfaceDirection.V) + ? refinedSurface.ControlPoints[0].Count - 1 + : refinedSurface.ControlPoints.Count - 1; + } + + return direction == SurfaceDirection.V + ? new NurbsCurve(refinedSurface.DegreeU, refinedSurface.KnotsU, CollectionHelpers.Transpose2DArray(refinedSurface.ControlPoints)[span]) + : new NurbsCurve(refinedSurface.DegreeV, refinedSurface.KnotsV, refinedSurface.ControlPoints[span]); + } + } +} \ No newline at end of file diff --git a/src/GShark/Operation/LegendreGaussData.cs b/src/GShark/Core/LegendreGaussData.cs similarity index 99% rename from src/GShark/Operation/LegendreGaussData.cs rename to src/GShark/Core/LegendreGaussData.cs index 52ea19cb..c53d76a4 100644 --- a/src/GShark/Operation/LegendreGaussData.cs +++ b/src/GShark/Core/LegendreGaussData.cs @@ -1,4 +1,4 @@ -namespace GShark.Operation +namespace GShark.Core { /// /// Contains the value of Gaussian Quadrature Weights and Abscissae.
diff --git a/src/GShark/Core/Trigonometry.cs b/src/GShark/Core/Trigonometry.cs index 167bf33a..d4bc0feb 100644 --- a/src/GShark/Core/Trigonometry.cs +++ b/src/GShark/Core/Trigonometry.cs @@ -106,5 +106,45 @@ public static int PointOrder(Point3 pt1, Point3 pt2, Point3 pt3) return (result < 0) ? 1 : 2; } + + /// + /// Finds the closest point on a segment.
+ /// It is not a segment like a line or ray but the part that compose the segment.
+ /// The segment is deconstruct in two points and two t values. + /// + /// Point to project. + /// First point of the segment. + /// Second point of the segment. + /// First t value of the segment. + /// Second t value of the segment. + /// Tuple with the closest point its corresponding tValue on the curve. + public static (double tValue, Point3 pt) ClosestPointToSegment(Point3 point, Point3 segmentPt0, Point3 segmentPt1, double valueT0, double valueT1) + { + Vector3 direction = segmentPt1 - segmentPt0; + double length = direction.Length; + + if (length < GSharkMath.Epsilon) + { + return (tValue: valueT0, pt: segmentPt0); + } + + Vector3 vecUnitized = direction.Unitize(); + Vector3 ptToSegPt0 = point - segmentPt0; + double dotResult = Vector3.DotProduct(ptToSegPt0, vecUnitized); + + if (dotResult < 0.0) + { + return (tValue: valueT0, pt: segmentPt0); + } + + if (dotResult > length) + { + return (tValue: valueT1, pt: segmentPt1); + } + + Point3 pointResult = segmentPt0 + (vecUnitized * dotResult); + double tValueResult = valueT0 + (valueT1 - valueT0) * dotResult / length; + return (tValue: tValueResult, pt: pointResult); + } } } diff --git a/src/GShark/Geometry/NurbsBase.cs b/src/GShark/Geometry/NurbsBase.cs index be20cda3..74628732 100644 --- a/src/GShark/Geometry/NurbsBase.cs +++ b/src/GShark/Geometry/NurbsBase.cs @@ -85,7 +85,7 @@ protected NurbsBase(int degree, KnotVector knots, List controlPoints) /// public KnotVector Knots { get; protected set; } - public double Length => Analyze.CurveLength(this); + public double Length => Analyze.Curve.Length(this); public Point3 StartPoint => PointAt(0.0); @@ -187,7 +187,7 @@ public Point3 PointAt(double t) /// public Point3 PointAtLength(double length) { - double parameter = Analyze.CurveParameterAtLength(this, length); + double parameter = Analyze.Curve.ParameterAtLength(this, length); return Evaluation.CurvePointAt(this, parameter); } @@ -272,7 +272,7 @@ public Vector3 CurvatureAt(double t) } List derivatives = Evaluation.RationalCurveDerivatives(this, t, 2); - return Analyze.Curvature(derivatives[1], derivatives[2]); + return Analyze.Curve.Curvature(derivatives[1], derivatives[2]); } /// @@ -298,7 +298,7 @@ public Plane PerpendicularFrameAt(double t) Vector3 normal = (derivatives[2].Length == 0.0) ? Vector3.PerpendicularTo(derivatives[1]) - : Analyze.Curvature(derivatives[1], derivatives[2]); + : Analyze.Curve.Curvature(derivatives[1], derivatives[2]); Vector3 yDir = Vector3.CrossProduct(derivatives[1], normal); return new Plane(derivatives[0], normal, yDir); @@ -323,7 +323,9 @@ public NurbsBase Reverse() /// public Point3 ClosestPoint(Point3 point) { - return Point4.PointDehomogenizer(Analyze.CurveClosestPoint(this, point, out _)); + double t = Analyze.Curve.ClosestParameter(this, point); + Point3 pointAt = Evaluation.CurvePointAt(this, t); + return Point4.PointDehomogenizer(pointAt); } /// @@ -331,7 +333,7 @@ public Point3 ClosestPoint(Point3 point) /// public double ClosestParameter(Point3 pt) { - return Analyze.CurveClosestParameter(this, pt); + return Analyze.Curve.ClosestParameter(this, pt); } /// @@ -351,7 +353,7 @@ public double ParameterAtLength(double segmentLength) return 1.0; } - return Analyze.CurveParameterAtLength(this, segmentLength); + return Analyze.Curve.ParameterAtLength(this, segmentLength); } /// @@ -369,7 +371,7 @@ public double LengthAt(double t) return Length; } - return Analyze.CurveLength(this, t); + return Analyze.Curve.Length(this, t); } /// @@ -530,7 +532,7 @@ public NurbsBase ReduceDegree(double tolerance = 10e-4) if (equalSegmentLengths) { List curves = Modify.Curve.DecomposeIntoBeziers(this); - List curveLengths = curves.Select(curve => Analyze.BezierCurveLength(curve)).ToList(); + List curveLengths = curves.Select(curve => Analyze.Curve.BezierLength(curve)).ToList(); double totalLength = curveLengths.Sum(); len = totalLength / Math.Ceiling(totalLength / maxSegmentLength); diff --git a/src/GShark/Geometry/NurbsSurface.cs b/src/GShark/Geometry/NurbsSurface.cs index 16a56ec7..dd092d13 100644 --- a/src/GShark/Geometry/NurbsSurface.cs +++ b/src/GShark/Geometry/NurbsSurface.cs @@ -251,7 +251,7 @@ public static NurbsSurface CreateRuledSurface(NurbsBase curveA, NurbsBase curveB /// The closest point on the surface. public Point3 ClosestPoint(Point3 point) { - var (u, v) = Analyze.SurfaceClosestParameter(this, point); + var (u, v) = Analyze.Surface.ClosestParameter(this, point); return new Point3(Evaluation.SurfacePointAt(this, u, v)); } @@ -260,7 +260,10 @@ public Point3 ClosestPoint(Point3 point) /// /// The point to test against. /// The U and V parameters of the surface that are closest to the test point. - public (double U, double V) ClosestParameters(Point3 point) => Analyze.SurfaceClosestParameter(this, point); + public (double U, double V) ClosestParameter(Point3 point) + { + return Analyze.Surface.ClosestParameter(this, point); + } /// /// Evaluate the surface at the given U and V parameters. @@ -298,6 +301,17 @@ public NurbsSurface[] SplitAt(double parameter, SplitDirection direction) return Sampling.Surface.Split(this, parameter, direction); } + /// + /// Extracts the isoparametric curves (isocurves) at the given parameter and surface direction. + /// + /// The parameter between 0.0 to 1.0 whether the isocurve will be extracted. + /// The U or V direction whether the isocurve will be extracted. + /// The isocurve extracted. + public NurbsCurve IsoCurve(double parameter, SurfaceDirection direction) + { + return Analyze.Surface.Isocurve(this, parameter, direction); + } + /// /// Transforms a NURBS surface with the given transformation matrix. /// diff --git a/src/GShark/Operation/Analyze.cs b/src/GShark/Operation/Analyze.cs deleted file mode 100644 index f6955f6e..00000000 --- a/src/GShark/Operation/Analyze.cs +++ /dev/null @@ -1,584 +0,0 @@ -using GShark.Core; -using GShark.ExtendedMethods; -using GShark.Geometry; -using System; -using System.Collections.Generic; -using System.Linq; -using GShark.Enumerations; - -namespace GShark.Operation -{ - /// - /// Contains methods for analyzing NURBS geometry. - /// - public class Analyze - { - /// - /// Computes the approximate length of a rational curve by gaussian quadrature - assumes a smooth curve. - /// - /// The curve object. - /// The parameter at which to approximate the length. - /// - /// The degree of gaussian quadrature to perform. - /// A higher number yields a more exact result, default set to 17. - /// - /// The approximate length. - public static double CurveLength(NurbsBase curve, double u = -1.0, int gaussDegIncrease = 17) - { - double uSet = u < 0.0 ? curve.Knots.Last() : u; - - List crvs = Modify.Curve.DecomposeIntoBeziers(curve); - double sum = 0.0; - - foreach (NurbsBase bezier in crvs) - { - if (!(bezier.Knots[0] + GSharkMath.Epsilon < uSet)) - { - break; - } - double param = Math.Min(bezier.Knots.Last(), uSet); - sum += BezierCurveLength(bezier, param, gaussDegIncrease); - } - - return sum; - } - - /// - /// Computes the curvature vector. - /// The vector has length equal to the radius of the curvature circle and with direction to the center of the circle. - /// https://github.com/mcneel/opennurbs/blob/484ba88836bbedff8fe0b9e574fcd6434b49c21c/opennurbs_math.cpp#L839 - /// - /// The first derivative. - /// The second derivative. - /// The curvature vector. - public static Vector3 Curvature(Vector3 derivative1, Vector3 derivative2) - { - double d1 = derivative1.Length; - - // If the first derivative exists and is zero the curvature is a zero vector. - if (d1 == 0.0) - { - return Vector3.Zero; - } - - Vector3 tangent = derivative1 / d1; - double d2DotTang = (derivative2 * -1) * tangent; - d1 = 1.0 / (d1 * d1); - Vector3 curvature = d1 * (derivative2 + d2DotTang * tangent); // usually identified as k. - - double curvatureLength = curvature.Length; - if (curvatureLength < 1.490116119385E-08) // SqrtEpsilon value that is used when comparing square roots. - { - throw new Exception("Curvature is infinite."); - } - - double radius = (curvature / (curvatureLength * curvatureLength)).Length; - return curvature.Unitize().Amplify(radius); - } - - /// - /// Computes the approximate length of a rational bezier curve by gaussian quadrature - assumes a smooth curve. - /// - /// The curve object. - /// The parameter at which to approximate the length. - /// - /// the degree of gaussian quadrature to perform. - /// A higher number yields a more exact result, default set to 17. - /// - /// The approximate length of a bezier. - public static double BezierCurveLength(NurbsBase curve, double u = -1.0, int gaussDegIncrease = 17) - { - double uSet = u < 0.0 ? curve.Knots.Last() : u; - double z = (uSet - curve.Knots[0]) / 2; - double sum = 0.0; - int gaussDegree = curve.Degree + gaussDegIncrease; - - for (int i = 0; i < gaussDegree; i++) - { - double cu = z * LegendreGaussData.tValues[gaussDegree][i] + z + curve.Knots[0]; - List tan = Evaluation.RationalCurveDerivatives(curve, cu); - - sum += LegendreGaussData.cValues[gaussDegree][i] * tan[1].Length; - } - - return z * sum; - } - - /// - /// Computes the curve parameter at a given length. - /// - /// The curve object. - /// The length to find the parameter. - /// If set less or equal 0.0, the tolerance used is 1e-10. - /// The parameter at the given length. - public static double BezierCurveParamAtLength(NurbsBase curve, double segmentLength, double tolerance) - { - if (segmentLength < 0) - { - return curve.Knots[0]; - } - - // We compute the whole length. - double curveLength = BezierCurveLength(curve); - - if (segmentLength > curveLength) - { - return curve.Knots[curve.Knots.Count - 1]; - } - - // Divide and conquer. - double setTolerance = tolerance <= 0.0 ? GSharkMath.Epsilon : tolerance; - - double startT = curve.Knots[0]; - double startLength = 0.0; - - double endT = curve.Knots[curve.Knots.Count - 1]; - double endLength = curveLength; - - while (endLength - startLength > setTolerance) - { - double midT = (startT + endT) / 2; - double midLength = BezierCurveLength(curve, midT); - - if (midLength > segmentLength) - { - endT = midT; - endLength = midLength; - } - else - { - startT = midT; - startLength = midLength; - } - } - - return (startT + endT) / 2; - } - - /// - /// Computes the closest point on a curve to a given point. - /// - /// The curve object. - /// Point to search from. - /// Parameter of local closest point. - /// The closest point on the curve. - public static Point3 CurveClosestPoint(NurbsBase curve, Point3 point, out double t) - { - t = CurveClosestParameter(curve, point); - return Evaluation.CurvePointAt(curve, t); - } - - /// - /// Computes the closest parameters on a curve to a given point.
- /// Corresponds to page 244 chapter six from The NURBS Book by Piegl and Tiller. - /// - /// The curve object. - /// Point to search from. - /// The closest parameter on the curve. - public static double CurveClosestParameter(NurbsBase curve, Point3 point) - { - double minimumDistance = double.PositiveInfinity; - double tParameter = 0D; - List ctrlPts = curve.ControlPointLocations; - - (List tValues, List pts) = Sampling.Curve.RegularSample(curve, ctrlPts.Count * curve.Degree); - - for (int i = 0; i < pts.Count - 1; i++) - { - double t0 = tValues[i]; - double t1 = tValues[i + 1]; - - Point3 pt0 = pts[i]; - Point3 pt1 = pts[i + 1]; - - var (tValue, pt) = ClosestPointToSegment(point, pt0, pt1, t0, t1); - double distance = pt.DistanceTo(point); - - if (!(distance < minimumDistance)) continue; - minimumDistance = distance; - tParameter = tValue; - } - - int maxIterations = 5; - int j = 0; - // Two zero tolerances can be used to indicate convergence: - double tol1 = GSharkMath.MaxTolerance; // a measure of Euclidean distance; - double tol2 = 0.0005; // a zero cosine measure. - double tVal0 = curve.Knots[0]; - double tVal1 = curve.Knots[curve.Knots.Count - 1]; - bool closedCurve = (ctrlPts[0] - ctrlPts[ctrlPts.Count - 1]).SquareLength < GSharkMath.Epsilon; - double Cu = tParameter; - - // To avoid infinite loop we limited the interaction. - while (j < maxIterations) - { - List e = Evaluation.RationalCurveDerivatives(curve, Cu, 2); - Vector3 diff = e[0] - new Vector3(point); // C(u) - P - - - // First condition, point coincidence: - // |C(u) - p| < e1 - double c1v = diff.Length; - bool c1 = c1v <= tol1; - - // Second condition, zero cosine: - // C'(u) * (C(u) - P) - // ------------------ < e2 - // |C'(u)| |C(u) - P| - double c2n = Vector3.DotProduct(e[1], diff); - double c2d = (e[1] * c1v).Length; - double c2v = c2n / c2d; - bool c2 = Math.Abs(c2v) <= tol2; - - // If at least one of these conditions is not satisfied, - // a new value, ui+l> is computed using the NewtonIterationCurve. - // Then two more conditions are checked. - if (c1 && c2) return Cu; - double ct = NewtonIterationCurve(Cu, e, diff); - - // Ensure that the parameter stays within the boundary of the curve. - if (ct < tVal0) ct = closedCurve ? tVal1 - (ct - tVal0) : tVal0; - if (ct > tVal1) ct = closedCurve ? tVal0 + (ct - tVal1) : tVal1; - - // the parameter does not change significantly, the point is off the end of the curve. - double c3v = (e[1] * (ct - Cu)).Length; - if (c3v < tol1) return Cu; - - Cu = ct; - j++; - } - - return Cu; - } - - /// - /// Computes the closest parameters on a surface to a given point.
- /// Corresponds to page 244 chapter six from The NURBS Book by Piegl and Tiller. - /// - /// The surface object. - /// Point to search from. - /// The closest parameter on the surface. - public static (double u, double v) SurfaceClosestParameter(NurbsSurface surface, Point3 point) - { - double minimumDistance = double.PositiveInfinity; - (double u, double v) selectedUV = (0.5, 0.5); - NurbsSurface splitSrf = surface; - double param = 0.5; - int maxIterations = 5; - - for (int i = 0; i < maxIterations; i++) - { - NurbsSurface[] surfaces = splitSrf.SplitAt(0.5, SplitDirection.Both); - Point3[] pts = surfaces.Select(s => s.PointAt(0.5, 0.5)).ToArray(); - double[] distanceBetweenPts = pts.Select(point.DistanceTo).ToArray(); - if (distanceBetweenPts.All(d => d > minimumDistance)) break; - - (double, double)[] srfUV = DefiningUV(selectedUV, param); - - for (int j = 0; j < distanceBetweenPts.Length; j++) - { - if (!(distanceBetweenPts[j] < minimumDistance)) continue; - minimumDistance = distanceBetweenPts[j]; - selectedUV = srfUV[j]; - splitSrf = surfaces[j]; - } - - param *= 0.5; - } - - int t = 0; - // Two zero tolerances can be used to indicate convergence: - double tol1 = GSharkMath.MaxTolerance; // a measure of Euclidean distance; - double tol2 = 0.0005; // a zero cosine measure. - double minU = surface.KnotsU[0]; - double maxU = surface.KnotsU.Last(); - double minV = surface.KnotsV[0]; - double maxV = surface.KnotsV.Last(); - bool closedDirectionU = surface.IsClosed(SurfaceDirection.U); - bool closedDirectionV = surface.IsClosed(SurfaceDirection.V); - - // To avoid infinite loop we limited the interaction. - while (t < maxIterations) - { - // Get derivatives. - var eval = Evaluation.RationalSurfaceDerivatives(surface, selectedUV.u, selectedUV.v, 2); - - // Convergence criteria: - // First condition, point coincidence: - // |S(u,v) - p| < e1 - Vector3 diff = eval[0, 0] - new Vector3(point); - double c1v = diff.Length; - bool c1 = c1v <= tol1; - - // Second condition, zero cosine: - // |Su(u,v)*(S(u,v) - P)| - // ---------------------- < e2 - // |Su(u,v)| |S(u,v) - P| - // - // |Sv(u,v)*(S(u,v) - P)| - // ---------------------- < e2 - // |Sv(u,v)| |S(u,v) - P| - double c2an = eval[1, 0] * diff; - double c2ad = eval[1, 0].Length * c1v; - - double c2bn = eval[0, 1] * diff; - double c2bd = eval[0, 1].Length * c1v; - - double c2av = c2an / c2ad; - double c2bv = c2bn / c2bd; - - bool c2a = c2av <= tol2; - bool c2b = c2bv <= tol2; - - // If all the criteria are satisfied we are done. - if (c1 && c2a && c2b) - { - return selectedUV; - } - - // Otherwise a new value ( Ui + 1, Vi + 1) is computed using Eq. 6.7 - var ct = NewtonIterationSurface(selectedUV, eval, diff); - - // Ensure the parameters stay in range (Ui+1 E [minU, maxU] and Vi+1 E [minV, maxV]). - if (ct.u < minU) - { - ct = (closedDirectionU) ? (maxU - (minU - ct.u), ct.v) : (minU, ct.v); - } - - if (ct.u > maxU) - { - ct = (closedDirectionU) ? (minU + (ct.u - maxU), ct.v) : (maxU, ct.v); - } - - if (ct.v < minV) - { - ct = (closedDirectionV) ? (ct.u, maxV - (minV - ct.v)) : (ct.u, minV); - } - - if (ct.v > maxV) - { - ct = (closedDirectionV) ? (ct.u, minV + (ct.v - maxV)) : (ct.u, maxV); - } - - // Parameters do not change significantly. - double c3u = (eval[1, 0] * (ct.u - selectedUV.u)).Length; - double c3v = (eval[0, 1] * (ct.v - selectedUV.v)).Length; - - if (c3u + c3v < tol1) - { - return selectedUV; - } - - selectedUV = ct; - t++; - } - - return selectedUV; - } - - /// - /// Newton iteration to minimize the distance between a point and a surface. - /// Corresponds to Eq. 6.5 at page 232 from The NURBS Book by Piegl and Tiller. - /// - /// The parameter uv obtained at the ith Newton iteration. - /// Derivatives of the surface identify as S(u,v). - /// Representing the difference from S(u,v) - P. - /// The minimized parameter. - private static (double u, double v) NewtonIterationSurface((double u, double v) uv, Vector3[,] derivatives, Vector3 r) - { - Vector3 Su = derivatives[1, 0]; - Vector3 Sv = derivatives[0, 1]; - - Vector3 Suu = derivatives[2, 0]; - Vector3 Svv = derivatives[0, 2]; - - Vector3 Suv = derivatives[1, 1]; - Vector3 Svu = derivatives[1, 1]; - - double f = Su * r; - double g = Sv * r; - - // Eq. 6.5 - Vector k = new Vector { -f, -g }; - - Matrix J = new Matrix - { - new List {Su * Su + Suu * r, Su * Sv + Suv * r}, - new List {Su * Sv + Svu * r, Sv * Sv + Svv * r} - }; - - // Eq. 6.6 - Matrix matrixLu = Matrix.Decompose(J, out int[] permutation); - Vector d = Matrix.Solve(matrixLu, permutation, k); - - // Eq. 6.7 - return (d[0] + uv.u, d[1] + uv.v); - } - - /// - /// Defines the U and V parameters for a surface split in both direction, subtracting or adding half of the input parameter based on the quadrant. - /// - private static (double u, double v)[] DefiningUV((double u, double v) surfaceUV, double parameter) - { - double halfParameter = parameter * 0.5; - var UV = new (double u, double v)[4] - { - (surfaceUV.u + halfParameter, surfaceUV.v - halfParameter), - (surfaceUV.u + halfParameter, surfaceUV.v + halfParameter), - (surfaceUV.u - halfParameter, surfaceUV.v - halfParameter), - (surfaceUV.u - halfParameter, surfaceUV.v + halfParameter) - }; - - return UV; - } - - /// - /// Newton iteration to minimize the distance between a point and a curve. - /// - /// The parameter obtained at the ith Newton iteration. - /// Point on curve identify as C'(u) - /// Representing the difference from C(u) - P. - /// The minimized parameter. - private static double NewtonIterationCurve(double u, List derivativePts, Vector3 difference) - { - // The distance from P to C(u) is minimum when f(u) = 0, whether P is on the curve or not. - // C'(u) * ( C(u) - P ) = 0 = f(u) - // C(u) is the curve, p is the point, * is a dot product - double f = Vector3.DotProduct(derivativePts[1], difference); - - // f' = C"(u) * ( C(u) - p ) + C'(u) * C'(u) - double s0 = Vector3.DotProduct(derivativePts[2], difference); - double s1 = Vector3.DotProduct(derivativePts[1], derivativePts[1]); - double df = s0 + s1; - - return u - f / df; - } - - /// - /// Approximates the parameter at a given length on a curve. - /// - /// The curve object. - /// The arc length for which to do the procedure. - /// If set less or equal 0.0, the tolerance used is 1e-10. - /// The parameter on the curve. - public static double CurveParameterAtLength(NurbsBase curve, double segmentLength, double tolerance = -1) - { - if (segmentLength < GSharkMath.Epsilon) return curve.Knots[0]; - if (Math.Abs(curve.Length - segmentLength) < GSharkMath.Epsilon) return curve.Knots[curve.Knots.Count - 1]; - - List curves = Modify.Curve.DecomposeIntoBeziers(curve); - int i = 0; - double curveLength = -GSharkMath.Epsilon; - double segmentLengthLeft = segmentLength; - - // Iterate through the curves consuming the bezier's, summing their length along the way. - while (curveLength < segmentLength && i < curves.Count) - { - double bezierLength = BezierCurveLength(curves[i]); - curveLength += bezierLength; - - if (segmentLength < curveLength + GSharkMath.Epsilon) - return BezierCurveParamAtLength(curves[i], segmentLengthLeft, tolerance); - i++; - segmentLengthLeft -= bezierLength; - } - - return -1; - } - - /// - /// Extracts the isoparametric curves (isocurves) at the given parameter and surface direction. - /// - /// The surface object to extract the isocurve. - /// The parameter between 0.0 to 1.0 whether the isocurve will be extracted. - /// The U or V direction whether the isocurve will be extracted. - /// The isocurve extracted. - public static NurbsBase Isocurve(NurbsSurface surface, double parameter, SurfaceDirection direction) - { - KnotVector knots = (direction == SurfaceDirection.V) ? surface.KnotsV : surface.KnotsU; - int degree = (direction == SurfaceDirection.V) ? surface.DegreeV : surface.DegreeU; - - Dictionary knotMultiplicity = knots.Multiplicities(); - // If the knotVector already exists in the array, don't make duplicates. - double knotKey = -1; - foreach (KeyValuePair keyValuePair in knotMultiplicity) - { - if (!(Math.Abs(parameter - keyValuePair.Key) < GSharkMath.Epsilon)) continue; - knotKey = keyValuePair.Key; - break; - } - - int knotToInsert = degree + 1; - if (knotKey >= 0) - { - knotToInsert = knotToInsert - knotMultiplicity[knotKey]; - } - - // Insert knots - NurbsSurface refinedSurface = surface; - if (knotToInsert > 0) - { - List knotsToInsert = CollectionHelpers.RepeatData(parameter, knotToInsert); - refinedSurface = KnotVector.Refine(surface, knotsToInsert, direction); - } - - // Obtain the correct index of control points to extract. - int span = knots.Span(degree, parameter); - - if (Math.Abs(parameter - knots[0]) < GSharkMath.Epsilon) - { - span = 0; - } - if (Math.Abs(parameter - knots.Last()) < GSharkMath.Epsilon) - { - span = (direction == SurfaceDirection.V) - ? refinedSurface.ControlPoints[0].Count - 1 - : refinedSurface.ControlPoints.Count - 1; - } - - return direction == SurfaceDirection.V - ? new NurbsCurve(refinedSurface.DegreeU, refinedSurface.KnotsU, CollectionHelpers.Transpose2DArray(refinedSurface.ControlPoints)[span]) - : new NurbsCurve(refinedSurface.DegreeV, refinedSurface.KnotsV, refinedSurface.ControlPoints[span]); - } - - /// - /// Finds the closest point on a segment.
- /// The segment is deconstruct in two points and two t values. - /// - /// Point to project. - /// First point of the segment. - /// Second point of the segment. - /// First t value of the segment. - /// Second t value of the segment. - /// Tuple with the closest point its corresponding tValue on the curve. - public static (double tValue, Point3 pt) ClosestPointToSegment(Point3 point, Point3 segmentPt0, - Point3 segmentPt1, double valueT0, double valueT1) - { - Vector3 direction = segmentPt1 - segmentPt0; - double length = direction.Length; - - if (length < GSharkMath.Epsilon) - { - return (tValue: valueT0, pt: segmentPt0); - } - - Vector3 vecUnitized = direction.Unitize(); - Vector3 ptToSegPt0 = point - segmentPt0; - double dotResult = Vector3.DotProduct(ptToSegPt0, vecUnitized); - - if (dotResult < 0.0) - { - return (tValue: valueT0, pt: segmentPt0); - } - - if (dotResult > length) - { - return (tValue: valueT1, pt: segmentPt1); - } - - Point3 pointResult = segmentPt0 + (vecUnitized * dotResult); - double tValueResult = valueT0 + (valueT1 - valueT0) * dotResult / length; - return (tValue: tValueResult, pt: pointResult); - } - } -} \ No newline at end of file diff --git a/src/GShark/Sampling/Curve.cs b/src/GShark/Sampling/Curve.cs index a9b50168..cf404f2f 100644 --- a/src/GShark/Sampling/Curve.cs +++ b/src/GShark/Sampling/Curve.cs @@ -27,7 +27,7 @@ public static class Curve /// A tuple define the t values where the curve is divided and the lengths between each division. internal static List ByCount(NurbsBase curve, int divisions) { - double approximatedLength = Analyze.CurveLength(curve); + double approximatedLength = Analyze.Curve.Length(curve); double arcLengthSeparation = approximatedLength / divisions; var divisionByLength = ByLength(curve, arcLengthSeparation); var tValues = divisionByLength.tValues; @@ -45,7 +45,7 @@ internal static List ByCount(NurbsBase curve, int divisions) internal static (List tValues, List lengths) ByLength(NurbsBase curve, double length) { List curves = Modify.Curve.DecomposeIntoBeziers(curve); - List curveLengths = curves.Select(NurbsBase => Analyze.BezierCurveLength(NurbsBase)).ToList(); + List curveLengths = curves.Select(NurbsBase => Analyze.Curve.BezierLength(NurbsBase)).ToList(); double totalLength = curveLengths.Sum(); List tValues = new List { curve.Knots[0] }; @@ -64,7 +64,7 @@ internal static (List tValues, List lengths) ByLength(NurbsBase while (segmentLength < sum + GSharkMath.Epsilon) { - double t = Analyze.BezierCurveParamAtLength(curves[i], segmentLength - sum2, GSharkMath.MinTolerance); + double t = Analyze.Curve.BezierParameterAtLength(curves[i], segmentLength - sum2, GSharkMath.MinTolerance); tValues.Add(t); divisionLengths.Add(segmentLength);