Dev/mibi/surface closest parameter

src/GShark.Test.XUnit/Geometry/NurbsSurfaceTests.cs

@@ -223,5 +223,24 @@ public void Split_Surface_Throws_An_Exception_If_Parameter_Is_Outside_The_Domain
             // Assert
             func.Should().Throw<ArgumentOutOfRangeException>();
         }
+
+        [Theory]
+        [InlineData(new double[] { 2.60009, 7.69754, 3.408162 }, new double[] { 2.5, 7, 5 })]
+        [InlineData(new double[] { 2.511373, 1.994265, 0.887211 }, new double[] { 2.5, 1.5, 2 })]
+        [InlineData(new double[] { 8.952827, 2.572942, 0.735217 }, new double[] { 9, 2.5, 1 })]
+        [InlineData(new double[] { 5.073733, 4.577509, 1.978153 }, new double[] { 5, 5, 1 })]
+        public void Returns_The_Closest_Point_On_The_Surface(double[] expectedPt, double[] testPt)
+        {
+            // Arrange
+            NurbsSurface surface = NurbsSurfaceCollection.SurfaceFromPoints();
+            Point3 pt = new Point3(testPt[0], testPt[1], testPt[2]);
+            Point3 expectedClosestPt = new Point3(expectedPt[0], expectedPt[1], expectedPt[2]);
+
+            // Act
+            Point3 closestPt = surface.ClosestPoint(pt);
+
+            // Assert
+            closestPt.DistanceTo(expectedClosestPt).Should().BeLessThan(GeoSharkMath.MaxTolerance);
+        }
     }
 }

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

@@ -124,5 +124,25 @@ public void RationalCurveParameterAtLength_Returns_Parameter_t_At_The_Given_Leng
             // 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(GeoSharkMath.MaxTolerance);
+            (closestParameter.v - expectedUV.v).Should().BeLessThan(GeoSharkMath.MaxTolerance);
+        }
     }
 }

src/GShark.Test.XUnit/VerbTests.cs

@@ -365,5 +365,95 @@ public void CurveKnotRefine2()
 
             _testOutput.WriteLine($"{after.controlPoints}");
         }
+
+        [Fact]
+        public void SurfaceClosestPt()
+        {
+            Array<double> pt1 = new Array<double>(new double[] { 0.0, 0.0, 0.0 });
+            Array<double> pt2 = new Array<double>(new double[] { 10.0, 0.0, 0.0 });
+            Array<double> pt3 = new Array<double>(new double[] { 10.0, 10.0, 2.0 });
+            Array<double> pt4 = new Array<double>(new double[] { 0.0, 10.0, 4.0 });
+            Array<double> pt5 = new Array<double>(new double[] { 5.0, 0.0, 0.0 });
+            Array<double> pt6 = new Array<double>(new double[] { 5.0, 10.0, 5.0 });
+
+            Array<object> pts = new Array<object>();
+            Array<object> c1 = new Array<object>();
+            c1.push(pt1);
+            c1.push(pt4);
+
+            Array<object> c2 = new Array<object>();
+            c2.push(pt5);
+            c2.push(pt6);
+
+            Array<object> c3 = new Array<object>();
+            c3.push(pt2);
+            c3.push(pt3);
+
+            pts.push(c1);
+            pts.push(c2);
+            pts.push(c3);
+
+            Array<object> weight = new Array<object>();
+            Array<object> w1 = new Array<object>();
+            w1.push(1);
+            w1.push(1);
+
+            Array<object> w2 = new Array<object>();
+            w2.push(1);
+            w2.push(2);
+
+            Array<object> w3 = new Array<object>();
+            w3.push(1);
+            w3.push(1);
+
+            weight.push(w1);
+            weight.push(w2);
+            weight.push(w3);
+
+            Array<double> knotsU = new Array<double>(new double[] { 0, 0, 0, 1, 1, 1 });
+            Array<double> knotsV = new Array<double>(new double[] { 0, 0, 1, 1 });
+
+            var surface = verb.geom.NurbsSurface.byKnotsControlPointsWeights(2, 1, knotsU, knotsV, pts, weight);
+            Array<object> derivatives = verb.eval.Eval.rationalSurfaceDerivatives(surface._data, 0.25, 0.25, 2);
+
+            var pt = new Array<double>(new double[] {2.5, 1.5, 2});
+
+            var rational =
+                verb.eval.Analyze.rationalSurfaceClosestParam(surface._data, pt);
+
+            var ptd = (Array<object>) derivatives.__a[0];
+            var diff = verb.core.Vec.sub((Array<double>)ptd.__a[0], pt);
+
+            var Su = (Array<object>)derivatives.__a[1];
+            var Sv = (Array<object>)derivatives.__a[0];
+
+            var Suu = (Array<object>)derivatives.__a[2];
+            var Svv = (Array<object>)derivatives.__a[0];
+
+            var Suv = (Array<object>)derivatives.__a[1];
+            var Svu = (Array<object>)derivatives.__a[1];
+
+            var f = verb.core.Vec.dot((Array<double>) Su.__a[0], diff);
+            var g = verb.core.Vec.dot((Array<double>)Sv.__a[1], diff);
+
+            var k = new Array<double>(new double[] { -f, -g });
+            var uv = new Array<double>(new double[] { 0.25, 0.25 });
+
+            var J00 = Vec.dot((Array<double>)Su.__a[0], (Array<double>)Su.__a[0]) + Vec.dot((Array<double>)Suu.__a[0], diff);
+            var J01 = Vec.dot((Array<double>)Su.__a[0], (Array<double>)Sv.__a[1]) + Vec.dot((Array<double>)Suv.__a[1], diff);
+            var J10 = Vec.dot((Array<double>)Su.__a[0], (Array<double>)Sv.__a[1]) + Vec.dot((Array<double>)Svu.__a[1], diff);
+            var J11 = Vec.dot((Array<double>)Sv.__a[1], (Array<double>)Sv.__a[1]) + Vec.dot((Array<double>)Svv.__a[2], diff);
+
+            Array<object> matrix = new Array<object>();
+            matrix.push(new Array<object>(new object[] { J00, J01 }));
+            matrix.push(new Array<object>(new object[] { J10, J11 }));
+
+            verb.core._Mat.LUDecomp matrixLU = Mat.LU(matrix);
+
+            var d = verb.core.Mat.solve(matrix, k);
+            var res = Vec.add(d, uv);
+
+            //_testOutput.WriteLine($"{{{rational}}}");
+        }
     }
 }

