Revolved Surface

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

@@ -2,11 +2,11 @@
 using GShark.Core;
 using GShark.Enumerations;
 using GShark.Geometry;
-using GShark.Operation;
 using GShark.Test.XUnit.Data;
 using System;
 using System.Collections.Generic;
 using System.Linq;
+using System.Xml.Serialization;
 using Xunit;
 using Xunit.Abstractions;
 
@@ -301,5 +301,134 @@ public void Returns_A_Ruled_Surface_Between_A_Polyline_And_A_Nurbs_Curve(double
             // Assert
             pointAt.EpsilonEquals(expectedPt, GSharkMath.MinTolerance).Should().BeTrue();
         }
+
+        [Fact]
+        public void It_Returns_A_Revolved_Surface_From_A_Line()
+        {
+            // Arrange
+            Ray axis = new Ray(Point3.Origin, Vector3.ZAxis);
+            Line profile = new Line(new Point3(1, 0, 0), new Point3(0, 0, 1));
+
+            List<List<Point3>> expectedPts0 = new List<List<Point3>>
+            {
+                new List<Point3>
+                    {
+                        new Point3(1, 0, 0),
+                        new Point3(0, 0, 1)
+                    },
+                new List<Point3>
+                {
+                    new Point3(1, 0.41421356237309526, 0),
+                    new Point3(0, 0, 1)
+                },
+                new List<Point3>
+                {
+                    new Point3(0.7071067811865476,0.7071067811865476,0),
+                    new Point3(0,0,1)
+                }
+            };
+            List<List<Point3>> expectedPts1 = new List<List<Point3>>
+            {
+                new List<Point3>
+                {
+                    new Point3(1, 0, 0),
+                    new Point3(0, 0, 1)
+                },
+                new List<Point3>
+                {
+                    new Point3(1, 1, 0),
+                    new Point3(0, 0, 1)
+                },
+                new List<Point3>
+                {
+                    new Point3(0,1,0),
+                    new Point3(0,0,1)
+                },
+                new List<Point3>
+                {
+                    new Point3(-1,1,0),
+                    new Point3(0,0,1)
+                },
+                new List<Point3>
+                {
+                    new Point3(-1,0,0),
+                    new Point3(0,0,1)
+                },
+                new List<Point3>
+                {
+                    new Point3(-1,-1,0),
+                    new Point3(0,0,1)
+                },
+                new List<Point3>
+                {
+                    new Point3(0,-1,0),
+                    new Point3(0,0,1)
+                },
+                new List<Point3>
+                {
+                    new Point3(1,-1,0),
+                    new Point3(0,0,1)
+                },
+                new List<Point3>
+                {
+                    new Point3(1, 0, 0),
+                    new Point3(0, 0, 1)
+                }
+            };
+
+            // Act
+            NurbsSurface revolvedSurface0 = NurbsSurface.CreateRevolvedSurface(profile, axis, Math.PI * 0.25);
+            NurbsSurface revolvedSurface1 = NurbsSurface.CreateRevolvedSurface(profile, axis, 2 * Math.PI);
+
+            // Assert
+            revolvedSurface0.ControlPointLocations
+                .Select((pts, i) => pts.SequenceEqual(expectedPts0[i]))
+                .All(res => res)
+                .Should().BeTrue();
+
+            var t = revolvedSurface1.ControlPointLocations
+                .Select((pts, i) => pts.SequenceEqual(expectedPts1[i]))
+                .All(res => res)
+                .Should().BeTrue();
+        }
+
+        [Fact]
+        public void It_Returns_A_Revolved_Surface_From_An_Arc()
+        {
+            // Arrange
+            Ray axis = new Ray(Point3.Origin, Vector3.ZAxis);
+            Arc profile = new Arc(Plane.PlaneZX, 1, new Interval(0, Math.PI * 0.5));
+
+            List<List<Point3>> expectedPts = new List<List<Point3>>
+            {
+                new List<Point3>
+                    {
+                        new Point3(0, 0, 1),
+                        new Point3(1, 0, 1),
+                        new Point3(1, 0, 0)
+                    },
+                new List<Point3>
+                {
+                    new Point3(0, 0, 1),
+                    new Point3(1, 0.41421356237309526, 1),
+                    new Point3(1, 0.41421356237309526, 0)
+                },
+                new List<Point3>
+                {
+                    new Point3(0,0,1),
+                    new Point3(0.7071067811865476,0.7071067811865476,1),
+                    new Point3(0.7071067811865476,0.7071067811865476,0)
+                }
+            };
+
+            // Act
+            NurbsSurface revolvedSurface = NurbsSurface.CreateRevolvedSurface(profile, axis, Math.PI * 0.25);
+
+            // Assert
+            revolvedSurface.ControlPointLocations
+                .Select((pts, i) => pts.SequenceEqual(expectedPts[i]))
+                .All(res => res)
+                .Should().BeTrue();
+        }
     }
 }

