diff --git a/TessellationAndVoxelizationGeometryLibrary/Constants.cs b/TessellationAndVoxelizationGeometryLibrary/Constants.cs index 5c570387e..ed76e952e 100644 --- a/TessellationAndVoxelizationGeometryLibrary/Constants.cs +++ b/TessellationAndVoxelizationGeometryLibrary/Constants.cs @@ -769,4 +769,32 @@ public enum BooleanOperationType XOR } + public enum QuadricType + { + Unknown, + Ellipsoid, + ImaginaryEllipsoid, + EllipticParaboloid, + HyperbolicParaboloid, + HyperboloidOneSheet, + HyperboloidTwoSheets, + EllipticCylinder, + ImaginaryEllipticCylinder, + ParabolicCylinder, + HyperbolicCylinder, + Cone, + ImaginaryCone, + IntersectingPlanes, + ImaginaryIntersectingPlanes, + ParallelPlanes, + ImaginaryParallelPlanes, + Plane, + Sphere, + Line, + Point, + NoPoint, + AllPoints, + Other + } + } \ No newline at end of file diff --git a/TessellationAndVoxelizationGeometryLibrary/Enclosure Operations/ConvexHull4D.cs b/TessellationAndVoxelizationGeometryLibrary/Enclosure Operations/ConvexHull4D.cs index e8f193c58..24d09f1c1 100644 --- a/TessellationAndVoxelizationGeometryLibrary/Enclosure Operations/ConvexHull4D.cs +++ b/TessellationAndVoxelizationGeometryLibrary/Enclosure Operations/ConvexHull4D.cs @@ -119,7 +119,7 @@ internal Vector4 GetNormal(bool tryToRepair) if (face == null) continue; var other = face.OwnedTetra == this ? face.OtherTetra : face.OwnedTetra; if (other == null) continue; - if (other.Normal.IsNull()) continue; + if (other.Normal.IsNegligible()) continue; normal += other.Normal; validNeighborCount++; } diff --git a/TessellationAndVoxelizationGeometryLibrary/Enclosure Operations/MinimumCircleCylinder.cs b/TessellationAndVoxelizationGeometryLibrary/Enclosure Operations/MinimumCircleCylinder.cs index f17dc98f2..1185bd2f3 100644 --- a/TessellationAndVoxelizationGeometryLibrary/Enclosure Operations/MinimumCircleCylinder.cs +++ b/TessellationAndVoxelizationGeometryLibrary/Enclosure Operations/MinimumCircleCylinder.cs @@ -184,7 +184,7 @@ private static Circle FindCircle(Vector2[] points) var minRadiusSqd = double.PositiveInfinity; if (Circle.CreateFrom3Points(points[0], points[1], points[2], out tempCircle) && (points[3] - tempCircle.Center).LengthSquared() <= tempCircle.RadiusSquared) - { // this one uses IsGreaterThanNonNegligible to prevent infinite cycling when more points are on the circle + { circle = tempCircle; minRadiusSqd = circle.RadiusSquared; } diff --git a/TessellationAndVoxelizationGeometryLibrary/Enclosure Operations/MinimumGaussSpherePlane.cs b/TessellationAndVoxelizationGeometryLibrary/Enclosure Operations/MinimumGaussSpherePlane.cs index 9a0f4af4c..0caad174d 100644 --- a/TessellationAndVoxelizationGeometryLibrary/Enclosure Operations/MinimumGaussSpherePlane.cs +++ b/TessellationAndVoxelizationGeometryLibrary/Enclosure Operations/MinimumGaussSpherePlane.cs @@ -142,7 +142,7 @@ private static Plane CreatePlaneFrom2Points(Vector3 p1, Vector3 p2, Vector3 orie { var center = (p1 + p2) / 2; var plane = new Plane(center, center); - if (plane.Normal.IsNull()) + if (plane.Normal.IsNegligible()) { // the two points midpoint happens to be the origin plane.Normal = orientingVector.Normalize(); plane.DistanceToOrigin = 0; @@ -154,6 +154,5 @@ private static Plane CreatePlaneFrom2Points(Vector3 p1, Vector3 p2, Vector3 orie } return plane; } - } } diff --git a/TessellationAndVoxelizationGeometryLibrary/Enclosure Operations/MinimumSphere.cs b/TessellationAndVoxelizationGeometryLibrary/Enclosure Operations/MinimumSphere.cs index 300bb5b18..21be46e97 100644 --- a/TessellationAndVoxelizationGeometryLibrary/Enclosure Operations/MinimumSphere.cs +++ b/TessellationAndVoxelizationGeometryLibrary/Enclosure Operations/MinimumSphere.cs @@ -117,79 +117,77 @@ private static sphereLight FirstSphere(Vector3[] points) // first check diametrically opposing points // 0,1 & check 2 & 3 var sphere = CreateFrom2Points(points[0], points[1]); - if ((points[2] - sphere.Center).LengthSquared() <= sphere.RadiusSquared - && (points[3] - sphere.Center).LengthSquared() <= sphere.RadiusSquared) + if (PointIsInSphere(points[2], sphere) && PointIsInSphere(points[3], sphere)) return sphere; // 0,2 & check 1 & 3 sphere = CreateFrom2Points(points[0], points[2]); - if ((points[1] - sphere.Center).LengthSquared() <= sphere.RadiusSquared - && (points[3] - sphere.Center).LengthSquared() <= sphere.RadiusSquared) + if (PointIsInSphere(points[1], sphere) && PointIsInSphere(points[3], sphere)) return sphere; // 0,3 & check 1 & 2 sphere = CreateFrom2Points(points[0], points[3]); - if ((points[1] - sphere.Center).LengthSquared() <= sphere.RadiusSquared - && (points[2] - sphere.Center).LengthSquared() <= sphere.RadiusSquared) + if (PointIsInSphere(points[1], sphere) && PointIsInSphere(points[2], sphere)) return sphere; // 1,2 & check 0 & 3 sphere = CreateFrom2Points(points[1], points[2]); - if ((points[0] - sphere.Center).LengthSquared() <= sphere.RadiusSquared - && (points[3] - sphere.Center).LengthSquared() <= sphere.RadiusSquared) + if (PointIsInSphere(points[0], sphere) && PointIsInSphere(points[3], sphere)) return sphere; // 1,3 & check 0 & 2 sphere = CreateFrom2Points(points[1], points[3]); - if ((points[0] - sphere.Center).LengthSquared() <= sphere.RadiusSquared - && (points[2] - sphere.Center).LengthSquared() <= sphere.RadiusSquared) + if (PointIsInSphere(points[0], sphere) && PointIsInSphere(points[2], sphere)) return sphere; // 2,3 & check 0 & 1 sphere = CreateFrom2Points(points[2], points[3]); - if ((points[0] - sphere.Center).LengthSquared() <= sphere.RadiusSquared - && (points[1] - sphere.Center).LengthSquared() <= sphere.RadiusSquared) + if (PointIsInSphere(points[0], sphere) && PointIsInSphere(points[1], sphere)) return sphere; + sphere = default; var minRadiusSqd = double.PositiveInfinity; // now check 3-point spheres. here we need to find the smallest sphere! this wasn't the // case for the above diametrically opposed points (since the two defining points of the sphere // are farthest apart, but we need to check this here as 3 nearly collinear points will make // a huge sphere // 0,1,2 & check 3 - sphere = CreateFrom3Points(points[0], points[1], points[2]); - if ((points[3] - sphere.Center).LengthSquared() <= sphere.RadiusSquared) + var tempSphere3 = CreateFrom3Points(points[0], points[1], points[2]); + if (PointIsInSphere(points[3], tempSphere3)) + { + sphere = tempSphere3; minRadiusSqd = sphere.RadiusSquared; + } // 0,1,3 & check 2 var swap3And2 = false; - var tempSphere = CreateFrom3Points(points[0], points[1], points[3]); - if ((points[2] - tempSphere.Center).LengthSquared() <= tempSphere.RadiusSquared - && tempSphere.RadiusSquared < minRadiusSqd) + tempSphere3 = CreateFrom3Points(points[0], points[1], points[3]); + if (PointIsInSphere(points[2], tempSphere3) + && tempSphere3.RadiusSquared < minRadiusSqd) { - sphere = tempSphere; + sphere = tempSphere3; swap3And2 = true; minRadiusSqd = sphere.RadiusSquared; } // 0,2,3 & check 1 var swap3And1 = false; - tempSphere = CreateFrom3Points(points[0], points[2], points[3]); - if ((points[1] - tempSphere.Center).LengthSquared() <= tempSphere.RadiusSquared - && tempSphere.RadiusSquared < minRadiusSqd) + tempSphere3 = CreateFrom3Points(points[0], points[2], points[3]); + if (PointIsInSphere(points[1], tempSphere3) + && tempSphere3.RadiusSquared < minRadiusSqd) { - sphere = tempSphere; + sphere = tempSphere3; swap3And1 = true; minRadiusSqd = sphere.RadiusSquared; } // 1,2,3 & check 0 var swap3And0 = false; - tempSphere = CreateFrom3Points(points[1], points[2], points[3]); - if ((points[0] - tempSphere.Center).LengthSquared() <= tempSphere.RadiusSquared - && tempSphere.RadiusSquared < minRadiusSqd) + tempSphere3 = CreateFrom3Points(points[1], points[2], points[3]); + if (PointIsInSphere(points[0], tempSphere3) + && tempSphere3.RadiusSquared < minRadiusSqd) { - sphere = tempSphere; + sphere = tempSphere3; swap3And0 = true; minRadiusSqd = sphere.RadiusSquared; } var fourPointIsBest = false; - tempSphere = CreateFrom4Points(points[0], points[1], points[2], points[3]); - if (tempSphere.RadiusSquared < minRadiusSqd) + var tempSphere4 = CreateFrom4Points(points[0], points[1], points[2], points[3]); + if (tempSphere4.RadiusSquared < minRadiusSqd) { - sphere = tempSphere; + sphere = tempSphere4; fourPointIsBest = true; } if (!fourPointIsBest) @@ -198,143 +196,133 @@ private static sphereLight FirstSphere(Vector3[] points) else if (swap3And1) Global.SwapItemsInList(3, 1, points); else if (swap3And2) Global.SwapItemsInList(3, 2, points); } + if (sphere.RadiusSquared == 0.0) + sphere = double.IsNaN(tempSphere3.RadiusSquared) || double.IsInfinity(tempSphere3.RadiusSquared) + ? tempSphere4 : tempSphere3; return sphere; } + private static sphereLight FindSphere(Vector3[] points) { + // first check diametrically opposing points: the good news is that the zeroth // point must be in the set of points defining the sphere. The bad news is that // we need to include cases up to 5 points! // 0,1 & check 2, 3 & 4 var sphere = CreateFrom2Points(points[0], points[1]); - if ((points[2] - sphere.Center).LengthSquared() <= sphere.RadiusSquared - && (points[3] - sphere.Center).LengthSquared() <= sphere.RadiusSquared - && (points[4] - sphere.Center).LengthSquared() <= sphere.RadiusSquared) + if (PointIsInSphere(points[2], sphere) && PointIsInSphere(points[3], sphere) && PointIsInSphere(points[4], sphere)) return sphere; // 0,2 & check 1, 3 & 4 sphere = CreateFrom2Points(points[0], points[2]); - if ((points[1] - sphere.Center).LengthSquared() <= sphere.RadiusSquared - && (points[3] - sphere.Center).LengthSquared() <= sphere.RadiusSquared - && (points[4] - sphere.Center).LengthSquared() <= sphere.RadiusSquared) + if (PointIsInSphere(points[1], sphere) && PointIsInSphere(points[3], sphere) && PointIsInSphere(points[4], sphere)) return sphere; // 0,3 & check 1, 2 & 4 sphere = CreateFrom2Points(points[0], points[3]); - if ((points[1] - sphere.Center).LengthSquared() <= sphere.RadiusSquared - && (points[2] - sphere.Center).LengthSquared() <= sphere.RadiusSquared - && (points[4] - sphere.Center).LengthSquared() <= sphere.RadiusSquared) + if (PointIsInSphere(points[1], sphere) && PointIsInSphere(points[2], sphere) && PointIsInSphere(points[4], sphere)) return sphere; // 0,4 & check 1, 2 & 3 sphere = CreateFrom2Points(points[0], points[4]); - if ((points[1] - sphere.Center).LengthSquared() <= sphere.RadiusSquared - && (points[2] - sphere.Center).LengthSquared() <= sphere.RadiusSquared - && (points[3] - sphere.Center).LengthSquared() <= sphere.RadiusSquared) + if (PointIsInSphere(points[1], sphere) && PointIsInSphere(points[2], sphere) && PointIsInSphere(points[3], sphere)) return sphere; + sphere = default; var minRadiusSqd = double.PositiveInfinity; // now check 3-point spheres. here we need to find the smallest sphere! this wasn't the // case for the above diametrically opposed points (since the two defining points of the sphere // are farthest apart, but we need to check this here as 3 nearly collinear points will make // a huge sphere // 0,1,2 & check 3 & 4 - sphere = CreateFrom3Points(points[0], points[1], points[2]); - if ((points[3] - sphere.Center).LengthSquared() <= sphere.RadiusSquared - && (points[4] - sphere.Center).LengthSquared() <= sphere.RadiusSquared) + var tempSphere3 = CreateFrom3Points(points[0], points[1], points[2]); + if (PointIsInSphere(points[3], tempSphere3) && PointIsInSphere(points[4], tempSphere3)) + { + sphere = tempSphere3; minRadiusSqd = sphere.RadiusSquared; + } // 0,1,3 & check 2 & 4 var swap3And2 = false; - var tempSphere = CreateFrom3Points(points[0], points[1], points[3]); - if ((points[2] - tempSphere.Center).LengthSquared() <= tempSphere.RadiusSquared - && (points[4] - tempSphere.Center).LengthSquared() <= tempSphere.RadiusSquared - && tempSphere.RadiusSquared < minRadiusSqd) + tempSphere3 = CreateFrom3Points(points[0], points[1], points[3]); + if (PointIsInSphere(points[2], tempSphere3) && PointIsInSphere(points[4], tempSphere3) + && tempSphere3.RadiusSquared < minRadiusSqd) { - sphere = tempSphere; + sphere = tempSphere3; swap3And2 = true; minRadiusSqd = sphere.RadiusSquared; } // 0,1,4 & check 2 & 3 var swap4And2 = false; - tempSphere = CreateFrom3Points(points[0], points[1], points[4]); - if ((points[2] - tempSphere.Center).LengthSquared() <= tempSphere.RadiusSquared - && (points[3] - tempSphere.Center).LengthSquared() <= tempSphere.RadiusSquared - && tempSphere.RadiusSquared < minRadiusSqd) + tempSphere3 = CreateFrom3Points(points[0], points[1], points[4]); + if (PointIsInSphere(points[2], tempSphere3) && PointIsInSphere(points[3], tempSphere3) + && tempSphere3.RadiusSquared < minRadiusSqd) { - sphere = tempSphere; + sphere = tempSphere3; swap4And2 = true; minRadiusSqd = sphere.RadiusSquared; } // 0,2,3 & check 1 & 4 var swap3And1 = false; - tempSphere = CreateFrom3Points(points[0], points[2], points[3]); - if ((points[1] - tempSphere.Center).LengthSquared() <= tempSphere.RadiusSquared - && (points[4] - tempSphere.Center).LengthSquared() <= tempSphere.RadiusSquared - && tempSphere.RadiusSquared < minRadiusSqd) + tempSphere3 = CreateFrom3Points(points[0], points[2], points[3]); + if (PointIsInSphere(points[1], tempSphere3) && PointIsInSphere(points[4], tempSphere3) + && tempSphere3.RadiusSquared < minRadiusSqd) { - sphere = tempSphere; + sphere = tempSphere3; swap3And1 = true; minRadiusSqd = sphere.RadiusSquared; } // 0,2,4 & check 1 & 3 var swap1With4 = false; - tempSphere = CreateFrom3Points(points[0], points[2], points[4]); - if ((points[1] - tempSphere.Center).LengthSquared() <= tempSphere.RadiusSquared - && (points[3] - tempSphere.Center).LengthSquared() <= tempSphere.RadiusSquared - && tempSphere.RadiusSquared < minRadiusSqd) + tempSphere3 = CreateFrom3Points(points[0], points[2], points[4]); + if (PointIsInSphere(points[1], tempSphere3) && PointIsInSphere(points[3], tempSphere3) + && tempSphere3.RadiusSquared < minRadiusSqd) { - sphere = tempSphere; + sphere = tempSphere3; swap1With4 = true; minRadiusSqd = sphere.RadiusSquared; } // 0,3,4 & check 1 & 2 var swap12With34 = false; - tempSphere = CreateFrom3Points(points[0], points[3], points[4]); - if ((points[1] - tempSphere.Center).LengthSquared() <= tempSphere.RadiusSquared - && (points[2] - tempSphere.Center).LengthSquared() <= tempSphere.RadiusSquared - && tempSphere.RadiusSquared < minRadiusSqd) + tempSphere3 = CreateFrom3Points(points[0], points[3], points[4]); + if (PointIsInSphere(points[1], tempSphere3) && PointIsInSphere(points[2], tempSphere3) + && tempSphere3.RadiusSquared < minRadiusSqd) { - sphere = tempSphere; + sphere = tempSphere3; swap12With34 = true; minRadiusSqd = sphere.RadiusSquared; } // now the 4-point spheres var fourPointIsBest = false; // 0,1,2,3 & check 4 - tempSphere = CreateFrom4Points(points[0], points[1], points[2], points[3]); - if (!(points[4] - tempSphere.Center).LengthSquared().IsGreaterThanNonNegligible(tempSphere.RadiusSquared) - // this one uses IsGreaterThanNonNegligible to prevent infinite cycling when more points are on the sphere - && tempSphere.RadiusSquared < minRadiusSqd) + var tempSphere4 = CreateFrom4Points(points[0], points[1], points[2], points[3]); + if (PointIsInSphere(points[4], tempSphere4) && tempSphere4.RadiusSquared < minRadiusSqd) { - sphere = tempSphere; + sphere = tempSphere4; fourPointIsBest = true; minRadiusSqd = sphere.RadiusSquared; } // 0,1,2,4 & check 3 var swap4And3 = false; - tempSphere = CreateFrom4Points(points[0], points[1], points[2], points[4]); - if ((points[3] - tempSphere.Center).LengthSquared() <= tempSphere.RadiusSquared - && tempSphere.RadiusSquared < minRadiusSqd) + tempSphere4 = CreateFrom4Points(points[0], points[1], points[2], points[4]); + if (PointIsInSphere(points[3], tempSphere4) && tempSphere4.RadiusSquared < minRadiusSqd) { - sphere = tempSphere; + sphere = tempSphere4; swap4And3 = fourPointIsBest = true; minRadiusSqd = sphere.RadiusSquared; } // 0,1,3,4 & check 2 var swap4With2 = false; - tempSphere = CreateFrom4Points(points[0], points[1], points[3], points[4]); - if ((points[2] - tempSphere.Center).LengthSquared() <= tempSphere.RadiusSquared - && tempSphere.RadiusSquared < minRadiusSqd) + tempSphere4 = CreateFrom4Points(points[0], points[1], points[3], points[4]); + if (PointIsInSphere(points[2], tempSphere4) && tempSphere4.RadiusSquared < minRadiusSqd) { - sphere = tempSphere; + sphere = tempSphere4; swap4With2 = fourPointIsBest = true; minRadiusSqd = sphere.RadiusSquared; } // 0,2,3,4 & check 1 var swap4With1 = false; - tempSphere = CreateFrom4Points(points[0], points[2], points[3], points[4]); - if ((points[1] - tempSphere.Center).LengthSquared() <= tempSphere.RadiusSquared - && tempSphere.RadiusSquared < minRadiusSqd) + tempSphere4 = CreateFrom4Points(points[0], points[2], points[3], points[4]); + if (PointIsInSphere(points[1], tempSphere4) && tempSphere4.RadiusSquared < minRadiusSqd) { - sphere = tempSphere; + sphere = tempSphere4; swap4With1 = fourPointIsBest = true; //minRadiusSqd = sphere.RadiusSquared; // don't need this anymore...no more checks @@ -359,9 +347,15 @@ private static sphereLight FindSphere(Vector3[] points) else if (swap4And2) Global.SwapItemsInList(4, 2, points); else if (swap3And2) Global.SwapItemsInList(3, 2, points); } + if (sphere.RadiusSquared == 0.0) + sphere = double.IsNaN(tempSphere3.RadiusSquared) || double.IsInfinity(tempSphere3.RadiusSquared) + ? tempSphere4 : tempSphere3; return sphere; } + private static bool PointIsInSphere(Vector3 point, sphereLight s) + => (point - s.Center).LengthSquared() <= 1.0001 * s.RadiusSquared; + private static sphereLight CreateFrom4Points(Vector3 p1, Vector3 p2, Vector3 p3, Vector3 p4) diff --git a/TessellationAndVoxelizationGeometryLibrary/Implicits and Primitives/Primitive Surfaces/GeneralQuadric.cs b/TessellationAndVoxelizationGeometryLibrary/Implicits and Primitives/Primitive Surfaces/GeneralQuadric.cs index a94d86562..be2507e53 100644 --- a/TessellationAndVoxelizationGeometryLibrary/Implicits and Primitives/Primitive Surfaces/GeneralQuadric.cs +++ b/TessellationAndVoxelizationGeometryLibrary/Implicits and Primitives/Primitive Surfaces/GeneralQuadric.cs @@ -11,11 +11,13 @@ // // // *********************************************************************** +using Microsoft.Extensions.Options; using StarMathLib; using System; using System.Collections.Generic; using System.Linq; using System.Runtime.CompilerServices; +using System.Security.Principal; using System.Text.Json.Serialization; namespace TVGL @@ -115,6 +117,36 @@ public double W init => w = value; } private double w; + + /// + /// Specifies the type of quadric represented. + /// + private QuadricType type = QuadricType.Unknown; + public QuadricType Type + { + get + { + if (type == QuadricType.Unknown) + type = SetQuadricType(XSqdCoeff, YSqdCoeff, ZSqdCoeff, XYCoeff, + XZCoeff, YZCoeff, XCoeff, YCoeff, ZCoeff, W); + return type; + } + set { type = value; } + } + + + /// + /// Represent quadric specific coefficients + /// For plane, a,b,c are plane normal components and d is the distance from the origin. + /// For parallel planes, a,b,c, and d describe the mid-plane while e is the distance between the two planes. + /// ... + /// + public double a; + public double b; + public double c; + public double d; + public double e; + /// /// Gets all the coefficients as an enumerable in the order listed above. public IEnumerable Coefficients @@ -284,13 +316,15 @@ public override bool PointIsInside(Vector3 x) /// A Vector3. public override Vector3 GetNormalAtPoint(Vector3 point) { - var gradient = GetGradient(point).Normalize(); + if (point.IsNull()) return Vector3.Null; + var gradient = GetGradient(point); - if (!gradient.Length().IsPositiveNonNegligible()) return GetNormalAtUmbilicPoint(point); + if (gradient.Length().IsNegligible()) + gradient = GetNormalAtStationaryPoint(point); if (IsPositive.GetValueOrDefault(true)) - return gradient; - else return -gradient; + return gradient.Normalize(); + else return -gradient.Normalize(); } /// @@ -314,23 +348,26 @@ public double QuadricValue(Vector3 point) + W); } - private Vector3 GetNormalAtUmbilicPoint(Vector3 point) + private Vector3 GetNormalAtStationaryPoint(Vector3 point) { - // At an umbilic point, the normal is zero, so we need to find the direction of maximum increase. + // At a stationary point, the normal is zero, so we need to find the direction of maximum increase. // We assume that the quadric field is not flat everywhere - // Return a normal at a point slightly offset from the umbilic point - Vector3 normal = Vector3.Zero; + // Return a normal at a point slightly offset from the stationary point double delta = 1E-6; // small offset distance int iters = 0; - while (!normal.Length().IsPositiveNonNegligible() && iters < 100) + Vector3 gradient = Vector3.Zero; + while (!gradient.Length().IsPositiveNonNegligible() && iters < 100) { var random = new Random(); - normal = GetNormalAtPoint(point + delta * new Vector3(random.NextDouble(), random.NextDouble(), random.NextDouble()).Normalize()); + gradient = GetGradient(point + delta * new Vector3(random.NextDouble(), random.NextDouble(), random.NextDouble())); iters++; } - return normal; + if (IsPositive.GetValueOrDefault(true)) + return gradient; + else return -gradient; } + private Vector3 GetNearbyPointOnQuadric(Vector3 anchor) { Vector3 normal = GetNormalAtPoint(anchor); @@ -343,13 +380,10 @@ private Vector3 GetNearbyPointOnQuadric(Vector3 anchor) while (!intersections.GetEnumerator().MoveNext() && iters++ < 100) { //This is the distance to the tangent point on the normal to the level set. - t = -normal.Dot(normal) / GetNormalAtPoint(normal).Dot(normal); + t = normal.Dot(normal) / ((new Vector3(XCoeff, YCoeff, ZCoeff) - GetNormalAtPoint(normal)).Dot(normal)); + var tangencyCheck = normal.Dot(GetNormalAtPoint(newAnchor + t * normal)); newAnchor = newAnchor + t * normal; normal = GetNormalAtPoint(newAnchor); - if (normal.Length() == 0) - { - normal = GetNormalAtUmbilicPoint(newAnchor); - } intersections = LineIntersection(newAnchor, normal); } if (!intersections.GetEnumerator().MoveNext()) return newAnchor; @@ -396,12 +430,13 @@ public override Vector3 ClosestPointOnSurfaceToPoint(Vector3 point) return closestPt; } + public override double DistanceToPoint(Vector3 point) { Vector3 closestPt = ClosestPointOnSurfaceToPoint(point); //scaling the quadric value by the norm of the normal vector to get the approximate distance locally - if (closestPt == Vector3.Null) return DistanceToPointQuick(point); + if (closestPt.IsNull()) return DistanceToPointQuick(point); return Math.Sign(QuadricValue(point)) * point.Distance(closestPt); } @@ -420,7 +455,7 @@ public double DistanceToPointQuick(Vector3 point) protected override void CalculateIsPositive() { if (Faces == null || !Faces.Any() || Area.IsNegligible()) return; - isPositive = LargestFace.Normal.Dot(LargestFace.Center - Center) > 0; + isPositive = LargestFace.Normal.Dot(LargestFace.Center - StationaryPoint) > 0; } protected override void SetPrimitiveLimits() @@ -433,33 +468,28 @@ protected override void SetPrimitiveLimits() /// /// Finds the points on the quadric that has a gradient in the same direction as /// the specified gradient vector. This input does not have to be normalized. - /// If onlyAlignedDir is true, then only the points where the gradient is in the - /// same direction as the input gradient will be returned. - /// If false, then points where the gradient is in the opposite direction will also be returned. /// - /// The method determines the anchor point and direction of the line using the quadric's - /// coefficients and the provided gradient vector. Ensure that the gradient vector is non-zero to avoid - /// undefined behavior. - /// The gradient vector that defines the direction of the line for which intersection points with the surface - /// are calculated. Must not be the zero vector. - /// An enumerable collection of points representing the intersection points of the line - /// with the quadric surface. The collection may be empty if there are no intersections. - public IEnumerable PointsWithGradient(Vector3 gradient, bool onlyAlignedDir = false) + /// + /// + /// + public bool PointsAtGivenGradient(Vector3 gradient, out Vector3 point) { var c = new Vector3(-XCoeff, -YCoeff, -ZCoeff); var gradMatrix = new Matrix3x3( 2 * XSqdCoeff, XYCoeff, XZCoeff, XYCoeff, 2 * YSqdCoeff, YZCoeff, XZCoeff, YZCoeff, 2 * ZSqdCoeff); - Matrix3x3.Invert(gradMatrix, out var invGradMatrix); + if (!Matrix3x3.Invert(gradMatrix, out var invGradMatrix)) + { + point = Vector3.Null; + return false; + } var anchor = invGradMatrix.Multiply(c); var dir = invGradMatrix.Multiply(gradient).Normalize(); - if (onlyAlignedDir) - return LineIntersection(anchor, dir) - .Where(x => GetGradient(x.intersection).Dot(gradient) > 0) - .Select(x => x.intersection); - else - return LineIntersection(anchor, dir).Select(x => x.intersection); + point = LineIntersection(anchor, dir) + .Where(x => GetGradient(x.intersection).Dot(gradient) > 0) + .Select(x => x.intersection).First(); + return true; } /// @@ -482,40 +512,47 @@ public IEnumerable PointsWithGradient(Vector3 gradient, bool onlyAligne + XYCoeff * anchor.X * anchor.Y + XZCoeff * anchor.X * anchor.Z + YZCoeff * anchor.Y * anchor.Z + XCoeff * anchor.X + YCoeff * anchor.Y + ZCoeff * anchor.Z + W; (var root1, var root2) = PolynomialSolve.Quadratic(a, b, c); - if (root1.IsRealNumber && root1.Real.IsPracticallySame(root2.Real)) + yield return (anchor + root1.Real * direction, root1.Real); + else if (root1.IsRealNumber) { - var t = 0.5 * (root1.Real + root2.Real); - yield return (anchor + t * direction, root1.Real); - yield break; - } - if (root1.IsRealNumber) yield return (anchor + root1.Real * direction, root1.Real); - if (root2.IsRealNumber) yield return (anchor + root2.Real * direction, root2.Real); + } + //else + // yield return (Vector3.Null, double.PositiveInfinity); + } + + public bool PracticallyOnSurface(Vector3 point, double tolerance = Constants.BaseTolerance) + { + var magGradient = GetGradient(point).Length(); + // Use max(gradient, 1) to prevent the tolerance from collapsing to zero + // near vertices/tips where the gradient vanishes (e.g. tip of a hyperboloid) + return QuadricValue(point).IsNegligible(tolerance * Math.Max(magGradient, 1.0)); } /// - /// The is the mathematical center of the quadric. It is not the centroid of the faces. - /// It is the location where the gradient of the quadric function is zero. + /// The is the mathematical center of the quadric is referred to as the stationary point. + /// It is not the centroid of the faces, but rather the point where the quadric function has + /// a zero gradient. /// It would be the center of an ellipsoid, but for other quadrics, it is ...? /// Desmos it up! /// /// [JsonIgnore] - public Vector3 Center + public Vector3 StationaryPoint { get { - if (center.IsNull()) + if (staionaryPoint.IsNull()) { Matrix3x3 coeffs = GetHessian(); // the center is defined where the gradient is zero for the quadric function // this produces a system of 3 equations with 3 unknowns var b = new Vector3(-XCoeff, -YCoeff, -ZCoeff); - center = coeffs.Solve(b); + staionaryPoint = coeffs.Solve(b); } - return center; + return staionaryPoint; } } @@ -527,14 +564,14 @@ private Vector3 GetGradient(Vector3 point) 2 * ZSqdCoeff * point.Z + XZCoeff * point.X + YZCoeff * point.Y + ZCoeff); } - private Matrix3x3 GetHessian() + public Matrix3x3 GetHessian() { return new Matrix3x3(2 * XSqdCoeff, XYCoeff, XZCoeff, XYCoeff, 2 * YSqdCoeff, YZCoeff, XZCoeff, YZCoeff, 2 * ZSqdCoeff); } - Vector3 center = Vector3.Null; + Vector3 staionaryPoint = Vector3.Null; [JsonIgnore] public Vector3 Axis1 @@ -578,6 +615,21 @@ public Vector3 Axis3 Vector3 axis3 = Vector3.Null; + + public void Negate() + { + xSqdCoeff = -xSqdCoeff; + ySqdCoeff = -ySqdCoeff; + zSqdCoeff = -zSqdCoeff; + xyCoeff = -xyCoeff; + xzCoeff = -xzCoeff; + yzCoeff = -yzCoeff; + xCoeff = -xCoeff; + yCoeff = -yCoeff; + zCoeff = -zCoeff; + w = -w; + } + [JsonIgnore] public override string KeyString => "Quadric|" + XSqdCoeff.ToString("F5") + "|" + YSqdCoeff.ToString("F5") + "|" + ZSqdCoeff.ToString("F5") + "|" + @@ -932,5 +984,145 @@ public static (Vector3 tangent, Vector3 normal, Vector3 center, double radius) C return (tangent, curvatureVector, point + radius * curvatureVector, radius); } + + public void RemoveSmallCoefficents(double tolerance = 1E-3) + { + this.xSqdCoeff = xSqdCoeff.IsNegligible(tolerance) ? 0 : xSqdCoeff; + this.ySqdCoeff = ySqdCoeff.IsNegligible(tolerance) ? 0 : ySqdCoeff; + this.zSqdCoeff = zSqdCoeff.IsNegligible(tolerance) ? 0 : zSqdCoeff; + this.xyCoeff = xyCoeff.IsNegligible(tolerance) ? 0 : xyCoeff; + this.xzCoeff = xzCoeff.IsNegligible(tolerance) ? 0 : xzCoeff; + this.yzCoeff = yzCoeff.IsNegligible(tolerance) ? 0 : yzCoeff; + this.xCoeff = xCoeff.IsNegligible(tolerance) ? 0 : xCoeff; + this.yCoeff = yCoeff.IsNegligible(tolerance) ? 0 : yCoeff; + this.zCoeff = zCoeff.IsNegligible(tolerance) ? 0 : zCoeff; + this.w = w.IsNegligible(tolerance) ? 0 : w; + } + + public static QuadricType SetQuadricType(GeneralQuadric q, double tol = 1E-12) + => SetQuadricType(q.XSqdCoeff, q.YSqdCoeff, q.ZSqdCoeff, q.XYCoeff, + q.XZCoeff, q.YZCoeff, q.XCoeff, q.YCoeff, q.ZCoeff, q.W, tol); + + public static QuadricType SetQuadricType(double XSqdCoeff, double YSqdCoeff, double ZSqdCoeff, double XYCoeff, + double XZCoeff, double YZCoeff, double XCoeff, double YCoeff, double ZCoeff, double W, double tol = 1E-12) + { + var A = new Matrix3x3(XSqdCoeff, XYCoeff / 2, XZCoeff / 2, XYCoeff / 2, YSqdCoeff / 2, + YZCoeff / 2, XZCoeff / 2, YZCoeff / 2, ZSqdCoeff); + var M = new Matrix4x4(XSqdCoeff, XYCoeff / 2, XZCoeff / 2, XCoeff / 2, + XYCoeff / 2, YSqdCoeff, YZCoeff / 2, YCoeff / 2, + XZCoeff / 2, YZCoeff / 2, ZSqdCoeff, ZCoeff / 2, + XCoeff / 2, YCoeff / 2, ZCoeff / 2, W); + + var AEigvals = A.GetEigenValues().Select(e => e.Real); + var MEigvals = M.GetEigenValues().Select(e => e.Real); + int numPosAEigvals = AEigvals.Count(e => e.IsPositiveNonNegligible(tol)); + int Arank = AEigvals.Count(e => !e.IsNegligible(tol)); + int Mrank = MEigvals.Count(e => !e.IsNegligible(tol)); + bool MDetSign = (MEigvals.Count(e => e < 0) % 2) == 0; + bool ANonZeroEigValsSameSign = AEigvals.Where(e => !e.IsNegligible(tol)).All(e => e > 0) || AEigvals.Where(e => !e.IsNegligible(tol)).All(e => e < 0); + bool MNonZeroEigValsSameSign = MEigvals.Where(e => !e.IsNegligible(tol)).All(e => e > 0) || MEigvals.Where(e => !e.IsNegligible(tol)).All(e => e < 0); + + if (Arank == 3) + { + if (Mrank == 4) + { + if (!MDetSign) + { + if (ANonZeroEigValsSameSign) + return QuadricType.Ellipsoid; + else return QuadricType.HyperboloidTwoSheets; + } + else + { + if (ANonZeroEigValsSameSign) return QuadricType.ImaginaryEllipsoid; + else return QuadricType.HyperboloidOneSheet; + } + } + + else if (Mrank == 3) + { + if (ANonZeroEigValsSameSign) return QuadricType.ImaginaryCone; + else return QuadricType.Cone; + } + } + + else if (Arank == 2) + { + if (Mrank == 4) + { + if (ANonZeroEigValsSameSign) return QuadricType.EllipticParaboloid; + else return QuadricType.HyperbolicParaboloid; + } + else if (Mrank == 3) + { + if (ANonZeroEigValsSameSign) + { + if (MNonZeroEigValsSameSign) return QuadricType.ImaginaryEllipticCylinder; + else return QuadricType.EllipticCylinder; + } + else return QuadricType.HyperbolicCylinder; + } + else if (Mrank == 2) + { + if (ANonZeroEigValsSameSign) return QuadricType.ImaginaryIntersectingPlanes; + else return QuadricType.IntersectingPlanes; + } + } + + else if (Arank == 1) + { + if (Mrank == 3) return QuadricType.ParabolicCylinder; + else if (Mrank == 2) + { + if (MNonZeroEigValsSameSign) return QuadricType.ImaginaryParallelPlanes; + else return QuadricType.ParallelPlanes; + } + else if (Mrank == 1) return QuadricType.Plane; + } + else + { + if (Mrank == 2) return QuadricType.Plane; + } + return QuadricType.Unknown; + } + + public void SetTypeSpecificCoefficients() + { + if (Type == QuadricType.Plane) + { + double normalizer = Math.Sqrt(XCoeff * XCoeff + YCoeff * YCoeff + ZCoeff * ZCoeff); + a = XCoeff / normalizer; + b = YCoeff / normalizer; + c = ZCoeff / normalizer; + d = W / normalizer; + } + if (Type == QuadricType.ParallelPlanes) + { + double normalizer = Math.Sqrt(Math.Abs(XSqdCoeff + YSqdCoeff + ZSqdCoeff)); + a = Math.Sqrt(Math.Abs(XSqdCoeff)) / normalizer; + d = XCoeff / (2 * a); + b = YCoeff / (2 * d); + c = ZCoeff / (2 * d); + e = 2 * Math.Sqrt((XCoeff * XCoeff) / (4 * XSqdCoeff) - W); + + } + } + + /// + /// Constructs and returns a transformation matrix representing the coefficients of the quadric surface in 3D + /// space. + /// + /// The returned matrix encodes the coefficients for the X, Y, and Z dimensions as + /// currently set in the instance. This matrix can be used for further geometric computations or transformations + /// involving the quadric surface. + /// A containing the coefficients for the quadric surface transformation. + public Matrix4x4 GetCoefficientMatrix() + { + return new Matrix4x4( + XSqdCoeff, 0.5 * XYCoeff, 0.5 * XZCoeff, 0.5 * XCoeff, + 0.5 * XYCoeff, YSqdCoeff, 0.5 * YZCoeff, 0.5 * YCoeff, + 0.5 * XZCoeff, 0.5 * YZCoeff, ZSqdCoeff, 0.5 * ZCoeff, + 0.5 * XCoeff, 0.5 * YCoeff, 0.5 * ZCoeff, W); + } } } \ No newline at end of file diff --git a/TessellationAndVoxelizationGeometryLibrary/Implicits and Primitives/Primitive Surfaces/PrimitiveSurface.cs b/TessellationAndVoxelizationGeometryLibrary/Implicits and Primitives/Primitive Surfaces/PrimitiveSurface.cs index 5c608287f..3f00acbf5 100644 --- a/TessellationAndVoxelizationGeometryLibrary/Implicits and Primitives/Primitive Surfaces/PrimitiveSurface.cs +++ b/TessellationAndVoxelizationGeometryLibrary/Implicits and Primitives/Primitive Surfaces/PrimitiveSurface.cs @@ -67,7 +67,7 @@ public void SetFacesAndVertices(IEnumerable faces, bool connectFac ResetFaceDependentValues(); ResetFaceDependentConnectivity(); if (faces == null) Faces = new HashSet(); - else Faces = new HashSet(faces); + else Faces = faces.ToHashSet(); FaceIndices = Faces.Select(f => f.IndexInList).ToArray(); if (connectFacesToPrimitive) foreach (var face in Faces) @@ -170,6 +170,7 @@ public virtual double CalculateMeanSquareError(IEnumerable points) foreach (var c in points) { var d = DistanceToPoint(c); + if (double.IsNaN(d)) ; mse += d * d; n++; } @@ -442,9 +443,10 @@ public HashSet OuterEdges private HashSet _outerEdges; /// - /// Defines the inner outer edges. + /// Defines the inner and outer edges. This recreates the two hashsets within the PrimitiveSurface + /// by going through the faces stored in the primitive. /// - internal void DefineInnerOuterEdges() + public void DefineInnerOuterEdges() { MiscFunctions.DefineInnerOuterEdges(Faces, out _innerEdges, out _outerEdges); } @@ -504,42 +506,42 @@ public void AddFace(TriangleFace face) } else OuterEdges.Add(e); } - else //basically, this is for cases where edges are not yet defined. - { - var faceVertexIndices = new List { face.A.IndexInList, face.B.IndexInList, face.B.IndexInList, face.C.IndexInList, face.C.IndexInList, face.A.IndexInList }; - //this is kind of hacky, but the faceVertexIndices don't need to be in order, so we can just add the same vertex twice. The if statement below will catch it. - var outerEdgesToRemove = new List(); - foreach (var outerEdge in OuterEdges) - { - var vertexIndex1 = outerEdge.From.IndexInList; - var vertexIndex2 = outerEdge.To.IndexInList; - if (faceVertexIndices.Contains(vertexIndex1) && faceVertexIndices.Contains(vertexIndex2)) - { - faceVertexIndices.Remove(vertexIndex1); - faceVertexIndices.Remove(vertexIndex2); - outerEdge.UpdateWithNewFace(face); - outerEdgesToRemove.Add(outerEdge); - InnerEdges.Add(outerEdge); - if (faceVertexIndices.Count == 0) break; - } - } - foreach (var edge in outerEdgesToRemove) - OuterEdges.Remove(edge); - while (faceVertexIndices.Count > 0) - { - var vIndex1 = faceVertexIndices[0]; - faceVertexIndices.RemoveAt(0); - var vIndex2 = faceVertexIndices.First(index => index != vIndex1); - faceVertexIndices.Remove(vIndex2); - var v1 = face.Vertices.FindIndex(v => v.IndexInList == vIndex1); - var v2 = face.Vertices.FindIndex(v => v.IndexInList == vIndex2); - var step = v2 - v1; - if (step < 0) step += 3; - if (step == 1) - OuterEdges.Add(new Edge(face.A, face.B, face, null, false)); - else OuterEdges.Add(new Edge(face.B, face.A, face, null, false)); - } - } + //else //basically, this is for cases where edges are not yet defined. + //{ + // var faceVertexIndices = new List { face.A.IndexInList, face.B.IndexInList, face.B.IndexInList, face.C.IndexInList, face.C.IndexInList, face.A.IndexInList }; + // //this is kind of hacky, but the faceVertexIndices don't need to be in order, so we can just add the same vertex twice. The if statement below will catch it. + // var outerEdgesToRemove = new List(); + // foreach (var outerEdge in OuterEdges) + // { + // var vertexIndex1 = outerEdge.From.IndexInList; + // var vertexIndex2 = outerEdge.To.IndexInList; + // if (faceVertexIndices.Contains(vertexIndex1) && faceVertexIndices.Contains(vertexIndex2)) + // { + // faceVertexIndices.Remove(vertexIndex1); + // faceVertexIndices.Remove(vertexIndex2); + // outerEdge.UpdateWithNewFace(face); + // outerEdgesToRemove.Add(outerEdge); + // InnerEdges.Add(outerEdge); + // if (faceVertexIndices.Count == 0) break; + // } + // } + // foreach (var edge in outerEdgesToRemove) + // OuterEdges.Remove(edge); + // while (faceVertexIndices.Count > 0) + // { + // var vIndex1 = faceVertexIndices[0]; + // faceVertexIndices.RemoveAt(0); + // var vIndex2 = faceVertexIndices.First(index => index != vIndex1); + // faceVertexIndices.Remove(vIndex2); + // var v1 = face.Vertices.FindIndex(v => v.IndexInList == vIndex1); + // var v2 = face.Vertices.FindIndex(v => v.IndexInList == vIndex2); + // var step = v2 - v1; + // if (step < 0) step += 3; + // if (step == 1) + // OuterEdges.Add(new Edge(face.A, face.B, face, null, false)); + // else OuterEdges.Add(new Edge(face.B, face.A, face, null, false)); + // } + //} Faces.Add(face); face.BelongsToPrimitive = this; } @@ -948,10 +950,15 @@ public PrimitiveSurface Copy(IEnumerable faces, bool doNotResetFac /// Copies the primitive surface including the faces and vertices (edges are not defined in the copy) /// /// - public PrimitiveSurface Copy() + public PrimitiveSurface Copy(bool copyTessellationElements = true) { var copy = (PrimitiveSurface)this.Clone(); - if (Faces != null && Faces.Count != 0) + copy.Faces = new HashSet(); + copy.Vertices = new HashSet(); + copy.InnerEdges = new HashSet(); + copy.OuterEdges = new HashSet(); + copy.IsPositive = IsPositive; + if (copyTessellationElements && Faces != null && Faces.Count != 0) { var i = 0; //NumberOfVertices = vertices.Count; @@ -975,16 +982,9 @@ public PrimitiveSurface Copy() i++; } copy.SetFacesAndVertices(newFaces, true); - if (Faces.SelectMany(f => f.Edges).Any()) - { - var tempSolid = new TessellatedSolid(copy.Faces, copy.Vertices, - TessellatedSolidBuildOptions.Minimal); - //Edges are defined, so we need to makes these. This is like the MakeEdgesIfNonExistent - var repair = new TessellationInspectAndRepair(tempSolid); - repair.CheckModelIntegrityPreBuild(); - repair.MakeEdges(); - } } + copy.ResetFaceDependentValues(); + return copy; } @@ -1007,15 +1007,5 @@ private protected string GetCommonKeyDetails() return key; } - - public Polygon GetAsPolygon() - { - var polygons = new List(); - var vertexLoops = OuterEdges.MakeEdgePaths(true, - new EdgePathLoopsAroundInputFaces(Faces)).Select(ep => ep.GetVertices().ToList()); - foreach (var loop in vertexLoops) - polygons.Add(new Polygon(TransformFrom3DTo2D(loop.Select(v => v.Coordinates), true))); - return polygons.CreateShallowPolygonTrees(false).First(); - } } } \ No newline at end of file diff --git a/TessellationAndVoxelizationGeometryLibrary/Implicits and Primitives/Primitive Surfaces/PrimitiveSurfaceExtensions.cs b/TessellationAndVoxelizationGeometryLibrary/Implicits and Primitives/Primitive Surfaces/PrimitiveSurfaceExtensions.cs index 4f75c25de..c661f0149 100644 --- a/TessellationAndVoxelizationGeometryLibrary/Implicits and Primitives/Primitive Surfaces/PrimitiveSurfaceExtensions.cs +++ b/TessellationAndVoxelizationGeometryLibrary/Implicits and Primitives/Primitive Surfaces/PrimitiveSurfaceExtensions.cs @@ -198,7 +198,7 @@ public static Vector3 GetAnchor(this PrimitiveSurface s) else if (s is Torus torus) return torus.Center; else if (s is Sphere sphere) return sphere.Center; else if (s is Capsule capsule) return 0.5 * (capsule.Anchor1 + capsule.Anchor2); - else if (s is GeneralQuadric gq) return gq.Center; + else if (s is GeneralQuadric gq) return gq.StationaryPoint; else return GetCenterOfMass(s.Faces); } @@ -410,7 +410,44 @@ private static bool PointASeesPointB(Vector2 pointA, Vector2 vA1, Vector2 vA2, V /// Thrown if an unexpected case occurs where two edges are removed from a triangle, which indicates an invalid /// tessellation state. public static void TrimTessellationToPositiveVertices(this PrimitiveSurface surface, - IDictionary vertexDistances, IDictionary edgeIntersections = null) + IDictionary vertexDistances) + { + var edgeIntersections = new Dictionary(); + + foreach (var edge in surface.InnerEdges.Concat(surface.OuterEdges)) + { + var fromDist = vertexDistances[edge.From]; + var toDist = vertexDistances[edge.To]; + if (fromDist < 0 && toDist < 0) + continue; // both negative, so skip + if (fromDist >= 0 && toDist >= 0) + continue; // both non-negative, so skip + // edge is crossing, so add the intersection to the dictionary for use in the main function below + edgeIntersections.Add(edge, edge.Length * (fromDist / (fromDist - toDist))); + } + var removeVertices = vertexDistances.Where(kvp => kvp.Value < 0).Select(kvp => kvp.Key).ToHashSet(); + TrimTessellationToPositiveVertices(surface, edgeIntersections, removeVertices); + } + + + /// + /// Trims the tessellation of the surface by removing faces and edges associated with vertices that have + /// negative distance values, retaining only those with positive or zero distances. New vertices, edges, + /// and faces are made at the boundarys where edges cross from negative to non-negative distances. + /// + /// This method modifies the input surface in place. After execution, only the portions + /// of the surface associated with non-negative vertex distances remain. The method is typically used to clip or + /// trim a mesh based on a scalar field or distance function. + /// The primitive surface whose tessellation will be modified to exclude faces and edges with negative vertex + /// distances. + /// A mapping of edges to their corresponding intersection values. This is used to + /// determine where new vertices should be created along edges that cross from negative to non-negative distances. Distances are + /// not normalized (vary from 0 to edge length) and start from the "from" vertex. + /// A set of vertices that have negative distance values and should be removed from the tessellation. + /// Thrown if an unexpected case occurs where two edges are removed from a triangle, which indicates an invalid + /// tessellation state. + public static void TrimTessellationToPositiveVertices(this PrimitiveSurface surface, + IDictionary edgeIntersections, HashSet removeVertices) { var edgesToRemove = new HashSet(); // we keep track of edges to remove to remove their // reference from the vertices later @@ -420,108 +457,123 @@ public static void TrimTessellationToPositiveVertices(this PrimitiveSurface surf // most are replaced by facesToAdd below (unless they are all negative) // we keep track of these to remove their references from the vertices later var facesToAdd = new List(); // added to the primitive surface with the existing in SetFacesAndVertices + + + // var fGroups = new List[4]; + // var fcolors = new Color[]{ new Color(KnownColors.LightGreen),new Color(KnownColors.Yellow), + // new Color(KnownColors.Orange),new Color(KnownColors.Red)}; + // for (int i = 0; i < 4; i++) + // fGroups[i] = new List(); + // foreach (var f in surface.Faces) + // { + // var count = f.Vertices.Sum(v => removeVertices.Contains(v) ? 1 : 0); + // fGroups[count].Add(f); + // } + // for (int i = 0; i < 4; i++) + // foreach (var f in fGroups[i]) + // f.Color = fcolors[i]; + // //Presenter.ShowAndHang(surface.Faces); + //Presenter.ShowAndHang( //[surface.InnerEdges.Select(e=>new[] { e.From.Coordinates, e.To.Coordinates } ), + // edgeIntersections.Keys.Select(e=>new[] { e.From.Coordinates, e.To.Coordinates } ),lineThicknesses: [3], + // colors: [new Color(KnownColors.Black)], + // faces:surface.Faces + // ); foreach (var face in surface.Faces) { - if (!vertexDistances.TryGetValue(face.A, out var aDist)) - aDist = double.MinValue; - if (!vertexDistances.TryGetValue(face.B, out var bDist)) - bDist = double.MinValue; - if (!vertexDistances.TryGetValue(face.C, out var cDist)) - cDist = double.MinValue; - if (aDist >= 0 && bDist >= 0 && cDist >= 0) - continue; // all positive, so keep face as is + var aRemoved = removeVertices.Contains(face.A); + var bRemoved = removeVertices.Contains(face.B); + var cRemoved = removeVertices.Contains(face.C); + if (!aRemoved && !bRemoved && !cRemoved) + continue; // all keep vertices, so skip to next facesToRemove.Add(face); - var thisPatchRemovedEdges = new List(); - if (aDist < 0) + if (aRemoved && bRemoved && cRemoved) + continue; // all remove vertices, so - yes- remove face (as above, but then leave) + //face.Color = new Color(KnownColors.Red); + //Presenter.ShowAndHang(surface.Faces.Take(40).Concat(new[] { face })); + Edge thisPatchRemovedEdge = null; + if (aRemoved) { - if (bDist < 0) // then edge AB is below and is to be removed + if (bRemoved) // then edge AB is below (both vertices are to be removed), so it is to be removed { edgesToRemove.Add(face.AB); - thisPatchRemovedEdges.Add(face.AB); + thisPatchRemovedEdge = face.AB; } else { // edge AB is crossing, and B is above, so replace if (alreadyModifiedEdges.Add(face.AB)) - ReplaceNegativeVertexOnEdge(face.AB, face.A, face.B, aDist, bDist, edgeIntersections); + ReplaceNegativeVertexOnEdge(face.AB, face.A, edgeIntersections); } - if (cDist < 0) + if (cRemoved) { edgesToRemove.Add(face.CA); - thisPatchRemovedEdges.Add(face.CA); + thisPatchRemovedEdge = face.CA; } else { // edge CA is crossing, and C is above, so replace if (alreadyModifiedEdges.Add(face.CA)) - ReplaceNegativeVertexOnEdge(face.CA, face.A, face.C, aDist, cDist, edgeIntersections); + ReplaceNegativeVertexOnEdge(face.CA, face.A, edgeIntersections); } } - if (bDist < 0) + if (bRemoved) { // edge AB already handled above if both are negative, just need to check // if A is positive and the edge is to be modified - if (aDist >= 0 && alreadyModifiedEdges.Add(face.AB)) - ReplaceNegativeVertexOnEdge(face.AB, face.B, face.A, bDist, aDist, edgeIntersections); - if (cDist < 0) + if (!aRemoved && alreadyModifiedEdges.Add(face.AB)) + ReplaceNegativeVertexOnEdge(face.AB, face.B, edgeIntersections); + if (cRemoved) { edgesToRemove.Add(face.BC); - thisPatchRemovedEdges.Add(face.BC); + thisPatchRemovedEdge = face.BC; } else { // edge BC is crossing, and C is above, so replace if (alreadyModifiedEdges.Add(face.BC)) - ReplaceNegativeVertexOnEdge(face.BC, face.B, face.C, bDist, cDist, edgeIntersections); + ReplaceNegativeVertexOnEdge(face.BC, face.B, edgeIntersections); } } - if (cDist < 0) + if (cRemoved) { // now the three edge removals would have happened in the above // conditions, so here just need to check if A or B are positive // and the edge is to be modified - if (aDist >= 0 && alreadyModifiedEdges.Add(face.CA)) - ReplaceNegativeVertexOnEdge(face.CA, face.C, face.A, cDist, aDist, edgeIntersections); - if (bDist >= 0 && alreadyModifiedEdges.Add(face.BC)) - ReplaceNegativeVertexOnEdge(face.BC, face.C, face.B, cDist, bDist, edgeIntersections); + if (!aRemoved && alreadyModifiedEdges.Add(face.CA)) + ReplaceNegativeVertexOnEdge(face.CA, face.C, edgeIntersections); + if (!bRemoved && alreadyModifiedEdges.Add(face.BC)) + ReplaceNegativeVertexOnEdge(face.BC, face.C, edgeIntersections); + //face.Color = new Color(KnownColors.DarkOliveGreen); + //Presenter.ShowAndHang(surface.Faces); } - if (thisPatchRemovedEdges.Count == 3) - continue; // entire face is removed - if (thisPatchRemovedEdges.Count == 2) - throw new Exception("This case should not happen - two edges removed from a triangle implies all three vertices are negative."); - if (thisPatchRemovedEdges.Count == 1) - { // one edge removed, so create one new face - var orderedFaceVerts = GetVerticesFromModifyEdges(face, thisPatchRemovedEdges[0], out var oldVertIndex); - var newFace = new TriangleFace(orderedFaceVerts.Take(3), true); - facesToAdd.Add(newFace); - if (oldVertIndex == 1) - newFace.AddEdge(new Edge(newFace.C, newFace.A, newFace, null, true)); - else - newFace.AddEdge(new Edge(newFace.B, newFace.C, newFace, null, true)); - if (face.AB != thisPatchRemovedEdges[0]) - newFace.AddEdge(face.AB); - if (face.BC != thisPatchRemovedEdges[0]) - newFace.AddEdge(face.BC); - if (face.CA != thisPatchRemovedEdges[0]) - newFace.AddEdge(face.CA); - } - else //if (thisPatchRemovedEdges.Count==0) + if (thisPatchRemovedEdge == null) { - // no edges removed, so must be two crossing edges, so create two new faces - // find the one edge that is kept + // no edges removed means that two new faces must be created find the one edge that is kept var orderedFaceVerts = GetVerticesFromModifyEdges(face, null, out var oldVertIndex); var newFace = new TriangleFace(orderedFaceVerts.Take(3), true); - newFace.AddEdge(face.AB); + facesToAdd.Add(newFace); + newFace.AddEdge(face.AB); // how do you know that face.AB is part of newFace? because orderedFaceVerts + // would have started with that + if (face.AB.OwnedFace == face) + face.AB.OwnedFace = newFace; + else face.AB.OtherFace = newFace; var innerEdge = new Edge(newFace.C, newFace.A, newFace, null, true); newFace.AddEdge(innerEdge); if (oldVertIndex == 0) { - newFace.AddEdge(new Edge(newFace.B, newFace.C, newFace, null, true)); + newFace.AddEdge(new Edge(newFace.B, newFace.C, newFace, null, true)); // this is the border edge + // so, no need to do any more here newFace = new TriangleFace(orderedFaceVerts.Skip(2), true); facesToAdd.Add(newFace); newFace.AddEdge(face.BC); + if (face.BC.OwnedFace == face) + face.BC.OwnedFace = newFace; + else face.BC.OtherFace = newFace; } else if (oldVertIndex == 1) { newFace.AddEdge(face.BC); + if (face.BC.OwnedFace == face) + face.BC.OwnedFace = newFace; + else face.BC.OtherFace = newFace; newFace = new TriangleFace(orderedFaceVerts.Skip(2), true); facesToAdd.Add(newFace); newFace.AddEdge(new Edge(newFace.A, newFace.B, newFace, null, true)); @@ -529,26 +581,57 @@ public static void TrimTessellationToPositiveVertices(this PrimitiveSurface surf else // oldVertIndex ==2 { newFace.AddEdge(face.BC); + if (face.BC.OwnedFace == face) + face.BC.OwnedFace = newFace; + else face.BC.OtherFace = newFace; newFace = new TriangleFace(orderedFaceVerts.Skip(2).Concat([orderedFaceVerts[0]]), true); facesToAdd.Add(newFace); newFace.AddEdge(new Edge(newFace.B, newFace.C, newFace, null, true)); } newFace.AddEdge(face.CA); + if (face.CA.OwnedFace == face) + face.CA.OwnedFace = newFace; + else face.CA.OtherFace = newFace; newFace.AddEdge(innerEdge); + innerEdge.OtherFace = newFace; + } + else + { // one edge removed, so create one new face + var orderedFaceVerts = GetVerticesFromModifyEdges(face, thisPatchRemovedEdge, out var keptVertex); + var newFace = new TriangleFace(orderedFaceVerts.Take(3), true); + facesToAdd.Add(newFace); + Edge newBorderEdge; + if (keptVertex == 1) // then vertex B (A is 0, C is 2) is the one that is kept, so edges AB and BC are + // clipped above but not delete. Therefore, we need a new CA edge + newBorderEdge = new Edge(newFace.C, newFace.A, newFace, null, true); + else // then keptvertex A or C so edges BC and CA are clipped above but not delete. But, (this is really + // confusing) for the new face, it is alway the second ede (BC) that needs to be recreated + newBorderEdge = new Edge(newFace.B, newFace.C, newFace, null, true); + newFace.AddEdge(newBorderEdge); + foreach (var e in face.Edges) + { + if (e != thisPatchRemovedEdge) + { + newFace.AddEdge(e); + if (e.OwnedFace == face) + e.OwnedFace = newFace; + else e.OtherFace = newFace; + } + } } - foreach (var f in facesToAdd) - if (f.BC == null) ; } // edgesToRemove ,facesToRemove,facesToAdd surface.Faces.ExceptWith(facesToRemove); surface.SetFacesAndVertices(surface.Faces.Concat(facesToAdd), true); surface.DefineInnerOuterEdges(); foreach (var removeFace in facesToRemove) - foreach (var v in surface.Vertices) + foreach (var v in removeFace.Vertices) v.Faces.Remove(removeFace); foreach (var removeEdge in edgesToRemove) - foreach (var v in surface.Vertices) - v.Edges.Remove(removeEdge); + { + removeEdge.To.Edges.Remove(removeEdge); + removeEdge.From.Edges.Remove(removeEdge); + } ReindexVertices(surface); } @@ -566,6 +649,14 @@ public static void ReindexVertices(this PrimitiveSurface surface) v.IndexInList = i++; } + /// + /// + /// + /// + /// + /// the highest vertex index that was not cut. If removedEdge is not null, then + /// this is the only vertex index that is not cut. + /// private static List GetVerticesFromModifyEdges(TriangleFace face, Edge removedEdge, out int oldVertIndex) { var result = new List(); @@ -591,14 +682,10 @@ private static List GetVerticesFromModifyEdges(TriangleFace face, Edge r return result; } - private static Vertex ReplaceNegativeVertexOnEdge(Edge edge, Vertex negV, Vertex posV, double negDist, double posDist, - IDictionary edgeIntersections) + + private static Vertex ReplaceNegativeVertexOnEdge(Edge edge, Vertex negV, IDictionary edgeIntersections) { - Vector3 position; - if (edgeIntersections == null) - position = posV.Coordinates + (negV.Coordinates - posV.Coordinates) * (posDist / (posDist - negDist)); - else - position = edge.From.Coordinates + edge.Vector * edgeIntersections[edge]; + var position = edge.From.Coordinates + edge.UnitVector * edgeIntersections[edge]; var newVertex = new Vertex(position); if (edge.To == negV) edge.To = newVertex; else edge.From = newVertex; @@ -1120,6 +1207,18 @@ public static TessellatedSolid Tessellate(this Torus torus, int numPointsInMajor #endregion + + + public static Polygon GetAs2DPolygon(this PrimitiveSurface surface) + { + var polygons = new List(); + var vertexLoops = surface.OuterEdges.MakeEdgePaths(true, + new EdgePathLoopsAroundInputFaces(surface.Faces)).Select(ep => ep.GetVertices().ToList()); + foreach (var loop in vertexLoops) + polygons.Add(new Polygon(surface.TransformFrom3DTo2D(loop.Select(v => v.Coordinates), true))); + return polygons.CreateShallowPolygonTrees(false).First(); + } + private static IEnumerable<(Vector3 intersection, double lineT)> GetPrimitiveAndLineIntersections(PrimitiveSurface surface, double xCoord, double yCoord, double zCoord) { diff --git a/TessellationAndVoxelizationGeometryLibrary/Implicits and Primitives/TwoDimensional/GeneralConicSection.cs b/TessellationAndVoxelizationGeometryLibrary/Implicits and Primitives/TwoDimensional/GeneralConicSection.cs index 8ab59d8f0..8bd76f3d1 100644 --- a/TessellationAndVoxelizationGeometryLibrary/Implicits and Primitives/TwoDimensional/GeneralConicSection.cs +++ b/TessellationAndVoxelizationGeometryLibrary/Implicits and Primitives/TwoDimensional/GeneralConicSection.cs @@ -50,27 +50,28 @@ public enum PrimitiveCurveType public struct GeneralConicSection : ICurve { /// - /// a + /// the coefficient multiplying x^2 /// public double A; /// - /// The b + /// the coefficient multiplying xy /// public double B; /// - /// The c + /// the coefficient multiplying y^2 /// public double C; /// /// The d + /// the coefficient multiplying x /// public double D; /// - /// The e + /// the coefficient multiplying y /// public double E; /// - /// The constant is zero + /// The constant is zero, otherwise the constant is 1 /// public bool ConstantIsZero; /// @@ -85,26 +86,68 @@ public struct GeneralConicSection : ICurve public PrimitiveCurveType CurveType { get; private set; } /// - /// Initializes a new instance of the struct. + /// Initializes a new instance of the GeneralConicSection class using the specified coefficients of the general + /// conic equation. /// - /// a. - /// The b. - /// The c. - /// The d. - /// The e. + /// The general conic equation is typically represented as: a·x² + b·xy + c·y² + d·x + + /// e·y + w = 0. If the constant term w is negligible, the coefficients are used as provided; otherwise, all + /// coefficients are normalized by dividing by w. + /// The coefficient of the x² term in the conic equation. + /// The coefficient of the xy term in the conic equation. + /// The coefficient of the y² term in the conic equation. + /// The coefficient of the x term in the conic equation. + /// The coefficient of the y term in the conic equation. /// if set to true [constant is zero]. public GeneralConicSection(double a, double b, double c, double d, double e, bool constantIsZero) { - var max = (new[] { Math.Abs(a), Math.Abs(b), Math.Abs(c), Math.Abs(d), Math.Abs(e) }).Max(); - A = Math.Abs(a / max) < conicTolerance ? 0 : a; - B = Math.Abs(b / max) < conicTolerance ? 0 : b; - C = Math.Abs(c / max) < conicTolerance ? 0 : c; - D = Math.Abs(d / max) < conicTolerance ? 0 : d; - E = Math.Abs(e / max) < conicTolerance ? 0 : e; + A = a; + B = b; + C = c; + D = d; + E = e; ConstantIsZero = constantIsZero; - CurveType = PrimitiveCurveType.StraightLine; + SetConicType(); + } + + /// + /// Initializes a new instance of the GeneralConicSection class using the specified coefficients of the general + /// conic equation. + /// + /// The general conic equation is typically represented as: a·x² + b·xy + c·y² + d·x + + /// e·y + w = 0. If the constant term w is negligible, the coefficients are used as provided; otherwise, all + /// coefficients are normalized by dividing by w. + /// The coefficient of the x² term in the conic equation. + /// The coefficient of the xy term in the conic equation. + /// The coefficient of the y² term in the conic equation. + /// The coefficient of the x term in the conic equation. + /// The coefficient of the y term in the conic equation. + /// The constant term in the conic equation. If negligible, the conic is considered to have a zero constant + /// term. + public GeneralConicSection(double a, double b, double c, double d, double e, double w) : this() + { + if (w.IsNegligible()) + { + A = a; + B = b; + C = c; + D = d; + E = e; + ConstantIsZero = true; + } + else + { + A = a / w; + B = b / w; + C = c / w; + D = d / w; + E = e / w; + ConstantIsZero = false; + } + SetConicType(); + } - #region SetConicType + private void SetConicType() + { if (A.IsNegligible() && B.IsNegligible() && C.IsNegligible()) { A = B = C = 0; @@ -127,9 +170,9 @@ public GeneralConicSection(double a, double b, double c, double d, double e, boo if (det > 0) CurveType = PrimitiveCurveType.Ellipse; else CurveType = PrimitiveCurveType.Hyperbola; } - #endregion } + /// /// Calculates at point. /// @@ -139,7 +182,8 @@ public double CalculateAtPoint(Vector2 point) { double x = point.X; double y = point.Y; - return A * x * x + B * x * y + C * y * y + D * x + E * y + 1; + var constant = ConstantIsZero ? 0.0 : 1.0; + return A * x * x + B * x * y + C * y * y + D * x + E * y + constant; } @@ -273,10 +317,7 @@ public static double DistancePointToConic(GeneralConicSection conic, T point, var rj = conic.E - conic.B * point[0] + 2 * conic.A * point[1]; var sj = -conic.D + conic.B * point[1] - 2 * conic.C * point[0]; var tj = conic.D * point[1] - conic.E * point[0]; - GeneralConicSection conicJ; - if (tj.IsNegligible()) - conicJ = new GeneralConicSection(aj, bj, cj, rj, sj, true); - else conicJ = new GeneralConicSection(aj / tj, bj / tj, cj / tj, rj / tj, sj / tj, false); + var conicJ = new GeneralConicSection(aj, bj, cj, rj, sj, tj); var minDistance = double.PositiveInfinity; pointOnCurve = Vector2.Null; foreach (var p in IntersectingConics(conic, conicJ)) @@ -290,7 +331,114 @@ public static double DistancePointToConic(GeneralConicSection conic, T point, } return minDistance; } + public static GeneralConicSection CreateFromQuadric(GeneralQuadric quadric, Plane plane) + { + var mTranspose = plane.AsTransformFromXYPlane.Transpose(); + var qMatrix = quadric.GetCoefficientMatrix(); + // In TVGL's row-vector convention, a 2D point p maps to 3D as x = p * M. + // Substituting into the quadric x * Q * x^T = 0 gives p * (M * Q * M^T) * p^T = 0. + var qNew = plane.AsTransformFromXYPlane * qMatrix * mTranspose; + var a = qNew.M11; + var b = 2 * qNew.M12; + var c = qNew.M22; + var d = 2 * qNew.M14; + var e = 2 * qNew.M24; + var w = qNew.M44; + + return new GeneralConicSection(a, b, c, d, e, w); + } + + /// + /// Finds the points on the conic that has a gradient in the same direction as + /// the specified gradient vector. This input does not have to be normalized. + /// + /// + /// + /// + public bool PointsAtGivenGradient(Vector2 gradient, out Vector2 point) + { + point = Vector2.Null; + // The gradient of f(x,y) = Ax² + Bxy + Cy² + Dx + Ey + const is: + // ∇f = [2Ax + By + D, Bx + 2Cy + E] + // For ∇f ∥ gradient, the 2D cross product ∇f × gradient = 0: + // (2Ax + By + D)*gy - (Bx + 2Cy + E)*gx = 0 + // This defines a line: alpha*x + beta*y + gamma = 0 + // Intersecting this line with the conic finds the desired points. + // Unlike inverting the gradient matrix, this works for all conic types + // including parabolas where the matrix is singular. + var gx = gradient.X; + var gy = gradient.Y; + var alpha = 2 * A * gy - B * gx; + var beta = B * gy - 2 * C * gx; + var gamma = D * gy - E * gx; + + Vector2 anchor, dir; + if (Math.Abs(alpha) >= Math.Abs(beta)) + { + if (alpha.IsNegligible()) return false; + anchor = new Vector2(-gamma / alpha, 0); + dir = new Vector2(-beta / alpha, 1); + } + else + { + anchor = new Vector2(0, -gamma / beta); + dir = new Vector2(1, -alpha / beta); + } + + foreach (var tuple in LineIntersection(anchor, dir)) + { + var intersection = tuple.intersection; + if (GetGradient(intersection).Dot(gradient) > 0) + { + point = intersection; + return true; + } + } + return false; + } + + private Vector2 GetGradient(Vector2 pt) + { + return new Vector2(2 * A * pt.X + B * pt.Y + D, 2 * C * pt.Y + B * pt.X + E); + } + + + /// + /// Returns the intersection points between this quadric and the given line. + /// + /// + /// + /// + public IEnumerable<(Vector2 intersection, double lineT)> LineIntersection(Vector2 anchor, Vector2 direction) + { + //put equation for line p = anchor + t * direction into the conic equation, + //we get a quadratic equation in terms of t. Solve it and we get the intersection points. + // so, here we just find the at^2 + bt + c = 0, where a, b, c are calculated as below + // x = anchor.X + t * direction.X + // y = anchor.Y + t * direction.Y + // A * x^2 + B * x * y + C * y^2 + D * x + E * y + constant = 0 + // collect t^2, t, and constant terms, we get + var constant = ConstantIsZero ? 0.0 : 1.0; + var a = A * direction.X * direction.X + B * direction.X * direction.Y + C * direction.Y * direction.Y; + var b = 2 * A * anchor.X * direction.X + B * (anchor.X * direction.Y + anchor.Y * direction.X) + 2 * C * anchor.Y * direction.Y + + D * direction.X + E * direction.Y; + var c = A * anchor.X * anchor.X + B * anchor.X * anchor.Y + C * anchor.Y * anchor.Y + D * anchor.X + E * anchor.Y + constant; + (var root1, var root2) = PolynomialSolve.Quadratic(a, b, c); + + if (!root1.IsRealNumber) + yield break; + if (root1.Real.IsPracticallySame(root2.Real)) + { + var t = 0.5 * (root1.Real + root2.Real); + yield return (anchor + t * direction, root1.Real); + yield break; + } + if (root1.IsRealNumber) + yield return (anchor + root1.Real * direction, root1.Real); + if (root2.IsRealNumber) + yield return (anchor + root2.Real * direction, root2.Real); + } /// /// Finds the zero to four points of two intersectings the conics. diff --git a/TessellationAndVoxelizationGeometryLibrary/MarchingCubes/MarchingCubes.Implicit.cs b/TessellationAndVoxelizationGeometryLibrary/MarchingCubes/MarchingCubes.Implicit.cs index f626627d1..d34162361 100644 --- a/TessellationAndVoxelizationGeometryLibrary/MarchingCubes/MarchingCubes.Implicit.cs +++ b/TessellationAndVoxelizationGeometryLibrary/MarchingCubes/MarchingCubes.Implicit.cs @@ -180,12 +180,20 @@ protected override void MakeFacesInCube(int xIndex, int yIndex, int zIndex) var faceVertices = new Vertex[3]; for (var i = 0; i < NumFacesTable[cubeType]; i++) { + var validFace = true; for (var j = 0; j < 3; j++) { var vertexIndex = FaceVertexIndicesTable[cubeType][3 * i + j]; - faceVertices[j] = EdgeVertex[vertexIndex]; + var vertex = EdgeVertex[vertexIndex]; + if (vertex.Coordinates.IsNull()) + { + validFace = false; + break; + } + faceVertices[j] = vertex; } - faces.Add(new TriangleFace(faceVertices)); + if (validFace) + faces.Add(new TriangleFace(faceVertices)); } } diff --git a/TessellationAndVoxelizationGeometryLibrary/Miscellaneous Functions/MiscFunctions.cs b/TessellationAndVoxelizationGeometryLibrary/Miscellaneous Functions/MiscFunctions.cs index a741d1a00..b7bc9fd85 100644 --- a/TessellationAndVoxelizationGeometryLibrary/Miscellaneous Functions/MiscFunctions.cs +++ b/TessellationAndVoxelizationGeometryLibrary/Miscellaneous Functions/MiscFunctions.cs @@ -254,7 +254,7 @@ public static (T minPoint, double dotDistance, int index) GetMinVertexDistanceAl } } return (minPoint, dotDistance, minPointIndex); - } + } /// /// Finds the vertices or points that are closest and furthest along a given direction vector in a single pass. @@ -432,6 +432,8 @@ private static (Vector2, double)[] GetPointDistances(this IList points, } return distances; } + #endregion Sort Along Direction + /// /// From a collection of faces, this method identifies and separates the edges into two sets: "inner" edges and "outer" edges. @@ -449,22 +451,112 @@ public static void DefineInnerOuterEdges(IEnumerable faces, out Ha { innerEdgeHash = new HashSet(); outerEdgeHash = new HashSet(); - if (faces != null && faces.Any()) + if (faces == null || !faces.Any()) + return; + if (faces.First().Edges.FirstOrDefault() == null) + DefineInnerOuterEdgesWithNoEdges(faces, innerEdgeHash, outerEdgeHash); + else + { foreach (var face in faces) { foreach (var edge in face.Edges) { if (innerEdgeHash.Contains(edge)) continue; - if (!outerEdgeHash.Contains(edge)) outerEdgeHash.Add(edge); - else + if (!outerEdgeHash.Add(edge)) { innerEdgeHash.Add(edge); outerEdgeHash.Remove(edge); } } } + } + } + + /// + /// From a collection of faces, this method identifies and separates the edges into two sets: "inner" edges and "outer" edges. + /// + /// The collection of triangle faces. + /// An output hash set containing the inner edges. An edge is considered "inner" if it is shared by two faces in the collection. + /// An output hash set containing the outer edges. An edge is considered "outer" or a "border" edge if it belongs to only one face in the collection. + /// + /// This method is useful for identifying the boundary of a patch of faces on a larger mesh, or for manifold checking. + /// It works by iterating through all edges of the provided faces. It uses a hash set to keep track of edges it has already seen. + /// If an edge is seen a second time, it's moved to the inner edge set. + /// Common search terms: "find boundary edges", "identify border edges", "inner vs outer edges". + /// + private static void DefineInnerOuterEdgesWithNoEdges(IEnumerable faces, HashSet innerEdgeHash, + HashSet outerEdgeHash) + { + var allVertices = new HashSet(); + var vertexIndices = new HashSet(); + var maxVIndex = int.MinValue; + var needToReIndex = false; + foreach (var f in faces) + { + foreach (var v in f.Vertices) + { + if (allVertices.Add(v)) //first time seeing this vertex + { + if (v.IndexInList < 0 || !vertexIndices.Add(v.IndexInList)) + needToReIndex = true; + else if (maxVIndex < v.IndexInList) + maxVIndex = v.IndexInList; + } + } + } + if (needToReIndex) + { + var k = 0; + foreach (var v in allVertices) + v.IndexInList = k++; + maxVIndex = k; + } + else maxVIndex++; + var edges = new Dictionary(); + var allEdgeIndices = new List(); + var doubleEdgeIndices = new HashSet(); + foreach (var f in faces) + { + var prevV = f.C; + var edge = f.CA; + var edgeEnumerator = f.Edges.GetEnumerator(); + foreach (var v in f.Vertices) + { + long edgeIndex = (v.IndexInList < prevV.IndexInList) + ? (long)prevV.IndexInList * (long)maxVIndex + v.IndexInList + : prevV.IndexInList + (long)maxVIndex * (long)v.IndexInList; + allEdgeIndices.Add(edgeIndex); + if (edge == null) + { + if (edges.TryGetValue(edgeIndex, out edge)) + { + doubleEdgeIndices.Add(edgeIndex); + edge.OtherFace = f; + } + else + edges[edgeIndex] = edge = new Edge(prevV, v, f, null, true, edgeIndex); + f.AddEdge(edge); + } + else + { + if (edges.ContainsKey(edgeIndex)) doubleEdgeIndices.Add(edgeIndex); + else edges.Add(edgeIndex, edge); + } + prevV = v; + edgeEnumerator.MoveNext(); + edge = edgeEnumerator.Current; + } + } + var ei = 0; + foreach (var i in allEdgeIndices) + { + var edge = edges[i]; + edge.IndexInList = ei++; + if (doubleEdgeIndices.Contains(i)) + innerEdgeHash.Add(edges[i]); + else outerEdgeHash.Add(edges[i]); + } } - #endregion Sort Along Direction #region Perimeter @@ -956,7 +1048,7 @@ public static IEnumerable GetMultipleSolids(this TessellatedSo unusedFaces = new HashSet(ts.Faces); } // now, the bigger job of walking through the faces to find groups - faceGroups.AddRange(GetContiguousFaceGroups(unusedFaces)); + faceGroups.AddRange(GetContiguousFaceGroups(unusedFaces).Select(fg => fg.ToList())); if (faceGroups.Count == 1 && faceGroups[0].Count == ts.NumberOfFaces) { @@ -1011,32 +1103,64 @@ public static IEnumerable GetMultipleSolids(this TessellatedSo } } - public static List> GetContiguousFaceGroups(this IEnumerable facesInput) - { - var unusedFaces = facesInput as HashSet ?? facesInput.ToHashSet(); - var faceGroups = new List>(); - while (unusedFaces.Any()) + + /// + /// Identifies and returns contiguous patches of triangle faces, grouping faces that are connected without + /// crossing specified stop edges or stop faces. + /// + /// This method is useful for segmenting a mesh into regions based on connectivity, with + /// the ability to specify custom boundaries. The returned patches are disjoint, and each face from the input is + /// included in at most one patch unless it is specified as a stop face. + /// The collection of triangle faces to be partitioned into patches. + /// An optional set of edges that act as boundaries; patches will not cross these edges. If null, no edge + /// boundaries are applied. + /// An optional set of faces that act as boundaries; patches will not include or cross these faces. If null, no + /// face boundaries are applied. + /// An enumerable collection of face patches, where each patch is a set of connected triangle faces. Each patch + /// contains faces that are contiguous and do not cross any stop edges or stop faces. + public static IEnumerable> GetContiguousFaceGroups(this IEnumerable faces, + HashSet stopEdges = null, HashSet stopFaces = null) + { + var remainingFaces = new HashSet(faces); + if (stopEdges != null) + stopEdges = new HashSet(stopEdges, stopEdges.Comparer); + else stopEdges = new HashSet(); + if (stopFaces != null) + stopFaces = new HashSet(stopFaces, stopFaces.Comparer); + else stopFaces = new HashSet(); + //Pick a start edge, then collect all adjacent faces on one side of the face, without crossing over significant edges. + //This collection of faces will be used to create a patch. + while (remainingFaces.Any()) { - var groupHash = new HashSet(); - var stack = new Stack(new[] { unusedFaces.First() }); + var startFace = remainingFaces.First(); + remainingFaces.Remove(startFace); + if (stopFaces.Contains(startFace)) continue; // this is redundant with the below check for the same, but + // check here says a split second and the creation of the next two collections. + var patch = new HashSet(); + var stack = new Stack(); + stack.Push(startFace); while (stack.Any()) { var face = stack.Pop(); - if (groupHash.Contains(face)) continue; - groupHash.Add(face); - unusedFaces.Remove(face); - foreach (var adjacentFace in face.AdjacentFaces) + if (stopFaces.Contains(face)) continue; + stopFaces.Add(face); + patch.Add(face); + foreach (var edge in face.Edges) { - if (adjacentFace == null) continue; //This is an error. Handle it in the error function. - if (!unusedFaces.Contains(adjacentFace)) continue; //Cannot assign the same face twice - stack.Push(adjacentFace); + if (stopEdges.Contains(edge)) continue;//Don't cross over significant edges + var otherFace = face == edge.OwnedFace ? edge.OtherFace : edge.OwnedFace; + if (remainingFaces.Contains(otherFace)) + { + stack.Push(otherFace); + remainingFaces.Remove(otherFace); + } } } - faceGroups.Add(groupHash.ToList()); + yield return patch; } - return faceGroups; } + public static List> GetContiguousFaceGroups(TessellatedSolid ts, List faceGroupsThatAreBodies, out HashSet unusedFaces) { @@ -1363,6 +1487,63 @@ public static Matrix4x4 TransformToXYPlane(this Vector3 direction, out Matrix4x4 return backTransform.Transpose(); } } + /// + /// Creates a transformation matrix that rotates a given 3D direction vector to align with the Z-axis (0, 0, 1), effectively creating a 2D projection plane. + /// + /// The 3D vector to be aligned with the Z-axis. This vector represents the normal of the desired 2D plane. + /// An output parameter for the inverse transformation matrix, which can rotate the Z-axis back to the original `direction`. + /// A tolerance used to determine if the input direction is close enough to a Cartesian axis to use a simpler, predefined transformation. + /// A 4x4 transformation matrix that projects points onto the XY plane. + /// + /// This is a core function for converting 3D geometry into a 2D representation for slicing, analysis, or visualization. The `backTransform` is essential for converting results from 2D calculations back into the original 3D coordinate system. + /// The method includes an optimization to "snap" to major Cartesian axes, which is faster and avoids potential floating-point precision issues. + /// Common search terms: "view matrix", "projection matrix", "align vector to axis", "lookat transform". + /// + public static Matrix4x4 TransformToXYPlaneMaybeBetter(this Vector3 direction, out Matrix4x4 backTransform, double tolerance = Constants.DefaultEqualityTolerance) + { + var closestCartesianDirection = SnapDirectionToCartesian(direction, out var withinTolerance, tolerance); + if (withinTolerance) + return TransformToXYPlane(closestCartesianDirection, out backTransform); + + var dirMagSquared = direction.LengthSquared(); + if (dirMagSquared.IsPracticallySame(0.0)) + throw new ArgumentException("The direction cannot be the zero vector."); + // Normalize the direction. If already unit length, this is a no-op multiplication by ~1. + var oneOverH = 1 / Math.Sqrt(dirMagSquared); + var dx = direction.X * oneOverH; + var dy = direction.Y * oneOverH; + var dz = direction.Z * oneOverH; + + // g = length of the projection of the unit direction onto the XZ plane. + // When g ≈ 0 the direction is along ±Y and the two-Givens-rotation approach + // (rotate in XZ then YZ) is singular. Handle that case with a simple 90° rotation. + var g = Math.Sqrt(dx * dx + dz * dz); + if (g < 1e-12) + { + // direction is essentially ±Y. Rotate 90° around the Z-axis so that Y maps to Z. + // For +Y: new X = old X, new Y = old Z, new Z = old Y (right-handed) + // For -Y: flip sign on Z row. + if (dy > 0) + { + backTransform = new Matrix4x4(1, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0); + return new Matrix4x4(1, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0); + } + else + { + backTransform = new Matrix4x4(1, 0, 0, 0, 0, -1, 0, -1, 0, 0, 0, 0); + return new Matrix4x4(1, 0, 0, 0, 0, -1, 0, -1, 0, 0, 0, 0); + } + } + var oneOverG = 1 / g; + var xOverG = dx * oneOverG; + var zOverG = dz * oneOverG; + // backTransform columns are the new X, Y, Z axes expressed in the original frame. + // Row 0 (new X axis): rotate in XZ plane to eliminate X component → (z/g, 0, -x/g) + // Row 1 (new Y axis): perpendicular to both → (-x*y/g, g, -z*y/g) + // Row 2 (new Z axis): the normalized direction itself → (dx, dy, dz) + backTransform = new Matrix4x4(zOverG, 0, -xOverG, -xOverG * dy, g, -zOverG * dy, dx, dy, dz, 0, 0, 0); + return backTransform.Transpose(); + } /// @@ -2533,7 +2714,7 @@ public static bool IsVertexInsideAABB(this Vertex vertexInQuestion, Solid solid, /// If true, on boundary is inside. /// A bool. public static bool IsPointInsideAABB(this Vector3 pointInQuestion, Solid solid, - bool onBoundaryIsInside = true) + bool onBoundaryIsInside = true) => IsPointInsideAABB(pointInQuestion, solid.Bounds[0], solid.Bounds[1], onBoundaryIsInside); /// @@ -2758,7 +2939,7 @@ public static Vector4 Unique3DLine(Vector3 anchor, Vector3 direction) var PolarAngle = Math.Acos(direction.Z * oneOverRadius); // when z_dir is 1, then angle is zero or pi var AzimuthAngle = Math.Atan2(direction.Y, direction.X); // regardless of length, we can find azimuth by tangent var inPlaneYDir = new Vector3(-direction.Y, direction.X, 0).Normalize(); // the y-dir will never have a z-component - it is like - // the latitude lines on a globe + // the latitude lines on a globe var inPlaneXDir = // here, x-dir determined by y cross z where z is the given direction // cross product is i = (y1z2 - y2z1), j = (x2z1 - x1z2), k = (x1y2 - x2y1) // x1 = -direction.Y, y1 = direction.X, z1=0 diff --git a/TessellationAndVoxelizationGeometryLibrary/Numerics/PolynomialSolve.cs b/TessellationAndVoxelizationGeometryLibrary/Numerics/PolynomialSolve.cs index eeaf9a5f1..adcb163e1 100644 --- a/TessellationAndVoxelizationGeometryLibrary/Numerics/PolynomialSolve.cs +++ b/TessellationAndVoxelizationGeometryLibrary/Numerics/PolynomialSolve.cs @@ -120,7 +120,7 @@ public static (ComplexNumber, ComplexNumber) Quadratic(double squaredCoeff, doub { if (squaredCoeff.IsNegligible()) { - return (new ComplexNumber(-constant / linearCoeff), ComplexNumber.NaN); + return (new ComplexNumber(-constant / linearCoeff), new ComplexNumber(-constant / linearCoeff)); } if ((constant / squaredCoeff).IsNegligible()) { @@ -128,20 +128,22 @@ public static (ComplexNumber, ComplexNumber) Quadratic(double squaredCoeff, doub } var oneOverDenom = 1 / (2 * squaredCoeff); var radicalTerm = linearCoeff * linearCoeff - 4 * squaredCoeff * constant; // more commonly known as b^2 - 4ac - if (radicalTerm < 0) // then imaginary roots + if (radicalTerm.IsNegativeNonNegligible()) // then imaginary roots { radicalTerm = Math.Sqrt(-radicalTerm); radicalTerm *= oneOverDenom; var negBTerm = -oneOverDenom * linearCoeff; return (new ComplexNumber(negBTerm, -radicalTerm), new ComplexNumber(negBTerm, radicalTerm)); } - else + if (radicalTerm.IsPositiveNonNegligible()) // then two distinct real roots { - radicalTerm = Math.Sqrt(radicalTerm); + radicalTerm = radicalTerm > 0 ? Math.Sqrt(radicalTerm) : 0; radicalTerm *= oneOverDenom; var negBTerm = oneOverDenom * linearCoeff; return (new ComplexNumber(radicalTerm - negBTerm), new ComplexNumber(-radicalTerm - negBTerm)); } + var singleRoot = new ComplexNumber(-linearCoeff * oneOverDenom, 0); + return (singleRoot, singleRoot); } /// diff --git a/TessellationAndVoxelizationGeometryLibrary/Polygon Operations/Polygon3D.cs b/TessellationAndVoxelizationGeometryLibrary/Polygon Operations/Polygon3D.cs new file mode 100644 index 000000000..9dff4fff7 --- /dev/null +++ b/TessellationAndVoxelizationGeometryLibrary/Polygon Operations/Polygon3D.cs @@ -0,0 +1,137 @@ +// *********************************************************************** +// Assembly : TessellationAndVoxelizationGeometryLibrary +// Author : matth +// Created : 04-03-2023 +// +// Last Modified By : matth +// Last Modified On : 04-03-2023 +// *********************************************************************** +// +// 2026 +// +// +// *********************************************************************** +using System.Collections.Generic; +using System.Linq; + + +namespace TVGL +{ + /// + /// Represents a 3D polygon, consisting of vertices and optional holes. + /// + public class Polygon3D + { + /// + /// Initializes a new instance of the class. + /// + public Polygon3D() { } + + /// + /// Initializes a new instance of the class with vertices, a closed flag, and an index. + /// + /// The vertices. + /// if set to true, the polygon is closed. + /// The index. + public Polygon3D(IEnumerable vertices, bool isClosed, int index = -1) + { + Vertices = vertices as IList ?? vertices.ToList(); + Index = index; + IsClosed = isClosed; + } + + /// + /// Initializes a new instance of the class with vertices, a closed flag, and holes. + /// + /// The vertices. + /// if set to true, the polygon is closed. + /// The holes in the polygon. + public Polygon3D(IEnumerable vertices, bool isClosed, IEnumerable> holes) + : this(vertices, isClosed) + { + if (holes is ICollection> collection) + { + var numHoles = collection.Count; + Holes = new IList[numHoles]; + var i = 0; + foreach (var hole in holes) + Holes[i++] = hole as IList ?? hole.ToArray(); + } + else + { + Holes = new List>(); + foreach (var hole in holes) + Holes.Add(hole as IList ?? hole.ToArray()); + } + } + + /// + /// Gets the primary outer vertices of the polygon. + /// + public IList Vertices { get; init; } + + /// + /// Gets or sets the holes in the polygon. + /// + public IList> Holes { get; set; } + + /// + /// Gets a value indicating whether this is closed. + /// + public bool IsClosed { get; } + + /// + /// Iterates over all paths starting with the outer perimeter and + /// then any holes. + /// + /// An enumeration of all paths. + public IEnumerable> AllPaths() + { + yield return Vertices; + if (Holes != null) + foreach (var hole in Holes) + yield return hole; + } + + /// + /// Gets the booleans for all the paths in the 3D polygon + /// starting with the outer one. All holes are considerd closed anyway, + /// so the output is always [t/f t t t ...] + /// + /// An enumeration of boolean closed states. + public IEnumerable AllIsClosed() + { + yield return IsClosed; + if (Holes != null) + foreach (var hole in Holes) + yield return true; + } + + /// + /// Creates a deep copy of this . + /// + /// A new object. + public Polygon3D Copy() + { + var copiedHoles = new Vector3[Holes.Count][]; + for (int i = 0; i < Holes.Count; i++) + copiedHoles[i] = Holes[i].ToArray(); + return new Polygon3D(Vertices.ToArray(), IsClosed, copiedHoles); + } + + /// + /// Gets the total number of paths (1 for the outer perimeter + number of holes). + /// + public int NumPaths => 1 + (Holes?.Count ?? 0); + + /// + /// Gets the total number of points across all paths. + /// + public int NumPoints => AllPaths().Sum(path => path.Count); + + /// + /// Gets or sets the index of the polygon. + /// + public int Index { get; set; } + } +} diff --git a/TessellationAndVoxelizationGeometryLibrary/Polygon Operations/PolygonOperations.Triangulate.cs b/TessellationAndVoxelizationGeometryLibrary/Polygon Operations/PolygonOperations.Triangulate.cs index df4371284..1e4543d71 100644 --- a/TessellationAndVoxelizationGeometryLibrary/Polygon Operations/PolygonOperations.Triangulate.cs +++ b/TessellationAndVoxelizationGeometryLibrary/Polygon Operations/PolygonOperations.Triangulate.cs @@ -40,7 +40,8 @@ public static partial class PolygonOperations /// if set to true [force to positive]. /// IEnumerable<Vertex[]> where each represents a triangular polygonal face. /// The vertices must all have a unique IndexInList value - vertexLoop - public static IEnumerable<(Vertex A, Vertex B, Vertex C)> Triangulate(this IEnumerable vertexLoop, Vector3 normal, bool forceToPositive = false) + public static IEnumerable<(Vertex A, Vertex B, Vertex C)> Triangulate(this IEnumerable vertexLoop, + Vector3 normal, bool forceToPositive = false, bool handleSelfIntersects = true, double suggestedAngle = 0.0) { var transform = normal.TransformToXYPlane(out _); var coords = new List(); @@ -54,7 +55,7 @@ public static partial class PolygonOperations } var polygon = new Polygon(coords); if (forceToPositive && !polygon.IsPositive) polygon.IsPositive = true; - foreach (var triangleIndices in polygon.TriangulateToIndices()) + foreach (var triangleIndices in polygon.TriangulateToIndices(handleSelfIntersects,suggestedAngle)) { if (indexToVertexDict[triangleIndices.A] != indexToVertexDict[triangleIndices.B] && indexToVertexDict[triangleIndices.B] != indexToVertexDict[triangleIndices.C] @@ -72,7 +73,8 @@ public static partial class PolygonOperations /// The normal direction. /// IEnumerable<Vertex[]> where each represents a triangular polygonal face. /// The vertices must all have a unique IndexInList value - vertexLoops - public static IEnumerable Triangulate(this IEnumerable> vertexLoops, Vector3 normal) + public static IEnumerable Triangulate(this IEnumerable> vertexLoops, Vector3 normal, + bool handleSelfIntersects = true, double suggestedAngle = 0.0) { var transform = normal.TransformToXYPlane(out _); var polygons = new List(); @@ -92,7 +94,7 @@ public static IEnumerable Triangulate(this IEnumerable> polygons = polygons.CreateShallowPolygonTrees(false); foreach (var polygon in polygons) { - foreach (var triangleIndices in polygon.TriangulateToIndices()) + foreach (var triangleIndices in polygon.TriangulateToIndices(handleSelfIntersects, suggestedAngle)) yield return new[] {indexToVertexDict[triangleIndices.A], indexToVertexDict[triangleIndices.B], indexToVertexDict[triangleIndices.C]}; @@ -105,9 +107,10 @@ public static IEnumerable Triangulate(this IEnumerable> /// /// The polygon. /// List<System.Int32[]>. - public static IEnumerable TriangulateToCoordinates(this Polygon polygon) + public static IEnumerable TriangulateToCoordinates(this Polygon polygon, bool handleSelfIntersects = true, + double suggestedAngle = 0.0) { - foreach (var triangle in polygon.Triangulate()) + foreach (var triangle in polygon.Triangulate(handleSelfIntersects, suggestedAngle)) yield return new[] { triangle[0].Coordinates, triangle[1].Coordinates, triangle[2].Coordinates }; } @@ -117,27 +120,28 @@ public static IEnumerable TriangulateToCoordinates(this Polygon polyg /// The polygon. /// if set to true [handle self intersects]. /// List<System.Int32[]>. - public static IEnumerable<(int A, int B, int C)> TriangulateToIndices(this Polygon polygon, bool handleSelfIntersects = true) + public static IEnumerable<(int A, int B, int C)> TriangulateToIndices(this Polygon polygon, bool handleSelfIntersects = true, + double suggestedAngle = 0.0) { var vertexIndices = new HashSet(); - var needToReIndex = false; - foreach (var v in polygon.AllPolygons.SelectMany(p => p.Vertices)) - { - if (vertexIndices.Contains(v.IndexInList)) - { - needToReIndex = true; - break; - } - vertexIndices.Add(v.IndexInList); - } + //var needToReIndex = false; + //foreach (var v in polygon.AllPolygons.SelectMany(p => p.Vertices)) + //{ + // if (vertexIndices.Contains(v.IndexInList)) + // { + // needToReIndex = true; + // break; + // } + // vertexIndices.Add(v.IndexInList); + //} var index = 0; - if (needToReIndex) - { + //if (needToReIndex) + //{ foreach (var subPolygon in polygon.AllPolygons) foreach (var vertex in subPolygon.Vertices) vertex.IndexInList = index++; - } - foreach (var triangle in polygon.Triangulate(handleSelfIntersects)) + //} + foreach (var triangle in polygon.Triangulate(handleSelfIntersects, suggestedAngle)) yield return (triangle[0].IndexInList, triangle[1].IndexInList, triangle[2].IndexInList); } /// @@ -148,7 +152,8 @@ public static IEnumerable TriangulateToCoordinates(this Polygon polyg /// List<System.Int32[]>. /// Triangulate Polygon requires a positive polygon. A negative one was provided. - polygon /// Unable to triangulate polygon. - public static List Triangulate(this Polygon polygon, bool handleSelfIntersects = true) + public static List Triangulate(this Polygon polygon, bool handleSelfIntersects = true, + double suggestedAngle =0.0) { if (polygon.Area.IsNegligible() || (polygon.IsConvex && !polygon.InnerPolygons.Any())) { @@ -176,7 +181,7 @@ public static List Triangulate(this Polygon polygon, bool handleSelf if (concaveEdge != null) { while (verts[0].IndexInList != concaveEdge.IndexInList) - verts = new[] { verts[1], verts[2], verts[3], verts[0] }; + verts = [verts[1], verts[2], verts[3], verts[0]]; } return new List { new[] { verts[0], verts[1], verts[2] }, new[] { verts[0], verts[2], verts[3] } }; } @@ -184,7 +189,7 @@ public static List Triangulate(this Polygon polygon, bool handleSelf { var triangleList = new List(); for (int i = 2; i < polygon.Vertices.Count; i++) - triangleList.Add(new[] { polygon.Vertices[0], polygon.Vertices[i - 1], polygon.Vertices[i] }); + triangleList.Add([polygon.Vertices[0], polygon.Vertices[i - 1], polygon.Vertices[i]]); return triangleList; } var triangleFaceList = new List(); @@ -193,7 +198,7 @@ public static List Triangulate(this Polygon polygon, bool handleSelf // in case this is a deep polygon tree - recurse down and solve for the inner positive polygons foreach (var hole in polygon.InnerPolygons) foreach (var smallInnerPolys in hole.InnerPolygons) - triangleFaceList.AddRange(Triangulate(smallInnerPolys)); + triangleFaceList.AddRange(Triangulate(smallInnerPolys,handleSelfIntersects,suggestedAngle)); var selfIntersections = polygon.GetSelfIntersections().Where(intersect => intersect.Relationship != SegmentRelationship.NoOverlap).ToList(); if (selfIntersections.Count > 0) @@ -216,7 +221,7 @@ public static List Triangulate(this Polygon polygon, bool handleSelf var attempts = 0; var random = new Random(1); var successful = false; - var angle = 0.0; + var angle = suggestedAngle; var localTriangleFaceList = new List(); do { @@ -227,7 +232,7 @@ public static List Triangulate(this Polygon polygon, bool handleSelf if (angle != 0) { var rotateMatrix = new Matrix3x3(c, s, -s, c, 0, 0); - polygon.Transform(rotateMatrix); //This destructive operation implies that polygon is mutable (i.e., polygons should not be used in multiple locations). + polygon.Transform(rotateMatrix); //This destructively alters polygon coordinates, but if used - we rotate back in 22 lines } try { diff --git a/TessellationAndVoxelizationGeometryLibrary/TessellatedSolid/ModifyTessellation/ComplexifyTessellation.cs b/TessellationAndVoxelizationGeometryLibrary/TessellatedSolid/ModifyTessellation/ComplexifyTessellation.cs index da8a35f17..8802325e6 100644 --- a/TessellationAndVoxelizationGeometryLibrary/TessellatedSolid/ModifyTessellation/ComplexifyTessellation.cs +++ b/TessellationAndVoxelizationGeometryLibrary/TessellatedSolid/ModifyTessellation/ComplexifyTessellation.cs @@ -14,7 +14,6 @@ using System; using System.Collections.Generic; using System.Linq; -using System.Net.Http.Headers; namespace TVGL { @@ -84,81 +83,106 @@ public static void Complexify(TessellatedSolid ts, int targetNumberOfFaces, doub /// When this method returns, contains the list of vertices that were added during the operation. /// When this method returns, contains the list of triangle faces that were added during the operation. /// The desired total number of faces in the mesh after the operation. Must be a non-negative integer. - /// The maximum allowable length for edges to be considered for subdivision. Must be a positive value. - public static void Complexify(IEnumerable edges, out List addedEdges, out List addedVertices, out List addedFaces, - int targetNumberOfFaces, double maxLength) + /// The maximum allowable length for edges to be considered for subdivision. Must be a positive value. + public static void Complexify(IEnumerable edges, out List addedEdges, out List addedVertices, + out List addedFaces, int targetNumberOfFaces, double maxSurfaceDeviation) { - var edgeQueue = new PriorityQueue(edges.Where(e => e.Length >= maxLength).Select(e => (e, e.Length)), - new ReverseSort()); + //var edgeLengthList = edges.OrderByDescending(e => e.Length).ToArray(); + var initEdgePlot = edges.Select(e => new[] { e.From.Coordinates, e.To.Coordinates }).ToArray(); + var edgeQueue = new PriorityQueue<(Edge, Vector3), double>(new ReverseSort()); + foreach (var edge in edges) + EnqueueEdgeAndFindNewPoint(edgeQueue, edge); addedEdges = new List(); addedVertices = new List(); addedFaces = new List(); var iterations = targetNumberOfFaces > 0 ? (int)Math.Ceiling(targetNumberOfFaces / 2.0) : targetNumberOfFaces; - while (iterations-- != 0 && edgeQueue.TryDequeue(out var edge, out var edgeLength)) + var edgeCounter = edgeQueue.Count; + while (iterations-- != 0 && edgeQueue.TryDequeue(out (Edge edge, Vector3 mpt) c, out _)) { - var origLeftFace = edge.OtherFace; - var origRightFace = edge.OwnedFace; - var leftFarVertex = origLeftFace?.OtherVertex(edge); - var rightFarVertex = origRightFace?.OtherVertex(edge); - var fromVertex = edge.From; - var toVertex = edge.To; - var primitive1 = origLeftFace?.BelongsToPrimitive; - var primitive2 = origRightFace == null ? primitive1 : origRightFace.BelongsToPrimitive; - if (primitive1 == null) primitive1 = primitive2; - var commonPrimitive = (primitive1 == primitive2) ? primitive1 : null; - var addedVertex = new Vertex(DetermineIntermediateVertexPosition(fromVertex.Coordinates, - toVertex.Coordinates, commonPrimitive)); + //var map = edgeLengthList.IndexOf(c.edge); + //Console.WriteLine(map); + //if (iterations % 1000 <= 0) + // Presenter.ShowAndHang([initEdgePlot, addedEdges.Select(e => new[] { e.From.Coordinates, e.To.Coordinates }), [[c.edge.From.Coordinates, c.mpt, c.edge.To.Coordinates]]], + // [false, false, false], colors: [new Color(KnownColors.LightGray), new Color(KnownColors.Blue), new Color(KnownColors.Red)]); + var origLeftFace = c.edge.OtherFace; + var origRightFace = c.edge.OwnedFace; + var leftFarVertex = origLeftFace?.OtherVertex(c.edge); + var rightFarVertex = origRightFace?.OtherVertex(c.edge); + var fromVertex = c.edge.From; + var toVertex = c.edge.To; + var addedVertex = new Vertex(c.mpt); // modify original faces with new intermediate vertex origLeftFace?.ReplaceVertex(toVertex, addedVertex); origRightFace?.ReplaceVertex(toVertex, addedVertex); - var newLeftFace = new TriangleFace(toVertex, addedVertex, leftFarVertex); - var newRightFace = new TriangleFace(addedVertex, toVertex, rightFarVertex); + TriangleFace newLeftFace = null, newRightFace = null; + if (leftFarVertex != null) + newLeftFace = new TriangleFace(toVertex, addedVertex, leftFarVertex) + { BelongsToPrimitive = origLeftFace.BelongsToPrimitive }; + if (rightFarVertex != null) + newRightFace = new TriangleFace(addedVertex, toVertex, rightFarVertex) + { BelongsToPrimitive = origRightFace.BelongsToPrimitive }; toVertex.Faces.Remove(origLeftFace); toVertex.Faces.Remove(origRightFace); - var inlineEdge = new Edge(addedVertex, toVertex, newRightFace, newLeftFace, true); - toVertex.Edges.Remove(edge); - edge.To = addedVertex; - addedVertex.Edges.Add(edge); - edge.Update(); - var newLeftEdge = new Edge(leftFarVertex, addedVertex, origLeftFace, newLeftFace, true); - var newRightEdge = new Edge(rightFarVertex, addedVertex, newRightFace, origRightFace, true); - origLeftFace.AddEdge(newLeftEdge); - origRightFace.AddEdge(newRightEdge); - var bottomEdge = toVertex.Edges.First(e => e.OtherVertex(toVertex) == leftFarVertex); - if (bottomEdge.OwnedFace == origLeftFace) - bottomEdge.OwnedFace = newLeftFace; - else bottomEdge.OtherFace = newLeftFace; - newLeftFace.AddEdge(bottomEdge); - bottomEdge.Update(); - - bottomEdge = toVertex.Edges.First(e => e.OtherVertex(toVertex) == rightFarVertex); - if (bottomEdge.OwnedFace == origRightFace) - bottomEdge.OwnedFace = newRightFace; - else bottomEdge.OtherFace = newRightFace; - newRightFace.AddEdge(bottomEdge); - bottomEdge.Update(); + var inlineEdge = new Edge(addedVertex, toVertex, newRightFace, newLeftFace, true, edgeCounter++); + toVertex.Edges.Remove(c.edge); + c.edge.To = addedVertex; + addedVertex.Edges.Add(c.edge); + c.edge.Update(); + Edge newLeftEdge = null, newRightEdge = null; + if (leftFarVertex != null) + { + newLeftEdge = new Edge(leftFarVertex, addedVertex, origLeftFace, newLeftFace, true, edgeCounter++); + origLeftFace.AddEdge(newLeftEdge); + var bottomEdge = toVertex.Edges.First(e => e.OtherVertex(toVertex) == leftFarVertex); + if (bottomEdge.OwnedFace == origLeftFace) + bottomEdge.OwnedFace = newLeftFace; + else bottomEdge.OtherFace = newLeftFace; + newLeftFace.AddEdge(bottomEdge); + bottomEdge.Update(); + } + if (rightFarVertex != null) + { + newRightEdge = new Edge(rightFarVertex, addedVertex, newRightFace, origRightFace, true, edgeCounter++); + origRightFace.AddEdge(newRightEdge); + var bottomEdge = toVertex.Edges.First(e => e.OtherVertex(toVertex) == rightFarVertex); + if (bottomEdge.OwnedFace == origRightFace) + bottomEdge.OwnedFace = newRightFace; + else bottomEdge.OtherFace = newRightFace; + newRightFace.AddEdge(bottomEdge); + bottomEdge.Update(); + } // need to re-add the edge. It was modified in the SplitEdge function (now, half the lenght), but // it may still be met by this criteria - if (edge.Length >= maxLength) - edgeQueue.Enqueue(edge, edge.Length); - if (inlineEdge.Length >= maxLength) - edgeQueue.Enqueue(inlineEdge, inlineEdge.Length); - if (newLeftEdge.Length >= maxLength) - edgeQueue.Enqueue(newLeftEdge, newLeftEdge.Length); - if (newRightEdge.Length >= maxLength) - edgeQueue.Enqueue(newRightEdge, newRightEdge.Length); + EnqueueEdgeAndFindNewPoint(edgeQueue, c.edge); + EnqueueEdgeAndFindNewPoint(edgeQueue, inlineEdge); + if (newLeftEdge != null) + EnqueueEdgeAndFindNewPoint(edgeQueue, newLeftEdge); + if (newRightEdge != null) + EnqueueEdgeAndFindNewPoint(edgeQueue, newRightEdge); - addedEdges.Add(inlineEdge); - addedEdges.Add(newLeftEdge); - addedEdges.Add(newRightEdge); - addedFaces.Add(newLeftFace); - addedFaces.Add(newRightFace); addedVertices.Add(addedVertex); + addedEdges.Add(inlineEdge); + if (newLeftEdge != null) addedEdges.Add(newLeftEdge); + if (newRightEdge != null) addedEdges.Add(newRightEdge); + if (newLeftFace != null) addedFaces.Add(newLeftFace); + if (newRightFace != null) addedFaces.Add(newRightFace); } } + + + private static void EnqueueEdgeAndFindNewPoint(PriorityQueue<(Edge, Vector3), double> edgeQueue, Edge edge) + { + var midPoint = DetermineIntermediateVertexPosition(edge); + var lineUnitVector = edge.UnitVector; + var t = lineUnitVector.Dot(midPoint - edge.From.Coordinates); + var pointOnLine = edge.From.Coordinates + t * lineUnitVector; + var distanceToSurf = (midPoint - pointOnLine).LengthSquared(); + var newLength = 2 * distanceToSurf + 2 * t * t + edge.Vector.LengthSquared() - 2 * t * edge.Length; + edgeQueue.Enqueue((edge, midPoint), newLength); + } } } \ No newline at end of file diff --git a/TessellationAndVoxelizationGeometryLibrary/TessellatedSolid/ModifyTessellation/DetermineIntermediateVertex.cs b/TessellationAndVoxelizationGeometryLibrary/TessellatedSolid/ModifyTessellation/DetermineIntermediateVertex.cs index f8c02b57f..dd9936b31 100644 --- a/TessellationAndVoxelizationGeometryLibrary/TessellatedSolid/ModifyTessellation/DetermineIntermediateVertex.cs +++ b/TessellationAndVoxelizationGeometryLibrary/TessellatedSolid/ModifyTessellation/DetermineIntermediateVertex.cs @@ -24,24 +24,26 @@ namespace TVGL /// public static partial class ModifyTessellation { - /// - /// Adjusts the position of kept vertex. - /// - /// The keep vertex. - /// The other vertex. - /// Vector3. + + internal static Vector3 DetermineIntermediateVertexPosition(Edge edge) + { + var fromVertex = edge.From; + var toVertex = edge.To; + var primitive1 = edge.OwnedFace?.BelongsToPrimitive; + var primitive2 = edge.OtherFace == null ? primitive1 : edge.OtherFace.BelongsToPrimitive; + if (primitive1 == null) primitive1 = primitive2; + if (primitive1 == primitive2) + return DetermineIntermediateVertexPosition(edge.From.Coordinates, edge.To.Coordinates, primitive1); + else + return DetermineIntermediateVertexPosition(edge.From.Coordinates, edge.To.Coordinates, primitive1, primitive2); + } + internal static Vector3 DetermineIntermediateVertexPosition(Vector3 pt1, Vector3 pt2, - PrimitiveSurface primitive = null) + PrimitiveSurface primitive) { if (primitive == null || primitive is Plane) // then average positions return 0.5 * (pt1 + pt2); //else we're going to project the average position onto the primitive surface - if (primitive is Sphere sphere) - // the direction from the center to the average of the two points is the direction to - // move out to the sphere surface. Note we can cheese the average by the Normalize - // function because the direction is all we care about, but we need to find the - // relative movements so, the points are each subtracted from the center - return sphere.Center + sphere.Radius * (pt1 + pt2 - 2 * sphere.Center).Normalize(); if (primitive is Cylinder cylinder) { var anchorToPt1 = pt1 - cylinder.Anchor; @@ -75,17 +77,76 @@ internal static Vector3 DetermineIntermediateVertexPosition(Vector3 pt1, Vector3 axisT1 * (pt2 - pt2DiskCenter)).Normalize(); } if (primitive is GeneralQuadric quadric) - { - var p1Normal = quadric.GetNormalAtPoint(pt1); - var p2Normal = quadric.GetNormalAtPoint(pt2); - var pCenter = 0.5 * (pt1 + pt2); - var planeNormal = (p1Normal.Cross(p2Normal)).Normalize(); - var outwardDir = (pt2 - pt1).Cross(planeNormal); - // the point we seek is farthest from the line connecting pt1 to pt2 - // we don't know exactly where it is along the line, but we know its - // normal should be inline with the outwardDir. - return quadric.PointsWithGradient(outwardDir).MinBy(q => q.DistanceSquared(pCenter)); + { // this one is a bit complicated. The obvious solution to divide the line in half and then project + // it to the quadric is not the best, but it is kept here as a fallback at the end. Instead, we first + // find the plane made by the normal of the midpoint and the original line segment (it would seem better + // to use the endpoints' normals, but this was not the case). From this plane, we slice to find the conic + // on that plane made by the quadric, then find the point on the curve between the end points that is farthest + // from the original line. This point has the gradient (normal of the conic) defined by the outwardDirc + var midpoint = 0.5 * (pt1 + pt2); + var prevVector = pt2 - pt1; + var midPtNormal = quadric.GetNormalAtPoint(midpoint); + var planenormal = midPtNormal.Cross(prevVector).Normalize(); + var outwarddir = prevVector.Cross(planenormal).Normalize(); + var plane = new Plane(planenormal.Dot(pt1), planenormal); + var midPoint2d = midpoint.ConvertTo2DCoordinates(planenormal, out _); + var outwarddir2d = outwarddir.ConvertTo2DCoordinates(planenormal, out _); + var conic = GeneralConicSection.CreateFromQuadric(quadric, plane); + var posSuccess = conic.PointsAtGivenGradient(outwarddir2d, out var posPt); + var negSuccess = conic.PointsAtGivenGradient(-outwarddir2d, out var negPt); + var oddAxis = GetOddPrincipalAxis(quadric); + var stationaryPt = quadric.StationaryPoint; + var returnPt = Vector3.Null; + if (oddAxis.Dot(midpoint - stationaryPt) < 0) + { + oddAxis = -oddAxis; // make sure the odd axis is pointing in the same direction as the midpoint. + } + if (posSuccess && (!negSuccess || posPt.DistanceSquared(midPoint2d) < negPt.DistanceSquared(midPoint2d))) + { + var posPt3 = posPt.ConvertTo3DLocation(plane.AsTransformFromXYPlane); + //if (oddAxis.IsNull() || oddAxis.Dot(posPt3 - stationaryPt) > 0) + returnPt = posPt3; + } + else if (negSuccess) + { + var negPt3 = negPt.ConvertTo3DLocation(plane.AsTransformFromXYPlane); + //if (oddAxis.IsNull() || oddAxis.Dot(negPt3 - stationaryPt) > 0) + returnPt = negPt3; + } + + // fallback plan 1 + var midPointProjections = quadric.LineIntersection(midpoint, outwarddir); + + if (midPointProjections.Any()) + { + // Fallback: project the midpoint onto the quadric along the surface normal + var bestT = double.PositiveInfinity; + Vector3 result = midpoint; + foreach (var (intersection, t) in midPointProjections) + { + if (Math.Abs(t) < Math.Abs(bestT)) + { + bestT = t; + result = intersection; + } + } + //if (oddAxis.IsNull() || oddAxis.Dot(result - stationaryPt) > 0) + if (returnPt.IsNull() || returnPt.DistanceSquared(midpoint) > result.DistanceSquared(midpoint)) + returnPt = result; + } + // fall back plan 2 + var result2 = new Plane(stationaryPt, oddAxis).LineIntersection(midpoint, outwarddir).First().intersection; + + if (returnPt.IsNull() || returnPt.DistanceSquared(midpoint) > result2.DistanceSquared(midpoint)) + returnPt = result2; + return returnPt; } + if (primitive is Sphere sphere) + // the direction from the center to the average of the two points is the direction to + // move out to the sphere surface. Note we can cheese the average by the Normalize + // function because the direction is all we care about, but we need to find the + // relative movements so, the points are each subtracted from the center + return sphere.Center + sphere.Radius * (pt1 + pt2 - 2 * sphere.Center).Normalize(); if (primitive is Torus torus) { var anchorToPt1 = pt1 - torus.Center; @@ -97,6 +158,32 @@ internal static Vector3 DetermineIntermediateVertexPosition(Vector3 pt1, Vector3 throw new NotImplementedException(); } + + internal static Vector3 GetOddPrincipalAxis(GeneralQuadric quadric) + { + // The principal axes are the axes along which the rate of change of the gradient is maximized or minimized. + // This can be found by solving for the eigenvectors of the Hessian matrix of the quadric function. + // In particular we want the odd principal axis, which is the one that is different in sign from the others. + // This is the same as the cone or hyperboloid axis. + Matrix3x3 coeffs = quadric.GetHessian(); + var eigenVals = coeffs.GetEigenValuesAndVectors(out var eigenVecs); + + for (int i = 0; i < eigenVals.Length; i++) + { + if (eigenVals[i].Real < 0 && eigenVals.Aggregate(1.0, (a, b) => a * b.Real) < 0 || + eigenVals[i].Real > 0 && eigenVals.Aggregate(1.0, (a, b) => a * b.Real) > 0) + return new Vector3(eigenVecs[i][0].Real, eigenVecs[i][1].Real, eigenVecs[i][2].Real); + } + return Vector3.Null; + } + + + + private static Vector3 DetermineIntermediateVertexPosition(Vector3 coordinates1, Vector3 coordinates2, PrimitiveSurface primitive1, PrimitiveSurface primitive2) + { + throw new NotImplementedException(); + } + /// /// Adjusts the position of kept vertex experimental. /// diff --git a/TessellationAndVoxelizationGeometryLibrary/TessellatedSolid/ModifyTessellation/ModifyTessellation.AddRemoveElements.cs b/TessellationAndVoxelizationGeometryLibrary/TessellatedSolid/ModifyTessellation/ModifyTessellation.AddRemoveElements.cs index 2be49ed5b..efff8fe44 100644 --- a/TessellationAndVoxelizationGeometryLibrary/TessellatedSolid/ModifyTessellation/ModifyTessellation.AddRemoveElements.cs +++ b/TessellationAndVoxelizationGeometryLibrary/TessellatedSolid/ModifyTessellation/ModifyTessellation.AddRemoveElements.cs @@ -180,7 +180,8 @@ public static bool MergeVertexAndKill3EdgesAnd2Faces(this Vertex removedVertex, removedFace1.BelongsToPrimitive == removedFace2.BelongsToPrimitive) keepVertex.Coordinates = DetermineIntermediateVertexPosition(keepVertex.Coordinates, removedVertex.Coordinates, removedFace1.BelongsToPrimitive); - keepVertex.Coordinates = DetermineIntermediateVertexPosition(keepVertex.Coordinates, removedVertex.Coordinates); + keepVertex.Coordinates = DetermineIntermediateVertexPosition(keepVertex.Coordinates, removedVertex.Coordinates, removedFace1.BelongsToPrimitive, + removedFace2.BelongsToPrimitive); foreach (var e in keepVertex.Edges) e.Update(); foreach (var f in keepVertex.Faces) diff --git a/TessellationAndVoxelizationGeometryLibrary/TessellatedSolid/ModifyTessellation/SimplifyTessellation.cs b/TessellationAndVoxelizationGeometryLibrary/TessellatedSolid/ModifyTessellation/SimplifyTessellation.cs index f1e92e499..a2ff28944 100644 --- a/TessellationAndVoxelizationGeometryLibrary/TessellatedSolid/ModifyTessellation/SimplifyTessellation.cs +++ b/TessellationAndVoxelizationGeometryLibrary/TessellatedSolid/ModifyTessellation/SimplifyTessellation.cs @@ -56,7 +56,7 @@ public static void SimplifyFlatPatches(this TessellatedSolid ts) foreach (var triangle in triangulatedList) { var newFace = new TriangleFace(triangle, flat.Normal); - if (newFace.Area.IsNegligible() && newFace.Normal.IsNull()) continue; + if (newFace.Area.IsNegligible() && newFace.Normal.IsNegligible()) continue; newFaces.Add(newFace); foreach (var fromVertex in newFace.Vertices) { diff --git a/TessellationAndVoxelizationGeometryLibrary/TessellatedSolid/ModifyTessellation/Single3DPolygonTriangulation.cs b/TessellationAndVoxelizationGeometryLibrary/TessellatedSolid/ModifyTessellation/Single3DPolygonTriangulation.cs index 5448b25f4..7cb0f53d9 100644 --- a/TessellationAndVoxelizationGeometryLibrary/TessellatedSolid/ModifyTessellation/Single3DPolygonTriangulation.cs +++ b/TessellationAndVoxelizationGeometryLibrary/TessellatedSolid/ModifyTessellation/Single3DPolygonTriangulation.cs @@ -577,11 +577,11 @@ private static double CalcObjFunction(double dotWeight, TriangulationLoop triang //return area; //return area * (2 - neighborNormal.Dot(triangle.Normal)); var dotPenalty = 0.0; - if (!neighborNormal1.IsNull()) + if (!neighborNormal1.IsNegligible()) dotPenalty += (1 - neighborNormal1.Dot(triangle.Normal)); - if (!neighborNormal2.IsNull()) + if (!neighborNormal2.IsNegligible()) dotPenalty += (1 - neighborNormal2.Dot(triangle.Normal)); - if (!neighborNormal3.IsNull()) + if (!neighborNormal3.IsNegligible()) dotPenalty += (1 - neighborNormal3.Dot(triangle.Normal)); foreach (var ead in triangle) { diff --git a/TessellationAndVoxelizationGeometryLibrary/TessellatedSolid/ModifyTessellation/TessellationInspectAndRepair.cs b/TessellationAndVoxelizationGeometryLibrary/TessellatedSolid/ModifyTessellation/TessellationInspectAndRepair.cs index ec2d9bf51..7aa152d1e 100644 --- a/TessellationAndVoxelizationGeometryLibrary/TessellatedSolid/ModifyTessellation/TessellationInspectAndRepair.cs +++ b/TessellationAndVoxelizationGeometryLibrary/TessellatedSolid/ModifyTessellation/TessellationInspectAndRepair.cs @@ -672,7 +672,7 @@ bool SetNegligibleAreaFaceNormals() { foreach (var adjFace in facesToCheck[i].AdjacentFaces) { - if (!adjFace._normal.IsNull()) + if (!adjFace._normal.IsNegligible()) { facesToCheck[i]._normal = adjFace._normal; facesToCheck.RemoveAt(i); diff --git a/TessellationAndVoxelizationGeometryLibrary/TessellatedSolid/TriangleFace.cs b/TessellationAndVoxelizationGeometryLibrary/TessellatedSolid/TriangleFace.cs index 1872d6576..b40dab12d 100644 --- a/TessellationAndVoxelizationGeometryLibrary/TessellatedSolid/TriangleFace.cs +++ b/TessellationAndVoxelizationGeometryLibrary/TessellatedSolid/TriangleFace.cs @@ -239,9 +239,13 @@ public override Vector3 Normal get { if (_normal.IsNull()) - _normal = ((B.Coordinates - A.Coordinates) - .Cross(C.Coordinates - A.Coordinates)) - .Normalize(); + { + _normal = (B.Coordinates - A.Coordinates) + .Cross(C.Coordinates - A.Coordinates); + if (_normal.LengthSquared().IsNegligible()) + _normal = Vector3.Zero; + else _normal.Normalize(); + } return _normal; } } @@ -306,12 +310,9 @@ public IEnumerable Edges { get { - if (AB != null) - yield return AB; - if (BC != null) - yield return BC; - if (CA != null) - yield return CA; + yield return AB; + yield return BC; + yield return CA; } }