diff --git a/Scripts/ConvexHullCalculator.cs b/Scripts/ConvexHullCalculator.cs index 9da0f75..f7c1e35 100644 --- a/Scripts/ConvexHullCalculator.cs +++ b/Scripts/ConvexHullCalculator.cs @@ -7,8 +7,8 @@ namespace GK { /// - /// A not very optimized implementation of the quickhull algorithm - /// for generating 3d convex hulls. + /// An implementation of the quickhull algorithm for generating + /// 3d convex hulls. /// /// The algorithm works like this: you start with an initial /// "seed" hull, that is just a simple tetrahedron made up of @@ -54,8 +54,8 @@ namespace GK { /// the cone of new faces. /// /// There are a number of things you could possibly optimize in - /// this version of the algorithm, but this is correct and - /// reasonably fast. Fast enough for my purposes, anyway. + /// this version of the algorithm, but this should be correct + /// and reasonably fast. /// /// A note: the code uses a right-handed coordinate system, where /// the cross-product uses the right-hand rule and the faces are @@ -81,7 +81,7 @@ public class ConvexHullCalculator { const int INSIDE = -1; /// - /// Epsilon value. + /// Epsilon value. /// const float EPSILON = 0.0001f; @@ -98,8 +98,6 @@ public class ConvexHullCalculator { /// Vertex2 and Vertex1), etc. /// /// Normal is (unsurprisingly) the normal of the triangle. - /// - /// TODO: SOA instead of AOS? /// struct Face { public int Vertex0; @@ -138,7 +136,7 @@ public bool Equals(Face other) { /// face. These are used in the openSet array. /// /// Point is the index of the point in the points array, - /// Face is the key of the face in the Key dictioanry, + /// Face is the key of the face in the Key dictionary, /// Distance is the distance from the face to the point. /// struct PointFace { @@ -188,8 +186,7 @@ struct HorizonEdge { /// the Faces EVER created in the algorithm, and skip the /// ones that have been marked as deleted. However, looping /// through a list is fairly fast, and it might be worth it - /// to avoid Dictionary overhead. I should probably test and - /// profile both ways. + /// to avoid Dictionary overhead. /// /// TODO test converting to a List instead. /// @@ -333,34 +330,21 @@ void Initialize(List points, bool splitVerts) { /// /// Create initial seed hull. - /// - /// The good way to do this is probably to find the extreme - /// points and create the seed hull from those, but I'm just - /// using the four first points for it, for now. Obviously - /// it can be optimized :) /// void GenerateInitialHull(List points) { - // TODO use extreme points to generate seed hull. I wonder - // how much difference that actually makes, you would - // imagine that even with a tiny seed hull, it would grow - // pretty quickly. Anyway, the rest should be the same, - // you only need to change how you find b0/b1/b2/b3 - - // TODO i'm a bit worried what happens if these points are - // too close to each other or if the fourth point is - // coplanar with the triangle. I should loop through the - // point set to find suitable points instead. - var b0 = 0; - var b1 = 1; - var b2 = 2; - var b3 = 3; + // Find points suitable for use as the seed hull. Some + // varieties of this algorithm pick extreme points here, + // but I'm not convinced you gain all that much from that. + // Currently what it does is just find the first four + // points that are not coplanar. + int b0, b1, b2, b3; + FindInitialHullIndices(points, out b0, out b1, out b2, out b3); var v0 = points[b0]; var v1 = points[b1]; var v2 = points[b2]; var v3 = points[b3]; - // TODO use epsilon? var above = Dot(v3 - v1, Cross(v1 - v0, v2 - v0)) > 0.0f; // Create the faces of the seed hull. You need to draw a @@ -468,6 +452,43 @@ void GenerateInitialHull(List points) { VerifyOpenSet(points); } + /// + /// Find four points in the point cloud that are not + /// coplanar for the seed hull + /// + void FindInitialHullIndices(List points, out int b0, out int b1, out int b2, out int b3) { + var count = points.Count; + + for (int i0 = 0; i0 < count - 3; i0++) { + for (int i1 = i0 + 1; i1 < count - 2; i1++) { + var p0 = points[i0]; + var p1 = points[i1]; + + if (AreCoincident(p0, p1)) continue; + + for (int i2 = i1 + 1; i2 < count - 1; i2++) { + var p2 = points[i2]; + + if (AreCollinear(p0, p1, p2)) continue; + + for (int i3 = i2 + 1; i3 < count - 0; i3++) { + var p3 = points[i3]; + + if(AreCoplanar(p0, p1, p2, p3)) continue; + + b0 = i0; + b1 = i1; + b2 = i2; + b3 = i3; + return; + } + } + } + } + + throw new System.ArgumentException("Can't generate hull, points are coplanar"); + } + /// /// Grow the hull. This method takes the current hull, and /// expands it to encompass the point in openSet with the @@ -923,6 +944,32 @@ static Vector3 Cross(Vector3 a, Vector3 b) { a.x*b.y - a.y*b.x); } + /// + /// Check if two points are coincident + /// + [MethodImpl(MethodImplOptions.AggressiveInlining)] + bool AreCoincident(Vector3 a, Vector3 b) { + return (a - b).magnitude <= EPSILON; + } + + /// + /// Check if three points are collinear + /// + [MethodImpl(MethodImplOptions.AggressiveInlining)] + bool AreCollinear(Vector3 a, Vector3 b, Vector3 c) { + return Cross(c - a, c - b).magnitude <= EPSILON; + } + + /// + /// Check if four points are coplanar + /// + [MethodImpl(MethodImplOptions.AggressiveInlining)] + bool AreCoplanar(Vector3 a, Vector3 b, Vector3 c, Vector3 d) { + var n1 = Cross(c - a, c - b); + var n2 = Cross(d - a, d - b); + + return AreCollinear(Vector3.zero, n1, n2); + } /// /// Method used for debugging, verifies that the openSet is