Dev/mibi/join curve

src/GShark.Test.XUnit/Operation/ModifyTests.cs

@@ -1,10 +1,12 @@
-using FluentAssertions;
+using System;
+using FluentAssertions;
 using GShark.Core;
 using GShark.Geometry;
 using GShark.Geometry.Interfaces;
 using GShark.Operation;
 using GShark.Test.XUnit.Data;
 using System.Collections.Generic;
+using Newtonsoft.Json.Bson;
 using Xunit;
 using Xunit.Abstractions;
 
@@ -222,5 +224,101 @@ public void It_Returns_A_Curve_Where_Degree_Is_Reduced_From_5_To_4()
             ptOnCurve0.DistanceTo(ptOnReducedDegreeCurve0).Should().BeLessThan(GeoSharkMath.MinTolerance);
             ptOnCurve1.DistanceTo(ptOnReducedDegreeCurve1).Should().BeLessThan(tolerance);
         }
+
+        [Fact]
+        public void JoinCurve_Throw_An_Exception_If_The_Number_Of_Curves_Is_Insufficient()
+        {
+            // Arrange
+            ICurve[] curves = {NurbsCurveCollection.NurbsCurvePlanarExample()};
+
+            // Act
+            Func<object> func = () => Modify.JoinCurve(curves);
+
+            // Assert
+            func.Should().Throw<Exception>();
+        }
+
+        [Fact]
+        public void JoinCurve_Throw_An_Exception_If_Curves_Are_Close_Enough_To_Be_Joined()
+        {
+            // Arrange
+            ICurve[] curves = { NurbsCurveCollection.NurbsCurvePlanarExample(), NurbsCurveCollection.NurbsCurveQuadratic3DBezier() };
+
+            // Act
+            Func<object> func = () => Modify.JoinCurve(curves);
+
+            // Assert
+            func.Should().Throw<Exception>();
+        }
+
+        [Fact]
+        public void Returns_A_Curve_Joining_Different_Types_Of_Curves()
+        {
+            // Arrange
+            int degree = 3;
+            List<Point3> pts = new List<Point3>
+            {
+                new Point3(0, 5, 5),
+                new Point3(0, 0, 0),
+                new Point3(5, 0, 0),
+                new Point3(5, 0, 5),
+                new Point3(5, 5, 5),
+                new Point3(5, 5, 0)
+            };
+
+            NurbsCurve curve = new NurbsCurve(pts, degree);
+            Line ln = new Line(new Point3(5, 5, 0), new Point3(5, 5, -2.5));
+            Arc arc = Arc.ByStartEndDirection(new Point3(5, 5, -2.5), new Point3(10, 5, -7.5), new Vector3(0,0,-1));
+            Point3 expectedPt1 = new Point3(3.34125, 0.005, 0.6125 );
+            Point3 expectedPt2 = new Point3(5, 2.363281, 4.570313 );
+            Point3 expectedPt3 = new Point3(5.351058, 5, -4.340474);
+            ICurve[] curves = {curve, ln, arc};
+
+            // Act
+            ICurve joinedCurve = Modify.JoinCurve(curves);
+            Point3 pt1 = joinedCurve.PointAt(0.1);
+            Point3 pt2 = joinedCurve.PointAt(0.25);
+            Point3 pt3 = joinedCurve.PointAt(0.75);
+
+            // Arrange
+            pt1.DistanceTo(expectedPt1).Should().BeLessThan(GeoSharkMath.MinTolerance);
+            pt2.DistanceTo(expectedPt2).Should().BeLessThan(GeoSharkMath.MinTolerance);
+            pt3.DistanceTo(expectedPt3).Should().BeLessThan(GeoSharkMath.MinTolerance);
+        }
+
+        [Fact]
+        public void Returns_A_Curve_Joining_Polylines_And_Lines()
+        {
+            // Arrange
+            List<Point3> pts = new List<Point3>
+            {
+                new Point3(0, 5, 5),
+                new Point3(0, 0, 0),
+                new Point3(5, 0, 0),
+                new Point3(5, 0, 5),
+                new Point3(5, 5, 5),
+                new Point3(5, 5, 0)
+            };
+
+            Polyline poly = new Polyline(pts);
+            Line ln = new Line(new Point3(5, 5, 0), new Point3(5, 5, -2.5));
+            Point3 expectedPt1 = new Point3(0, 2.0, 2.0);
+            Point3 expectedPt2 = new Point3(2.5, 0, 0);
+            Point3 expectedPt3 = new Point3(5, 5.0, 4.0);
+            ICurve[] curves = { poly, ln };
+
+            // Act
+            ICurve joinedCurve = Modify.JoinCurve(curves);
+            Point3 pt1 = joinedCurve.PointAt(0.1);
+            Point3 pt2 = joinedCurve.PointAt(0.25);
+            Point3 pt3 = joinedCurve.PointAt(0.70);
+
+            // Arrange
+            joinedCurve.Degree.Should().Be(1);
+            pt1.DistanceTo(expectedPt1).Should().BeLessThan(GeoSharkMath.MinTolerance);
+            pt2.DistanceTo(expectedPt2).Should().BeLessThan(GeoSharkMath.MinTolerance);
+            pt3.DistanceTo(expectedPt3).Should().BeLessThan(GeoSharkMath.MinTolerance);
+        }
     }
 }