src/GShark.Test.XUnit/VerbTests.cs

@@ -1,5 +1,6 @@
 using verb.core;
 using verb.eval;
+using verb.geom;
 using Xunit;
 using Xunit.Abstractions;
 
@@ -468,5 +469,31 @@ public void ElevateDegree()
                 _testOutput.WriteLine($"{elevatePoly.controlPoints[i]}");
             }
         }
+
+        [Fact]
+        public void RevolvedSurface()
+        {
+            var axis = new Array<double>(new double[] { 0, 0, 1 });
+            var pt = new Array<double>(new double[] { 0, 0, 1 });
+            var xaxis = new Array<double>(new double[] { 1, 0, 0 });
+            var center = new Array<double>(new double[] { 0, 0, 0 });
+
+            var arc = new Arc(center, axis, xaxis, 1, 0.0, Math.PI);
+
+            Array<object> pts = new Array<object>();
+
+            pts.push(new Array<double>(new double[] { 1, 0, 0, 1 }));
+            pts.push(new Array<double>(new double[] { 0, 0, 1, 1 }));
+            Array<double> knots = new Array<double>(new double[] { 0.0, 0.0, 1.0, 1.0 });
+            var profile = new verb.core.NurbsCurveData(1, knots, pts);
+
+            var comps = verb.eval.Make.revolvedSurface(arc._data, center, axis, 0.25 * Math.PI);
+
+            for (int i = 0; i < comps.controlPoints.length; i++)
+            {
+                var h = verb.eval.Eval.dehomogenize1d((Array<object>) comps.controlPoints.__a[i]);
+                _testOutput.WriteLine($"{h}");
+            }
+        }
     }
 }

src/GShark/Geometry/NurbsSurface.cs

@@ -1,5 +1,6 @@
 using GShark.Core;
 using GShark.Enumerations;
+using GShark.ExtendedMethods;
 using GShark.Interfaces;
 using System;
 using System.Collections.Generic;
@@ -234,6 +235,131 @@ public static NurbsSurface CreateRuledSurface(NurbsBase curveA, NurbsBase curveB
                 new List<List<Point4>> { curves[0].ControlPoints, curves[1].ControlPoints });
         }
 
