Fix collapsed line bug, remove representative line, add tolerance

demo/demo_helper.cpp

@@ -222,11 +222,11 @@ void DemoConfiguration::displaySettings()
     std::cout << "[Source Size: " << scaling << "]";
     std::cout << "[Source Vertices: " << vertexCount0 << " and " << vertexCount1 << "]";
     std::cout << "[GlobalScaling:" << globalScaling << "]";
+    std::cout << "[Tolerance:" << tolerance << "]";
     std::cout << std::endl;
 
     std::cout << "  [Early stop: " << getStatusString(detectApertureOnly) << "]";
     std::cout << "[Silhouette: " << getStatusString(silhouetteOptimisation) << "]";
-    std::cout << "[Middle line: " << getStatusString(representativeLineSampling) << "]";
     std::cout << "[Normalization: " << getStatusString(normalization) << "]" << std::endl;
 
 #if EXACT_ARITHMETIC            
@@ -249,15 +249,14 @@ void DemoConfiguration::writeConfig(const std::string& filename)
     output << "vertexCount1 = " << vertexCount1 << std::endl;
     output << "silhouetteOptimisation  = " << silhouetteOptimisation << std::endl;
     output << "detectApertureOnly  = " << detectApertureOnly << std::endl;
-    output << "representativeLineSampling  = " << representativeLineSampling << std::endl;
     output << "normalization  = " << normalization << std::endl;
     output << "scaling  = " << scaling << std::endl;
     output << "phi  = " << phi << std::endl;
     output << "eta = " << eta << std::endl;
     output << "sceneIndex = " << sceneIndex << std::endl;
     output << "globalScaling = " << globalScaling << std::endl;
     output << "precisionType = " << precisionType << std::endl;
-
+    output << "tolerance = " << tolerance << std::endl;
     output.close();
 }
 
@@ -278,15 +277,14 @@ void DemoConfiguration::readConfig(const std::string& filename)
         else if (tokens[0] == "vertexCount1") { vertexCount1 = atoi(tokens[2].c_str()); }
         else if (tokens[0] == "silhouetteOptimisation") { silhouetteOptimisation = atoi(tokens[2].c_str()); }
         else if (tokens[0] == "detectApertureOnly") { detectApertureOnly = atoi(tokens[2].c_str()); }
-        else if (tokens[0] == "representativeLineSampling") { representativeLineSampling = atoi(tokens[2].c_str()); }
         else if (tokens[0] == "normalization") { normalization = atoi(tokens[2].c_str()); }
         else if (tokens[0] == "scaling") { scaling = (float)atof(tokens[2].c_str()); }
         else if (tokens[0] == "phi") { phi = (float)atof(tokens[2].c_str()); }
         else if (tokens[0] == "eta") { eta = (float)atof(tokens[2].c_str()); }
         else if (tokens[0] == "sceneIndex") { sceneIndex = atoi(tokens[2].c_str()); }
         else if (tokens[0] == "globalScaling") { globalScaling = (float)atof(tokens[2].c_str()); }
         else if (tokens[0] == "precisionType") { precisionType = (VisibilityExactQueryConfiguration::PrecisionType)atoi(tokens[2].c_str()); }
-
+        else if (tokens[0] == "tolerance") { tolerance =  (double)atof(tokens[2].c_str()); }
         else { V_ASSERT(0); }
     }
     input.close();
@@ -295,7 +293,7 @@ void DemoConfiguration::readConfig(const std::string& filename)
 
 void DemoConfiguration::displaySummary()
 {
-    std::cout << "VisibilityTest [Scaling: " << globalScaling << ", v0: " << vertexCount0 << ", vv1: " << vertexCount1 << "; phi:" << phi << "; precision: "<< precisionType<<"] ";
+    std::cout << "VisibilityTest [Scaling: " << globalScaling << ", v0: " << vertexCount0 << ", vv1: " << vertexCount1 << "; phi:" << phi << "; precision: "<< precisionType<<"; tolerance:" << tolerance << "] ";
 }
 
 #if EMBREE

demo/demo_helper.h

@@ -34,7 +34,6 @@ namespace visilibDemo
         size_t vertexCount0 = 3;
         size_t vertexCount1 = 3;
         bool silhouetteOptimisation = true;
-        bool representativeLineSampling = true;
         VisibilityExactQueryConfiguration::PrecisionType precisionType = VisibilityExactQueryConfiguration::DOUBLE;
         bool detectApertureOnly = false;
         bool normalization = true;
