From 88443cada841979c36b81225bb7db19944aef523 Mon Sep 17 00:00:00 2001 From: Enguerrand DE SMET Date: Tue, 28 Jul 2020 18:56:46 +0200 Subject: [PATCH 01/26] [fuseCut] add struct for advanced geometry intersection --- src/aliceVision/fuseCut/DelaunayGraphCut.hpp | 50 ++++++++++++++++++++ 1 file changed, 50 insertions(+) diff --git a/src/aliceVision/fuseCut/DelaunayGraphCut.hpp b/src/aliceVision/fuseCut/DelaunayGraphCut.hpp index 10a61afd4c..105e712f99 100644 --- a/src/aliceVision/fuseCut/DelaunayGraphCut.hpp +++ b/src/aliceVision/fuseCut/DelaunayGraphCut.hpp @@ -77,6 +77,56 @@ class DelaunayGraphCut VertexIndex localVertexIndex = GEO::NO_VERTEX; }; + struct Edge + { + Edge() = default; + Edge(VertexIndex v0_, VertexIndex v1_) + : v0{v0_} + , v1{v1_} + {} + + VertexIndex v0 = GEO::NO_VERTEX; + VertexIndex v1 = GEO::NO_VERTEX; + }; + + enum class EGeometryType + { + Vertex, + Edge, + Facet, + None + }; + + struct GeometryIntersection + { + EGeometryType type = EGeometryType::None; + union + { + Facet facet; + VertexIndex vertexIndex; + Edge edge; + }; + GeometryIntersection() {} + explicit GeometryIntersection(const Facet& f) + : facet{f} + , type{EGeometryType::Facet} + {} + explicit GeometryIntersection(const VertexIndex& v) + : vertexIndex{v} + , type{EGeometryType::Vertex} + {} + explicit GeometryIntersection(const Edge& e) + : edge{e} + , type{EGeometryType::Edge} + {} + }; + + enum class EDirection + { + toTheCam, + behindThePoint + }; + mvsUtils::MultiViewParams* mp; GEO::Delaunay_var _tetrahedralization; From 8e44731238798860937423e12976cd1492523fc2 Mon Sep 17 00:00:00 2001 From: Enguerrand DE SMET Date: Tue, 28 Jul 2020 19:14:30 +0200 Subject: [PATCH 02/26] [fuseCut] delaunayGraphCut: Add a function to get neighbouring cells by the edge --- src/aliceVision/fuseCut/DelaunayGraphCut.cpp | 11 +++++++++++ src/aliceVision/fuseCut/DelaunayGraphCut.hpp | 10 +++++++++- 2 files changed, 20 insertions(+), 1 deletion(-) diff --git a/src/aliceVision/fuseCut/DelaunayGraphCut.cpp b/src/aliceVision/fuseCut/DelaunayGraphCut.cpp index e3ecaf247d..6fbd24e1cd 100644 --- a/src/aliceVision/fuseCut/DelaunayGraphCut.cpp +++ b/src/aliceVision/fuseCut/DelaunayGraphCut.cpp @@ -445,6 +445,17 @@ void DelaunayGraphCut::saveDh(const std::string& fileNameDh, const std::string& mvsUtils::printfElapsedTime(t1); } +std::vector DelaunayGraphCut::getNeighboringCellsByEdge(const Edge& e) const +{ + const std::vector& v0ci = getNeighboringCellsByVertexIndex(e.v0); + const std::vector& v1ci = getNeighboringCellsByVertexIndex(e.v1); + + std::vector neighboringCells; + std::set_intersection(v0ci.begin(), v0ci.end(), v1ci.begin(), v1ci.end(), std::back_inserter(neighboringCells)); + + return neighboringCells; +} + void DelaunayGraphCut::initVertices() { ALICEVISION_LOG_DEBUG("initVertices ...\n"); diff --git a/src/aliceVision/fuseCut/DelaunayGraphCut.hpp b/src/aliceVision/fuseCut/DelaunayGraphCut.hpp index 105e712f99..b8862ae342 100644 --- a/src/aliceVision/fuseCut/DelaunayGraphCut.hpp +++ b/src/aliceVision/fuseCut/DelaunayGraphCut.hpp @@ -310,13 +310,21 @@ class DelaunayGraphCut * @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 + * @return a vector of neighboring cell indices */ inline const std::vector& getNeighboringCellsByVertexIndex(VertexIndex vi) const { return _neighboringCellsPerVertex.at(vi); } + /** + * @brief Retrieves the global indexes of neighboring cells around one edge. + * + * @param e the concerned edge + * @return a vector of neighboring cell indices + */ + std::vector getNeighboringCellsByEdge(const Edge& e) const; + void initVertices(); void computeDelaunay(); void initCells(); From d43ed9703f71ea6c3fe577a96f434da3023cfb61 Mon Sep 17 00:00:00 2001 From: Enguerrand DE SMET Date: Tue, 28 Jul 2020 19:22:35 +0200 Subject: [PATCH 03/26] [fuseCut] DelaunayGraphCut: fix the order during the insertion of camera's points --- src/aliceVision/fuseCut/DelaunayGraphCut.cpp | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/aliceVision/fuseCut/DelaunayGraphCut.cpp b/src/aliceVision/fuseCut/DelaunayGraphCut.cpp index 6fbd24e1cd..7c8d111b49 100644 --- a/src/aliceVision/fuseCut/DelaunayGraphCut.cpp +++ b/src/aliceVision/fuseCut/DelaunayGraphCut.cpp @@ -2195,12 +2195,6 @@ void DelaunayGraphCut::createDensePointCloud(Point3d hexah[8], const StaticVecto 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); - // add points from depth maps if(depthMapsFuseParams != nullptr) fuseFromDepthMaps(cams, hexah, *depthMapsFuseParams); @@ -2209,6 +2203,12 @@ void DelaunayGraphCut::createDensePointCloud(Point3d hexah[8], const StaticVecto if(sfmData != nullptr) addPointsFromSfM(hexah, cams, *sfmData); + // add points for cam centers + addPointsFromCameraCenters(cams, minDist); + + // add 6 points to prevent singularities + addPointsToPreventSingularities(hexah, minDist); + const int nGridHelperVolumePointsDim = mp->userParams.get("LargeScale.nGridHelperVolumePointsDim", 10); // add volume points to prevent singularities From 0e0049c016d8c52f3c214e05dadc9606e6cd6f4d Mon Sep 17 00:00:00 2001 From: Enguerrand DE SMET Date: Tue, 28 Jul 2020 19:23:30 +0200 Subject: [PATCH 04/26] [fuseCut] DelaunayGraphCut: Add progress Bar --- src/aliceVision/fuseCut/DelaunayGraphCut.cpp | 12 +++++++++++- 1 file changed, 11 insertions(+), 1 deletion(-) diff --git a/src/aliceVision/fuseCut/DelaunayGraphCut.cpp b/src/aliceVision/fuseCut/DelaunayGraphCut.cpp index 7c8d111b49..6e91793941 100644 --- a/src/aliceVision/fuseCut/DelaunayGraphCut.cpp +++ b/src/aliceVision/fuseCut/DelaunayGraphCut.cpp @@ -28,6 +28,7 @@ #include #include +#include namespace aliceVision { namespace fuseCut { @@ -1573,10 +1574,19 @@ void DelaunayGraphCut::fillGraph(bool fixesSigma, float nPixelSizeBehind, int avCams = 0; int nAvCams = 0; + boost::progress_display progressBar(std::min(size_t(100), verticesRandIds.size()), std::cout, "fillGraphPartPtRc\n"); + size_t progressStep = verticesRandIds.size() / 100; + progressStep = std::max(size_t(1), progressStep); #pragma omp parallel for reduction(+:avStepsFront,aAvStepsFront,avStepsBehind,nAvStepsBehind,avCams,nAvCams) for(int i = 0; i < verticesRandIds.size(); i++) { - int vertexIndex = verticesRandIds[i]; + if(i % progressStep == 0) + { +#pragma omp critical + ++progressBar; + } + + const int vertexIndex = verticesRandIds[i]; const GC_vertexInfo& v = _verticesAttr[vertexIndex]; if(v.isReal() && (v.nrc > 0)) From 91d5f9272cca50a6423f91f9548e8026e21f27ed Mon Sep 17 00:00:00 2001 From: Enguerrand DE SMET Date: Thu, 30 Jul 2020 15:59:50 +0200 Subject: [PATCH 05/26] [fuseCut] delaunayGraphCut: add geometry intersections functions Add new functions to handle intersections with new geometry struct. Allows us to deal with the specific cases of edges and vertices --- src/aliceVision/fuseCut/DelaunayGraphCut.cpp | 156 +++++++++++++++++++ src/aliceVision/fuseCut/DelaunayGraphCut.hpp | 30 ++++ 2 files changed, 186 insertions(+) diff --git a/src/aliceVision/fuseCut/DelaunayGraphCut.cpp b/src/aliceVision/fuseCut/DelaunayGraphCut.cpp index 6e91793941..ff4d5aedae 100644 --- a/src/aliceVision/fuseCut/DelaunayGraphCut.cpp +++ b/src/aliceVision/fuseCut/DelaunayGraphCut.cpp @@ -1370,6 +1370,162 @@ DelaunayGraphCut::Facet DelaunayGraphCut::getFacetFromVertexOnTheRayToTheCam(Ver return Facet(); } + +DelaunayGraphCut::GeometryIntersection +DelaunayGraphCut::intersectNextGeom(const DelaunayGraphCut::GeometryIntersection& inGeometry, + const Point3d& originPt, + const Point3d& dirVect, Point3d& intersectPt, + const float epsilon, const Point3d& lastIntersectPt) const +{ + + switch (inGeometry.type) + { + case EGeometryType::Facet: + { + const CellIndex tetrahedronIndex = inGeometry.facet.cellIndex; + + // Test all facets of the tetrahedron using i as localVertexIndex to define next intersectionFacet + for (int i = 0; i < 4; ++i) + { + // Because we can't intersect with incoming facet (same localVertexIndex) + if (i == inGeometry.facet.localVertexIndex) + continue; + + const Facet intersectionFacet(tetrahedronIndex, i); + + const GeometryIntersection result = rayIntersectTriangle(originPt, dirVect, intersectionFacet, intersectPt, epsilon, &lastIntersectPt); + if (result.type != EGeometryType::None) + return result; + } + } + break; + + case EGeometryType::Vertex: + { + for (CellIndex adjCellIndex : getNeighboringCellsByVertexIndex(inGeometry.vertexIndex)) + { + if (isInfiniteCell(adjCellIndex)) + continue; + + // Get local vertex index + const VertexIndex localVertexIndex = _tetrahedralization->index(adjCellIndex, inGeometry.vertexIndex); + + // Define the facet to intersect + const Facet facet(adjCellIndex, localVertexIndex); + + const GeometryIntersection result = rayIntersectTriangle(originPt, dirVect, facet, intersectPt, epsilon, &lastIntersectPt); + if (result.type != EGeometryType::None) + return result; + } + } + break; + + case EGeometryType::Edge: + { + GeometryIntersection result; + + for (CellIndex adjCellIndex : getNeighboringCellsByEdge(inGeometry.edge)) + { + if (isInfiniteCell(adjCellIndex)) + continue; + // Local vertices indices + const VertexIndex lvi0 = _tetrahedralization->index(adjCellIndex, inGeometry.edge.v0); + const VertexIndex lvi1 = _tetrahedralization->index(adjCellIndex, inGeometry.edge.v1); + + // The two facets that do not touch this edge need to be tested + const std::array opositeFacets{ {{adjCellIndex, lvi0}, {adjCellIndex, lvi1}} }; + + for (const Facet& facet : opositeFacets) + { + const GeometryIntersection result = rayIntersectTriangle(originPt, dirVect, facet, intersectPt, epsilon, &lastIntersectPt); + if (result.type != EGeometryType::None) + return result; + } + } + } + break; + + case EGeometryType::None: + throw std::runtime_error("[intersectNextGeom] intersection with input none geometry should not happen."); + break; + } + return GeometryIntersection(); +} + +DelaunayGraphCut::GeometryIntersection DelaunayGraphCut::rayIntersectTriangle(const Point3d& originPt, + const Point3d& DirVec, + const DelaunayGraphCut::Facet& facet, + Point3d& intersectPt, + const float epsilon, const Point3d* lastIntersectPt) const +{ + const VertexIndex AvertexIndex = getVertexIndex(facet, 0); + const VertexIndex BvertexIndex = getVertexIndex(facet, 1); + const VertexIndex CvertexIndex = getVertexIndex(facet, 2); + + const Point3d* A = &_verticesCoords[AvertexIndex]; + const Point3d* B = &_verticesCoords[BvertexIndex]; + const Point3d* C = &_verticesCoords[CvertexIndex]; + + Point3d tempIntersectPt; + const Point2d triangleUv = getLineTriangleIntersectBarycCoords(&tempIntersectPt, A, B, C, &originPt, &DirVec); + + if (!isnormal(tempIntersectPt.x) || !isnormal(tempIntersectPt.y) || !isnormal(tempIntersectPt.z)) + { + ALICEVISION_LOG_DEBUG("Invalid/notNormal intersection point found during rayIntersectTriangle."); + } + + const float u = triangleUv.x; // A to C + const float v = triangleUv.y; // A to B + + // If we find invalid uv coordinate + if (!std::isfinite(u) && !std::isfinite(v)) + return GeometryIntersection(); + + // Ouside the triangle with epsilon margin + if (u < -epsilon || v < -epsilon || (u + v) >(1 + epsilon)) + return GeometryIntersection(); + + // In case intersectPt is provided, check if intersectPt is in front of lastIntersectionPt + // in the DirVec direction to ensure that we are moving forward in the right direction + if (lastIntersectPt != nullptr && dot(DirVec, (tempIntersectPt - *lastIntersectPt).normalize()) <= 0) + return GeometryIntersection(); + + // Change intersection point only if tempIntersectionPt is in the right direction (mean we intersect something) + intersectPt = tempIntersectPt; + + if (v < epsilon) // along A C edge + { + if (u < epsilon) + { + intersectPt = *A; + return GeometryIntersection(AvertexIndex); // vertex A + } + if (u > 1.0f - epsilon) + { + intersectPt = *C; + return GeometryIntersection(CvertexIndex); // vertex C + } + + return GeometryIntersection(Edge(AvertexIndex, CvertexIndex)); // edge AC + } + + if (u < epsilon) // along A B edge + { + if (v > 1.0f - epsilon) + { + intersectPt = *B; + return GeometryIntersection(BvertexIndex); // vertex B + } + + return GeometryIntersection(Edge(AvertexIndex, BvertexIndex)); // edge AB + } + + if (u + v > 1.0f - epsilon) + return GeometryIntersection(Edge(BvertexIndex, CvertexIndex)); // edge BC + + return GeometryIntersection(facet); +} + DelaunayGraphCut::Facet DelaunayGraphCut::getFirstFacetOnTheRayFromCamToThePoint(int cam, const Point3d& p, Point3d& intersectPt) const { int cam_vi = _camsVertexes[cam]; diff --git a/src/aliceVision/fuseCut/DelaunayGraphCut.hpp b/src/aliceVision/fuseCut/DelaunayGraphCut.hpp index b8862ae342..1c6f767acc 100644 --- a/src/aliceVision/fuseCut/DelaunayGraphCut.hpp +++ b/src/aliceVision/fuseCut/DelaunayGraphCut.hpp @@ -360,6 +360,36 @@ class DelaunayGraphCut Facet getFacetFromVertexOnTheRayToTheCam(VertexIndex globalVertexIndex, int cam, bool nearestFarest) const; Facet getFirstFacetOnTheRayFromCamToThePoint(int cam, const Point3d& p, Point3d& intersectPt) const; + /** + * @brief Function that returns the next geometry intersected by the ray. + * The function handles different cases, whether we come from an edge, a facet or a vertex. + * + * @param inGeometry the geometry we come from + * @param originPt ray origin point + * @param dirVect ray direction + * @param intersectPt a reference that will store the computed intersection point for the intersected geometry + * @param epsilon used to define the boundary when we have to consider either a collision with an edge/vertex or a facet + * @param lastIntersectPt constant reference to the last intersection point used to test the direction. + * @return + */ + GeometryIntersection intersectNextGeom(const GeometryIntersection& inGeometry, const Point3d& originPt, + const Point3d& dirVect, Point3d& intersectPt, const float epsilon, const Point3d& lastIntersectPt) const; + + /** + * @brief Function that returns the next geometry intersected by the ray on a given facet or None if there are no intersected geometry. + * The function distinguishes the intersections cases using epsilon. + * + * @param originPt ray origin point + * @param DirVec ray direction + * @param facet the given facet to intersect with + * @param intersectPt a reference that will store the computed intersection point for the next intersecting geometry + * @param epsilon used to define the boundary when we have to consider either a collision with an edge/vertex or a facet + * @param lastIntersectPt pointer to the last intersection point used to test the direction (if not nulllptr) + * @return + */ + GeometryIntersection rayIntersectTriangle(const Point3d& originPt, const Point3d& DirVec, const Facet& facet, + Point3d& intersectPt, const float epsilon, const Point3d* lastIntersectPt = nullptr) const; + float distFcn(float maxDist, float dist, float distFcnHeight) const; inline double conj(double val) const { return val; } From 4fe8bbb2e919a6b2943c2529d3f435bf653a17bc Mon Sep 17 00:00:00 2001 From: Enguerrand DE SMET Date: Thu, 30 Jul 2020 16:35:07 +0200 Subject: [PATCH 06/26] [fuseCut] DelaunayGraphCut: new implementation using geometry (CRITIC) --- src/aliceVision/fuseCut/DelaunayGraphCut.cpp | 374 ++++++++++--------- 1 file changed, 201 insertions(+), 173 deletions(-) diff --git a/src/aliceVision/fuseCut/DelaunayGraphCut.cpp b/src/aliceVision/fuseCut/DelaunayGraphCut.cpp index ff4d5aedae..8d1186bb09 100644 --- a/src/aliceVision/fuseCut/DelaunayGraphCut.cpp +++ b/src/aliceVision/fuseCut/DelaunayGraphCut.cpp @@ -1783,152 +1783,167 @@ void DelaunayGraphCut::fillGraphPartPtRc(int& out_nstepsFront, int& out_nstepsBe 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; - out_nstepsBehind = 0; - - int maxint = 1000000; // std::numeric_limits::std::max() + const int maxint = 1000000; // std::numeric_limits::std::max() const Point3d& originPt = _verticesCoords[vertexIndex]; const float pixSize = mp->getCamPixelSize(originPt, cam); float maxDist = nPixelSizeBehind * pixSize; if(fixesSigma) - { maxDist = nPixelSizeBehind; - } assert(cam >= 0); assert(cam < mp->ncams); - // printf("%f %f %f\n",originPt.x,originPt.y,originPt.z); - int nsteps = 0; if(fillOut) { + // Initialisation + GeometryIntersection geometry(vertexIndex); // Starting on global vertex index + Point3d intersectPt = originPt; + const EDirection dir = EDirection::toTheCam; + const Point3d dirVect = (mp->CArr[cam] - originPt).normalize() * (dir == EDirection::toTheCam ? 1.0 : -1.0); + out_nstepsFront = 0; - // true here mean nearest - const bool nearestFarest = true; - // 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 = facet.cellIndex != GEO::NO_CELL; - while(ok) + Facet lastIntersectedFacet; + // As long as we find a next geometry + do { - { -#pragma OMP_ATOMIC_UPDATE - _cellsAttr[facet.cellIndex].emptinessScore += weight; - } + // Keep previous informations + const GeometryIntersection previousGeometry = geometry; + const Point3d lastIntersectPt = intersectPt; ++out_nstepsFront; - ++nsteps; - - Facet outFacet; - 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, facet, outFacet, nearestFarest, intersectPt)) + geometry = intersectNextGeom(previousGeometry, originPt, dirVect, intersectPt, 1.0e-3, lastIntersectPt); + + if (geometry.type == EGeometryType::None) { - ok = false; + // throw std::runtime_error("[Error]: fillGraphPartPtRc toTheCam, cause: geometry cannot be found."); + ALICEVISION_LOG_TRACE("[Error]: fillGraphPartPtRc toTheCam, cause: geometry cannot be found."); + break; } - else + + if(geometry.type == EGeometryType::Facet) { + lastIntersectedFacet = geometry.facet; + // Vote for the cell between the previous vertex/edge and the intersected facet + // Because if there is an edge or a vertex, no vote has been made, we have to vote here for the right facet once we find it + if(previousGeometry.type != EGeometryType::Facet) { - const float dist = distFcn(maxDist, (originPt - p).size(), distFcnHeight); -#pragma OMP_ATOMIC_UPDATE - _cellsAttr[outFacet.cellIndex].gEdgeVisWeight[outFacet.localVertexIndex] += weight * dist; + _cellsAttr[geometry.facet.cellIndex].emptinessScore += weight; } + const float dist = distFcn(maxDist, (originPt - intersectPt).size(), distFcnHeight); +#pragma OMP_ATOMIC_UPDATE + _cellsAttr[geometry.facet.cellIndex].gEdgeVisWeight[geometry.facet.localVertexIndex] += weight * dist; + // Take the mirror facet to iterate over the next cell - facet = mirrorFacet(outFacet); + const Facet mFacet = mirrorFacet(geometry.facet); + if (isInvalidOrInfiniteCell(mFacet.cellIndex)) + { + // throw std::runtime_error("[Error]: fillGraphPartPtRc toTheCam, cause: invalidOrInfinite miror facet."); + ALICEVISION_LOG_TRACE("[Error]: fillGraphPartPtRc toTheCam, cause: invalidOrInfinite miror facet."); + break; + } + geometry.facet = mFacet; - if(facet.cellIndex == GEO::NO_CELL) - ok = false; - p = intersectPt; + // Vote for the new cell after taking the mirror facet +#pragma OMP_ATOMIC_UPDATE + _cellsAttr[geometry.facet.cellIndex].emptinessScore += weight; } - } - // get the outer tetrahedron of camera c for the ray to p = the last tetrahedron - if(facet.cellIndex != GEO::NO_CELL) + // Break only when we reach our camera vertex + } while(!(geometry.type == EGeometryType::Vertex && (mp->CArr[cam] - intersectPt).size() < 1.0e-3)); + + // Vote for the last intersected facet (close to the cam) + if (lastIntersectedFacet.cellIndex != GEO::NO_CELL) { #pragma OMP_ATOMIC_WRITE - _cellsAttr[facet.cellIndex].cellSWeight = (float)maxint; + _cellsAttr[lastIntersectedFacet.cellIndex].cellSWeight = (float)maxint; } } - { + // Initialisation + GeometryIntersection geometry(vertexIndex); // Starting on global vertex index + Point3d intersectPt = originPt; + const EDirection dir = EDirection::behindThePoint; + const Point3d dirVect = (mp->CArr[cam] - originPt).normalize() * (dir == EDirection::toTheCam ? 1 : -1); + out_nstepsBehind = 0; - // get the tetrahedron next to point p on the ray from c - // False here mean farest - const bool nearestFarest = false; - Facet facet = getFacetFromVertexOnTheRayToTheCam(vertexIndex, cam, nearestFarest); - Facet outFacet; - if(facet.cellIndex != GEO::NO_CELL) + bool firstIteration = true; + Facet lastIntersectedFacet; + // As long as we find a next geometry + do { -#pragma OMP_ATOMIC_UPDATE - _cellsAttr[facet.cellIndex].on += weight; - } + // Keep previous informations + const GeometryIntersection previousGeometry = geometry; + const Point3d lastIntersectPt = intersectPt; - Point3d p = originPt; // HAS TO BE HERE !!! + ++out_nstepsBehind; - bool ok = (facet.cellIndex != GEO::NO_CELL); - while(ok) - { - GC_cellInfo& c = _cellsAttr[facet.cellIndex]; + // Vote for the first facet found (only once) + if(firstIteration && geometry.type == EGeometryType::Facet) { #pragma OMP_ATOMIC_UPDATE - c.fullnessScore += weight; + _cellsAttr[geometry.facet.cellIndex].on += weight; + firstIteration = false; } + geometry = intersectNextGeom(previousGeometry, originPt, dirVect, intersectPt, 1.0e-3, lastIntersectPt); - ++out_nstepsBehind; - ++nsteps; - - Point3d intersectPt; - - // 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, facet, outFacet, nearestFarest, intersectPt) || - ((originPt - p).size() >= maxDist)) + if(geometry.type == EGeometryType::None) { - ok = false; + // If we come from a facet, the next intersection must exist (even if the mirror facet is invalid, which is verified after taking mirror facet) + if(previousGeometry.type == EGeometryType::Facet) + { + // throw std::runtime_error("[DelaunayGraphCut::fillGraphPartPtRc] first loop: Geometry cannot be found."); + ALICEVISION_LOG_TRACE("[Error]: fillGraphPartPtRc behindThePoint, cause: None geometry but previous is Facet."); + } + // Break if we reach the end of the tetrahedralization volume + break; } - else + + if(geometry.type == EGeometryType::Facet) { - // float dist = 1.0f; - // 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 + lastIntersectedFacet = geometry.facet; + // Vote for the cell between the previous vertex/edge and the intersected facet + // Because if there is an edge or a vertex, no vote has been made, we have to vote here for the right facet once we find it + if (previousGeometry.type != EGeometryType::Facet) + { +#pragma OMP_ATOMIC_UPDATE + _cellsAttr[geometry.facet.cellIndex].fullnessScore += weight; + } + const float dist = distFcn(maxDist, (originPt - intersectPt).size(), distFcnHeight); +#pragma OMP_ATOMIC_UPDATE + _cellsAttr[geometry.facet.cellIndex].gEdgeVisWeight[geometry.facet.localVertexIndex] += weight * dist; + // Take the mirror facet to iterate over the next cell - facet = mirrorFacet(outFacet); - if(facet.cellIndex == GEO::NO_CELL) + const Facet mFacet = mirrorFacet(geometry.facet); + if (isInvalidOrInfiniteCell(mFacet.cellIndex)) { - ok = false; + // Break if we reach the end of the tetrahedralization volume (mirror facet cannot be found) + break; } - else - { - const float dist = distFcn(maxDist, (originPt - p).size(), distFcnHeight); + geometry.facet = mFacet; + + // Vote for the new cell after taking the mirror facet #pragma OMP_ATOMIC_UPDATE - _cellsAttr[facet.cellIndex].gEdgeVisWeight[facet.localVertexIndex] += weight * dist; - } - p = intersectPt; + _cellsAttr[geometry.facet.cellIndex].fullnessScore += weight; } - } - // cv: is the tetrahedron in distance 2*sigma behind the point p in the direction of the camera c (called Lcp in the paper) + // While we are within the surface margin + } while((originPt - intersectPt).size() < maxDist); - if(facet.cellIndex != GEO::NO_CELL) + // cv: is the tetrahedron in distance 2*sigma behind the point p in the direction of the camera c (called Lcp in the paper) Vote for the last found facet + // Vote for the last intersected facet (farthest from the camera) + if (lastIntersectedFacet.cellIndex != GEO::NO_CELL) { -#pragma OMP_ATOMIC_UPDATE - _cellsAttr[facet.cellIndex].cellTWeight += weight; +#pragma OMP_ATOMIC_WRITE + _cellsAttr[lastIntersectedFacet.cellIndex].cellTWeight += weight; } } } @@ -1977,26 +1992,15 @@ void DelaunayGraphCut::forceTedgesByGradientIJCV(bool fixesSigma, float nPixelSi for(int i = 0; i < verticesRandIds.size(); ++i) { const int vertexIndex = verticesRandIds[i]; - GC_vertexInfo& v = _verticesAttr[vertexIndex]; + const GC_vertexInfo& v = _verticesAttr[vertexIndex]; if(v.isVirtual()) continue; const Point3d& originPt = _verticesCoords[vertexIndex]; - for(int c = 0; c < v.cams.size(); ++c) + // For each camera that has visibility over the vertex v (vertexIndex) + for(const int cam : v.cams) { - int nstepsFront = 0; - int nstepsBehind = 0; - - int cam = v.cams[c]; - float maxDist = 0.0f; - if(fixesSigma) - { - maxDist = nPixelSizeBehind; - } - else - { - maxDist = nPixelSizeBehind * mp->getCamPixelSize(originPt, cam); - } + const float maxDist = nPixelSizeBehind * (fixesSigma ? 1.0f : mp->getCamPixelSize(originPt, cam)); float minJump = 10000000.0f; float minSilent = 10000000.0f; @@ -2005,91 +2009,117 @@ void DelaunayGraphCut::forceTedgesByGradientIJCV(bool fixesSigma, float nPixelSi float midSilent = 10000000.0f; { - // True here mean nearest - const bool nearestFarest = true; - Facet facet = getFacetFromVertexOnTheRayToTheCam(vertexIndex, cam, nearestFarest); - Point3d p = originPt; // HAS TO BE HERE !!! - bool ok = (facet.cellIndex != GEO::NO_CELL); - while(ok) + // Initialisation + GeometryIntersection geometry(vertexIndex); // Starting on global vertex index + Point3d intersectPt = originPt; + const EDirection dir = EDirection::toTheCam; + const Point3d dirVect = (mp->CArr[cam] - originPt).normalize() * (dir == EDirection::toTheCam ? 1.0 : -1.0); + + int nstepsFront = 0; + + // As long as we find a next geometry + do { + // Keep previous informations + const GeometryIntersection previousGeometry = geometry; + const Point3d lastIntersectPt = intersectPt; + ++nstepsFront; - const GC_cellInfo& c = _cellsAttr[facet.cellIndex]; - if((p - originPt).size() > nsigmaFrontSilentPart * maxDist) // (p-originPt).size() > 2 * sigma - { - minJump = std::min(minJump, c.emptinessScore); - maxJump = std::max(maxJump, c.emptinessScore); - } - else - { - minSilent = std::min(minSilent, c.emptinessScore); - maxSilent = std::max(maxSilent, c.emptinessScore); - } + geometry = intersectNextGeom(previousGeometry, originPt, dirVect, intersectPt, 1.0e-3, lastIntersectPt); - Facet outFacet; - Point3d intersectPt; - // 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, facet, outFacet, nearestFarest, intersectPt)) + if (geometry.type == EGeometryType::None) { - ok = false; + // throw std::runtime_error("[DelaunayGraphCut::fillGraphPartPtRc] first loop: Geometry cannot be found."); + ALICEVISION_LOG_TRACE("[Error]: forceTedgesByGradientIJCV toTheCam, cause: geometry cannot be found."); + break; } - else + + if(geometry.type == EGeometryType::Facet) { + const GC_cellInfo& c = _cellsAttr[geometry.facet.cellIndex]; + if((intersectPt - originPt).size() > nsigmaFrontSilentPart * maxDist) // (p-originPt).size() > 2 * sigma + { + minJump = std::min(minJump, c.emptinessScore); + maxJump = std::max(maxJump, c.emptinessScore); + } + else + { + minSilent = std::min(minSilent, c.emptinessScore); + maxSilent = std::max(maxSilent, c.emptinessScore); + } + // Take the mirror facet to iterate over the next cell - facet = mirrorFacet(outFacet); - if(facet.cellIndex == GEO::NO_CELL) - ok = false; - p = intersectPt; + const Facet mFacet = mirrorFacet(geometry.facet); + if (isInvalidOrInfiniteCell(mFacet.cellIndex)) + { + // throw std::runtime_error("[DelaunayGraphCut::fillGraphPartPtRc] first loop: Geometry cannot be found."); + ALICEVISION_LOG_TRACE("[Error]: forceTedgesByGradientIJCV toTheCam, cause: invalidOrInfinite miror facet."); + break; + } + geometry.facet = mFacet; } - } - } + // Break only when we reach our camera vertex or our distance condition is unsatisfied + } while(!(geometry.type == EGeometryType::Vertex && (mp->CArr[cam] - intersectPt).size() < 1.0e-3) + && (intersectPt - originPt).size() <= (nsigmaJumpPart + nsigmaFrontSilentPart) * maxDist); + + avStepsFront += nstepsFront; + aAvStepsFront += 1; + } { - // False here mean farest - const bool nearestFarest = false; - Facet facet = getFacetFromVertexOnTheRayToTheCam(vertexIndex, cam, nearestFarest); - Point3d p = originPt; // HAS TO BE HERE !!! - bool ok = (facet.cellIndex != GEO::NO_CELL); - if(ok) - { - midSilent = _cellsAttr[facet.cellIndex].emptinessScore; - } + // Initialisation + GeometryIntersection geometry(vertexIndex); + Point3d intersectPt = originPt; + const EDirection dir = EDirection::behindThePoint; + const Point3d dirVect = (mp->CArr[cam] - originPt).normalize() * (dir == EDirection::toTheCam ? 1 : -1); - while(ok) + int nstepsBehind = 0; + + Facet lastIntersectedFacet; + // As long as we find a next geometry + do { - nstepsBehind++; - const GC_cellInfo& c = _cellsAttr[facet.cellIndex]; + // Keep previous informations + const GeometryIntersection previousGeometry = geometry; + const Point3d lastIntersectPt = intersectPt; - minSilent = std::min(minSilent, c.emptinessScore); - maxSilent = std::max(maxSilent, c.emptinessScore); + ++nstepsBehind; - Facet outFacet; - Point3d intersectPt; + geometry = intersectNextGeom(previousGeometry, originPt, dirVect, intersectPt, 1.0e-3, lastIntersectPt); - // 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 - originPt).size() > nsigmaBackSilentPart * maxDist) || // (p-originPt).size() > 2 * sigma - !rayCellIntersection(mp->CArr[cam], p, facet, outFacet, nearestFarest, intersectPt)) + if(geometry.type == EGeometryType::None) { - ok = false; + lastIntersectedFacet = geometry.facet; + // If we come from a facet, the next intersection must exist (even if the mirror facet is invalid, which is verified later) + if (previousGeometry.type == EGeometryType::Facet) + { + // throw std::runtime_error("[Error]: forceTedgesByGradientIJCV behindThePoint, cause: geometry cannot be found."); + ALICEVISION_LOG_TRACE("[Error]: forceTedgesByGradientIJCV behindThePoint, cause: geometry cannot be found."); + } + // Break if we reach the end of the tetrahedralization volume + break; } - else + + if(geometry.type == EGeometryType::Facet) { + const GC_cellInfo& c = _cellsAttr[geometry.facet.cellIndex]; + minSilent = std::min(minSilent, c.emptinessScore); + maxSilent = std::max(maxSilent, c.emptinessScore); + // Take the mirror facet to iterate over the next cell - facet = mirrorFacet(outFacet); - if(facet.cellIndex == GEO::NO_CELL) - ok = false; - p = intersectPt; + const Facet mFacet = mirrorFacet(geometry.facet); + if (isInvalidOrInfiniteCell(mFacet.cellIndex)) + { + // Break if we reach the end of the tetrahedralization volume (mirror facet cannot be found) + break; + } + geometry.facet = mFacet; } - } - if(facet.cellIndex != GEO::NO_CELL) + } while((intersectPt - originPt).size() <= nsigmaBackSilentPart * maxDist); + + if (lastIntersectedFacet.cellIndex != GEO::NO_CELL) { // Equation 6 in paper // (g / B) < k_rel @@ -2112,15 +2142,13 @@ void DelaunayGraphCut::forceTedgesByGradientIJCV(bool fixesSigma, float nPixelSi //(maxSilent-minSilent Date: Thu, 30 Jul 2020 16:42:21 +0200 Subject: [PATCH 07/26] [fuseCut] DelaunayGraphCut: remove old implementation --- src/aliceVision/fuseCut/DelaunayGraphCut.cpp | 135 ------------------- src/aliceVision/fuseCut/DelaunayGraphCut.hpp | 6 - 2 files changed, 141 deletions(-) diff --git a/src/aliceVision/fuseCut/DelaunayGraphCut.cpp b/src/aliceVision/fuseCut/DelaunayGraphCut.cpp index 8d1186bb09..365faa380e 100644 --- a/src/aliceVision/fuseCut/DelaunayGraphCut.cpp +++ b/src/aliceVision/fuseCut/DelaunayGraphCut.cpp @@ -1271,106 +1271,6 @@ void DelaunayGraphCut::removeSmallSegs(int minSegSize) initVertices(); } -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(); - double currentDist = (camCenter - p).size(); - - // go thru all triangles - bool existsTriOnRay = false; - 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, 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 - }}; - GEO::index_t oppositeVertexIndex = GEO::NO_VERTEX; - Point3d intersectPt; - double intersectDist; - // Test all facets of the tetrahedron - for(int i = 0; i < facetsPoints.size(); ++i) - { - - 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 - // behind the point (from the camera) if nearestFarest == false - if(nearestFarest ? (intersectDist < currentDist) : (intersectDist > currentDist)) - { - oppositeVertexIndex = i; - existsTriOnRay = true; - currentDist = intersectDist; - outIntersectPt = intersectPt; - } - } - } - - if(!existsTriOnRay) - return false; - - outFacet.cellIndex = tetrahedronIndex; // cell handle - outFacet.localVertexIndex = oppositeVertexIndex; // vertex index that is not in the intersected facet - return true; -} - -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: " << globalVertexIndex); - } - const Point3d& p = _verticesCoords[globalVertexIndex]; - - double currentDist = (mp->CArr[cam] - p).size(); // initialize currentDist to the distance from 3d point p to camera center - - for(CellIndex adjCellIndex : getNeighboringCellsByVertexIndex(globalVertexIndex)) // GEOGRAM: set_stores_cicl(true) required - { - if(isInfiniteCell(adjCellIndex)) - continue; - - // Get local vertex index - const VertexIndex localVertexIndex = _tetrahedralization->index(adjCellIndex, globalVertexIndex); - - // Define facet - const Facet facet(adjCellIndex, localVertexIndex); - - const std::array facetPoints = getFacetsPoints(facet); - - Point3d intersectPt; - 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)) - { - return facet; - } - } - } - return Facet(); -} - - DelaunayGraphCut::GeometryIntersection DelaunayGraphCut::intersectNextGeom(const DelaunayGraphCut::GeometryIntersection& inGeometry, const Point3d& originPt, @@ -1526,41 +1426,6 @@ DelaunayGraphCut::GeometryIntersection DelaunayGraphCut::rayIntersectTriangle(co return GeometryIntersection(facet); } -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; - Facet farestFacet; - - for(CellIndex adjCellIndex : getNeighboringCellsByVertexIndex(cam_vi)) // GEOGRAM: set_stores_cicl(true) required - { - if(isInfiniteCell(adjCellIndex)) - continue; - - // Get local vertex index - const VertexIndex localVertexIndex = _tetrahedralization->index(adjCellIndex, cam_vi); - - // Define facet - const Facet facet(adjCellIndex, localVertexIndex); - - const std::array facetPoints = getFacetsPoints(facet); - - const Point3d lineVect = (camBehind - mp->CArr[cam]).normalize(); - if(isLineInTriangle(&intersectPt, facetPoints[0], facetPoints[1], facetPoints[2], &mp->CArr[cam], &lineVect)) - { - const float dist = orientedPointPlaneDistance(intersectPt, camBehind, dir); - if(dist > maxdist) - { - maxdist = dist; - farestFacet = facet; - } - } - } - return farestFacet; -} - float DelaunayGraphCut::distFcn(float maxDist, float dist, float distFcnHeight) const { // distFcnHeight ... 0 for distFcn == 1 for all dist, 0.1 distFcn std::min diff --git a/src/aliceVision/fuseCut/DelaunayGraphCut.hpp b/src/aliceVision/fuseCut/DelaunayGraphCut.hpp index 1c6f767acc..4fd81677a2 100644 --- a/src/aliceVision/fuseCut/DelaunayGraphCut.hpp +++ b/src/aliceVision/fuseCut/DelaunayGraphCut.hpp @@ -354,12 +354,6 @@ class DelaunayGraphCut void computeVerticesSegSize(bool allPoints, float alpha = 0.0f); void removeSmallSegs(int minSegSize); - 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; - Facet getFirstFacetOnTheRayFromCamToThePoint(int cam, const Point3d& p, Point3d& intersectPt) const; - /** * @brief Function that returns the next geometry intersected by the ray. * The function handles different cases, whether we come from an edge, a facet or a vertex. From 5f70f8f4f316f00b4ff37e387642d09762ee4fc7 Mon Sep 17 00:00:00 2001 From: Enguerrand DE SMET Date: Thu, 30 Jul 2020 16:45:29 +0200 Subject: [PATCH 08/26] [fuseCut] clean --- src/aliceVision/fuseCut/DelaunayGraphCut.cpp | 7 +++---- src/aliceVision/fuseCut/DelaunayGraphCut.hpp | 6 +++--- src/software/pipeline/main_meshing.cpp | 2 +- 3 files changed, 7 insertions(+), 8 deletions(-) diff --git a/src/aliceVision/fuseCut/DelaunayGraphCut.cpp b/src/aliceVision/fuseCut/DelaunayGraphCut.cpp index 365faa380e..5b48ab1bab 100644 --- a/src/aliceVision/fuseCut/DelaunayGraphCut.cpp +++ b/src/aliceVision/fuseCut/DelaunayGraphCut.cpp @@ -24,7 +24,9 @@ #include #include + #include +#include #include #include @@ -493,8 +495,6 @@ void DelaunayGraphCut::computeDelaunay() void DelaunayGraphCut::initCells() { - ALICEVISION_LOG_DEBUG("initCells ...\n"); - _cellsAttr.resize(_tetrahedralization->nb_cells()); // or nb_finite_cells() if keeps_infinite() ALICEVISION_LOG_INFO(_cellsAttr.size() << " cells created by tetrahedralization."); @@ -513,7 +513,7 @@ void DelaunayGraphCut::initCells() } } - ALICEVISION_LOG_DEBUG("initCells done\n"); + ALICEVISION_LOG_DEBUG("initCells [" << _tetrahedralization->nb_cells() << "] done"); } void DelaunayGraphCut::displayStatistics() @@ -1635,7 +1635,6 @@ void DelaunayGraphCut::fillGraph(bool fixesSigma, float nPixelSizeBehind, nAvCams += 1; } } - 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)); diff --git a/src/aliceVision/fuseCut/DelaunayGraphCut.hpp b/src/aliceVision/fuseCut/DelaunayGraphCut.hpp index 4fd81677a2..e21f493caf 100644 --- a/src/aliceVision/fuseCut/DelaunayGraphCut.hpp +++ b/src/aliceVision/fuseCut/DelaunayGraphCut.hpp @@ -264,11 +264,11 @@ class DelaunayGraphCut std::map> neighboringCellsPerVertexTmp; int coutInvalidVertices = 0; - for(CellIndex ci = 0, nbCells = _tetrahedralization->nb_cells(); ci < nbCells; ++ci) + for (CellIndex ci = 0; ci < _tetrahedralization->nb_cells(); ++ci) { for(VertexIndex k = 0; k < 4; ++k) { - CellIndex vi = _tetrahedralization->cell_vertex(ci, k); + const VertexIndex vi = _tetrahedralization->cell_vertex(ci, k); if(vi == GEO::NO_VERTEX || vi >= _verticesCoords.size()) { ++coutInvalidVertices; @@ -298,7 +298,7 @@ class DelaunayGraphCut * @param lvi * @return global index of the lvi'th neighboring cell */ - CellIndex vertexToCells(VertexIndex vi, int lvi) const + inline CellIndex vertexToCells(VertexIndex vi, int lvi) const { const std::vector& localCells = _neighboringCellsPerVertex.at(vi); if(lvi >= localCells.size()) diff --git a/src/software/pipeline/main_meshing.cpp b/src/software/pipeline/main_meshing.cpp index 03fca301bd..95a8288471 100644 --- a/src/software/pipeline/main_meshing.cpp +++ b/src/software/pipeline/main_meshing.cpp @@ -455,7 +455,7 @@ int aliceVision_main(int argc, char* argv[]) if (boundingBox.isInitialized()) boundingBox.toHexahedron(&hexah[0]); - else if(meshingFromDepthMaps && !estimateSpaceFromSfM) + else if(meshingFromDepthMaps && (!estimateSpaceFromSfM || sfmData.getLandmarks().empty())) fs.divideSpaceFromDepthMaps(&hexah[0], minPixSize); else fs.divideSpaceFromSfM(sfmData, &hexah[0], estimateSpaceMinObservations, estimateSpaceMinObservationAngle); From 4fc538388b180eb3b73f131da1409696dc063383 Mon Sep 17 00:00:00 2001 From: Enguerrand DE SMET Date: Thu, 30 Jul 2020 16:50:07 +0200 Subject: [PATCH 09/26] [FuseCut] DelaunayGraphCut: Add geometries history - add a history structure for debugging purposes --- src/aliceVision/fuseCut/DelaunayGraphCut.cpp | 56 ++++++++++++++++++++ src/aliceVision/fuseCut/DelaunayGraphCut.hpp | 18 +++++++ 2 files changed, 74 insertions(+) diff --git a/src/aliceVision/fuseCut/DelaunayGraphCut.cpp b/src/aliceVision/fuseCut/DelaunayGraphCut.cpp index 5b48ab1bab..2b2ed18a3a 100644 --- a/src/aliceVision/fuseCut/DelaunayGraphCut.cpp +++ b/src/aliceVision/fuseCut/DelaunayGraphCut.cpp @@ -4,6 +4,9 @@ // v. 2.0. If a copy of the MPL was not distributed with this file, // You can obtain one at https://mozilla.org/MPL/2.0/. +// To analyse history of intersected geometries during vote +// #define ALICEVISION_DEBUG_VOTE + #include "DelaunayGraphCut.hpp" // #include #include @@ -395,6 +398,17 @@ void createVerticesWithVisibilities(const StaticVector& cams, std::vector

CArr[cam] - originPt).normalize() * (dir == EDirection::toTheCam ? 1.0 : -1.0); +#ifdef ALICEVISION_DEBUG_VOTE + IntersectionHistory history(mp->CArr[cam], originPt, dirVect); +#endif out_nstepsFront = 0; Facet lastIntersectedFacet; @@ -1678,12 +1695,18 @@ void DelaunayGraphCut::fillGraphPartPtRc(int& out_nstepsFront, int& out_nstepsBe const GeometryIntersection previousGeometry = geometry; const Point3d lastIntersectPt = intersectPt; +#ifdef ALICEVISION_DEBUG_VOTE + history.append(geometry, intersectPt); +#endif ++out_nstepsFront; geometry = intersectNextGeom(previousGeometry, originPt, dirVect, intersectPt, 1.0e-3, lastIntersectPt); if (geometry.type == EGeometryType::None) { +#ifdef ALICEVISION_DEBUG_VOTE + // exportBackPropagationMesh("fillGraph_ToCam_typeNone", history.geometries, originPt, mp->CArr[cam]); +#endif // throw std::runtime_error("[Error]: fillGraphPartPtRc toTheCam, cause: geometry cannot be found."); ALICEVISION_LOG_TRACE("[Error]: fillGraphPartPtRc toTheCam, cause: geometry cannot be found."); break; @@ -1707,6 +1730,9 @@ void DelaunayGraphCut::fillGraphPartPtRc(int& out_nstepsFront, int& out_nstepsBe const Facet mFacet = mirrorFacet(geometry.facet); if (isInvalidOrInfiniteCell(mFacet.cellIndex)) { +#ifdef ALICEVISION_DEBUG_VOTE + // exportBackPropagationMesh("fillGraph_ToCam_invalidMirorFacet", history.geometries, originPt, mp->CArr[cam]); +#endif // throw std::runtime_error("[Error]: fillGraphPartPtRc toTheCam, cause: invalidOrInfinite miror facet."); ALICEVISION_LOG_TRACE("[Error]: fillGraphPartPtRc toTheCam, cause: invalidOrInfinite miror facet."); break; @@ -1735,6 +1761,9 @@ void DelaunayGraphCut::fillGraphPartPtRc(int& out_nstepsFront, int& out_nstepsBe const EDirection dir = EDirection::behindThePoint; const Point3d dirVect = (mp->CArr[cam] - originPt).normalize() * (dir == EDirection::toTheCam ? 1 : -1); +#ifdef ALICEVISION_DEBUG_VOTE + IntersectionHistory history(mp->CArr[cam], originPt, dirVect); +#endif out_nstepsBehind = 0; bool firstIteration = true; @@ -1746,6 +1775,9 @@ void DelaunayGraphCut::fillGraphPartPtRc(int& out_nstepsFront, int& out_nstepsBe const GeometryIntersection previousGeometry = geometry; const Point3d lastIntersectPt = intersectPt; +#ifdef ALICEVISION_DEBUG_VOTE + history.append(geometry, intersectPt); +#endif ++out_nstepsBehind; // Vote for the first facet found (only once) @@ -1762,6 +1794,9 @@ void DelaunayGraphCut::fillGraphPartPtRc(int& out_nstepsFront, int& out_nstepsBe // If we come from a facet, the next intersection must exist (even if the mirror facet is invalid, which is verified after taking mirror facet) if(previousGeometry.type == EGeometryType::Facet) { +#ifdef ALICEVISION_DEBUG_VOTE + // exportBackPropagationMesh("fillGraph_behindThePoint_NoneButPreviousIsFacet", history.geometries, originPt, mp->CArr[cam]); +#endif // throw std::runtime_error("[DelaunayGraphCut::fillGraphPartPtRc] first loop: Geometry cannot be found."); ALICEVISION_LOG_TRACE("[Error]: fillGraphPartPtRc behindThePoint, cause: None geometry but previous is Facet."); } @@ -1879,6 +1914,9 @@ void DelaunayGraphCut::forceTedgesByGradientIJCV(bool fixesSigma, float nPixelSi const EDirection dir = EDirection::toTheCam; const Point3d dirVect = (mp->CArr[cam] - originPt).normalize() * (dir == EDirection::toTheCam ? 1.0 : -1.0); +#ifdef ALICEVISION_DEBUG_VOTE + IntersectionHistory history(mp->CArr[cam], originPt, dirVect); +#endif int nstepsFront = 0; // As long as we find a next geometry @@ -1888,12 +1926,18 @@ void DelaunayGraphCut::forceTedgesByGradientIJCV(bool fixesSigma, float nPixelSi const GeometryIntersection previousGeometry = geometry; const Point3d lastIntersectPt = intersectPt; +#ifdef ALICEVISION_DEBUG_VOTE + history.append(geometry, intersectPt); +#endif ++nstepsFront; geometry = intersectNextGeom(previousGeometry, originPt, dirVect, intersectPt, 1.0e-3, lastIntersectPt); if (geometry.type == EGeometryType::None) { +#ifdef ALICEVISION_DEBUG_VOTE + // exportBackPropagationMesh("forceTedgesByGradientIJCV_ToCam_typeNone", history.geometries, originPt, mp->CArr[cam]); +#endif // throw std::runtime_error("[DelaunayGraphCut::fillGraphPartPtRc] first loop: Geometry cannot be found."); ALICEVISION_LOG_TRACE("[Error]: forceTedgesByGradientIJCV toTheCam, cause: geometry cannot be found."); break; @@ -1917,6 +1961,9 @@ void DelaunayGraphCut::forceTedgesByGradientIJCV(bool fixesSigma, float nPixelSi const Facet mFacet = mirrorFacet(geometry.facet); if (isInvalidOrInfiniteCell(mFacet.cellIndex)) { +#ifdef ALICEVISION_DEBUG_VOTE + // exportBackPropagationMesh("forceTedgesByGradientIJCV_ToCam_invalidMirorFacet", history.geometries, originPt, mp->CArr[cam]); +#endif // throw std::runtime_error("[DelaunayGraphCut::fillGraphPartPtRc] first loop: Geometry cannot be found."); ALICEVISION_LOG_TRACE("[Error]: forceTedgesByGradientIJCV toTheCam, cause: invalidOrInfinite miror facet."); break; @@ -1938,6 +1985,9 @@ void DelaunayGraphCut::forceTedgesByGradientIJCV(bool fixesSigma, float nPixelSi const EDirection dir = EDirection::behindThePoint; const Point3d dirVect = (mp->CArr[cam] - originPt).normalize() * (dir == EDirection::toTheCam ? 1 : -1); +#ifdef ALICEVISION_DEBUG_VOTE + IntersectionHistory history(mp->CArr[cam], originPt, dirVect); +#endif int nstepsBehind = 0; Facet lastIntersectedFacet; @@ -1948,6 +1998,9 @@ void DelaunayGraphCut::forceTedgesByGradientIJCV(bool fixesSigma, float nPixelSi const GeometryIntersection previousGeometry = geometry; const Point3d lastIntersectPt = intersectPt; +#ifdef ALICEVISION_DEBUG_VOTE + history.append(geometry, intersectPt); +#endif ++nstepsBehind; geometry = intersectNextGeom(previousGeometry, originPt, dirVect, intersectPt, 1.0e-3, lastIntersectPt); @@ -1958,6 +2011,9 @@ void DelaunayGraphCut::forceTedgesByGradientIJCV(bool fixesSigma, float nPixelSi // If we come from a facet, the next intersection must exist (even if the mirror facet is invalid, which is verified later) if (previousGeometry.type == EGeometryType::Facet) { +#ifdef ALICEVISION_DEBUG_VOTE + // exportBackPropagationMesh("forceTedgesByGradientIJCV_behindThePoint_typeNone", history.geometries, originPt, mp->CArr[cam]); +#endif // throw std::runtime_error("[Error]: forceTedgesByGradientIJCV behindThePoint, cause: geometry cannot be found."); ALICEVISION_LOG_TRACE("[Error]: forceTedgesByGradientIJCV behindThePoint, cause: geometry cannot be found."); } diff --git a/src/aliceVision/fuseCut/DelaunayGraphCut.hpp b/src/aliceVision/fuseCut/DelaunayGraphCut.hpp index e21f493caf..d3bfb9d1a9 100644 --- a/src/aliceVision/fuseCut/DelaunayGraphCut.hpp +++ b/src/aliceVision/fuseCut/DelaunayGraphCut.hpp @@ -127,6 +127,24 @@ class DelaunayGraphCut behindThePoint }; + struct IntersectionHistory + { + size_t steps = 0; + Point3d cam; + Point3d originPt; + Point3d dirVect; + + std::vector geometries; + std::vector intersectPts; + std::vector vecToCam; + std::vector distToCam; + std::vector angleToCam; + IntersectionHistory(const Point3d& c, const Point3d& oPt, const Point3d& diV) : cam{c}, originPt{oPt}, dirVect{diV} {} + + void append(const GeometryIntersection& geom, const Point3d& intersectPt); + }; + + mvsUtils::MultiViewParams* mp; GEO::Delaunay_var _tetrahedralization; From 08123bac99fa56953e9d8e025af0da879c76b8ef Mon Sep 17 00:00:00 2001 From: Enguerrand DE SMET Date: Thu, 30 Jul 2020 17:00:57 +0200 Subject: [PATCH 10/26] [FuseCut] DelaunayGraphCut: update mesh export for debugging - small clean and new optional variables for mesh export - add new functions to export meshes for debugging --- src/aliceVision/fuseCut/DelaunayGraphCut.cpp | 170 ++++++++++++++----- src/aliceVision/fuseCut/DelaunayGraphCut.hpp | 5 +- 2 files changed, 134 insertions(+), 41 deletions(-) diff --git a/src/aliceVision/fuseCut/DelaunayGraphCut.cpp b/src/aliceVision/fuseCut/DelaunayGraphCut.cpp index 2b2ed18a3a..206efc3368 100644 --- a/src/aliceVision/fuseCut/DelaunayGraphCut.cpp +++ b/src/aliceVision/fuseCut/DelaunayGraphCut.cpp @@ -2705,11 +2705,11 @@ mesh::Mesh* DelaunayGraphCut::createMesh(bool filterHelperPointsTriangles) return me; } -mesh::Mesh* DelaunayGraphCut::createTetrahedralMesh() const +mesh::Mesh* DelaunayGraphCut::createTetrahedralMesh(bool filter, const float& downscaleFactor, const std::function getScore) const { ALICEVISION_LOG_INFO("Create mesh of the tetrahedralization."); - ALICEVISION_LOG_INFO("# vertixes: " << _verticesCoords.size()); + ALICEVISION_LOG_INFO("# vertices: " << _verticesCoords.size()); if(_cellsAttr.empty()) return nullptr; @@ -2742,7 +2742,7 @@ mesh::Mesh* DelaunayGraphCut::createTetrahedralMesh() const displayAcc("acc_camSize", acc_camSize); } - float max_score = 1.0; + float maxScore = 1.0; { Accumulator acc_cellScore; Accumulator acc_cellSWeight; @@ -2751,10 +2751,10 @@ mesh::Mesh* DelaunayGraphCut::createTetrahedralMesh() const Accumulator acc_fullnessScore; Accumulator acc_emptinessScore; Accumulator acc_on; + Accumulator acc_selectedScore; - for(CellIndex ci = 0; ci < _cellsAttr.size(); ++ci) + for(const GC_cellInfo& cellAttr : _cellsAttr) { - const GC_cellInfo& cellAttr = _cellsAttr[ci]; acc_cellScore(cellAttr.cellSWeight - cellAttr.cellTWeight); acc_cellSWeight(cellAttr.cellSWeight); acc_cellTWeight(cellAttr.cellTWeight); @@ -2767,6 +2767,8 @@ mesh::Mesh* DelaunayGraphCut::createTetrahedralMesh() const acc_fullnessScore(cellAttr.fullnessScore); acc_emptinessScore(cellAttr.emptinessScore); acc_on(cellAttr.on); + + acc_selectedScore(getScore(cellAttr)); } displayAcc("cellScore", acc_cellScore); @@ -2776,9 +2778,17 @@ mesh::Mesh* DelaunayGraphCut::createTetrahedralMesh() const displayAcc("fullnessScore", acc_fullnessScore); displayAcc("emptinessScore", acc_emptinessScore); displayAcc("on", acc_on); - max_score = 4.0f * extract::mean(acc_emptinessScore); + + displayAcc("selected", acc_selectedScore); + maxScore = 4.0f * extract::mean(acc_selectedScore); } - ALICEVISION_LOG_DEBUG("createTetrahedralMesh: max_score: " << max_score); + + + ALICEVISION_LOG_DEBUG("createTetrahedralMesh: maxScore: " << maxScore); + + // Prevent division by zero (getRGBFromJetColorMap) + if (maxScore == 0) + maxScore = 1; // loop over all facets for(CellIndex ci = 0; ci < _cellsAttr.size(); ++ci) @@ -2791,13 +2801,13 @@ mesh::Mesh* DelaunayGraphCut::createTetrahedralMesh() const { const VertexIndex vi = _tetrahedralization->cell_vertex(ci, k); const GC_vertexInfo& vertexAttr = _verticesAttr[vi]; - if(vertexAttr.cams.size() <= 3) + if(filter && vertexAttr.cams.size() <= 3) { ++weakVertex; } const Facet f1(ci, k); const Facet f2 = mirrorFacet(f1); - if(isInvalidOrInfiniteCell(f2.cellIndex)) + if(filter && isInvalidOrInfiniteCell(f2.cellIndex)) { invalid = true; continue; @@ -2805,12 +2815,12 @@ mesh::Mesh* DelaunayGraphCut::createTetrahedralMesh() const pointscellCenter = pointscellCenter + _verticesCoords[vi]; } pointscellCenter = pointscellCenter / 4.0; - if(invalid) + if(filter && invalid) { - // If one of the facet is invalid, discard the tetrahedron. + // If one of the mirror facet is invalid, discard the tetrahedron. continue; } - if (weakVertex >= 3) + if (filter && (weakVertex >= 3)) { // If the tetrahedron mainly composed of weak vertices, skip it. // So keep tetrahedra that are linked to at least 3 good vertices. @@ -2819,8 +2829,8 @@ mesh::Mesh* DelaunayGraphCut::createTetrahedralMesh() const } const GC_cellInfo& cellAttr = _cellsAttr[ci]; - const float score = cellAttr.emptinessScore; // cellAttr.cellSWeight - cellAttr.cellTWeight; - if(score < (max_score / 1000.0f)) + const float score = getScore(cellAttr); // cellAttr.cellSWeight - cellAttr.cellTWeight; + if(filter && (score < (maxScore / 1000.0f))) { // Skip too low score cells continue; @@ -2840,44 +2850,43 @@ mesh::Mesh* DelaunayGraphCut::createTetrahedralMesh() const { points[k] = _verticesCoords[vertices[k]]; // Downscale cell for visibility - points[k] = pointscellCenter + ((points[k] - pointscellCenter) * 0.95); + points[k] = pointscellCenter + ((points[k] - pointscellCenter) * downscaleFactor); } const Facet f2 = mirrorFacet(f1); + bool clockwise = false; + //// 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(); + if (!isInvalidOrInfiniteCell(f2.cellIndex)) + { + const Point3d D1 = _verticesCoords[getOppositeVertexIndex(f1)]; + const Point3d D2 = _verticesCoords[getOppositeVertexIndex(f2)]; - const double dd1 = orientedPointPlaneDistance(D1, points[0], N); - const double dd2 = orientedPointPlaneDistance(D2, points[0], N); + const Point3d N = cross((points[1] - points[0]).normalize(), (points[2] - points[0]).normalize()).normalize(); - bool clockwise = false; - if(dd1 == 0.0) - { - if(dd2 == 0.0) + const double dd1 = orientedPointPlaneDistance(D1, points[0], N); + const double dd2 = orientedPointPlaneDistance(D2, points[0], N); + if (dd1 == 0.0) { - ALICEVISION_LOG_WARNING("createMesh: bad triangle orientation."); - } - if(dd2 > 0.0) - { - clockwise = true; + if (dd2 == 0.0) + { + ALICEVISION_LOG_WARNING("createMesh: bad triangle orientation."); + } + if (dd2 > 0.0) + { + clockwise = true; + } } - } - else - { - if(dd1 < 0.0) + else { - clockwise = true; + if (dd1 < 0.0) + { + clockwise = true; + } } } - const rgb color = getRGBFromJetColorMap(score / max_score); + const rgb color = getRGBFromJetColorMap(score / maxScore); const std::size_t vertexBaseIndex = me->pts.size(); for(const Point3d& p: points) { @@ -2910,6 +2919,87 @@ mesh::Mesh* DelaunayGraphCut::createTetrahedralMesh() const return me; } +void DelaunayGraphCut::exportDebugMesh(const std::string& filename, const Point3d& fromPt, const Point3d& toPt) +{ + std::unique_ptr mesh(createTetrahedralMesh(false, 0.999f)); + std::unique_ptr meshf(createTetrahedralMesh(true, 0.999f)); + + // hack add direction vector from fromPt to toPt + { + const Point3d dirVec = fromPt - toPt; + const Point3d deltaVec = (dirVec * 0.1f); + + mesh->pts.push_back(toPt - deltaVec); + mesh->colors().push_back(rgb(255, 0, 0)); + mesh->pts.push_back(fromPt + deltaVec); + mesh->colors().push_back(rgb(0, 0, 255)); + mesh->pts.push_back(fromPt + deltaVec * 0.001f); + mesh->colors().push_back(rgb(0, 0, 255)); + + meshf->pts.push_back(toPt - deltaVec); + meshf->colors().push_back(rgb(255, 0, 0)); + meshf->pts.push_back(fromPt + deltaVec); + meshf->colors().push_back(rgb(0, 0, 255)); + meshf->pts.push_back(fromPt + deltaVec * 0.001f); + meshf->colors().push_back(rgb(0, 0, 255)); + } + { + mesh::Mesh::triangle t; + t.alive = true; + t.v[0] = meshf->pts.size() - 3; + t.v[1] = meshf->pts.size() - 2; + t.v[2] = meshf->pts.size() - 1; + meshf->tris.push_back(t); + } + { + mesh::Mesh::triangle t; + t.alive = true; + t.v[0] = mesh->pts.size() - 3; + t.v[1] = mesh->pts.size() - 2; + t.v[2] = mesh->pts.size() - 1; + mesh->tris.push_back(t); + } + + const std::string tempDirPath = boost::filesystem::temp_directory_path().generic_string(); + mesh->saveToObj(tempDirPath + "/" + filename + ".obj"); + meshf->saveToObj(tempDirPath + "/" + filename + "Filtered.obj"); +} + +void DelaunayGraphCut::exportFullScoreMeshs() +{ + std::unique_ptr meshEmptiness(createTetrahedralMesh(false, 0.999f, [](const GC_cellInfo& c) { return c.emptinessScore; })); + std::unique_ptr meshFullness(createTetrahedralMesh(false, 0.999f, [](const GC_cellInfo& c) { return c.fullnessScore; })); + std::unique_ptr meshSWeight(createTetrahedralMesh(false, 0.999f, [](const GC_cellInfo& c) { return c.cellSWeight; })); + std::unique_ptr meshTWeight(createTetrahedralMesh(false, 0.999f, [](const GC_cellInfo& c) { return c.cellTWeight; })); + std::unique_ptr meshOn(createTetrahedralMesh(false, 0.999f, [](const GC_cellInfo& c) { return c.on; })); + + const std::string tempDirPath = boost::filesystem::temp_directory_path().generic_string(); + + meshEmptiness->saveToObj(tempDirPath + "/meshEmptiness.obj"); + meshFullness->saveToObj(tempDirPath + "/meshFullness.obj"); + meshSWeight->saveToObj(tempDirPath + "/meshSWeight.obj"); + meshTWeight->saveToObj(tempDirPath + "/meshTWeight.obj"); + meshOn->saveToObj(tempDirPath + "/meshOn.obj"); +} + +void DelaunayGraphCut::exportBackPropagationMesh(const std::string& filename, std::vector& intersectedGeom, const Point3d& fromPt, const Point3d& toPt) +{ + // Clean _cellsAttr emptinessScore + for (int i = 0; i < _cellsAttr.size(); ++i) + _cellsAttr[i].emptinessScore = 0.0f; + + // Vote only for listed intersected geom facets + for (size_t i = 0; i < intersectedGeom.size(); i++) + { + const GeometryIntersection& geo = intersectedGeom[i]; + + if (geo.type == EGeometryType::Facet) + _cellsAttr[geo.facet.cellIndex].emptinessScore += i; + } + + exportDebugMesh(filename + "_BackProp", fromPt, toPt); +} + void DelaunayGraphCut::segmentFullOrFree(bool full, StaticVector** 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 d3bfb9d1a9..3f42142f0a 100644 --- a/src/aliceVision/fuseCut/DelaunayGraphCut.hpp +++ b/src/aliceVision/fuseCut/DelaunayGraphCut.hpp @@ -449,7 +449,10 @@ class DelaunayGraphCut void leaveLargestFullSegmentOnly(); mesh::Mesh* createMesh(bool filterHelperPointsTriangles = true); - mesh::Mesh* createTetrahedralMesh() const; + mesh::Mesh* createTetrahedralMesh(bool filter = true, const float& downscaleFactor = 0.95f, const std::function getScore = [](const GC_cellInfo& c) { return c.emptinessScore; }) const; + void exportDebugMesh(const std::string& filename, const Point3d& fromPt, const Point3d& toPt); + void exportFullScoreMeshs(); + void exportBackPropagationMesh(const std::string& filename, std::vector& intersectedGeom, const Point3d& fromPt, const Point3d& toPt); }; } // namespace fuseCut From 3c244267afa0dd8e52ff8b338f9020733a42ee8b Mon Sep 17 00:00:00 2001 From: Enguerrand DE SMET Date: Fri, 31 Jul 2020 16:37:55 +0200 Subject: [PATCH 11/26] [fuseCut] DelaunayGraphCut: add unit test --- src/aliceVision/fuseCut/CMakeLists.txt | 9 + .../fuseCut/DelaunayGraphCut_test.cpp | 174 ++++++++++++++++++ 2 files changed, 183 insertions(+) create mode 100644 src/aliceVision/fuseCut/DelaunayGraphCut_test.cpp diff --git a/src/aliceVision/fuseCut/CMakeLists.txt b/src/aliceVision/fuseCut/CMakeLists.txt index f222fb2d5f..100541ed5b 100644 --- a/src/aliceVision/fuseCut/CMakeLists.txt +++ b/src/aliceVision/fuseCut/CMakeLists.txt @@ -38,3 +38,12 @@ alicevision_add_library(aliceVision_fuseCut nanoflann Boost::boost ) + +# Unit tests +alicevision_add_test(DelaunayGraphCut_test.cpp + NAME "fuseCut_delaunayGraphCut" + LINKS aliceVision_fuseCut + aliceVision_sfm + aliceVision_multiview + aliceVision_multiview_test_data +) \ No newline at end of file diff --git a/src/aliceVision/fuseCut/DelaunayGraphCut_test.cpp b/src/aliceVision/fuseCut/DelaunayGraphCut_test.cpp new file mode 100644 index 0000000000..992fab08b8 --- /dev/null +++ b/src/aliceVision/fuseCut/DelaunayGraphCut_test.cpp @@ -0,0 +1,174 @@ +// This file is part of the AliceVision project. +// Copyright (c) 2020 AliceVision contributors. +// This Source Code Form is subject to the terms of the Mozilla Public License, +// v. 2.0. If a copy of the MPL was not distributed with this file, +// You can obtain one at https://mozilla.org/MPL/2.0/. + +#include +#include +#include + +#include + +#include + +#define BOOST_TEST_MODULE fuseCut + +#include +#include + +using namespace aliceVision; +using namespace aliceVision::fuseCut; +using namespace aliceVision::sfmData; + +SfMData generateSfm(const NViewDatasetConfigurator& config, const size_t size = 3, camera::EINTRINSIC eintrinsic = camera::EINTRINSIC::PINHOLE_CAMERA_RADIAL3); + +BOOST_AUTO_TEST_CASE(fuseCut_delaunayGraphCut) +{ + system::Logger::get()->setLogLevel(system::EVerboseLevel::Trace); + + const NViewDatasetConfigurator config(1000, 1000, 500, 500, 1, 0); + SfMData sfmData = generateSfm(config, 6); + + mvsUtils::MultiViewParams mp(sfmData, "", "", "", false); + + mp.userParams.put("LargeScale.universePercentile", 0.999); + mp.userParams.put("delaunaycut.forceTEdgeDelta", 0.1f); + mp.userParams.put("delaunaycut.seed", 1); + + std::array hexah; + + Fuser fs(&mp); + const size_t minObservations = 2; + const float minObservationsAngle = 0.01f; + fs.divideSpaceFromSfM(sfmData, &hexah[0], minObservations, minObservationsAngle); + + StaticVector cams; + cams.resize(mp.getNbCameras()); + for (int i = 0; i < cams.size(); ++i) + cams[i] = i; + + const std::string tempDirPath = boost::filesystem::temp_directory_path().generic_string(); + + DelaunayGraphCut delaunayGC(&mp); + ALICEVISION_LOG_TRACE("Creating dense point cloud witout support pts."); + + // delaunayGC.createDensePointCloud(&hexah[0], cams, &sfmData, nullptr); + const float minDist = (hexah[0] - hexah[1]).size() / 1000.0f; + // add points for cam centers + delaunayGC.addPointsFromCameraCenters(cams, minDist); + // add points from sfm + delaunayGC.addPointsFromSfM(&hexah[0], cams, sfmData); + + ALICEVISION_LOG_TRACE("Generated pts:"); + for (size_t i = 0; i < delaunayGC._verticesCoords.size(); i++) + ALICEVISION_LOG_TRACE("[" << i << "]: " << delaunayGC._verticesCoords[i].x << ", " << delaunayGC._verticesCoords[i].y << ", " << delaunayGC._verticesCoords[i].z); + + // delaunayGC.createGraphCut(&hexah[0], cams, tempDirPath + "/", tempDirPath + "/SpaceCamsTracks/", false); + delaunayGC.initVertices(); + delaunayGC.computeDelaunay(); + delaunayGC.displayStatistics(); + delaunayGC.computeVerticesSegSize(true, 0.0f); + delaunayGC.voteFullEmptyScore(cams, tempDirPath); + delaunayGC.reconstructGC(&hexah[0]); + + ALICEVISION_LOG_TRACE("CreateGraphCut Done."); +} + +/** + * @brief Generate syntesize dataset with succesion of n(size) alignaed regular thetraedron and two camera on the last thetrahedron. + * + * @param size + * @param eintrinsic + * @return + */ +SfMData generateSfm(const NViewDatasetConfigurator& config, const size_t size, camera::EINTRINSIC eintrinsic) +{ + assert(size > 0); + + SfMData sfm_data; + + const double heightT = std::sqrt(0.75); + const double centerMinT = 0.5 / (2.0 * heightT); + + // Y axis UP + std::vector camsPts; + camsPts.push_back({ static_cast(size) + 10.0, 0.0, 0.0 }); + camsPts.push_back({ static_cast(size) - 0.5 + 10.0, 0.0, heightT - 0.0}); + + const size_t ptsSize = 5 * size - 3; + // Construct our pts (Y axis UP) + Mat3X matPts(3, ptsSize); + + // Bottom down + for (size_t i = 0; i < size; i++) + matPts.col(i) << i, 0.0, 0.0; + + // Bottom up + for (size_t i = 0; i < size - 1; i++) + matPts.col(size + i) << (i + 0.5), 0.0, heightT; + + // Top + for (size_t i = 0; i < 2 * size - 1; i++) + matPts.col(2 * size - 1 + i) << (i / 2.0 + 0.5), heightT, (i % 2 == 0 ? centerMinT : heightT - centerMinT); + + // Center + for (size_t i = 0; i < size - 1; i ++) + matPts.col(4 * size - 2 + i) << (i + (1 + (i % 2)) * 0.5), centerMinT + (i % 2) * 0.1, (i % 2 == 0 ? centerMinT : heightT - centerMinT); // center + + // jitters all pts + matPts += Mat3X::Random(3, ptsSize) * 0.05; + + // Same internals parameters for all cams + Mat3 internalParams; + internalParams << config._fx, 0, config._cx, + 0, config._fy, config._cy, + 0, 0, 1; + + // 1. Views + for (int i = 0; i < camsPts.size(); ++i) + { + const IndexT id_view = i, id_pose = i, id_intrinsic = 0; // shared intrinsics + sfm_data.views[i] = std::make_shared("", id_view, id_intrinsic, id_pose, config._cx * 2, config._cy * 2); + } + + // 2. Poses + // and compute projected pts for cam + std::vector projectedPtsPerCam; + for (size_t i = 0; i < camsPts.size(); i++) + { + Mat34 P; + // All cam look to the first thetraherdron center + const Vec3 camCenter = camsPts[i]; + const Mat3 rotationMat = LookAt(matPts.col(0) - camCenter); + const Vec3 translation = -rotationMat * camCenter; + P_from_KRt(internalParams, rotationMat, translation, &P); + projectedPtsPerCam.push_back(project(P, matPts)); + + geometry::Pose3 pose(rotationMat, camCenter); + sfm_data.setPose(*sfm_data.views.at(i), CameraPose(pose)); + } + + // 3. Intrinsic data (shared, so only one camera intrinsic is defined) + { + const unsigned int w = config._cx * 2; + const unsigned int h = config._cy * 2; + sfm_data.intrinsics[0] = createIntrinsic(eintrinsic, w, h, config._fx, config._cx, config._cy); + } + + // 4. Landmarks + const double unknownScale = 0.0; + for (size_t i = 0; i < ptsSize; ++i) + { + Landmark landmark; + landmark.X = matPts.col(i); + for (int j = 0; j < camsPts.size(); ++j) + { + const Vec2 pt = projectedPtsPerCam[j].col(i); + landmark.observations[j] = Observation(pt, i, unknownScale); + } + sfm_data.structure[i] = landmark; + } + + return sfm_data; +} \ No newline at end of file From e3f77e330521a91e1c4584d82fcc580975682430 Mon Sep 17 00:00:00 2001 From: Enguerrand DE SMET Date: Tue, 4 Aug 2020 10:39:58 +0200 Subject: [PATCH 12/26] [fuseCut] DelaunayGraphCut: Use an adaptive epsilon factor for the intersection of geometries --- src/aliceVision/fuseCut/DelaunayGraphCut.cpp | 11 ++++++----- src/aliceVision/fuseCut/DelaunayGraphCut.hpp | 10 ++++++---- 2 files changed, 12 insertions(+), 9 deletions(-) diff --git a/src/aliceVision/fuseCut/DelaunayGraphCut.cpp b/src/aliceVision/fuseCut/DelaunayGraphCut.cpp index 206efc3368..00f14d13a6 100644 --- a/src/aliceVision/fuseCut/DelaunayGraphCut.cpp +++ b/src/aliceVision/fuseCut/DelaunayGraphCut.cpp @@ -1289,7 +1289,7 @@ DelaunayGraphCut::GeometryIntersection DelaunayGraphCut::intersectNextGeom(const DelaunayGraphCut::GeometryIntersection& inGeometry, const Point3d& originPt, const Point3d& dirVect, Point3d& intersectPt, - const float epsilon, const Point3d& lastIntersectPt) const + const double epsilonFactor, const Point3d& lastIntersectPt) const { switch (inGeometry.type) @@ -1307,7 +1307,7 @@ DelaunayGraphCut::intersectNextGeom(const DelaunayGraphCut::GeometryIntersection const Facet intersectionFacet(tetrahedronIndex, i); - const GeometryIntersection result = rayIntersectTriangle(originPt, dirVect, intersectionFacet, intersectPt, epsilon, &lastIntersectPt); + const GeometryIntersection result = rayIntersectTriangle(originPt, dirVect, intersectionFacet, intersectPt, epsilonFactor, &lastIntersectPt); if (result.type != EGeometryType::None) return result; } @@ -1327,7 +1327,7 @@ DelaunayGraphCut::intersectNextGeom(const DelaunayGraphCut::GeometryIntersection // Define the facet to intersect const Facet facet(adjCellIndex, localVertexIndex); - const GeometryIntersection result = rayIntersectTriangle(originPt, dirVect, facet, intersectPt, epsilon, &lastIntersectPt); + const GeometryIntersection result = rayIntersectTriangle(originPt, dirVect, facet, intersectPt, epsilonFactor, &lastIntersectPt); if (result.type != EGeometryType::None) return result; } @@ -1351,7 +1351,7 @@ DelaunayGraphCut::intersectNextGeom(const DelaunayGraphCut::GeometryIntersection for (const Facet& facet : opositeFacets) { - const GeometryIntersection result = rayIntersectTriangle(originPt, dirVect, facet, intersectPt, epsilon, &lastIntersectPt); + const GeometryIntersection result = rayIntersectTriangle(originPt, dirVect, facet, intersectPt, epsilonFactor, &lastIntersectPt); if (result.type != EGeometryType::None) return result; } @@ -1370,7 +1370,7 @@ DelaunayGraphCut::GeometryIntersection DelaunayGraphCut::rayIntersectTriangle(co const Point3d& DirVec, const DelaunayGraphCut::Facet& facet, Point3d& intersectPt, - const float epsilon, const Point3d* lastIntersectPt) const + const double epsilonFactor, const Point3d* lastIntersectPt) const { const VertexIndex AvertexIndex = getVertexIndex(facet, 0); const VertexIndex BvertexIndex = getVertexIndex(facet, 1); @@ -1380,6 +1380,7 @@ DelaunayGraphCut::GeometryIntersection DelaunayGraphCut::rayIntersectTriangle(co const Point3d* B = &_verticesCoords[BvertexIndex]; const Point3d* C = &_verticesCoords[CvertexIndex]; + const double epsilon = std::min(std::min((*A - *B).size(), (*B - *C).size()), (*A - *C).size()) * epsilonFactor; Point3d tempIntersectPt; const Point2d triangleUv = getLineTriangleIntersectBarycCoords(&tempIntersectPt, A, B, C, &originPt, &DirVec); diff --git a/src/aliceVision/fuseCut/DelaunayGraphCut.hpp b/src/aliceVision/fuseCut/DelaunayGraphCut.hpp index 3f42142f0a..0fff4c3583 100644 --- a/src/aliceVision/fuseCut/DelaunayGraphCut.hpp +++ b/src/aliceVision/fuseCut/DelaunayGraphCut.hpp @@ -380,12 +380,13 @@ class DelaunayGraphCut * @param originPt ray origin point * @param dirVect ray direction * @param intersectPt a reference that will store the computed intersection point for the intersected geometry - * @param epsilon used to define the boundary when we have to consider either a collision with an edge/vertex or a facet + * @param epsilonFactor a multiplicative factor on the smaller side of the facet used to define the boundary when we + * have to consider either a collision with an edge/vertex or a facet. * @param lastIntersectPt constant reference to the last intersection point used to test the direction. * @return */ GeometryIntersection intersectNextGeom(const GeometryIntersection& inGeometry, const Point3d& originPt, - const Point3d& dirVect, Point3d& intersectPt, const float epsilon, const Point3d& lastIntersectPt) const; + const Point3d& dirVect, Point3d& intersectPt, const double epsilonFactor, const Point3d& lastIntersectPt) const; /** * @brief Function that returns the next geometry intersected by the ray on a given facet or None if there are no intersected geometry. @@ -395,12 +396,13 @@ class DelaunayGraphCut * @param DirVec ray direction * @param facet the given facet to intersect with * @param intersectPt a reference that will store the computed intersection point for the next intersecting geometry - * @param epsilon used to define the boundary when we have to consider either a collision with an edge/vertex or a facet + * @param epsilonFactor a multiplicative factor on the smaller side of the facet used to define the boundary when we + * have to consider either a collision with an edge/vertex or a facet. * @param lastIntersectPt pointer to the last intersection point used to test the direction (if not nulllptr) * @return */ GeometryIntersection rayIntersectTriangle(const Point3d& originPt, const Point3d& DirVec, const Facet& facet, - Point3d& intersectPt, const float epsilon, const Point3d* lastIntersectPt = nullptr) const; + Point3d& intersectPt, const double epsilonFactor, const Point3d* lastIntersectPt = nullptr) const; float distFcn(float maxDist, float dist, float distFcnHeight) const; From 0c5fb6051a0a68882fad5d25c165ce63bd3840de Mon Sep 17 00:00:00 2001 From: Enguerrand DE SMET Date: Wed, 5 Aug 2020 11:23:46 +0200 Subject: [PATCH 13/26] [fuseCut] DelaunayGraphCut: Remove unused EDirection struct --- src/aliceVision/fuseCut/DelaunayGraphCut.cpp | 16 ++++++++-------- src/aliceVision/fuseCut/DelaunayGraphCut.hpp | 6 ------ 2 files changed, 8 insertions(+), 14 deletions(-) diff --git a/src/aliceVision/fuseCut/DelaunayGraphCut.cpp b/src/aliceVision/fuseCut/DelaunayGraphCut.cpp index 00f14d13a6..aa98cd8aea 100644 --- a/src/aliceVision/fuseCut/DelaunayGraphCut.cpp +++ b/src/aliceVision/fuseCut/DelaunayGraphCut.cpp @@ -1680,8 +1680,8 @@ void DelaunayGraphCut::fillGraphPartPtRc(int& out_nstepsFront, int& out_nstepsBe // Initialisation GeometryIntersection geometry(vertexIndex); // Starting on global vertex index Point3d intersectPt = originPt; - const EDirection dir = EDirection::toTheCam; - const Point3d dirVect = (mp->CArr[cam] - originPt).normalize() * (dir == EDirection::toTheCam ? 1.0 : -1.0); + // toTheCam + const Point3d dirVect = (mp->CArr[cam] - originPt).normalize(); #ifdef ALICEVISION_DEBUG_VOTE IntersectionHistory history(mp->CArr[cam], originPt, dirVect); @@ -1759,8 +1759,8 @@ void DelaunayGraphCut::fillGraphPartPtRc(int& out_nstepsFront, int& out_nstepsBe // Initialisation GeometryIntersection geometry(vertexIndex); // Starting on global vertex index Point3d intersectPt = originPt; - const EDirection dir = EDirection::behindThePoint; - const Point3d dirVect = (mp->CArr[cam] - originPt).normalize() * (dir == EDirection::toTheCam ? 1 : -1); + // behindThePoint + const Point3d dirVect = (originPt - mp->CArr[cam]).normalize(); #ifdef ALICEVISION_DEBUG_VOTE IntersectionHistory history(mp->CArr[cam], originPt, dirVect); @@ -1912,8 +1912,8 @@ void DelaunayGraphCut::forceTedgesByGradientIJCV(bool fixesSigma, float nPixelSi // Initialisation GeometryIntersection geometry(vertexIndex); // Starting on global vertex index Point3d intersectPt = originPt; - const EDirection dir = EDirection::toTheCam; - const Point3d dirVect = (mp->CArr[cam] - originPt).normalize() * (dir == EDirection::toTheCam ? 1.0 : -1.0); + // toTheCam + const Point3d dirVect = (mp->CArr[cam] - originPt).normalize(); #ifdef ALICEVISION_DEBUG_VOTE IntersectionHistory history(mp->CArr[cam], originPt, dirVect); @@ -1983,8 +1983,8 @@ void DelaunayGraphCut::forceTedgesByGradientIJCV(bool fixesSigma, float nPixelSi // Initialisation GeometryIntersection geometry(vertexIndex); Point3d intersectPt = originPt; - const EDirection dir = EDirection::behindThePoint; - const Point3d dirVect = (mp->CArr[cam] - originPt).normalize() * (dir == EDirection::toTheCam ? 1 : -1); + // behindThePoint + const Point3d dirVect = (originPt - mp->CArr[cam]).normalize(); #ifdef ALICEVISION_DEBUG_VOTE IntersectionHistory history(mp->CArr[cam], originPt, dirVect); diff --git a/src/aliceVision/fuseCut/DelaunayGraphCut.hpp b/src/aliceVision/fuseCut/DelaunayGraphCut.hpp index 0fff4c3583..01f24d159d 100644 --- a/src/aliceVision/fuseCut/DelaunayGraphCut.hpp +++ b/src/aliceVision/fuseCut/DelaunayGraphCut.hpp @@ -121,12 +121,6 @@ class DelaunayGraphCut {} }; - enum class EDirection - { - toTheCam, - behindThePoint - }; - struct IntersectionHistory { size_t steps = 0; From 30f8e92a56fdfad44a499822c2c634200837a860 Mon Sep 17 00:00:00 2001 From: Enguerrand DE SMET Date: Mon, 3 Aug 2020 12:36:27 +0200 Subject: [PATCH 14/26] [fuseCut] DelaunayGraphCut: export csv data for debugging purposes --- src/aliceVision/fuseCut/DelaunayGraphCut.cpp | 42 ++++++++++++++++++++ src/aliceVision/fuseCut/DelaunayGraphCut.hpp | 1 + 2 files changed, 43 insertions(+) diff --git a/src/aliceVision/fuseCut/DelaunayGraphCut.cpp b/src/aliceVision/fuseCut/DelaunayGraphCut.cpp index aa98cd8aea..6e1a24f5eb 100644 --- a/src/aliceVision/fuseCut/DelaunayGraphCut.cpp +++ b/src/aliceVision/fuseCut/DelaunayGraphCut.cpp @@ -3001,6 +3001,48 @@ void DelaunayGraphCut::exportBackPropagationMesh(const std::string& filename, st exportDebugMesh(filename + "_BackProp", fromPt, toPt); } +void DelaunayGraphCut::writeScoreInCsv(const std::string& filePath, const size_t& sizeLimit) +{ + assert(boost::filesystem::path(filePath).extension().string() == std::string(".csv")); + + const unsigned int seed = (unsigned int)mp->userParams.get("delaunaycut.seed", 0); + std::mt19937 generator(seed != 0 ? seed : std::random_device{}()); + + std::vector idx(_cellsAttr.size()); + std::iota(idx.begin(), idx.end(), 0); + std::shuffle(idx.begin(), idx.end(), generator); + + std::ofstream csv(filePath); + const char sep = ','; // separator + csv << "fullnessScore" << sep << + "emptinessScore" << sep << + "cellSWeight" << sep << + "cellTWeight" << sep << + "on" << sep << + "gEdgeVisWeight0" << sep << + "gEdgeVisWeight1" << sep << + "gEdgeVisWeight2" << sep << + "gEdgeVisWeight3" << '\n'; + const size_t size = sizeLimit > 0 ? std::min(sizeLimit, _cellsAttr.size()) : _cellsAttr.size(); + for (size_t i = 0; i < size; ++i) + { + const GC_cellInfo& cellAttr = _cellsAttr[idx.back()]; + idx.pop_back(); + csv << cellAttr.fullnessScore << sep << + cellAttr.emptinessScore << sep << + cellAttr.cellSWeight << sep << + cellAttr.cellTWeight << sep << + cellAttr.on << sep << + cellAttr.gEdgeVisWeight[0] << sep << + cellAttr.gEdgeVisWeight[1] << sep << + cellAttr.gEdgeVisWeight[2] << sep << + cellAttr.gEdgeVisWeight[3] << '\n'; + } + + csv.close(); + ALICEVISION_LOG_INFO("Csv exported: " << filePath); +} + void DelaunayGraphCut::segmentFullOrFree(bool full, StaticVector** 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 01f24d159d..4a53bef301 100644 --- a/src/aliceVision/fuseCut/DelaunayGraphCut.hpp +++ b/src/aliceVision/fuseCut/DelaunayGraphCut.hpp @@ -449,6 +449,7 @@ class DelaunayGraphCut void exportDebugMesh(const std::string& filename, const Point3d& fromPt, const Point3d& toPt); void exportFullScoreMeshs(); void exportBackPropagationMesh(const std::string& filename, std::vector& intersectedGeom, const Point3d& fromPt, const Point3d& toPt); + void writeScoreInCsv(const std::string& filePath, const size_t& sizeLimit = 1000); }; } // namespace fuseCut From 1992601047b9a8f3d4def7f4b0305324e3007bf5 Mon Sep 17 00:00:00 2001 From: Enguerrand DE SMET Date: Tue, 4 Aug 2020 15:52:02 +0200 Subject: [PATCH 15/26] [fuseCut] DelaunayGraphCut: renaming and more detailed log message for fillGraph and forceTedges --- src/aliceVision/fuseCut/DelaunayGraphCut.cpp | 156 +++++++++++++------ src/aliceVision/fuseCut/DelaunayGraphCut.hpp | 20 ++- 2 files changed, 125 insertions(+), 51 deletions(-) diff --git a/src/aliceVision/fuseCut/DelaunayGraphCut.cpp b/src/aliceVision/fuseCut/DelaunayGraphCut.cpp index 6e1a24f5eb..600e412c1c 100644 --- a/src/aliceVision/fuseCut/DelaunayGraphCut.cpp +++ b/src/aliceVision/fuseCut/DelaunayGraphCut.cpp @@ -1603,17 +1603,23 @@ void DelaunayGraphCut::fillGraph(bool fixesSigma, float nPixelSizeBehind, const unsigned int seed = (unsigned int)mp->userParams.get("delaunaycut.seed", 0); const std::vector verticesRandIds = mvsUtils::createRandomArrayOfIntegers(_verticesAttr.size(), seed); - int64_t avStepsFront = 0; - int64_t aAvStepsFront = 0; - int64_t avStepsBehind = 0; - int64_t nAvStepsBehind = 0; - int avCams = 0; - int nAvCams = 0; + int64_t totalStepsFront = 0; + int64_t totalRayFront = 0; + int64_t totalStepsBehind = 0; + int64_t totalRayBehind = 0; + + size_t totalCamHaveVisibilityOnVertex = 0; + size_t totalOfVertex = 0; + + size_t totalIsRealNrc = 0; + + GeometriesCount totalGeometriesIntersectedFrontCount; + GeometriesCount totalGeometriesIntersectedBehindCount; boost::progress_display progressBar(std::min(size_t(100), verticesRandIds.size()), std::cout, "fillGraphPartPtRc\n"); size_t progressStep = verticesRandIds.size() / 100; progressStep = std::max(size_t(1), progressStep); -#pragma omp parallel for reduction(+:avStepsFront,aAvStepsFront,avStepsBehind,nAvStepsBehind,avCams,nAvCams) +#pragma omp parallel for reduction(+:totalStepsFront,totalRayFront,totalStepsBehind,totalRayBehind,totalCamHaveVisibilityOnVertex,totalOfVertex) for(int i = 0; i < verticesRandIds.size(); i++) { if(i % progressStep == 0) @@ -1627,6 +1633,7 @@ void DelaunayGraphCut::fillGraph(bool fixesSigma, float nPixelSizeBehind, if(v.isReal() && (v.nrc > 0)) { + ++totalIsRealNrc; // "weight" is called alpha(p) in the paper const float weight = weightFcn((float)v.nrc, labatutWeights, v.getNbCameras()); // number of cameras @@ -1635,30 +1642,40 @@ void DelaunayGraphCut::fillGraph(bool fixesSigma, float nPixelSizeBehind, assert(v.cams[c] >= 0); assert(v.cams[c] < mp->ncams); - int nstepsFront = 0; - int nstepsBehind = 0; - fillGraphPartPtRc(nstepsFront, nstepsBehind, vertexIndex, v.cams[c], weight, fixesSigma, nPixelSizeBehind, + int stepsFront = 0; + int stepsBehind = 0; + GeometriesCount geometriesIntersectedFrontCount; + GeometriesCount geometriesIntersectedBehindCount; + fillGraphPartPtRc(stepsFront, stepsBehind, geometriesIntersectedFrontCount, geometriesIntersectedBehindCount, vertexIndex, v.cams[c], weight, fixesSigma, nPixelSizeBehind, fillOut, distFcnHeight); - avStepsFront += nstepsFront; - aAvStepsFront += 1; - avStepsBehind += nstepsBehind; - nAvStepsBehind += 1; + totalStepsFront += stepsFront; + totalRayFront += 1; + totalStepsBehind += stepsBehind; + totalRayBehind += 1; + + totalGeometriesIntersectedFrontCount += geometriesIntersectedFrontCount; + totalGeometriesIntersectedBehindCount += geometriesIntersectedBehindCount; } // for c - avCams += v.cams.size(); - nAvCams += 1; + totalCamHaveVisibilityOnVertex += v.cams.size(); + totalOfVertex += 1; } } - 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)); - ALICEVISION_LOG_DEBUG("avCams = " << mvsUtils::num2str(avCams) << " // " << mvsUtils::num2str(nAvCams)); + ALICEVISION_LOG_DEBUG("_verticesAttr.size(): " << _verticesAttr.size() << "(" << verticesRandIds.size() << ")"); + ALICEVISION_LOG_DEBUG("totalIsRealNrc: " << totalIsRealNrc); + ALICEVISION_LOG_DEBUG("totalStepsFront//totalRayFront = " << totalStepsFront << " // " << totalRayFront); + ALICEVISION_LOG_DEBUG("totalStepsBehind//totalRayBehind = " << totalStepsBehind << " // " << totalRayBehind); + ALICEVISION_LOG_DEBUG("totalCamHaveVisibilityOnVertex//totalOfVertex = " << totalCamHaveVisibilityOnVertex << " // " << totalOfVertex); + + ALICEVISION_LOG_DEBUG(" - Geometries Intersected count - \n"); + ALICEVISION_LOG_DEBUG("Front: edges " << totalGeometriesIntersectedFrontCount.edges << ", vertices: " << totalGeometriesIntersectedFrontCount.vertices << ", facets: " << totalGeometriesIntersectedFrontCount.facets); + ALICEVISION_LOG_DEBUG("Behind: edges " << totalGeometriesIntersectedBehindCount.edges << ", vertices: " << totalGeometriesIntersectedBehindCount.vertices << ", facets: " << totalGeometriesIntersectedBehindCount.facets); mvsUtils::printfElapsedTime(t1, "s-t graph weights computed : "); } -void DelaunayGraphCut::fillGraphPartPtRc(int& out_nstepsFront, int& out_nstepsBehind, int vertexIndex, int cam, +void DelaunayGraphCut::fillGraphPartPtRc(int& outTotalStepsFront, int& outTotalStepsBehind, GeometriesCount& outFrontCount, GeometriesCount& outBehindCount, int vertexIndex, int cam, float weight, bool fixesSigma, float nPixelSizeBehind, bool fillOut, float distFcnHeight) // fixesSigma=true nPixelSizeBehind=2*spaceSteps allPoints=1 behind=0 fillOut=1 distFcnHeight=0 { @@ -1673,8 +1690,6 @@ void DelaunayGraphCut::fillGraphPartPtRc(int& out_nstepsFront, int& out_nstepsBe assert(cam >= 0); assert(cam < mp->ncams); - int nsteps = 0; - if(fillOut) { // Initialisation @@ -1686,8 +1701,7 @@ void DelaunayGraphCut::fillGraphPartPtRc(int& out_nstepsFront, int& out_nstepsBe #ifdef ALICEVISION_DEBUG_VOTE IntersectionHistory history(mp->CArr[cam], originPt, dirVect); #endif - out_nstepsFront = 0; - + outTotalStepsFront = 0; Facet lastIntersectedFacet; // As long as we find a next geometry do @@ -1699,7 +1713,7 @@ void DelaunayGraphCut::fillGraphPartPtRc(int& out_nstepsFront, int& out_nstepsBe #ifdef ALICEVISION_DEBUG_VOTE history.append(geometry, intersectPt); #endif - ++out_nstepsFront; + ++outTotalStepsFront; geometry = intersectNextGeom(previousGeometry, originPt, dirVect, intersectPt, 1.0e-3, lastIntersectPt); @@ -1715,7 +1729,7 @@ void DelaunayGraphCut::fillGraphPartPtRc(int& out_nstepsFront, int& out_nstepsBe if(geometry.type == EGeometryType::Facet) { - lastIntersectedFacet = geometry.facet; + ++outFrontCount.facets; // Vote for the cell between the previous vertex/edge and the intersected facet // Because if there is an edge or a vertex, no vote has been made, we have to vote here for the right facet once we find it if(previousGeometry.type != EGeometryType::Facet) @@ -1744,6 +1758,14 @@ void DelaunayGraphCut::fillGraphPartPtRc(int& out_nstepsFront, int& out_nstepsBe #pragma OMP_ATOMIC_UPDATE _cellsAttr[geometry.facet.cellIndex].emptinessScore += weight; } + else if (geometry.type == EGeometryType::Vertex) + { + ++outFrontCount.vertices; + } + else if (geometry.type == EGeometryType::Edge) + { + ++outFrontCount.edges; + } // Break only when we reach our camera vertex } while(!(geometry.type == EGeometryType::Vertex && (mp->CArr[cam] - intersectPt).size() < 1.0e-3)); @@ -1765,7 +1787,7 @@ void DelaunayGraphCut::fillGraphPartPtRc(int& out_nstepsFront, int& out_nstepsBe #ifdef ALICEVISION_DEBUG_VOTE IntersectionHistory history(mp->CArr[cam], originPt, dirVect); #endif - out_nstepsBehind = 0; + outTotalStepsBehind = 0; bool firstIteration = true; Facet lastIntersectedFacet; @@ -1779,10 +1801,10 @@ void DelaunayGraphCut::fillGraphPartPtRc(int& out_nstepsFront, int& out_nstepsBe #ifdef ALICEVISION_DEBUG_VOTE history.append(geometry, intersectPt); #endif - ++out_nstepsBehind; + ++outTotalStepsBehind; // Vote for the first facet found (only once) - if(firstIteration && geometry.type == EGeometryType::Facet) + if (firstIteration && geometry.type == EGeometryType::Facet) { #pragma OMP_ATOMIC_UPDATE _cellsAttr[geometry.facet.cellIndex].on += weight; @@ -1807,6 +1829,7 @@ void DelaunayGraphCut::fillGraphPartPtRc(int& out_nstepsFront, int& out_nstepsBe if(geometry.type == EGeometryType::Facet) { + ++outBehindCount.facets; lastIntersectedFacet = geometry.facet; // Vote for the cell between the previous vertex/edge and the intersected facet @@ -1834,6 +1857,14 @@ void DelaunayGraphCut::fillGraphPartPtRc(int& out_nstepsFront, int& out_nstepsBe #pragma OMP_ATOMIC_UPDATE _cellsAttr[geometry.facet.cellIndex].fullnessScore += weight; } + else if (geometry.type == EGeometryType::Vertex) + { + ++outBehindCount.vertices; + } + else if (geometry.type == EGeometryType::Edge) + { + ++outBehindCount.edges; + } // While we are within the surface margin } while((originPt - intersectPt).size() < maxDist); @@ -1871,7 +1902,6 @@ void DelaunayGraphCut::forceTedgesByGradientIJCV(bool fixesSigma, float nPixelSi float nsigmaBackSilentPart = (float)mp->userParams.get("delaunaycut.nsigmaBackSilentPart", 2.0f); ALICEVISION_LOG_DEBUG("nsigmaBackSilentPart: " << nsigmaBackSilentPart); - for(GC_cellInfo& c: _cellsAttr) { c.on = 0.0f; @@ -1883,12 +1913,20 @@ void DelaunayGraphCut::forceTedgesByGradientIJCV(bool fixesSigma, float nPixelSi const unsigned int seed = (unsigned int)mp->userParams.get("delaunaycut.seed", 0); const std::vector verticesRandIds = mvsUtils::createRandomArrayOfIntegers(_verticesAttr.size(), seed); - int64_t avStepsFront = 0; - int64_t aAvStepsFront = 0; - int64_t avStepsBehind = 0; - int64_t nAvStepsBehind = 0; + size_t totalStepsFront = 0; + size_t totalRayFront = 0; + size_t totalStepsBehind = 0; + size_t totalRayBehind = 0; -#pragma omp parallel for reduction(+:avStepsFront,aAvStepsFront,avStepsBehind,nAvStepsBehind) + size_t totalCamHaveVisibilityOnVertex = 0; + size_t totalOfVertex = 0; + + size_t totalVertexIsVirtual = 0; + + GeometriesCount totalGeometriesIntersectedFrontCount; + GeometriesCount totalGeometriesIntersectedBehindCount; + +#pragma omp parallel for reduction(+:totalStepsFront,totalRayFront,totalStepsBehind,totalRayBehind) for(int i = 0; i < verticesRandIds.size(); ++i) { const int vertexIndex = verticesRandIds[i]; @@ -1896,6 +1934,7 @@ void DelaunayGraphCut::forceTedgesByGradientIJCV(bool fixesSigma, float nPixelSi if(v.isVirtual()) continue; + ++totalVertexIsVirtual; const Point3d& originPt = _verticesCoords[vertexIndex]; // For each camera that has visibility over the vertex v (vertexIndex) for(const int cam : v.cams) @@ -1918,8 +1957,6 @@ void DelaunayGraphCut::forceTedgesByGradientIJCV(bool fixesSigma, float nPixelSi #ifdef ALICEVISION_DEBUG_VOTE IntersectionHistory history(mp->CArr[cam], originPt, dirVect); #endif - int nstepsFront = 0; - // As long as we find a next geometry do { @@ -1930,7 +1967,7 @@ void DelaunayGraphCut::forceTedgesByGradientIJCV(bool fixesSigma, float nPixelSi #ifdef ALICEVISION_DEBUG_VOTE history.append(geometry, intersectPt); #endif - ++nstepsFront; + ++totalStepsFront; geometry = intersectNextGeom(previousGeometry, originPt, dirVect, intersectPt, 1.0e-3, lastIntersectPt); @@ -1946,6 +1983,7 @@ void DelaunayGraphCut::forceTedgesByGradientIJCV(bool fixesSigma, float nPixelSi if(geometry.type == EGeometryType::Facet) { + ++totalGeometriesIntersectedFrontCount.facets; const GC_cellInfo& c = _cellsAttr[geometry.facet.cellIndex]; if((intersectPt - originPt).size() > nsigmaFrontSilentPart * maxDist) // (p-originPt).size() > 2 * sigma { @@ -1971,13 +2009,20 @@ void DelaunayGraphCut::forceTedgesByGradientIJCV(bool fixesSigma, float nPixelSi } geometry.facet = mFacet; } + else if (geometry.type == EGeometryType::Vertex) + { + ++totalGeometriesIntersectedFrontCount.vertices; + } + else if (geometry.type == EGeometryType::Edge) + { + ++totalGeometriesIntersectedFrontCount.edges; + } // Break only when we reach our camera vertex or our distance condition is unsatisfied } while(!(geometry.type == EGeometryType::Vertex && (mp->CArr[cam] - intersectPt).size() < 1.0e-3) && (intersectPt - originPt).size() <= (nsigmaJumpPart + nsigmaFrontSilentPart) * maxDist); - avStepsFront += nstepsFront; - aAvStepsFront += 1; + ++totalRayFront; } { // Initialisation @@ -1989,7 +2034,6 @@ void DelaunayGraphCut::forceTedgesByGradientIJCV(bool fixesSigma, float nPixelSi #ifdef ALICEVISION_DEBUG_VOTE IntersectionHistory history(mp->CArr[cam], originPt, dirVect); #endif - int nstepsBehind = 0; Facet lastIntersectedFacet; // As long as we find a next geometry @@ -2002,7 +2046,7 @@ void DelaunayGraphCut::forceTedgesByGradientIJCV(bool fixesSigma, float nPixelSi #ifdef ALICEVISION_DEBUG_VOTE history.append(geometry, intersectPt); #endif - ++nstepsBehind; + ++totalStepsBehind; geometry = intersectNextGeom(previousGeometry, originPt, dirVect, intersectPt, 1.0e-3, lastIntersectPt); @@ -2024,6 +2068,8 @@ void DelaunayGraphCut::forceTedgesByGradientIJCV(bool fixesSigma, float nPixelSi if(geometry.type == EGeometryType::Facet) { + ++totalGeometriesIntersectedFrontCount.facets; + const GC_cellInfo& c = _cellsAttr[geometry.facet.cellIndex]; minSilent = std::min(minSilent, c.emptinessScore); maxSilent = std::max(maxSilent, c.emptinessScore); @@ -2037,6 +2083,14 @@ void DelaunayGraphCut::forceTedgesByGradientIJCV(bool fixesSigma, float nPixelSi } geometry.facet = mFacet; } + else if (geometry.type == EGeometryType::Vertex) + { + ++totalGeometriesIntersectedFrontCount.vertices; + } + else if (geometry.type == EGeometryType::Edge) + { + ++totalGeometriesIntersectedFrontCount.edges; + } } while((intersectPt - originPt).size() <= nsigmaBackSilentPart * maxDist); @@ -2066,11 +2120,11 @@ void DelaunayGraphCut::forceTedgesByGradientIJCV(bool fixesSigma, float nPixelSi _cellsAttr[lastIntersectedFacet.cellIndex].on += (maxJump - midSilent); } } - - avStepsBehind += nstepsBehind; - nAvStepsBehind += 1; + ++totalRayBehind; } } + totalCamHaveVisibilityOnVertex += v.cams.size(); + totalOfVertex += 1; } for(GC_cellInfo& c: _cellsAttr) @@ -2082,10 +2136,12 @@ void DelaunayGraphCut::forceTedgesByGradientIJCV(bool fixesSigma, float nPixelSi // c.cellTWeight = std::max(c.cellTWeight,fit->info().on); } - { - ALICEVISION_LOG_DEBUG("avStepsFront = " << avStepsFront << " // " << aAvStepsFront); - ALICEVISION_LOG_DEBUG("avStepsBehind = " << avStepsBehind << " // " << nAvStepsBehind); - } + ALICEVISION_LOG_DEBUG("_verticesAttr.size(): " << _verticesAttr.size() << "(" << verticesRandIds.size() << ")"); + ALICEVISION_LOG_DEBUG("totalVertexIsVirtual: " << totalVertexIsVirtual); + ALICEVISION_LOG_DEBUG("totalStepsFront//totalRayFront = " << totalStepsFront << " // " << totalRayFront); + ALICEVISION_LOG_DEBUG("totalStepsBehind//totalRayBehind = " << totalStepsBehind << " // " << totalRayBehind); + ALICEVISION_LOG_DEBUG("totalCamHaveVisibilityOnVertex//totalOfVertex = " << totalCamHaveVisibilityOnVertex << " // " << totalOfVertex); + mvsUtils::printfElapsedTime(t2, "t-edges forced: "); } diff --git a/src/aliceVision/fuseCut/DelaunayGraphCut.hpp b/src/aliceVision/fuseCut/DelaunayGraphCut.hpp index 4a53bef301..2fb95773cd 100644 --- a/src/aliceVision/fuseCut/DelaunayGraphCut.hpp +++ b/src/aliceVision/fuseCut/DelaunayGraphCut.hpp @@ -138,6 +138,24 @@ class DelaunayGraphCut void append(const GeometryIntersection& geom, const Point3d& intersectPt); }; + /** + * @brief Used for debug purposes to store count about geometries intersected during fillGraph and forceTedges. + */ + struct GeometriesCount + { + size_t facets = 0; + size_t vertices = 0; + size_t edges = 0; + + GeometriesCount& operator+=(const GeometriesCount& gc) + { + edges += gc.edges; + vertices += gc.vertices; + facets += gc.facets; + + return *this; + } + }; mvsUtils::MultiViewParams* mp; @@ -410,7 +428,7 @@ class DelaunayGraphCut 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, + void fillGraphPartPtRc(int& out_nstepsFront, int& out_nstepsBehind, GeometriesCount& outFrontCount, GeometriesCount& outBehindCount, int vertexIndex, int cam, float weight, bool fixesSigma, float nPixelSizeBehind, bool fillOut, float distFcnHeight); From ca5864676c657d2f883b2674370acc043501f670 Mon Sep 17 00:00:00 2001 From: Enguerrand DE SMET Date: Wed, 5 Aug 2020 19:37:48 +0200 Subject: [PATCH 16/26] [fuseCut] DelaunayGraphCut: generalize getNeighboringCells for geometries --- src/aliceVision/fuseCut/DelaunayGraphCut.cpp | 28 ++++++++++++++++++++ src/aliceVision/fuseCut/DelaunayGraphCut.hpp | 17 ++++++++++++ 2 files changed, 45 insertions(+) diff --git a/src/aliceVision/fuseCut/DelaunayGraphCut.cpp b/src/aliceVision/fuseCut/DelaunayGraphCut.cpp index 600e412c1c..c14969496d 100644 --- a/src/aliceVision/fuseCut/DelaunayGraphCut.cpp +++ b/src/aliceVision/fuseCut/DelaunayGraphCut.cpp @@ -462,6 +462,34 @@ void DelaunayGraphCut::saveDh(const std::string& fileNameDh, const std::string& mvsUtils::printfElapsedTime(t1); } + +std::vector DelaunayGraphCut::getNeighboringCellsByGeometry(const GeometryIntersection& g) const +{ + switch (g.type) + { + case EGeometryType::Edge: + return getNeighboringCellsByEdge(g.edge); + case EGeometryType::Vertex: + return getNeighboringCellsByVertexIndex(g.vertexIndex); + case EGeometryType::Facet: + return getNeighboringCellsByFacet(g.facet); + case EGeometryType::None: + throw std::runtime_error("[error] getNeighboringCellsByGeometry: an undefined/None geometry has no neighboring cells."); + } +} + +std::vector DelaunayGraphCut::getNeighboringCellsByFacet(const Facet& f) const +{ + std::vector neighboringCells; + neighboringCells.push_back(f.cellIndex); + + const Facet mFacet = mirrorFacet(f); + if(!isInvalidOrInfiniteCell(mFacet.cellIndex)) + neighboringCells.push_back(mFacet.cellIndex); + + return neighboringCells; +} + std::vector DelaunayGraphCut::getNeighboringCellsByEdge(const Edge& e) const { const std::vector& v0ci = getNeighboringCellsByVertexIndex(e.v0); diff --git a/src/aliceVision/fuseCut/DelaunayGraphCut.hpp b/src/aliceVision/fuseCut/DelaunayGraphCut.hpp index 2fb95773cd..8c903e57ca 100644 --- a/src/aliceVision/fuseCut/DelaunayGraphCut.hpp +++ b/src/aliceVision/fuseCut/DelaunayGraphCut.hpp @@ -336,6 +336,23 @@ class DelaunayGraphCut return localCells[lvi]; } + + /** + * @brief Retrieves the global indexes of neighboring cells around a geometry. + * + * @param g the concerned geometry + * @return a vector of neighboring cell indices + */ + std::vector getNeighboringCellsByGeometry(const GeometryIntersection& g) const; + + /** + * @brief Retrieves the two global indexes of neighboring cells using a facet. + * + * @param f the concerned facet + * @return a vector of neighboring cell indices + */ + std::vector getNeighboringCellsByFacet(const Facet& f) const; + /** * @brief Retrieves the global indexes of neighboring cells using the global index of a vertex. * From 6b7218f788b78956ba189441884b2330d3597d67 Mon Sep 17 00:00:00 2001 From: Enguerrand DE SMET Date: Wed, 5 Aug 2020 19:39:44 +0200 Subject: [PATCH 17/26] [fuseCut] DelaunayGraphCut: add equality operators for few struct --- src/aliceVision/fuseCut/DelaunayGraphCut.hpp | 39 +++++++++++++++++--- 1 file changed, 34 insertions(+), 5 deletions(-) diff --git a/src/aliceVision/fuseCut/DelaunayGraphCut.hpp b/src/aliceVision/fuseCut/DelaunayGraphCut.hpp index 8c903e57ca..361cf21653 100644 --- a/src/aliceVision/fuseCut/DelaunayGraphCut.hpp +++ b/src/aliceVision/fuseCut/DelaunayGraphCut.hpp @@ -66,27 +66,38 @@ class DelaunayGraphCut struct Facet { + CellIndex cellIndex = GEO::NO_CELL; + /// local opposite vertex index + VertexIndex localVertexIndex = GEO::NO_VERTEX; + Facet(){} Facet(CellIndex ci, VertexIndex lvi) : cellIndex(ci) , localVertexIndex(lvi) {} - CellIndex cellIndex = GEO::NO_CELL; - /// local opposite vertex index - VertexIndex localVertexIndex = GEO::NO_VERTEX; + bool operator==(const Facet& f) const + { + return cellIndex == f.cellIndex && localVertexIndex == f.localVertexIndex; + } }; struct Edge { + VertexIndex v0 = GEO::NO_VERTEX; + VertexIndex v1 = GEO::NO_VERTEX; + Edge() = default; Edge(VertexIndex v0_, VertexIndex v1_) : v0{v0_} , v1{v1_} {} - VertexIndex v0 = GEO::NO_VERTEX; - VertexIndex v1 = GEO::NO_VERTEX; + bool operator==(const Edge& e) const + { + return v0 == e.v0 && v1 == e.v1; + } + }; enum class EGeometryType @@ -119,6 +130,24 @@ class DelaunayGraphCut : edge{e} , type{EGeometryType::Edge} {} + + bool operator==(const GeometryIntersection& g) const + { + if (type != g.type) + return false; + + switch (type) + { + case EGeometryType::Vertex: + return vertexIndex == g.vertexIndex; + case EGeometryType::Edge: + return edge == g.edge; + case EGeometryType::Facet: + return facet == g.facet; + case EGeometryType::None: + return true; + } + } }; struct IntersectionHistory From 64db42f99efa627368b6455dfc59f1c02cb7d6b0 Mon Sep 17 00:00:00 2001 From: Enguerrand DE SMET Date: Thu, 6 Aug 2020 18:58:43 +0200 Subject: [PATCH 18/26] [fuseCut] DelaunayGraphCut: fix new implementation --- src/aliceVision/fuseCut/DelaunayGraphCut.cpp | 222 ++++++++++++------- 1 file changed, 143 insertions(+), 79 deletions(-) diff --git a/src/aliceVision/fuseCut/DelaunayGraphCut.cpp b/src/aliceVision/fuseCut/DelaunayGraphCut.cpp index c14969496d..5a5a53bfbc 100644 --- a/src/aliceVision/fuseCut/DelaunayGraphCut.cpp +++ b/src/aliceVision/fuseCut/DelaunayGraphCut.cpp @@ -1731,8 +1731,8 @@ void DelaunayGraphCut::fillGraphPartPtRc(int& outTotalStepsFront, int& outTotalS #endif outTotalStepsFront = 0; Facet lastIntersectedFacet; - // As long as we find a next geometry - do + // Break only when we reach our camera vertex (as long as we find a next geometry) + while (geometry.type != EGeometryType::Vertex || (mp->CArr[cam] - intersectPt).size() >= 1.0e-3) { // Keep previous informations const GeometryIntersection previousGeometry = geometry; @@ -1744,7 +1744,7 @@ void DelaunayGraphCut::fillGraphPartPtRc(int& outTotalStepsFront, int& outTotalS ++outTotalStepsFront; geometry = intersectNextGeom(previousGeometry, originPt, dirVect, intersectPt, 1.0e-3, lastIntersectPt); - + if (geometry.type == EGeometryType::None) { #ifdef ALICEVISION_DEBUG_VOTE @@ -1755,22 +1755,23 @@ void DelaunayGraphCut::fillGraphPartPtRc(int& outTotalStepsFront, int& outTotalS break; } - if(geometry.type == EGeometryType::Facet) + if (geometry.type == EGeometryType::Facet) { ++outFrontCount.facets; - // Vote for the cell between the previous vertex/edge and the intersected facet - // Because if there is an edge or a vertex, no vote has been made, we have to vote here for the right facet once we find it - if(previousGeometry.type != EGeometryType::Facet) { +#pragma OMP_ATOMIC_UPDATE _cellsAttr[geometry.facet.cellIndex].emptinessScore += weight; } - const float dist = distFcn(maxDist, (originPt - intersectPt).size(), distFcnHeight); + { + const float dist = distFcn(maxDist, (originPt - lastIntersectPt).size(), distFcnHeight); #pragma OMP_ATOMIC_UPDATE - _cellsAttr[geometry.facet.cellIndex].gEdgeVisWeight[geometry.facet.localVertexIndex] += weight * dist; + _cellsAttr[geometry.facet.cellIndex].gEdgeVisWeight[geometry.facet.localVertexIndex] += weight * dist; + } // Take the mirror facet to iterate over the next cell const Facet mFacet = mirrorFacet(geometry.facet); + geometry.facet = mFacet; if (isInvalidOrInfiniteCell(mFacet.cellIndex)) { #ifdef ALICEVISION_DEBUG_VOTE @@ -1780,23 +1781,28 @@ void DelaunayGraphCut::fillGraphPartPtRc(int& outTotalStepsFront, int& outTotalS ALICEVISION_LOG_TRACE("[Error]: fillGraphPartPtRc toTheCam, cause: invalidOrInfinite miror facet."); break; } - geometry.facet = mFacet; - - // Vote for the new cell after taking the mirror facet -#pragma OMP_ATOMIC_UPDATE - _cellsAttr[geometry.facet.cellIndex].emptinessScore += weight; + lastIntersectedFacet = mFacet; } - else if (geometry.type == EGeometryType::Vertex) - { - ++outFrontCount.vertices; - } - else if (geometry.type == EGeometryType::Edge) + else { - ++outFrontCount.edges; - } + // We have just intersected a vertex or edge. + // These geometries do not have a cellIndex, so we use the previousGeometry to retrieve the cell between the previous geometry and the current one. + if (previousGeometry.type == EGeometryType::Facet) + { +#pragma OMP_ATOMIC_UPDATE + _cellsAttr[previousGeometry.facet.cellIndex].emptinessScore += weight; + } - // Break only when we reach our camera vertex - } while(!(geometry.type == EGeometryType::Vertex && (mp->CArr[cam] - intersectPt).size() < 1.0e-3)); + if (geometry.type == EGeometryType::Vertex) + { + ++outFrontCount.vertices; + } + else if (geometry.type == EGeometryType::Edge) + { + ++outFrontCount.edges; + } + } + } // Vote for the last intersected facet (close to the cam) if (lastIntersectedFacet.cellIndex != GEO::NO_CELL) @@ -1819,8 +1825,8 @@ void DelaunayGraphCut::fillGraphPartPtRc(int& outTotalStepsFront, int& outTotalS bool firstIteration = true; Facet lastIntersectedFacet; - // As long as we find a next geometry - do + // While we are within the surface margin (as long as we find a next geometry) + while ((originPt - intersectPt).size() < maxDist) { // Keep previous informations const GeometryIntersection previousGeometry = geometry; @@ -1831,19 +1837,12 @@ void DelaunayGraphCut::fillGraphPartPtRc(int& outTotalStepsFront, int& outTotalS #endif ++outTotalStepsBehind; - // Vote for the first facet found (only once) - if (firstIteration && geometry.type == EGeometryType::Facet) - { -#pragma OMP_ATOMIC_UPDATE - _cellsAttr[geometry.facet.cellIndex].on += weight; - firstIteration = false; - } geometry = intersectNextGeom(previousGeometry, originPt, dirVect, intersectPt, 1.0e-3, lastIntersectPt); - if(geometry.type == EGeometryType::None) + if (geometry.type == EGeometryType::None) { // If we come from a facet, the next intersection must exist (even if the mirror facet is invalid, which is verified after taking mirror facet) - if(previousGeometry.type == EGeometryType::Facet) + if (previousGeometry.type == EGeometryType::Facet) { #ifdef ALICEVISION_DEBUG_VOTE // exportBackPropagationMesh("fillGraph_behindThePoint_NoneButPreviousIsFacet", history.geometries, originPt, mp->CArr[cam]); @@ -1855,53 +1854,86 @@ void DelaunayGraphCut::fillGraphPartPtRc(int& outTotalStepsFront, int& outTotalS break; } - if(geometry.type == EGeometryType::Facet) + if (geometry.type == EGeometryType::Facet) { ++outBehindCount.facets; - lastIntersectedFacet = geometry.facet; - // Vote for the cell between the previous vertex/edge and the intersected facet - // Because if there is an edge or a vertex, no vote has been made, we have to vote here for the right facet once we find it - if (previousGeometry.type != EGeometryType::Facet) + // Vote for the first cell found (only once) + if (firstIteration) { #pragma OMP_ATOMIC_UPDATE - _cellsAttr[geometry.facet.cellIndex].fullnessScore += weight; + _cellsAttr[geometry.facet.cellIndex].on += weight; + firstIteration = false; } - const float dist = distFcn(maxDist, (originPt - intersectPt).size(), distFcnHeight); + { #pragma OMP_ATOMIC_UPDATE - _cellsAttr[geometry.facet.cellIndex].gEdgeVisWeight[geometry.facet.localVertexIndex] += weight * dist; - + _cellsAttr[geometry.facet.cellIndex].fullnessScore += weight; + } + // Take the mirror facet to iterate over the next cell const Facet mFacet = mirrorFacet(geometry.facet); + lastIntersectedFacet = mFacet; + geometry.facet = mFacet; if (isInvalidOrInfiniteCell(mFacet.cellIndex)) { // Break if we reach the end of the tetrahedralization volume (mirror facet cannot be found) break; } - geometry.facet = mFacet; - // Vote for the new cell after taking the mirror facet + { + const float dist = distFcn(maxDist, (originPt - lastIntersectPt).size(), distFcnHeight); #pragma OMP_ATOMIC_UPDATE - _cellsAttr[geometry.facet.cellIndex].fullnessScore += weight; - } - else if (geometry.type == EGeometryType::Vertex) - { - ++outBehindCount.vertices; + _cellsAttr[geometry.facet.cellIndex].gEdgeVisWeight[geometry.facet.localVertexIndex] += weight * dist; + } } - else if (geometry.type == EGeometryType::Edge) + else { - ++outBehindCount.edges; - } + // Vote for the first cell found (only once) + // if we come from an edge or vertex to an other we have to vote for the first intersected cell. + if (firstIteration) + { + if (previousGeometry.type != EGeometryType::Vertex) + throw std::runtime_error("[error] The firstIteration vote could only happen during for the first cell when we come from the first vertex."); + // the information of first intersected cell can only be found by taking intersection of neighbouring cells for both geometries + const std::vector previousNeighbouring = getNeighboringCellsByVertexIndex(previousGeometry.vertexIndex); + const std::vector currentNeigbouring = getNeighboringCellsByGeometry(geometry); - // While we are within the surface margin - } while((originPt - intersectPt).size() < maxDist); + std::vector neighboringCells; + std::set_intersection(previousNeighbouring.begin(), previousNeighbouring.end(), currentNeigbouring.begin(), currentNeigbouring.end(), std::back_inserter(neighboringCells)); + + for (const CellIndex& ci : neighboringCells) + { +#pragma OMP_ATOMIC_UPDATE + _cellsAttr[neighboringCells[0]].on += weight; + } + firstIteration = false; + } + + // We have just intersected a vertex or edge. + // These geometries do not have a cellIndex, so we use the previousGeometry to retrieve the cell between the previous geometry and the current one. + if (previousGeometry.type == EGeometryType::Facet) + { +#pragma OMP_ATOMIC_UPDATE + _cellsAttr[previousGeometry.facet.cellIndex].fullnessScore += weight; + } + + if (geometry.type == EGeometryType::Vertex) + { + ++outBehindCount.vertices; + } + else if (geometry.type == EGeometryType::Edge) + { + ++outBehindCount.edges; + } + } + } // cv: is the tetrahedron in distance 2*sigma behind the point p in the direction of the camera c (called Lcp in the paper) Vote for the last found facet // Vote for the last intersected facet (farthest from the camera) if (lastIntersectedFacet.cellIndex != GEO::NO_CELL) { -#pragma OMP_ATOMIC_WRITE +#pragma OMP_ATOMIC_UPDATE _cellsAttr[lastIntersectedFacet.cellIndex].cellTWeight += weight; } } @@ -1927,6 +1959,7 @@ void DelaunayGraphCut::forceTedgesByGradientIJCV(bool fixesSigma, float nPixelSi float nsigmaFrontSilentPart = (float)mp->userParams.get("delaunaycut.nsigmaFrontSilentPart", 2.0f); ALICEVISION_LOG_DEBUG("nsigmaFrontSilentPart: " << nsigmaFrontSilentPart); + // This parameter allows to enlage the surface margin behind the point float nsigmaBackSilentPart = (float)mp->userParams.get("delaunaycut.nsigmaBackSilentPart", 2.0f); ALICEVISION_LOG_DEBUG("nsigmaBackSilentPart: " << nsigmaBackSilentPart); @@ -1986,11 +2019,14 @@ void DelaunayGraphCut::forceTedgesByGradientIJCV(bool fixesSigma, float nPixelSi IntersectionHistory history(mp->CArr[cam], originPt, dirVect); #endif // As long as we find a next geometry - do + Point3d lastIntersectPt = originPt; + // Iterate on geometries in the direction of camera's vertex within margin defined by maxDist (as long as we find a next geometry) + while ((geometry.type != EGeometryType::Vertex || (mp->CArr[cam] - intersectPt).size() < 1.0e-3) // We reach our camera vertex + && (lastIntersectPt - originPt).size() <= (nsigmaJumpPart + nsigmaFrontSilentPart) * maxDist) // We to far from the originPt { // Keep previous informations const GeometryIntersection previousGeometry = geometry; - const Point3d lastIntersectPt = intersectPt; + lastIntersectPt = intersectPt; #ifdef ALICEVISION_DEBUG_VOTE history.append(geometry, intersectPt); @@ -2009,11 +2045,11 @@ void DelaunayGraphCut::forceTedgesByGradientIJCV(bool fixesSigma, float nPixelSi break; } - if(geometry.type == EGeometryType::Facet) + if (geometry.type == EGeometryType::Facet) { ++totalGeometriesIntersectedFrontCount.facets; const GC_cellInfo& c = _cellsAttr[geometry.facet.cellIndex]; - if((intersectPt - originPt).size() > nsigmaFrontSilentPart * maxDist) // (p-originPt).size() > 2 * sigma + if ((lastIntersectPt - originPt).size() > nsigmaFrontSilentPart * maxDist) // (p-originPt).size() > 2 * sigma { minJump = std::min(minJump, c.emptinessScore); maxJump = std::max(maxJump, c.emptinessScore); @@ -2045,11 +2081,7 @@ void DelaunayGraphCut::forceTedgesByGradientIJCV(bool fixesSigma, float nPixelSi { ++totalGeometriesIntersectedFrontCount.edges; } - - // Break only when we reach our camera vertex or our distance condition is unsatisfied - } while(!(geometry.type == EGeometryType::Vertex && (mp->CArr[cam] - intersectPt).size() < 1.0e-3) - && (intersectPt - originPt).size() <= (nsigmaJumpPart + nsigmaFrontSilentPart) * maxDist); - + } ++totalRayFront; } { @@ -2064,12 +2096,15 @@ void DelaunayGraphCut::forceTedgesByGradientIJCV(bool fixesSigma, float nPixelSi #endif Facet lastIntersectedFacet; - // As long as we find a next geometry - do + bool firstIteration = true; + Point3d lastIntersectPt = originPt; + + // While we are within the surface margin defined by maxDist (as long as we find a next geometry) + while ((lastIntersectPt - originPt).size() <= nsigmaBackSilentPart * maxDist) { // Keep previous informations const GeometryIntersection previousGeometry = geometry; - const Point3d lastIntersectPt = intersectPt; + lastIntersectPt = intersectPt; #ifdef ALICEVISION_DEBUG_VOTE history.append(geometry, intersectPt); @@ -2080,7 +2115,6 @@ void DelaunayGraphCut::forceTedgesByGradientIJCV(bool fixesSigma, float nPixelSi if(geometry.type == EGeometryType::None) { - lastIntersectedFacet = geometry.facet; // If we come from a facet, the next intersection must exist (even if the mirror facet is invalid, which is verified later) if (previousGeometry.type == EGeometryType::Facet) { @@ -2098,29 +2132,59 @@ void DelaunayGraphCut::forceTedgesByGradientIJCV(bool fixesSigma, float nPixelSi { ++totalGeometriesIntersectedFrontCount.facets; + // Vote for the first cell found (only once) + if (firstIteration) + { + midSilent = _cellsAttr[geometry.facet.cellIndex].emptinessScore; + firstIteration = false; + } + const GC_cellInfo& c = _cellsAttr[geometry.facet.cellIndex]; minSilent = std::min(minSilent, c.emptinessScore); maxSilent = std::max(maxSilent, c.emptinessScore); // Take the mirror facet to iterate over the next cell const Facet mFacet = mirrorFacet(geometry.facet); + lastIntersectedFacet = mFacet; + geometry.facet = mFacet; if (isInvalidOrInfiniteCell(mFacet.cellIndex)) { // Break if we reach the end of the tetrahedralization volume (mirror facet cannot be found) break; } - geometry.facet = mFacet; - } - else if (geometry.type == EGeometryType::Vertex) - { - ++totalGeometriesIntersectedFrontCount.vertices; - } - else if (geometry.type == EGeometryType::Edge) - { - ++totalGeometriesIntersectedFrontCount.edges; } + else { + + // Vote for the first cell found (only once) + // if we come from an edge or vertex to an other we have to vote for the first intersected cell. + if (firstIteration) + { + if (previousGeometry.type != EGeometryType::Vertex) + throw std::runtime_error("[error] The firstIteration vote could only happen during for the first cell when we come from the first vertex."); + // the information of first intersected cell can only be found by taking intersection of neighbouring cells for both geometries + const std::vector previousNeighbouring = getNeighboringCellsByVertexIndex(previousGeometry.vertexIndex); + const std::vector currentNeigbouring = getNeighboringCellsByGeometry(geometry); - } while((intersectPt - originPt).size() <= nsigmaBackSilentPart * maxDist); + std::vector neighboringCells; + std::set_intersection(previousNeighbouring.begin(), previousNeighbouring.end(), currentNeigbouring.begin(), currentNeigbouring.end(), std::back_inserter(neighboringCells)); + + for (const CellIndex& ci : neighboringCells) + { + midSilent = _cellsAttr[geometry.facet.cellIndex].emptinessScore; + } + firstIteration = false; + } + + if (geometry.type == EGeometryType::Vertex) + { + ++totalGeometriesIntersectedFrontCount.vertices; + } + else if (geometry.type == EGeometryType::Edge) + { + ++totalGeometriesIntersectedFrontCount.edges; + } + } + } if (lastIntersectedFacet.cellIndex != GEO::NO_CELL) { From e26fa797b8894fdf0574cc773073b4fa8b66dd2e Mon Sep 17 00:00:00 2001 From: Enguerrand DE SMET Date: Fri, 7 Aug 2020 12:46:08 +0200 Subject: [PATCH 19/26] [mvsData] getLineTriangleIntersectBarycCoords: use double instead of float reduce the number of failure cases --- src/aliceVision/mvsData/geometry.cpp | 66 ++++++++++++++-------------- 1 file changed, 33 insertions(+), 33 deletions(-) diff --git a/src/aliceVision/mvsData/geometry.cpp b/src/aliceVision/mvsData/geometry.cpp index e388501a81..c503375f3a 100644 --- a/src/aliceVision/mvsData/geometry.cpp +++ b/src/aliceVision/mvsData/geometry.cpp @@ -263,42 +263,42 @@ void rotPointAroundVect(Point3d* out, const Point3d* X, const Point3d* vect, con Point2d getLineTriangleIntersectBarycCoords(Point3d* P, const Point3d* A, const Point3d* B, const Point3d* C, const Point3d* linePoint, const Point3d* lineVect) { - float A_x = A->x; - float A_y = A->y; - float A_z = A->z; - float linePoint_x = linePoint->x; - float linePoint_y = linePoint->y; - float linePoint_z = linePoint->z; - float lineVect_x = lineVect->x; - float lineVect_y = lineVect->y; - float lineVect_z = lineVect->z; - float v0_x = C->x - A_x; - float v0_y = C->y - A_y; - float v0_z = C->z - A_z; - float v1_x = B->x - A_x; - float v1_y = B->y - A_y; - float v1_z = B->z - A_z; - float _n_x = v0_y * v1_z - v0_z * v1_y; - float _n_y = v0_z * v1_x - v0_x * v1_z; - float _n_z = v0_x * v1_y - v0_y * v1_x; - float k = ((A_x * _n_x + A_y * _n_y + A_z * _n_z) - (_n_x * linePoint_x + _n_y * linePoint_y + _n_z * linePoint_z)) / (_n_x * lineVect_x + _n_y * lineVect_y + _n_z * lineVect_z); - float P_x = linePoint_x + lineVect_x * k; - float P_y = linePoint_y + lineVect_y * k; - float P_z = linePoint_z + lineVect_z * k; + const double A_x = A->x; + const double A_y = A->y; + const double A_z = A->z; + const double linePoint_x = linePoint->x; + const double linePoint_y = linePoint->y; + const double linePoint_z = linePoint->z; + const double lineVect_x = lineVect->x; + const double lineVect_y = lineVect->y; + const double lineVect_z = lineVect->z; + const double v0_x = C->x - A_x; + const double v0_y = C->y - A_y; + const double v0_z = C->z - A_z; + const double v1_x = B->x - A_x; + const double v1_y = B->y - A_y; + const double v1_z = B->z - A_z; + const double _n_x = v0_y * v1_z - v0_z * v1_y; + const double _n_y = v0_z * v1_x - v0_x * v1_z; + const double _n_z = v0_x * v1_y - v0_y * v1_x; + const double k = ((A_x * _n_x + A_y * _n_y + A_z * _n_z) - (_n_x * linePoint_x + _n_y * linePoint_y + _n_z * linePoint_z)) / (_n_x * lineVect_x + _n_y * lineVect_y + _n_z * lineVect_z); + const double P_x = linePoint_x + lineVect_x * k; + const double P_y = linePoint_y + lineVect_y * k; + const double P_z = linePoint_z + lineVect_z * k; // Compute vectors - float v2_x = P_x - A_x; - float v2_y = P_y - A_y; - float v2_z = P_z - A_z; + const double v2_x = P_x - A_x; + const double v2_y = P_y - A_y; + const double v2_z = P_z - A_z; // Compute dot products - float dot00 = (v0_x * v0_x + v0_y * v0_y + v0_z * v0_z); - float dot01 = (v0_x * v1_x + v0_y * v1_y + v0_z * v1_z); - float dot02 = (v0_x * v2_x + v0_y * v2_y + v0_z * v2_z); - float dot11 = (v1_x * v1_x + v1_y * v1_y + v1_z * v1_z); - float dot12 = (v1_x * v2_x + v1_y * v2_y + v1_z * v2_z); + const double dot00 = (v0_x * v0_x + v0_y * v0_y + v0_z * v0_z); + const double dot01 = (v0_x * v1_x + v0_y * v1_y + v0_z * v1_z); + const double dot02 = (v0_x * v2_x + v0_y * v2_y + v0_z * v2_z); + const double dot11 = (v1_x * v1_x + v1_y * v1_y + v1_z * v1_z); + const double dot12 = (v1_x * v2_x + v1_y * v2_y + v1_z * v2_z); // Compute barycentric coordinates - float invDenom = 1.0 / (dot00 * dot11 - dot01 * dot01); - float u = (dot11 * dot02 - dot01 * dot12) * invDenom; - float v = (dot00 * dot12 - dot01 * dot02) * invDenom; + const double invDenom = 1.0 / (dot00 * dot11 - dot01 * dot01); + const double u = (dot11 * dot02 - dot01 * dot12) * invDenom; + const double v = (dot00 * dot12 - dot01 * dot02) * invDenom; P->x = P_x; P->y = P_y; P->z = P_z; From 432964cfd50401ad3cddbb4dca39eabeff7f9d8e Mon Sep 17 00:00:00 2001 From: Enguerrand DE SMET Date: Fri, 7 Aug 2020 13:05:52 +0200 Subject: [PATCH 20/26] [fuseCut] DelaunayGraphCut: handle ambigious case in rayIntersectTriangle - Add ambiguityEpsilon to handle cases where the direction test is ambiguous - Use double instead of float - renaming epsilon => marginEpsilon --- src/aliceVision/fuseCut/DelaunayGraphCut.cpp | 105 +++++++++++++++---- src/aliceVision/fuseCut/DelaunayGraphCut.hpp | 3 +- 2 files changed, 85 insertions(+), 23 deletions(-) diff --git a/src/aliceVision/fuseCut/DelaunayGraphCut.cpp b/src/aliceVision/fuseCut/DelaunayGraphCut.cpp index 5a5a53bfbc..baefc7d41d 100644 --- a/src/aliceVision/fuseCut/DelaunayGraphCut.cpp +++ b/src/aliceVision/fuseCut/DelaunayGraphCut.cpp @@ -1319,6 +1319,8 @@ DelaunayGraphCut::intersectNextGeom(const DelaunayGraphCut::GeometryIntersection const Point3d& dirVect, Point3d& intersectPt, const double epsilonFactor, const Point3d& lastIntersectPt) const { + GeometryIntersection bestMatch; + Point3d bestMatchIntersectPt; switch (inGeometry.type) { @@ -1334,10 +1336,21 @@ DelaunayGraphCut::intersectNextGeom(const DelaunayGraphCut::GeometryIntersection continue; const Facet intersectionFacet(tetrahedronIndex, i); + bool ambiguous = false; - const GeometryIntersection result = rayIntersectTriangle(originPt, dirVect, intersectionFacet, intersectPt, epsilonFactor, &lastIntersectPt); + const GeometryIntersection result = rayIntersectTriangle(originPt, dirVect, intersectionFacet, intersectPt, epsilonFactor, ambiguous, &lastIntersectPt); if (result.type != EGeometryType::None) - return result; + { + if (!ambiguous) + return result; + + // Retrieve the best intersected point (farthest from origin point) + if (bestMatch.type == EGeometryType::None || (originPt - intersectPt).size() > (originPt - bestMatchIntersectPt).size()) + { + bestMatchIntersectPt = intersectPt; + bestMatch = result; + } + } } } break; @@ -1354,10 +1367,21 @@ DelaunayGraphCut::intersectNextGeom(const DelaunayGraphCut::GeometryIntersection // Define the facet to intersect const Facet facet(adjCellIndex, localVertexIndex); + bool ambiguous = false; - const GeometryIntersection result = rayIntersectTriangle(originPt, dirVect, facet, intersectPt, epsilonFactor, &lastIntersectPt); + const GeometryIntersection result = rayIntersectTriangle(originPt, dirVect, facet, intersectPt, epsilonFactor, ambiguous, &lastIntersectPt); if (result.type != EGeometryType::None) - return result; + { + if (!ambiguous) + return result; + + // Retrieve the best intersected point (farthest from origin point) + if (bestMatch.type == EGeometryType::None || (originPt - intersectPt).size() > (originPt - bestMatchIntersectPt).size()) + { + bestMatchIntersectPt = intersectPt; + bestMatch = result; + } + } } } break; @@ -1379,9 +1403,20 @@ DelaunayGraphCut::intersectNextGeom(const DelaunayGraphCut::GeometryIntersection for (const Facet& facet : opositeFacets) { - const GeometryIntersection result = rayIntersectTriangle(originPt, dirVect, facet, intersectPt, epsilonFactor, &lastIntersectPt); + bool ambiguous = false; + const GeometryIntersection result = rayIntersectTriangle(originPt, dirVect, facet, intersectPt, epsilonFactor, ambiguous, &lastIntersectPt); if (result.type != EGeometryType::None) - return result; + { + if (!ambiguous) + return result; + + // Retrieve the best intersected point (farthest from origin point) + if (bestMatch.type == EGeometryType::None || (originPt - intersectPt).size() > (originPt - bestMatchIntersectPt).size()) + { + bestMatchIntersectPt = intersectPt; + bestMatch = result; + } + } } } } @@ -1391,15 +1426,19 @@ DelaunayGraphCut::intersectNextGeom(const DelaunayGraphCut::GeometryIntersection throw std::runtime_error("[intersectNextGeom] intersection with input none geometry should not happen."); break; } - return GeometryIntersection(); + + intersectPt = bestMatchIntersectPt; + return bestMatch; } DelaunayGraphCut::GeometryIntersection DelaunayGraphCut::rayIntersectTriangle(const Point3d& originPt, const Point3d& DirVec, const DelaunayGraphCut::Facet& facet, Point3d& intersectPt, - const double epsilonFactor, const Point3d* lastIntersectPt) const + const double epsilonFactor, bool &ambiguous, const Point3d* lastIntersectPt) const { + ambiguous = false; + const VertexIndex AvertexIndex = getVertexIndex(facet, 0); const VertexIndex BvertexIndex = getVertexIndex(facet, 1); const VertexIndex CvertexIndex = getVertexIndex(facet, 2); @@ -1408,42 +1447,64 @@ DelaunayGraphCut::GeometryIntersection DelaunayGraphCut::rayIntersectTriangle(co const Point3d* B = &_verticesCoords[BvertexIndex]; const Point3d* C = &_verticesCoords[CvertexIndex]; - const double epsilon = std::min(std::min((*A - *B).size(), (*B - *C).size()), (*A - *C).size()) * epsilonFactor; + const double ABSize = (*A - *B).size(); + const double BCSize = (*B - *C).size(); + const double ACSize = (*A - *C).size(); + + const double marginEpsilon = std::min(std::min(ABSize, BCSize), ACSize) * epsilonFactor; + const double ambiguityEpsilon = (ABSize + BCSize + ACSize) / 3.0 * 1.0e-5; + Point3d tempIntersectPt; const Point2d triangleUv = getLineTriangleIntersectBarycCoords(&tempIntersectPt, A, B, C, &originPt, &DirVec); if (!isnormal(tempIntersectPt.x) || !isnormal(tempIntersectPt.y) || !isnormal(tempIntersectPt.z)) { - ALICEVISION_LOG_DEBUG("Invalid/notNormal intersection point found during rayIntersectTriangle."); + // This is not suppose to happen in real life, we log a warning instead of raising an exeption if we face a border case + ALICEVISION_LOG_WARNING("Invalid/notNormal intersection point found during rayIntersectTriangle."); + return GeometryIntersection(); } - const float u = triangleUv.x; // A to C - const float v = triangleUv.y; // A to B + const double u = triangleUv.x; // A to C + const double v = triangleUv.y; // A to B // If we find invalid uv coordinate if (!std::isfinite(u) && !std::isfinite(v)) return GeometryIntersection(); - // Ouside the triangle with epsilon margin - if (u < -epsilon || v < -epsilon || (u + v) >(1 + epsilon)) + // Ouside the triangle with marginEpsilon margin + if (u < -marginEpsilon || v < -marginEpsilon || (u + v) >(1 + marginEpsilon)) return GeometryIntersection(); // In case intersectPt is provided, check if intersectPt is in front of lastIntersectionPt // in the DirVec direction to ensure that we are moving forward in the right direction - if (lastIntersectPt != nullptr && dot(DirVec, (tempIntersectPt - *lastIntersectPt).normalize()) <= 0) - return GeometryIntersection(); + if (lastIntersectPt != nullptr) { + const Point3d diff = tempIntersectPt - *lastIntersectPt; + const double dotValue = dot(DirVec, diff.normalize()); + + if (diff.size() < ambiguityEpsilon) + { + ambiguous = true; + } + else + { + if (dotValue <= 0) + { + return GeometryIntersection(); + } + } + } // Change intersection point only if tempIntersectionPt is in the right direction (mean we intersect something) intersectPt = tempIntersectPt; - if (v < epsilon) // along A C edge + if (v < marginEpsilon) // along A C edge { - if (u < epsilon) + if (u < marginEpsilon) { intersectPt = *A; return GeometryIntersection(AvertexIndex); // vertex A } - if (u > 1.0f - epsilon) + if (u > 1.0 - marginEpsilon) { intersectPt = *C; return GeometryIntersection(CvertexIndex); // vertex C @@ -1452,9 +1513,9 @@ DelaunayGraphCut::GeometryIntersection DelaunayGraphCut::rayIntersectTriangle(co return GeometryIntersection(Edge(AvertexIndex, CvertexIndex)); // edge AC } - if (u < epsilon) // along A B edge + if (u < marginEpsilon) // along A B edge { - if (v > 1.0f - epsilon) + if (v > 1.0 - marginEpsilon) { intersectPt = *B; return GeometryIntersection(BvertexIndex); // vertex B @@ -1463,7 +1524,7 @@ DelaunayGraphCut::GeometryIntersection DelaunayGraphCut::rayIntersectTriangle(co return GeometryIntersection(Edge(AvertexIndex, BvertexIndex)); // edge AB } - if (u + v > 1.0f - epsilon) + if (u + v > 1.0 - marginEpsilon) return GeometryIntersection(Edge(BvertexIndex, CvertexIndex)); // edge BC return GeometryIntersection(facet); diff --git a/src/aliceVision/fuseCut/DelaunayGraphCut.hpp b/src/aliceVision/fuseCut/DelaunayGraphCut.hpp index 361cf21653..0219594bd7 100644 --- a/src/aliceVision/fuseCut/DelaunayGraphCut.hpp +++ b/src/aliceVision/fuseCut/DelaunayGraphCut.hpp @@ -456,11 +456,12 @@ class DelaunayGraphCut * @param intersectPt a reference that will store the computed intersection point for the next intersecting geometry * @param epsilonFactor a multiplicative factor on the smaller side of the facet used to define the boundary when we * have to consider either a collision with an edge/vertex or a facet. + * @param ambiguous boolean used to know if our intersection is ambiguous or not * @param lastIntersectPt pointer to the last intersection point used to test the direction (if not nulllptr) * @return */ GeometryIntersection rayIntersectTriangle(const Point3d& originPt, const Point3d& DirVec, const Facet& facet, - Point3d& intersectPt, const double epsilonFactor, const Point3d* lastIntersectPt = nullptr) const; + Point3d& intersectPt, const double epsilonFactor, bool& ambiguous, const Point3d* lastIntersectPt = nullptr) const; float distFcn(float maxDist, float dist, float distFcnHeight) const; From 24117a9ce02a4b34e664fe2ec991c3fe37a6a6ca Mon Sep 17 00:00:00 2001 From: Enguerrand DE SMET Date: Fri, 7 Aug 2020 15:46:18 +0200 Subject: [PATCH 21/26] [fuseCut] DelaunayGraphCut : define an adapted value for marginEpsilon --- src/aliceVision/fuseCut/DelaunayGraphCut.cpp | 11 +++++++---- 1 file changed, 7 insertions(+), 4 deletions(-) diff --git a/src/aliceVision/fuseCut/DelaunayGraphCut.cpp b/src/aliceVision/fuseCut/DelaunayGraphCut.cpp index baefc7d41d..66d58b6d43 100644 --- a/src/aliceVision/fuseCut/DelaunayGraphCut.cpp +++ b/src/aliceVision/fuseCut/DelaunayGraphCut.cpp @@ -1769,6 +1769,7 @@ void DelaunayGraphCut::fillGraphPartPtRc(int& outTotalStepsFront, int& outTotalS bool fillOut, float distFcnHeight) // fixesSigma=true nPixelSizeBehind=2*spaceSteps allPoints=1 behind=0 fillOut=1 distFcnHeight=0 { const int maxint = 1000000; // std::numeric_limits::std::max() + const double marginEpsilonFactor = 1.0e-4; const Point3d& originPt = _verticesCoords[vertexIndex]; const float pixSize = mp->getCamPixelSize(originPt, cam); @@ -1804,7 +1805,7 @@ void DelaunayGraphCut::fillGraphPartPtRc(int& outTotalStepsFront, int& outTotalS #endif ++outTotalStepsFront; - geometry = intersectNextGeom(previousGeometry, originPt, dirVect, intersectPt, 1.0e-3, lastIntersectPt); + geometry = intersectNextGeom(previousGeometry, originPt, dirVect, intersectPt, marginEpsilonFactor, lastIntersectPt); if (geometry.type == EGeometryType::None) { @@ -1898,7 +1899,7 @@ void DelaunayGraphCut::fillGraphPartPtRc(int& outTotalStepsFront, int& outTotalS #endif ++outTotalStepsBehind; - geometry = intersectNextGeom(previousGeometry, originPt, dirVect, intersectPt, 1.0e-3, lastIntersectPt); + geometry = intersectNextGeom(previousGeometry, originPt, dirVect, intersectPt, marginEpsilonFactor, lastIntersectPt); if (geometry.type == EGeometryType::None) { @@ -2031,6 +2032,8 @@ void DelaunayGraphCut::forceTedgesByGradientIJCV(bool fixesSigma, float nPixelSi // c.out = c.gEdgeVisWeight[0] + c.gEdgeVisWeight[1] + c.gEdgeVisWeight[2] + c.gEdgeVisWeight[3]; } + const double marginEpsilonFactor = 1.0e-4; + // choose random order to prevent waiting const unsigned int seed = (unsigned int)mp->userParams.get("delaunaycut.seed", 0); const std::vector verticesRandIds = mvsUtils::createRandomArrayOfIntegers(_verticesAttr.size(), seed); @@ -2094,7 +2097,7 @@ void DelaunayGraphCut::forceTedgesByGradientIJCV(bool fixesSigma, float nPixelSi #endif ++totalStepsFront; - geometry = intersectNextGeom(previousGeometry, originPt, dirVect, intersectPt, 1.0e-3, lastIntersectPt); + geometry = intersectNextGeom(previousGeometry, originPt, dirVect, intersectPt, marginEpsilonFactor, lastIntersectPt); if (geometry.type == EGeometryType::None) { @@ -2172,7 +2175,7 @@ void DelaunayGraphCut::forceTedgesByGradientIJCV(bool fixesSigma, float nPixelSi #endif ++totalStepsBehind; - geometry = intersectNextGeom(previousGeometry, originPt, dirVect, intersectPt, 1.0e-3, lastIntersectPt); + geometry = intersectNextGeom(previousGeometry, originPt, dirVect, intersectPt, marginEpsilonFactor, lastIntersectPt); if(geometry.type == EGeometryType::None) { From 87903cb31ca2ae4adff90d0a1304241b1a6dca04 Mon Sep 17 00:00:00 2001 From: Enguerrand DE SMET Date: Fri, 7 Aug 2020 15:58:59 +0200 Subject: [PATCH 22/26] [FuseCut] DelaunayGraphCut: cleaner debug messages --- src/aliceVision/fuseCut/DelaunayGraphCut.cpp | 50 ++++++++++++++------ 1 file changed, 35 insertions(+), 15 deletions(-) diff --git a/src/aliceVision/fuseCut/DelaunayGraphCut.cpp b/src/aliceVision/fuseCut/DelaunayGraphCut.cpp index 66d58b6d43..a7abbfa233 100644 --- a/src/aliceVision/fuseCut/DelaunayGraphCut.cpp +++ b/src/aliceVision/fuseCut/DelaunayGraphCut.cpp @@ -1812,11 +1812,23 @@ void DelaunayGraphCut::fillGraphPartPtRc(int& outTotalStepsFront, int& outTotalS #ifdef ALICEVISION_DEBUG_VOTE // exportBackPropagationMesh("fillGraph_ToCam_typeNone", history.geometries, originPt, mp->CArr[cam]); #endif - // throw std::runtime_error("[Error]: fillGraphPartPtRc toTheCam, cause: geometry cannot be found."); - ALICEVISION_LOG_TRACE("[Error]: fillGraphPartPtRc toTheCam, cause: geometry cannot be found."); + ALICEVISION_LOG_DEBUG("[Error]: fillGraph(toTheCam) cause: geometry cannot be found."); break; } +#ifdef ALICEVISION_DEBUG_VOTE + { + const auto end = history.geometries.end(); + auto it = std::find(history.geometries.begin(), end, geometry); + if (it != end) + { + // exportBackPropagationMesh("fillGraph_ToCam_alreadyIntersected", history.geometries, originPt, mp->CArr[cam]); + ALICEVISION_LOG_DEBUG("[Error]: fillGraph(toTheCam) cause: intersected geometry has already been intersected."); + break; + } + } +#endif + if (geometry.type == EGeometryType::Facet) { ++outFrontCount.facets; @@ -1839,8 +1851,7 @@ void DelaunayGraphCut::fillGraphPartPtRc(int& outTotalStepsFront, int& outTotalS #ifdef ALICEVISION_DEBUG_VOTE // exportBackPropagationMesh("fillGraph_ToCam_invalidMirorFacet", history.geometries, originPt, mp->CArr[cam]); #endif - // throw std::runtime_error("[Error]: fillGraphPartPtRc toTheCam, cause: invalidOrInfinite miror facet."); - ALICEVISION_LOG_TRACE("[Error]: fillGraphPartPtRc toTheCam, cause: invalidOrInfinite miror facet."); + ALICEVISION_LOG_DEBUG("[Error]: fillGraph(toTheCam) cause: invalidOrInfinite miror facet."); break; } lastIntersectedFacet = mFacet; @@ -1909,8 +1920,7 @@ void DelaunayGraphCut::fillGraphPartPtRc(int& outTotalStepsFront, int& outTotalS #ifdef ALICEVISION_DEBUG_VOTE // exportBackPropagationMesh("fillGraph_behindThePoint_NoneButPreviousIsFacet", history.geometries, originPt, mp->CArr[cam]); #endif - // throw std::runtime_error("[DelaunayGraphCut::fillGraphPartPtRc] first loop: Geometry cannot be found."); - ALICEVISION_LOG_TRACE("[Error]: fillGraphPartPtRc behindThePoint, cause: None geometry but previous is Facet."); + ALICEVISION_LOG_DEBUG("[Error]: fillGraph(behindThePoint) cause: None geometry but previous is Facet."); } // Break if we reach the end of the tetrahedralization volume break; @@ -2102,13 +2112,25 @@ void DelaunayGraphCut::forceTedgesByGradientIJCV(bool fixesSigma, float nPixelSi if (geometry.type == EGeometryType::None) { #ifdef ALICEVISION_DEBUG_VOTE - // exportBackPropagationMesh("forceTedgesByGradientIJCV_ToCam_typeNone", history.geometries, originPt, mp->CArr[cam]); + // exportBackPropagationMesh("forceTedges_ToCam_typeNone", history.geometries, originPt, mp->CArr[cam]); #endif - // throw std::runtime_error("[DelaunayGraphCut::fillGraphPartPtRc] first loop: Geometry cannot be found."); - ALICEVISION_LOG_TRACE("[Error]: forceTedgesByGradientIJCV toTheCam, cause: geometry cannot be found."); + ALICEVISION_LOG_DEBUG("[Error]: forceTedges(toTheCam) cause: geometry cannot be found."); break; } +#ifdef ALICEVISION_DEBUG_VOTE + { + const auto end = history.geometries.end(); + auto it = std::find(history.geometries.begin(), end, geometry); + if (it != end) + { + // exportBackPropagationMesh("forceTedges_ToCam_alreadyIntersected", history.geometries, originPt, mp->CArr[cam]); + ALICEVISION_LOG_DEBUG("[Error]: forceTedges(toTheCam) cause: intersected geometry has already been intersected."); + break; + } + } +#endif + if (geometry.type == EGeometryType::Facet) { ++totalGeometriesIntersectedFrontCount.facets; @@ -2129,10 +2151,9 @@ void DelaunayGraphCut::forceTedgesByGradientIJCV(bool fixesSigma, float nPixelSi if (isInvalidOrInfiniteCell(mFacet.cellIndex)) { #ifdef ALICEVISION_DEBUG_VOTE - // exportBackPropagationMesh("forceTedgesByGradientIJCV_ToCam_invalidMirorFacet", history.geometries, originPt, mp->CArr[cam]); + // exportBackPropagationMesh("forceTedges_ToCam_invalidMirorFacet", history.geometries, originPt, mp->CArr[cam]); #endif - // throw std::runtime_error("[DelaunayGraphCut::fillGraphPartPtRc] first loop: Geometry cannot be found."); - ALICEVISION_LOG_TRACE("[Error]: forceTedgesByGradientIJCV toTheCam, cause: invalidOrInfinite miror facet."); + ALICEVISION_LOG_DEBUG("[Error]: forceTedges(toTheCam) cause: invalidOrInfinite miror facet."); break; } geometry.facet = mFacet; @@ -2183,10 +2204,9 @@ void DelaunayGraphCut::forceTedgesByGradientIJCV(bool fixesSigma, float nPixelSi if (previousGeometry.type == EGeometryType::Facet) { #ifdef ALICEVISION_DEBUG_VOTE - // exportBackPropagationMesh("forceTedgesByGradientIJCV_behindThePoint_typeNone", history.geometries, originPt, mp->CArr[cam]); + // exportBackPropagationMesh("forceTedges_behindThePoint_NoneButPreviousIsFacet", history.geometries, originPt, mp->CArr[cam]); #endif - // throw std::runtime_error("[Error]: forceTedgesByGradientIJCV behindThePoint, cause: geometry cannot be found."); - ALICEVISION_LOG_TRACE("[Error]: forceTedgesByGradientIJCV behindThePoint, cause: geometry cannot be found."); + ALICEVISION_LOG_DEBUG("[Error]: forceTedges(behindThePoint) cause: None geometry but previous is Facet."); } // Break if we reach the end of the tetrahedralization volume break; From c620efc7ac1ef3b97a56637a329ec38f63b072ce Mon Sep 17 00:00:00 2001 From: Enguerrand DE SMET Date: Fri, 7 Aug 2020 18:08:07 +0200 Subject: [PATCH 23/26] [FuseCut] DelaunayGraphCut: add GeometryIntersection operators --- src/aliceVision/fuseCut/DelaunayGraphCut.cpp | 20 ++++++++++++++++++++ src/aliceVision/fuseCut/DelaunayGraphCut.hpp | 7 +++++++ 2 files changed, 27 insertions(+) diff --git a/src/aliceVision/fuseCut/DelaunayGraphCut.cpp b/src/aliceVision/fuseCut/DelaunayGraphCut.cpp index a7abbfa233..b75fb1d02b 100644 --- a/src/aliceVision/fuseCut/DelaunayGraphCut.cpp +++ b/src/aliceVision/fuseCut/DelaunayGraphCut.cpp @@ -3469,5 +3469,25 @@ void DelaunayGraphCut::leaveLargestFullSegmentOnly() ALICEVISION_LOG_DEBUG("Largest full segment only done."); } +std::ostream& operator<<(std::ostream& stream, const DelaunayGraphCut::EGeometryType type) +{ + switch (type) + { + case DelaunayGraphCut::EGeometryType::Facet: + stream << "Facet"; + break; + case DelaunayGraphCut::EGeometryType::Vertex: + stream << "Vertex"; + break; + case DelaunayGraphCut::EGeometryType::Edge: + stream << "Edge"; + break; + case DelaunayGraphCut::EGeometryType::None: + stream << "None"; + break; + } + return stream; +} + } // namespace fuseCut } // namespace aliceVision diff --git a/src/aliceVision/fuseCut/DelaunayGraphCut.hpp b/src/aliceVision/fuseCut/DelaunayGraphCut.hpp index 0219594bd7..cb2f2aa661 100644 --- a/src/aliceVision/fuseCut/DelaunayGraphCut.hpp +++ b/src/aliceVision/fuseCut/DelaunayGraphCut.hpp @@ -148,6 +148,10 @@ class DelaunayGraphCut return true; } } + bool operator!=(const GeometryIntersection& g) const + { + return !(*this == g); + } }; struct IntersectionHistory @@ -517,5 +521,8 @@ class DelaunayGraphCut void writeScoreInCsv(const std::string& filePath, const size_t& sizeLimit = 1000); }; + +std::ostream& operator<<(std::ostream& stream, const DelaunayGraphCut::EGeometryType type); + } // namespace fuseCut } // namespace aliceVision From 2760cbe44b3154695442ca41dfa1aee8339e0c1b Mon Sep 17 00:00:00 2001 From: Enguerrand DE SMET Date: Fri, 7 Aug 2020 19:02:07 +0200 Subject: [PATCH 24/26] [software] meshing: add log info for boundingBox size --- src/software/pipeline/main_meshing.cpp | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/src/software/pipeline/main_meshing.cpp b/src/software/pipeline/main_meshing.cpp index 95a8288471..b6b9c9d9d1 100644 --- a/src/software/pipeline/main_meshing.cpp +++ b/src/software/pipeline/main_meshing.cpp @@ -460,6 +460,14 @@ int aliceVision_main(int argc, char* argv[]) else fs.divideSpaceFromSfM(sfmData, &hexah[0], estimateSpaceMinObservations, estimateSpaceMinObservationAngle); + { + const double length = hexah[0].x - hexah[1].x; + const double width = hexah[0].y - hexah[3].y; + const double height = hexah[0].z - hexah[4].z; + + ALICEVISION_LOG_INFO("bounding Box : length: " << length << ", width: " << width << ", height: " << height); + } + StaticVector cams; if(meshingFromDepthMaps) { From 5b2930a4a8a71563dcb34ef99b1ed770b14c137d Mon Sep 17 00:00:00 2001 From: Enguerrand DE SMET Date: Fri, 7 Aug 2020 19:03:22 +0200 Subject: [PATCH 25/26] [FuseCut] DelaunayGraphCut: function documentation - forceTedgesByGradientIJCV --- src/aliceVision/fuseCut/DelaunayGraphCut.hpp | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/src/aliceVision/fuseCut/DelaunayGraphCut.hpp b/src/aliceVision/fuseCut/DelaunayGraphCut.hpp index cb2f2aa661..fd8276b151 100644 --- a/src/aliceVision/fuseCut/DelaunayGraphCut.hpp +++ b/src/aliceVision/fuseCut/DelaunayGraphCut.hpp @@ -483,6 +483,12 @@ class DelaunayGraphCut bool fixesSigma, float nPixelSizeBehind, bool fillOut, float distFcnHeight); + /** + * @brief Estimate the cells property "on" based on the analysis of the visibility of neigbouring cells. + * + * @param fixesSigma Use constant sigma for all points or take into account the pixelSize + * @param nPixelSizeBehind Used to define the surface margin + */ void forceTedgesByGradientIJCV(bool fixesSigma, float nPixelSizeBehind); int setIsOnSurface(); From 5e545ae1fa6eead4d046f2dd79e51d610666d2eb Mon Sep 17 00:00:00 2001 From: Fabien Castan Date: Wed, 20 Jan 2021 20:58:56 +0100 Subject: [PATCH 26/26] [fuseCut] build fix: missing return of non-void function --- src/aliceVision/fuseCut/DelaunayGraphCut.cpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/aliceVision/fuseCut/DelaunayGraphCut.cpp b/src/aliceVision/fuseCut/DelaunayGraphCut.cpp index b75fb1d02b..8ef8e09cbe 100644 --- a/src/aliceVision/fuseCut/DelaunayGraphCut.cpp +++ b/src/aliceVision/fuseCut/DelaunayGraphCut.cpp @@ -474,8 +474,9 @@ std::vector DelaunayGraphCut::getNeighboringCellsBy case EGeometryType::Facet: return getNeighboringCellsByFacet(g.facet); case EGeometryType::None: - throw std::runtime_error("[error] getNeighboringCellsByGeometry: an undefined/None geometry has no neighboring cells."); + break; } + throw std::runtime_error("[error] getNeighboringCellsByGeometry: an undefined/None geometry has no neighboring cells."); } std::vector DelaunayGraphCut::getNeighboringCellsByFacet(const Facet& f) const