+

src/GShark/Core/LinearAlgebra.cs

@@ -118,22 +118,24 @@ public static Vector3 GetRotationAxis(Transform transform)
         /// <summary>
         /// Gets the orientation between tree points in the plane.<br/>
         /// The order can be: collinear (result 0), clockwise (result 1), counterclockwise (result 2)<br/>
-        /// https://www.geeksforgeeks.org/orientation-3-ordered-points/
+        /// https://math.stackexchange.com/questions/2386810/orientation-of-three-points-in-3d-space
         /// </summary>
         /// <param name="pt1">First point.</param>
         /// <param name="pt2">Second point.</param>
         /// <param name="pt3">Third point.</param>
         /// <returns>The result expressed as a value between 0 and 2.</returns>
         public static int Orientation(Point3 pt1, Point3 pt2, Point3 pt3)
         {
-            double result = (pt2[1] - pt1[1]) * (pt3[0] - pt2[0]) - (pt2[0] - pt1[0]) * (pt3[1] - pt2[1]);
+            Plane pl = new Plane(pt1, pt2, pt3);
+            Vector3 n = Vector3.CrossProduct(pt2 - pt1, pt3 - pt1);
+            double result = Vector3.DotProduct(pl.ZAxis, n.Unitize());
 
             if (Math.Abs(result) < GeoSharkMath.Epsilon)
             {
                 return 0;
             }
 
-            return (result > 0) ? 1 : 2;
+            return (result < 0) ? 1 : 2;
         }
 
         /// <summary>

src/GShark/Geometry/Arc.cs