@@ -43,6 +42,7 @@ namespace visilibDemo
         float eta = 0;
         int   sceneIndex = 2;
         float globalScaling = 1;
+        double tolerance = -1;
 #if EMBREE
         bool embree = false;
 #endif

demo/demo_main.cpp

@@ -92,8 +92,9 @@ namespace visilibDemo
             config.silhouetteOptimization = mDemoConfiguration.silhouetteOptimisation;
             config.hyperSphereNormalization = mDemoConfiguration.normalization;
             config.precision = mDemoConfiguration.precisionType;
-            config.representativeLineSampling = mDemoConfiguration.representativeLineSampling;
+
             config.detectApertureOnly = mDemoConfiguration.detectApertureOnly;
+            config.tolerance = mDemoConfiguration.tolerance;
 #if EMBREE 
             config.useEmbree = mDemoConfiguration.embree;
 #endif
@@ -155,7 +156,6 @@ namespace visilibDemo
 #if EMBREE
             std::cout << "  g: enable/disable embree ray tracing" << std::endl;
 #endif
-            std::cout << "  r: enable/disable representative line sampling strategy" << std::endl;
             std::cout << "  f: enable/disable detect aperture only" << std::endl;
 
             std::cout << "  n: enable/disable Plucker normalization" << std::endl;
@@ -164,6 +164,8 @@ namespace visilibDemo
             std::cout << "  x: change scene " << std::endl;
             std::cout << "  +/-: increase/decrease scaling of query polygons" << std::endl;
             std::cout << "  1/2: increase/decrease number of vertices of query polygons" << std::endl;
+            std::cout << "  /-*: increase/decrease global scaling" << std::endl;
+            std::cout << "  [/]: increase/decrease scaling" << std::endl;
 
             std::cout << "  w: write config" << std::endl;
             std::cout << "  o: open config" << std::endl;
@@ -231,6 +233,17 @@ namespace visilibDemo
 
                 break;
 
+            case '[':                
+                mDemoConfiguration.tolerance = mDemoConfiguration.tolerance == -1 ? 1e-8 : mDemoConfiguration.tolerance/2;
+                forceDisplay = true;                                
+
+                break;
+
+            case ']':
+                mDemoConfiguration.tolerance = mDemoConfiguration.tolerance == -1 ? 1e-8 : mDemoConfiguration.tolerance*2;
+                forceDisplay = true;                
+                break;
+
             case 'h':
 
                 writeHelp();
@@ -254,10 +267,6 @@ namespace visilibDemo
                 initScene(mDemoConfiguration.sceneIndex);
                 forceDisplay = true;
                 break;
-            case 'r':
-                mDemoConfiguration.representativeLineSampling = !mDemoConfiguration.representativeLineSampling;
-                forceDisplay = true;
-                break;
 
 
 #ifdef EXACT_ARITHMETIC

test/src/test_visibility.cpp

@@ -41,7 +41,6 @@ bool testConfiguration(const DemoConfiguration& configuration, VisibilityResult
     config.silhouetteOptimization = true;
     config.hyperSphereNormalization = configuration.normalization;
     config.precision = configuration.precisionType;
-    config.representativeLineSampling = true;
     config.detectApertureOnly = false;
 
     int s = 2;

visilib/math_combinatorial.h

@@ -33,6 +33,7 @@ namespace visilib
     public:
         static bool haveAtLeastNCommonFacets(const std::vector<size_t>& aFacetsDescription1, const std::vector<size_t>& aFacetsDescription2, size_t n = 3);
         static void initFacets(const std::vector<size_t>& aFacetsDescription1, const std::vector<size_t>& aFacetsDescription2, size_t anHyperplane, std::vector<size_t>& result);
+        static void initFacets(const std::vector<size_t>& aFacetsDescription1, const std::vector<size_t>& aFacetsDescription2, std::vector<size_t>& result);
         static bool hasFacet(const std::vector<size_t>& facets, size_t aFace);
     };
 
@@ -76,6 +77,42 @@ namespace visilib
         V_ASSERT(std::is_sorted(result.begin(), result.end()));
     }
 
