From e67cca0048acc427cc5929e85c5ecb350a205f0a Mon Sep 17 00:00:00 2001 From: Fabien Castan Date: Thu, 18 Feb 2021 13:15:00 +0100 Subject: [PATCH] [fuseCut] option to densify vertex candidates and more robust filtering of helper points - option to densify vertex candidates - more robust filtering of helper points - debug export of tetrahedralization as dense mesh - more argument on meshing cmd line --- src/aliceVision/fuseCut/DelaunayGraphCut.cpp | 367 ++++++++++++------- src/aliceVision/fuseCut/DelaunayGraphCut.hpp | 35 +- src/software/pipeline/main_meshing.cpp | 47 ++- 3 files changed, 305 insertions(+), 144 deletions(-) diff --git a/src/aliceVision/fuseCut/DelaunayGraphCut.cpp b/src/aliceVision/fuseCut/DelaunayGraphCut.cpp index c3d8826f19..91eb856b5a 100644 --- a/src/aliceVision/fuseCut/DelaunayGraphCut.cpp +++ b/src/aliceVision/fuseCut/DelaunayGraphCut.cpp @@ -824,12 +824,59 @@ void DelaunayGraphCut::addPointsToPreventSingularities(const Point3d voxel[8], f ALICEVISION_LOG_WARNING("Add " << addedPoints << " points to prevent singularities"); } -void DelaunayGraphCut::addHelperPoints(int nGridHelperVolumePointsDim, const Point3d voxel[8], float minDist) +void DelaunayGraphCut::densifyWithHelperPoints(int nbFront, int nbBack, double scale) { - if(nGridHelperVolumePointsDim <= 0) + if(nbFront <= 0 && nbBack <= 0) return; - int ns = nGridHelperVolumePointsDim; + const std::size_t nbInputVertices = _verticesCoords.size(); + + std::vector newHelperPoints; + newHelperPoints.reserve((nbFront + nbBack) * nbInputVertices); + + for(std::size_t vi = 0; vi < nbInputVertices; ++vi) + { + const Point3d& v = _verticesCoords[vi]; + const GC_vertexInfo& vAttr = _verticesAttr[vi]; + + if(vAttr.cams.empty() || vAttr.pixSize <= std::numeric_limits::epsilon()) + continue; + + Point3d mainCamDir; + for(int camId: vAttr.cams) + { + const Point3d& cam = mp->CArr[camId]; + const Point3d d = (cam - v).normalize(); + mainCamDir += d; + } + mainCamDir /= double(vAttr.cams.size()); + mainCamDir = mainCamDir.normalize() * vAttr.pixSize; + + for(int iFront = 1; iFront < nbFront + 1; ++iFront) + newHelperPoints.push_back(v + mainCamDir * iFront * scale); + for(int iBack = 1; iBack < nbBack + 1; ++iBack) + newHelperPoints.push_back(v - mainCamDir * iBack * scale); + } + _verticesCoords.resize(nbInputVertices + newHelperPoints.size()); + _verticesAttr.resize(nbInputVertices + newHelperPoints.size()); + for(std::size_t vi = 0; vi < newHelperPoints.size(); ++vi) + { + _verticesCoords[nbInputVertices + vi] = newHelperPoints[vi]; + // GC_vertexInfo& vAttr = _verticesAttr[nbInputVertices + vi]; + // Keep vertexInfo with default/empty values, so they will be removed at the end as other helper points + // vAttr.nrc = 0; + } + + ALICEVISION_LOG_WARNING("Densify the " << nbInputVertices << " vertices with " << newHelperPoints.size() + << " new helper points."); +} + +void DelaunayGraphCut::addGridHelperPoints(int helperPointsGridSize, const Point3d voxel[8], float minDist) +{ + if(helperPointsGridSize <= 0) + return; + + int ns = helperPointsGridSize; float md = 1.0f / 500.0f; Point3d vx = (voxel[1] - voxel[0]); Point3d vy = (voxel[3] - voxel[0]); @@ -879,6 +926,9 @@ void DelaunayGraphCut::addHelperPoints(int nGridHelperVolumePointsDim, const Poi void DelaunayGraphCut::addMaskHelperPoints(const Point3d voxel[8], const StaticVector& cams, const FuseParams& params) { + if(params.maskHelperPointsWeight <= 0.0) + return; + ALICEVISION_LOG_INFO("Add Mask Helper Points."); Point3d inflatedVoxel[8]; @@ -889,17 +939,9 @@ void DelaunayGraphCut::addMaskHelperPoints(const Point3d voxel[8], const StaticV { nbPixels += imgParams.size; } - int step = std::floor(std::sqrt(double(nbPixels) / double(params.maxInputPoints))); - step = std::max(step, params.minStep); - std::size_t realMaxVertices = 0; - std::vector startIndex(mp->getNbCameras(), 0); - for(int i = 0; i < mp->getNbCameras(); ++i) - { - const auto& imgParams = mp->getImageParams(i); - startIndex[i] = realMaxVertices; - realMaxVertices += std::ceil(imgParams.width / step) * std::ceil(imgParams.height / step); - } - + // int step = std::floor(std::sqrt(double(nbPixels) / double(params.maxInputPoints))); + // step = std::max(step, params.minStep); + const int step = 1; int nbAddedPoints = 0; ALICEVISION_LOG_INFO("Load depth maps and add points."); @@ -925,7 +967,6 @@ void DelaunayGraphCut::addMaskHelperPoints(const Point3d voxel[8], const StaticV { for(int sx = 0; sx < sxMax; ++sx) { - int index = startIndex[c] + sy * sxMax + sx; float bestScore = 0; int bestX = 0; @@ -936,12 +977,13 @@ void DelaunayGraphCut::addMaskHelperPoints(const Point3d voxel[8], const StaticV { const std::size_t index = y * width + x; const float depth = depthMap[index]; - - if(depth >= -1.0f) + + // -2 means that the pixels should be masked-out with mask helper points + if(depth > -1.5f) continue; int nbValidDepth = 0; - const int kernelSize = 1; + const int kernelSize = params.maskBorderSize; for(int ly = std::max(y - kernelSize, 0), lyMax = std::min(y + kernelSize, height - 1); ly < lyMax; ++ly) { @@ -983,7 +1025,7 @@ void DelaunayGraphCut::addMaskHelperPoints(const Point3d voxel[8], const StaticV } } GC_vertexInfo newv; - newv.nrc = 1; + newv.nrc = params.maskHelperPointsWeight; newv.pixSize = 0.0f; newv.cams.push_back_distinct(c); @@ -1345,7 +1387,7 @@ void DelaunayGraphCut::computeVerticesSegSize(std::vector& out_segments, { int rc = v.getCamera(0); - // go thru all neighbour points + // go through all the neighbouring points GEO::vector adjVertices; _tetrahedralization->get_neighbors(vi, adjVertices); @@ -1795,8 +1837,9 @@ float DelaunayGraphCut::weightFcn(float nrc, bool labatutWeights, int /*ncams*/ return weight; } -void DelaunayGraphCut::fillGraph(double nPixelSizeBehind, - bool labatutWeights, bool fillOut, float distFcnHeight) // nPixelSizeBehind=2*spaceSteps allPoints=1 behind=0 labatutWeights=0 fillOut=1 distFcnHeight=0 +void DelaunayGraphCut::fillGraph(double nPixelSizeBehind, bool labatutWeights, bool fillOut, float distFcnHeight, + float fullWeight) // nPixelSizeBehind=2*spaceSteps allPoints=1 behind=0 + // labatutWeights=0 fillOut=1 distFcnHeight=0 { ALICEVISION_LOG_INFO("Computing s-t graph weights."); long t1 = clock(); @@ -1847,7 +1890,7 @@ void DelaunayGraphCut::fillGraph(double nPixelSizeBehind, const int vertexIndex = verticesRandIds[i]; const GC_vertexInfo& v = _verticesAttr[vertexIndex]; - if(v.isReal() && (v.nrc > 0)) + if(v.isReal()) { ++totalIsRealNrc; // "weight" is called alpha(p) in the paper @@ -1862,7 +1905,9 @@ void DelaunayGraphCut::fillGraph(double nPixelSizeBehind, int stepsBehind = 0; GeometriesCount geometriesIntersectedFrontCount; GeometriesCount geometriesIntersectedBehindCount; - fillGraphPartPtRc(stepsFront, stepsBehind, geometriesIntersectedFrontCount, geometriesIntersectedBehindCount, vertexIndex, v.cams[c], weight, nPixelSizeBehind, + fillGraphPartPtRc(stepsFront, stepsBehind, geometriesIntersectedFrontCount, + geometriesIntersectedBehindCount, vertexIndex, v.cams[c], weight, fullWeight, + nPixelSizeBehind, fillOut, distFcnHeight); totalStepsFront += stepsFront; @@ -1905,8 +1950,9 @@ void DelaunayGraphCut::fillGraph(double nPixelSizeBehind, mvsUtils::printfElapsedTime(t1, "s-t graph weights computed : "); } -void DelaunayGraphCut::fillGraphPartPtRc(int& outTotalStepsFront, int& outTotalStepsBehind, GeometriesCount& outFrontCount, GeometriesCount& outBehindCount, int vertexIndex, int cam, - float weight, double nPixelSizeBehind, +void DelaunayGraphCut::fillGraphPartPtRc( + int& outTotalStepsFront, int& outTotalStepsBehind, GeometriesCount& outFrontCount, GeometriesCount& outBehindCount, + int vertexIndex, int cam, float weight, float fullWeight, double nPixelSizeBehind, bool fillOut, float distFcnHeight) // nPixelSizeBehind=2*spaceSteps allPoints=1 behind=0 fillOut=1 distFcnHeight=0 { const int maxint = 1000000; // std::numeric_limits::std::max() @@ -1919,7 +1965,7 @@ void DelaunayGraphCut::fillGraphPartPtRc(int& outTotalStepsFront, int& outTotalS assert(cam >= 0); assert(cam < mp->ncams); - if(fillOut) + if(fillOut) // EMPTY part { // Initialisation GeometryIntersection geometry(vertexIndex); // Starting on global vertex index @@ -2079,8 +2125,9 @@ void DelaunayGraphCut::fillGraphPartPtRc(int& outTotalStepsFront, int& outTotalS // } } - if(pixSize > 0.0f) // fillIn + if(pixSize > 0.0f) // fillIn FULL part { + const float fWeight = fullWeight * weight; // Initialisation GeometryIntersection geometry(vertexIndex); // Starting on global vertex index Point3d intersectPt = originPt; @@ -2135,13 +2182,13 @@ void DelaunayGraphCut::fillGraphPartPtRc(int& outTotalStepsFront, int& outTotalS if (firstIteration) { #pragma OMP_ATOMIC_UPDATE - _cellsAttr[geometry.facet.cellIndex].on += weight; + _cellsAttr[geometry.facet.cellIndex].on += fWeight; firstIteration = false; } { #pragma OMP_ATOMIC_UPDATE - _cellsAttr[geometry.facet.cellIndex].fullnessScore += weight; + _cellsAttr[geometry.facet.cellIndex].fullnessScore += fWeight; } // Take the mirror facet to iterate over the next cell @@ -2157,7 +2204,8 @@ void DelaunayGraphCut::fillGraphPartPtRc(int& outTotalStepsFront, int& outTotalS { 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] += + fWeight * dist; } if(previousGeometry.type == EGeometryType::Facet && outBehindCount.facets > 1000) { @@ -2186,7 +2234,7 @@ void DelaunayGraphCut::fillGraphPartPtRc(int& outTotalStepsFront, int& outTotalS for (const CellIndex& ci : neighboringCells) { #pragma OMP_ATOMIC_UPDATE - _cellsAttr[neighboringCells[0]].on += weight; + _cellsAttr[neighboringCells[0]].on += fWeight; } firstIteration = false; } @@ -2196,7 +2244,7 @@ void DelaunayGraphCut::fillGraphPartPtRc(int& outTotalStepsFront, int& outTotalS if (previousGeometry.type == EGeometryType::Facet) { #pragma OMP_ATOMIC_UPDATE - _cellsAttr[previousGeometry.facet.cellIndex].fullnessScore += weight; + _cellsAttr[previousGeometry.facet.cellIndex].fullnessScore += fWeight; } if (geometry.type == EGeometryType::Vertex) @@ -2225,7 +2273,7 @@ void DelaunayGraphCut::fillGraphPartPtRc(int& outTotalStepsFront, int& outTotalS if (lastIntersectedFacet.cellIndex != GEO::NO_CELL) { #pragma OMP_ATOMIC_UPDATE - _cellsAttr[lastIntersectedFacet.cellIndex].cellTWeight += weight; + _cellsAttr[lastIntersectedFacet.cellIndex].cellTWeight += fWeight; } } } @@ -3008,6 +3056,10 @@ void DelaunayGraphCut::createDensePointCloud(const Point3d hexah[8], const Stati ALICEVISION_LOG_INFO("Creating dense point cloud."); float minDist = hexah ? (hexah[0] - hexah[1]).size() / 1000.0f : 0.00001f; + const int helperPointsGridSize = mp->userParams.get("LargeScale.helperPointsGridSize", 10); + const int densifyNbFront = mp->userParams.get("LargeScale.densifyNbFront", 0); + const int densifyNbBack = mp->userParams.get("LargeScale.densifyNbBack", 0); + const double densifyScale = mp->userParams.get("LargeScale.densifyScale", 1.0); // add points from depth maps if(depthMapsFuseParams != nullptr) @@ -3017,22 +3069,22 @@ void DelaunayGraphCut::createDensePointCloud(const Point3d hexah[8], const Stati if(sfmData != nullptr) addPointsFromSfM(hexah, cams, *sfmData); - const int nGridHelperVolumePointsDim = mp->userParams.get("LargeScale.nGridHelperVolumePointsDim", 10); - // add points for cam centers addPointsFromCameraCenters(cams, minDist); // add 6 points to prevent singularities addPointsToPreventSingularities(hexah, minDist); + densifyWithHelperPoints(densifyNbFront, densifyNbBack, densifyScale); + // add volume points to prevent singularities { Point3d hexahExt[8]; mvsUtils::inflateHexahedron(hexah, hexahExt, 1.1); - addHelperPoints(nGridHelperVolumePointsDim, hexahExt, minDist); + addGridHelperPoints(helperPointsGridSize, hexahExt, minDist); // add point for shape from silhouette - if(mp->userParams.get("delaunaycut.addMaskHelperPoints", false)) + if(depthMapsFuseParams != nullptr) addMaskHelperPoints(hexahExt, cams, *depthMapsFuseParams); } @@ -3042,7 +3094,9 @@ void DelaunayGraphCut::createDensePointCloud(const Point3d hexah[8], const Stati ALICEVISION_LOG_WARNING("Final dense point cloud: " << _verticesCoords.size() << " points."); } -void DelaunayGraphCut::createGraphCut(const Point3d hexah[8], const StaticVector& cams, const std::string& folderName, const std::string& tmpCamsPtsFolderName, bool removeSmallSegments) +void DelaunayGraphCut::createGraphCut(const Point3d hexah[8], const StaticVector& cams, + const std::string& folderName, const std::string& tmpCamsPtsFolderName, + bool removeSmallSegments, bool exportDebugTetrahedralization) { // Create tetrahedralization computeDelaunay(); @@ -3058,32 +3112,8 @@ void DelaunayGraphCut::createGraphCut(const Point3d hexah[8], const StaticVector voteFullEmptyScore(cams, folderName); - if(false) - { - std::unique_ptr mesh(createTetrahedralMesh(true, 0.99f, [](const GC_cellInfo& c) { return c.cellSWeight; })); - mesh->saveToObj(folderName + "tetrahedralMesh_SWeight.obj"); - } - - if(false) - { - { - std::unique_ptr meshf(createTetrahedralMesh(false, 0.99f, [](const fuseCut::GC_cellInfo& c) { return c.emptinessScore; })); - meshf->saveToObj(folderName + "tetrahedralMesh_emptiness.obj"); - } - { - std::unique_ptr meshf(createTetrahedralMesh(false, 0.99f, [](const fuseCut::GC_cellInfo& c) { return c.fullnessScore; })); - meshf->saveToObj(folderName + "tetrahedralMesh_fullness.obj"); - } - /* - { - std::unique_ptr meshf(createTetrahedralMesh(false, 0.9f, [](const fuseCut::GC_cellInfo& c) { return c.fullnessScore - c.emptinessScore; })); - meshf->saveToObj(folderName + "tetrahedralMesh_emptiness-fullness.obj"); - } - { - std::unique_ptr meshf(createTetrahedralMesh(false, 0.9f, [](const fuseCut::GC_cellInfo& c) { return c.cellSWeight - c.cellTWeight; })); - meshf->saveToObj(folderName + "tetrahedralMesh_SmT.obj"); - }*/ - } + if(exportDebugTetrahedralization) + exportFullScoreMeshs(folderName, ""); maxflow(); } @@ -3208,12 +3238,13 @@ void DelaunayGraphCut::voteFullEmptyScore(const StaticVector& cams, const s long t1; // TODO FACA: nPixelSizeBehind 2 or 4 by default? - const double nPixelSizeBehind = mp->userParams.get("delaunaycut.nPixelSizeBehind", 4.0f); // sigma value + const double nPixelSizeBehind = mp->userParams.get("delaunaycut.nPixelSizeBehind", 4.0); // sigma value + const float fullWeight = float(mp->userParams.get("delaunaycut.fullWeight", 1.0)); ALICEVISION_LOG_INFO("nPixelSizeBehind: " << nPixelSizeBehind); // 0 for distFcn equals 1 all the time - const float distFcnHeight = (float)mp->userParams.get("delaunaycut.distFcnHeight", 0.0f); + const float distFcnHeight = (float)mp->userParams.get("delaunaycut.distFcnHeight", 0.0); const bool labatutCFG09 = mp->userParams.get("global.LabatutCFG09", false); // jancosekIJCV: "Exploiting Visibility Information in Surface Reconstruction to Preserve Weakly Supported Surfaces", Michal Jancosek and Tomas Pajdla, 2014 @@ -3226,7 +3257,7 @@ void DelaunayGraphCut::voteFullEmptyScore(const StaticVector& cams, const s displayCellsStats(); // compute weights on edge between tetrahedra - fillGraph(nPixelSizeBehind, false, true, distFcnHeight); + fillGraph(nPixelSizeBehind, false, true, distFcnHeight, fullWeight); displayCellsStats(); @@ -3257,14 +3288,112 @@ void DelaunayGraphCut::voteFullEmptyScore(const StaticVector& cams, const s if(labatutCFG09) { ALICEVISION_LOG_INFO("Labatut CFG 2009 method:"); - fillGraph(nPixelSizeBehind, true, true, distFcnHeight); + fillGraph(nPixelSizeBehind, true, true, distFcnHeight, fullWeight); if(saveTemporaryBinFiles) saveDhInfo(folderName + "delaunayTriangulationInfoInit.bin"); } } -mesh::Mesh* DelaunayGraphCut::createMesh(bool filterHelperPointsTriangles) +void DelaunayGraphCut::filterLargeHelperPoints(std::vector& out_reliableVertices, + const std::vector& vertexIsOnSurface, int maxSegSize) +{ + ALICEVISION_LOG_DEBUG("DelaunayGraphCut::filterLargeHelperPoints"); + out_reliableVertices.clear(); + + // Do not filter helper points if maxSegSize is negative/infinit + if(maxSegSize < 0) + return; + + out_reliableVertices.resize(_verticesAttr.size(), true); + + std::vector edges; + edges.reserve(_verticesAttr.size()); + + // Create edges for all connected helper points + for(VertexIndex vi = 0; vi < _verticesAttr.size(); ++vi) + { + const GC_vertexInfo& v = _verticesAttr[vi]; + const Point3d& p = _verticesCoords[vi]; + if(v.nrc <= 0 && vertexIsOnSurface[vi]) + { + // go through all the neighbouring points + GEO::vector adjVertices; + _tetrahedralization->get_neighbors(vi, adjVertices); + + for(VertexIndex nvi : adjVertices) + { + // ignore itself + if(vi == nvi) + continue; + // avoid duplicates + if(vi < nvi) + continue; + // ignore points not on the surface + if(!vertexIsOnSurface[vi]) + continue; + // ignore valid vertices + const GC_vertexInfo& nv = _verticesAttr[nvi]; + if(nv.nrc > 0) + continue; + // Declare a new edge between 2 helper points both on the selected surface + edges.emplace_back(vi, nvi); + } + } + } + + Universe u(_verticesAttr.size()); + + if (maxSegSize > 0) + { + int s = (int)edges.size(); + // Fuse all edges collected to be merged + for(int i = 0; i < s; i++) + { + int a = u.find(edges[i].x); + int b = u.find(edges[i].y); + if(a != b) + { + u.join(a, b); + } + } + } + + // Last loop over vertices to update segId + for(int vi = 0; vi < _verticesAttr.size(); ++vi) + { + if(!vertexIsOnSurface[vi]) + { + // Point is not on the surface + out_reliableVertices[vi] = false; + continue; + } + const GC_vertexInfo& v = _verticesAttr[vi]; + if(v.nrc > 0) + { + // This is not a helper point, so it is reliable. + out_reliableVertices[vi] = true; + continue; + } + + if(maxSegSize > 0) + { + // It is an helper point, so it is reliable only if it is a small group. + const int sigId = u.find(vi); + const int segSize = u.elts[sigId].size; + out_reliableVertices[vi] = segSize <= maxSegSize; + } + else + { + // It is an helper point and should be removed. + out_reliableVertices[vi] = false; + } + } + + ALICEVISION_LOG_DEBUG("DelaunayGraphCut::filterLargeHelperPoints done."); +} + +mesh::Mesh* DelaunayGraphCut::createMesh(int maxNbConnectedHelperPoints) { ALICEVISION_LOG_INFO("Extract mesh from Graph Cut."); @@ -3287,50 +3416,7 @@ mesh::Mesh* DelaunayGraphCut::createMesh(bool filterHelperPointsTriangles) } std::vector reliableVertices; - if(filterHelperPointsTriangles) - { - // Some vertices have not been created by depth maps but - // have been created for the tetrahedralization like - // camera centers, helper points, etc. - // These points have no visibility information (nrc == 0). - // We want to remove these fake points, but we don't want to create holes. - // So if the vertex is alone in the middle of valid points, we want to keep it. - // - // Algo: For each surface vertex without visibility, - // we check the neighbor surface vertices. If there is another one - // without visibility, we declare it unreliable. - reliableVertices.resize(_verticesCoords.size()); - for(VertexIndex vi = 0; vi < _verticesCoords.size(); ++vi) - { - if(!vertexIsOnSurface[vi]) - { - // this vertex is not on the surface, so no interest to spend time here - reliableVertices[vi] = false; - } - else if(_verticesAttr[vi].nrc > 0) - { - // this is a valid point without ambiguity - reliableVertices[vi] = true; - } - else - { - // this vertex has no visibility, check if it is connected to other weak vertices - reliableVertices[vi] = true; // reliable by default - GEO::vector neighbors; - _tetrahedralization->get_neighbors(vi, neighbors); - for(GEO::index_t nvi: neighbors) - { - if(vertexIsOnSurface[nvi] && (_verticesAttr[nvi].nrc == 0)) - { - // this vertex has no visibility and is connected to another - // surface vertex without visibility, so we declare it unreliable. - reliableVertices[vi] = false; - break; - } - } - } - } - } + filterLargeHelperPoints(reliableVertices, vertexIsOnSurface, maxNbConnectedHelperPoints); me->tris = StaticVector(); me->tris.reserve(nbSurfaceFacets); @@ -3366,7 +3452,7 @@ mesh::Mesh* DelaunayGraphCut::createMesh(bool filterHelperPointsTriangles) vertices[1] = getVertexIndex(f1, 1); vertices[2] = getVertexIndex(f1, 2); - if(filterHelperPointsTriangles) + if(!reliableVertices.empty()) { // We skip triangles if it contains one unreliable vertex. bool invalidTriangle = false; @@ -3707,21 +3793,44 @@ void DelaunayGraphCut::exportDebugMesh(const std::string& filename, const Point3 meshf->saveToObj(tempDirPath + "/" + filename + "Filtered.obj"); } -void DelaunayGraphCut::exportFullScoreMeshs() +void DelaunayGraphCut::exportFullScoreMeshs(const std::string& outputFolder, const std::string& name) const { - 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"); + const std::string nameExt = (name.empty() ? "" : "_" + name) + ".obj"; + { + std::unique_ptr meshEmptiness( + createTetrahedralMesh(false, 0.999f, [](const GC_cellInfo& c) { return c.emptinessScore; })); + meshEmptiness->saveToObj(outputFolder + "/mesh_emptiness" + nameExt); + } + { + std::unique_ptr meshFullness( + createTetrahedralMesh(false, 0.999f, [](const GC_cellInfo& c) { return c.fullnessScore; })); + meshFullness->saveToObj(outputFolder + "/mesh_fullness" + nameExt); + } + { + std::unique_ptr meshSWeight( + createTetrahedralMesh(false, 0.999f, [](const GC_cellInfo& c) { return c.cellSWeight; })); + meshSWeight->saveToObj(outputFolder + "/mesh_sWeight" + nameExt); + } + { + std::unique_ptr meshTWeight( + createTetrahedralMesh(false, 0.999f, [](const GC_cellInfo& c) { return c.cellTWeight; })); + meshTWeight->saveToObj(outputFolder + "/mesh_tWeight" + nameExt); + } + { + std::unique_ptr meshOn( + createTetrahedralMesh(false, 0.999f, [](const GC_cellInfo& c) { return c.on; })); + meshOn->saveToObj(outputFolder + "/mesh_on" + nameExt); + } + { + std::unique_ptr mesh(createTetrahedralMesh( + false, 0.99f, [](const fuseCut::GC_cellInfo& c) { return c.fullnessScore - c.emptinessScore; })); + mesh->saveToObj(outputFolder + "/mesh_fullness-emptiness" + nameExt); + } + { + std::unique_ptr mesh(createTetrahedralMesh( + false, 0.99f, [](const fuseCut::GC_cellInfo& c) { return c.cellSWeight - c.cellTWeight; })); + mesh->saveToObj(outputFolder + "/mesh_s-t" + nameExt); + } } void DelaunayGraphCut::exportBackPropagationMesh(const std::string& filename, std::vector& intersectedGeom, const Point3d& fromPt, const Point3d& toPt) diff --git a/src/aliceVision/fuseCut/DelaunayGraphCut.hpp b/src/aliceVision/fuseCut/DelaunayGraphCut.hpp index cf3f1449d0..11adb852a3 100644 --- a/src/aliceVision/fuseCut/DelaunayGraphCut.hpp +++ b/src/aliceVision/fuseCut/DelaunayGraphCut.hpp @@ -55,6 +55,9 @@ struct FuseParams float simGaussianSize = 10.0f; double minAngleThreshold = 0.1; bool refineFuse = true; + // Weight for helper points from mask. Do not create helper points if zero. + float maskHelperPointsWeight = 0.0; + int maskBorderSize = 1; }; @@ -434,10 +437,12 @@ class DelaunayGraphCut void addPointsFromCameraCenters(const StaticVector& cams, float minDist); void addPointsToPreventSingularities(const Point3d Voxel[8], float minDist); + void densifyWithHelperPoints(int nbFront, int nbBack, double scale); + /** * @brief Add volume points to prevent singularities */ - void addHelperPoints(int nGridHelperVolumePointsDim, const Point3d Voxel[8], float minDist); + void addGridHelperPoints(int helperPointsGridSize, const Point3d Voxel[8], float minDist); void addMaskHelperPoints(const Point3d voxel[8], const StaticVector& cams, const FuseParams& params); @@ -495,10 +500,10 @@ class DelaunayGraphCut float weightFcn(float nrc, bool labatutWeights, int ncams); - virtual void fillGraph(double nPixelSizeBehind, bool labatutWeights, - bool fillOut, float distFcnHeight = 0.0f); + void fillGraph(double nPixelSizeBehind, bool labatutWeights, bool fillOut, float distFcnHeight, + float fullWeight); void fillGraphPartPtRc(int& out_nstepsFront, int& out_nstepsBehind, GeometriesCount& outFrontCount, GeometriesCount& outBehindCount, int vertexIndex, int cam, float weight, - double nPixelSizeBehind, bool fillOut, float distFcnHeight); + float fullWeight, double nPixelSizeBehind, bool fillOut, float distFcnHeight); /** * @brief Estimate the cells property "on" based on the analysis of the visibility of neigbouring cells. @@ -517,7 +522,8 @@ class DelaunayGraphCut void createDensePointCloud(const Point3d hexah[8], const StaticVector& cams, const sfmData::SfMData* sfmData, const FuseParams* depthMapsFuseParams); - void createGraphCut(const Point3d hexah[8], const StaticVector& cams, const std::string& folderName, const std::string& tmpCamsPtsFolderName, bool removeSmallSegments); + void createGraphCut(const Point3d hexah[8], const StaticVector& cams, const std::string& folderName, + const std::string& tmpCamsPtsFolderName, bool removeSmallSegments, bool exportDebugTetrahedralization); /** * @brief Invert full/empty status of cells if they represent a too small group after labelling. @@ -538,12 +544,27 @@ class DelaunayGraphCut int removeDust(int minSegSize); void leaveLargestFullSegmentOnly(); - mesh::Mesh* createMesh(bool filterHelperPointsTriangles = true); + /** + * Some vertices have not been created by depth maps but + * have been created for the tetrahedralization like + * camera centers, helper points, etc. + * These points have no visibility information (nrc == 0). + * We want to remove these fake points, but we do not want to create holes. + * So if the vertex is alone in the middle of valid points, we want to keep it. + */ + void filterLargeHelperPoints(std::vector& out_reliableVertices, const std::vector& vertexIsOnSurface, int maxSegSize); + + /** + * @brief Create a mesh from the tetrahedral scores + * @param[in] maxNbConnectedHelperPoints: maximum number of connected helper points before we remove the group. 0 means that we remove all helper points. -1 means that we do not filter helper points at all. + */ + mesh::Mesh* createMesh(int maxNbConnectedHelperPoints); mesh::Mesh* createTetrahedralMesh(bool filter = true, const float& downscaleFactor = 0.95f, const std::function getScore = [](const GC_cellInfo& c) { return c.emptinessScore; }) const; void displayCellsStats() const; void exportDebugMesh(const std::string& filename, const Point3d& fromPt, const Point3d& toPt); - void exportFullScoreMeshs(); + void exportFullScoreMeshs(const std::string& outputFolder, const std::string& name) const; + 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); }; diff --git a/src/software/pipeline/main_meshing.cpp b/src/software/pipeline/main_meshing.cpp index 8b4d1ea8dc..23443fa360 100644 --- a/src/software/pipeline/main_meshing.cpp +++ b/src/software/pipeline/main_meshing.cpp @@ -265,7 +265,6 @@ int aliceVision_main(int argc, char* argv[]) bool meshingFromDepthMaps = true; bool estimateSpaceFromSfM = true; bool addLandmarksToTheDensePointCloud = false; - bool addMaskHelperPoints = false; bool saveRawDensePointCloud = false; bool colorizeOutput = false; bool voteFilteringForWeaklySupportedSurfaces = true; @@ -276,7 +275,15 @@ int aliceVision_main(int argc, char* argv[]) BoundingBox boundingBox; fuseCut::FuseParams fuseParams; + + int helperPointsGridSize = 10; + int densifyNbFront = 0; + int densifyNbBack = 0; + double densifyScale = 1.0; double nPixelSizeBehind = 4.0; + double fullWeight = 1.0; + bool exportDebugTetrahedralization = false; + int maxNbConnectedHelperPoints = 50; po::options_description allParams("AliceVision meshing"); @@ -345,10 +352,22 @@ int aliceVision_main(int argc, char* argv[]) "minAngleThreshold") ("refineFuse", po::value(&fuseParams.refineFuse)->default_value(fuseParams.refineFuse), "refineFuse") + ("helperPointsGridSize", po::value(&helperPointsGridSize)->default_value(helperPointsGridSize), + "Helper points grid size.") + ("densifyNbFront", po::value(&densifyNbFront)->default_value(densifyNbFront), + "Number of points in front of the vertices to densify the scene.") + ("densifyNbBack", po::value(&densifyNbBack)->default_value(densifyNbBack), + "Number of points behind the vertices to densify the scene.") + ("densifyScale", po::value(&densifyScale)->default_value(densifyScale), + "Scale between points used to densify the scene.") + ("maskHelperPointsWeight", po::value(&fuseParams.maskHelperPointsWeight)->default_value(fuseParams.maskHelperPointsWeight), + "Mask helper points weight. Zero to disable it.") + ("maskBorderSize", po::value(&fuseParams.maskBorderSize)->default_value(fuseParams.maskBorderSize), + "How many pixels on mask borders? 1 by default.") ("nPixelSizeBehind", po::value(&nPixelSizeBehind)->default_value(nPixelSizeBehind), "Number of pixel size units to vote behind the vertex with FULL status.") - ("addMaskHelperPoints", po::value(&addMaskHelperPoints)->default_value(addMaskHelperPoints), - "Add Helper points on the outline of the depth maps masks.") + ("fullWeight", po::value(&fullWeight)->default_value(fullWeight), + "Weighting of the FULL cells.") ("saveRawDensePointCloud", po::value(&saveRawDensePointCloud)->default_value(saveRawDensePointCloud), "Save dense point cloud before cut and filtering.") ("voteFilteringForWeaklySupportedSurfaces", po::value(&voteFilteringForWeaklySupportedSurfaces)->default_value(voteFilteringForWeaklySupportedSurfaces), @@ -358,9 +377,13 @@ int aliceVision_main(int argc, char* argv[]) ("minSolidAngleRatio", po::value(&minSolidAngleRatio)->default_value(minSolidAngleRatio), "Filter cells status on surface around vertices to improve smoothness using solid angle ratio between full/empty parts.") ("nbSolidAngleFilteringIterations", po::value(&nbSolidAngleFilteringIterations)->default_value(nbSolidAngleFilteringIterations), - "Number of iterations to filter the status cells based on solid angle ratio.") + "Number of iterations to filter the status cells based on solid angle ratio.") + ("maxNbConnectedHelperPoints", po::value(&maxNbConnectedHelperPoints)->default_value(maxNbConnectedHelperPoints), + "Maximum number of connected helper points before we remove them.") + ("exportDebugTetrahedralization", po::value(&exportDebugTetrahedralization)->default_value(exportDebugTetrahedralization), + "Export debug cells score as tetrahedral mesh. WARNING: could create huge meshes, only use on very small datasets.") ("seed", po::value(&seed)->default_value(seed), - "Seed used in random processes. (0 to use a random seed)."); + "Seed used in random processes. (0 to use a random seed)."); po::options_description logParams("Log parameters"); logParams.add_options() @@ -432,13 +455,18 @@ int aliceVision_main(int argc, char* argv[]) mvsUtils::MultiViewParams mp(sfmData, "", "", depthMapsFolder, meshingFromDepthMaps); mp.userParams.put("LargeScale.universePercentile", universePercentile); + mp.userParams.put("LargeScale.helperPointsGridSize", helperPointsGridSize); + mp.userParams.put("LargeScale.densifyNbFront", densifyNbFront); + mp.userParams.put("LargeScale.densifyNbBack", densifyNbBack); + mp.userParams.put("LargeScale.densifyScale", densifyScale); + mp.userParams.put("delaunaycut.seed", seed); mp.userParams.put("delaunaycut.nPixelSizeBehind", nPixelSizeBehind); + mp.userParams.put("delaunaycut.fullWeight", fullWeight); mp.userParams.put("delaunaycut.voteFilteringForWeaklySupportedSurfaces", voteFilteringForWeaklySupportedSurfaces); mp.userParams.put("hallucinationsFiltering.invertTetrahedronBasedOnNeighborsNbIterations", invertTetrahedronBasedOnNeighborsNbIterations); mp.userParams.put("hallucinationsFiltering.minSolidAngleRatio", minSolidAngleRatio); mp.userParams.put("hallucinationsFiltering.nbSolidAngleFilteringIterations", nbSolidAngleFilteringIterations); - mp.userParams.put("delaunaycut.addMaskHelperPoints", addMaskHelperPoints); int ocTreeDim = mp.userParams.get("LargeScale.gridLevel0", 1024); const auto baseDir = mp.userParams.get("LargeScale.baseDirName", "root01024"); @@ -518,10 +546,13 @@ int aliceVision_main(int argc, char* argv[]) sfmDataIO::Save(densePointCloud, (outDirectory/"densePointCloud_raw.abc").string(), sfmDataIO::ESfMData::ALL_DENSE); } - delaunayGC.createGraphCut(&hexah[0], cams, outDirectory.string()+"/", outDirectory.string()+"/SpaceCamsTracks/", false); + delaunayGC.createGraphCut(&hexah[0], cams, outDirectory.string() + "/", + outDirectory.string() + "/SpaceCamsTracks/", false, + exportDebugTetrahedralization); delaunayGC.graphCutPostProcessing(&hexah[0], outDirectory.string()+"/"); - mesh = delaunayGC.createMesh(); + + mesh = delaunayGC.createMesh(maxNbConnectedHelperPoints); delaunayGC.createPtsCams(ptsCams); mesh::meshPostProcessing(mesh, ptsCams, mp, outDirectory.string()+"/", nullptr, &hexah[0]);