[fuseCut] adjustments to the tetrahedral vote

src/aliceVision/fuseCut/DelaunayGraphCut.hpp

@@ -97,7 +97,11 @@ class DelaunayGraphCut
         {
             return v0 == e.v0 && v1 == e.v1;
         }
-
+        bool isSameUndirectionalEdge(const Edge& e) const
+        {
+            return (v0 == e.v0 && v1 == e.v1) ||
+                   (v0 == e.v1 && v1 == e.v0);
+        }
     };
 
     enum class EGeometryType
@@ -145,8 +149,9 @@ class DelaunayGraphCut
             case EGeometryType::Facet:
                 return facet == g.facet;
             case EGeometryType::None:
-                return true;
+                break;
             }
+            return true;
         }
         bool operator!=(const GeometryIntersection& g) const
         {
@@ -185,7 +190,13 @@ class DelaunayGraphCut
             edges += gc.edges;
             vertices += gc.vertices;
             facets += gc.facets;
-
+            return *this;
+        }
+        GeometriesCount& operator/=(const size_t v)
+        {
+            edges /= v;
+            vertices /= v;
+            facets /= v;
             return *this;
         }
     };
@@ -431,8 +442,14 @@ class DelaunayGraphCut
 
     void fuseFromDepthMaps(const StaticVector<int>& cams, const Point3d voxel[8], const FuseParams& params);
 
-    void computeVerticesSegSize(bool allPoints, float alpha = 0.0f);
-    void removeSmallSegs(int minSegSize);
+    /**
+     * @brief Compute connected segments size
+     * @param[out] out_segments
+     * @param[in] useVertex: sub-set of vertices to compute
+     * @param[in] alpha
+     */
+    void computeVerticesSegSize(std::vector<GC_Seg>& out_segments, const std::vector<bool>& useVertex, float alpha = 0.0f);
+    void removeSmallSegs(const std::vector<GC_Seg>& segments, int minSegSize);
 
     /**
      * @brief Function that returns the next geometry intersected by the ray.
@@ -491,36 +508,34 @@ class DelaunayGraphCut
      */
     void forceTedgesByGradientIJCV(bool fixesSigma, float nPixelSizeBehind);
 
-    int setIsOnSurface();
+    int computeIsOnSurface(std::vector<bool>& vertexIsOnSurface) const;
 
     void addToInfiniteSw(float sW);
 
-    void freeUnwantedFullCells(const Point3d* hexah);
-
-    void reconstructGC(const Point3d* hexah);
-
     void maxflow();
 
     void voteFullEmptyScore(const StaticVector<int>& cams, const std::string& folderName);
 
-    void createDensePointCloud(Point3d hexah[8], const StaticVector<int>& cams, const sfmData::SfMData* sfmData, const FuseParams* depthMapsFuseParams);
+    void createDensePointCloud(const Point3d hexah[8], const StaticVector<int>& cams, const sfmData::SfMData* sfmData, const FuseParams* depthMapsFuseParams);
 
-    void createGraphCut(Point3d hexah[8], const StaticVector<int>& cams, const std::string& folderName, const std::string& tmpCamsPtsFolderName, bool removeSmallSegments);
+    void createGraphCut(const Point3d hexah[8], const StaticVector<int>& cams, const std::string& folderName, const std::string& tmpCamsPtsFolderName, bool removeSmallSegments);
 
     /**
      * @brief Invert full/empty status of cells if they represent a too small group after labelling.
      */
     void invertFullStatusForSmallLabels();
 
-    void graphCutPostProcessing();
+    void graphCutPostProcessing(const Point3d hexah[8], const std::string& folderName);
 
-    void segmentFullOrFree(bool full, StaticVector<int>** inColors, int& nsegments);
+    void segmentFullOrFree(bool full, StaticVector<int>& out_fullSegsColor, int& nsegments);
     int removeBubbles();
     int removeDust(int minSegSize);
     void leaveLargestFullSegmentOnly();
 
     mesh::Mesh* createMesh(bool filterHelperPointsTriangles = true);
     mesh::Mesh* createTetrahedralMesh(bool filter = true, const float& downscaleFactor = 0.95f, const std::function<float(const GC_cellInfo&)> getScore = [](const GC_cellInfo& c) { return c.emptinessScore; }) const;
+
+    void displayCellsStats() const;
     void exportDebugMesh(const std::string& filename, const Point3d& fromPt, const Point3d& toPt);
     void exportFullScoreMeshs();
     void exportBackPropagationMesh(const std::string& filename, std::vector<GeometryIntersection>& intersectedGeom, const Point3d& fromPt, const Point3d& toPt);