src/GShark/Geometry/NurbsSurface.cs

@@ -88,15 +88,32 @@ internal NurbsSurface(int degreeU, int degreeV, KnotVector knotsU, KnotVector kn
         public List<List<double>> Weights { get; }
 
         /// <summary>
-        /// A 2D collection of points, the vertical U direction increases from bottom to top, the V direction from left to right.
+        /// A 2D collection of points, U direction increases from left to right, the V direction from bottom to top.
         /// </summary>
         public List<List<Point3>> LocationPoints { get; }
 
         /// <summary>
-        /// A 2d collection of control points, the vertical U direction increases from bottom to top, the V direction from left to right.
+        /// A 2d collection of control points, U direction increases from left to right, the V direction from bottom to top.
         /// </summary>
         internal List<List<Point4>> ControlPoints { get; }
 
+        /// <summary>
+        /// Checks if a NURBS surface is closed.<br/>
+        /// A surface is closed if the first points and the lasts in a direction are coincident.
+        /// </summary>
+        /// <returns>True if the curve is closed.</returns>
+        // ToDo: Add test using lofts.
+        public bool IsClosed(SurfaceDirection direction)
+        {
+            KnotVector knot = (direction == SurfaceDirection.U) ? KnotsU : KnotsV;
+            int degree = (direction == SurfaceDirection.U) ? DegreeU : DegreeV;
+
+            if (!knot.IsClamped(degree)) return knot.IsKnotVectorPeriodic(degree);
+
+            var pts2d = (direction == SurfaceDirection.U) ? LocationPoints : Sets.Reverse2DMatrixData(LocationPoints);
+            return pts2d.All(pts => pts[0].DistanceTo(pts.Last()) < GeoSharkMath.Epsilon);
+        }
+
         /// <summary>
         /// Constructs a NURBS surface from four corners are expected in counter-clockwise order.<br/>
         /// The surface is defined with degree 1.
@@ -144,6 +161,24 @@ public static NurbsSurface CreateFromPoints(int degreeU, int degreeV, List<List<
         /// <returns>A evaluated point.</returns>
         public Point3 PointAt(double u, double v) => new Point3(Evaluation.SurfacePointAt(this, u, v));
 
+        /// <summary>
+        /// Computes the point on the surface that is closest to the test point.
+        /// </summary>
+        /// <param name="point">The point to test against.</param>
+        /// <returns>The closest point on the surface.</returns>
+        public Point3 ClosestPoint(Point3 point)
+        {
+            var param = Analyze.SurfaceClosestParameter(this, point);
+            return new Point3(Evaluation.SurfacePointAt(this, param.u, param.v));
+        }
+
+        /// <summary>
+        /// Computes the U and V parameters of the surface that is closest to the test point.
+        /// </summary>
+        /// <param name="point">The point to test against.</param>
+        /// <returns>The U and V parameters of the surface that are closest to the test point.</returns>
+        public (double U, double V) ClosestParameters(Point3 point) => Analyze.SurfaceClosestParameter(this, point);
+
         /// <summary>
         /// Evaluate the surface at the given U and V parameters.
         /// </summary>

src/GShark/Geometry/Vector3.cs

@@ -692,8 +692,8 @@ public int IsParallelTo(Vector3 other, double angleTolerance = GeoSharkMath.Angl
         /// <summary>
         /// Computes the dot product between two vectors.
         /// </summary>
-        /// <param name="a">The first vector.</param>
-        /// <param name="b">The second vector with which compute the dot product.</param>
+        /// <param name="vector1">The first vector.</param>
+        /// <param name="vector2">The second vector with which compute the dot product.</param>
         /// <returns>The dot product.</returns>
         public static double DotProduct(Vector3 vector1, Vector3 vector2)
         {