Changeset 9945 for branches/HeuristicLab.Analysis.AlgorithmBehavior/MIConvexHull

Timestamp:

09/10/13 21:06:55 (11 years ago)

Author:

ascheibe

Message:

#1886

• added new convex hull algorithm based on SMO/SVM
• added unit test and a visualization
• updated MIConvexHull

Location:

branches/HeuristicLab.Analysis.AlgorithmBehavior

Files:

• . (modified) (1 prop)
• MIConvexHull/Constants.cs (added)
• MIConvexHull/ConvexHull/ConvexFaceInternal.cs (modified) (1 diff)
• MIConvexHull/ConvexHull/ConvexHullInternal.cs (modified) (18 diffs)
• MIConvexHull/ConvexHull/VertexSort.cs (added)

branches/HeuristicLab.Analysis.AlgorithmBehavior/MIConvexHull/ConvexHull/ConvexFaceInternal.cs

@@ -40 +40 @@
             return x.Index.CompareTo(y.Index);
         }
+    }
+
+    /// <summary>
+    /// For deferred face addition.
+    /// </summary>
+    sealed class DeferredFace
+    {
+        /// <summary>
+        /// The faces.
+        /// </summary>
+        public ConvexFaceInternal Face, Pivot, OldFace;
+
+        /// <summary>
+        /// The indices.
+        /// </summary>
+        public int FaceIndex, PivotIndex;
     }
 

branches/HeuristicLab.Analysis.AlgorithmBehavior/MIConvexHull/ConvexHull/ConvexHullInternal.cs

@@ -42 +42 @@
         VertexBuffer BeyondBuffer;
         List<ConvexFaceInternal> AffectedFaceBuffer;
+        List<DeferredFace> ConeFaceBuffer;
+        Stack<DeferredFace> ConeFaceStack;
+        HashSet<VertexWrap> SingularVertices;
 
         const int ConnectorTableSize = 2017;
@@ -69 +72 @@
             EmptyBuffer = new VertexBuffer();
             AffectedFaceBuffer = new List<ConvexFaceInternal>();
+            ConeFaceBuffer = new List<DeferredFace>();
+            ConeFaceStack = new Stack<DeferredFace>();
+            SingularVertices = new HashSet<VertexWrap>();
             BeyondBuffer = new VertexBuffer();
 
@@ -130 +136 @@
         /// <param name="face"></param>
         /// <returns></returns>
-        private void CalculateFacePlane(ConvexFaceInternal face)
+        private bool CalculateFacePlane(ConvexFaceInternal face)
         {
             var vertices = face.Vertices;
@@ -136 +142 @@
             FindNormalVector(vertices, normal);
 
-            if (double.IsNaN(normal[0])) ThrowSingular();
+            if (double.IsNaN(normal[0]))
+            {
+                return false;
+                //ThrowSingular();
+            }
 
             double offset = 0.0;
@@ -157 +167 @@
             }
             else face.IsNormalFlipped = false;
+
+            return true;
         }
         
@@ -318 +330 @@
 
         /// <summary>
+        /// Creates a new deferred face.
+        /// </summary>
+        /// <param name="face"></param>
+        /// <param name="faceIndex"></param>
+        /// <param name="pivot"></param>
+        /// <param name="pivotIndex"></param>
+        /// <param name="oldFace"></param>
+        /// <returns></returns>
+        DeferredFace GetDeferredFace(ConvexFaceInternal face, int faceIndex, ConvexFaceInternal pivot, int pivotIndex, ConvexFaceInternal oldFace)
+        {
+            var ret = ConeFaceStack.Count > 0 ? ConeFaceStack.Pop() : new DeferredFace();
+            
+            ret.Face = face;
+            ret.FaceIndex = faceIndex;
+            ret.Pivot = pivot;
+            ret.PivotIndex = pivotIndex;
+            ret.OldFace = oldFace;
+
+            return ret;
+        }
+
+        /// <summary>
         /// Connect faces using a connector.
         /// </summary>
@@ -326 +360 @@
             var list = ConnectorTable[index];
 