@@ -529,6 +544,10 @@ class DelaunayGraphCut
 
 
 std::ostream& operator<<(std::ostream& stream, const DelaunayGraphCut::EGeometryType type);
+std::ostream& operator<<(std::ostream& stream, const DelaunayGraphCut::Facet& facet);
+std::ostream& operator<<(std::ostream& stream, const DelaunayGraphCut::Edge& edge);
+std::ostream& operator<<(std::ostream& stream, const DelaunayGraphCut::GeometryIntersection& intersection);
+std::ostream& operator<<(std::ostream& stream, const DelaunayGraphCut::GeometriesCount& count);
 
 } // namespace fuseCut
 } // namespace aliceVision

src/aliceVision/fuseCut/DelaunayGraphCut_test.cpp

@@ -64,13 +64,15 @@ BOOST_AUTO_TEST_CASE(fuseCut_delaunayGraphCut)
     for (size_t i = 0; i < delaunayGC._verticesCoords.size(); i++)
         ALICEVISION_LOG_TRACE("[" << i << "]: " << delaunayGC._verticesCoords[i].x << ", " << delaunayGC._verticesCoords[i].y << ", " << delaunayGC._verticesCoords[i].z);
 
-    // delaunayGC.createGraphCut(&hexah[0], cams, tempDirPath + "/", tempDirPath + "/SpaceCamsTracks/", false);
+    delaunayGC.createGraphCut(&hexah[0], cams, tempDirPath + "/", tempDirPath + "/SpaceCamsTracks/", false);
+    /*
     delaunayGC.initVertices();
     delaunayGC.computeDelaunay();
     delaunayGC.displayStatistics();
     delaunayGC.computeVerticesSegSize(true, 0.0f);
     delaunayGC.voteFullEmptyScore(cams, tempDirPath);
     delaunayGC.reconstructGC(&hexah[0]);
+    */
 
     ALICEVISION_LOG_TRACE("CreateGraphCut Done.");
 }
@@ -171,4 +173,4 @@ SfMData generateSfm(const NViewDatasetConfigurator& config, const size_t size, c
     }
 
     return sfm_data;
-}
+}

src/aliceVision/fuseCut/delaunayGraphCutTypes.hpp

@@ -55,14 +55,17 @@ struct GC_cellInfo
     }
 };
 