+    /** @brief Initialize a new facets description result from the two input facets description aFacetsDescription1 an aFacetsDescription1 of two given vertices
+    The resulting facets description is a sorted list of facets, containing the intersection set of the two input facets descriptions list and the additional facet anHyperplane.
+    Remark: the precondition is that the input facets list are sorted .*/
+
+    inline void MathCombinatorial::initFacets(const std::vector<size_t>& aFacetsDescription1, const std::vector<size_t>& aFacetsDescription2, std::vector<size_t>& result)
+    {
+        int i = 0, j = 0;
+        int n1 = aFacetsDescription1.size();
+        int n2 = aFacetsDescription2.size();
+
+        while (i < n1 && j < n2)
+        {
+            if (aFacetsDescription1[i] < aFacetsDescription2[j])
+            {
+                aFacetsDescription1[i++];
+            }
+            else if (aFacetsDescription1[i] > aFacetsDescription2[j])
+            {
+                result.push_back(aFacetsDescription2[j++]);
+            }
+            else
+            { 
+                i++;j++;                
+            }
+        }
+
+        while (i < n1)
+        {
+            result.push_back(aFacetsDescription1[i++]);
+        }
+        while (j < n2)
+        {
+            result.push_back(aFacetsDescription2[j++]);
+        }
+    }
+
     /** @brief Determine if the sorted facets description contains aFace. Behaviour is only defined if aFacetsDescription is sorted.*/
 
     inline bool MathCombinatorial::hasFacet(const std::vector<size_t> & aFacetsDescription, size_t aFace)

visilib/plucker_polytope.h

@@ -117,7 +117,16 @@ namespace visilib
         @param S: the tolerance against which the length of the edges is tested
         */
 
-        template<class S> bool hasAllEdgesCollapsed(PluckerPolyhedron<P>* polyhedron, S Tolerance);
+        template<class S> bool hasSomeEdgesCollapsed(PluckerPolyhedron<P>* polyhedron, S Tolerance);
+
+        /** @brief Remove the collasped edge (the length of all edges are below the given tolerance).
+        
+        @param polyhedron: the polyhedron containing all the vertices definition
+        @param S: the tolerance against which the length of the edges is tested
+        */
+
+        template<class S> bool removeCollapsedEdges(PluckerPolyhedron<P>* polyhedron, S Tolerance);
+
 
         /** @brief Compute the extremal stabbing lines of the polytope*/
         template<class S> void computeExtremalStabbingLines(PluckerPolyhedron<P>* polyhedron, S Tolerance);
@@ -267,15 +276,13 @@ namespace visilib
     template<class P>
     void PluckerPolytope<P>::outputProperties(std::ostream& o, PluckerPolyhedron<P>* polyhedron)
     {
-        return;
-        o << "Polyhedron Lines:" << polyhedron->getLinesCount() << std::endl;
-
         o << "Polytope ESL: " << mExtremalStabbingLines.size() << std::endl;
         o << "Polytope Edges: " << mEdges.size() << std::endl;
         o << "Polytope Vertices: " << mVertices.size() << std::endl;
 
         if (polyhedron != nullptr)
         {
+            o << "Polyhedron Lines:" << polyhedron->getLinesCount() << std::endl;
             std::set<size_t> myFacets;
             getFacets(myFacets, polyhedron);
 
@@ -302,6 +309,7 @@ namespace visilib
                 o << std::endl;
             }
 
+        }
             o << "Polytope edges facets: " << std::endl;
 
             int line = 0;
@@ -310,24 +318,126 @@ namespace visilib
                 o << " e[" << iter->first << "," << iter->second << "]";
             }
             o << std::endl;