-            int count = 0;
             for (var current = list.First; current != null; current = current.Next)
             {
@@ -347 +380 @@
         /// Removes the faces "covered" by the current vertex and adds the newly created ones.
         /// </summary>
-        private void CreateCone()
-        {
-            var oldFaces = AffectedFaceBuffer;
-
+        private bool CreateCone()
+        {
             var currentVertexIndex = CurrentVertex.Index;
-
-            for (int fIndex = 0; fIndex < oldFaces.Count; fIndex++)
-            {
-                var oldFace = oldFaces[fIndex];
+            ConeFaceBuffer.Clear();
+
+            for (int fIndex = 0; fIndex < AffectedFaceBuffer.Count; fIndex++)
+            {
+                var oldFace = AffectedFaceBuffer[fIndex];
 
                 // Find the faces that need to be updated
@@ -368 +400 @@
                         ++updateCount;
                     }
-                }
-
-                // Recycle the face for future use
-                if (updateCount == 0)
-                {
-                    // If the face is present in the unprocessed list, remove it 
-                    UnprocessedFaces.Remove(oldFace);
-
-                    RecycleFace(oldFace);
-                    RecycledFaceStack.Push(oldFace);
                 }
 
@@ -402 +424 @@
                     VertexWrap[] vertices;
 
-                    // Recycle the oldFace
-                    if (i == updateCount - 1)
-                    {
-                        RecycleFace(oldFace);
-                        newFace = oldFace;
-                        vertices = newFace.Vertices;
-                        oldVertexIndex = vertices[forbidden].Index;
-                    }
-                    else // Pop a face from the recycled stack or create a new one
-                    {
-                        newFace = GetNewFace();
-                        vertices = newFace.Vertices;                        
-                        for (int j = 0; j < Dimension; j++) vertices[j] = oldFace.Vertices[j];
-                        oldVertexIndex = vertices[forbidden].Index;
-                    }
+                    newFace = GetNewFace();
+                    vertices = newFace.Vertices;                        
+                    for (int j = 0; j < Dimension; j++) vertices[j] = oldFace.Vertices[j];
+                    oldVertexIndex = vertices[forbidden].Index;
 
                     int orderedPivotIndex;
@@ -450 +461 @@
                     vertices[orderedPivotIndex] = CurrentVertex;
 
-                    CalculateFacePlane(newFace);
-                    newFace.AdjacentFaces[orderedPivotIndex] = adjacentFace;
-                    adjacentFace.AdjacentFaces[oldFaceAdjacentIndex] = newFace;
-
-                    // let there be a connection.
-                    for (int j = 0; j < Dimension; j++)
+                    if (!CalculateFacePlane(newFace))
                     {
-                        if (j == orderedPivotIndex) continue;
-                        var connector = GetNewConnector();
-                        connector.Update(newFace, j, Dimension);
-                        ConnectFace(connector);
+                        return false;
                     }
-                    
-                    // This could slightly help...
-                    if (adjacentFace.VerticesBeyond.Count < oldFace.VerticesBeyond.Count)
-                    {
-                        FindBeyondVertices(newFace, adjacentFace.VerticesBeyond, oldFace.VerticesBeyond);
-                    }
-                    else
-                    {
-                        FindBeyondVertices(newFace, oldFace.VerticesBeyond, adjacentFace.VerticesBeyond);
-                    }
-
-                    // This face will definitely lie on the hull
-                    if (newFace.VerticesBeyond.Count == 0)
-                    {
-                        ConvexFaces.Add(newFace);
-                        UnprocessedFaces.Remove(newFace);
-                        EmptyBufferStack.Push(newFace.VerticesBeyond);
-                        newFace.VerticesBeyond = EmptyBuffer;
-                    }
-                    else // Add the face to the list
-                    {
-                        UnprocessedFaces.Add(newFace);
-                    }
-                }
+
+                    ConeFaceBuffer.Add(GetDeferredFace(newFace, orderedPivotIndex, adjacentFace, oldFaceAdjacentIndex, oldFace));
+                }
+            }
+            
+            return true;
+        }
+
+        /// <summary>
+        /// Commits a cone and adds a vertex to the convex hull.
+        /// </summary>
+        void CommitCone()
+        {
+            // Add the current vertex.
+            ConvexHull.Add(CurrentVertex);
+            
+            // Fill the adjacency.
+            for (int i = 0; i < ConeFaceBuffer.Count; i++)
+            {
+                var face = ConeFaceBuffer[i];
+
+                var newFace = face.Face;
+                var adjacentFace = face.Pivot;
+                var oldFace = face.OldFace;
+                var orderedPivotIndex = face.FaceIndex;
+
+                newFace.AdjacentFaces[orderedPivotIndex] = adjacentFace;
+                adjacentFace.AdjacentFaces[face.PivotIndex] = newFace;
+
+                // let there be a connection.
+                for (int j = 0; j < Dimension; j++)
+                {
+                    if (j == orderedPivotIndex) continue;
+                    var connector = GetNewConnector();
+                    connector.Update(newFace, j, Dimension);
+                    ConnectFace(connector);
+                }
+
+                // This could slightly help...
+                if (adjacentFace.VerticesBeyond.Count < oldFace.VerticesBeyond.Count)
+                {
+                    FindBeyondVertices(newFace, adjacentFace.VerticesBeyond, oldFace.VerticesBeyond);
+                }
+                else
+                {
+                    FindBeyondVertices(newFace, oldFace.VerticesBeyond, adjacentFace.VerticesBeyond);
+                }
+
+                // This face will definitely lie on the hull
+                if (newFace.VerticesBeyond.Count == 0)
+                {
+                    ConvexFaces.Add(newFace);
+                    UnprocessedFaces.Remove(newFace);
+                    EmptyBufferStack.Push(newFace.VerticesBeyond);
+                    newFace.VerticesBeyond = EmptyBuffer;
+                }
+                else // Add the face to the list
+                {
+                    UnprocessedFaces.Add(newFace);
+                }
+
+                // recycle the object.
+                ConeFaceStack.Push(face);
+            }
+
+            // Recycle the affected faces.
+            for (int fIndex = 0; fIndex < AffectedFaceBuffer.Count; fIndex++)
+            {
+                var face = AffectedFaceBuffer[fIndex];
+                UnprocessedFaces.Remove(face);
+                RecycleFace(face);
+                RecycledFaceStack.Push(face);
             }
         }
