From 9f285e8b6c440587d35e97f59cb0439c5bbce25e Mon Sep 17 00:00:00 2001
From: Enguerrand DE SMET <desmet.enguerrand@gmail.com>
Date: Mon, 6 Jul 2020 17:24:59 +0200
Subject: [PATCH 01/28] [fuseCut] refactoring neighbourCellToTheCamOnTheRay

refactoring nearestNeighCellToTheCamOnTheRay and farestNeighCellToTheCamOnTheRay  functions into a common function with a boolean parameter.
---
 src/aliceVision/fuseCut/DelaunayGraphCut.cpp | 12 +++---
 src/aliceVision/fuseCut/DelaunayGraphCut.hpp | 41 ++++++++------------
 2 files changed, 22 insertions(+), 31 deletions(-)

diff --git a/src/aliceVision/fuseCut/DelaunayGraphCut.cpp b/src/aliceVision/fuseCut/DelaunayGraphCut.cpp
index 999d12eace..bf408b209d 100644
--- a/src/aliceVision/fuseCut/DelaunayGraphCut.cpp
+++ b/src/aliceVision/fuseCut/DelaunayGraphCut.cpp
@@ -1734,7 +1734,7 @@ void DelaunayGraphCut::fillGraphPartPtRc(int& out_nstepsFront, int& out_nstepsBe
             Facet f1, f2;
             Point3d lpi;
             // find cell which is nearest to the cam and which is intersected with cam-p ray
-            if(!nearestNeighCellToTheCamOnTheRay(mp->CArr[cam], p, ci, f1, f2, lpi))
+            if(!nearestNeighbourCellToTheCamOnTheRay(mp->CArr[cam], p, ci, f1, f2, lpi))
             {
                 ok = false;
             }
@@ -1797,7 +1797,7 @@ void DelaunayGraphCut::fillGraphPartPtRc(int& out_nstepsFront, int& out_nstepsBe
             Point3d pold = p;
             Point3d lpi;
             // find cell which is farest to the cam and which intersect cam-p ray
-            if((!farestNeighCellToTheCamOnTheRay(mp->CArr[cam], p, ci, f1, f2, lpi)) ||
+            if((!farestNeighbourCellToTheCamOnTheRay(mp->CArr[cam], p, ci, f1, f2, lpi)) ||
                ((po - pold).size() >= maxDist) || (!allPoints))
             {
                 ok = false;
@@ -1897,7 +1897,7 @@ void DelaunayGraphCut::forceTedgesByGradientCVPR11(bool fixesSigma, float nPixel
                     // find cell which is farest to the cam and which is
                     // intersected with cam-p
                     // ray
-                    if((!farestNeighCellToTheCamOnTheRay(mp->CArr[cam], p, ci, ff1, ff2, lpi)) ||
+                    if((!farestNeighbourCellToTheCamOnTheRay(mp->CArr[cam], p, ci, ff1, ff2, lpi)) ||
                        ((po - pold).size() >= maxDist))
                     {
                         ok = false;
@@ -2031,7 +2031,7 @@ void DelaunayGraphCut::forceTedgesByGradientIJCV(bool fixesSigma, float nPixelSi
                     Point3d lpi;
                     // find cell which is nearest to the cam and which intersect cam-p ray
                     if(((p - po).size() > (nsigmaJumpPart + nsigmaFrontSilentPart) * maxDist) || // (2 + 2) * sigma
-                       (!nearestNeighCellToTheCamOnTheRay(mp->CArr[cam], p, ci, f1, f2, lpi)))
+                       (!nearestNeighbourCellToTheCamOnTheRay(mp->CArr[cam], p, ci, f1, f2, lpi)))
                     {
                         ok = false;
                     }
@@ -2065,7 +2065,7 @@ void DelaunayGraphCut::forceTedgesByGradientIJCV(bool fixesSigma, float nPixelSi
                     Point3d lpi;
                     // find cell which is farest to the cam and which intersect cam-p ray
                     if(((p - po).size() > nsigmaBackSilentPart * maxDist) || // (p-po).size() > 2 * sigma
-                       (!farestNeighCellToTheCamOnTheRay(mp->CArr[cam], p, ci, f1, f2, lpi)))
+                       (!farestNeighbourCellToTheCamOnTheRay(mp->CArr[cam], p, ci, f1, f2, lpi)))
                     {
                         ok = false;
                     }
@@ -2209,7 +2209,7 @@ void DelaunayGraphCut::updateGraphFromTmpPtsCamsHexahRC(int rc, Point3d hexah[8]
                     Point3d lpi;
                     // find cell which is nearest to the pont and which is
                     // intersected with point-p ray
-                    if(!farestNeighCellToTheCamOnTheRay(camBehind, p, tmp_ci, f1, f2, lpi))
+                    if(!farestNeighbourCellToTheCamOnTheRay(camBehind, p, tmp_ci, f1, f2, lpi))
                     {
                         ok = false;
                     }
diff --git a/src/aliceVision/fuseCut/DelaunayGraphCut.hpp b/src/aliceVision/fuseCut/DelaunayGraphCut.hpp
index ba94f92053..90a21f3cfe 100644
--- a/src/aliceVision/fuseCut/DelaunayGraphCut.hpp
+++ b/src/aliceVision/fuseCut/DelaunayGraphCut.hpp
@@ -345,44 +345,35 @@ class DelaunayGraphCut
     mesh::Mesh* createMesh(bool filterHelperPointsTriangles = true);
 };
 
-
-inline bool DelaunayGraphCut::nearestNeighCellToTheCamOnTheRay(const Point3d& camC, Point3d& out_p, int tetrahedron, Facet& out_f1,
-                                                      Facet& out_f2, Point3d& out_lpi) const
-{
-    out_f1.cellIndex = GEO::NO_CELL;
-    out_f1.localVertexIndex = GEO::NO_VERTEX;
-    out_f2.cellIndex = GEO::NO_CELL;
-    out_f2.localVertexIndex = GEO::NO_VERTEX;
-    out_lpi = out_p;
-
-    if(rayCellIntersection(camC, out_p, tetrahedron, out_f1, true, out_lpi) == true)
-    {
-        out_f2 = mirrorFacet(out_f1);
-        out_p = out_lpi;
-        return true;
-    }
-    return false;
-}
-
-inline bool DelaunayGraphCut::farestNeighCellToTheCamOnTheRay(Point3d& camC, Point3d& p, int tetrahedron, Facet& f1,
-                                                     Facet& f2, Point3d& lpi) const
+inline bool DelaunayGraphCut::neighbourCellToTheCamOnTheRay(const Point3d& camC, Point3d& p, int tetrahedron, bool nearestFarest, 
+                                                    Facet& f1, Facet& f2, Point3d& outIntersectPt) const
 {
     f1.cellIndex = GEO::NO_CELL;
     f1.localVertexIndex = GEO::NO_VERTEX;
     f2.cellIndex = GEO::NO_CELL;
     f2.localVertexIndex = GEO::NO_VERTEX;
-    lpi = p;
+    outIntersectPt = p;
 
-    if(rayCellIntersection(camC, p, tetrahedron, f1, false, lpi) == true)
+    if(rayCellIntersection(camC, p, tetrahedron, f1, nearestFarest, outIntersectPt))
     {
         f2 = mirrorFacet(f1);
-        p = lpi;
+        p = outIntersectPt;
         return true;
     }
-
     return false;
 }
 
+inline bool DelaunayGraphCut::nearestNeighbourCellToTheCamOnTheRay(const Point3d& camC, Point3d& p, int tetrahedron, Facet& f1, Facet& f2, 
+                                                            Point3d& outIntersectPt) const
+{
+    return neighbourCellToTheCamOnTheRay(camC, p, tetrahedron, true, f1, f2, outIntersectPt);
+}
+inline bool DelaunayGraphCut::farestNeighbourCellToTheCamOnTheRay(const Point3d& camC, Point3d& p, int tetrahedron, Facet& f1,  Facet& f2, 
+                                                            Point3d& outIntersectPt) const
+{
+    return neighbourCellToTheCamOnTheRay(camC, p, tetrahedron, false, f1, f2, outIntersectPt);
+}
+
 inline DelaunayGraphCut::Facet DelaunayGraphCut::getFacetInFrontVertexOnTheRayToTheCam(int vertexIndex,
                                                                                int cam) const
 {

From 2f90040d95d3b7a84651ea6147636e8e6e98ed18 Mon Sep 17 00:00:00 2001
From: Enguerrand DE SMET <desmet.enguerrand@gmail.com>
Date: Mon, 6 Jul 2020 18:07:28 +0200
Subject: [PATCH 02/28] [fuseCut] renaming and use appropriate variable type
 GEO::index_t

---
 src/aliceVision/fuseCut/DelaunayGraphCut.cpp | 74 ++++++++++----------
 src/aliceVision/fuseCut/DelaunayGraphCut.hpp | 18 +++--
 2 files changed, 48 insertions(+), 44 deletions(-)

diff --git a/src/aliceVision/fuseCut/DelaunayGraphCut.cpp b/src/aliceVision/fuseCut/DelaunayGraphCut.cpp
index bf408b209d..4f3ba0108c 100644
--- a/src/aliceVision/fuseCut/DelaunayGraphCut.cpp
+++ b/src/aliceVision/fuseCut/DelaunayGraphCut.cpp
@@ -1308,28 +1308,28 @@ void DelaunayGraphCut::removeSmallSegs(int minSegSize)
     initVertices();
 }
 
-bool DelaunayGraphCut::rayCellIntersection(const Point3d& camC, const Point3d& p, int tetrahedron, Facet& out_facet,
-                                        bool nearestFarest, Point3d& out_nlpi) const
+bool DelaunayGraphCut::rayCellIntersection(const Point3d& camCenter, const Point3d& p, GEO::index_t tetrahedronIndex, Facet& outFacet,
+                                           bool nearestFarest, Point3d& outIntersectPt) const
 {
-    out_nlpi = p; // important
-    out_facet.cellIndex = -1;
-    out_facet.localVertexIndex = -1;
+    outIntersectPt = p; // important
+    outFacet.cellIndex = GEO::NO_CELL;
+    outFacet.localVertexIndex = GEO::NO_VERTEX;
 
-    if(isInfiniteCell(tetrahedron))
+    if(isInfiniteCell(tetrahedronIndex))
     {
         return false;
     }
 
     Point3d linePoint = p;
-    Point3d lineVect = (camC - p).normalize();
-    double mind = (camC - p).size();
+    Point3d lineVect = (camCenter - p).normalize();
+    double currentDist = (camCenter - p).size();
 
     // go thru all triangles
     bool existsTriOnRay = false;
-    const Point3d* A = &(_verticesCoords[_tetrahedralization->cell_vertex(tetrahedron, 0)]);
-    const Point3d* B = &(_verticesCoords[_tetrahedralization->cell_vertex(tetrahedron, 1)]);
-    const Point3d* C = &(_verticesCoords[_tetrahedralization->cell_vertex(tetrahedron, 2)]);
-    const Point3d* D = &(_verticesCoords[_tetrahedralization->cell_vertex(tetrahedron, 3)]);
+    const Point3d* A = &(_verticesCoords[_tetrahedralization->cell_vertex(tetrahedronIndex, 0)]);
+    const Point3d* B = &(_verticesCoords[_tetrahedralization->cell_vertex(tetrahedronIndex, 1)]);
+    const Point3d* C = &(_verticesCoords[_tetrahedralization->cell_vertex(tetrahedronIndex, 2)]);
+    const Point3d* D = &(_verticesCoords[_tetrahedralization->cell_vertex(tetrahedronIndex, 3)]);
 
     // All the facets of the tetrahedron
     std::array<std::array<const Point3d*, 3>, 4> facets {{
@@ -1339,34 +1339,34 @@ bool DelaunayGraphCut::rayCellIntersection(const Point3d& camC, const Point3d& p
         {A, B, C}  // opposite vertex D, index 3
     }};
     int oppositeVertexIndex = -1;
-    Point3d lpi;
-    double dist;
+    Point3d intersectPt;
+    double intersectDist;
     // Test all facets of the tetrahedron
     for(int i = 0; i < 4; ++i)
     {
         const auto& facet= facets[i];
-        if(isLineInTriangle(&lpi, facet[0], facet[1], facet[2], &linePoint, &lineVect))
+        if(isLineInTriangle(&intersectPt, facet[0], facet[1], facet[2], &linePoint, &lineVect))
         {
-            dist = (camC - lpi).size();
+            intersectDist = (camCenter - intersectPt).size();
             if(nearestFarest)
             {
-                if(dist < mind) // between the camera and the point
+                if(intersectDist < currentDist) // between the camera and the point
                 {
                     oppositeVertexIndex = i;
                     existsTriOnRay = true;
-                    mind = dist;
-                    out_nlpi = lpi;
+                    currentDist = intersectDist;
+                    outIntersectPt = intersectPt;
                     //break;
                 }
             }
             else
             {
-                if(dist > mind) // behind the point (from the camera)
+                if(intersectDist > currentDist) // behind the point (from the camera)
                 {
                     oppositeVertexIndex = i;
                     existsTriOnRay = true;
-                    mind = dist;
-                    out_nlpi = lpi;
+                    currentDist = intersectDist;
+                    outIntersectPt = intersectPt;
                     //break;
                 }
             }
@@ -1376,8 +1376,8 @@ bool DelaunayGraphCut::rayCellIntersection(const Point3d& camC, const Point3d& p
     if(!existsTriOnRay)
         return false;
 
-    out_facet.cellIndex = tetrahedron; // cell handle
-    out_facet.localVertexIndex = oppositeVertexIndex; // vertex index that is not in the intersected facet
+    outFacet.cellIndex = tetrahedronIndex;           // cell handle
+    outFacet.localVertexIndex = oppositeVertexIndex; // vertex index that is not in the intersected facet
     return true;
 }
 
@@ -1388,8 +1388,8 @@ DelaunayGraphCut::Facet DelaunayGraphCut::getFacetInFrontVertexOnTheRayToThePoin
 
     double minDist = (ptt - p).size(); // initialize minDist to the distance from 3d point p to camera center c
     Facet nearestFacet;
-    nearestFacet.cellIndex = -1;
-    nearestFacet.localVertexIndex = -1;
+    nearestFacet.cellIndex = GEO::NO_CELL;
+    nearestFacet.localVertexIndex = GEO::NO_VERTEX;
 
     for(int k = 0; true; ++k)
     {
@@ -1431,28 +1431,28 @@ DelaunayGraphCut::Facet DelaunayGraphCut::getFacetBehindVertexOnTheRayToTheCam(V
         if(isInfiniteCell(adjCellIndex))
             continue;
 
-        Facet f1;
-        Point3d lpi;
-        if(rayCellIntersection(mp->CArr[cam], p, adjCellIndex, f1, false, lpi) == true)
+        Facet outFacet;
+        Point3d intersectionPoint;
+        if(rayCellIntersection(mp->CArr[cam], p, adjCellIndex, outFacet, false, intersectionPoint))
         {
             // if it is after the point p (along the axis from the camera)
-            if((mp->CArr[cam] - lpi).size() > maxDist)
+            if((mp->CArr[cam] - intersectionPoint).size() > maxDist)
             {
-                farestFacet = f1;
-                maxDist = (mp->CArr[cam] - lpi).size();
+                farestFacet = outFacet;
+                maxDist = (mp->CArr[cam] - intersectionPoint).size();
             }
         }
     }
     return farestFacet;
 }
 
-int DelaunayGraphCut::getFirstCellOnTheRayFromCamToThePoint(int cam, Point3d& p, Point3d& lpi) const
+GEO::index_t DelaunayGraphCut::getFirstCellOnTheRayFromCamToThePoint(int cam, Point3d& p, Point3d& intersectPoint) const
 {
     int cam_vi = _camsVertexes[cam];
     Point3d camBehind = mp->CArr[cam] + (mp->CArr[cam] - p);
     Point3d dir = (p - mp->CArr[cam]).normalize();
     float maxdist = 0.0f;
-    int farestCell = -1;
+    GEO::index_t farestCell = GEO::NO_CELL;
 
     for(int k = 0; true; ++k)
     {
@@ -1462,10 +1462,10 @@ int DelaunayGraphCut::getFirstCellOnTheRayFromCamToThePoint(int cam, Point3d& p,
         if(isInfiniteCell(adjCellIndex))
             continue;
 
-        Facet f1;
-        if(rayCellIntersection(camBehind, mp->CArr[cam], adjCellIndex, f1, false, lpi) == true)
+        Facet outFacet;
+        if(rayCellIntersection(camBehind, mp->CArr[cam], adjCellIndex, outFacet, false, intersectPoint))
         {
-            float dist = orientedPointPlaneDistance(lpi, camBehind, dir);
+            float dist = orientedPointPlaneDistance(intersectPoint, camBehind, dir);
             if(dist > maxdist)
             {
                 maxdist = dist;
diff --git a/src/aliceVision/fuseCut/DelaunayGraphCut.hpp b/src/aliceVision/fuseCut/DelaunayGraphCut.hpp
index 90a21f3cfe..08c2e258fc 100644
--- a/src/aliceVision/fuseCut/DelaunayGraphCut.hpp
+++ b/src/aliceVision/fuseCut/DelaunayGraphCut.hpp
@@ -267,16 +267,20 @@ class DelaunayGraphCut
     void computeVerticesSegSize(bool allPoints, float alpha = 0.0f);
     void removeSmallSegs(int minSegSize);
 
-    bool rayCellIntersection(const Point3d &camC, const Point3d &p, int c, Facet& out_facet, bool nearestFarest,
-                            Point3d& out_nlpi) const;
-    inline bool nearestNeighCellToTheCamOnTheRay(const Point3d &camC, Point3d& out_p, int tetrahedron, Facet &out_f1, Facet &f2,
-                                          Point3d& out_lpi) const;
-    inline bool farestNeighCellToTheCamOnTheRay(Point3d& camC, Point3d& p, int tetrahedron, Facet &f1, Facet &f2,
-                                         Point3d& lpi) const;
+    bool rayCellIntersection(const Point3d& camCenter, const Point3d& p, GEO::index_t tetrahedronIndex, Facet& outFacet,
+                             bool nearestFarest, Point3d& outIntersectPt) const;
+
+    inline bool neighbourCellToTheCamOnTheRay(const Point3d& camC, Point3d& p, int tetrahedron, bool nearestFarest, 
+                                            Facet& f1, Facet& f2, Point3d& outIntersectPt) const;
+    inline bool nearestNeighbourCellToTheCamOnTheRay(const Point3d& camC, Point3d& p, int tetrahedron, Facet& f1, Facet& f2,
+                                                     Point3d& outIntersectPt) const;
+    inline bool farestNeighbourCellToTheCamOnTheRay(const Point3d& camC, Point3d& p, int tetrahedron, Facet& f1, Facet& f2,
+                                                    Point3d& outIntersectPt) const;
+
     inline Facet getFacetInFrontVertexOnTheRayToTheCam(int vertexIndex, int cam) const;
     Facet getFacetInFrontVertexOnTheRayToThePoint3d(VertexIndex vi, Point3d& ptt) const;
     Facet getFacetBehindVertexOnTheRayToTheCam(VertexIndex vi, int cam) const;
-    int getFirstCellOnTheRayFromCamToThePoint(int cam, Point3d& p, Point3d& lpi) const;
+    GEO::index_t getFirstCellOnTheRayFromCamToThePoint(int cam, Point3d& p, Point3d& lpi) const;
 
     float distFcn(float maxDist, float dist, float distFcnHeight) const;
 

From eb6e5203fa019e2cbbe7b8c56952847964ec747d Mon Sep 17 00:00:00 2001
From: Enguerrand DE SMET <desmet.enguerrand@gmail.com>
Date: Mon, 6 Jul 2020 18:48:15 +0200
Subject: [PATCH 03/28] [fuseCut] Use a constant input and move the assignment
 to a higher level

- remove unnecessary double assignment :
    f1.cellIndex = GEO::NO_CELL;
    f1.localVertexIndex = GEO::NO_VERTEX;
- move 'p = outIntersectPt' assigment in a higher level in loop with 'p = lpi;' each time the functions nearestNeighbourCellToTheCamOnTheRay and farestNeighbourCellToTheCamOnTheRay are used
---
 src/aliceVision/fuseCut/DelaunayGraphCut.cpp |  6 ++++++
 src/aliceVision/fuseCut/DelaunayGraphCut.hpp | 15 ++++++---------
 2 files changed, 12 insertions(+), 9 deletions(-)

diff --git a/src/aliceVision/fuseCut/DelaunayGraphCut.cpp b/src/aliceVision/fuseCut/DelaunayGraphCut.cpp
index 4f3ba0108c..47dc9e9f2f 100644
--- a/src/aliceVision/fuseCut/DelaunayGraphCut.cpp
+++ b/src/aliceVision/fuseCut/DelaunayGraphCut.cpp
@@ -1750,6 +1750,7 @@ void DelaunayGraphCut::fillGraphPartPtRc(int& out_nstepsFront, int& out_nstepsBe
                 if(f2.cellIndex == GEO::NO_CELL)
                     ok = false;
                 ci = f2.cellIndex;
+                p = lpi;
                 lastFinite = f2.cellIndex;
             }
         }
@@ -1822,6 +1823,7 @@ void DelaunayGraphCut::fillGraphPartPtRc(int& out_nstepsFront, int& out_nstepsBe
                     _cellsAttr[f2.cellIndex].gEdgeVisWeight[f2.localVertexIndex] += weight * dist;
                 }
                 ci = f2.cellIndex;
+                p = lpi;
             }
         }
 
@@ -1907,6 +1909,7 @@ void DelaunayGraphCut::forceTedgesByGradientCVPR11(bool fixesSigma, float nPixel
                         if(ff2.cellIndex == GEO::NO_CELL)
                             ok = false;
                         ci = ff2.cellIndex;
+                        p = lpi;
                         f1 = ff1;
                         f2 = ff2;
                     }
@@ -2040,6 +2043,7 @@ void DelaunayGraphCut::forceTedgesByGradientIJCV(bool fixesSigma, float nPixelSi
                         if(f2.cellIndex == GEO::NO_CELL)
                             ok = false;
                         ci = f2.cellIndex;
+                        p = lpi;
                     }
                 }
             }
@@ -2074,6 +2078,7 @@ void DelaunayGraphCut::forceTedgesByGradientIJCV(bool fixesSigma, float nPixelSi
                         if(f2.cellIndex == GEO::NO_CELL)
                             ok = false;
                         ci = f2.cellIndex;
+                        p = lpi;
                     }
                 }
 
@@ -2225,6 +2230,7 @@ void DelaunayGraphCut::updateGraphFromTmpPtsCamsHexahRC(int rc, Point3d hexah[8]
                             _cellsAttr[f2.cellIndex].gEdgeVisWeight[f2.localVertexIndex] += weight;
                         }
                         tmp_ci = f2.cellIndex;
+                        p = lpi;
                         ++nwup;
                     }
                 }
diff --git a/src/aliceVision/fuseCut/DelaunayGraphCut.hpp b/src/aliceVision/fuseCut/DelaunayGraphCut.hpp
index 08c2e258fc..17fa047ad9 100644
--- a/src/aliceVision/fuseCut/DelaunayGraphCut.hpp
+++ b/src/aliceVision/fuseCut/DelaunayGraphCut.hpp
@@ -270,11 +270,11 @@ class DelaunayGraphCut
     bool rayCellIntersection(const Point3d& camCenter, const Point3d& p, GEO::index_t tetrahedronIndex, Facet& outFacet,
                              bool nearestFarest, Point3d& outIntersectPt) const;
 
-    inline bool neighbourCellToTheCamOnTheRay(const Point3d& camC, Point3d& p, int tetrahedron, bool nearestFarest, 
+    inline bool neighbourCellToTheCamOnTheRay(const Point3d& camC, const Point3d& p, int tetrahedron, bool nearestFarest, 
                                             Facet& f1, Facet& f2, Point3d& outIntersectPt) const;
-    inline bool nearestNeighbourCellToTheCamOnTheRay(const Point3d& camC, Point3d& p, int tetrahedron, Facet& f1, Facet& f2,
+    inline bool nearestNeighbourCellToTheCamOnTheRay(const Point3d& camC, const Point3d& p, int tetrahedron, Facet& f1, Facet& f2,
                                                      Point3d& outIntersectPt) const;
-    inline bool farestNeighbourCellToTheCamOnTheRay(const Point3d& camC, Point3d& p, int tetrahedron, Facet& f1, Facet& f2,
+    inline bool farestNeighbourCellToTheCamOnTheRay(const Point3d& camC, const Point3d& p, int tetrahedron, Facet& f1, Facet& f2,
                                                     Point3d& outIntersectPt) const;
 
     inline Facet getFacetInFrontVertexOnTheRayToTheCam(int vertexIndex, int cam) const;
@@ -349,11 +349,9 @@ class DelaunayGraphCut
     mesh::Mesh* createMesh(bool filterHelperPointsTriangles = true);
 };
 
-inline bool DelaunayGraphCut::neighbourCellToTheCamOnTheRay(const Point3d& camC, Point3d& p, int tetrahedron, bool nearestFarest, 
+inline bool DelaunayGraphCut::neighbourCellToTheCamOnTheRay(const Point3d& camC, const Point3d& p, int tetrahedron, bool nearestFarest, 
                                                     Facet& f1, Facet& f2, Point3d& outIntersectPt) const
 {
-    f1.cellIndex = GEO::NO_CELL;
-    f1.localVertexIndex = GEO::NO_VERTEX;
     f2.cellIndex = GEO::NO_CELL;
     f2.localVertexIndex = GEO::NO_VERTEX;
     outIntersectPt = p;
@@ -361,18 +359,17 @@ inline bool DelaunayGraphCut::neighbourCellToTheCamOnTheRay(const Point3d& camC,
     if(rayCellIntersection(camC, p, tetrahedron, f1, nearestFarest, outIntersectPt))
     {
         f2 = mirrorFacet(f1);
-        p = outIntersectPt;
         return true;
     }
     return false;
 }
 
-inline bool DelaunayGraphCut::nearestNeighbourCellToTheCamOnTheRay(const Point3d& camC, Point3d& p, int tetrahedron, Facet& f1, Facet& f2, 
+inline bool DelaunayGraphCut::nearestNeighbourCellToTheCamOnTheRay(const Point3d& camC, const Point3d& p, int tetrahedron, Facet& f1, Facet& f2, 
                                                             Point3d& outIntersectPt) const
 {
     return neighbourCellToTheCamOnTheRay(camC, p, tetrahedron, true, f1, f2, outIntersectPt);
 }
-inline bool DelaunayGraphCut::farestNeighbourCellToTheCamOnTheRay(const Point3d& camC, Point3d& p, int tetrahedron, Facet& f1,  Facet& f2, 
+inline bool DelaunayGraphCut::farestNeighbourCellToTheCamOnTheRay(const Point3d& camC, const Point3d& p, int tetrahedron, Facet& f1,  Facet& f2, 
                                                             Point3d& outIntersectPt) const
 {
     return neighbourCellToTheCamOnTheRay(camC, p, tetrahedron, false, f1, f2, outIntersectPt);

From 244fdfd48c7e710cc3f501736ad90dacb439c3ec Mon Sep 17 00:00:00 2001
From: Enguerrand DE SMET <desmet.enguerrand@gmail.com>
Date: Tue, 7 Jul 2020 11:00:25 +0200
Subject: [PATCH 04/28] [fuseCut] rayCellIntersection: group duplicate codes
 using ternary operator

---
 src/aliceVision/fuseCut/DelaunayGraphCut.cpp | 27 +++++---------------
 1 file changed, 7 insertions(+), 20 deletions(-)

diff --git a/src/aliceVision/fuseCut/DelaunayGraphCut.cpp b/src/aliceVision/fuseCut/DelaunayGraphCut.cpp
index 47dc9e9f2f..ae2fd36259 100644
--- a/src/aliceVision/fuseCut/DelaunayGraphCut.cpp
+++ b/src/aliceVision/fuseCut/DelaunayGraphCut.cpp
@@ -1348,27 +1348,14 @@ bool DelaunayGraphCut::rayCellIntersection(const Point3d& camCenter, const Point
         if(isLineInTriangle(&intersectPt, facet[0], facet[1], facet[2], &linePoint, &lineVect))
         {
             intersectDist = (camCenter - intersectPt).size();
-            if(nearestFarest)
+            // between the camera and the point if nearestFarest == true
+            // behind the point (from the camera) if nearestFarest == false
+            if(nearestFarest ? (intersectDist < currentDist) : (intersectDist > currentDist))
             {
-                if(intersectDist < currentDist) // between the camera and the point
-                {
-                    oppositeVertexIndex = i;
-                    existsTriOnRay = true;
-                    currentDist = intersectDist;
-                    outIntersectPt = intersectPt;
-                    //break;
-                }
-            }
-            else
-            {
-                if(intersectDist > currentDist) // behind the point (from the camera)
-                {
-                    oppositeVertexIndex = i;
-                    existsTriOnRay = true;
-                    currentDist = intersectDist;
-                    outIntersectPt = intersectPt;
-                    //break;
-                }
+                oppositeVertexIndex = i;
+                existsTriOnRay = true;
+                currentDist = intersectDist;
+                outIntersectPt = intersectPt;
             }
         }
     }

From 28e9233d968176910b0d99f82b0ce4965e1cf9b7 Mon Sep 17 00:00:00 2001
From: Enguerrand DE SMET <desmet.enguerrand@gmail.com>
Date: Tue, 7 Jul 2020 17:21:09 +0200
Subject: [PATCH 05/28] [fuseCut]  DelaunayGraphCut: remove unnecessary
 neighbourCellToTheCamOnTheRay function (CRITIC)

Remove the unnecessary neighborCellToTheCamOnTheRay function by using rayCellIntersection directly in loops and use mirrorFacet when it' s necessary to ensure loop recursion and get the next cell.
---
 src/aliceVision/fuseCut/DelaunayGraphCut.cpp | 109 ++++++++++++-------
 src/aliceVision/fuseCut/DelaunayGraphCut.hpp |  33 ------
 2 files changed, 67 insertions(+), 75 deletions(-)

diff --git a/src/aliceVision/fuseCut/DelaunayGraphCut.cpp b/src/aliceVision/fuseCut/DelaunayGraphCut.cpp
index ae2fd36259..1c99864a4b 100644
--- a/src/aliceVision/fuseCut/DelaunayGraphCut.cpp
+++ b/src/aliceVision/fuseCut/DelaunayGraphCut.cpp
@@ -1718,10 +1718,14 @@ void DelaunayGraphCut::fillGraphPartPtRc(int& out_nstepsFront, int& out_nstepsBe
             ++nsteps;
 
             Point3d pold = p;
-            Facet f1, f2;
+            Facet outFacet;
             Point3d lpi;
-            // find cell which is nearest to the cam and which is intersected with cam-p ray
-            if(!nearestNeighbourCellToTheCamOnTheRay(mp->CArr[cam], p, ci, f1, f2, lpi))
+
+            // Intersection with the next facet in the current tetrahedron (ci) in order to find the cell nearest to the
+            // cam which is intersected with cam-p ray
+            // true here mean nearest
+            const bool nearestFarest = true;
+            if(!rayCellIntersection(mp->CArr[cam], p, ci, outFacet, nearestFarest, lpi))
             {
                 ok = false;
             }
@@ -1731,14 +1735,16 @@ void DelaunayGraphCut::fillGraphPartPtRc(int& out_nstepsFront, int& out_nstepsBe
 
                 {
 #pragma OMP_ATOMIC_UPDATE
-                    _cellsAttr[f1.cellIndex].gEdgeVisWeight[f1.localVertexIndex] += weight * dist;
+                    _cellsAttr[outFacet.cellIndex].gEdgeVisWeight[outFacet.localVertexIndex] += weight * dist;
                 }
 
-                if(f2.cellIndex == GEO::NO_CELL)
+                // Take the mirror facet to iterate over the next cell
+                ci = mirrorFacet(outFacet).cellIndex;
+
+                if(ci == GEO::NO_CELL)
                     ok = false;
-                ci = f2.cellIndex;
                 p = lpi;
-                lastFinite = f2.cellIndex;
+                lastFinite = ci;
             }
         }
 
@@ -1754,7 +1760,7 @@ void DelaunayGraphCut::fillGraphPartPtRc(int& out_nstepsFront, int& out_nstepsBe
         out_nstepsBehind = 0;
         // get the tetrahedron next to point p on the ray from c
         Facet f1 = getFacetBehindVertexOnTheRayToTheCam(vertexIndex, cam);
-        Facet f2;
+        Facet outFacet;
 
         CellIndex ci = f1.cellIndex;
         if(ci != GEO::NO_CELL)
@@ -1784,8 +1790,12 @@ void DelaunayGraphCut::fillGraphPartPtRc(int& out_nstepsFront, int& out_nstepsBe
 
             Point3d pold = p;
             Point3d lpi;
-            // find cell which is farest to the cam and which intersect cam-p ray
-            if((!farestNeighbourCellToTheCamOnTheRay(mp->CArr[cam], p, ci, f1, f2, lpi)) ||
+
+            // Intersection with the next facet in the current tetrahedron (ci) in order to find the cell farest to the
+            // cam which is intersected with cam-p ray
+            // False here mean farest
+            const bool nearestFarest = false;
+            if(!rayCellIntersection(mp->CArr[cam], p, ci, outFacet, nearestFarest, lpi) ||
                ((po - pold).size() >= maxDist) || (!allPoints))
             {
                 ok = false;
@@ -1800,16 +1810,18 @@ void DelaunayGraphCut::fillGraphPartPtRc(int& out_nstepsFront, int& out_nstepsBe
 
                 float dist = distFcn(maxDist, (po - pold).size(), distFcnHeight);
 
-                if(f2.cellIndex == GEO::NO_CELL)
+                // Take the mirror facet to iterate over the next cell
+                const Facet mFacet = mirrorFacet(outFacet);
+                ci = mFacet.cellIndex;
+                if(ci == GEO::NO_CELL)
                 {
                     ok = false;
                 }
                 else
                 {
 #pragma OMP_ATOMIC_UPDATE
-                    _cellsAttr[f2.cellIndex].gEdgeVisWeight[f2.localVertexIndex] += weight * dist;
+                    _cellsAttr[ci].gEdgeVisWeight[mFacet.localVertexIndex] += weight * dist;
                 }
-                ci = f2.cellIndex;
                 p = lpi;
             }
         }
@@ -1861,14 +1873,12 @@ void DelaunayGraphCut::forceTedgesByGradientCVPR11(bool fixesSigma, float nPixel
             Facet fFirst = getFacetInFrontVertexOnTheRayToTheCam(vi, cam);
 
             // get the tetrahedron next to point p on the ray from c
-            Facet f2;
             Facet f1 = getFacetBehindVertexOnTheRayToTheCam(vi, cam);
 
             if((fFirst.cellIndex != GEO::NO_CELL) && (f1.cellIndex != GEO::NO_CELL) && (!isInfiniteCell(f1.cellIndex)))
             {
                 float eFirst = _cellsAttr[fFirst.cellIndex].out;
 
-                f2 = mirrorFacet(f1);
                 CellIndex ci = f1.cellIndex;
                 Point3d p = po; // HAS TO BE HERE !!!
                 float maxDist = nPixelSizeBehind * mp->getCamPixelSize(p, cam);
@@ -1882,29 +1892,30 @@ void DelaunayGraphCut::forceTedgesByGradientCVPR11(bool fixesSigma, float nPixel
                 {
                     Point3d pold = p;
                     Point3d lpi;
-                    Facet ff1, ff2;
-                    // find cell which is farest to the cam and which is
-                    // intersected with cam-p
-                    // ray
-                    if((!farestNeighbourCellToTheCamOnTheRay(mp->CArr[cam], p, ci, ff1, ff2, lpi)) ||
+                    Facet outFacet;
+
+                    // Intersection with the next facet in the current tetrahedron (ci) in order to find the cell farest
+                    // to the cam which is intersected with cam-p ray
+                    // False here mean farest
+                    const bool nearestFarest = false;
+                    if(!rayCellIntersection(mp->CArr[cam], p, ci, outFacet, nearestFarest, lpi) ||
                        ((po - pold).size() >= maxDist))
                     {
                         ok = false;
                     }
                     else
                     {
-                        if(ff2.cellIndex == GEO::NO_CELL)
+                        // Take the mirror facet to iterate over the next cell
+                        ci = mirrorFacet(outFacet).cellIndex;
+                        if(ci == GEO::NO_CELL)
                             ok = false;
-                        ci = ff2.cellIndex;
                         p = lpi;
-                        f1 = ff1;
-                        f2 = ff2;
                     }
                 }
 
                 if(ci != GEO::NO_CELL)
                 {
-                    float eLast = _cellsAttr[f2.cellIndex].out;
+                    float eLast = _cellsAttr[ci].out;
                     if((eFirst > eLast) && (eFirst < beta) && (eLast / eFirst < delta))
                     {
 #pragma OMP_ATOMIC_UPDATE
@@ -2017,19 +2028,23 @@ void DelaunayGraphCut::forceTedgesByGradientIJCV(bool fixesSigma, float nPixelSi
                         maxSilent = std::max(maxSilent, c.out);
                     }
 
-                    Facet f1, f2;
+                    Facet outFacet;
                     Point3d lpi;
-                    // find cell which is nearest to the cam and which intersect cam-p ray
+                    // Intersection with the next facet in the current tetrahedron (ci) in order to find the cell nearest
+                    // to the cam which is intersected with cam-p ray
+                    // True here mean nearest
+                    const bool nearestFarest = true;
                     if(((p - po).size() > (nsigmaJumpPart + nsigmaFrontSilentPart) * maxDist) || // (2 + 2) * sigma
-                       (!nearestNeighbourCellToTheCamOnTheRay(mp->CArr[cam], p, ci, f1, f2, lpi)))
+                       !rayCellIntersection(mp->CArr[cam], p, ci, outFacet, nearestFarest, lpi))
                     {
                         ok = false;
                     }
                     else
                     {
-                        if(f2.cellIndex == GEO::NO_CELL)
+                        // Take the mirror facet to iterate over the next cell
+                        ci = mirrorFacet(outFacet).cellIndex;
+                        if(ci == GEO::NO_CELL)
                             ok = false;
-                        ci = f2.cellIndex;
                         p = lpi;
                     }
                 }
@@ -2052,19 +2067,24 @@ void DelaunayGraphCut::forceTedgesByGradientIJCV(bool fixesSigma, float nPixelSi
                     minSilent = std::min(minSilent, c.out);
                     maxSilent = std::max(maxSilent, c.out);
 
-                    Facet f1, f2;
+                    Facet outFacet;
                     Point3d lpi;
-                    // find cell which is farest to the cam and which intersect cam-p ray
+
+                    // Intersection with the next facet in the current tetrahedron (ci) in order to find the cell farest
+                    // to the cam which is intersected with cam-p ray
+                    // False here mean farest
+                    const bool nearestFarest = false;
                     if(((p - po).size() > nsigmaBackSilentPart * maxDist) || // (p-po).size() > 2 * sigma
-                       (!farestNeighbourCellToTheCamOnTheRay(mp->CArr[cam], p, ci, f1, f2, lpi)))
+                       !rayCellIntersection(mp->CArr[cam], p, ci, outFacet, nearestFarest, lpi))
                     {
                         ok = false;
                     }
                     else
                     {
-                        if(f2.cellIndex == GEO::NO_CELL)
+                        // Take the mirror facet to iterate over the next cell
+                        ci = mirrorFacet(outFacet).cellIndex;
+                        if(ci == GEO::NO_CELL)
                             ok = false;
-                        ci = f2.cellIndex;
                         p = lpi;
                     }
                 }
@@ -2197,26 +2217,31 @@ void DelaunayGraphCut::updateGraphFromTmpPtsCamsHexahRC(int rc, Point3d hexah[8]
                         _cellsAttr[tmp_ci].out += weight;
                     }
 
-                    Facet f1, f2;
+                    Facet outFacet;
                     Point3d lpi;
-                    // find cell which is nearest to the pont and which is
-                    // intersected with point-p ray
-                    if(!farestNeighbourCellToTheCamOnTheRay(camBehind, p, tmp_ci, f1, f2, lpi))
+
+                    // Intersection with the next facet in the current tetrahedron (ci) in order to find the cell farest
+                    // to the cam which is intersected with cam-p ray
+                    // False here mean farest
+                    const bool nearestFarest = false;
+                    if(!rayCellIntersection(camBehind, p, tmp_ci, outFacet, nearestFarest, lpi))
                     {
                         ok = false;
                     }
                     else
                     {
-                        if(f2.cellIndex == GEO::NO_CELL)
+                        // Take the mirror facet to iterate over the next cell
+                        const Facet mFacet = mirrorFacet(outFacet);
+                        tmp_ci = mFacet.cellIndex;
+                        if(tmp_ci == GEO::NO_CELL)
                         {
                             ok = false;
                         }
                         else
                         {
 #pragma OMP_ATOMIC_UPDATE
-                            _cellsAttr[f2.cellIndex].gEdgeVisWeight[f2.localVertexIndex] += weight;
+                            _cellsAttr[tmp_ci].gEdgeVisWeight[mFacet.localVertexIndex] += weight;
                         }
-                        tmp_ci = f2.cellIndex;
                         p = lpi;
                         ++nwup;
                     }
diff --git a/src/aliceVision/fuseCut/DelaunayGraphCut.hpp b/src/aliceVision/fuseCut/DelaunayGraphCut.hpp
index 17fa047ad9..4d8e2d6bb9 100644
--- a/src/aliceVision/fuseCut/DelaunayGraphCut.hpp
+++ b/src/aliceVision/fuseCut/DelaunayGraphCut.hpp
@@ -270,13 +270,6 @@ class DelaunayGraphCut
     bool rayCellIntersection(const Point3d& camCenter, const Point3d& p, GEO::index_t tetrahedronIndex, Facet& outFacet,
                              bool nearestFarest, Point3d& outIntersectPt) const;
 
-    inline bool neighbourCellToTheCamOnTheRay(const Point3d& camC, const Point3d& p, int tetrahedron, bool nearestFarest, 
-                                            Facet& f1, Facet& f2, Point3d& outIntersectPt) const;
-    inline bool nearestNeighbourCellToTheCamOnTheRay(const Point3d& camC, const Point3d& p, int tetrahedron, Facet& f1, Facet& f2,
-                                                     Point3d& outIntersectPt) const;
-    inline bool farestNeighbourCellToTheCamOnTheRay(const Point3d& camC, const Point3d& p, int tetrahedron, Facet& f1, Facet& f2,
-                                                    Point3d& outIntersectPt) const;
-
     inline Facet getFacetInFrontVertexOnTheRayToTheCam(int vertexIndex, int cam) const;
     Facet getFacetInFrontVertexOnTheRayToThePoint3d(VertexIndex vi, Point3d& ptt) const;
     Facet getFacetBehindVertexOnTheRayToTheCam(VertexIndex vi, int cam) const;
@@ -349,32 +342,6 @@ class DelaunayGraphCut
     mesh::Mesh* createMesh(bool filterHelperPointsTriangles = true);
 };
 
-inline bool DelaunayGraphCut::neighbourCellToTheCamOnTheRay(const Point3d& camC, const Point3d& p, int tetrahedron, bool nearestFarest, 
-                                                    Facet& f1, Facet& f2, Point3d& outIntersectPt) const
-{
-    f2.cellIndex = GEO::NO_CELL;
-    f2.localVertexIndex = GEO::NO_VERTEX;
-    outIntersectPt = p;
-
-    if(rayCellIntersection(camC, p, tetrahedron, f1, nearestFarest, outIntersectPt))
-    {
-        f2 = mirrorFacet(f1);
-        return true;
-    }
-    return false;
-}
-
-inline bool DelaunayGraphCut::nearestNeighbourCellToTheCamOnTheRay(const Point3d& camC, const Point3d& p, int tetrahedron, Facet& f1, Facet& f2, 
-                                                            Point3d& outIntersectPt) const
-{
-    return neighbourCellToTheCamOnTheRay(camC, p, tetrahedron, true, f1, f2, outIntersectPt);
-}
-inline bool DelaunayGraphCut::farestNeighbourCellToTheCamOnTheRay(const Point3d& camC, const Point3d& p, int tetrahedron, Facet& f1,  Facet& f2, 
-                                                            Point3d& outIntersectPt) const
-{
-    return neighbourCellToTheCamOnTheRay(camC, p, tetrahedron, false, f1, f2, outIntersectPt);
-}
-
 inline DelaunayGraphCut::Facet DelaunayGraphCut::getFacetInFrontVertexOnTheRayToTheCam(int vertexIndex,
                                                                                int cam) const
 {

From 596b3346a7ea74e288a015a059409f981bd8225a Mon Sep 17 00:00:00 2001
From: Enguerrand DE SMET <desmet.enguerrand@gmail.com>
Date: Tue, 7 Jul 2020 17:56:12 +0200
Subject: [PATCH 06/28] [fuseCut] DelaunayGraphCut: more meaningful variable
 names

po -> orignPt
lip -> intersectPt
---
 src/aliceVision/fuseCut/DelaunayGraphCut.cpp | 87 ++++++++++----------
 1 file changed, 43 insertions(+), 44 deletions(-)

diff --git a/src/aliceVision/fuseCut/DelaunayGraphCut.cpp b/src/aliceVision/fuseCut/DelaunayGraphCut.cpp
index 1c99864a4b..ce02844095 100644
--- a/src/aliceVision/fuseCut/DelaunayGraphCut.cpp
+++ b/src/aliceVision/fuseCut/DelaunayGraphCut.cpp
@@ -1387,14 +1387,14 @@ DelaunayGraphCut::Facet DelaunayGraphCut::getFacetInFrontVertexOnTheRayToThePoin
         if(isInfiniteCell(adjCellIndex))
             continue;
         Facet f1;
-        Point3d lpi;
-        if(rayCellIntersection(ptt, p, adjCellIndex, f1, true, lpi) == true)
+        Point3d intersectPt;
+        if(rayCellIntersection(ptt, p, adjCellIndex, f1, true, intersectPt) == true)
         {
             // if it is inbetween the camera and the point
-            if((ptt - lpi).size() < minDist)
+            if((ptt - intersectPt).size() < minDist)
             {
                 nearestFacet = f1;
-                minDist = (ptt - lpi).size();
+                minDist = (ptt - intersectPt).size();
             }
             // TODO FACA: maybe we can break and remove minDist?
         }
@@ -1683,8 +1683,8 @@ void DelaunayGraphCut::fillGraphPartPtRc(int& out_nstepsFront, int& out_nstepsBe
 
     int maxint = 1000000; // std::numeric_limits<int>::std::max()
 
-    const Point3d& po = _verticesCoords[vertexIndex];
-    const float pixSize = mp->getCamPixelSize(po, cam);
+    const Point3d& originPt = _verticesCoords[vertexIndex];
+    const float pixSize = mp->getCamPixelSize(originPt, cam);
     float maxDist = nPixelSizeBehind * pixSize;
     if(fixesSigma)
     {
@@ -1694,7 +1694,7 @@ void DelaunayGraphCut::fillGraphPartPtRc(int& out_nstepsFront, int& out_nstepsBe
     assert(cam >= 0);
     assert(cam < mp->ncams);
 
-    // printf("%f %f %f\n",po.x,po.y,po.z);
+    // printf("%f %f %f\n",originPt.x,originPt.y,originPt.z);
 
     int nsteps = 0;
 
@@ -1704,7 +1704,7 @@ void DelaunayGraphCut::fillGraphPartPtRc(int& out_nstepsFront, int& out_nstepsBe
         // tetrahedron connected to the point p and which intersect the ray from camera c to point p
         CellIndex ci = getFacetInFrontVertexOnTheRayToTheCam(vertexIndex, cam).cellIndex;
 
-        Point3d p = po;
+        Point3d p = originPt;
         CellIndex lastFinite = GEO::NO_CELL;
         bool ok = ci != GEO::NO_CELL;
         while(ok)
@@ -1719,20 +1719,19 @@ void DelaunayGraphCut::fillGraphPartPtRc(int& out_nstepsFront, int& out_nstepsBe
 
             Point3d pold = p;
             Facet outFacet;
-            Point3d lpi;
+            Point3d intersectPt;
 
             // Intersection with the next facet in the current tetrahedron (ci) in order to find the cell nearest to the
             // cam which is intersected with cam-p ray
             // true here mean nearest
             const bool nearestFarest = true;
-            if(!rayCellIntersection(mp->CArr[cam], p, ci, outFacet, nearestFarest, lpi))
+            if(!rayCellIntersection(mp->CArr[cam], p, ci, outFacet, nearestFarest, intersectPt))
             {
                 ok = false;
             }
             else
             {
-                float dist = distFcn(maxDist, (po - pold).size(), distFcnHeight);
-
+                float dist = distFcn(maxDist, (originPt - pold).size(), distFcnHeight);
                 {
 #pragma OMP_ATOMIC_UPDATE
                     _cellsAttr[outFacet.cellIndex].gEdgeVisWeight[outFacet.localVertexIndex] += weight * dist;
@@ -1743,7 +1742,7 @@ void DelaunayGraphCut::fillGraphPartPtRc(int& out_nstepsFront, int& out_nstepsBe
 
                 if(ci == GEO::NO_CELL)
                     ok = false;
-                p = lpi;
+                p = intersectPt;
                 lastFinite = ci;
             }
         }
@@ -1769,7 +1768,7 @@ void DelaunayGraphCut::fillGraphPartPtRc(int& out_nstepsFront, int& out_nstepsBe
             _cellsAttr[ci].on += weight;
         }
 
-        Point3d p = po; // HAS TO BE HERE !!!
+        Point3d p = originPt; // HAS TO BE HERE !!!
 
         bool ok = (ci != GEO::NO_CELL) && allPoints;
         while(ok)
@@ -1789,26 +1788,26 @@ void DelaunayGraphCut::fillGraphPartPtRc(int& out_nstepsFront, int& out_nstepsBe
             ++nsteps;
 
             Point3d pold = p;
-            Point3d lpi;
+            Point3d intersectPt;
 
             // Intersection with the next facet in the current tetrahedron (ci) in order to find the cell farest to the
             // cam which is intersected with cam-p ray
             // False here mean farest
             const bool nearestFarest = false;
-            if(!rayCellIntersection(mp->CArr[cam], p, ci, outFacet, nearestFarest, lpi) ||
-               ((po - pold).size() >= maxDist) || (!allPoints))
+            if(!rayCellIntersection(mp->CArr[cam], p, ci, outFacet, nearestFarest, intersectPt) ||
+               ((originPt - pold).size() >= maxDist) || (!allPoints))
             {
                 ok = false;
             }
             else
             {
                 // float dist = 1.0f;
-                // float dist = distFcn(maxDist,(po-pold).size()); // not sure if it is OK ... TODO
+                // float dist = distFcn(maxDist,(originPt-pold).size()); // not sure if it is OK ... TODO
                 // check ... with, without, and so on ....
                 // because labatutCFG09 with 32 gives much better result than nrc
                 // but when using just nrc and not using distFcn then the result is the same as labatutCGF09
 
-                float dist = distFcn(maxDist, (po - pold).size(), distFcnHeight);
+                float dist = distFcn(maxDist, (originPt - pold).size(), distFcnHeight);
 
                 // Take the mirror facet to iterate over the next cell
                 const Facet mFacet = mirrorFacet(outFacet);
@@ -1822,7 +1821,7 @@ void DelaunayGraphCut::fillGraphPartPtRc(int& out_nstepsFront, int& out_nstepsBe
 #pragma OMP_ATOMIC_UPDATE
                     _cellsAttr[ci].gEdgeVisWeight[mFacet.localVertexIndex] += weight * dist;
                 }
-                p = lpi;
+                p = intersectPt;
             }
         }
 
@@ -1865,7 +1864,7 @@ void DelaunayGraphCut::forceTedgesByGradientCVPR11(bool fixesSigma, float nPixel
         if(v.isVirtual())
             continue;
 
-        const Point3d& po = _verticesCoords[vi];
+        const Point3d& originPt = _verticesCoords[vi];
         for(int c = 0; c < v.cams.size(); ++c)
         {
             int cam = v.cams[c];
@@ -1880,7 +1879,7 @@ void DelaunayGraphCut::forceTedgesByGradientCVPR11(bool fixesSigma, float nPixel
                 float eFirst = _cellsAttr[fFirst.cellIndex].out;
 
                 CellIndex ci = f1.cellIndex;
-                Point3d p = po; // HAS TO BE HERE !!!
+                Point3d p = originPt; // HAS TO BE HERE !!!
                 float maxDist = nPixelSizeBehind * mp->getCamPixelSize(p, cam);
                 if(fixesSigma)
                 {
@@ -1891,15 +1890,15 @@ void DelaunayGraphCut::forceTedgesByGradientCVPR11(bool fixesSigma, float nPixel
                 while(ok)
                 {
                     Point3d pold = p;
-                    Point3d lpi;
+                    Point3d intersectPt;
                     Facet outFacet;
 
                     // Intersection with the next facet in the current tetrahedron (ci) in order to find the cell farest
                     // to the cam which is intersected with cam-p ray
                     // False here mean farest
                     const bool nearestFarest = false;
-                    if(!rayCellIntersection(mp->CArr[cam], p, ci, outFacet, nearestFarest, lpi) ||
-                       ((po - pold).size() >= maxDist))
+                    if(!rayCellIntersection(mp->CArr[cam], p, ci, outFacet, nearestFarest, intersectPt) ||
+                       ((originPt - pold).size() >= maxDist))
                     {
                         ok = false;
                     }
@@ -1909,7 +1908,7 @@ void DelaunayGraphCut::forceTedgesByGradientCVPR11(bool fixesSigma, float nPixel
                         ci = mirrorFacet(outFacet).cellIndex;
                         if(ci == GEO::NO_CELL)
                             ok = false;
-                        p = lpi;
+                        p = intersectPt;
                     }
                 }
 
@@ -1985,7 +1984,7 @@ void DelaunayGraphCut::forceTedgesByGradientIJCV(bool fixesSigma, float nPixelSi
         if(v.isVirtual())
             continue;
 
-        const Point3d& po = _verticesCoords[vi];
+        const Point3d& originPt = _verticesCoords[vi];
         for(int c = 0; c < v.cams.size(); ++c)
         {
             int nstepsFront = 0;
@@ -1999,7 +1998,7 @@ void DelaunayGraphCut::forceTedgesByGradientIJCV(bool fixesSigma, float nPixelSi
             }
             else
             {
-                maxDist = nPixelSizeBehind * mp->getCamPixelSize(po, cam);
+                maxDist = nPixelSizeBehind * mp->getCamPixelSize(originPt, cam);
             }
 
             float minJump = 10000000.0f;
@@ -2010,14 +2009,14 @@ void DelaunayGraphCut::forceTedgesByGradientIJCV(bool fixesSigma, float nPixelSi
 
             {
                 CellIndex ci = getFacetInFrontVertexOnTheRayToTheCam(vi, cam).cellIndex;
-                Point3d p = po; // HAS TO BE HERE !!!
+                Point3d p = originPt; // HAS TO BE HERE !!!
                 bool ok = (ci != GEO::NO_CELL);
                 while(ok)
                 {
                     ++nstepsFront;
 
                     const GC_cellInfo& c = _cellsAttr[ci];
-                    if((p - po).size() > nsigmaFrontSilentPart * maxDist) // (p-po).size() > 2 * sigma
+                    if((p - originPt).size() > nsigmaFrontSilentPart * maxDist) // (p-originPt).size() > 2 * sigma
                     {
                         minJump = std::min(minJump, c.out);
                         maxJump = std::max(maxJump, c.out);
@@ -2029,13 +2028,13 @@ void DelaunayGraphCut::forceTedgesByGradientIJCV(bool fixesSigma, float nPixelSi
                     }
 
                     Facet outFacet;
-                    Point3d lpi;
+                    Point3d intersectPt;
                     // Intersection with the next facet in the current tetrahedron (ci) in order to find the cell nearest
                     // to the cam which is intersected with cam-p ray
                     // True here mean nearest
                     const bool nearestFarest = true;
-                    if(((p - po).size() > (nsigmaJumpPart + nsigmaFrontSilentPart) * maxDist) || // (2 + 2) * sigma
-                       !rayCellIntersection(mp->CArr[cam], p, ci, outFacet, nearestFarest, lpi))
+                    if(((p - originPt).size() > (nsigmaJumpPart + nsigmaFrontSilentPart) * maxDist) || // (2 + 2) * sigma
+                       !rayCellIntersection(mp->CArr[cam], p, ci, outFacet, nearestFarest, intersectPt))
                     {
                         ok = false;
                     }
@@ -2045,14 +2044,14 @@ void DelaunayGraphCut::forceTedgesByGradientIJCV(bool fixesSigma, float nPixelSi
                         ci = mirrorFacet(outFacet).cellIndex;
                         if(ci == GEO::NO_CELL)
                             ok = false;
-                        p = lpi;
+                        p = intersectPt;
                     }
                 }
             }
 
             {
                 CellIndex ci = getFacetBehindVertexOnTheRayToTheCam(vi, cam).cellIndex; // T1
-                Point3d p = po; // HAS TO BE HERE !!!
+                Point3d p = originPt; // HAS TO BE HERE !!!
                 bool ok = (ci != GEO::NO_CELL);
                 if(ok)
                 {
@@ -2068,14 +2067,14 @@ void DelaunayGraphCut::forceTedgesByGradientIJCV(bool fixesSigma, float nPixelSi
                     maxSilent = std::max(maxSilent, c.out);
 
                     Facet outFacet;
-                    Point3d lpi;
+                    Point3d intersectPt;
 
                     // Intersection with the next facet in the current tetrahedron (ci) in order to find the cell farest
                     // to the cam which is intersected with cam-p ray
                     // False here mean farest
                     const bool nearestFarest = false;
-                    if(((p - po).size() > nsigmaBackSilentPart * maxDist) || // (p-po).size() > 2 * sigma
-                       !rayCellIntersection(mp->CArr[cam], p, ci, outFacet, nearestFarest, lpi))
+                    if(((p - originPt).size() > nsigmaBackSilentPart * maxDist) || // (p-originPt).size() > 2 * sigma
+                       !rayCellIntersection(mp->CArr[cam], p, ci, outFacet, nearestFarest, intersectPt))
                     {
                         ok = false;
                     }
@@ -2085,7 +2084,7 @@ void DelaunayGraphCut::forceTedgesByGradientIJCV(bool fixesSigma, float nPixelSi
                         ci = mirrorFacet(outFacet).cellIndex;
                         if(ci == GEO::NO_CELL)
                             ok = false;
-                        p = lpi;
+                        p = intersectPt;
                     }
                 }
 
@@ -2198,9 +2197,9 @@ void DelaunayGraphCut::updateGraphFromTmpPtsCamsHexahRC(int rc, Point3d hexah[8]
            ((!doFilterOctreeTracks) || ((doFilterOctreeTracks) && (pnt.ncams >= minNumOfConsistentCams))))
         {
             nin++;
-            Point3d lpi = pt;
+            Point3d intersectPt = pt;
             Point3d camBehind = mp->CArr[rc] + (mp->CArr[rc] - pt);
-            CellIndex ci = getFirstCellOnTheRayFromCamToThePoint(rc, pt, lpi);
+            CellIndex ci = getFirstCellOnTheRayFromCamToThePoint(rc, pt, intersectPt);
             if(ci != GEO::NO_CELL)
             {
                 // update weights on the sp-cam half line
@@ -2218,13 +2217,13 @@ void DelaunayGraphCut::updateGraphFromTmpPtsCamsHexahRC(int rc, Point3d hexah[8]
                     }
 
                     Facet outFacet;
-                    Point3d lpi;
+                    Point3d intersectPt;
 
                     // Intersection with the next facet in the current tetrahedron (ci) in order to find the cell farest
                     // to the cam which is intersected with cam-p ray
                     // False here mean farest
                     const bool nearestFarest = false;
-                    if(!rayCellIntersection(camBehind, p, tmp_ci, outFacet, nearestFarest, lpi))
+                    if(!rayCellIntersection(camBehind, p, tmp_ci, outFacet, nearestFarest, intersectPt))
                     {
                         ok = false;
                     }
@@ -2242,7 +2241,7 @@ void DelaunayGraphCut::updateGraphFromTmpPtsCamsHexahRC(int rc, Point3d hexah[8]
 #pragma OMP_ATOMIC_UPDATE
                             _cellsAttr[tmp_ci].gEdgeVisWeight[mFacet.localVertexIndex] += weight;
                         }
-                        p = lpi;
+                        p = intersectPt;
                         ++nwup;
                     }
                 }

From 0933d84063a65a0acf2209549964a4757cec2017 Mon Sep 17 00:00:00 2001
From: Enguerrand DE SMET <desmet.enguerrand@gmail.com>
Date: Tue, 7 Jul 2020 18:44:59 +0200
Subject: [PATCH 07/28] [fuseCut] DelaunayGraphCut: remove and rename
 unnecessary function

Remove and rename unnecessary getFacetInFrontVertexOnTheRayToThePoint3d function using directly getFacetInFrontVertexOnTheRayToTheCam function.
---
 src/aliceVision/fuseCut/DelaunayGraphCut.cpp | 26 +++++++++++---------
 src/aliceVision/fuseCut/DelaunayGraphCut.hpp | 20 +--------------
 2 files changed, 16 insertions(+), 30 deletions(-)

diff --git a/src/aliceVision/fuseCut/DelaunayGraphCut.cpp b/src/aliceVision/fuseCut/DelaunayGraphCut.cpp
index ce02844095..9229267d40 100644
--- a/src/aliceVision/fuseCut/DelaunayGraphCut.cpp
+++ b/src/aliceVision/fuseCut/DelaunayGraphCut.cpp
@@ -1368,33 +1368,37 @@ bool DelaunayGraphCut::rayCellIntersection(const Point3d& camCenter, const Point
     return true;
 }
 
-DelaunayGraphCut::Facet DelaunayGraphCut::getFacetInFrontVertexOnTheRayToThePoint3d(VertexIndex vi,
-                                                                                   Point3d& ptt) const
+DelaunayGraphCut::Facet DelaunayGraphCut::getFacetInFrontVertexOnTheRayToTheCam(VertexIndex vertexIndex, int cam) const
 {
-    const Point3d& p = _verticesCoords[vi];
+    if((cam < 0) || (cam >= mp->ncams))
+    {
+        ALICEVISION_LOG_WARNING("Bad camId, cam: " << cam << ", ptid: " << vertexIndex);
+    }
+    const Point3d& p = _verticesCoords[vertexIndex];
 
-    double minDist = (ptt - p).size(); // initialize minDist to the distance from 3d point p to camera center c
+    double minDist = (mp->CArr[cam] - p).size(); // initialize minDist to the distance from 3d point p to camera center
     Facet nearestFacet;
     nearestFacet.cellIndex = GEO::NO_CELL;
     nearestFacet.localVertexIndex = GEO::NO_VERTEX;
 
     for(int k = 0; true; ++k)
     {
-        CellIndex adjCellIndex = vertexToCells(vi, k); // GEOGRAM: set_stores_cicl(true) required
+        CellIndex adjCellIndex = vertexToCells(vertexIndex, k); // GEOGRAM: set_stores_cicl(true) required
         if(adjCellIndex == GEO::NO_CELL) // last one
             break;
-
         if(isInfiniteCell(adjCellIndex))
             continue;
-        Facet f1;
+
+        Facet outFacet;
         Point3d intersectPt;
-        if(rayCellIntersection(ptt, p, adjCellIndex, f1, true, intersectPt) == true)
+        if(rayCellIntersection(mp->CArr[cam], p, adjCellIndex, outFacet, true, intersectPt) == true)
         {
+            const double intersectDist = (mp->CArr[cam] - intersectPt).size();
             // if it is inbetween the camera and the point
-            if((ptt - intersectPt).size() < minDist)
+            if(intersectDist < minDist)
             {
-                nearestFacet = f1;
-                minDist = (ptt - intersectPt).size();
+                nearestFacet = outFacet;
+                minDist = intersectDist;
             }
             // TODO FACA: maybe we can break and remove minDist?
         }
diff --git a/src/aliceVision/fuseCut/DelaunayGraphCut.hpp b/src/aliceVision/fuseCut/DelaunayGraphCut.hpp
index 4d8e2d6bb9..ac72f2a026 100644
--- a/src/aliceVision/fuseCut/DelaunayGraphCut.hpp
+++ b/src/aliceVision/fuseCut/DelaunayGraphCut.hpp
@@ -270,8 +270,7 @@ class DelaunayGraphCut
     bool rayCellIntersection(const Point3d& camCenter, const Point3d& p, GEO::index_t tetrahedronIndex, Facet& outFacet,
                              bool nearestFarest, Point3d& outIntersectPt) const;
 
-    inline Facet getFacetInFrontVertexOnTheRayToTheCam(int vertexIndex, int cam) const;
-    Facet getFacetInFrontVertexOnTheRayToThePoint3d(VertexIndex vi, Point3d& ptt) const;
+    Facet getFacetInFrontVertexOnTheRayToTheCam(VertexIndex vertexIndex, int cam) const;
     Facet getFacetBehindVertexOnTheRayToTheCam(VertexIndex vi, int cam) const;
     GEO::index_t getFirstCellOnTheRayFromCamToThePoint(int cam, Point3d& p, Point3d& lpi) const;
 
@@ -342,22 +341,5 @@ class DelaunayGraphCut
     mesh::Mesh* createMesh(bool filterHelperPointsTriangles = true);
 };
 
-inline DelaunayGraphCut::Facet DelaunayGraphCut::getFacetInFrontVertexOnTheRayToTheCam(int vertexIndex,
-                                                                               int cam) const
-{
-    // if (btest) {
-    //	printf("cam %i pt %f %f %f\n",cam,p.x,p.y,p.z);
-    //	printf("ptid %i\n",vp->info().id);
-    //	printf("campos %f %f
-    //%f\n",mp->CArr[cam].x,mp->CArr[cam].y,mp->CArr[cam].z);
-    //};
-    if((cam < 0) || (cam >= mp->ncams))
-    {
-        ALICEVISION_LOG_WARNING("Bad camId, cam: " << cam << ", ptid: " << vertexIndex);
-    }
-
-    return getFacetInFrontVertexOnTheRayToThePoint3d(vertexIndex, mp->CArr[cam]);
-}
-
 } // namespace fuseCut
 } // namespace aliceVision

From 0755c86a8a8b98420898bb90bc13eda49c1bc28f Mon Sep 17 00:00:00 2001
From: Enguerrand DE SMET <desmet.enguerrand@gmail.com>
Date: Tue, 7 Jul 2020 19:34:36 +0200
Subject: [PATCH 08/28] [fuseCut] DelaunayGraphCut:  rename and group duplicate
 codes

about getFacetInFrontVertexOnTheRayToTheCam and getFacetBehindVertexOnTheRayToTheCam which use the same code except one boolean.
---
 src/aliceVision/fuseCut/DelaunayGraphCut.cpp | 76 ++++++++------------
 src/aliceVision/fuseCut/DelaunayGraphCut.hpp |  3 +-
 2 files changed, 29 insertions(+), 50 deletions(-)

diff --git a/src/aliceVision/fuseCut/DelaunayGraphCut.cpp b/src/aliceVision/fuseCut/DelaunayGraphCut.cpp
index 9229267d40..4940bfa651 100644
--- a/src/aliceVision/fuseCut/DelaunayGraphCut.cpp
+++ b/src/aliceVision/fuseCut/DelaunayGraphCut.cpp
@@ -1368,7 +1368,7 @@ bool DelaunayGraphCut::rayCellIntersection(const Point3d& camCenter, const Point
     return true;
 }
 
-DelaunayGraphCut::Facet DelaunayGraphCut::getFacetInFrontVertexOnTheRayToTheCam(VertexIndex vertexIndex, int cam) const
+DelaunayGraphCut::Facet DelaunayGraphCut::getFacetFromVertexOnTheRayToTheCam(VertexIndex vertexIndex, int cam, bool nearestFarest) const
 {
     if((cam < 0) || (cam >= mp->ncams))
     {
@@ -1376,10 +1376,8 @@ DelaunayGraphCut::Facet DelaunayGraphCut::getFacetInFrontVertexOnTheRayToTheCam(
     }
     const Point3d& p = _verticesCoords[vertexIndex];
 
-    double minDist = (mp->CArr[cam] - p).size(); // initialize minDist to the distance from 3d point p to camera center
-    Facet nearestFacet;
-    nearestFacet.cellIndex = GEO::NO_CELL;
-    nearestFacet.localVertexIndex = GEO::NO_VERTEX;
+    double currentDist = (mp->CArr[cam] - p).size(); // initialize currentDist to the distance from 3d point p to camera center
+    Facet facet;
 
     for(int k = 0; true; ++k)
     {
@@ -1390,51 +1388,22 @@ DelaunayGraphCut::Facet DelaunayGraphCut::getFacetInFrontVertexOnTheRayToTheCam(
             continue;
 
         Facet outFacet;
+
         Point3d intersectPt;
-        if(rayCellIntersection(mp->CArr[cam], p, adjCellIndex, outFacet, true, intersectPt) == true)
+        if(rayCellIntersection(mp->CArr[cam], p, adjCellIndex, outFacet, nearestFarest, intersectPt) == true)
         {
             const double intersectDist = (mp->CArr[cam] - intersectPt).size();
-            // if it is inbetween the camera and the point
-            if(intersectDist < minDist)
+            // if it's between the camera and the point if nearestFarest == true
+            // if it's behind the point (from the camera) if nearestFarest == false
+            if(nearestFarest ? (intersectDist < currentDist) : (intersectDist > currentDist))
             {
-                nearestFacet = outFacet;
-                minDist = intersectDist;
+                facet = outFacet;
+                currentDist = intersectDist;
             }
             // TODO FACA: maybe we can break and remove minDist?
         }
     }
-    return nearestFacet;
-}
-
-DelaunayGraphCut::Facet DelaunayGraphCut::getFacetBehindVertexOnTheRayToTheCam(VertexIndex vi,
-                                                                              int cam) const
-{
-    const Point3d& p = _verticesCoords[vi];
-
-    double maxDist = (mp->CArr[cam] - p).size();
-    Facet farestFacet;
-
-    for(int k = 0; true; ++k)
-    {
-        CellIndex adjCellIndex = vertexToCells(vi, k); // GEOGRAM: set_stores_cicl(true) required
-        if(adjCellIndex == GEO::NO_CELL) // last one
-            break;
-        if(isInfiniteCell(adjCellIndex))
-            continue;
-
-        Facet outFacet;
-        Point3d intersectionPoint;
-        if(rayCellIntersection(mp->CArr[cam], p, adjCellIndex, outFacet, false, intersectionPoint))
-        {
-            // if it is after the point p (along the axis from the camera)
-            if((mp->CArr[cam] - intersectionPoint).size() > maxDist)
-            {
-                farestFacet = outFacet;
-                maxDist = (mp->CArr[cam] - intersectionPoint).size();
-            }
-        }
-    }
-    return farestFacet;
+    return facet;
 }
 
 GEO::index_t DelaunayGraphCut::getFirstCellOnTheRayFromCamToThePoint(int cam, Point3d& p, Point3d& intersectPoint) const
@@ -1705,8 +1674,10 @@ void DelaunayGraphCut::fillGraphPartPtRc(int& out_nstepsFront, int& out_nstepsBe
     if(fillOut)
     {
         out_nstepsFront = 0;
+        // true here mean nearest
+        const bool nearestFarest = true;
         // tetrahedron connected to the point p and which intersect the ray from camera c to point p
-        CellIndex ci = getFacetInFrontVertexOnTheRayToTheCam(vertexIndex, cam).cellIndex;
+        CellIndex ci = getFacetFromVertexOnTheRayToTheCam(vertexIndex, cam, nearestFarest).cellIndex;
 
         Point3d p = originPt;
         CellIndex lastFinite = GEO::NO_CELL;
@@ -1762,7 +1733,9 @@ void DelaunayGraphCut::fillGraphPartPtRc(int& out_nstepsFront, int& out_nstepsBe
     {
         out_nstepsBehind = 0;
         // get the tetrahedron next to point p on the ray from c
-        Facet f1 = getFacetBehindVertexOnTheRayToTheCam(vertexIndex, cam);
+        // False here mean farest
+        const bool nearestFarest = false;
+        Facet f1 = getFacetFromVertexOnTheRayToTheCam(vertexIndex, cam, nearestFarest);
         Facet outFacet;
 
         CellIndex ci = f1.cellIndex;
@@ -1873,10 +1846,13 @@ void DelaunayGraphCut::forceTedgesByGradientCVPR11(bool fixesSigma, float nPixel
         {
             int cam = v.cams[c];
 
-            Facet fFirst = getFacetInFrontVertexOnTheRayToTheCam(vi, cam);
+            // False here mean nearest
+            const bool nearestFarest = true;
+            Facet fFirst = getFacetFromVertexOnTheRayToTheCam(vi, cam, nearestFarest);
 
             // get the tetrahedron next to point p on the ray from c
-            Facet f1 = getFacetBehindVertexOnTheRayToTheCam(vi, cam);
+            // False here mean farest
+            Facet f1 = getFacetFromVertexOnTheRayToTheCam(vi, cam, false);
 
             if((fFirst.cellIndex != GEO::NO_CELL) && (f1.cellIndex != GEO::NO_CELL) && (!isInfiniteCell(f1.cellIndex)))
             {
@@ -2012,7 +1988,9 @@ void DelaunayGraphCut::forceTedgesByGradientIJCV(bool fixesSigma, float nPixelSi
             float midSilent = 10000000.0f;
 
             {
-                CellIndex ci = getFacetInFrontVertexOnTheRayToTheCam(vi, cam).cellIndex;
+                // True here mean nearest
+                const bool nearestFarest = true;
+                CellIndex ci = getFacetFromVertexOnTheRayToTheCam(vi, cam, nearestFarest).cellIndex;
                 Point3d p = originPt; // HAS TO BE HERE !!!
                 bool ok = (ci != GEO::NO_CELL);
                 while(ok)
@@ -2054,7 +2032,9 @@ void DelaunayGraphCut::forceTedgesByGradientIJCV(bool fixesSigma, float nPixelSi
             }
 
             {
-                CellIndex ci = getFacetBehindVertexOnTheRayToTheCam(vi, cam).cellIndex; // T1
+                // False here mean farest
+                const bool nearestFarest = false;
+                CellIndex ci = getFacetFromVertexOnTheRayToTheCam(vi, cam, nearestFarest).cellIndex; // T1
                 Point3d p = originPt; // HAS TO BE HERE !!!
                 bool ok = (ci != GEO::NO_CELL);
                 if(ok)
diff --git a/src/aliceVision/fuseCut/DelaunayGraphCut.hpp b/src/aliceVision/fuseCut/DelaunayGraphCut.hpp
index ac72f2a026..2f92642c0e 100644
--- a/src/aliceVision/fuseCut/DelaunayGraphCut.hpp
+++ b/src/aliceVision/fuseCut/DelaunayGraphCut.hpp
@@ -270,8 +270,7 @@ class DelaunayGraphCut
     bool rayCellIntersection(const Point3d& camCenter, const Point3d& p, GEO::index_t tetrahedronIndex, Facet& outFacet,
                              bool nearestFarest, Point3d& outIntersectPt) const;
 
-    Facet getFacetInFrontVertexOnTheRayToTheCam(VertexIndex vertexIndex, int cam) const;
-    Facet getFacetBehindVertexOnTheRayToTheCam(VertexIndex vi, int cam) const;
+    Facet getFacetFromVertexOnTheRayToTheCam(VertexIndex vertexIndex, int cam, bool nearestFarest) const;
     GEO::index_t getFirstCellOnTheRayFromCamToThePoint(int cam, Point3d& p, Point3d& lpi) const;
 
     float distFcn(float maxDist, float dist, float distFcnHeight) const;

From cbdbf0aaeae8e12b439605950596e75238b6586f Mon Sep 17 00:00:00 2001
From: Enguerrand DE SMET <desmet.enguerrand@gmail.com>
Date: Wed, 8 Jul 2020 08:35:36 +0200
Subject: [PATCH 09/28] [fusecut] DelaunayGraphCut: clean and fix wrong
 functions signatures

---
 src/aliceVision/fuseCut/DelaunayGraphCut.cpp |  2 +-
 src/aliceVision/fuseCut/DelaunayGraphCut.hpp | 14 ++------------
 2 files changed, 3 insertions(+), 13 deletions(-)

diff --git a/src/aliceVision/fuseCut/DelaunayGraphCut.cpp b/src/aliceVision/fuseCut/DelaunayGraphCut.cpp
index 4940bfa651..cb93e0eaec 100644
--- a/src/aliceVision/fuseCut/DelaunayGraphCut.cpp
+++ b/src/aliceVision/fuseCut/DelaunayGraphCut.cpp
@@ -2654,7 +2654,7 @@ void DelaunayGraphCut::maxflow()
 }
 
 void DelaunayGraphCut::reconstructExpetiments(const StaticVector<int>& cams, const std::string& folderName,
-                                            bool update, Point3d* hexahInflated, const std::string& tmpCamsPtsFolderName,
+                                            bool update, Point3d hexahInflated[8], const std::string& tmpCamsPtsFolderName,
                                             const Point3d& spaceSteps)
 {
     ALICEVISION_LOG_INFO("DelaunayGraphCut::reconstructExpetiments");
diff --git a/src/aliceVision/fuseCut/DelaunayGraphCut.hpp b/src/aliceVision/fuseCut/DelaunayGraphCut.hpp
index 2f92642c0e..e265e0e6db 100644
--- a/src/aliceVision/fuseCut/DelaunayGraphCut.hpp
+++ b/src/aliceVision/fuseCut/DelaunayGraphCut.hpp
@@ -251,19 +251,16 @@ class DelaunayGraphCut
 
     void addPointsFromSfM(const Point3d hexah[8], const StaticVector<int>& cams, const sfmData::SfMData& sfmData);
     void addPointsFromCameraCenters(const StaticVector<int>& cams, float minDist);
-    void addPointsToPreventSingularities(const Point3d Voxel[], float minDist);
+    void addPointsToPreventSingularities(const Point3d Voxel[8], float minDist);
 
     /**
      * @brief Add volume points to prevent singularities
      */
-    void addHelperPoints(int nGridHelperVolumePointsDim, const Point3d Voxel[], float minDist);
+    void addHelperPoints(int nGridHelperVolumePointsDim, const Point3d Voxel[8], float minDist);
 
     void fuseFromDepthMaps(const StaticVector<int>& cams, const Point3d voxel[8], const FuseParams& params);
     void loadPrecomputedDensePoints(const StaticVector<int>* voxelsIds, const Point3d voxel[8], VoxelsGrid* ls);
 
-    void createTetrahedralizationFromDepthMapsCamsVoxel(const StaticVector<int>& allCams,
-                                                   StaticVector<int>* voxelsIds, Point3d Voxel[8], VoxelsGrid* ls);
-
     void computeVerticesSegSize(bool allPoints, float alpha = 0.0f);
     void removeSmallSegs(int minSegSize);
 
@@ -280,7 +277,6 @@ class DelaunayGraphCut
     double maxEdgeLength() const;
     Point3d cellCircumScribedSphereCentre(CellIndex ci) const;
     double getFaceWeight(const Facet &f1) const;
-    float weightFromSim(float sim);
 
     float weightFcn(float nrc, bool labatutWeights, int ncams);
 
@@ -326,12 +322,6 @@ class DelaunayGraphCut
 
     void graphCutPostProcessing();
 
-    void clearAllPointsInFreeSpace();
-    void clearAllPointsNotOnSurface();
-    void addNewPointsToOccupiedSpace();
-    void clearOutAddIn();
-    StaticVector<int>* getNearestTrisFromMeshTris(mesh::Mesh* otherMesh);
-
     void segmentFullOrFree(bool full, StaticVector<int>** inColors, int& nsegments);
     int removeBubbles();
     int removeDust(int minSegSize);

From 9873b13efb89fdaf6ddfe5ddd84ab1f2cd9103d0 Mon Sep 17 00:00:00 2001
From: Enguerrand DE SMET <desmet.enguerrand@gmail.com>
Date: Wed, 8 Jul 2020 08:52:30 +0200
Subject: [PATCH 10/28] [fuseCut] DelaunayGraphCut:  unnecessary  variables

Delete the unnecessary pold variable and move the dist variable.
---
 src/aliceVision/fuseCut/DelaunayGraphCut.cpp | 8 +++-----
 1 file changed, 3 insertions(+), 5 deletions(-)

diff --git a/src/aliceVision/fuseCut/DelaunayGraphCut.cpp b/src/aliceVision/fuseCut/DelaunayGraphCut.cpp
index cb93e0eaec..0c1fe20dd1 100644
--- a/src/aliceVision/fuseCut/DelaunayGraphCut.cpp
+++ b/src/aliceVision/fuseCut/DelaunayGraphCut.cpp
@@ -1692,7 +1692,6 @@ void DelaunayGraphCut::fillGraphPartPtRc(int& out_nstepsFront, int& out_nstepsBe
             ++out_nstepsFront;
             ++nsteps;
 
-            Point3d pold = p;
             Facet outFacet;
             Point3d intersectPt;
 
@@ -1706,8 +1705,8 @@ void DelaunayGraphCut::fillGraphPartPtRc(int& out_nstepsFront, int& out_nstepsBe
             }
             else
             {
-                float dist = distFcn(maxDist, (originPt - pold).size(), distFcnHeight);
                 {
+                    const float dist = distFcn(maxDist, (originPt - p).size(), distFcnHeight);
 #pragma OMP_ATOMIC_UPDATE
                     _cellsAttr[outFacet.cellIndex].gEdgeVisWeight[outFacet.localVertexIndex] += weight * dist;
                 }
@@ -1764,7 +1763,6 @@ void DelaunayGraphCut::fillGraphPartPtRc(int& out_nstepsFront, int& out_nstepsBe
             ++out_nstepsBehind;
             ++nsteps;
 
-            Point3d pold = p;
             Point3d intersectPt;
 
             // Intersection with the next facet in the current tetrahedron (ci) in order to find the cell farest to the
@@ -1772,7 +1770,7 @@ void DelaunayGraphCut::fillGraphPartPtRc(int& out_nstepsFront, int& out_nstepsBe
             // False here mean farest
             const bool nearestFarest = false;
             if(!rayCellIntersection(mp->CArr[cam], p, ci, outFacet, nearestFarest, intersectPt) ||
-               ((originPt - pold).size() >= maxDist) || (!allPoints))
+               ((originPt - p).size() >= maxDist) || (!allPoints))
             {
                 ok = false;
             }
@@ -1784,7 +1782,6 @@ void DelaunayGraphCut::fillGraphPartPtRc(int& out_nstepsFront, int& out_nstepsBe
                 // because labatutCFG09 with 32 gives much better result than nrc
                 // but when using just nrc and not using distFcn then the result is the same as labatutCGF09
 
-                float dist = distFcn(maxDist, (originPt - pold).size(), distFcnHeight);
 
                 // Take the mirror facet to iterate over the next cell
                 const Facet mFacet = mirrorFacet(outFacet);
@@ -1795,6 +1792,7 @@ void DelaunayGraphCut::fillGraphPartPtRc(int& out_nstepsFront, int& out_nstepsBe
                 }
                 else
                 {
+                    const float dist = distFcn(maxDist, (originPt - p).size(), distFcnHeight);
 #pragma OMP_ATOMIC_UPDATE
                     _cellsAttr[ci].gEdgeVisWeight[mFacet.localVertexIndex] += weight * dist;
                 }

From 4e196776989cca40a65e8fe4ef79257c05900b6f Mon Sep 17 00:00:00 2001
From: Enguerrand DE SMET <desmet.enguerrand@gmail.com>
Date: Wed, 8 Jul 2020 10:02:54 +0200
Subject: [PATCH 11/28] [fuseCut] DelaunayGraphCut: more meaningful loop

Add a new function getNeighboringCellsByVertexIndex used to loop over neighboring cells in a more meaningful way than using a loop 'for(int k = 0; true; ++k)'.
---
 src/aliceVision/fuseCut/DelaunayGraphCut.cpp | 15 +++------------
 src/aliceVision/fuseCut/DelaunayGraphCut.hpp | 11 +++++++++++
 2 files changed, 14 insertions(+), 12 deletions(-)

diff --git a/src/aliceVision/fuseCut/DelaunayGraphCut.cpp b/src/aliceVision/fuseCut/DelaunayGraphCut.cpp
index 0c1fe20dd1..4e771187be 100644
--- a/src/aliceVision/fuseCut/DelaunayGraphCut.cpp
+++ b/src/aliceVision/fuseCut/DelaunayGraphCut.cpp
@@ -1379,11 +1379,8 @@ DelaunayGraphCut::Facet DelaunayGraphCut::getFacetFromVertexOnTheRayToTheCam(Ver
     double currentDist = (mp->CArr[cam] - p).size(); // initialize currentDist to the distance from 3d point p to camera center
     Facet facet;
 
-    for(int k = 0; true; ++k)
+    for(CellIndex adjCellIndex : getNeighboringCellsByVertexIndex(vertexIndex)) // GEOGRAM: set_stores_cicl(true) required
     {
-        CellIndex adjCellIndex = vertexToCells(vertexIndex, k); // GEOGRAM: set_stores_cicl(true) required
-        if(adjCellIndex == GEO::NO_CELL) // last one
-            break;
         if(isInfiniteCell(adjCellIndex))
             continue;
 
@@ -1414,11 +1411,8 @@ GEO::index_t DelaunayGraphCut::getFirstCellOnTheRayFromCamToThePoint(int cam, Po
     float maxdist = 0.0f;
     GEO::index_t farestCell = GEO::NO_CELL;
 
-    for(int k = 0; true; ++k)
+    for(CellIndex adjCellIndex : getNeighboringCellsByVertexIndex(cam_vi)) // GEOGRAM: set_stores_cicl(true) required
     {
-        CellIndex adjCellIndex = vertexToCells(cam_vi, k); // GEOGRAM: set_stores_cicl(true) required
-        if(adjCellIndex == GEO::NO_CELL) // last one
-            break;
         if(isInfiniteCell(adjCellIndex))
             continue;
 
@@ -2335,11 +2329,8 @@ void DelaunayGraphCut::freeUnwantedFullCells(const Point3d* hexah)
         VertexIndex cam_vi = _camsVertexes[rc];
         if(cam_vi == GEO::NO_VERTEX)
             continue;
-        for(int k = 0; true; ++k)
+        for(CellIndex adjCellIndex : getNeighboringCellsByVertexIndex(cam_vi)) // GEOGRAM: set_stores_cicl(true) required
         {
-            CellIndex adjCellIndex = vertexToCells(cam_vi, k); // GEOGRAM: set_stores_cicl(true) required
-            if(adjCellIndex == GEO::NO_CELL) // last one
-                break;
             if(isInfiniteCell(adjCellIndex))
                 continue;
 
diff --git a/src/aliceVision/fuseCut/DelaunayGraphCut.hpp b/src/aliceVision/fuseCut/DelaunayGraphCut.hpp
index e265e0e6db..2b6e693579 100644
--- a/src/aliceVision/fuseCut/DelaunayGraphCut.hpp
+++ b/src/aliceVision/fuseCut/DelaunayGraphCut.hpp
@@ -234,6 +234,17 @@ class DelaunayGraphCut
         return localCells[lvi];
     }
 
+    /**
+     * @brief Retrieves the global indexes of neighboring cells using the global index of a vertex.
+     * 
+     * @param vi the global vertexIndex
+     * @return a vector of neighboring cell indexes
+     */
+    inline const std::vector<CellIndex>& getNeighboringCellsByVertexIndex(VertexIndex vi) const
+    {
+        return _neighboringCellsPerVertex.at(vi);
+    }
+
     void initVertices();
     void computeDelaunay();
     void initCells();

From df1fd8907c772292a66ca7cbccb660f6b9dd080a Mon Sep 17 00:00:00 2001
From: Enguerrand DE SMET <desmet.enguerrand@gmail.com>
Date: Wed, 8 Jul 2020 10:35:02 +0200
Subject: [PATCH 12/28] [fuseCut] DelaunayGraphCut: add function documentation

---
 src/aliceVision/fuseCut/DelaunayGraphCut.hpp | 17 +++++++++++++++--
 1 file changed, 15 insertions(+), 2 deletions(-)

diff --git a/src/aliceVision/fuseCut/DelaunayGraphCut.hpp b/src/aliceVision/fuseCut/DelaunayGraphCut.hpp
index 2b6e693579..6ffc079391 100644
--- a/src/aliceVision/fuseCut/DelaunayGraphCut.hpp
+++ b/src/aliceVision/fuseCut/DelaunayGraphCut.hpp
@@ -97,12 +97,25 @@ class DelaunayGraphCut
     DelaunayGraphCut(mvsUtils::MultiViewParams* _mp);
     virtual ~DelaunayGraphCut();
 
-    /// Get absolute opposite vertex index
+    /**
+     * @brief Retrieve the global vertex index of the localVertexIndex of the facet.
+     * 
+     * @param f the facet
+     * @return the global vertex index
+     */
     inline VertexIndex getOppositeVertexIndex(const Facet& f) const
     {
         return _tetrahedralization->cell_vertex(f.cellIndex, f.localVertexIndex);
     }
-    /// Get absolute vertex index
+
+    /**
+     * @brief Retrieve the global vertex index of a vertex from a facet and an relative index 
+     * compared to the localVertexIndex of the facet.
+     * 
+     * @param f the facet
+     * @param i the relative index (relative to the localVertexIndex of the facet)
+     * @return the global vertex index
+     */
     inline VertexIndex getVertexIndex(const Facet& f, int i) const
     {
         return _tetrahedralization->cell_vertex(f.cellIndex, ((f.localVertexIndex + i + 1) % 4));

From eecc4dcc3061f129db649262b53ff9e51fabee86 Mon Sep 17 00:00:00 2001
From: Enguerrand DE SMET <desmet.enguerrand@gmail.com>
Date: Wed, 8 Jul 2020 12:57:23 +0200
Subject: [PATCH 13/28] [fuseCut] DelaunayGraphCut: Use directly
 isLineInTriangle (CRITIC)

Use directly isLineInTriangle in getFacetFromVertexOnTheRayToTheCam instead of rayCellIntersection because we come directly from a vertex and only the opposite facet of an adjacent cell can be intersected.

Stop the collision tests as soon as an intersection with a facet is found.
The configuration with convex geometry of the adjacent
 cells(tetrahedrons) to a vertex means that only two intersections can be found with a ray crossing this vertex.
When the intersection verifying the distance condition is found, it can be returned without iterating on the other cells.
---
 src/aliceVision/fuseCut/DelaunayGraphCut.cpp | 26 +++++++++++---------
 src/aliceVision/fuseCut/DelaunayGraphCut.hpp |  9 ++++++-
 2 files changed, 23 insertions(+), 12 deletions(-)

diff --git a/src/aliceVision/fuseCut/DelaunayGraphCut.cpp b/src/aliceVision/fuseCut/DelaunayGraphCut.cpp
index 4e771187be..d30f56161f 100644
--- a/src/aliceVision/fuseCut/DelaunayGraphCut.cpp
+++ b/src/aliceVision/fuseCut/DelaunayGraphCut.cpp
@@ -1368,39 +1368,43 @@ bool DelaunayGraphCut::rayCellIntersection(const Point3d& camCenter, const Point
     return true;
 }
 
-DelaunayGraphCut::Facet DelaunayGraphCut::getFacetFromVertexOnTheRayToTheCam(VertexIndex vertexIndex, int cam, bool nearestFarest) const
+DelaunayGraphCut::Facet DelaunayGraphCut::getFacetFromVertexOnTheRayToTheCam(VertexIndex globalVertexIndex, int cam, bool nearestFarest) const
 {
     if((cam < 0) || (cam >= mp->ncams))
     {
-        ALICEVISION_LOG_WARNING("Bad camId, cam: " << cam << ", ptid: " << vertexIndex);
+        ALICEVISION_LOG_WARNING("Bad camId, cam: " << cam << ", ptid: " << globalVertexIndex);
     }
-    const Point3d& p = _verticesCoords[vertexIndex];
+    const Point3d& p = _verticesCoords[globalVertexIndex];
 
     double currentDist = (mp->CArr[cam] - p).size(); // initialize currentDist to the distance from 3d point p to camera center
-    Facet facet;
 
-    for(CellIndex adjCellIndex : getNeighboringCellsByVertexIndex(vertexIndex)) // GEOGRAM: set_stores_cicl(true) required
+    for(CellIndex adjCellIndex : getNeighboringCellsByVertexIndex(globalVertexIndex)) // GEOGRAM: set_stores_cicl(true) required
     {
         if(isInfiniteCell(adjCellIndex))
             continue;
 
-        Facet outFacet;
+        // Get local vertex index
+        const VertexIndex localVertexIndex = _tetrahedralization->index(adjCellIndex, globalVertexIndex);
+
+        // Define facet
+        const Facet facet(adjCellIndex, localVertexIndex);
+
+        const std::array<const Point3d*, 3> facetPoints = getFacetsPoints(facet);
 
         Point3d intersectPt;
-        if(rayCellIntersection(mp->CArr[cam], p, adjCellIndex, outFacet, nearestFarest, intersectPt) == true)
+        const Point3d lineVect = (mp->CArr[cam] - p).normalize();
+        if(isLineInTriangle(&intersectPt, facetPoints[0], facetPoints[1], facetPoints[2], &p, &lineVect))
         {
             const double intersectDist = (mp->CArr[cam] - intersectPt).size();
             // if it's between the camera and the point if nearestFarest == true
             // if it's behind the point (from the camera) if nearestFarest == false
             if(nearestFarest ? (intersectDist < currentDist) : (intersectDist > currentDist))
             {
-                facet = outFacet;
-                currentDist = intersectDist;
+                return facet;
             }
-            // TODO FACA: maybe we can break and remove minDist?
         }
     }
-    return facet;
+    return Facet();
 }
 
 GEO::index_t DelaunayGraphCut::getFirstCellOnTheRayFromCamToThePoint(int cam, Point3d& p, Point3d& intersectPoint) const
diff --git a/src/aliceVision/fuseCut/DelaunayGraphCut.hpp b/src/aliceVision/fuseCut/DelaunayGraphCut.hpp
index 6ffc079391..c6ece5a7ed 100644
--- a/src/aliceVision/fuseCut/DelaunayGraphCut.hpp
+++ b/src/aliceVision/fuseCut/DelaunayGraphCut.hpp
@@ -121,6 +121,13 @@ class DelaunayGraphCut
         return _tetrahedralization->cell_vertex(f.cellIndex, ((f.localVertexIndex + i + 1) % 4));
     }
 
+    inline const std::array<const Point3d*, 3> getFacetsPoints(const Facet& f) const
+    {
+        return {&(_verticesCoords[getVertexIndex(f, 0)]),
+                &(_verticesCoords[getVertexIndex(f, 1)]),
+                &(_verticesCoords[getVertexIndex(f, 2)])};
+    }
+
     inline std::size_t getNbVertices() const
     {
         return _verticesAttr.size();
@@ -291,7 +298,7 @@ class DelaunayGraphCut
     bool rayCellIntersection(const Point3d& camCenter, const Point3d& p, GEO::index_t tetrahedronIndex, Facet& outFacet,
                              bool nearestFarest, Point3d& outIntersectPt) const;
 
-    Facet getFacetFromVertexOnTheRayToTheCam(VertexIndex vertexIndex, int cam, bool nearestFarest) const;
+    Facet getFacetFromVertexOnTheRayToTheCam(VertexIndex globalVertexIndex, int cam, bool nearestFarest) const;
     GEO::index_t getFirstCellOnTheRayFromCamToThePoint(int cam, Point3d& p, Point3d& lpi) const;
 
     float distFcn(float maxDist, float dist, float distFcnHeight) const;

From 2f13ccd102e09ddbce3dd606095f18899cd71179 Mon Sep 17 00:00:00 2001
From: Enguerrand DE SMET <desmet.enguerrand@gmail.com>
Date: Wed, 8 Jul 2020 15:51:59 +0200
Subject: [PATCH 14/28] [fuseCut] DelaunayGraphCut: remove unnecessary
 lastFinite variable

---
 src/aliceVision/fuseCut/DelaunayGraphCut.cpp | 6 ++----
 1 file changed, 2 insertions(+), 4 deletions(-)

diff --git a/src/aliceVision/fuseCut/DelaunayGraphCut.cpp b/src/aliceVision/fuseCut/DelaunayGraphCut.cpp
index d30f56161f..5e655c8339 100644
--- a/src/aliceVision/fuseCut/DelaunayGraphCut.cpp
+++ b/src/aliceVision/fuseCut/DelaunayGraphCut.cpp
@@ -1678,7 +1678,6 @@ void DelaunayGraphCut::fillGraphPartPtRc(int& out_nstepsFront, int& out_nstepsBe
         CellIndex ci = getFacetFromVertexOnTheRayToTheCam(vertexIndex, cam, nearestFarest).cellIndex;
 
         Point3d p = originPt;
-        CellIndex lastFinite = GEO::NO_CELL;
         bool ok = ci != GEO::NO_CELL;
         while(ok)
         {
@@ -1715,15 +1714,14 @@ void DelaunayGraphCut::fillGraphPartPtRc(int& out_nstepsFront, int& out_nstepsBe
                 if(ci == GEO::NO_CELL)
                     ok = false;
                 p = intersectPt;
-                lastFinite = ci;
             }
         }
 
         // get the outer tetrahedron of camera c for the ray to p = the last tetrahedron
-        if(lastFinite != GEO::NO_CELL)
+        if(ci != GEO::NO_CELL)
         {
 #pragma OMP_ATOMIC_WRITE
-            _cellsAttr[lastFinite].cellSWeight = (float)maxint;
+            _cellsAttr[ci].cellSWeight = (float)maxint;
         }
     }
 

From 5ebdfc39d466a74d6e9b837bcc7c2db08dc86887 Mon Sep 17 00:00:00 2001
From: Enguerrand DE SMET <desmet.enguerrand@gmail.com>
Date: Fri, 10 Jul 2020 17:45:15 +0200
Subject: [PATCH 15/28] [fuseCut] DelaunayGraphCut: Use directly
 isLineInTriangle (CRITIC)

Use directly isLineInTriangle in getFirstFacetOnTheRayFromCamToThePoint instead of rayCellIntersection .
---
 src/aliceVision/fuseCut/DelaunayGraphCut.cpp | 18 +++++++++++++-----
 src/aliceVision/fuseCut/DelaunayGraphCut.hpp |  2 +-
 2 files changed, 14 insertions(+), 6 deletions(-)

diff --git a/src/aliceVision/fuseCut/DelaunayGraphCut.cpp b/src/aliceVision/fuseCut/DelaunayGraphCut.cpp
index 5e655c8339..02500c71c5 100644
--- a/src/aliceVision/fuseCut/DelaunayGraphCut.cpp
+++ b/src/aliceVision/fuseCut/DelaunayGraphCut.cpp
@@ -1407,7 +1407,7 @@ DelaunayGraphCut::Facet DelaunayGraphCut::getFacetFromVertexOnTheRayToTheCam(Ver
     return Facet();
 }
 
-GEO::index_t DelaunayGraphCut::getFirstCellOnTheRayFromCamToThePoint(int cam, Point3d& p, Point3d& intersectPoint) const
+GEO::index_t DelaunayGraphCut::getFirstCellOnTheRayFromCamToThePoint(int cam, const Point3d& p, Point3d& intersectPt) const
 {
     int cam_vi = _camsVertexes[cam];
     Point3d camBehind = mp->CArr[cam] + (mp->CArr[cam] - p);
@@ -1420,14 +1420,22 @@ GEO::index_t DelaunayGraphCut::getFirstCellOnTheRayFromCamToThePoint(int cam, Po
         if(isInfiniteCell(adjCellIndex))
             continue;
 
-        Facet outFacet;
-        if(rayCellIntersection(camBehind, mp->CArr[cam], adjCellIndex, outFacet, false, intersectPoint))
+        // Get local vertex index
+        const VertexIndex localVertexIndex = _tetrahedralization->index(adjCellIndex, cam_vi);
+
+        // Define facet
+        const Facet facet(adjCellIndex, localVertexIndex);
+
+        const std::array<const Point3d*, 3> facetPoints = getFacetsPoints(facet);
+
+        const Point3d lineVect = (camBehind - mp->CArr[cam]).normalize();
+        if(isLineInTriangle(&intersectPt, facetPoints[0], facetPoints[1], facetPoints[2], &mp->CArr[cam], &lineVect))
         {
-            float dist = orientedPointPlaneDistance(intersectPoint, camBehind, dir);
+            const float dist = orientedPointPlaneDistance(intersectPt, camBehind, dir);
             if(dist > maxdist)
             {
                 maxdist = dist;
-                farestCell = adjCellIndex;
+                farestCell = facet.cellIndex;
             }
         }
     }
diff --git a/src/aliceVision/fuseCut/DelaunayGraphCut.hpp b/src/aliceVision/fuseCut/DelaunayGraphCut.hpp
index c6ece5a7ed..b407658ce6 100644
--- a/src/aliceVision/fuseCut/DelaunayGraphCut.hpp
+++ b/src/aliceVision/fuseCut/DelaunayGraphCut.hpp
@@ -299,7 +299,7 @@ class DelaunayGraphCut
                              bool nearestFarest, Point3d& outIntersectPt) const;
 
     Facet getFacetFromVertexOnTheRayToTheCam(VertexIndex globalVertexIndex, int cam, bool nearestFarest) const;
-    GEO::index_t getFirstCellOnTheRayFromCamToThePoint(int cam, Point3d& p, Point3d& lpi) const;
+    GEO::index_t getFirstCellOnTheRayFromCamToThePoint(int cam, const Point3d& p, Point3d& intersectPt) const;
 
     float distFcn(float maxDist, float dist, float distFcnHeight) const;
 

From 777abdcaa985654ec41b2b36ca4245e0116c0246 Mon Sep 17 00:00:00 2001
From: Enguerrand DE SMET <desmet.enguerrand@gmail.com>
Date: Fri, 10 Jul 2020 18:13:57 +0200
Subject: [PATCH 16/28] [fuseCut] DelaunayGraphCut: Consider previously
 intersected facet

Consider the previously intersected facet in rayCellIntersection instead of CellIndex in order to avoid later testing the collision with the previously intersected facet. (performance improvement)
---
 src/aliceVision/fuseCut/DelaunayGraphCut.cpp | 138 +++++++++----------
 src/aliceVision/fuseCut/DelaunayGraphCut.hpp |   4 +-
 2 files changed, 68 insertions(+), 74 deletions(-)

diff --git a/src/aliceVision/fuseCut/DelaunayGraphCut.cpp b/src/aliceVision/fuseCut/DelaunayGraphCut.cpp
index 02500c71c5..eb2bab31fa 100644
--- a/src/aliceVision/fuseCut/DelaunayGraphCut.cpp
+++ b/src/aliceVision/fuseCut/DelaunayGraphCut.cpp
@@ -1308,17 +1308,16 @@ void DelaunayGraphCut::removeSmallSegs(int minSegSize)
     initVertices();
 }
 
-bool DelaunayGraphCut::rayCellIntersection(const Point3d& camCenter, const Point3d& p, GEO::index_t tetrahedronIndex, Facet& outFacet,
+bool DelaunayGraphCut::rayCellIntersection(const Point3d& camCenter, const Point3d& p, const Facet& inFacet, Facet& outFacet,
                                            bool nearestFarest, Point3d& outIntersectPt) const
 {
     outIntersectPt = p; // important
     outFacet.cellIndex = GEO::NO_CELL;
     outFacet.localVertexIndex = GEO::NO_VERTEX;
 
+    const CellIndex tetrahedronIndex = inFacet.cellIndex;
     if(isInfiniteCell(tetrahedronIndex))
-    {
         return false;
-    }
 
     Point3d linePoint = p;
     Point3d lineVect = (camCenter - p).normalize();
@@ -1332,20 +1331,21 @@ bool DelaunayGraphCut::rayCellIntersection(const Point3d& camCenter, const Point
     const Point3d* D = &(_verticesCoords[_tetrahedralization->cell_vertex(tetrahedronIndex, 3)]);
 
     // All the facets of the tetrahedron
-    std::array<std::array<const Point3d*, 3>, 4> facets {{
+    std::array<std::array<const Point3d*, 3>, 4> facetsPoints {{
         {B, C, D}, // opposite vertex A, index 0
         {A, C, D}, // opposite vertex B, index 1
         {A, B, D}, // opposite vertex C, index 2
         {A, B, C}  // opposite vertex D, index 3
     }};
-    int oppositeVertexIndex = -1;
+    GEO::index_t oppositeVertexIndex = GEO::NO_VERTEX;
     Point3d intersectPt;
     double intersectDist;
     // Test all facets of the tetrahedron
-    for(int i = 0; i < 4; ++i)
+    for(int i = 0; i < facetsPoints.size(); ++i)
     {
-        const auto& facet= facets[i];
-        if(isLineInTriangle(&intersectPt, facet[0], facet[1], facet[2], &linePoint, &lineVect))
+
+        const auto& facetPoints = facetsPoints[i];
+        if(isLineInTriangle(&intersectPt, facetPoints[0], facetPoints[1], facetPoints[2], &linePoint, &lineVect))
         {
             intersectDist = (camCenter - intersectPt).size();
             // between the camera and the point if nearestFarest == true
@@ -1407,13 +1407,13 @@ DelaunayGraphCut::Facet DelaunayGraphCut::getFacetFromVertexOnTheRayToTheCam(Ver
     return Facet();
 }
 
-GEO::index_t DelaunayGraphCut::getFirstCellOnTheRayFromCamToThePoint(int cam, const Point3d& p, Point3d& intersectPt) const
+DelaunayGraphCut::Facet DelaunayGraphCut::getFirstFacetOnTheRayFromCamToThePoint(int cam, const Point3d& p, Point3d& intersectPt) const
 {
     int cam_vi = _camsVertexes[cam];
     Point3d camBehind = mp->CArr[cam] + (mp->CArr[cam] - p);
     Point3d dir = (p - mp->CArr[cam]).normalize();
     float maxdist = 0.0f;
-    GEO::index_t farestCell = GEO::NO_CELL;
+    Facet farestFacet;
 
     for(CellIndex adjCellIndex : getNeighboringCellsByVertexIndex(cam_vi)) // GEOGRAM: set_stores_cicl(true) required
     {
@@ -1435,12 +1435,11 @@ GEO::index_t DelaunayGraphCut::getFirstCellOnTheRayFromCamToThePoint(int cam, co
             if(dist > maxdist)
             {
                 maxdist = dist;
-                farestCell = facet.cellIndex;
+                farestFacet = facet;
             }
         }
     }
-
-    return farestCell;
+    return farestFacet;
 }
 
 float DelaunayGraphCut::distFcn(float maxDist, float dist, float distFcnHeight) const
@@ -1682,16 +1681,16 @@ void DelaunayGraphCut::fillGraphPartPtRc(int& out_nstepsFront, int& out_nstepsBe
         out_nstepsFront = 0;
         // true here mean nearest
         const bool nearestFarest = true;
-        // tetrahedron connected to the point p and which intersect the ray from camera c to point p
-        CellIndex ci = getFacetFromVertexOnTheRayToTheCam(vertexIndex, cam, nearestFarest).cellIndex;
+        // tetrahedron connected to the point define by the vertexIndex and which intersect the ray from camera
+        Facet facet = getFacetFromVertexOnTheRayToTheCam(vertexIndex, cam, nearestFarest);
 
         Point3d p = originPt;
-        bool ok = ci != GEO::NO_CELL;
+        bool ok = facet.cellIndex != GEO::NO_CELL;
         while(ok)
         {
             {
 #pragma OMP_ATOMIC_UPDATE
-                _cellsAttr[ci].out += weight;
+                _cellsAttr[facet.cellIndex].out += weight;
             }
 
             ++out_nstepsFront;
@@ -1704,7 +1703,7 @@ void DelaunayGraphCut::fillGraphPartPtRc(int& out_nstepsFront, int& out_nstepsBe
             // cam which is intersected with cam-p ray
             // true here mean nearest
             const bool nearestFarest = true;
-            if(!rayCellIntersection(mp->CArr[cam], p, ci, outFacet, nearestFarest, intersectPt))
+            if(!rayCellIntersection(mp->CArr[cam], p, facet, outFacet, nearestFarest, intersectPt))
             {
                 ok = false;
             }
@@ -1717,19 +1716,19 @@ void DelaunayGraphCut::fillGraphPartPtRc(int& out_nstepsFront, int& out_nstepsBe
                 }
 
                 // Take the mirror facet to iterate over the next cell
-                ci = mirrorFacet(outFacet).cellIndex;
+                facet = mirrorFacet(outFacet);
 
-                if(ci == GEO::NO_CELL)
+                if(facet.cellIndex == GEO::NO_CELL)
                     ok = false;
                 p = intersectPt;
             }
         }
 
         // get the outer tetrahedron of camera c for the ray to p = the last tetrahedron
-        if(ci != GEO::NO_CELL)
+        if(facet.cellIndex != GEO::NO_CELL)
         {
 #pragma OMP_ATOMIC_WRITE
-            _cellsAttr[ci].cellSWeight = (float)maxint;
+            _cellsAttr[facet.cellIndex].cellSWeight = (float)maxint;
         }
     }
 
@@ -1738,22 +1737,21 @@ void DelaunayGraphCut::fillGraphPartPtRc(int& out_nstepsFront, int& out_nstepsBe
         // get the tetrahedron next to point p on the ray from c
         // False here mean farest
         const bool nearestFarest = false;
-        Facet f1 = getFacetFromVertexOnTheRayToTheCam(vertexIndex, cam, nearestFarest);
+        Facet facet = getFacetFromVertexOnTheRayToTheCam(vertexIndex, cam, nearestFarest);
         Facet outFacet;
 
-        CellIndex ci = f1.cellIndex;
-        if(ci != GEO::NO_CELL)
+        if(facet.cellIndex != GEO::NO_CELL)
         {
 #pragma OMP_ATOMIC_UPDATE
-            _cellsAttr[ci].on += weight;
+            _cellsAttr[facet.cellIndex].on += weight;
         }
 
         Point3d p = originPt; // HAS TO BE HERE !!!
 
-        bool ok = (ci != GEO::NO_CELL) && allPoints;
+        bool ok = (facet.cellIndex != GEO::NO_CELL) && allPoints;
         while(ok)
         {
-            GC_cellInfo& c = _cellsAttr[ci];
+            GC_cellInfo& c = _cellsAttr[facet.cellIndex];
             {
                 if(behind)
                 {
@@ -1769,11 +1767,11 @@ void DelaunayGraphCut::fillGraphPartPtRc(int& out_nstepsFront, int& out_nstepsBe
 
             Point3d intersectPt;
 
-            // Intersection with the next facet in the current tetrahedron (ci) in order to find the cell farest to the
+            // Intersection with the next facet in the current tetrahedron (facet.cellIndex) in order to find the cell farest to the
             // cam which is intersected with cam-p ray
             // False here mean farest
             const bool nearestFarest = false;
-            if(!rayCellIntersection(mp->CArr[cam], p, ci, outFacet, nearestFarest, intersectPt) ||
+            if(!rayCellIntersection(mp->CArr[cam], p, facet, outFacet, nearestFarest, intersectPt) ||
                ((originPt - p).size() >= maxDist) || (!allPoints))
             {
                 ok = false;
@@ -1788,9 +1786,8 @@ void DelaunayGraphCut::fillGraphPartPtRc(int& out_nstepsFront, int& out_nstepsBe
 
 
                 // Take the mirror facet to iterate over the next cell
-                const Facet mFacet = mirrorFacet(outFacet);
-                ci = mFacet.cellIndex;
-                if(ci == GEO::NO_CELL)
+                facet = mirrorFacet(outFacet);
+                if(facet.cellIndex == GEO::NO_CELL)
                 {
                     ok = false;
                 }
@@ -1798,7 +1795,7 @@ void DelaunayGraphCut::fillGraphPartPtRc(int& out_nstepsFront, int& out_nstepsBe
                 {
                     const float dist = distFcn(maxDist, (originPt - p).size(), distFcnHeight);
 #pragma OMP_ATOMIC_UPDATE
-                    _cellsAttr[ci].gEdgeVisWeight[mFacet.localVertexIndex] += weight * dist;
+                    _cellsAttr[facet.cellIndex].gEdgeVisWeight[facet.localVertexIndex] += weight * dist;
                 }
                 p = intersectPt;
             }
@@ -1807,10 +1804,10 @@ void DelaunayGraphCut::fillGraphPartPtRc(int& out_nstepsFront, int& out_nstepsBe
         // cv: is the tetrahedron in distance 2*sigma behind the point p in the direction of the camera c (called Lcp in the paper)
         if(!behind)
         {
-            if(ci != GEO::NO_CELL)
+            if(facet.cellIndex != GEO::NO_CELL)
             {
 #pragma OMP_ATOMIC_UPDATE
-                _cellsAttr[ci].cellTWeight += weight;
+                _cellsAttr[facet.cellIndex].cellTWeight += weight;
             }
         }
     }
@@ -1854,13 +1851,12 @@ void DelaunayGraphCut::forceTedgesByGradientCVPR11(bool fixesSigma, float nPixel
 
             // get the tetrahedron next to point p on the ray from c
             // False here mean farest
-            Facet f1 = getFacetFromVertexOnTheRayToTheCam(vi, cam, false);
+            Facet facet = getFacetFromVertexOnTheRayToTheCam(vi, cam, false);
 
-            if((fFirst.cellIndex != GEO::NO_CELL) && (f1.cellIndex != GEO::NO_CELL) && (!isInfiniteCell(f1.cellIndex)))
+            if((fFirst.cellIndex != GEO::NO_CELL) && (facet.cellIndex != GEO::NO_CELL) && (!isInfiniteCell(facet.cellIndex)))
             {
                 float eFirst = _cellsAttr[fFirst.cellIndex].out;
 
-                CellIndex ci = f1.cellIndex;
                 Point3d p = originPt; // HAS TO BE HERE !!!
                 float maxDist = nPixelSizeBehind * mp->getCamPixelSize(p, cam);
                 if(fixesSigma)
@@ -1868,7 +1864,7 @@ void DelaunayGraphCut::forceTedgesByGradientCVPR11(bool fixesSigma, float nPixel
                     maxDist = nPixelSizeBehind;
                 }
 
-                bool ok = (ci != GEO::NO_CELL);
+                bool ok = (facet.cellIndex != GEO::NO_CELL);
                 while(ok)
                 {
                     Point3d pold = p;
@@ -1879,7 +1875,7 @@ void DelaunayGraphCut::forceTedgesByGradientCVPR11(bool fixesSigma, float nPixel
                     // to the cam which is intersected with cam-p ray
                     // False here mean farest
                     const bool nearestFarest = false;
-                    if(!rayCellIntersection(mp->CArr[cam], p, ci, outFacet, nearestFarest, intersectPt) ||
+                    if(!rayCellIntersection(mp->CArr[cam], p, facet, outFacet, nearestFarest, intersectPt) ||
                        ((originPt - pold).size() >= maxDist))
                     {
                         ok = false;
@@ -1887,20 +1883,20 @@ void DelaunayGraphCut::forceTedgesByGradientCVPR11(bool fixesSigma, float nPixel
                     else
                     {
                         // Take the mirror facet to iterate over the next cell
-                        ci = mirrorFacet(outFacet).cellIndex;
-                        if(ci == GEO::NO_CELL)
+                        facet = mirrorFacet(outFacet);
+                        if(facet.cellIndex == GEO::NO_CELL)
                             ok = false;
                         p = intersectPt;
                     }
                 }
 
-                if(ci != GEO::NO_CELL)
+                if(facet.cellIndex != GEO::NO_CELL)
                 {
-                    float eLast = _cellsAttr[ci].out;
+                    const float eLast = _cellsAttr[facet.cellIndex].out;
                     if((eFirst > eLast) && (eFirst < beta) && (eLast / eFirst < delta))
                     {
 #pragma OMP_ATOMIC_UPDATE
-                        _cellsAttr[ci].on += (eFirst - eLast);
+                        _cellsAttr[facet.cellIndex].on += (eFirst - eLast);
                     }
                 }
             }
@@ -1992,14 +1988,14 @@ void DelaunayGraphCut::forceTedgesByGradientIJCV(bool fixesSigma, float nPixelSi
             {
                 // True here mean nearest
                 const bool nearestFarest = true;
-                CellIndex ci = getFacetFromVertexOnTheRayToTheCam(vi, cam, nearestFarest).cellIndex;
+                Facet facet = getFacetFromVertexOnTheRayToTheCam(vi, cam, nearestFarest);
                 Point3d p = originPt; // HAS TO BE HERE !!!
-                bool ok = (ci != GEO::NO_CELL);
+                bool ok = (facet.cellIndex != GEO::NO_CELL);
                 while(ok)
                 {
                     ++nstepsFront;
 
-                    const GC_cellInfo& c = _cellsAttr[ci];
+                    const GC_cellInfo& c = _cellsAttr[facet.cellIndex];
                     if((p - originPt).size() > nsigmaFrontSilentPart * maxDist) // (p-originPt).size() > 2 * sigma
                     {
                         minJump = std::min(minJump, c.out);
@@ -2013,20 +2009,20 @@ void DelaunayGraphCut::forceTedgesByGradientIJCV(bool fixesSigma, float nPixelSi
 
                     Facet outFacet;
                     Point3d intersectPt;
-                    // Intersection with the next facet in the current tetrahedron (ci) in order to find the cell nearest
+                    // Intersection with the next facet in the current tetrahedron (facet.cellIndex) in order to find the cell nearest
                     // to the cam which is intersected with cam-p ray
                     // True here mean nearest
                     const bool nearestFarest = true;
                     if(((p - originPt).size() > (nsigmaJumpPart + nsigmaFrontSilentPart) * maxDist) || // (2 + 2) * sigma
-                       !rayCellIntersection(mp->CArr[cam], p, ci, outFacet, nearestFarest, intersectPt))
+                       !rayCellIntersection(mp->CArr[cam], p, facet, outFacet, nearestFarest, intersectPt))
                     {
                         ok = false;
                     }
                     else
                     {
                         // Take the mirror facet to iterate over the next cell
-                        ci = mirrorFacet(outFacet).cellIndex;
-                        if(ci == GEO::NO_CELL)
+                        facet = mirrorFacet(outFacet);
+                        if(facet.cellIndex == GEO::NO_CELL)
                             ok = false;
                         p = intersectPt;
                     }
@@ -2036,18 +2032,18 @@ void DelaunayGraphCut::forceTedgesByGradientIJCV(bool fixesSigma, float nPixelSi
             {
                 // False here mean farest
                 const bool nearestFarest = false;
-                CellIndex ci = getFacetFromVertexOnTheRayToTheCam(vi, cam, nearestFarest).cellIndex; // T1
+                Facet facet = getFacetFromVertexOnTheRayToTheCam(vi, cam, nearestFarest);
                 Point3d p = originPt; // HAS TO BE HERE !!!
-                bool ok = (ci != GEO::NO_CELL);
+                bool ok = (facet.cellIndex != GEO::NO_CELL);
                 if(ok)
                 {
-                    midSilent = _cellsAttr[ci].out;
+                    midSilent = _cellsAttr[facet.cellIndex].out;
                 }
 
                 while(ok)
                 {
                     nstepsBehind++;
-                    const GC_cellInfo& c = _cellsAttr[ci];
+                    const GC_cellInfo& c = _cellsAttr[facet.cellIndex];
 
                     minSilent = std::min(minSilent, c.out);
                     maxSilent = std::max(maxSilent, c.out);
@@ -2060,21 +2056,21 @@ void DelaunayGraphCut::forceTedgesByGradientIJCV(bool fixesSigma, float nPixelSi
                     // False here mean farest
                     const bool nearestFarest = false;
                     if(((p - originPt).size() > nsigmaBackSilentPart * maxDist) || // (p-originPt).size() > 2 * sigma
-                       !rayCellIntersection(mp->CArr[cam], p, ci, outFacet, nearestFarest, intersectPt))
+                       !rayCellIntersection(mp->CArr[cam], p, facet, outFacet, nearestFarest, intersectPt))
                     {
                         ok = false;
                     }
                     else
                     {
                         // Take the mirror facet to iterate over the next cell
-                        ci = mirrorFacet(outFacet).cellIndex;
-                        if(ci == GEO::NO_CELL)
+                        facet = mirrorFacet(outFacet);
+                        if(facet.cellIndex == GEO::NO_CELL)
                             ok = false;
                         p = intersectPt;
                     }
                 }
 
-                if(ci != GEO::NO_CELL)
+                if(facet.cellIndex != GEO::NO_CELL)
                 {
                     // Equation 6 in paper
                     //   (g / B) < k_rel
@@ -2097,7 +2093,7 @@ void DelaunayGraphCut::forceTedgesByGradientIJCV(bool fixesSigma, float nPixelSi
                         //(maxSilent-minSilent<maxSilentPartRange))
                     {
 #pragma OMP_ATOMIC_UPDATE
-                        _cellsAttr[ci].on += (maxJump - midSilent);
+                        _cellsAttr[facet.cellIndex].on += (maxJump - midSilent);
                     }
                 }
             }
@@ -2185,21 +2181,20 @@ void DelaunayGraphCut::updateGraphFromTmpPtsCamsHexahRC(int rc, Point3d hexah[8]
             nin++;
             Point3d intersectPt = pt;
             Point3d camBehind = mp->CArr[rc] + (mp->CArr[rc] - pt);
-            CellIndex ci = getFirstCellOnTheRayFromCamToThePoint(rc, pt, intersectPt);
-            if(ci != GEO::NO_CELL)
+            Facet facet = getFirstFacetOnTheRayFromCamToThePoint(rc, pt, intersectPt);
+            if(facet.cellIndex != GEO::NO_CELL)
             {
                 // update weights on the sp-cam half line
-                CellIndex tmp_ci = ci;
                 Point3d p = mp->CArr[rc];
                 // _cellsAttr[ci].cellSWeight = (float)maxint;
 
-                bool ok = tmp_ci != GEO::NO_CELL;
+                bool ok = facet.cellIndex != GEO::NO_CELL;
 
                 while(ok)
                 {
                     {
 #pragma OMP_ATOMIC_UPDATE
-                        _cellsAttr[tmp_ci].out += weight;
+                        _cellsAttr[facet.cellIndex].out += weight;
                     }
 
                     Facet outFacet;
@@ -2209,23 +2204,22 @@ void DelaunayGraphCut::updateGraphFromTmpPtsCamsHexahRC(int rc, Point3d hexah[8]
                     // to the cam which is intersected with cam-p ray
                     // False here mean farest
                     const bool nearestFarest = false;
-                    if(!rayCellIntersection(camBehind, p, tmp_ci, outFacet, nearestFarest, intersectPt))
+                    if(!rayCellIntersection(camBehind, p, facet, outFacet, nearestFarest, intersectPt))
                     {
                         ok = false;
                     }
                     else
                     {
                         // Take the mirror facet to iterate over the next cell
-                        const Facet mFacet = mirrorFacet(outFacet);
-                        tmp_ci = mFacet.cellIndex;
-                        if(tmp_ci == GEO::NO_CELL)
+                        facet = mirrorFacet(outFacet);
+                        if(facet.cellIndex == GEO::NO_CELL)
                         {
                             ok = false;
                         }
                         else
                         {
 #pragma OMP_ATOMIC_UPDATE
-                            _cellsAttr[tmp_ci].gEdgeVisWeight[mFacet.localVertexIndex] += weight;
+                            _cellsAttr[facet.cellIndex].gEdgeVisWeight[facet.localVertexIndex] += weight;
                         }
                         p = intersectPt;
                         ++nwup;
diff --git a/src/aliceVision/fuseCut/DelaunayGraphCut.hpp b/src/aliceVision/fuseCut/DelaunayGraphCut.hpp
index b407658ce6..a8b6a1a269 100644
--- a/src/aliceVision/fuseCut/DelaunayGraphCut.hpp
+++ b/src/aliceVision/fuseCut/DelaunayGraphCut.hpp
@@ -295,11 +295,11 @@ class DelaunayGraphCut
     void computeVerticesSegSize(bool allPoints, float alpha = 0.0f);
     void removeSmallSegs(int minSegSize);
 
-    bool rayCellIntersection(const Point3d& camCenter, const Point3d& p, GEO::index_t tetrahedronIndex, Facet& outFacet,
+    bool rayCellIntersection(const Point3d& camCenter, const Point3d& p, const Facet& inFacet, Facet& outFacet,
                              bool nearestFarest, Point3d& outIntersectPt) const;
 
     Facet getFacetFromVertexOnTheRayToTheCam(VertexIndex globalVertexIndex, int cam, bool nearestFarest) const;
-    GEO::index_t getFirstCellOnTheRayFromCamToThePoint(int cam, const Point3d& p, Point3d& intersectPt) const;
+    Facet getFirstFacetOnTheRayFromCamToThePoint(int cam, const Point3d& p, Point3d& intersectPt) const;
 
     float distFcn(float maxDist, float dist, float distFcnHeight) const;
 

From 005c02e3c8df566fd732d64e3ed85a723a687170 Mon Sep 17 00:00:00 2001
From: Enguerrand DE SMET <desmet.enguerrand@gmail.com>
Date: Mon, 13 Jul 2020 13:58:30 +0200
Subject: [PATCH 17/28] [fuseCut] Remove regular grid option

---
 src/aliceVision/fuseCut/DelaunayGraphCut.cpp  | 250 +-----------------
 src/aliceVision/fuseCut/DelaunayGraphCut.hpp  |  11 +-
 .../fuseCut/ReconstructionPlan.cpp            |  67 -----
 .../fuseCut/ReconstructionPlan.hpp            |   1 -
 src/software/pipeline/main_meshing.cpp        | 102 +------
 5 files changed, 12 insertions(+), 419 deletions(-)

diff --git a/src/aliceVision/fuseCut/DelaunayGraphCut.cpp b/src/aliceVision/fuseCut/DelaunayGraphCut.cpp
index eb2bab31fa..bac47b0870 100644
--- a/src/aliceVision/fuseCut/DelaunayGraphCut.cpp
+++ b/src/aliceVision/fuseCut/DelaunayGraphCut.cpp
@@ -1128,74 +1128,6 @@ void DelaunayGraphCut::fuseFromDepthMaps(const StaticVector<int>& cams, const Po
     ALICEVISION_LOG_INFO("fuseFromDepthMaps done: " << _verticesCoords.size() << " points created.");
 }
 
-void DelaunayGraphCut::loadPrecomputedDensePoints(const StaticVector<int>* voxelsIds, const Point3d voxel[8], VoxelsGrid* ls)
-{
-    Point3d cgpt;
-    int ncgpt = 0;
-
-    bool doFilterOctreeTracks = mp->userParams.get<bool>("LargeScale.doFilterOctreeTracks", true);
-    int minNumOfConsistentCams = mp->userParams.get<int>("filter.minNumOfConsistentCams", 2);
-
-    // add points from voxel
-    for(int i = 0; i < voxelsIds->size(); i++)
-    {
-        ALICEVISION_LOG_INFO("Add points from voxel " << i << " of " << voxelsIds->size() << ".");
-
-        std::string folderName = ls->getVoxelFolderName((*voxelsIds)[i]);
-        std::string fileNameTracksPts = folderName + "tracksGridPts.bin";
-        std::string fileNameTracksPtsCams = folderName + "tracksGridPtsCams.bin";
-
-        if(mvsUtils::FileExists(fileNameTracksPts))
-        {
-            StaticVector<Point3d>* tracksPoints = loadArrayFromFile<Point3d>(fileNameTracksPts);
-            StaticVector<StaticVector<Pixel>*>* tracksPointsCams =
-                loadArrayOfArraysFromFile<Pixel>(fileNameTracksPtsCams);
-
-            long t1 = mvsUtils::initEstimate();
-            for(int j = 0; j < tracksPoints->size(); j++)
-            {
-                Point3d tp = (*tracksPoints)[j];
-                StaticVector<Pixel>* cams = (*tracksPointsCams)[j];
-                if((mvsUtils::isPointInHexahedron(tp, voxel)) && (cams != nullptr) &&
-                   ((!doFilterOctreeTracks) || ((doFilterOctreeTracks) && (cams->size() >= minNumOfConsistentCams))))
-                {
-                    cgpt = cgpt + tp;
-                    ncgpt++;
-                    _verticesCoords.push_back(tp);
-
-                    GC_vertexInfo newv;
-                    newv.nrc = 0;
-                    newv.segSize = 0;
-                    newv.segId = -1;
-
-                    newv.cams.reserve(cams->size());
-
-                    for(int c = 0; c < cams->size(); c++)
-                    {
-                        int rc = (*cams)[c].x;
-                        int nptsrc = (*cams)[c].y;
-                        newv.cams.push_back(rc);
-                        newv.nrc += nptsrc;
-                    }
-
-                    _verticesAttr.push_back(newv);
-                }
-
-                mvsUtils::printfEstimate(j, tracksPoints->size(), t1);
-            } // for j
-            mvsUtils::finishEstimate();
-
-            // delete randIds;
-
-            delete tracksPoints;
-            deleteArrayOfArrays<Pixel>(&tracksPointsCams);
-        } // if fileexists
-    }
-
-    ALICEVISION_LOG_INFO("Dense points loaded.\n");
-}
-
-
 void DelaunayGraphCut::computeVerticesSegSize(bool allPoints, float alpha) // allPoints=true, alpha=0
 {
     ALICEVISION_LOG_DEBUG("creating universe");
@@ -2123,121 +2055,6 @@ void DelaunayGraphCut::forceTedgesByGradientIJCV(bool fixesSigma, float nPixelSi
     mvsUtils::printfElapsedTime(t2, "t-edges forced: ");
 }
 
-void DelaunayGraphCut::updateGraphFromTmpPtsCamsHexah(const StaticVector<int>& incams, Point3d hexah[8],
-                                                    std::string tmpCamsPtsFolderName, bool labatutWeights,
-                                                    float distFcnHeight)
-{
-    ALICEVISION_LOG_INFO("Updating: LSC.");
-
-#pragma omp parallel for
-    for(int c = 0; c < incams.size(); c++)
-    {
-        int rc = incams[c];
-        std::string camPtsFileName;
-        camPtsFileName = tmpCamsPtsFolderName + "camPtsGrid_" + mvsUtils::num2strFourDecimal(rc) + ".bin";
-        ALICEVISION_LOG_INFO("File exist? " << camPtsFileName);
-        if(mvsUtils::FileExists(camPtsFileName))
-        {
-            updateGraphFromTmpPtsCamsHexahRC(rc, hexah, tmpCamsPtsFolderName, labatutWeights, distFcnHeight);
-        }
-    } // for c
-}
-
-void DelaunayGraphCut::updateGraphFromTmpPtsCamsHexahRC(int rc, Point3d hexah[8], std::string tmpCamsPtsFolderName,
-                                                      bool labatutWeights, float  /*distFcnHeight*/)
-{
-    ALICEVISION_LOG_INFO("DelaunayGraphCut::updateGraphFromTmpPtsCamsHexahRC: rc: " << rc);
-
-    // fill edges
-    int nin = 0;
-    int nout = 0;
-    int cnull = 0;
-    int nwup = 0;
-
-    std::string camPtsFileName;
-    camPtsFileName = tmpCamsPtsFolderName + "camPtsGrid_" + mvsUtils::num2strFourDecimal(rc) + ".bin";
-
-    bool doFilterOctreeTracks = mp->userParams.get<bool>("LargeScale.doFilterOctreeTracks", true);
-    int minNumOfConsistentCams = mp->userParams.get<int>("filter.minNumOfConsistentCams", 2);
-
-    /////////////////////////////////////////////////////////////////////////////////
-    FILE* f = fopen(camPtsFileName.c_str(), "rb");
-    while(feof(f) == 0)
-    {
-        GC_camVertexInfo pnt;
-        pnt.freadinfo(f);
-        // printf("%f, %i, %f, %f,
-        // %f\n",pnt.sim,pnt.nrc,pnt.point.x,pnt.point.y,pnt.point.z);
-        Point3d pt = pnt.point;
-
-        float weight = weightFcn((float)pnt.nrc, labatutWeights, (float)pnt.ncams);
-
-        nout++;
-
-        StaticVector<Point3d>* lshi = mvsUtils::lineSegmentHexahedronIntersection(pt, mp->CArr[rc], hexah);
-        if((!mvsUtils::isPointInHexahedron(pt, hexah)) && (lshi->size() >= 1) && (feof(f) == 0) &&
-           ((!doFilterOctreeTracks) || ((doFilterOctreeTracks) && (pnt.ncams >= minNumOfConsistentCams))))
-        {
-            nin++;
-            Point3d intersectPt = pt;
-            Point3d camBehind = mp->CArr[rc] + (mp->CArr[rc] - pt);
-            Facet facet = getFirstFacetOnTheRayFromCamToThePoint(rc, pt, intersectPt);
-            if(facet.cellIndex != GEO::NO_CELL)
-            {
-                // update weights on the sp-cam half line
-                Point3d p = mp->CArr[rc];
-                // _cellsAttr[ci].cellSWeight = (float)maxint;
-
-                bool ok = facet.cellIndex != GEO::NO_CELL;
-
-                while(ok)
-                {
-                    {
-#pragma OMP_ATOMIC_UPDATE
-                        _cellsAttr[facet.cellIndex].out += weight;
-                    }
-
-                    Facet outFacet;
-                    Point3d intersectPt;
-
-                    // Intersection with the next facet in the current tetrahedron (ci) in order to find the cell farest
-                    // to the cam which is intersected with cam-p ray
-                    // False here mean farest
-                    const bool nearestFarest = false;
-                    if(!rayCellIntersection(camBehind, p, facet, outFacet, nearestFarest, intersectPt))
-                    {
-                        ok = false;
-                    }
-                    else
-                    {
-                        // Take the mirror facet to iterate over the next cell
-                        facet = mirrorFacet(outFacet);
-                        if(facet.cellIndex == GEO::NO_CELL)
-                        {
-                            ok = false;
-                        }
-                        else
-                        {
-#pragma OMP_ATOMIC_UPDATE
-                            _cellsAttr[facet.cellIndex].gEdgeVisWeight[facet.localVertexIndex] += weight;
-                        }
-                        p = intersectPt;
-                        ++nwup;
-                    }
-                }
-            }
-            else
-            {
-                ++cnull;
-            }
-        }
-        delete lshi;
-    }
-    fclose(f);
-
-    // ALICEVISION_LOG_DEBUG("in: " << nin, ", cnull:" << cnull << " nwup: " << nwup << ", out: " << nout);
-}
-
 int DelaunayGraphCut::setIsOnSurface()
 {
     // set is on surface
@@ -2451,31 +2268,6 @@ void DelaunayGraphCut::invertFullStatusForSmallLabels()
     delete buff;
 }
 
-void DelaunayGraphCut::createDensePointCloudFromPrecomputedDensePoints(Point3d hexah[8], const StaticVector<int>& cams, StaticVector<int>* voxelsIds, VoxelsGrid* ls)
-{
-  // Load tracks
-  ALICEVISION_LOG_INFO("Creating delaunay tetrahedralization from depth maps voxel");
-
-  float minDist = hexah ? (hexah[0] - hexah[1]).size() / 1000.0f : 0.00001f;
-
-  // add points for cam centers
-  addPointsFromCameraCenters(cams, minDist);
-
-  // add 6 points to prevent singularities
-  addPointsToPreventSingularities(hexah, minDist);
-
-  loadPrecomputedDensePoints(voxelsIds, hexah, ls);
-
-  // initialize random seed
-  srand(time(nullptr));
-
-  int nGridHelperVolumePointsDim = mp->userParams.get<int>("LargeScale.nGridHelperVolumePointsDim", 10);
-
-  // add volume points to prevent singularities
-  addHelperPoints(nGridHelperVolumePointsDim, hexah, minDist);
-}
-
-
 void DelaunayGraphCut::createDensePointCloud(Point3d hexah[8], const StaticVector<int>& cams, const sfmData::SfMData* sfmData, const FuseParams* depthMapsFuseParams)
 {
   assert(sfmData != nullptr || depthMapsFuseParams != nullptr);
@@ -2507,7 +2299,7 @@ void DelaunayGraphCut::createDensePointCloud(Point3d hexah[8], const StaticVecto
   addHelperPoints(nGridHelperVolumePointsDim, hexah, minDist);
 }
 
-void DelaunayGraphCut::createGraphCut(Point3d hexah[8], const StaticVector<int>& cams, VoxelsGrid* ls, const std::string& folderName, const std::string& tmpCamsPtsFolderName, bool removeSmallSegments, const Point3d& spaceSteps)
+void DelaunayGraphCut::createGraphCut(Point3d hexah[8], const StaticVector<int>& cams, const std::string& folderName, const std::string& tmpCamsPtsFolderName, bool removeSmallSegments, const Point3d& spaceSteps)
 {
   initVertices();
 
@@ -2520,9 +2312,7 @@ void DelaunayGraphCut::createGraphCut(Point3d hexah[8], const StaticVector<int>&
   if(removeSmallSegments) // false
     removeSmallSegs(2500); // TODO FACA: to decide
 
-  const bool updateLSC = ls ? mp->userParams.get<bool>("LargeScale.updateLSC", true) : false;
-
-  reconstructExpetiments(cams, folderName, updateLSC,
+  reconstructExpetiments(cams, folderName,
                          hexah, tmpCamsPtsFolderName,
                          spaceSteps);
 
@@ -2647,7 +2437,7 @@ void DelaunayGraphCut::maxflow()
 }
 
 void DelaunayGraphCut::reconstructExpetiments(const StaticVector<int>& cams, const std::string& folderName,
-                                            bool update, Point3d hexahInflated[8], const std::string& tmpCamsPtsFolderName,
+                                            Point3d hexahInflated[8], const std::string& tmpCamsPtsFolderName,
                                             const Point3d& spaceSteps)
 {
     ALICEVISION_LOG_INFO("DelaunayGraphCut::reconstructExpetiments");
@@ -2655,12 +2445,8 @@ void DelaunayGraphCut::reconstructExpetiments(const StaticVector<int>& cams, con
 
     long t1;
 
-    bool fixesSigma = (update && spaceSteps.size() != 0.0);
-    float sigma = (float)mp->userParams.get<double>("delaunaycut.sigma", 4.0f); // TODO FACA: 2 or 4?
-    if(fixesSigma)
-        sigma *= spaceSteps.size();
+    const float sigma = (float)mp->userParams.get<double>("delaunaycut.sigma", 4.0f); // TODO FACA: 2 or 4?
 
-    ALICEVISION_LOG_INFO("fixesSigma: " << fixesSigma);
     ALICEVISION_LOG_INFO("sigma: " << sigma);
 
     // 0 for distFcn equals 1 all the time
@@ -2681,14 +2467,8 @@ void DelaunayGraphCut::reconstructExpetiments(const StaticVector<int>& cams, con
 
         ALICEVISION_LOG_INFO("Jancosek CVPR 2011 method ( delta*100 = " << static_cast<int>(delta * 100.0f) << "):");
 
-        fillGraph(fixesSigma, sigma, true, false, false, true, distFcnHeight);
+        fillGraph(false, sigma, true, false, false, true, distFcnHeight);
 
-        if(update)
-        {
-            t1 = clock();
-            updateGraphFromTmpPtsCamsHexah(cams, hexahInflated, tmpCamsPtsFolderName, false, distFcnHeight);
-            mvsUtils::printfElapsedTime(t1);
-        }
         addToInfiniteSw((float)maxint);
 
         if(saveTemporaryBinFiles)
@@ -2696,7 +2476,7 @@ void DelaunayGraphCut::reconstructExpetiments(const StaticVector<int>& cams, con
 
         if((delta > 0.0f) && (delta < 1.0f))
         {
-            forceTedgesByGradientCVPR11(fixesSigma, sigma);
+            forceTedgesByGradientCVPR11(false, sigma);
         }
 
         if(saveTemporaryBinFiles)
@@ -2715,14 +2495,8 @@ void DelaunayGraphCut::reconstructExpetiments(const StaticVector<int>& cams, con
         ALICEVISION_LOG_INFO("Jancosek IJCV method ( delta*100 = " << static_cast<int>(delta * 100.0f) << " ): ");
 
         // compute weights on edge between tetrahedra
-        fillGraph(fixesSigma, sigma, true, false, false, true, distFcnHeight);
+        fillGraph(false, sigma, true, false, false, true, distFcnHeight);
 
-        if(update) // true by default
-        {
-            t1 = clock();
-            updateGraphFromTmpPtsCamsHexah(cams, hexahInflated, tmpCamsPtsFolderName, false, distFcnHeight);
-            mvsUtils::printfElapsedTime(t1);
-        }
         addToInfiniteSw((float)maxint);
 
         if(saveTemporaryBinFiles)
@@ -2730,7 +2504,7 @@ void DelaunayGraphCut::reconstructExpetiments(const StaticVector<int>& cams, con
 
         if((delta > 0.0f) && (delta < 1.0f))
         {
-            forceTedgesByGradientIJCV(fixesSigma, sigma);
+            forceTedgesByGradientIJCV(false, sigma);
         }
 
         if(saveTemporaryBinFiles)
@@ -2745,13 +2519,7 @@ void DelaunayGraphCut::reconstructExpetiments(const StaticVector<int>& cams, con
     if(labatutCFG09)
     {
         ALICEVISION_LOG_INFO("Labatut CFG 2009 method:");
-        fillGraph(fixesSigma, sigma, true, false, true, true, distFcnHeight);
-        if(update)
-        {
-            t1 = clock();
-            updateGraphFromTmpPtsCamsHexah(cams, hexahInflated, tmpCamsPtsFolderName, distFcnHeight != 0.0f);
-            mvsUtils::printfElapsedTime(t1);
-        }
+        fillGraph(false, sigma, true, false, true, true, distFcnHeight);
 
         if(saveTemporaryBinFiles)
             saveDhInfo(folderName + "delaunayTriangulationInfoInit.bin");
diff --git a/src/aliceVision/fuseCut/DelaunayGraphCut.hpp b/src/aliceVision/fuseCut/DelaunayGraphCut.hpp
index a8b6a1a269..23a97ad08e 100644
--- a/src/aliceVision/fuseCut/DelaunayGraphCut.hpp
+++ b/src/aliceVision/fuseCut/DelaunayGraphCut.hpp
@@ -290,7 +290,6 @@ class DelaunayGraphCut
     void addHelperPoints(int nGridHelperVolumePointsDim, const Point3d Voxel[8], float minDist);
 
     void fuseFromDepthMaps(const StaticVector<int>& cams, const Point3d voxel[8], const FuseParams& params);
-    void loadPrecomputedDensePoints(const StaticVector<int>* voxelsIds, const Point3d voxel[8], VoxelsGrid* ls);
 
     void computeVerticesSegSize(bool allPoints, float alpha = 0.0f);
     void removeSmallSegs(int minSegSize);
@@ -320,11 +319,6 @@ class DelaunayGraphCut
     void forceTedgesByGradientCVPR11(bool fixesSigma, float nPixelSizeBehind);
     void forceTedgesByGradientIJCV(bool fixesSigma, float nPixelSizeBehind);
 
-    void updateGraphFromTmpPtsCamsHexah(const StaticVector<int>& incams, Point3d hexah[8], std::string tmpCamsPtsFolderName,
-                                        bool labatutWeights, float distFcnHeight = 0.0f);
-    void updateGraphFromTmpPtsCamsHexahRC(int rc, Point3d hexah[8], std::string tmpCamsPtsFolderName,
-                                          bool labatutWeights, float distFcnHeight);
-
     int setIsOnSurface();
 
     void addToInfiniteSw(float sW);
@@ -336,14 +330,13 @@ class DelaunayGraphCut
     void maxflow();
 
     void reconstructExpetiments(const StaticVector<int>& cams, const std::string& folderName,
-                                bool update, Point3d hexahInflated[8], const std::string& tmpCamsPtsFolderName,
+                                Point3d hexahInflated[8], const std::string& tmpCamsPtsFolderName,
                                 const Point3d& spaceSteps);
 
 
     void createDensePointCloud(Point3d hexah[8], const StaticVector<int>& cams, const sfmData::SfMData* sfmData, const FuseParams* depthMapsFuseParams);
-    void createDensePointCloudFromPrecomputedDensePoints(Point3d hexah[8], const StaticVector<int>& cams, StaticVector<int>* voxelsIds, VoxelsGrid* ls);
 
-    void createGraphCut(Point3d hexah[8], const StaticVector<int>& cams, VoxelsGrid* ls, const std::string& folderName, const std::string& tmpCamsPtsFolderName,
+    void createGraphCut(Point3d hexah[8], const StaticVector<int>& cams, const std::string& folderName, const std::string& tmpCamsPtsFolderName,
                         bool removeSmallSegments, const Point3d& spaceSteps);
 
     /**
diff --git a/src/aliceVision/fuseCut/ReconstructionPlan.cpp b/src/aliceVision/fuseCut/ReconstructionPlan.cpp
index cb8bf8d72f..9dc7cd4f56 100644
--- a/src/aliceVision/fuseCut/ReconstructionPlan.cpp
+++ b/src/aliceVision/fuseCut/ReconstructionPlan.cpp
@@ -248,73 +248,6 @@ void ReconstructionPlan::getHexahedronForID(float dist, int id, Point3d* out)
     mvsUtils::inflateHexahedron(&(*voxels)[id * 8], out, dist);
 }
 
-void reconstructSpaceAccordingToVoxelsArray(const std::string& voxelsArrayFileName, LargeScale* ls)
-{
-    StaticVector<Point3d>* voxelsArray = loadArrayFromFile<Point3d>(voxelsArrayFileName);
-
-    ReconstructionPlan* rp =
-        new ReconstructionPlan(ls->dimensions, &ls->space[0], ls->mp, ls->spaceVoxelsFolderName);
-
-    StaticVector<Point3d>* hexahsToExcludeFromResultingMesh = new StaticVector<Point3d>();
-    hexahsToExcludeFromResultingMesh->reserve(voxelsArray->size());
-
-    for(int i = 0; i < voxelsArray->size() / 8; i++)
-    {
-        ALICEVISION_LOG_INFO("Reconstructing " << (voxelsArray->size() / 8) << "-th Voxel of " <<  i << ".");
-
-        const std::string folderName = ls->getReconstructionVoxelFolder(i);
-        bfs::create_directory(folderName);
-
-        const std::string meshBinFilepath = folderName + "mesh.bin";
-        if(!mvsUtils::FileExists(meshBinFilepath))
-        {
-            StaticVector<int>* voxelsIds = rp->voxelsIdsIntersectingHexah(&(*voxelsArray)[i * 8]);
-            DelaunayGraphCut delaunayGC(ls->mp);
-            Point3d* hexah = &(*voxelsArray)[i * 8];
-
-            StaticVector<int> cams = ls->mp->findCamsWhichIntersectsHexahedron(hexah);
-
-            if(cams.empty())
-                throw std::logic_error("No camera to make the reconstruction");
-
-            delaunayGC.createDensePointCloudFromPrecomputedDensePoints(hexah, cams, voxelsIds, (VoxelsGrid*)rp);
-            delaunayGC.createGraphCut(hexah, cams,(VoxelsGrid*)rp, folderName, ls->getSpaceCamsTracksDir(), false, ls->getSpaceSteps());
-            delete voxelsIds;
-
-            // Save mesh as .bin and .obj
-            mesh::Mesh* mesh = delaunayGC.createMesh();
-            StaticVector<StaticVector<int>> ptsCams;
-            delaunayGC.createPtsCams(ptsCams);
-
-            mesh::meshPostProcessing(mesh, ptsCams, *ls->mp, folderName, hexahsToExcludeFromResultingMesh, hexah);
-            mesh->saveToBin(folderName + "mesh.bin");
-            mesh->saveToObj(folderName + "mesh.obj");
-
-            saveArrayOfArraysToFile<int>(folderName + "meshPtsCamsFromDGC.bin", ptsCams);
-
-            delete mesh;
-        }
-
-        /*
-        if(doComputeColoredMeshes)
-        {
-            std::string resultFolderName = folderName + "/";
-            computeColoredMesh(resultFolderName, ls);
-        }
-        */
-        Point3d hexahThin[8];
-        mvsUtils::inflateHexahedron(&(*voxelsArray)[i * 8], hexahThin, 0.9);
-        for(int k = 0; k < 8; k++)
-        {
-            hexahsToExcludeFromResultingMesh->push_back(hexahThin[k]);
-        }
-        delete hexahsToExcludeFromResultingMesh;
-    }
-    delete rp;
-    delete voxelsArray;
-}
-
-
 StaticVector<StaticVector<int>*>* loadLargeScalePtsCams(const std::vector<std::string>& recsDirs)
 {
     StaticVector<StaticVector<int>*>* ptsCamsFromDct = new StaticVector<StaticVector<int>*>();
diff --git a/src/aliceVision/fuseCut/ReconstructionPlan.hpp b/src/aliceVision/fuseCut/ReconstructionPlan.hpp
index c50ca132bb..23c0d80880 100644
--- a/src/aliceVision/fuseCut/ReconstructionPlan.hpp
+++ b/src/aliceVision/fuseCut/ReconstructionPlan.hpp
@@ -33,7 +33,6 @@ class ReconstructionPlan : public VoxelsGrid
 };
 
 void reconstructAccordingToOptimalReconstructionPlan(int gl, LargeScale* ls);
-void reconstructSpaceAccordingToVoxelsArray(const std::string& voxelsArrayFileName, LargeScale* ls);
 mesh::Mesh* joinMeshes(const std::vector<std::string>& recsDirs, StaticVector<Point3d>* voxelsArray, LargeScale* ls);
 mesh::Mesh* joinMeshes(int gl, LargeScale* ls);
 mesh::Mesh* joinMeshes(const std::string& voxelsArrayFileName, LargeScale* ls);
diff --git a/src/software/pipeline/main_meshing.cpp b/src/software/pipeline/main_meshing.cpp
index 8f19bfb0c7..a49ad978fe 100644
--- a/src/software/pipeline/main_meshing.cpp
+++ b/src/software/pipeline/main_meshing.cpp
@@ -63,15 +63,12 @@ enum ERepartitionMode
 {
     eRepartitionUndefined = 0,
     eRepartitionMultiResolution = 1,
-    eRepartitionRegularGrid = 2,
 };
 
 ERepartitionMode ERepartitionMode_stringToEnum(const std::string& s)
 {
     if(s == "multiResolution")
         return eRepartitionMultiResolution;
-    if(s == "regularGrid")
-        return eRepartitionRegularGrid;
     return eRepartitionUndefined;
 }
 
@@ -307,103 +304,6 @@ int aliceVision_main(int argc, char* argv[])
 
     switch(repartitionMode)
     {
-        case eRepartitionRegularGrid:
-        {
-            switch(partitioningMode)
-            {
-                case ePartitioningAuto:
-                {
-                    ALICEVISION_LOG_INFO("Meshing mode: regular Grid, partitioning: auto.");
-                    fuseCut::LargeScale lsbase(&mp, tmpDirectory.string() + "/");
-                    lsbase.generateSpace(maxPtsPerVoxel, ocTreeDim, true);
-                    std::string voxelsArrayFileName = lsbase.spaceFolderName + "hexahsToReconstruct.bin";
-                    StaticVector<Point3d>* voxelsArray = nullptr;
-                    if(fs::exists(voxelsArrayFileName))
-                    {
-                        // If already computed reload it.
-                        ALICEVISION_LOG_INFO("Voxels array already computed, reload from file: " << voxelsArrayFileName);
-                        voxelsArray = loadArrayFromFile<Point3d>(voxelsArrayFileName);
-                    }
-                    else
-                    {
-                        ALICEVISION_LOG_INFO("Compute voxels array.");
-                        fuseCut::ReconstructionPlan rp(lsbase.dimensions, &lsbase.space[0], lsbase.mp, lsbase.spaceVoxelsFolderName);
-                        voxelsArray = rp.computeReconstructionPlanBinSearch(fuseParams.maxPoints);
-                        saveArrayToFile<Point3d>(voxelsArrayFileName, voxelsArray);
-                    }
-                    fuseCut::reconstructSpaceAccordingToVoxelsArray(voxelsArrayFileName, &lsbase);
-                    // Join meshes and ptsCams
-                    mesh = fuseCut::joinMeshes(voxelsArrayFileName, &lsbase);
-                    fuseCut::loadLargeScalePtsCams(lsbase.getRecsDirs(voxelsArray), ptsCams);
-                    break;
-                }
-                case ePartitioningSingleBlock:
-                {
-                    ALICEVISION_LOG_INFO("Meshing mode: regular Grid, partitioning: single block.");
-                    fuseCut::LargeScale ls0(&mp, tmpDirectory.string() + "/");
-                    ls0.generateSpace(maxPtsPerVoxel, ocTreeDim, true);
-                    unsigned long ntracks = std::numeric_limits<unsigned long>::max();
-                    while(ntracks > fuseParams.maxPoints)
-                    {
-                        fs::path dirName = outDirectory/("LargeScaleMaxPts" + mvsUtils::num2strFourDecimal(ocTreeDim));
-                        fuseCut::LargeScale* ls = ls0.cloneSpaceIfDoesNotExists(ocTreeDim, dirName.string() + "/");
-                        fuseCut::VoxelsGrid vg(ls->dimensions, &ls->space[0], ls->mp, ls->spaceVoxelsFolderName);
-                        ntracks = vg.getNTracks();
-                        delete ls;
-                        ALICEVISION_LOG_INFO("Number of track candidates: " << ntracks);
-                        if(ntracks > fuseParams.maxPoints)
-                        {
-                            ALICEVISION_LOG_INFO("ocTreeDim: " << ocTreeDim);
-                            double t = (double)ntracks / (double)fuseParams.maxPoints;
-                            ALICEVISION_LOG_INFO("downsample: " << ((t < 2.0) ? "slow" : "fast"));
-                            ocTreeDim = (t < 2.0) ? ocTreeDim-100 : ocTreeDim*0.5;
-                        }
-                    }
-                    ALICEVISION_LOG_INFO("Number of tracks: " << ntracks);
-                    ALICEVISION_LOG_INFO("ocTreeDim: " << ocTreeDim);
-                    fs::path dirName = outDirectory/("LargeScaleMaxPts" + mvsUtils::num2strFourDecimal(ocTreeDim));
-                    fuseCut::LargeScale lsbase(&mp, dirName.string()+"/");
-                    lsbase.loadSpaceFromFile();
-                    fuseCut::ReconstructionPlan rp(lsbase.dimensions, &lsbase.space[0], lsbase.mp, lsbase.spaceVoxelsFolderName);
-
-                    StaticVector<int> voxelNeighs;
-                    voxelNeighs.resize(rp.voxels->size() / 8);
-                    ALICEVISION_LOG_INFO("voxelNeighs.size(): " << voxelNeighs.size());
-                    for(int i = 0; i < voxelNeighs.size(); ++i)
-                        voxelNeighs[i] = i;
-
-                    Point3d* hexah = &lsbase.space[0];
-
-                    StaticVector<int> cams;
-                    if(hexah)
-                    {
-                      cams = mp.findCamsWhichIntersectsHexahedron(hexah);
-                    }
-                    else
-                    {
-                      cams.resize(mp.getNbCameras());
-                      for(int i = 0; i < cams.size(); ++i)
-                          cams[i] = i;
-                    }
-
-                    if(cams.size() < 1)
-                        throw std::logic_error("No camera to make the reconstruction");
-
-                    fuseCut::DelaunayGraphCut delaunayGC(&mp);
-                    delaunayGC.createDensePointCloudFromPrecomputedDensePoints(hexah, cams, &voxelNeighs, (fuseCut::VoxelsGrid*)&rp);
-                    delaunayGC.createGraphCut(hexah, cams, (fuseCut::VoxelsGrid*)&rp, outDirectory.string()+"/", lsbase.getSpaceCamsTracksDir(), false, lsbase.getSpaceSteps());
-                    delaunayGC.graphCutPostProcessing();
-                    mesh = delaunayGC.createMesh();
-                    delaunayGC.createPtsCams(ptsCams);
-                    mesh::meshPostProcessing(mesh, ptsCams, mp, outDirectory.string()+"/", nullptr, hexah);
-                    break;
-                }
-                case ePartitioningUndefined:
-                default:
-                    throw std::invalid_argument("Partitioning mode is not defined");
-            }
-            break;
-        }
         case eRepartitionMultiResolution:
         {
             switch(partitioningMode)
@@ -476,7 +376,7 @@ int aliceVision_main(int argc, char* argv[])
                       sfmDataIO::Save(densePointCloud, (outDirectory/"densePointCloud_raw.abc").string(), sfmDataIO::ESfMData::ALL_DENSE);
                     }
 
-                    delaunayGC.createGraphCut(&hexah[0], cams, nullptr, outDirectory.string()+"/", outDirectory.string()+"/SpaceCamsTracks/", false, spaceSteps);
+                    delaunayGC.createGraphCut(&hexah[0], cams, outDirectory.string()+"/", outDirectory.string()+"/SpaceCamsTracks/", false, spaceSteps);
                     delaunayGC.graphCutPostProcessing();
                     mesh = delaunayGC.createMesh();
                     delaunayGC.createPtsCams(ptsCams);

From 35da7a072df6376f6fb0eb6d605cbc36e6135cea Mon Sep 17 00:00:00 2001
From: Enguerrand DE SMET <desmet.enguerrand@gmail.com>
Date: Mon, 13 Jul 2020 14:39:16 +0200
Subject: [PATCH 18/28] [fuseCut] DelaunayGrapyCut: renaming emptiness and
 fullness score

---
 src/aliceVision/fuseCut/DelaunayGraphCut.cpp  | 30 +++++++++----------
 .../fuseCut/delaunayGraphCutTypes.hpp         | 12 ++++----
 2 files changed, 21 insertions(+), 21 deletions(-)

diff --git a/src/aliceVision/fuseCut/DelaunayGraphCut.cpp b/src/aliceVision/fuseCut/DelaunayGraphCut.cpp
index bac47b0870..3137f86398 100644
--- a/src/aliceVision/fuseCut/DelaunayGraphCut.cpp
+++ b/src/aliceVision/fuseCut/DelaunayGraphCut.cpp
@@ -489,8 +489,8 @@ void DelaunayGraphCut::initCells()
         c.cellSWeight = 0.0f;
         c.cellTWeight = 0.0f;
         c.on = 0.0f;
-        c.in = 0.0f;
-        c.out = 0.0f;
+        c.fullnessScore = 0.0f;
+        c.emptinessScore = 0.0f;
         for(int s = 0; s < 4; ++s)
         {
             c.gEdgeVisWeight[s] = 0.0f; // weights for the 4 faces of the tetrahedron
@@ -1523,8 +1523,8 @@ void DelaunayGraphCut::fillGraph(bool fixesSigma, float nPixelSizeBehind, bool a
     {
         c.cellSWeight = 0.0f;
         c.cellTWeight = 0.0f;
-        c.in = 0.0f;
-        c.out = 0.0f;
+        c.fullnessScore = 0.0f;
+        c.emptinessScore = 0.0f;
         c.on = 0.0f;
         for(int s = 0; s < 4; s++)
         {
@@ -1622,7 +1622,7 @@ void DelaunayGraphCut::fillGraphPartPtRc(int& out_nstepsFront, int& out_nstepsBe
         {
             {
 #pragma OMP_ATOMIC_UPDATE
-                _cellsAttr[facet.cellIndex].out += weight;
+                _cellsAttr[facet.cellIndex].emptinessScore += weight;
             }
 
             ++out_nstepsFront;
@@ -1691,7 +1691,7 @@ void DelaunayGraphCut::fillGraphPartPtRc(int& out_nstepsFront, int& out_nstepsBe
                     c.cellTWeight += weight;
                 }
 #pragma OMP_ATOMIC_UPDATE
-                c.in += weight;
+                c.fullnessScore += weight;
             }
 
             ++out_nstepsBehind;
@@ -1787,7 +1787,7 @@ void DelaunayGraphCut::forceTedgesByGradientCVPR11(bool fixesSigma, float nPixel
 
             if((fFirst.cellIndex != GEO::NO_CELL) && (facet.cellIndex != GEO::NO_CELL) && (!isInfiniteCell(facet.cellIndex)))
             {
-                float eFirst = _cellsAttr[fFirst.cellIndex].out;
+                float eFirst = _cellsAttr[fFirst.cellIndex].emptinessScore;
 
                 Point3d p = originPt; // HAS TO BE HERE !!!
                 float maxDist = nPixelSizeBehind * mp->getCamPixelSize(p, cam);
@@ -1824,7 +1824,7 @@ void DelaunayGraphCut::forceTedgesByGradientCVPR11(bool fixesSigma, float nPixel
 
                 if(facet.cellIndex != GEO::NO_CELL)
                 {
-                    const float eLast = _cellsAttr[facet.cellIndex].out;
+                    const float eLast = _cellsAttr[facet.cellIndex].emptinessScore;
                     if((eFirst > eLast) && (eFirst < beta) && (eLast / eFirst < delta))
                     {
 #pragma OMP_ATOMIC_UPDATE
@@ -1930,13 +1930,13 @@ void DelaunayGraphCut::forceTedgesByGradientIJCV(bool fixesSigma, float nPixelSi
                     const GC_cellInfo& c = _cellsAttr[facet.cellIndex];
                     if((p - originPt).size() > nsigmaFrontSilentPart * maxDist) // (p-originPt).size() > 2 * sigma
                     {
-                        minJump = std::min(minJump, c.out);
-                        maxJump = std::max(maxJump, c.out);
+                        minJump = std::min(minJump, c.emptinessScore);
+                        maxJump = std::max(maxJump, c.emptinessScore);
                     }
                     else
                     {
-                        minSilent = std::min(minSilent, c.out);
-                        maxSilent = std::max(maxSilent, c.out);
+                        minSilent = std::min(minSilent, c.emptinessScore);
+                        maxSilent = std::max(maxSilent, c.emptinessScore);
                     }
 
                     Facet outFacet;
@@ -1969,7 +1969,7 @@ void DelaunayGraphCut::forceTedgesByGradientIJCV(bool fixesSigma, float nPixelSi
                 bool ok = (facet.cellIndex != GEO::NO_CELL);
                 if(ok)
                 {
-                    midSilent = _cellsAttr[facet.cellIndex].out;
+                    midSilent = _cellsAttr[facet.cellIndex].emptinessScore;
                 }
 
                 while(ok)
@@ -1977,8 +1977,8 @@ void DelaunayGraphCut::forceTedgesByGradientIJCV(bool fixesSigma, float nPixelSi
                     nstepsBehind++;
                     const GC_cellInfo& c = _cellsAttr[facet.cellIndex];
 
-                    minSilent = std::min(minSilent, c.out);
-                    maxSilent = std::max(maxSilent, c.out);
+                    minSilent = std::min(minSilent, c.emptinessScore);
+                    maxSilent = std::max(maxSilent, c.emptinessScore);
 
                     Facet outFacet;
                     Point3d intersectPt;
diff --git a/src/aliceVision/fuseCut/delaunayGraphCutTypes.hpp b/src/aliceVision/fuseCut/delaunayGraphCutTypes.hpp
index be0e50dbf4..1a39b5ebe2 100644
--- a/src/aliceVision/fuseCut/delaunayGraphCutTypes.hpp
+++ b/src/aliceVision/fuseCut/delaunayGraphCutTypes.hpp
@@ -24,9 +24,9 @@ struct GC_cellInfo
     /// score for emptiness along each egde/facet
     std::array<float, 4> gEdgeVisWeight{{0.0f, 0.0f, 0.0f, 0.0f}};
     /// fullness score: sum of all weights for fullness (just after the point p)
-    float in = 0.0f;
+    float fullnessScore = 0.0f;
     /// emptiness score: sum of all weights for emptiness (before the point p)
-    float out = 0.0f;
+    float emptinessScore = 0.0f;
     /// first full tetrahedron score: sum of weights for T1 (tetrahedron just after the point p)
     float on = 0.0f;
 
@@ -34,8 +34,8 @@ struct GC_cellInfo
     {
         fwrite(&cellSWeight, sizeof(float), 1, f);
         fwrite(&cellTWeight, sizeof(float), 1, f);
-        fwrite(&in, sizeof(float), 1, f);
-        fwrite(&out, sizeof(float), 1, f);
+        fwrite(&fullnessScore, sizeof(float), 1, f);
+        fwrite(&emptinessScore, sizeof(float), 1, f);
         fwrite(&on, sizeof(float), 1, f);
 
         // fwrite(gEdgePhotoWeight,sizeof(float),4,f);
@@ -46,8 +46,8 @@ struct GC_cellInfo
     {
         fread(&cellSWeight, sizeof(float), 1, f);
         fread(&cellTWeight, sizeof(float), 1, f);
-        fread(&in, sizeof(float), 1, f);
-        fread(&out, sizeof(float), 1, f);
+        fread(&fullnessScore, sizeof(float), 1, f);
+        fread(&emptinessScore, sizeof(float), 1, f);
         fread(&on, sizeof(float), 1, f);
 
         // fread(gEdgePhotoWeight,sizeof(float),4,f);

From b8c94500ee43d6564028f87d0ce15bb78b9be756 Mon Sep 17 00:00:00 2001
From: Enguerrand DE SMET <desmet.enguerrand@gmail.com>
Date: Mon, 13 Jul 2020 16:39:51 +0200
Subject: [PATCH 19/28] [fuseCut] DelaunayGraphCut: remove unnecessary behind
 option

---
 src/aliceVision/fuseCut/DelaunayGraphCut.cpp | 25 +++++++-------------
 src/aliceVision/fuseCut/DelaunayGraphCut.hpp |  4 ++--
 2 files changed, 11 insertions(+), 18 deletions(-)

diff --git a/src/aliceVision/fuseCut/DelaunayGraphCut.cpp b/src/aliceVision/fuseCut/DelaunayGraphCut.cpp
index 3137f86398..b5f9c2fcd7 100644
--- a/src/aliceVision/fuseCut/DelaunayGraphCut.cpp
+++ b/src/aliceVision/fuseCut/DelaunayGraphCut.cpp
@@ -1510,7 +1510,7 @@ float DelaunayGraphCut::weightFcn(float nrc, bool labatutWeights, int  /*ncams*/
     return weight;
 }
 
-void DelaunayGraphCut::fillGraph(bool fixesSigma, float nPixelSizeBehind, bool allPoints, bool behind,
+void DelaunayGraphCut::fillGraph(bool fixesSigma, float nPixelSizeBehind, bool allPoints,
                                bool labatutWeights, bool fillOut, float distFcnHeight) // fixesSigma=true nPixelSizeBehind=2*spaceSteps allPoints=1 behind=0 labatutWeights=0 fillOut=1 distFcnHeight=0
 {
     ALICEVISION_LOG_INFO("Computing s-t graph weights.");
@@ -1561,7 +1561,7 @@ void DelaunayGraphCut::fillGraph(bool fixesSigma, float nPixelSizeBehind, bool a
                 int nstepsFront = 0;
                 int nstepsBehind = 0;
                 fillGraphPartPtRc(nstepsFront, nstepsBehind, iV, v.cams[c], weight, fixesSigma, nPixelSizeBehind,
-                                  allPoints, behind, fillOut, distFcnHeight);
+                                  allPoints, fillOut, distFcnHeight);
 
                 avStepsFront += nstepsFront;
                 aAvStepsFront += 1;
@@ -1586,7 +1586,7 @@ void DelaunayGraphCut::fillGraph(bool fixesSigma, float nPixelSizeBehind, bool a
 
 void DelaunayGraphCut::fillGraphPartPtRc(int& out_nstepsFront, int& out_nstepsBehind, int vertexIndex, int cam,
                                        float weight, bool fixesSigma, float nPixelSizeBehind, bool allPoints,
-                                       bool behind, bool fillOut, float distFcnHeight)  // fixesSigma=true nPixelSizeBehind=2*spaceSteps allPoints=1 behind=0 fillOut=1 distFcnHeight=0
+                                       bool fillOut, float distFcnHeight)  // fixesSigma=true nPixelSizeBehind=2*spaceSteps allPoints=1 behind=0 fillOut=1 distFcnHeight=0
 {
     out_nstepsFront = 0;
     out_nstepsBehind = 0;
@@ -1685,11 +1685,6 @@ void DelaunayGraphCut::fillGraphPartPtRc(int& out_nstepsFront, int& out_nstepsBe
         {
             GC_cellInfo& c = _cellsAttr[facet.cellIndex];
             {
-                if(behind)
-                {
-#pragma OMP_ATOMIC_UPDATE
-                    c.cellTWeight += weight;
-                }
 #pragma OMP_ATOMIC_UPDATE
                 c.fullnessScore += weight;
             }
@@ -1734,13 +1729,11 @@ void DelaunayGraphCut::fillGraphPartPtRc(int& out_nstepsFront, int& out_nstepsBe
         }
 
         // cv: is the tetrahedron in distance 2*sigma behind the point p in the direction of the camera c (called Lcp in the paper)
-        if(!behind)
+
+        if(facet.cellIndex != GEO::NO_CELL)
         {
-            if(facet.cellIndex != GEO::NO_CELL)
-            {
 #pragma OMP_ATOMIC_UPDATE
-                _cellsAttr[facet.cellIndex].cellTWeight += weight;
-            }
+            _cellsAttr[facet.cellIndex].cellTWeight += weight;
         }
     }
 }
@@ -2467,7 +2460,7 @@ void DelaunayGraphCut::reconstructExpetiments(const StaticVector<int>& cams, con
 
         ALICEVISION_LOG_INFO("Jancosek CVPR 2011 method ( delta*100 = " << static_cast<int>(delta * 100.0f) << "):");
 
-        fillGraph(false, sigma, true, false, false, true, distFcnHeight);
+        fillGraph(false, sigma, true, false, true, distFcnHeight);
 
         addToInfiniteSw((float)maxint);
 
@@ -2495,7 +2488,7 @@ void DelaunayGraphCut::reconstructExpetiments(const StaticVector<int>& cams, con
         ALICEVISION_LOG_INFO("Jancosek IJCV method ( delta*100 = " << static_cast<int>(delta * 100.0f) << " ): ");
 
         // compute weights on edge between tetrahedra
-        fillGraph(false, sigma, true, false, false, true, distFcnHeight);
+        fillGraph(false, sigma, true, false, true, distFcnHeight);
 
         addToInfiniteSw((float)maxint);
 
@@ -2519,7 +2512,7 @@ void DelaunayGraphCut::reconstructExpetiments(const StaticVector<int>& cams, con
     if(labatutCFG09)
     {
         ALICEVISION_LOG_INFO("Labatut CFG 2009 method:");
-        fillGraph(false, sigma, true, false, true, true, distFcnHeight);
+        fillGraph(false, sigma, true, true, true, distFcnHeight);
 
         if(saveTemporaryBinFiles)
             saveDhInfo(folderName + "delaunayTriangulationInfoInit.bin");
diff --git a/src/aliceVision/fuseCut/DelaunayGraphCut.hpp b/src/aliceVision/fuseCut/DelaunayGraphCut.hpp
index 23a97ad08e..a3428ce4f8 100644
--- a/src/aliceVision/fuseCut/DelaunayGraphCut.hpp
+++ b/src/aliceVision/fuseCut/DelaunayGraphCut.hpp
@@ -310,10 +310,10 @@ class DelaunayGraphCut
 
     float weightFcn(float nrc, bool labatutWeights, int ncams);
 
-    virtual void fillGraph(bool fixesSigma, float nPixelSizeBehind, bool allPoints, bool behind, bool labatutWeights,
+    virtual void fillGraph(bool fixesSigma, float nPixelSizeBehind, bool allPoints, bool labatutWeights,
                            bool fillOut, float distFcnHeight = 0.0f);
     void fillGraphPartPtRc(int& out_nstepsFront, int& out_nstepsBehind, int vertexIndex, int cam, float weight,
-                           bool fixesSigma, float nPixelSizeBehind, bool allPoints, bool behind, bool fillOut,
+                           bool fixesSigma, float nPixelSizeBehind, bool allPoint, bool fillOut,
                            float distFcnHeight);
 
     void forceTedgesByGradientCVPR11(bool fixesSigma, float nPixelSizeBehind);

From 8ff0611a5baec36431e2555183db16a296628991 Mon Sep 17 00:00:00 2001
From: Enguerrand DE SMET <desmet.enguerrand@gmail.com>
Date: Mon, 13 Jul 2020 17:18:14 +0200
Subject: [PATCH 20/28] [fuseCut] DelaunayGraphCut: expose and rename
 forceTEdgeDelta

---
 src/aliceVision/fuseCut/DelaunayGraphCut.cpp | 22 ++++++++++----------
 src/software/pipeline/main_meshing.cpp       |  6 +++++-
 2 files changed, 16 insertions(+), 12 deletions(-)

diff --git a/src/aliceVision/fuseCut/DelaunayGraphCut.cpp b/src/aliceVision/fuseCut/DelaunayGraphCut.cpp
index b5f9c2fcd7..f38d4eea09 100644
--- a/src/aliceVision/fuseCut/DelaunayGraphCut.cpp
+++ b/src/aliceVision/fuseCut/DelaunayGraphCut.cpp
@@ -1743,8 +1743,8 @@ void DelaunayGraphCut::forceTedgesByGradientCVPR11(bool fixesSigma, float nPixel
     ALICEVISION_LOG_INFO("Forcing t-edges.");
     long t2 = clock();
 
-    float delta = (float)mp->userParams.get<double>("delaunaycut.delta", 0.1f);
-    ALICEVISION_LOG_INFO("delta: " << delta);
+    float forceTEdgeDelta = (float)mp->userParams.get<double>("delaunaycut.forceTEdgeDelta", 0.1f);
+    ALICEVISION_LOG_INFO("forceTEdgeDelta: " << forceTEdgeDelta);
 
     float beta = (float)mp->userParams.get<double>("delaunaycut.beta", 1000.0f);
     ALICEVISION_LOG_INFO("beta: " << beta);
@@ -1818,7 +1818,7 @@ void DelaunayGraphCut::forceTedgesByGradientCVPR11(bool fixesSigma, float nPixel
                 if(facet.cellIndex != GEO::NO_CELL)
                 {
                     const float eLast = _cellsAttr[facet.cellIndex].emptinessScore;
-                    if((eFirst > eLast) && (eFirst < beta) && (eLast / eFirst < delta))
+                    if((eFirst > eLast) && (eFirst < beta) && (eLast / eFirst < forceTEdgeDelta))
                     {
 #pragma OMP_ATOMIC_UPDATE
                         _cellsAttr[facet.cellIndex].on += (eFirst - eLast);
@@ -1845,8 +1845,8 @@ void DelaunayGraphCut::forceTedgesByGradientIJCV(bool fixesSigma, float nPixelSi
     ALICEVISION_LOG_INFO("Forcing t-edges");
     long t2 = clock();
 
-    float delta = (float)mp->userParams.get<double>("delaunaycut.delta", 0.1f);
-    ALICEVISION_LOG_DEBUG("delta: " << delta);
+    float forceTEdgeDelta = (float)mp->userParams.get<double>("delaunaycut.forceTEdgeDelta", 0.1f);
+    ALICEVISION_LOG_DEBUG("forceTEdgeDelta: " << forceTEdgeDelta);
 
     float minJumpPartRange = (float)mp->userParams.get<double>("delaunaycut.minJumpPartRange", 10000.0f);
     ALICEVISION_LOG_DEBUG("minJumpPartRange: " << minJumpPartRange);
@@ -2012,7 +2012,7 @@ void DelaunayGraphCut::forceTedgesByGradientIJCV(bool fixesSigma, float nPixelSi
                     // maxSilent: max score of emptiness for the tetrahedron around the point p (+/- 2*sigma around p)
 
                     if(
-                       (midSilent / maxJump < delta) && // (g / B) < k_rel              //// k_rel=0.1
+                       (midSilent / maxJump < forceTEdgeDelta) && // (g / B) < k_rel              //// k_rel=0.1
                        (maxJump - midSilent > minJumpPartRange) && // (B - g) > k_abs   //// k_abs=10000 // 1000 in the paper
                        (maxSilent < maxSilentPartRange)) // g < k_outl                  //// k_outl=100  // 400 in the paper
                         //(maxSilent-minSilent<maxSilentPartRange))
@@ -2456,9 +2456,9 @@ void DelaunayGraphCut::reconstructExpetiments(const StaticVector<int>& cams, con
 
     if(jancosekCVPR11)
     {
-        float delta = (float)mp->userParams.get<double>("delaunaycut.delta", 0.1f);
+        float delta = (float)mp->userParams.get<double>("delaunaycut.forceTEdgeDelta", 0.1f);
 
-        ALICEVISION_LOG_INFO("Jancosek CVPR 2011 method ( delta*100 = " << static_cast<int>(delta * 100.0f) << "):");
+        ALICEVISION_LOG_INFO("Jancosek CVPR 2011 method ( forceTEdgeDelta*100 = " << static_cast<int>(delta * 100.0f) << "):");
 
         fillGraph(false, sigma, true, false, true, distFcnHeight);
 
@@ -2483,9 +2483,9 @@ void DelaunayGraphCut::reconstructExpetiments(const StaticVector<int>& cams, con
 
     if(jancosekIJCV) // true by default
     {
-        float delta = (float)mp->userParams.get<double>("delaunaycut.delta", 0.1f);
+        float forceTEdgeDelta = (float)mp->userParams.get<double>("delaunaycut.forceTEdgeDelta", 0.1f);
 
-        ALICEVISION_LOG_INFO("Jancosek IJCV method ( delta*100 = " << static_cast<int>(delta * 100.0f) << " ): ");
+        ALICEVISION_LOG_INFO("Jancosek IJCV method ( forceTEdgeDelta*100 = " << static_cast<int>(forceTEdgeDelta * 100.0f) << " ): ");
 
         // compute weights on edge between tetrahedra
         fillGraph(false, sigma, true, false, true, distFcnHeight);
@@ -2495,7 +2495,7 @@ void DelaunayGraphCut::reconstructExpetiments(const StaticVector<int>& cams, con
         if(saveTemporaryBinFiles)
             saveDhInfo(folderName + "delaunayTriangulationInfoInit.bin");
 
-        if((delta > 0.0f) && (delta < 1.0f))
+        if((forceTEdgeDelta > 0.0f) && (forceTEdgeDelta < 1.0f))
         {
             forceTedgesByGradientIJCV(false, sigma);
         }
diff --git a/src/software/pipeline/main_meshing.cpp b/src/software/pipeline/main_meshing.cpp
index a49ad978fe..4f7c4594e0 100644
--- a/src/software/pipeline/main_meshing.cpp
+++ b/src/software/pipeline/main_meshing.cpp
@@ -147,6 +147,7 @@ int aliceVision_main(int argc, char* argv[])
     bool addLandmarksToTheDensePointCloud = false;
     bool saveRawDensePointCloud = false;
     bool colorizeOutput = false;
+    float forceTEdgeDelta = 0.1f;
 
     fuseCut::FuseParams fuseParams;
 
@@ -214,7 +215,9 @@ int aliceVision_main(int argc, char* argv[])
         ("refineFuse", po::value<bool>(&fuseParams.refineFuse)->default_value(fuseParams.refineFuse),
             "refineFuse")
         ("saveRawDensePointCloud", po::value<bool>(&saveRawDensePointCloud)->default_value(saveRawDensePointCloud),
-            "Save dense point cloud before cut and filtering.");
+            "Save dense point cloud before cut and filtering.")
+        ("forceTEdgeDelta", po::value<float>(&forceTEdgeDelta)->default_value(forceTEdgeDelta),
+            "0 to disable force T edge in graphcut. Threshold for emptiness/fullness variation.");
 
     po::options_description logParams("Log parameters");
     logParams.add_options()
@@ -286,6 +289,7 @@ int aliceVision_main(int argc, char* argv[])
     mvsUtils::MultiViewParams mp(sfmData, "", "", depthMapsFolder, meshingFromDepthMaps);
 
     mp.userParams.put("LargeScale.universePercentile", universePercentile);
+    mp.userParams.put("delaunaycut.forceTEdgeDelta", forceTEdgeDelta);
 
     int ocTreeDim = mp.userParams.get<int>("LargeScale.gridLevel0", 1024);
     const auto baseDir = mp.userParams.get<std::string>("LargeScale.baseDirName", "root01024");

From 72f5746a84765339d140b2990e3e41acea30cbd0 Mon Sep 17 00:00:00 2001
From: Enguerrand DE SMET <desmet.enguerrand@gmail.com>
Date: Wed, 15 Jul 2020 16:05:07 +0200
Subject: [PATCH 21/28] [fuseCut] Fix random generation and expose seed

- Rewrite createRandomArrayOfIntegers using stl.
- Expose seed in meshing params using MultiViewParams.
---
 src/aliceVision/fuseCut/DelaunayGraphCut.cpp | 57 ++++++++++----------
 src/aliceVision/fuseCut/Fuser.cpp            |  5 +-
 src/aliceVision/mvsUtils/common.cpp          | 56 ++++---------------
 src/aliceVision/mvsUtils/common.hpp          |  2 +-
 src/software/pipeline/main_meshing.cpp       |  6 ++-
 5 files changed, 45 insertions(+), 81 deletions(-)

diff --git a/src/aliceVision/fuseCut/DelaunayGraphCut.cpp b/src/aliceVision/fuseCut/DelaunayGraphCut.cpp
index f38d4eea09..fcef1aeaae 100644
--- a/src/aliceVision/fuseCut/DelaunayGraphCut.cpp
+++ b/src/aliceVision/fuseCut/DelaunayGraphCut.cpp
@@ -24,6 +24,7 @@
 
 #include <boost/filesystem.hpp>
 #include <boost/filesystem/operations.hpp>
+#include <random>
 
 namespace aliceVision {
 namespace fuseCut {
@@ -794,6 +795,10 @@ void DelaunayGraphCut::addHelperPoints(int nGridHelperVolumePointsDim, const Poi
 
     float maxSize = 2.0f * (O - voxel[0]).size();
     Point3d CG = (voxel[0] + voxel[1] + voxel[2] + voxel[3] + voxel[4] + voxel[5] + voxel[6] + voxel[7]) / 8.0f;
+    
+    const unsigned int seed = (unsigned int)mp->userParams.get<unsigned int>("delaunaycut.seed", 0);
+    std::mt19937 generator(seed != 0 ? seed : std::random_device{}());
+    auto rand = std::bind(std::uniform_real_distribution<float>{0.0, 1.0}, generator);
 
     for(int x = 0; x <= ns; x++)
     {
@@ -803,7 +808,7 @@ void DelaunayGraphCut::addHelperPoints(int nGridHelperVolumePointsDim, const Poi
             {
                 Point3d pt = voxel[0] + vx * ((float)x / (float)ns) + vy * ((float)y / (float)ns) +
                              vz * ((float)z / (float)ns);
-                pt = pt + (CG - pt).normalize() * (maxSize * ((float)rand() / (float)RAND_MAX));
+                pt = pt + (CG - pt).normalize() * (maxSize * rand());
 
                 Point3d p(pt.x, pt.y, pt.z);
                 GEO::index_t vi = locateNearestVertex(p);
@@ -1533,7 +1538,8 @@ void DelaunayGraphCut::fillGraph(bool fixesSigma, float nPixelSizeBehind, bool a
     }
 
     // choose random order to prevent waiting
-    StaticVector<int>* vetexesToProcessIdsRand = mvsUtils::createRandomArrayOfIntegers(_verticesAttr.size());
+    const unsigned int seed = (unsigned int)mp->userParams.get<unsigned int>("delaunaycut.seed", 0);
+    const std::vector<int> verticesRandIds = mvsUtils::createRandomArrayOfIntegers(_verticesAttr.size(), seed);
 
     int64_t avStepsFront = 0;
     int64_t aAvStepsFront = 0;
@@ -1543,10 +1549,10 @@ void DelaunayGraphCut::fillGraph(bool fixesSigma, float nPixelSizeBehind, bool a
     int nAvCams = 0;
 
 #pragma omp parallel for reduction(+:avStepsFront,aAvStepsFront,avStepsBehind,nAvStepsBehind,avCams,nAvCams)
-    for(int i = 0; i < vetexesToProcessIdsRand->size(); i++)
+    for(int i = 0; i < verticesRandIds.size(); i++)
     {
-        int iV = (*vetexesToProcessIdsRand)[i];
-        const GC_vertexInfo& v = _verticesAttr[iV];
+        int vertexIndex = verticesRandIds[i];
+        const GC_vertexInfo& v = _verticesAttr[vertexIndex];
 
         if(v.isReal() && (allPoints || v.isOnSurface) && (v.nrc > 0))
         {
@@ -1560,7 +1566,7 @@ void DelaunayGraphCut::fillGraph(bool fixesSigma, float nPixelSizeBehind, bool a
 
                 int nstepsFront = 0;
                 int nstepsBehind = 0;
-                fillGraphPartPtRc(nstepsFront, nstepsBehind, iV, v.cams[c], weight, fixesSigma, nPixelSizeBehind,
+                fillGraphPartPtRc(nstepsFront, nstepsBehind, vertexIndex, v.cams[c], weight, fixesSigma, nPixelSizeBehind,
                                   allPoints, fillOut, distFcnHeight);
 
                 avStepsFront += nstepsFront;
@@ -1574,8 +1580,6 @@ void DelaunayGraphCut::fillGraph(bool fixesSigma, float nPixelSizeBehind, bool a
         }
     }
 
-    delete vetexesToProcessIdsRand;
-
     ALICEVISION_LOG_DEBUG("avStepsFront " << avStepsFront);
     ALICEVISION_LOG_DEBUG("avStepsFront = " << mvsUtils::num2str(avStepsFront) << " // " << mvsUtils::num2str(aAvStepsFront));
     ALICEVISION_LOG_DEBUG("avStepsBehind = " << mvsUtils::num2str(avStepsBehind) << " // " << mvsUtils::num2str(nAvStepsBehind));
@@ -1755,28 +1759,29 @@ void DelaunayGraphCut::forceTedgesByGradientCVPR11(bool fixesSigma, float nPixel
     }
 
     // choose random order to prevent waiting
-    StaticVector<int>* vetexesToProcessIdsRand = mvsUtils::createRandomArrayOfIntegers(_verticesAttr.size());
+    const unsigned int seed = (unsigned int)mp->userParams.get<unsigned int>("delaunaycut.seed", 0);
+    const std::vector<int> verticesRandIds = mvsUtils::createRandomArrayOfIntegers(_verticesAttr.size(), seed);
 
 #pragma omp parallel for
-    for(int i = 0; i < vetexesToProcessIdsRand->size(); ++i)
+    for(int i = 0; i < verticesRandIds.size(); ++i)
     {
-        int vi = (*vetexesToProcessIdsRand)[i];
-        GC_vertexInfo& v = _verticesAttr[vi];
+        const int vertexIndex = verticesRandIds[i];
+        GC_vertexInfo& v = _verticesAttr[vertexIndex];
         if(v.isVirtual())
             continue;
 
-        const Point3d& originPt = _verticesCoords[vi];
+        const Point3d& originPt = _verticesCoords[vertexIndex];
         for(int c = 0; c < v.cams.size(); ++c)
         {
             int cam = v.cams[c];
 
             // False here mean nearest
             const bool nearestFarest = true;
-            Facet fFirst = getFacetFromVertexOnTheRayToTheCam(vi, cam, nearestFarest);
+            Facet fFirst = getFacetFromVertexOnTheRayToTheCam(vertexIndex, cam, nearestFarest);
 
             // get the tetrahedron next to point p on the ray from c
             // False here mean farest
-            Facet facet = getFacetFromVertexOnTheRayToTheCam(vi, cam, false);
+            Facet facet = getFacetFromVertexOnTheRayToTheCam(vertexIndex, cam, false);
 
             if((fFirst.cellIndex != GEO::NO_CELL) && (facet.cellIndex != GEO::NO_CELL) && (!isInfiniteCell(facet.cellIndex)))
             {
@@ -1830,8 +1835,6 @@ void DelaunayGraphCut::forceTedgesByGradientCVPR11(bool fixesSigma, float nPixel
 
     } // for i
 
-    delete vetexesToProcessIdsRand;
-
     for(GC_cellInfo& c: _cellsAttr)
     {
         c.cellTWeight = std::max(c.cellTWeight, std::min(1000000.0f, std::max(1.0f, c.cellTWeight) * c.on));
@@ -1872,7 +1875,8 @@ void DelaunayGraphCut::forceTedgesByGradientIJCV(bool fixesSigma, float nPixelSi
     }
 
     // choose random order to prevent waiting
-    StaticVector<int>* vetexesToProcessIdsRand = mvsUtils::createRandomArrayOfIntegers(_verticesAttr.size());
+    const unsigned int seed = (unsigned int)mp->userParams.get<unsigned int>("delaunaycut.seed", 0);
+    const std::vector<int> verticesRandIds = mvsUtils::createRandomArrayOfIntegers(_verticesAttr.size(), seed);
 
     int64_t avStepsFront = 0;
     int64_t aAvStepsFront = 0;
@@ -1880,14 +1884,14 @@ void DelaunayGraphCut::forceTedgesByGradientIJCV(bool fixesSigma, float nPixelSi
     int64_t nAvStepsBehind = 0;
 
 #pragma omp parallel for reduction(+:avStepsFront,aAvStepsFront,avStepsBehind,nAvStepsBehind)
-    for(int i = 0; i < vetexesToProcessIdsRand->size(); ++i)
+    for(int i = 0; i < verticesRandIds.size(); ++i)
     {
-        int vi = (*vetexesToProcessIdsRand)[i];
-        GC_vertexInfo& v = _verticesAttr[vi];
+        const int vertexIndex = verticesRandIds[i];
+        GC_vertexInfo& v = _verticesAttr[vertexIndex];
         if(v.isVirtual())
             continue;
 
-        const Point3d& originPt = _verticesCoords[vi];
+        const Point3d& originPt = _verticesCoords[vertexIndex];
         for(int c = 0; c < v.cams.size(); ++c)
         {
             int nstepsFront = 0;
@@ -1913,7 +1917,7 @@ void DelaunayGraphCut::forceTedgesByGradientIJCV(bool fixesSigma, float nPixelSi
             {
                 // True here mean nearest
                 const bool nearestFarest = true;
-                Facet facet = getFacetFromVertexOnTheRayToTheCam(vi, cam, nearestFarest);
+                Facet facet = getFacetFromVertexOnTheRayToTheCam(vertexIndex, cam, nearestFarest);
                 Point3d p = originPt; // HAS TO BE HERE !!!
                 bool ok = (facet.cellIndex != GEO::NO_CELL);
                 while(ok)
@@ -1957,7 +1961,7 @@ void DelaunayGraphCut::forceTedgesByGradientIJCV(bool fixesSigma, float nPixelSi
             {
                 // False here mean farest
                 const bool nearestFarest = false;
-                Facet facet = getFacetFromVertexOnTheRayToTheCam(vi, cam, nearestFarest);
+                Facet facet = getFacetFromVertexOnTheRayToTheCam(vertexIndex, cam, nearestFarest);
                 Point3d p = originPt; // HAS TO BE HERE !!!
                 bool ok = (facet.cellIndex != GEO::NO_CELL);
                 if(ok)
@@ -2030,8 +2034,6 @@ void DelaunayGraphCut::forceTedgesByGradientIJCV(bool fixesSigma, float nPixelSi
         }
     }
 
-    delete vetexesToProcessIdsRand;
-
     for(GC_cellInfo& c: _cellsAttr)
     {
         float w = std::max(1.0f, c.cellTWeight);
@@ -2283,9 +2285,6 @@ void DelaunayGraphCut::createDensePointCloud(Point3d hexah[8], const StaticVecto
   if(sfmData != nullptr)
     addPointsFromSfM(hexah, cams, *sfmData);
 
-  // initialize random seed
-  srand(time(nullptr));
-
   const int nGridHelperVolumePointsDim = mp->userParams.get<int>("LargeScale.nGridHelperVolumePointsDim", 10);
 
   // add volume points to prevent singularities
diff --git a/src/aliceVision/fuseCut/Fuser.cpp b/src/aliceVision/fuseCut/Fuser.cpp
index 6d550b3068..0fa242cb8e 100644
--- a/src/aliceVision/fuseCut/Fuser.cpp
+++ b/src/aliceVision/fuseCut/Fuser.cpp
@@ -766,7 +766,7 @@ std::string generateTempPtsSimsFiles(std::string tmpDir, mvsUtils::MultiViewPara
                 }
 
                 int nnoisePts = ((percNoisePts / 100.0f) * (float)(idsAlive->size()));
-                StaticVector<int>* randIdsAlive = mvsUtils::createRandomArrayOfIntegers(idsAlive->size());
+                const std::vector<int> randIdsAlive = mvsUtils::createRandomArrayOfIntegers(idsAlive->size());
 
                 srand(time(nullptr));
 
@@ -775,7 +775,7 @@ std::string generateTempPtsSimsFiles(std::string tmpDir, mvsUtils::MultiViewPara
                     for(int x = 0; x < w; ++x)
                     {
                         int id = y * w + x;
-                        int i = (*idsAlive)[(*randIdsAlive)[id]];
+                        int i = (*idsAlive)[randIdsAlive[id]];
                         double depth = depthMap[i];
 
                         double sim = simMap[i];
@@ -823,7 +823,6 @@ std::string generateTempPtsSimsFiles(std::string tmpDir, mvsUtils::MultiViewPara
                     mvsUtils::printfElapsedTime(t1);
 
                 delete idsAlive;
-                delete randIdsAlive;
             }
             else
             {
diff --git a/src/aliceVision/mvsUtils/common.cpp b/src/aliceVision/mvsUtils/common.cpp
index 046701ae18..88d4adfaa0 100644
--- a/src/aliceVision/mvsUtils/common.cpp
+++ b/src/aliceVision/mvsUtils/common.cpp
@@ -14,6 +14,9 @@
 #include <aliceVision/mvsData/OrientedPoint.hpp>
 #include <aliceVision/mvsData/Pixel.hpp>
 
+#include <random>
+#include <numeric>
+
 namespace aliceVision {
 namespace mvsUtils {
 
@@ -727,57 +730,16 @@ StaticVector<Point3d>* computeVoxels(const Point3d* space, const Voxel& dimensio
     return voxels;
 }
 
-StaticVector<int>* createRandomArrayOfIntegers(int n)
+std::vector<int> createRandomArrayOfIntegers(const int size, const unsigned int seed)
 {
-    /* initialize random seed: */
-    srand(time(nullptr));
-
-    StaticVector<int>* tracksPointsRandomIds = new StaticVector<int>();
-    tracksPointsRandomIds->reserve(n);
-
-    for(int j = 0; j < n; j++)
-    {
-        tracksPointsRandomIds->push_back(j);
-    }
-
-    for(int j = 0; j < n - 1; j++)
-    {
-        int rid = rand() % (n - j);
-
-        /*
-        if ((j+rid<0)||(j+rid>=tracksPoints->size())) {
-                printf("WANRING rid ot of limits %i, 0 to %i !!!! \n",j+rid,tracksPoints->size());
-        };
-        */
-
-        int v = (*tracksPointsRandomIds)[j + rid];
-        (*tracksPointsRandomIds)[j + rid] = (*tracksPointsRandomIds)[j];
-        (*tracksPointsRandomIds)[j] = v;
-    }
-
-    // test
-    /*
-    {
-            StaticVectorBool *tracksPointsRandomIdsB = new StaticVectorBool(n);
-            tracksPointsRandomIdsB->resize_with(n,false);
-            for (int k=0;k<n;k++) {
-                    int j = (*tracksPointsRandomIds)[k];
-                    (*tracksPointsRandomIdsB)[j] = true;
-            };
-
-
-            for (int j=0;j<n;j++) {
-                    if ((*tracksPointsRandomIdsB)[j]==false) {
-                            printf("WANRING  ((*tracksPointsRandomIdsB)[j]==false) !!!! \n");
-                    };
-            };
+    std::mt19937 generator(seed != 0 ? seed : std::random_device{}());
 
+    std::vector<int> v(size);
+    std::iota(v.begin(), v.end(), 0);
 
-            delete tracksPointsRandomIdsB;
-    };
-    */
+    std::shuffle(v.begin(), v.end(), generator);
 
-    return tracksPointsRandomIds;
+    return v;
 }
 
 int findNSubstrsInString(const std::string& str, const std::string& val)
diff --git a/src/aliceVision/mvsUtils/common.hpp b/src/aliceVision/mvsUtils/common.hpp
index 62e457fb82..ce5b6c77c0 100644
--- a/src/aliceVision/mvsUtils/common.hpp
+++ b/src/aliceVision/mvsUtils/common.hpp
@@ -62,7 +62,7 @@ StaticVector<StaticVector<Pixel>*>* convertObjectsCamsToCamsObjects(const MultiV
 int computeStep(MultiViewParams* mp, int scale, int maxWidth, int maxHeight);
 
 StaticVector<Point3d>* computeVoxels(const Point3d* space, const Voxel& dimensions);
-StaticVector<int>* createRandomArrayOfIntegers(int n);
+std::vector<int> createRandomArrayOfIntegers(const int size, const unsigned int seed = 0);
 
 int findNSubstrsInString(const std::string& str, const std::string& val);
 std::string num2str(int num);
diff --git a/src/software/pipeline/main_meshing.cpp b/src/software/pipeline/main_meshing.cpp
index 4f7c4594e0..1ff7138656 100644
--- a/src/software/pipeline/main_meshing.cpp
+++ b/src/software/pipeline/main_meshing.cpp
@@ -148,6 +148,7 @@ int aliceVision_main(int argc, char* argv[])
     bool saveRawDensePointCloud = false;
     bool colorizeOutput = false;
     float forceTEdgeDelta = 0.1f;
+    unsigned int seed = 0;
 
     fuseCut::FuseParams fuseParams;
 
@@ -217,7 +218,9 @@ int aliceVision_main(int argc, char* argv[])
         ("saveRawDensePointCloud", po::value<bool>(&saveRawDensePointCloud)->default_value(saveRawDensePointCloud),
             "Save dense point cloud before cut and filtering.")
         ("forceTEdgeDelta", po::value<float>(&forceTEdgeDelta)->default_value(forceTEdgeDelta),
-            "0 to disable force T edge in graphcut. Threshold for emptiness/fullness variation.");
+            "0 to disable force T edge in graphcut. Threshold for emptiness/fullness variation.")
+        ("seed", po::value<unsigned int>(&seed)->default_value(seed),
+         "Seed used in random processes. (0 to use a random seed)."); 
 
     po::options_description logParams("Log parameters");
     logParams.add_options()
@@ -290,6 +293,7 @@ int aliceVision_main(int argc, char* argv[])
 
     mp.userParams.put("LargeScale.universePercentile", universePercentile);
     mp.userParams.put("delaunaycut.forceTEdgeDelta", forceTEdgeDelta);
+    mp.userParams.put("delaunaycut.seed", seed);
 
     int ocTreeDim = mp.userParams.get<int>("LargeScale.gridLevel0", 1024);
     const auto baseDir = mp.userParams.get<std::string>("LargeScale.baseDirName", "root01024");

From aa2e28cd3a058cc79b9cd24b0087689f6246b0e5 Mon Sep 17 00:00:00 2001
From: Enguerrand DE SMET <desmet.enguerrand@gmail.com>
Date: Thu, 16 Jul 2020 15:14:27 +0200
Subject: [PATCH 22/28] [fuseCut] DelaunayGraphCut: remove unused allPoints
 variable

variable allPoints always true
---
 src/aliceVision/fuseCut/DelaunayGraphCut.cpp | 18 +++++++++---------
 src/aliceVision/fuseCut/DelaunayGraphCut.hpp |  4 ++--
 2 files changed, 11 insertions(+), 11 deletions(-)

diff --git a/src/aliceVision/fuseCut/DelaunayGraphCut.cpp b/src/aliceVision/fuseCut/DelaunayGraphCut.cpp
index fcef1aeaae..2578dd7ad6 100644
--- a/src/aliceVision/fuseCut/DelaunayGraphCut.cpp
+++ b/src/aliceVision/fuseCut/DelaunayGraphCut.cpp
@@ -1515,7 +1515,7 @@ float DelaunayGraphCut::weightFcn(float nrc, bool labatutWeights, int  /*ncams*/
     return weight;
 }
 
-void DelaunayGraphCut::fillGraph(bool fixesSigma, float nPixelSizeBehind, bool allPoints,
+void DelaunayGraphCut::fillGraph(bool fixesSigma, float nPixelSizeBehind,
                                bool labatutWeights, bool fillOut, float distFcnHeight) // fixesSigma=true nPixelSizeBehind=2*spaceSteps allPoints=1 behind=0 labatutWeights=0 fillOut=1 distFcnHeight=0
 {
     ALICEVISION_LOG_INFO("Computing s-t graph weights.");
@@ -1554,7 +1554,7 @@ void DelaunayGraphCut::fillGraph(bool fixesSigma, float nPixelSizeBehind, bool a
         int vertexIndex = verticesRandIds[i];
         const GC_vertexInfo& v = _verticesAttr[vertexIndex];
 
-        if(v.isReal() && (allPoints || v.isOnSurface) && (v.nrc > 0))
+        if(v.isReal() && (v.nrc > 0))
         {
             for(int c = 0; c < v.cams.size(); c++)
             {
@@ -1567,7 +1567,7 @@ void DelaunayGraphCut::fillGraph(bool fixesSigma, float nPixelSizeBehind, bool a
                 int nstepsFront = 0;
                 int nstepsBehind = 0;
                 fillGraphPartPtRc(nstepsFront, nstepsBehind, vertexIndex, v.cams[c], weight, fixesSigma, nPixelSizeBehind,
-                                  allPoints, fillOut, distFcnHeight);
+                                  fillOut, distFcnHeight);
 
                 avStepsFront += nstepsFront;
                 aAvStepsFront += 1;
@@ -1589,7 +1589,7 @@ void DelaunayGraphCut::fillGraph(bool fixesSigma, float nPixelSizeBehind, bool a
 }
 
 void DelaunayGraphCut::fillGraphPartPtRc(int& out_nstepsFront, int& out_nstepsBehind, int vertexIndex, int cam,
-                                       float weight, bool fixesSigma, float nPixelSizeBehind, bool allPoints,
+                                       float weight, bool fixesSigma, float nPixelSizeBehind,
                                        bool fillOut, float distFcnHeight)  // fixesSigma=true nPixelSizeBehind=2*spaceSteps allPoints=1 behind=0 fillOut=1 distFcnHeight=0
 {
     out_nstepsFront = 0;
@@ -1684,7 +1684,7 @@ void DelaunayGraphCut::fillGraphPartPtRc(int& out_nstepsFront, int& out_nstepsBe
 
         Point3d p = originPt; // HAS TO BE HERE !!!
 
-        bool ok = (facet.cellIndex != GEO::NO_CELL) && allPoints;
+        bool ok = (facet.cellIndex != GEO::NO_CELL);
         while(ok)
         {
             GC_cellInfo& c = _cellsAttr[facet.cellIndex];
@@ -1703,7 +1703,7 @@ void DelaunayGraphCut::fillGraphPartPtRc(int& out_nstepsFront, int& out_nstepsBe
             // False here mean farest
             const bool nearestFarest = false;
             if(!rayCellIntersection(mp->CArr[cam], p, facet, outFacet, nearestFarest, intersectPt) ||
-               ((originPt - p).size() >= maxDist) || (!allPoints))
+               ((originPt - p).size() >= maxDist))
             {
                 ok = false;
             }
@@ -2459,7 +2459,7 @@ void DelaunayGraphCut::reconstructExpetiments(const StaticVector<int>& cams, con
 
         ALICEVISION_LOG_INFO("Jancosek CVPR 2011 method ( forceTEdgeDelta*100 = " << static_cast<int>(delta * 100.0f) << "):");
 
-        fillGraph(false, sigma, true, false, true, distFcnHeight);
+        fillGraph(false, sigma, false, true, distFcnHeight);
 
         addToInfiniteSw((float)maxint);
 
@@ -2487,7 +2487,7 @@ void DelaunayGraphCut::reconstructExpetiments(const StaticVector<int>& cams, con
         ALICEVISION_LOG_INFO("Jancosek IJCV method ( forceTEdgeDelta*100 = " << static_cast<int>(forceTEdgeDelta * 100.0f) << " ): ");
 
         // compute weights on edge between tetrahedra
-        fillGraph(false, sigma, true, false, true, distFcnHeight);
+        fillGraph(false, sigma, false, true, distFcnHeight);
 
         addToInfiniteSw((float)maxint);
 
@@ -2511,7 +2511,7 @@ void DelaunayGraphCut::reconstructExpetiments(const StaticVector<int>& cams, con
     if(labatutCFG09)
     {
         ALICEVISION_LOG_INFO("Labatut CFG 2009 method:");
-        fillGraph(false, sigma, true, true, true, distFcnHeight);
+        fillGraph(false, sigma, true, true, distFcnHeight);
 
         if(saveTemporaryBinFiles)
             saveDhInfo(folderName + "delaunayTriangulationInfoInit.bin");
diff --git a/src/aliceVision/fuseCut/DelaunayGraphCut.hpp b/src/aliceVision/fuseCut/DelaunayGraphCut.hpp
index a3428ce4f8..68acdc5782 100644
--- a/src/aliceVision/fuseCut/DelaunayGraphCut.hpp
+++ b/src/aliceVision/fuseCut/DelaunayGraphCut.hpp
@@ -310,10 +310,10 @@ class DelaunayGraphCut
 
     float weightFcn(float nrc, bool labatutWeights, int ncams);
 
-    virtual void fillGraph(bool fixesSigma, float nPixelSizeBehind, bool allPoints, bool labatutWeights,
+    virtual void fillGraph(bool fixesSigma, float nPixelSizeBehind, bool labatutWeights,
                            bool fillOut, float distFcnHeight = 0.0f);
     void fillGraphPartPtRc(int& out_nstepsFront, int& out_nstepsBehind, int vertexIndex, int cam, float weight,
-                           bool fixesSigma, float nPixelSizeBehind, bool allPoint, bool fillOut,
+                           bool fixesSigma, float nPixelSizeBehind, bool fillOut,
                            float distFcnHeight);
 
     void forceTedgesByGradientCVPR11(bool fixesSigma, float nPixelSizeBehind);

From fc9934c7be6a533f14d9dfc066e68a33217cc2d9 Mon Sep 17 00:00:00 2001
From: Fabien Castan <fabcastan@gmail.com>
Date: Thu, 16 Jul 2020 14:50:26 +0200
Subject: [PATCH 23/28] [fuseCut] add debug method to export a mesh
 representation of the tetrahedralization vote

---
 src/aliceVision/fuseCut/DelaunayGraphCut.cpp | 254 +++++++++++++++++--
 src/aliceVision/fuseCut/DelaunayGraphCut.hpp |   8 +-
 2 files changed, 233 insertions(+), 29 deletions(-)

diff --git a/src/aliceVision/fuseCut/DelaunayGraphCut.cpp b/src/aliceVision/fuseCut/DelaunayGraphCut.cpp
index 2578dd7ad6..215d8c6b66 100644
--- a/src/aliceVision/fuseCut/DelaunayGraphCut.cpp
+++ b/src/aliceVision/fuseCut/DelaunayGraphCut.cpp
@@ -26,6 +26,9 @@
 #include <boost/filesystem/operations.hpp>
 #include <random>
 
+#include <boost/accumulators/accumulators.hpp>
+#include <boost/accumulators/statistics.hpp>
+
 namespace aliceVision {
 namespace fuseCut {
 
@@ -1556,11 +1559,11 @@ void DelaunayGraphCut::fillGraph(bool fixesSigma, float nPixelSizeBehind,
 
         if(v.isReal() && (v.nrc > 0))
         {
+            // "weight" is called alpha(p) in the paper
+            const float weight = weightFcn((float)v.nrc, labatutWeights, v.getNbCameras()); // number of cameras
+
             for(int c = 0; c < v.cams.size(); c++)
             {
-                // "weight" is called alpha(p) in the paper
-                float weight = weightFcn((float)v.nrc, labatutWeights, v.getNbCameras()); // number of cameras
-
                 assert(v.cams[c] >= 0);
                 assert(v.cams[c] < mp->ncams);
 
@@ -2304,9 +2307,20 @@ void DelaunayGraphCut::createGraphCut(Point3d hexah[8], const StaticVector<int>&
   if(removeSmallSegments) // false
     removeSmallSegs(2500); // TODO FACA: to decide
 
-  reconstructExpetiments(cams, folderName,
-                         hexah, tmpCamsPtsFolderName,
-                         spaceSteps);
+  voteFullEmptyScore(cams, folderName, tmpCamsPtsFolderName, spaceSteps);
+
+  if(saveTemporaryBinFiles)
+  {
+      std::unique_ptr<mesh::Mesh> mesh(createTetrahedralMesh());
+      mesh->saveToObj(folderName + "tetrahedralMesh.obj");
+  }
+
+  reconstructGC(hexah);
+
+  if(saveTemporaryBinFiles)
+  {
+      saveDhInfo(folderName + "delaunayTriangulationInfoAfterHallRemoving.bin");
+  }
 
   if(mp->userParams.get<bool>("LargeScale.saveDelaunayTriangulation", false))
   {
@@ -2428,11 +2442,11 @@ void DelaunayGraphCut::maxflow()
     ALICEVISION_LOG_INFO("Maxflow: done.");
 }
 
-void DelaunayGraphCut::reconstructExpetiments(const StaticVector<int>& cams, const std::string& folderName,
-                                            Point3d hexahInflated[8], const std::string& tmpCamsPtsFolderName,
-                                            const Point3d& spaceSteps)
+void DelaunayGraphCut::voteFullEmptyScore(const StaticVector<int>& cams, const std::string& folderName,
+                                          const std::string& tmpCamsPtsFolderName,
+                                          const Point3d& spaceSteps)
 {
-    ALICEVISION_LOG_INFO("DelaunayGraphCut::reconstructExpetiments");
+    ALICEVISION_LOG_INFO("DelaunayGraphCut::voteFullEmptyScore");
     int maxint = 1000000.0f;
 
     long t1;
@@ -2473,11 +2487,6 @@ void DelaunayGraphCut::reconstructExpetiments(const StaticVector<int>& cams, con
 
         if(saveTemporaryBinFiles)
             saveDhInfo(folderName + "delaunayTriangulationInfoAfterForce.bin");
-
-        reconstructGC(hexahInflated);
-
-        if(saveTemporaryBinFiles)
-            saveDhInfo(folderName + "delaunayTriangulationInfoAfterHallRemoving.bin");
     }
 
     if(jancosekIJCV) // true by default
@@ -2501,11 +2510,6 @@ void DelaunayGraphCut::reconstructExpetiments(const StaticVector<int>& cams, con
 
         if(saveTemporaryBinFiles)
             saveDhInfo(folderName + "delaunayTriangulationInfoAfterForce.bin");
-
-        reconstructGC(hexahInflated);
-
-        if(saveTemporaryBinFiles)
-            saveDhInfo(folderName + "delaunayTriangulationInfoAfterHallRemoving.bin");
     }
 
     if(labatutCFG09)
@@ -2515,11 +2519,6 @@ void DelaunayGraphCut::reconstructExpetiments(const StaticVector<int>& cams, con
 
         if(saveTemporaryBinFiles)
             saveDhInfo(folderName + "delaunayTriangulationInfoInit.bin");
-
-        reconstructGC(hexahInflated);
-
-        if(saveTemporaryBinFiles)
-            saveDhInfo(folderName + "delaunayTriangulationInfoAfterHallRemoving.bin");
     }
 }
 
@@ -2696,6 +2695,211 @@ mesh::Mesh* DelaunayGraphCut::createMesh(bool filterHelperPointsTriangles)
     return me;
 }
 
+mesh::Mesh* DelaunayGraphCut::createTetrahedralMesh() const
+{
+    ALICEVISION_LOG_INFO("Create mesh of the tetrahedralization.");
+
+    ALICEVISION_LOG_INFO("# vertixes: " << _verticesCoords.size());
+    if(_cellsAttr.empty())
+        return nullptr;
+
+    mesh::Mesh* me = new mesh::Mesh();
+
+    // TODO: copy only surface points and remap visibilities
+    me->pts = StaticVector<Point3d>();
+    me->pts.reserve(10 * _verticesCoords.size());
+
+    me->tris = StaticVector<mesh::Mesh::triangle>();
+    me->tris.reserve(_verticesCoords.size());
+
+    using namespace boost::accumulators;
+    using Accumulator = accumulator_set<float, stats<tag::min, tag::max, tag::median, tag::mean>>;
+    auto displayAcc = [](const std::string& name, const Accumulator& acc) {
+        ALICEVISION_LOG_INFO(" [" << name << "]"
+                                  << " min: " << extract::min(acc) << " max: " << extract::max(acc)
+                                  << " mean: " << extract::mean(acc) << " median: " << extract::median(acc));
+    };
+    {
+        Accumulator acc_nrc;
+        Accumulator acc_camSize;
+        for(VertexIndex vi = 0; vi < _verticesAttr.size(); ++vi)
+        {
+            const GC_vertexInfo& vertexAttr = _verticesAttr[vi];
+            acc_nrc(vertexAttr.nrc);
+            acc_camSize(vertexAttr.cams.size());
+        }
+        displayAcc("acc_nrc", acc_nrc);
+        displayAcc("acc_camSize", acc_camSize);
+    }
+
+    float max_score = 1.0;
+    {
+        Accumulator acc_cellScore;
+        Accumulator acc_cellSWeight;
+        Accumulator acc_cellTWeight;
+        Accumulator acc_gEdgeVisWeight;
+        Accumulator acc_fullnessScore;
+        Accumulator acc_emptinessScore;
+        Accumulator acc_on;
+
+        for(CellIndex ci = 0; ci < _cellsAttr.size(); ++ci)
+        {
+            const GC_cellInfo& cellAttr = _cellsAttr[ci];
+            acc_cellScore(cellAttr.cellSWeight - cellAttr.cellTWeight);
+            acc_cellSWeight(cellAttr.cellSWeight);
+            acc_cellTWeight(cellAttr.cellTWeight);
+
+            acc_gEdgeVisWeight(cellAttr.gEdgeVisWeight[0]);
+            acc_gEdgeVisWeight(cellAttr.gEdgeVisWeight[1]);
+            acc_gEdgeVisWeight(cellAttr.gEdgeVisWeight[2]);
+            acc_gEdgeVisWeight(cellAttr.gEdgeVisWeight[3]);
+
+            acc_fullnessScore(cellAttr.fullnessScore);
+            acc_emptinessScore(cellAttr.emptinessScore);
+            acc_on(cellAttr.on);
+        }
+        
+        displayAcc("cellScore", acc_cellScore);
+        displayAcc("cellSWeight", acc_cellSWeight);
+        displayAcc("cellTWeight", acc_cellTWeight);
+        displayAcc("gEdgeVisWeight", acc_gEdgeVisWeight);
+        displayAcc("fullnessScore", acc_fullnessScore);
+        displayAcc("emptinessScore", acc_emptinessScore);
+        displayAcc("on", acc_on);
+        max_score = 4.0f * extract::mean(acc_emptinessScore);
+    }
+    ALICEVISION_LOG_DEBUG("createTetrahedralMesh: max_score: " << max_score);
+
+    // loop over all facets
+    for(CellIndex ci = 0; ci < _cellsAttr.size(); ++ci)
+    {
+        Point3d pointscellCenter(0.0, 0.0, 0.0);
+        {
+            bool invalid = false;
+            int weakVertex = 0;
+            for(VertexIndex k = 0; k < 4; ++k)
+            {
+                const VertexIndex vi = _tetrahedralization->cell_vertex(ci, k);
+                const GC_vertexInfo& vertexAttr = _verticesAttr[vi];
+                if(vertexAttr.cams.size() <= 3)
+                {
+                    ++weakVertex;
+                }
+                const Facet f1(ci, k);
+                const Facet f2 = mirrorFacet(f1);
+                if(isInvalidOrInfiniteCell(f2.cellIndex))
+                {
+                    invalid = true;
+                    continue;
+                }
+                pointscellCenter = pointscellCenter + _verticesCoords[vi];
+            }
+            pointscellCenter = pointscellCenter / 4.0;
+            if(invalid)
+            {
+                // If one of the facet is invalid, discard the tetrahedron.
+                continue;
+            }
+            if (weakVertex >= 3)
+            {
+                // If the tetrahedron mainly composed of weak vertices, skip it.
+                // So keep tetrahedra that are linked to at least 3 good vertices.
+                continue;
+            }
+        }
+
+        const GC_cellInfo& cellAttr = _cellsAttr[ci];
+        const float score = cellAttr.emptinessScore; // cellAttr.cellSWeight - cellAttr.cellTWeight;
+        if(score < (max_score / 1000.0f))
+        {
+            // Skip too low score cells
+            continue;
+        }
+
+        for(VertexIndex k = 0; k < 4; ++k)
+        {
+            const Facet f1(ci, k);
+
+            VertexIndex vertices[3];
+            vertices[0] = getVertexIndex(f1, 0);
+            vertices[1] = getVertexIndex(f1, 1);
+            vertices[2] = getVertexIndex(f1, 2);
+
+            Point3d points[3];
+            for(int k = 0; k < 3; ++k)
+            {
+                points[k] = _verticesCoords[vertices[k]];
+                // Downscale cell for visibility
+                points[k] = pointscellCenter + ((points[k] - pointscellCenter) * 0.95);
+            }
+
+            const Facet f2 = mirrorFacet(f1);
+            //// do not need to test again: already filtered before
+            // if(isInvalidOrInfiniteCell(f2.cellIndex))
+            //     continue;
+
+            const Point3d D1 = _verticesCoords[getOppositeVertexIndex(f1)];
+            const Point3d D2 = _verticesCoords[getOppositeVertexIndex(f2)];
+
+            const Point3d N =
+                cross((points[1] - points[0]).normalize(), (points[2] - points[0]).normalize()).normalize();
+
+            const double dd1 = orientedPointPlaneDistance(D1, points[0], N);
+            const double dd2 = orientedPointPlaneDistance(D2, points[0], N);
+
+            bool clockwise = false;
+            if(dd1 == 0.0)
+            {
+                if(dd2 == 0.0)
+                {
+                    ALICEVISION_LOG_WARNING("createMesh: bad triangle orientation.");
+                }
+                if(dd2 > 0.0)
+                {
+                    clockwise = true;
+                }
+            }
+            else
+            {
+                if(dd1 < 0.0)
+                {
+                    clockwise = true;
+                }
+            }
+
+            const rgb color = getRGBFromJetColorMap(score / max_score);
+            const std::size_t vertexBaseIndex = me->pts.size();
+            for(const Point3d& p: points)
+            {
+                me->pts.push_back(p);
+                me->colors().push_back(color);
+            }
+
+            if(clockwise)
+            {
+                mesh::Mesh::triangle t;
+                t.alive = true;
+                t.v[0] = vertexBaseIndex;
+                t.v[1] = vertexBaseIndex + 1;
+                t.v[2] = vertexBaseIndex + 2;
+                me->tris.push_back(t);
+            }
+            else
+            {
+                mesh::Mesh::triangle t;
+                t.alive = true;
+                t.v[0] = vertexBaseIndex;
+                t.v[1] = vertexBaseIndex + 2;
+                t.v[2] = vertexBaseIndex + 1;
+                me->tris.push_back(t);
+            }
+        }
+    }
+
+    ALICEVISION_LOG_INFO("Extract mesh from Graph Cut done.");
+    return me;
+}
+
 void DelaunayGraphCut::segmentFullOrFree(bool full, StaticVector<int>** out_fullSegsColor, int& out_nsegments)
 {
     ALICEVISION_LOG_DEBUG("segmentFullOrFree: segmenting connected space.");
diff --git a/src/aliceVision/fuseCut/DelaunayGraphCut.hpp b/src/aliceVision/fuseCut/DelaunayGraphCut.hpp
index 68acdc5782..e8d9f66b10 100644
--- a/src/aliceVision/fuseCut/DelaunayGraphCut.hpp
+++ b/src/aliceVision/fuseCut/DelaunayGraphCut.hpp
@@ -329,10 +329,9 @@ class DelaunayGraphCut
 
     void maxflow();
 
-    void reconstructExpetiments(const StaticVector<int>& cams, const std::string& folderName,
-                                Point3d hexahInflated[8], const std::string& tmpCamsPtsFolderName,
-                                const Point3d& spaceSteps);
-
+    void voteFullEmptyScore(const StaticVector<int>& cams, const std::string& folderName,
+                            const std::string& tmpCamsPtsFolderName,
+                            const Point3d& spaceSteps);
 
     void createDensePointCloud(Point3d hexah[8], const StaticVector<int>& cams, const sfmData::SfMData* sfmData, const FuseParams* depthMapsFuseParams);
 
@@ -352,6 +351,7 @@ class DelaunayGraphCut
     void leaveLargestFullSegmentOnly();
 
     mesh::Mesh* createMesh(bool filterHelperPointsTriangles = true);
+    mesh::Mesh* createTetrahedralMesh() const;
 };
 
 } // namespace fuseCut

From b5ee9fe7e8deba05f1226bbfb45d5c76bb7e2ce6 Mon Sep 17 00:00:00 2001
From: Fabien Castan <fabcastan@gmail.com>
Date: Thu, 16 Jul 2020 14:51:30 +0200
Subject: [PATCH 24/28] [mvsData] fix jetColorMap on unsigned char

---
 src/aliceVision/mvsData/jetColorMap.cpp | 36 +++++++++----------------
 1 file changed, 13 insertions(+), 23 deletions(-)

diff --git a/src/aliceVision/mvsData/jetColorMap.cpp b/src/aliceVision/mvsData/jetColorMap.cpp
index 9b1a622d35..9db9c3c4a3 100644
--- a/src/aliceVision/mvsData/jetColorMap.cpp
+++ b/src/aliceVision/mvsData/jetColorMap.cpp
@@ -30,35 +30,17 @@ static float jetb[64] = {0.5625, 0.6250, 0.6875, 0.7500, 0.8125, 0.8750, 0.9375,
                          0,      0,      0,      0,      0,      0,      0,      0,      0};
 
 
-rgb getRGBFromJetColorMap(float value)
-{
-    if(value <= 0.0f)
-        return rgb(0, 0, 0);
-    if(value >= 1.0f)
-        return rgb(1, 1, 1);
-    float idx_f = value * 63.0f;
-    float fractA, fractB, integral;
-    fractB = std::modf(idx_f, &integral);
-    fractA = 1.0f - fractB;
-    int idx = static_cast<int>(integral);
-    rgb c;
-    c.r = static_cast<unsigned char>((jetr[idx] * fractA + jetr[idx + 1] * fractB) * 255.0f);
-    c.g = static_cast<unsigned char>((jetg[idx] * fractA + jetg[idx + 1] * fractB) * 255.0f);
-    c.b = static_cast<unsigned char>((jetb[idx] * fractA + jetb[idx + 1] * fractB) * 255.0f);
-    return c;
-}
-
 Color getColorFromJetColorMap(float value)
 {
     if(value <= 0.0f)
         return Color(0, 0, 0);
     if(value >= 1.0f)
         return Color(1.0f, 1.0f, 1.0f);
-    float idx_f = value * 63.0f;
-    float fractA, fractB, integral;
-    fractB = std::modf(idx_f, &integral);
-    fractA = 1.0f - fractB;
-    int idx = static_cast<int>(integral);
+    const float idx_f = value * 63.0f;
+    float integral;
+    const float fractB = std::modf(idx_f, &integral);
+    const float fractA = 1.0f - fractB;
+    const int idx = static_cast<int>(integral);
     Color c;
     c.r = jetr[idx] * fractA + jetr[idx + 1] * fractB;
     c.g = jetg[idx] * fractA + jetg[idx + 1] * fractB;
@@ -66,4 +48,12 @@ Color getColorFromJetColorMap(float value)
     return c;
 }
 
+rgb getRGBFromJetColorMap(float value)
+{
+    const Color color = getColorFromJetColorMap(value);
+    return {static_cast<unsigned char>(color.r * 255.0f),
+            static_cast<unsigned char>(color.g * 255.0f),
+            static_cast<unsigned char>(color.b * 255.0f)};
+}
+
 } // namespace aliceVision

From 891a1f0f788efdb77307654b08480ecbacafa46c Mon Sep 17 00:00:00 2001
From: Fabien Castan <fabcastan@gmail.com>
Date: Thu, 16 Jul 2020 15:30:13 +0200
Subject: [PATCH 25/28] [fuseCut] remove unless code

---
 src/aliceVision/fuseCut/DelaunayGraphCut.cpp | 25 +++++++++-----------
 src/aliceVision/fuseCut/DelaunayGraphCut.hpp |  7 ++----
 src/software/pipeline/main_meshing.cpp       | 23 +-----------------
 3 files changed, 14 insertions(+), 41 deletions(-)

diff --git a/src/aliceVision/fuseCut/DelaunayGraphCut.cpp b/src/aliceVision/fuseCut/DelaunayGraphCut.cpp
index 215d8c6b66..f47914b27a 100644
--- a/src/aliceVision/fuseCut/DelaunayGraphCut.cpp
+++ b/src/aliceVision/fuseCut/DelaunayGraphCut.cpp
@@ -1138,7 +1138,7 @@ void DelaunayGraphCut::fuseFromDepthMaps(const StaticVector<int>& cams, const Po
 
 void DelaunayGraphCut::computeVerticesSegSize(bool allPoints, float alpha) // allPoints=true, alpha=0
 {
-    ALICEVISION_LOG_DEBUG("creating universe");
+    ALICEVISION_LOG_DEBUG("DelaunayGraphCut::computeVerticesSegSize");
     int scalePS = mp->userParams.get<int>("global.scalePS", 1);
     int step = mp->userParams.get<int>("global.step", 1);
     float pointToJoinPixSizeDist = (float)mp->userParams.get<double>("delaunaycut.pointToJoinPixSizeDist", 2.0) *
@@ -1188,17 +1188,17 @@ void DelaunayGraphCut::computeVerticesSegSize(bool allPoints, float alpha) // al
     }
     // mvsUtils::finishEstimate();
 
-    Universe* u = new Universe(_verticesAttr.size());
+    Universe u(_verticesAttr.size());
 
     // t1 = mvsUtils::initEstimate();
     int s = (int)edges.size(); // Fuse all edges collected to be merged
     for(int i = 0; i < s; i++)
     {
-        int a = u->find(edges[i].x);
-        int b = u->find(edges[i].y);
+        int a = u.find(edges[i].x);
+        int b = u.find(edges[i].y);
         if(a != b) // TODO FACA: Are they not different in all cases?
         {
-            u->join(a, b);
+            u.join(a, b);
         }
         // mvsUtils::printfEstimate(i, s, t1);
     }
@@ -1211,13 +1211,12 @@ void DelaunayGraphCut::computeVerticesSegSize(bool allPoints, float alpha) // al
         if(v.isVirtual())
             continue;
 
-        int a = u->find(vi);
-        v.segSize = u->elts[a].size;
+        int a = u.find(vi);
+        v.segSize = u.elts[a].size;
         v.segId = a;
     }
 
-    delete u;
-    ALICEVISION_LOG_DEBUG("creating universe done.");
+    ALICEVISION_LOG_DEBUG("DelaunayGraphCut::computeVerticesSegSize done.");
 }
 
 void DelaunayGraphCut::removeSmallSegs(int minSegSize)
@@ -2294,7 +2293,7 @@ void DelaunayGraphCut::createDensePointCloud(Point3d hexah[8], const StaticVecto
   addHelperPoints(nGridHelperVolumePointsDim, hexah, minDist);
 }
 
-void DelaunayGraphCut::createGraphCut(Point3d hexah[8], const StaticVector<int>& cams, const std::string& folderName, const std::string& tmpCamsPtsFolderName, bool removeSmallSegments, const Point3d& spaceSteps)
+void DelaunayGraphCut::createGraphCut(Point3d hexah[8], const StaticVector<int>& cams, const std::string& folderName, const std::string& tmpCamsPtsFolderName, bool removeSmallSegments)
 {
   initVertices();
 
@@ -2307,7 +2306,7 @@ void DelaunayGraphCut::createGraphCut(Point3d hexah[8], const StaticVector<int>&
   if(removeSmallSegments) // false
     removeSmallSegs(2500); // TODO FACA: to decide
 
-  voteFullEmptyScore(cams, folderName, tmpCamsPtsFolderName, spaceSteps);
+  voteFullEmptyScore(cams, folderName);
 
   if(saveTemporaryBinFiles)
   {
@@ -2442,9 +2441,7 @@ void DelaunayGraphCut::maxflow()
     ALICEVISION_LOG_INFO("Maxflow: done.");
 }
 
-void DelaunayGraphCut::voteFullEmptyScore(const StaticVector<int>& cams, const std::string& folderName,
-                                          const std::string& tmpCamsPtsFolderName,
-                                          const Point3d& spaceSteps)
+void DelaunayGraphCut::voteFullEmptyScore(const StaticVector<int>& cams, const std::string& folderName)
 {
     ALICEVISION_LOG_INFO("DelaunayGraphCut::voteFullEmptyScore");
     int maxint = 1000000.0f;
diff --git a/src/aliceVision/fuseCut/DelaunayGraphCut.hpp b/src/aliceVision/fuseCut/DelaunayGraphCut.hpp
index e8d9f66b10..22053d6aa5 100644
--- a/src/aliceVision/fuseCut/DelaunayGraphCut.hpp
+++ b/src/aliceVision/fuseCut/DelaunayGraphCut.hpp
@@ -329,14 +329,11 @@ class DelaunayGraphCut
 
     void maxflow();
 
-    void voteFullEmptyScore(const StaticVector<int>& cams, const std::string& folderName,
-                            const std::string& tmpCamsPtsFolderName,
-                            const Point3d& spaceSteps);
+    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 createGraphCut(Point3d hexah[8], const StaticVector<int>& cams, const std::string& folderName, const std::string& tmpCamsPtsFolderName,
-                        bool removeSmallSegments, const Point3d& spaceSteps);
+    void createGraphCut(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.
diff --git a/src/software/pipeline/main_meshing.cpp b/src/software/pipeline/main_meshing.cpp
index 1ff7138656..428469c9dd 100644
--- a/src/software/pipeline/main_meshing.cpp
+++ b/src/software/pipeline/main_meshing.cpp
@@ -333,27 +333,6 @@ int aliceVision_main(int argc, char* argv[])
                     else
                       fs.divideSpaceFromSfM(sfmData, &hexah[0], estimateSpaceMinObservations, estimateSpaceMinObservationAngle);
 
-                    Voxel dimensions = fs.estimateDimensions(&hexah[0], &hexah[0], 0, ocTreeDim, (meshingFromDepthMaps && !estimateSpaceFromSfM) ? nullptr : &sfmData);
-                    StaticVector<Point3d>* voxels = mvsUtils::computeVoxels(&hexah[0], dimensions);
-
-                    StaticVector<int> voxelNeighs;
-                    voxelNeighs.resize(voxels->size() / 8);
-                    ALICEVISION_LOG_INFO("voxelNeighs.size(): " << voxelNeighs.size());
-
-                    for(int i = 0; i < voxelNeighs.size(); ++i)
-                        voxelNeighs[i] = i;
-
-                    Point3d spaceSteps;
-                    {
-                        Point3d vx = hexah[1] - hexah[0];
-                        Point3d vy = hexah[3] - hexah[0];
-                        Point3d vz = hexah[4] - hexah[0];
-                        spaceSteps.x = (vx.size() / (double)dimensions.x) / (double)ocTreeDim;
-                        spaceSteps.y = (vy.size() / (double)dimensions.y) / (double)ocTreeDim;
-                        spaceSteps.z = (vz.size() / (double)dimensions.z) / (double)ocTreeDim;
-                    }
-                    delete voxels;
-
                     StaticVector<int> cams;
                     if(meshingFromDepthMaps)
                     {
@@ -384,7 +363,7 @@ int aliceVision_main(int argc, char* argv[])
                       sfmDataIO::Save(densePointCloud, (outDirectory/"densePointCloud_raw.abc").string(), sfmDataIO::ESfMData::ALL_DENSE);
                     }
 
-                    delaunayGC.createGraphCut(&hexah[0], cams, outDirectory.string()+"/", outDirectory.string()+"/SpaceCamsTracks/", false, spaceSteps);
+                    delaunayGC.createGraphCut(&hexah[0], cams, outDirectory.string()+"/", outDirectory.string()+"/SpaceCamsTracks/", false);
                     delaunayGC.graphCutPostProcessing();
                     mesh = delaunayGC.createMesh();
                     delaunayGC.createPtsCams(ptsCams);

From 1d48150ea0a69580161512f1d289eede9a60d573 Mon Sep 17 00:00:00 2001
From: Fabien Castan <fabcastan@gmail.com>
Date: Thu, 16 Jul 2020 20:40:19 +0200
Subject: [PATCH 26/28] [mvsUtils] do not fail when meshing an SfM without the
 source images

readImageMetadata throw an error if the path does not exist.
---
 src/aliceVision/mvsUtils/MultiViewParams.cpp | 26 ++++++++++++--------
 src/aliceVision/mvsUtils/MultiViewParams.hpp |  2 +-
 2 files changed, 17 insertions(+), 11 deletions(-)

diff --git a/src/aliceVision/mvsUtils/MultiViewParams.cpp b/src/aliceVision/mvsUtils/MultiViewParams.cpp
index 8489798190..d6ae29fe56 100644
--- a/src/aliceVision/mvsUtils/MultiViewParams.cpp
+++ b/src/aliceVision/mvsUtils/MultiViewParams.cpp
@@ -117,10 +117,16 @@ MultiViewParams::MultiViewParams(const sfmData::SfMData& sfmData,
         const ImageParams& imgParams = _imagesParams.at(i);
 
         oiio::ParamValueList metadata;
-        imageIO::readImageMetadata(imgParams.path, metadata);
-
-        const auto scaleIt = metadata.find("AliceVision:downscale");
-        const auto pIt = metadata.find("AliceVision:P");
+        oiio::ParamValueList::const_iterator scaleIt = metadata.end();
+        oiio::ParamValueList::const_iterator pIt = metadata.end();
+        
+        const bool fileExists = fs::exists(imgParams.path);
+        if(fileExists)
+        {
+            imageIO::readImageMetadata(imgParams.path, metadata);
+            scaleIt = metadata.find("AliceVision:downscale");
+            pIt = metadata.find("AliceVision:P");
+        }
 
         // find image scale information
         if(scaleIt != metadata.end() && scaleIt->type() == oiio::TypeDesc::INT)
@@ -128,7 +134,7 @@ MultiViewParams::MultiViewParams(const sfmData::SfMData& sfmData,
             // use aliceVision image metadata
             _imagesScale.at(i) = scaleIt->get_int();
         }
-        else
+        else if(fileExists)
         {
             // use image dimension
             int w, h, channels;
@@ -172,16 +178,16 @@ MultiViewParams::MultiViewParams(const sfmData::SfMData& sfmData,
             const std::string fileNameP = getFileNameFromIndex(this, i, EFileType::P);
             const std::string fileNameD = getFileNameFromIndex(this, i, EFileType::D);
 
-            if(fs::exists(fileNameP), fs::exists(fileNameD))
+            if(fs::exists(fileNameP) && fs::exists(fileNameD))
             {
-              ALICEVISION_LOG_DEBUG("Reading view " << getViewId(i) << " projection matrix from file '" << fileNameP << "'.");
+                ALICEVISION_LOG_DEBUG("Reading view " << getViewId(i) << " projection matrix from file '" << fileNameP << "'.");
 
-              loadMatricesFromTxtFile(i, fileNameP, fileNameD);
+                loadMatricesFromTxtFile(i, fileNameP, fileNameD);
             }
             else
             {
-              ALICEVISION_LOG_DEBUG("Reading view " << getViewId(i) << " projection matrix from SfMData.");
-              loadMatricesFromSfM(i);
+                ALICEVISION_LOG_DEBUG("Reading view " << getViewId(i) << " projection matrix from SfMData.");
+                loadMatricesFromSfM(i);
             }
         }
 
diff --git a/src/aliceVision/mvsUtils/MultiViewParams.hpp b/src/aliceVision/mvsUtils/MultiViewParams.hpp
index 7d85bf3b22..03e18e3c85 100644
--- a/src/aliceVision/mvsUtils/MultiViewParams.hpp
+++ b/src/aliceVision/mvsUtils/MultiViewParams.hpp
@@ -333,7 +333,7 @@ class MultiViewParams
         CArr.resize(ncams);
         iCamArr.resize(ncams);
         FocK1K2Arr.resize(ncams);
-        _imagesScale.resize(ncams);
+        _imagesScale.resize(ncams, 1);
     }
 };
 

From d970a82699cbed98208c9f8c9c6ad23eded2ec78 Mon Sep 17 00:00:00 2001
From: Enguerrand DE SMET <desmet.enguerrand@gmail.com>
Date: Mon, 20 Jul 2020 10:37:09 +0200
Subject: [PATCH 27/28] [fuseCut] delaunayGraphCut: remove the unused
 JancosekCVPR11 method

---
 src/aliceVision/fuseCut/DelaunayGraphCut.cpp  | 124 ------------------
 src/aliceVision/fuseCut/DelaunayGraphCut.hpp  |   1 -
 .../fuseCut/ReconstructionPlan.cpp            |   4 -
 3 files changed, 129 deletions(-)

diff --git a/src/aliceVision/fuseCut/DelaunayGraphCut.cpp b/src/aliceVision/fuseCut/DelaunayGraphCut.cpp
index f47914b27a..0d51ba517d 100644
--- a/src/aliceVision/fuseCut/DelaunayGraphCut.cpp
+++ b/src/aliceVision/fuseCut/DelaunayGraphCut.cpp
@@ -1744,107 +1744,6 @@ void DelaunayGraphCut::fillGraphPartPtRc(int& out_nstepsFront, int& out_nstepsBe
     }
 }
 
-void DelaunayGraphCut::forceTedgesByGradientCVPR11(bool fixesSigma, float nPixelSizeBehind)
-{
-    ALICEVISION_LOG_INFO("Forcing t-edges.");
-    long t2 = clock();
-
-    float forceTEdgeDelta = (float)mp->userParams.get<double>("delaunaycut.forceTEdgeDelta", 0.1f);
-    ALICEVISION_LOG_INFO("forceTEdgeDelta: " << forceTEdgeDelta);
-
-    float beta = (float)mp->userParams.get<double>("delaunaycut.beta", 1000.0f);
-    ALICEVISION_LOG_INFO("beta: " << beta);
-
-    for(GC_cellInfo& c: _cellsAttr)
-    {
-        c.on = 0.0f;
-    }
-
-    // choose random order to prevent waiting
-    const unsigned int seed = (unsigned int)mp->userParams.get<unsigned int>("delaunaycut.seed", 0);
-    const std::vector<int> verticesRandIds = mvsUtils::createRandomArrayOfIntegers(_verticesAttr.size(), seed);
-
-#pragma omp parallel for
-    for(int i = 0; i < verticesRandIds.size(); ++i)
-    {
-        const int vertexIndex = verticesRandIds[i];
-        GC_vertexInfo& v = _verticesAttr[vertexIndex];
-        if(v.isVirtual())
-            continue;
-
-        const Point3d& originPt = _verticesCoords[vertexIndex];
-        for(int c = 0; c < v.cams.size(); ++c)
-        {
-            int cam = v.cams[c];
-
-            // False here mean nearest
-            const bool nearestFarest = true;
-            Facet fFirst = getFacetFromVertexOnTheRayToTheCam(vertexIndex, cam, nearestFarest);
-
-            // get the tetrahedron next to point p on the ray from c
-            // False here mean farest
-            Facet facet = getFacetFromVertexOnTheRayToTheCam(vertexIndex, cam, false);
-
-            if((fFirst.cellIndex != GEO::NO_CELL) && (facet.cellIndex != GEO::NO_CELL) && (!isInfiniteCell(facet.cellIndex)))
-            {
-                float eFirst = _cellsAttr[fFirst.cellIndex].emptinessScore;
-
-                Point3d p = originPt; // HAS TO BE HERE !!!
-                float maxDist = nPixelSizeBehind * mp->getCamPixelSize(p, cam);
-                if(fixesSigma)
-                {
-                    maxDist = nPixelSizeBehind;
-                }
-
-                bool ok = (facet.cellIndex != GEO::NO_CELL);
-                while(ok)
-                {
-                    Point3d pold = p;
-                    Point3d intersectPt;
-                    Facet outFacet;
-
-                    // Intersection with the next facet in the current tetrahedron (ci) in order to find the cell farest
-                    // to the cam which is intersected with cam-p ray
-                    // False here mean farest
-                    const bool nearestFarest = false;
-                    if(!rayCellIntersection(mp->CArr[cam], p, facet, outFacet, nearestFarest, intersectPt) ||
-                       ((originPt - pold).size() >= maxDist))
-                    {
-                        ok = false;
-                    }
-                    else
-                    {
-                        // Take the mirror facet to iterate over the next cell
-                        facet = mirrorFacet(outFacet);
-                        if(facet.cellIndex == GEO::NO_CELL)
-                            ok = false;
-                        p = intersectPt;
-                    }
-                }
-
-                if(facet.cellIndex != GEO::NO_CELL)
-                {
-                    const float eLast = _cellsAttr[facet.cellIndex].emptinessScore;
-                    if((eFirst > eLast) && (eFirst < beta) && (eLast / eFirst < forceTEdgeDelta))
-                    {
-#pragma OMP_ATOMIC_UPDATE
-                        _cellsAttr[facet.cellIndex].on += (eFirst - eLast);
-                    }
-                }
-            }
-
-        } // for c
-
-    } // for i
-
-    for(GC_cellInfo& c: _cellsAttr)
-    {
-        c.cellTWeight = std::max(c.cellTWeight, std::min(1000000.0f, std::max(1.0f, c.cellTWeight) * c.on));
-    }
-
-    mvsUtils::printfElapsedTime(t2, "t-edges forced: ");
-}
-
 void DelaunayGraphCut::forceTedgesByGradientIJCV(bool fixesSigma, float nPixelSizeBehind)
 {
     ALICEVISION_LOG_INFO("Forcing t-edges");
@@ -2460,32 +2359,9 @@ void DelaunayGraphCut::voteFullEmptyScore(const StaticVector<int>& cams, const s
     ////////////////////////////////////////////////////////////////////////////////////////////////////////
 
     bool labatutCFG09 = mp->userParams.get<bool>("global.LabatutCFG09", false);
-    bool jancosekCVPR11 = mp->userParams.get<bool>("global.JancosekCVPR11", false);
     // jancosekIJCV: "Exploiting Visibility Information in Surface Reconstruction to Preserve Weakly Supported Surfaces", Michal Jancosek and Tomas Pajdla, 2014
     bool jancosekIJCV = mp->userParams.get<bool>("global.JancosekIJCV", true);
 
-    if(jancosekCVPR11)
-    {
-        float delta = (float)mp->userParams.get<double>("delaunaycut.forceTEdgeDelta", 0.1f);
-
-        ALICEVISION_LOG_INFO("Jancosek CVPR 2011 method ( forceTEdgeDelta*100 = " << static_cast<int>(delta * 100.0f) << "):");
-
-        fillGraph(false, sigma, false, true, distFcnHeight);
-
-        addToInfiniteSw((float)maxint);
-
-        if(saveTemporaryBinFiles)
-            saveDhInfo(folderName + "delaunayTriangulationInfoInit.bin");
-
-        if((delta > 0.0f) && (delta < 1.0f))
-        {
-            forceTedgesByGradientCVPR11(false, sigma);
-        }
-
-        if(saveTemporaryBinFiles)
-            saveDhInfo(folderName + "delaunayTriangulationInfoAfterForce.bin");
-    }
-
     if(jancosekIJCV) // true by default
     {
         float forceTEdgeDelta = (float)mp->userParams.get<double>("delaunaycut.forceTEdgeDelta", 0.1f);
diff --git a/src/aliceVision/fuseCut/DelaunayGraphCut.hpp b/src/aliceVision/fuseCut/DelaunayGraphCut.hpp
index 22053d6aa5..c9333df241 100644
--- a/src/aliceVision/fuseCut/DelaunayGraphCut.hpp
+++ b/src/aliceVision/fuseCut/DelaunayGraphCut.hpp
@@ -316,7 +316,6 @@ class DelaunayGraphCut
                            bool fixesSigma, float nPixelSizeBehind, bool fillOut,
                            float distFcnHeight);
 
-    void forceTedgesByGradientCVPR11(bool fixesSigma, float nPixelSizeBehind);
     void forceTedgesByGradientIJCV(bool fixesSigma, float nPixelSizeBehind);
 
     int setIsOnSurface();
diff --git a/src/aliceVision/fuseCut/ReconstructionPlan.cpp b/src/aliceVision/fuseCut/ReconstructionPlan.cpp
index 9dc7cd4f56..f52c054722 100644
--- a/src/aliceVision/fuseCut/ReconstructionPlan.cpp
+++ b/src/aliceVision/fuseCut/ReconstructionPlan.cpp
@@ -432,10 +432,6 @@ mesh::Mesh* joinMeshes(int gl, LargeScale* ls)
         {
             subFolderName = "LabatutCFG09";
         }
-        if(ls->mp->userParams.get<bool>("global.JancosekCVPR11", true))
-        {
-            subFolderName = "JancosekCVPR11";
-        }
     }
     subFolderName = subFolderName + "/";
 

From de2a440375f78b352d1ff6064d941323df6481c1 Mon Sep 17 00:00:00 2001
From: Fabien Castan <fabcastan@gmail.com>
Date: Wed, 29 Jul 2020 12:28:29 +0200
Subject: [PATCH 28/28] [fuseCut] add margin in divideSpaceFromSfM

The function "estimateDimensions" is no more called, and it was
indirectly adding a small margin.
---
 src/aliceVision/fuseCut/Fuser.cpp | 24 ++++++++++++++++++------
 1 file changed, 18 insertions(+), 6 deletions(-)

diff --git a/src/aliceVision/fuseCut/Fuser.cpp b/src/aliceVision/fuseCut/Fuser.cpp
index 0fa242cb8e..1fc6da0b30 100644
--- a/src/aliceVision/fuseCut/Fuser.cpp
+++ b/src/aliceVision/fuseCut/Fuser.cpp
@@ -629,12 +629,24 @@ void Fuser::divideSpaceFromSfM(const sfmData::SfMData& sfmData, Point3d* hexah,
     accMaxZ(z);
   }
 
-  const double xMin = quantile(accMinX, quantile_probability = 1.0 - percentile);
-  const double yMin = quantile(accMinY, quantile_probability = 1.0 - percentile);
-  const double zMin = quantile(accMinZ, quantile_probability = 1.0 - percentile);
-  const double xMax = quantile(accMaxX, quantile_probability = percentile);
-  const double yMax = quantile(accMaxY, quantile_probability = percentile);
-  const double zMax = quantile(accMaxZ, quantile_probability = percentile);
+  // Remove a percentile of the observations (to remove unstable points)
+  double xMin = quantile(accMinX, quantile_probability = 1.0 - percentile);
+  double yMin = quantile(accMinY, quantile_probability = 1.0 - percentile);
+  double zMin = quantile(accMinZ, quantile_probability = 1.0 - percentile);
+  double xMax = quantile(accMaxX, quantile_probability = percentile);
+  double yMax = quantile(accMaxY, quantile_probability = percentile);
+  double zMax = quantile(accMaxZ, quantile_probability = percentile);
+
+  // Add a margin on the result
+  const double xMargin = (xMax - xMin) * 0.05;
+  const double yMargin = (yMax - yMin) * 0.05;
+  const double zMargin = (zMax - zMin) * 0.05;
+  xMin -= xMargin;
+  yMin -= yMargin;
+  zMin -= zMargin;
+  xMax += xMargin;
+  yMax += yMargin;
+  zMax += zMargin;
 
   hexah[0] = Point3d(xMax, yMax, zMax);
   hexah[1] = Point3d(xMin, yMax, zMax);