Handle parallel-vector cross-product cases in MathUtil plane-angle ca…

…lcs. remove Frame3f asserts that were wrong. Support closed loops in TubeGenerator. improvements to MeshIsoCurves.

curve/DGraph3Util.cs

@@ -90,8 +90,9 @@ public static Curves ExtractCurves(DGraph3 graph)
                         path.Closed = true;
                         c.Loops.Add(path);
                         // need to mark incoming edge as used...but is it valid now?
-                        Util.gDevAssert(eid != int.MaxValue);
-                        used.Add(eid);
+                        //Util.gDevAssert(eid != int.MaxValue);
+                        if ( eid != int.MaxValue )
+                            used.Add(eid);
 
                     } else {
                         c.Paths.Add(path);

math/Frame3f.cs

@@ -138,12 +138,10 @@ public Frame3f Scaled(Vector3f scale)
 
         public void Rotate(Quaternionf q)
         {
-            Debug.Assert(rotation.w != 0);      // catch un-initialized quaternions
             rotation = q * rotation;
         }
         public Frame3f Rotated(Quaternionf q)
         {
-            Debug.Assert(rotation.w != 0);
             return new Frame3f(this.origin, q * this.rotation);
         }
         public Frame3f Rotated(float fAngle, int nAxis)
@@ -164,27 +162,23 @@ public void RotateAroundAxes(Quaternionf q)
 
         public void RotateAround(Vector3f point, Quaternionf q)
         {
-            Debug.Assert(rotation.w != 0);
             Vector3f dv = q * (origin - point);
             rotation = q * rotation;
             origin = point + dv;
         }
         public Frame3f RotatedAround(Vector3f point, Quaternionf q)
         {
-            Debug.Assert(rotation.w != 0);
             Vector3f dv = q * (this.origin - point);
             return new Frame3f(point + dv, q * this.rotation);
         }
 
         public void AlignAxis(int nAxis, Vector3f vTo)
         {
-            Debug.Assert(rotation.w != 0);
             Quaternionf rot = Quaternionf.FromTo(GetAxis(nAxis), vTo);
             Rotate(rot);
         }
         public void ConstrainedAlignAxis(int nAxis, Vector3f vTo, Vector3f vAround)
         {
-            Debug.Assert(rotation.w != 0);
             Vector3f axis = GetAxis(nAxis);
             float fAngle = MathUtil.PlaneAngleSignedD(axis, vTo, vAround);
             Quaternionf rot = Quaternionf.AxisAngleD(vAround, fAngle);

math/MathUtil.cs

@@ -246,7 +246,11 @@ public static float PlaneAngleSignedD(Vector3f vFrom, Vector3f vTo, int nPlaneNo
             vFrom[nPlaneNormalIdx] = vTo[nPlaneNormalIdx] = 0.0f;
             vFrom.Normalize();
             vTo.Normalize();
-            float fSign = Math.Sign(vFrom.Cross(vTo)[nPlaneNormalIdx]);
+            Vector3f c = vFrom.Cross(vTo);
+            if (c.LengthSquared < MathUtil.ZeroTolerancef) {        // vectors are parallel
+                return vFrom.Dot(vTo) < 0 ? 180.0f : 0;
+            }
+            float fSign = Math.Sign(c[nPlaneNormalIdx]);
             float fAngle = fSign * Vector3f.AngleD(vFrom, vTo);
             return fAngle;
         }
@@ -255,7 +259,11 @@ public static double PlaneAngleSignedD(Vector3d vFrom, Vector3d vTo, int nPlaneN
             vFrom[nPlaneNormalIdx] = vTo[nPlaneNormalIdx] = 0.0;
             vFrom.Normalize();
             vTo.Normalize();
-            double fSign = Math.Sign(vFrom.Cross(vTo)[nPlaneNormalIdx]);
+            Vector3d c = vFrom.Cross(vTo);
+            if (c.LengthSquared < MathUtil.ZeroTolerance) {        // vectors are parallel
+                return vFrom.Dot(vTo) < 0 ? 180.0 : 0;
+            }
+            double fSign = Math.Sign(c[nPlaneNormalIdx]);
             double fAngle = fSign * Vector3d.AngleD(vFrom, vTo);
             return fAngle;
         }
@@ -267,6 +275,9 @@ public static float PlaneAngleSignedD(Vector3f vFrom, Vector3f vTo, Vector3f pla
             vFrom.Normalize();
             vTo.Normalize();
             Vector3f c = Vector3f.Cross(vFrom, vTo);
+            if (c.LengthSquared < MathUtil.ZeroTolerancef) {        // vectors are parallel
+                return vFrom.Dot(vTo) < 0 ? 180.0f : 0;
+            }
             float fSign = Math.Sign(Vector3f.Dot(c, planeN));
             float fAngle = fSign * Vector3f.AngleD(vFrom, vTo);
             return fAngle;
@@ -278,6 +289,9 @@ public static double PlaneAngleSignedD(Vector3d vFrom, Vector3d vTo, Vector3d pl
             vFrom.Normalize();
             vTo.Normalize();
             Vector3d c = Vector3d.Cross(vFrom, vTo);
+            if (c.LengthSquared < MathUtil.ZeroTolerance) {        // vectors are parallel
+                return vFrom.Dot(vTo) < 0 ? 180.0 : 0;
+            }
             double fSign = Math.Sign(Vector3d.Dot(c, planeN));
             double fAngle = fSign * Vector3d.AngleD(vFrom, vTo);
             return fAngle;

mesh_generators/GenCylGenerators.cs

@@ -10,6 +10,7 @@ public class TubeGenerator : MeshGenerator
         public List<Vector3d> Vertices;
         public Polygon2d Polygon;
         public bool Capped = true;
+        public bool ClosedLoop = false;
 
         // [TODO] Frame3d ??
         public Frame3f Frame = Frame3f.Identity;
@@ -29,15 +30,16 @@ override public MeshGenerator Generate()
             int nRings = Vertices.Count;
             int nRingSize = (NoSharedVertices) ? Slices + 1 : Slices;
             int nCapVertices = (NoSharedVertices) ? Slices + 1 : 1;
-            if (Capped == false)
+            if (Capped == false || ClosedLoop == true)
                 nCapVertices = 0;
 
             vertices = new VectorArray3d(nRings * nRingSize + 2 * nCapVertices);
             uv = new VectorArray2f(vertices.Count);
             normals = new VectorArray3f(vertices.Count);
 
-            int nSpanTris = (Vertices.Count - 1) * (2 * Slices);
-            int nCapTris = (Capped) ? 2 * Slices : 0;
+            int quad_strips = ClosedLoop ? (nRings) : (nRings-1);
+            int nSpanTris = quad_strips * (2 * Slices);
+            int nCapTris = (Capped && ClosedLoop == false) ? 2 * Slices : 0;
             triangles = new IndexArray3i(nSpanTris + nCapTris);
 
             Frame3f fCur = new Frame3f(Frame);
@@ -78,9 +80,12 @@ override public MeshGenerator Generate()
 
             // generate triangles
             int ti = 0;
-            for (int ri = 0; ri < nRings-1; ++ri) {
+            int nStop = (ClosedLoop) ? nRings : (nRings - 1);
+            for (int ri = 0; ri < nStop; ++ri) {
                 int r0 = ri * nRingSize;
                 int r1 = r0 + nRingSize;
+                if (ClosedLoop && ri == nStop - 1)
+                    r1 = 0;
                 for (int k = 0; k < nRingSize - 1; ++k) {
                     triangles.Set(ti++, r0 + k, r0 + k + 1, r1 + k + 1, Clockwise);
                     triangles.Set(ti++, r0 + k, r1 + k + 1, r1 + k, Clockwise);
@@ -91,7 +96,7 @@ override public MeshGenerator Generate()
                 }
             }
 
-            if (Capped) {
+            if (Capped && ClosedLoop == false) {
                 // add endcap verts
                 int nBottomC = nRings * nRingSize;
                 vertices[nBottomC] = fStart.Origin;

mesh_ops/MeshIsoCurves.cs

@@ -10,6 +10,11 @@ public class MeshIsoCurves
         public DMesh3 Mesh;
         public Func<Vector3d, double> ValueF = null;
 
+        /// <summary>
+        /// Optional - provide value at vertex (which may be precomputed)
+        /// </summary>
+        public Func<int, double> VertexValueF = null;
+
         public DGraph3 Graph = null;
 
         public enum TriangleCase
@@ -64,7 +69,11 @@ protected void compute_full(IEnumerable<int> Triangles)
 
                 Vector3dTuple3 tv = new Vector3dTuple3();
                 Mesh.GetTriVertices(tid, ref tv.V0, ref tv.V1, ref tv.V2);
-                Vector3d f = new Vector3d(ValueF(tv.V0), ValueF(tv.V1), ValueF(tv.V2));
+                Index3i triVerts = Mesh.GetTriangle(tid);
+
+                Vector3d f = (VertexValueF != null) ?
+                    new Vector3d(VertexValueF(triVerts.a), VertexValueF(triVerts.b), VertexValueF(triVerts.c))
+                    : new Vector3d(ValueF(tv.V0), ValueF(tv.V1), ValueF(tv.V2));
 
                 // round f to 0 within epsilon?
 
@@ -73,7 +82,6 @@ protected void compute_full(IEnumerable<int> Triangles)
                 if (f.x > 0 && f.y > 0 && f.z > 0)
                     continue;
 
-                Index3i triVerts = Mesh.GetTriangle(tid);
                 Index3i triEdges = Mesh.GetTriEdges(tid);
 
                 if (f.x * f.y * f.z == 0) {
@@ -125,10 +133,15 @@ protected void compute_full(IEnumerable<int> Triangles)
                     int e1 = (cross_verts.c == int.MinValue) ? 1 : 2;
                     int ev0 = cross_verts[e0];
                     int ev1 = cross_verts[e1];
-                    Util.gDevAssert(ev0 != int.MinValue && ev1 != int.MinValue);
-                    int graph_eid = Graph.AppendEdge(ev0, ev1, (int)TriangleCase.EdgeEdge);
-                    if (WantGraphEdgeInfo)
-                        add_edge_edge(graph_eid, tid, new Index2i(triEdges[e0], triEdges[e1]));
+                    // [RMS] if function is garbage, we can end up w/ case where both crossings
+                    //   happen at same vertex, even though values are not the same (eg if
+                    //   some values are double.MaxValue). We will just fail in these cases.
+                    if (ev0 != ev1) {
+                        Util.gDevAssert(ev0 != int.MinValue && ev1 != int.MinValue);
+                        int graph_eid = Graph.AppendEdge(ev0, ev1, (int)TriangleCase.EdgeEdge);
+                        if (WantGraphEdgeInfo)
+                            add_edge_edge(graph_eid, tid, new Index2i(triEdges[e0], triEdges[e1]));
+                    }
                 }
             }
 
@@ -184,9 +197,11 @@ Vector3d find_crossing(Vector3d a, Vector3d b, double fA, double fB)
             double t = 0.5;
             if (fA < fB) {
                 t = (0 - fA) / (fB - fA);
+                t = MathUtil.Clamp(t, 0, 1);
                 return (1 - t) * a + (t) * b;
             } else if ( fB < fA ) {
                 t = (0 - fB) / (fA - fB);
+                t = MathUtil.Clamp(t, 0, 1);
                 return (1 - t) * b + (t) * a;
             } else
                 return a;

queries/MeshQueries.cs

@@ -37,6 +37,19 @@ public static Frame3f NearestPointFrame(DMesh3 mesh, ISpatial spatial, Vector3d
         }
 
 
+        /// <summary>
+        /// Find distance from point to mesh
+        /// Returns interpolated vertex-normal frame if available, otherwise tri-normal frame.
+        /// </summary>
+        public static double NearestPointDistance(DMesh3 mesh, ISpatial spatial, Vector3d queryPoint, double maxDist = double.MaxValue)
+        {
+            int tid = spatial.FindNearestTriangle(queryPoint, maxDist);
+            if (tid == DMesh3.InvalidID)
+                return double.MaxValue;
+            return Math.Sqrt(TriangleDistance(mesh, tid, queryPoint).DistanceSquared);
+        }
+
+
         /// <summary>
         /// convenience function to construct a IntrRay3Triangle3 object for a mesh triangle
         /// </summary>

spatial/MeshSignedDistanceGrid.cs

@@ -24,6 +24,15 @@ namespace g3
     ///    that "narrow" if you have large triangles on a diagonal to grid axes
     /// 
     /// 
+    /// Potential optimizations:
+    ///  - Often we have a spatial data structure that would allow faster computation of the
+    ///    narrow-band distances (which become quite expensive if we want a wider band!)
+    ///    Not clear how to take advantage of this though. Perhaps we could have a binary
+    ///    grid that, in first pass, we set bits inside triangle bboxes to 1? Or perhaps
+    ///    same as current code, but we use spatial-dist, and so for each ijk we only compute once?
+    ///    (then have to test for computed value at each cell of each triangle...)
+    ///    
+    /// 
     /// This code is based on the C++ implementation found at https://github.com/christopherbatty/SDFGen
     /// Original license was public domain. 
     /// Permission granted by Christopher Batty to include C# port under Boost license.