+        /// <summary>
+        /// Creates a surface of revolution through an arbitrary angle, and axis.
+        /// <em>Corresponds the algorithm A8.1 of The NURBS Book by Piegl and Tiller.</em>
+        /// </summary>
+        /// <param name="curveProfile">Profile curve.</param>
+        /// <param name="axis">Revolution axis.</param>
+        /// <param name="rotationAngle">Angle in radiance.</param>
+        /// <returns>The revolution surface.</returns>
+        public static NurbsSurface CreateRevolvedSurface(NurbsBase curveProfile, Ray axis, double rotationAngle)
+        {
+            // if angle is less than 90.
+            int arcCount = 1;
+            KnotVector knotsU = Vector.Zero1d(6).ToKnot();
+
+            if (rotationAngle <= Math.PI && rotationAngle > Math.PI / 2)
+            {
+                arcCount = 2;
+                knotsU[3] = knotsU[4] = 0.5;
+            }
+
+            if (rotationAngle <= 3 * Math.PI / 2 && rotationAngle > Math.PI)
+            {
+                arcCount = 3;
+                knotsU = Vector.Zero1d(6 + 2 * (arcCount - 1)).ToKnot();
+                knotsU[3] = knotsU[4] = (double)1 / 3;
+                knotsU[5] = knotsU[6] = (double)2 / 3;
+            }
+
+            if (rotationAngle <= 4 * Math.PI && rotationAngle > 3 * Math.PI / 2)
+            {
+                arcCount = 4;
+                knotsU = Vector.Zero1d(6 + 2 * (arcCount - 1)).ToKnot();
+                knotsU[3] = knotsU[4] = (double)1 / 4;
+                knotsU[5] = knotsU[6] = (double)1 / 2;
+                knotsU[7] = knotsU[8] = (double)3 / 4;
+            }
+
+            // load start and end knots.
+            int t = 3 + 2 * (arcCount - 1);
+            for (int i = 0; i < 3; i++, t++)
+            {
+                knotsU[i] = 0.0;
+                knotsU[t] = 1.0;
+            }
+
+            // some initialization.
+            double divideAngle = rotationAngle / arcCount;
+            int n = 2 * arcCount;
+            double wm = divideAngle / 2; // is the base angle.
+
+            // initialize the sines and cosines only once.
+            double angle = 0.0;
+            double[] sines = new double[arcCount + 1];
+            double[] cosines = new double[arcCount + 1];
+            for (int i = 1; i <= arcCount; i++)
+            {
+                angle += divideAngle;
+                sines[i] = Math.Sin(angle);
+                cosines[i] = Math.Cos(angle);
+            }
+
+            // loop and compute each u row of control points and weights.
+            List<List<Point4>> controlPts = new List<List<Point4>>();
+            for (int r = 0; r < 2 * arcCount + 1; r++)
+            {
+                List<Point4> temp = CollectionHelpers.RepeatData(Point4.Zero, curveProfile.ControlPoints.Count);
+                controlPts.Add(temp);
+            }
+
+            for (int j = 0; j < curveProfile.ControlPointLocations.Count; j++)
+            {
+                Point3 ptO = axis.ClosestPoint(curveProfile.ControlPointLocations[j]);
+                Vector3 vectorX = curveProfile.ControlPointLocations[j] - ptO;
+                double radius = vectorX.Length; // the radius at that length.
+                Vector3 vectorY = Vector3.CrossProduct(axis.Direction, vectorX);
+
+                if (radius > GSharkMath.Epsilon)
+                {
+                    vectorX *= (1 / radius);
+                    vectorY *= (1 / radius);
+                }
+
+                // initialize the first control points and weights.
+                Point3 pt0 = curveProfile.ControlPointLocations[j];
+                controlPts[0][j] = new Point4(pt0, curveProfile.Weights[j]);
+
+                Vector3 tangent0 = vectorY;
+                int index = 0;
+
+                for (int i = 1; i <= arcCount; i++)
+                {
+                    // rotated generatrix point.
+                    Point3 pt2 = (Math.Abs(radius) < GSharkMath.Epsilon)
+                        ? ptO
+                        : ptO + (vectorX * (cosines[i] * radius) + vectorY * (sines[i] * radius));
+
+                    controlPts[index + 2][j] = new Point4(pt2, curveProfile.Weights[j]);
+
+                    // construct the vector tangent to the rotation.
+                    Vector3 rotationTangent = vectorX * (-1 * sines[i]) + vectorY * cosines[i];
+
+                    // construct the next control point.
+                    if (Math.Abs(radius) < GSharkMath.Epsilon)
+                    {
+                        controlPts[index + 1][j] = ptO;
+                    }
+                    else
+                    {
+                        Line ln0 = new Line(pt0, tangent0, tangent0.Length);
+                        Line ln1 = new Line(pt2, rotationTangent, rotationTangent.Length);
+                        Intersection.Intersect.LineLine(ln0, ln1, out Point3 intersectionPt, out _, out _, out _);
+                        controlPts[index + 1][j] = new Point4(intersectionPt, wm * curveProfile.Weights[j]);
+                    }
+
+                    index += 2;
+                    if (i >= arcCount) continue;
+                    pt0 = pt2;
+                    tangent0 = rotationTangent;
+                }
+
+            }
+
+            return new NurbsSurface(2, curveProfile.Degree, knotsU, curveProfile.Knots, controlPts.Select(pts => pts.ToList()).ToList());
+        }
+
         /// <summary>
         /// Evaluates a point at a given U and V parameters.
         /// </summary>

src/GShark/Geometry/Point3.cs

@@ -373,16 +373,15 @@ public double this[int i]
         /// <returns>true if obj is a Point3 and has the same coordinates as this; otherwise false.</returns>
         public override bool Equals(object obj)
         {
-            return (obj is Point3 && this == (Point3)obj);
+            return obj is Point3 point3 && this == point3;
         }
 
         /// <summary>
         /// Check that all values in other are within epsilon of the values in this
         /// </summary>
         /// <param name="other"></param>
         /// <param name="epsilon"></param>
-        /// <returns></returns>
-        //ToDo Using EpsilonEquals everywhere. Perhaps should be moved into actual Equals method of classes.
+        /// <returns>True if the two points have the same coordinates as this; otherwise false.</returns>
         public bool EpsilonEquals(Point3 other, double epsilon)
         {
             return Math.Abs(X - other.X) <= epsilon &&