@@ -698 +749 @@
         void UpdateCenter()
         {
-            for (int i = 0; i < Dimension; i++) Center[i] *= (ConvexHull.Count - 1);
-            double f = 1.0 / ConvexHull.Count;
+            var count = ConvexHull.Count + 1;
+            for (int i = 0; i < Dimension; i++) Center[i] *= (count - 1);
+            double f = 1.0 / count;
             for (int i = 0; i < Dimension; i++) Center[i] = f * (Center[i] + CurrentVertex.PositionData[i]);
+        }
+
+        /// <summary>
+        /// Removes the last vertex from the center.
+        /// </summary>
+        void RollbackCenter()
+        {
+            var count = ConvexHull.Count + 1;
+            for (int i = 0; i < Dimension; i++) Center[i] *= count;
+            double f = 1.0 / (count - 1);
+            for (int i = 0; i < Dimension; i++) Center[i] = f * (Center[i] - CurrentVertex.PositionData[i]);
         }
 
@@ -715 +778 @@
             {
                 CurrentVertex = vertex;
+                // update center must be called before adding the vertex.
+                UpdateCenter();
                 ConvexHull.Add(CurrentVertex);
-                UpdateCenter();
                 InputVertices.Remove(vertex);
 
@@ -767 +831 @@
             for (int i = 2; i <= Dimension; i++)
             {
-                double maximum = 0.0001;
+                double maximum = 0.000001;
                 VertexWrap maxPoint = null;
                 for (int j = 0; j < extremes.Count; j++)
@@ -817 +881 @@
         double GetSimplexVolume(VertexWrap pivot, List<VertexWrap> initialPoints)
         {
-            var dim = initialPoints.Count;
-            var m = nDMatrix;
-
-            for (int i = 0; i < dim; i++)
-            {
-                var pts = initialPoints[i];
-                for (int j = 0; j < dim; j++) m[i, j] = pts.PositionData[j] - pivot.PositionData[j];
-            }
-
-            return Math.Abs(StarMath.determinantDestructive(m, dim));
+            var initPtsNum = initialPoints.Count;
+            //var m = nDMatrix;
+            var sum = 0.0;
+        
+            for (int i = 0; i < initPtsNum; i++)
+            {
+                var initPt = initialPoints[i];
+                for (int j = 0; j < Dimension; j++)
+                    sum += (initPt.PositionData[j] - pivot.PositionData[j]) 
+                        * (initPt.PositionData[j] - pivot.PositionData[j]);
+            }
+
+            return sum;
         }
 
@@ -878 +945 @@
 
         /// <summary>
+        /// Handles singular vertex.
+        /// </summary>
+        void HandleSingular()
+        {
+            RollbackCenter();
+            SingularVertices.Add(CurrentVertex);
+
+            // This means that all the affected faces must be on the hull and that all their "vertices beyond" are singular.
+            for (int fIndex = 0; fIndex < AffectedFaceBuffer.Count; fIndex++)
+            {
+                var face = AffectedFaceBuffer[fIndex];
+                var vb = face.VerticesBeyond;
+                for (int i = 0; i < vb.Count; i++)
+                {
+                    SingularVertices.Add(vb[i]);
+                }
+
+                ConvexFaces.Add(face);
+                UnprocessedFaces.Remove(face);
+                EmptyBufferStack.Push(face.VerticesBeyond);
+                face.VerticesBeyond = EmptyBuffer;
+            }
+        }
+
+        /// <summary>
         /// Fins the convex hull.
         /// </summary>
@@ -890 +982 @@
                 var currentFace = UnprocessedFaces.First;
                 CurrentVertex = currentFace.FurthestVertex;
-                ConvexHull.Add(CurrentVertex);
+                                                
                 UpdateCenter();
 
@@ -897 +989 @@
 
                 // Create the cone from the currentVertex and the affected faces horizon.
-                CreateCone();
+                if (!SingularVertices.Contains(CurrentVertex) && CreateCone()) CommitCone();
+                else HandleSingular();
 
                 // Need to reset the tags