+
+        }
+    template<class P> template<class S>
+    bool PluckerPolytope<P>::removeCollapsedEdges(PluckerPolyhedron<P>* polyhedron, S tolerance)
+    {
+        std::vector<int> myEdgesTable;
+        std::vector<int> myMergeTable;
+
+        for (auto iter = mEdges.begin(); iter != mEdges.end();iter++)
+        {
+            size_t i1 = iter->first;
+            size_t i2 = iter->second;                
+
+            const P& v1 = polyhedron->get(i1);
+            const P& v2 = polyhedron->get(i2);
+
+
+            if (MathPredicates::isEdgeCollapsed(v1, v2, tolerance))
+            {
+                std::cout << "collapsed edge dected: " << v1-v2 << "; tolerance: "<< 0.5 *tolerance<<std::endl;
+                myMergeTable.push_back(i1);                
+                myMergeTable.push_back(i2);
+            }
+            else
+            {
+                myEdgesTable.push_back(i1);
+                myEdgesTable.push_back(i2);
+            }
+        }
+
+        if (myMergeTable.empty())
+        {
+            return false;            
         }
+
+        std::cout << "EDGE TO BE MERGED DETECTED!" << std::endl;
+        outputProperties(std::cout,polyhedron);
+
+        std::cout << "myMergeTable table:";
+        for (int i=0; i<myMergeTable.size();i++)
+        {
+            std::cout<<myMergeTable[i]<<" ";
+        }
+        std::cout<<std::endl;
+        std::cout << "myEdgesTable table:";
+        for (int i=0; i<myEdgesTable.size();i++)
+        {
+            std::cout<<myEdgesTable[i]<<" ";
+        }
+        std::cout<<std::endl;
+        for (int i=0; i<myMergeTable.size();i+=2)
+        {
+            int merge_left  = myMergeTable[i];
+            int merge_right = myMergeTable[i+1];
+            if (merge_left!=merge_right)
+            {
+                for (int j=0; j<myEdgesTable.size();j++)
+                {
+                    if (myEdgesTable[j] == merge_right)
+                    {
+                        myEdgesTable[j] = merge_left;
+                    }
+                }
+
+                for (int k=i;k<myMergeTable.size();k++)
+                {
+                    if (myMergeTable[k] == merge_right)
+                    {
+                        myMergeTable[k] = merge_left;
+                    }
+                }
+
+                const std::vector<size_t>& facetsI1 = polyhedron->getFacetsDescription(merge_left);
+                const std::vector<size_t>& facetsI2 = polyhedron->getFacetsDescription(merge_right);
+
+                std::vector<size_t> myFacets;
+                MathCombinatorial::initFacets(facetsI1, 
+                                              facetsI2,                                               
+                                              myFacets);
+                polyhedron->initFacetsDescription(merge_left, myFacets);
+            }
+        }
+
+        std::cout << "myMergeTable table:";
+        for (int i=0; i<myMergeTable.size();i++)
+        {
+            std::cout<<myMergeTable[i]<<" ";
+        }
+        std::cout<<std::endl;
+        std::cout << "myEdgesTable table:";
+        for (int i=0; i<myEdgesTable.size();i++)
+        {
+            std::cout<<myEdgesTable[i]<<" ";
+        }
+        std::cout<<std::endl;   
+
+        outputProperties(std::cout,polyhedron);
+        mEdges.clear();
+        for (int i=0; i<myEdgesTable.size();i+=2)
+        {
+            addEdge(myEdgesTable[i], myEdgesTable[i+1],polyhedron);
+        }
+        return true;
     }
 
     template<class P> template<class S>
-    bool PluckerPolytope<P>::hasAllEdgesCollapsed(PluckerPolyhedron<P>* polyhedron, S tolerance)
+    bool PluckerPolytope<P>::hasSomeEdgesCollapsed(PluckerPolyhedron<P>* polyhedron, S tolerance)
     {
         for (auto iter = mEdges.begin(); iter != mEdges.end(); iter++)
         {
             const P& v1 = polyhedron->get(iter->first);
             const P& v2 = polyhedron->get(iter->second);
 
-            if (!MathPredicates::isEdgeCollapsed(v1, v2, tolerance))
+            if (MathPredicates::isEdgeCollapsed(v1, v2, tolerance))
             {
-                return false;
+                return true;
             }
         }
 
-        return true;
+        return false;
     }
 
     template<class P> template<class S>
@@ -421,9 +531,9 @@ namespace visilib
             isValid = false;
         }
 
-        if (hasAllEdgesCollapsed(polyhedron, tolerance))
+        if (hasSomeEdgesCollapsed(polyhedron, tolerance))
         {
-            std::cerr << "   IsValid error: hasAllEdgesCollapesd" << std::endl;
+            std::cerr << "   IsValid error: hasSomeEdgesCollapsed" << std::endl;
             outputProperties(std::cout);
 
             isValid = false;
@@ -452,6 +562,7 @@ namespace visilib
                     std::cerr << "   IsValid error: polytope vertex " << *iter << " is not on the facet " << *i << std::endl;
                     std::cerr << "        P:" << polyhedron->get(*iter) << std::endl;
                     std::cerr << "        Facet:" << polyhedron->get(*i) << std::endl;
+                    std::cerr << "        Facet.dot(P):" <<  polyhedron->get(*iter).dot(polyhedron->get(*i)) << std::endl;
                     isValid = false;
                 }
             }