@@ -57,6 +57,19 @@ public Arc(Plane plane, double radius, double angleRadians)
         /// <param name="pt3">End point of the arc.</param>
         public Arc(Point3 pt1, Point3 pt2, Point3 pt3)
         {
+            if (!pt1.IsValid)
+            {
+                throw new Exception("The first point is not valid.");
+            }
+            if (!pt2.IsValid)
+            {
+                throw new Exception("The second point is not valid.");
+            }
+            if (!pt3.IsValid)
+            {
+                throw new Exception("The third point is not valid.");
+            }
+
             Point3 center = Trigonometry.PointAtEqualDistanceFromThreePoints(pt1, pt2, pt3);
             Vector3 normal = Vector3.ZAxis.PerpendicularTo(pt1, pt2, pt3);
             Vector3 xDir = pt1 - center;
@@ -180,11 +193,29 @@ public BoundingBox BoundingBox
         /// <returns>An arc.</returns>
         public static Arc ByStartEndDirection(Point3 ptStart, Point3 ptEnd, Vector3 dir)
         {
+            if (!ptStart.IsValid)
+            {
+                throw new Exception("The first point is not valid.");
+            }
+            if (!ptEnd.IsValid)
+            {
+                throw new Exception("The second point is not valid.");
+            }
+            if (!dir.IsValid)
+            {
+                throw new Exception("The tangent is not valid.");
+            }
+
             Vector3 vec0 = dir.Unitize();
             Vector3 vec1 = (ptEnd - ptStart).Unitize();
-            if (vec1.Length == 0.0)
+            if (vec1.Length.Equals(0.0))
+            {
+                throw new Exception("Start and End point of the arc are coincident. Enable to create an arc");
+            }
+
+            if (vec0.IsParallelTo(vec1) != 0)
             {
-                throw new Exception("Points must not be coincident.");
+                throw new Exception("Tangent is parallel with the endpoints. Enable to create an arc");
             }
 
             Vector3 vec2 = (vec0 + vec1).Unitize();
@@ -293,19 +324,19 @@ private void ToNurbsCurve()
 
             // Number of arcs.
             double piNum = 0.5 * Math.PI;
-            if (Angle <= piNum)
+            if ((Angle - piNum) <= GeoSharkMath.Epsilon)
             {
                 numberOfArc = 1;
                 ctrPts = new Point3[3];
                 weights = new double[3];
             }
-            else if (Angle <= piNum * 2)
+            else if ((Angle - piNum * 2) <= GeoSharkMath.Epsilon)
             {
                 numberOfArc = 2;
                 ctrPts = new Point3[5];
                 weights = new double[5];
             }
-            else if (Angle <= piNum * 3)
+            else if ((Angle - piNum * 3) <= GeoSharkMath.Epsilon)
             {
                 numberOfArc = 3;
                 ctrPts = new Point3[7];

src/GShark/Operation/Modify.cs

@@ -691,5 +691,62 @@ public static ICurve ReduceDegree(ICurve curve, double tolerance = 10e-4)
 
             return new NurbsCurve(p - 1, Uh, Pw);
         }
+
+        /// <summary>
+        /// Joins all the curves, if it is possible.
+        /// </summary>
+        /// <param name="curves">Curves to join.</param>
+        /// <returns>A curve that is the result of joining all the curves.</returns>
+        public static ICurve JoinCurves(IList<ICurve> curves)
+        {
+            if (curves == null)
+            {
+                throw new Exception("The set of curves is empty.");
+            }
+
+            if (curves.Count <= 1)
+            {
+                throw new Exception("Insufficient curves for join operation.");
+            }
+
+            for (int i = 0; i < curves.Count - 1; i++)
+            {
+                if (curves[i].ControlPoints.Last().DistanceTo(curves[i + 1].ControlPoints[0]) > GeoSharkMath.MinTolerance)
+                {
+                    throw new Exception($"Curve at {i} and curve at {i + 1} don't touch each other.");
+                }
+            }
+
+            // Extract the biggest degree between the curves.
+            int finalDegree = curves.Max(c => c.Degree);
+
+            // Homogenized degree curves.
+            IEnumerable<ICurve> homogenizedCurves = curves.Select(curve => curve.Degree != finalDegree ? ElevateDegree(curve, finalDegree) : curve);
+
+            // Join curves.
+            List<double> joinedKnots = new List<double>();
+            List<Point4> joinedControlPts = new List<Point4>();
+            double endDomain = 0;
+
+            foreach (ICurve curve in homogenizedCurves)
+            {
+                if (joinedKnots.Count == 0)
+                {
+                    joinedKnots.AddRange(curve.Knots.Take(curve.Knots.Count - (finalDegree + 1)));
+                    joinedControlPts.AddRange(curve.ControlPoints);
+                }
+                else
+                {
+                    joinedKnots.AddRange(curve.Knots.Take(curve.Knots.Count - (finalDegree + 1)).Skip(1).Select(k => k + endDomain));
+                    joinedControlPts.AddRange(curve.ControlPoints.Skip(1));
+                }
+
+                endDomain += curve.Knots.Last();
+            }
+
+            // Appending the last knot to the end.
+            joinedKnots.AddRange(Sets.RepeatData(endDomain, finalDegree + 1));
+            return new NurbsCurve(finalDegree, joinedKnots.ToKnot().Normalize(), joinedControlPts);
+        }
     }
 }