+struct GC_Seg
+{
+    int segSize = 0;
+    int segId = -1;
+};
+
 struct GC_vertexInfo
 {
     float pixSize = 0.0f;
     /// Number of cameras which have contributed to the refinement of the vertex position, so nrc >= cams.size().
     int nrc = 0;
-    int segSize = 0;
-    int segId = -1;
-    bool isOnSurface = false;
     /// All cameras having a visibility of this vertex. Some of them may not have contributed to the vertex position
     StaticVector<int> cams;
 
@@ -86,9 +89,6 @@ struct GC_vertexInfo
     {
         fwrite(&pixSize, sizeof(float), 1, f);
         fwrite(&nrc, sizeof(int), 1, f);
-        fwrite(&segSize, sizeof(int), 1, f);
-        fwrite(&segId, sizeof(int), 1, f);
-        fwrite(&isOnSurface, sizeof(bool), 1, f);
         int n = cams.size();
         fwrite(&n, sizeof(int), 1, f);
         if(n > 0)
@@ -101,9 +101,6 @@ struct GC_vertexInfo
     {
         fread(&pixSize, sizeof(float), 1, f);
         fread(&nrc, sizeof(int), 1, f);
-        fread(&segSize, sizeof(int), 1, f);
-        fread(&segId, sizeof(int), 1, f);
-        fread(&isOnSurface, sizeof(bool), 1, f);
         int n;
         fread(&n, sizeof(int), 1, f);
         if(n > 0)

src/aliceVision/mesh/MeshAnalyze.cpp

@@ -150,7 +150,7 @@ void MeshAnalyze::getVertexPrincipalCurvatures(double Kh, double Kg, double& K1,
     K2 = Kh - temp;
 }
 
-bool MeshAnalyze::applyLaplacianOperator(int ptId, StaticVector<Point3d>& ptsToApplyLaplacianOp, Point3d& ln)
+bool MeshAnalyze::applyLaplacianOperator(int ptId, const StaticVector<Point3d>& ptsToApplyLaplacianOp, Point3d& ln)
 {
     StaticVector<int>& ptNeighPtsOrdered = ptsNeighPtsOrdered[ptId];
     if(ptNeighPtsOrdered.empty())
@@ -181,12 +181,6 @@ bool MeshAnalyze::applyLaplacianOperator(int ptId, StaticVector<Point3d>& ptsToA
         return false;
     }
 
-    if(std::isnan(d) || std::isnan(n.x) || std::isnan(n.y) || std::isnan(n.z)) // check if is not NaN
-    {
-        ALICEVISION_LOG_WARNING("MeshAnalyze::applyLaplacianOperator: nan");
-        return false;
-    }
-
     return true;
 }
 
@@ -198,9 +192,10 @@ bool MeshAnalyze::getLaplacianSmoothingVector(int ptId, Point3d& ln)
 }
 
 // kobbelt kampagna 98 Interactive Multi-Resolution Modeling on Arbitrary Meshes
-// page 5 - U1 - laplacian is obtained when apply to origina pts , U2 - bi-laplacian is obtained when apply to laplacian
-// pts
-bool MeshAnalyze::getBiLaplacianSmoothingVector(int ptId, StaticVector<Point3d>& ptsLaplacian, Point3d& tp)
+// page 5:
+// U1 - laplacian is obtained when apply to original points,
+// U2 - bi-laplacian is obtained when apply to laplacian points
+bool MeshAnalyze::getBiLaplacianSmoothingVector(int ptId, const StaticVector<Point3d>& ptsLaplacian, Point3d& tp)
 {
     if(applyLaplacianOperator(ptId, ptsLaplacian, tp))
     {

src/aliceVision/mesh/MeshAnalyze.hpp

@@ -25,9 +25,9 @@ class MeshAnalyze : public MeshClean
     int getVertexIdInTriangleForPtId(int ptId, int triId);
     bool getVertexMeanCurvatureNormal(int ptId, Point3d& Kh);
     void getVertexPrincipalCurvatures(double Kh, double Kg, double& K1, double& K2);
-    bool applyLaplacianOperator(int ptId, StaticVector<Point3d>& ptsToApplyLaplacianOp, Point3d& ln);
+    bool applyLaplacianOperator(int ptId, const StaticVector<Point3d>& ptsToApplyLaplacianOp, Point3d& ln);
     bool getLaplacianSmoothingVector(int ptId, Point3d& ln);
-    bool getBiLaplacianSmoothingVector(int ptId, StaticVector<Point3d>& ptsLaplacian, Point3d& tp);
+    bool getBiLaplacianSmoothingVector(int ptId, const StaticVector<Point3d>& ptsLaplacian, Point3d& tp);
     bool getMeanCurvAndLaplacianSmoothing(int ptId, Point3d& F, float epsilon);
     bool getVertexSurfaceNormal(int ptId, Point3d& N);
 };

src/aliceVision/mvsData/Point2d.hpp

@@ -7,6 +7,7 @@
 #pragma once
 
 #include <cmath>
+#include <ostream>
 
 namespace aliceVision {
 
@@ -90,4 +91,10 @@ inline double dot(const Point2d& p1, const Point2d& p2)
     return p1.x * p2.x + p1.y * p2.y;
 }
 
+inline std::ostream& operator<<(std::ostream& stream, const Point2d& p)
+{
+    stream << p.x << "," << p.y;
+    return stream;
+}
+
 } // namespace aliceVision

src/aliceVision/mvsData/Point3d.hpp

@@ -8,6 +8,7 @@
 
 #include <string>
 #include <cmath>
+#include <ostream>
 
 namespace aliceVision {
 
@@ -107,7 +108,7 @@ class Point3d
 
     friend double dot(const Point3d& p1, const Point3d& p2);
     friend Point3d cross(const Point3d& a, const Point3d& b);
-    friend Point3d proj(Point3d& e, Point3d& a);
+    friend Point3d proj(const Point3d& e, const Point3d& a);
 };
 
 inline double dist(const Point3d& p1, const Point3d& p2)
@@ -130,9 +131,15 @@ inline Point3d cross(const Point3d& a, const Point3d& b)
     return vc;
 }
 
-inline Point3d proj(Point3d& e, Point3d& a)
+inline Point3d proj(const Point3d& e, const Point3d& a)
 {
     return e * (dot(e, a) / dot(e, e));
 }
 
+inline std::ostream& operator<<(std::ostream& stream, const Point3d& p)
+{
+    stream << p.x << "," << p.y << "," << p.z;
+    return stream;
+}
+
 } // namespace aliceVision

src/aliceVision/mvsData/Point4d.hpp

@@ -6,6 +6,8 @@
 
 #pragma once
 
+#include <ostream>
+
 namespace aliceVision {
 
 class Point4d
@@ -99,4 +101,10 @@ class Point4d
     }
 };
 
+inline std::ostream& operator<<(std::ostream& stream, const Point4d& p)
+{
+    stream << p.x << "," << p.y << "," << p.z << "," << p.w;
+    return stream;
+}
+
 } // namespace aliceVision