diff --git a/apps/DensifyPointCloud/DensifyPointCloud.cpp b/apps/DensifyPointCloud/DensifyPointCloud.cpp index 757debefe..82b7d5daa 100644 --- a/apps/DensifyPointCloud/DensifyPointCloud.cpp +++ b/apps/DensifyPointCloud/DensifyPointCloud.cpp @@ -131,6 +131,7 @@ bool Application::Initialize(size_t argc, LPCTSTR* argv) unsigned nEstimationGeometricIters; unsigned nEstimateColors; unsigned nEstimateNormals; + unsigned nFuseFilter; unsigned nOptimize; int nIgnoreMaskLabel; bool bRemoveDmaps; @@ -146,7 +147,7 @@ bool Application::Initialize(size_t argc, LPCTSTR* argv) ("min-resolution", boost::program_options::value(&nMinResolution)->default_value(640), "do not scale images lower than this resolution") ("sub-resolution-levels", boost::program_options::value(&nSubResolutionLevels)->default_value(2), "number of patch-match sub-resolution iterations (0 - disabled)") ("number-views", boost::program_options::value(&nNumViews)->default_value(nNumViewsDefault), "number of views used for depth-map estimation (0 - all neighbor views available)") - ("number-views-fuse", boost::program_options::value(&nMinViewsFuse)->default_value(3), "minimum number of images that agrees with an estimate during fusion in order to consider it inlier (<2 - only merge depth-maps)") + ("number-views-fuse", boost::program_options::value(&nMinViewsFuse)->default_value(2), "minimum number of images that agrees with an estimate during fusion in order to consider it inlier (<2 - only merge depth-maps)") ("ignore-mask-label", boost::program_options::value(&nIgnoreMaskLabel)->default_value(-1), "label value to ignore in the image mask, stored in the MVS scene or next to each image with '.mask.png' extension (<0 - disabled)") ("mask-path", boost::program_options::value(&OPT::strMaskPath), "path to folder containing mask images with '.mask.png' extension") ("iters", boost::program_options::value(&nEstimationIters)->default_value(numIters), "number of patch-match iterations") @@ -157,6 +158,7 @@ bool Application::Initialize(size_t argc, LPCTSTR* argv) ("sub-scene-area", boost::program_options::value(&OPT::fMaxSubsceneArea)->default_value(0.f), "split the scene in sub-scenes such that each sub-scene surface does not exceed the given maximum sampling area (0 - disabled)") ("sample-mesh", boost::program_options::value(&OPT::fSampleMesh)->default_value(0.f), "uniformly samples points on a mesh (0 - disabled, <0 - number of points, >0 - sample density per square unit)") ("fusion-mode", boost::program_options::value(&OPT::nFusionMode)->default_value(0), "depth-maps fusion mode (-2 - fuse disparity-maps, -1 - export disparity-maps only, 0 - depth-maps & fusion, 1 - export depth-maps only)") + ("fusion-filter", boost::program_options::value(&nFuseFilter)->default_value(2), "filter used to fuse the depth-maps (0 - merge, 1 - fuse, 2 - dense-fuse)") ("postprocess-dmaps", boost::program_options::value(&nOptimize)->default_value(7), "flags used to filter the depth-maps after estimation (0 - disabled, 1 - remove-speckles, 2 - fill-gaps, 4 - adjust-filter)") ("filter-point-cloud", boost::program_options::value(&OPT::thFilterPointCloud)->default_value(0), "filter dense point-cloud based on visibility (0 - disabled)") ("export-number-views", boost::program_options::value(&OPT::nExportNumViews)->default_value(0), "export points with >= number of views (0 - disabled, <0 - save MVS project too)") @@ -248,6 +250,7 @@ bool Application::Initialize(size_t argc, LPCTSTR* argv) OPTDENSE::nEstimationGeometricIters = nEstimationGeometricIters; OPTDENSE::nEstimateColors = nEstimateColors; OPTDENSE::nEstimateNormals = nEstimateNormals; + OPTDENSE::nFuseFilter = nFuseFilter; OPTDENSE::nOptimize = nOptimize; OPTDENSE::nIgnoreMaskLabel = nIgnoreMaskLabel; OPTDENSE::bRemoveDmaps = bRemoveDmaps; diff --git a/apps/Tests/Tests.cpp b/apps/Tests/Tests.cpp index f0f8365b8..b8af6437d 100644 --- a/apps/Tests/Tests.cpp +++ b/apps/Tests/Tests.cpp @@ -82,21 +82,21 @@ bool PipelineTest(bool verbose=false) } OPTDENSE::init(); OPTDENSE::bRemoveDmaps = true; - if (!scene.DenseReconstruction() || scene.pointcloud.GetSize() < 200000u) { + if (!scene.DenseReconstruction() || scene.pointcloud.GetSize() < 50000u) { VERBOSE("ERROR: TestDataset failed estimating dense point cloud!"); return false; } if (verbose) scene.pointcloud.Save(MAKE_PATH("scene_dense.ply")); - if (!scene.ReconstructMesh() || scene.mesh.faces.size() < 75000u) { + if (!scene.ReconstructMesh() || scene.mesh.faces.size() < 40000u) { VERBOSE("ERROR: TestDataset failed reconstructing the mesh!"); return false; } if (verbose) scene.mesh.Save(MAKE_PATH("scene_dense_mesh.ply")); - constexpr float decimate = 0.5f; + constexpr float decimate = 0.7f; scene.mesh.Clean(decimate); - if (!ISINSIDE(scene.mesh.faces.size(), 35000u, 45000u)) { + if (!ISINSIDE(scene.mesh.faces.size(), 25000u, 45000u)) { VERBOSE("ERROR: TestDataset failed cleaning the mesh!"); return false; } diff --git a/libs/Common/Types.h b/libs/Common/Types.h index a0923b2b7..a34499a62 100644 --- a/libs/Common/Types.h +++ b/libs/Common/Types.h @@ -1298,9 +1298,15 @@ class TPoint2 : public cv::Point_ inline const TYPE* ptr() const { return &x; } inline TYPE* ptr() { return &x; } + // iterator base access to enable range-based for loops + inline const TYPE* begin() const { return &x; } + inline const TYPE* end() const { return &x+3; } + // 1D element access - inline const TYPE& operator [](size_t i) const { ASSERT(i>=0 && i<2); return ptr()[i]; } - inline TYPE& operator [](size_t i) { ASSERT(i>=0 && i<2); return ptr()[i]; } + inline const TYPE& operator ()(int i) const { ASSERT(i>=0 && i<2); return ptr()[i]; } + inline TYPE& operator ()(int i) { ASSERT(i>=0 && i<2); return ptr()[i]; } + inline const TYPE& operator [](int i) const { ASSERT(i>=0 && i<2); return ptr()[i]; } + inline TYPE& operator [](int i) { ASSERT(i>=0 && i<2); return ptr()[i]; } // Access point as Size equivalent inline operator const Size& () const { return *((const Size*)this); } @@ -1391,9 +1397,15 @@ class TPoint3 : public cv::Point3_ inline const TYPE* ptr() const { return &x; } inline TYPE* ptr() { return &x; } + // iterator base access to enable range-based for loops + inline const TYPE* begin() const { return &x; } + inline const TYPE* end() const { return &x+3; } + // 1D element access - inline const TYPE& operator [](BYTE i) const { ASSERT(i<3); return ptr()[i]; } - inline TYPE& operator [](BYTE i) { ASSERT(i<3); return ptr()[i]; } + inline const TYPE& operator ()(int i) const { ASSERT(i>=0 && i<3); return ptr()[i]; } + inline TYPE& operator ()(int i) { ASSERT(i>=0 && i<3); return ptr()[i]; } + inline const TYPE& operator [](int i) const { ASSERT(i>=0 && i<3); return ptr()[i]; } + inline TYPE& operator [](int i) { ASSERT(i>=0 && i<3); return ptr()[i]; } // Access point as vector equivalent inline operator const Vec& () const { return *reinterpret_cast(this); } diff --git a/libs/MVS/DepthMap.cpp b/libs/MVS/DepthMap.cpp index db39d19f5..5e4a7eca1 100644 --- a/libs/MVS/DepthMap.cpp +++ b/libs/MVS/DepthMap.cpp @@ -72,11 +72,15 @@ MDEFVAR_OPTDENSE_uint32(nMinResolution, "Min Resolution", "Do not scale images l DEFVAR_OPTDENSE_uint32(nSubResolutionLevels, "SubResolution levels", "Number of lower resolution levels to estimate the depth and normals", "2") DEFVAR_OPTDENSE_uint32(nMinViews, "Min Views", "minimum number of agreeing views to validate a depth", "2") MDEFVAR_OPTDENSE_uint32(nMaxViews, "Max Views", "maximum number of neighbor images used to compute the depth-map for the reference image", "12") -DEFVAR_OPTDENSE_uint32(nMinViewsFuse, "Min Views Fuse", "minimum number of images that agrees with an estimate during fusion in order to consider it inlier (<2 - only merge depth-maps)", "2") +DEFVAR_OPTDENSE_uint32(nMinViewsFuse, "Min Views Fuse", "minimum number of images that agrees with an estimate during fusion in order to consider it inlier", "2") +MDEFVAR_OPTDENSE_uint32(nMaxViewsFuse, "Max Views Fuse", "maximum number of neighbor depth-maps used during fusion", "32") DEFVAR_OPTDENSE_uint32(nMinViewsFilter, "Min Views Filter", "minimum number of images that agrees with an estimate in order to consider it inlier", "2") MDEFVAR_OPTDENSE_uint32(nMinViewsFilterAdjust, "Min Views Filter Adjust", "minimum number of images that agrees with an estimate in order to consider it inlier (0 - disabled)", "1") MDEFVAR_OPTDENSE_uint32(nMinViewsTrustPoint, "Min Views Trust Point", "min-number of views so that the point is considered for approximating the depth-maps (<2 - random initialization)", "2") MDEFVAR_OPTDENSE_uint32(nNumViews, "Num Views", "Number of views used for depth-map estimation (0 - all views available)", "0", "4", "8") +MDEFVAR_OPTDENSE_uint32(nMinPixelsFuse, "Min Pixels Fuse", "minimum number of depth-estimates that agree during fusion in order to consider it (multiple pixels can be from the same depth-map)", "5") +MDEFVAR_OPTDENSE_uint32(nMaxPointsFuse, "Max Points Fuse", "maximum number of pixels to fuse into a single point", "1000") +MDEFVAR_OPTDENSE_uint32(nMaxFuseDepth, "Max Fuse Depth", "maximum depth in fusion graph traversal", "100") MDEFVAR_OPTDENSE_uint32(nPointInsideROI, "Point Inside ROI", "consider a point shared only if inside ROI when estimating the neighbor views (0 - ignore ROI, 1 - weight more ROI points, 2 - consider only ROI points)", "1") MDEFVAR_OPTDENSE_bool(bFilterAdjust, "Filter Adjust", "adjust depth estimates during filtering", "1") MDEFVAR_OPTDENSE_bool(bAddCorners, "Add Corners", "add support points at image corners with nearest neighbor disparities", "0") @@ -89,6 +93,7 @@ MDEFVAR_OPTDENSE_float(fMinAngle, "Min Angle", "Min angle for accepting the dept MDEFVAR_OPTDENSE_float(fOptimAngle, "Optim Angle", "Optimal angle for computing the depth triangulation", "12.0") MDEFVAR_OPTDENSE_float(fMaxAngle, "Max Angle", "Max angle for accepting the depth triangulation", "65.0") MDEFVAR_OPTDENSE_float(fDescriptorMinMagnitudeThreshold, "Descriptor Min Magnitude Threshold", "minimum patch texture variance accepted when matching two patches (0 - disabled)", "0.02") // 0.02: pixels with patch texture variance below 0.0004 (0.02^2) will be removed from depthmap; 0.12: patch texture variance below 0.02 (0.12^2) is considered texture-less +MDEFVAR_OPTDENSE_float(fDepthReprojectionErrorThreshold, "Depth Reprojection Error Threshold", "maximum relative difference between measured and depth projected pixel", "1.2") MDEFVAR_OPTDENSE_float(fDepthDiffThreshold, "Depth Diff Threshold", "maximum variance allowed for the depths during refinement", "0.01") MDEFVAR_OPTDENSE_float(fNormalDiffThreshold, "Normal Diff Threshold", "maximum variance allowed for the normal during fusion (degrees)", "25") MDEFVAR_OPTDENSE_float(fPairwiseMul, "Pairwise Mul", "pairwise cost scale to match the unary cost", "0.3") @@ -98,8 +103,9 @@ MDEFVAR_OPTDENSE_uint32(nSpeckleSize, "Speckle Size", "maximal size of a speckle MDEFVAR_OPTDENSE_uint32(nIpolGapSize, "Interpolate Gap Size", "interpolate small gaps (left<->right, top<->bottom)", "7") MDEFVAR_OPTDENSE_int32(nIgnoreMaskLabel, "Ignore Mask Label", "label id used during ignore mask filter (<0 - disabled)", "-1") DEFVAR_OPTDENSE_uint32(nOptimize, "Optimize", "should we filter the extracted depth-maps?", "7") // see DepthFlags +DEFVAR_OPTDENSE_uint32(nFuseFilter, "Fuse Filter", "how to fuse the depth-maps into one dense point-cloud?", "2", "0", "1") // see FuseMode MDEFVAR_OPTDENSE_uint32(nEstimateColors, "Estimate Colors", "should we estimate the colors for the dense point-cloud?", "2", "0", "1") -MDEFVAR_OPTDENSE_uint32(nEstimateNormals, "Estimate Normals", "should we estimate the normals for the dense point-cloud?", "0", "1", "2") +MDEFVAR_OPTDENSE_uint32(nEstimateNormals, "Estimate Normals", "should we estimate the normals for the dense point-cloud?", "2", "0", "1") MDEFVAR_OPTDENSE_float(fNCCThresholdKeep, "NCC Threshold Keep", "Maximum 1-NCC score accepted for a match", "0.9", "0.5") DEFVAR_OPTDENSE_uint32(nEstimationIters, "Estimation Iters", "Number of patch-match iterations", "3") DEFVAR_OPTDENSE_uint32(nEstimationGeometricIters, "Estimation Geometric Iters", "Number of geometric consistent patch-match iterations (0 - disabled)", "2") diff --git a/libs/MVS/DepthMap.h b/libs/MVS/DepthMap.h index a5fabf2c8..52ed84cd5 100644 --- a/libs/MVS/DepthMap.h +++ b/libs/MVS/DepthMap.h @@ -90,6 +90,11 @@ enum DepthFlags { ADJUST_FILTER = (1 << 2), OPTIMIZE = (REMOVE_SPECKLES|FILL_GAPS) }; +enum FuseMode { + FUSE_NOFILTER = 0, + FUSE_FILTER, + FUSE_DENSEFILTER, +}; extern unsigned nResolutionLevel; extern unsigned nMaxResolution; extern unsigned nMinResolution; @@ -97,10 +102,14 @@ extern unsigned nSubResolutionLevels; extern unsigned nMinViews; extern unsigned nMaxViews; extern unsigned nMinViewsFuse; +extern unsigned nMaxViewsFuse; extern unsigned nMinViewsFilter; extern unsigned nMinViewsFilterAdjust; extern unsigned nMinViewsTrustPoint; extern unsigned nNumViews; +extern unsigned nMinPixelsFuse; +extern unsigned nMaxPointsFuse; +extern unsigned nMaxFuseDepth; extern unsigned nPointInsideROI; extern bool bFilterAdjust; extern bool bAddCorners; @@ -113,6 +122,7 @@ extern float fMinAngle; extern float fOptimAngle; extern float fMaxAngle; extern float fDescriptorMinMagnitudeThreshold; +extern float fDepthReprojectionErrorThreshold; extern float fDepthDiffThreshold; extern float fNormalDiffThreshold; extern float fPairwiseMul; @@ -122,6 +132,7 @@ extern unsigned nSpeckleSize; extern unsigned nIpolGapSize; extern int nIgnoreMaskLabel; extern unsigned nOptimize; +extern unsigned nFuseFilter; extern unsigned nEstimateColors; extern unsigned nEstimateNormals; extern float fNCCThresholdKeep; diff --git a/libs/MVS/SceneDensify.cpp b/libs/MVS/SceneDensify.cpp index 157259331..9a1f1bf1b 100644 --- a/libs/MVS/SceneDensify.cpp +++ b/libs/MVS/SceneDensify.cpp @@ -1295,7 +1295,7 @@ void DepthMapsData::MergeDepthMaps(PointCloud& pointcloud, bool bEstimateColor, // compute available memory to be used for depth-data caching // - numDMapsReserveFusion: maximum number of depth-maps for which to reserve memory for fusion -size_t GetAvailableMemory(const DepthDataArr& arrDepthData, const std::unordered_set& fusedDMaps, IIndex numDMapsReserveFusion, size_t currentCacheMemory = 0) +size_t GetAvailableMemory(const DepthDataArr& arrDepthData, const BoolArr& fusedDMaps, IIndex numDMapsReserveFusion, size_t currentCacheMemory = 0) { size_t resolution(0); IIndex numDMaps(0); @@ -1303,7 +1303,7 @@ size_t GetAvailableMemory(const DepthDataArr& arrDepthData, const std::unordered const DepthData& depthData = arrDepthData[idxImage]; if (!depthData.IsValid()) continue; - if (fusedDMaps.count(idxImage)) + if (fusedDMaps[idxImage]) continue; resolution += depthData.size.area(); if (++numDMaps >= numDMapsReserveFusion) @@ -1325,7 +1325,7 @@ size_t GetAvailableMemory(const DepthDataArr& arrDepthData, const std::unordered } // GetAvailableMemory // finds the best depth-map to fuse next that maximizes the number of neighbors already in cache -std::tuple FetchBestNextDMapIndex(const DepthDataArr& arrDepthData, DMapCache& cacheDMaps, std::unordered_set& fusedDMaps) { +std::tuple FetchBestNextDMapIndex(const DepthDataArr& arrDepthData, const DMapCache& cacheDMaps, const BoolArr& fusedDMaps) { const IIndexArr cachedImages = cacheDMaps.GetCachedImageIndices(true); IIndex bestImageIdx = NO_ID; unsigned bestImageScore = 0, bestImageSize = std::numeric_limits::max(); @@ -1333,7 +1333,7 @@ std::tuple FetchBestNextDMapIndex(const DepthDataA const DepthData& depthData = arrDepthData[idxImage]; if (!depthData.IsValid()) continue; - if (fusedDMaps.count(idxImage)) + if (fusedDMaps[idxImage]) continue; ASSERT(!depthData.neighbors.empty()); IIndexArr cachedNeighbors; @@ -1379,15 +1379,15 @@ void DepthMapsData::FuseDepthMaps(PointCloud& pointcloud, bool bEstimateColor, b typedef SEACAVE::cList ProjsArr; // fuse all depth-maps, processing the best connected images first - const unsigned nMinViewsFuse(MINF(OPTDENSE::nMinViewsFuse, scene.images.size())); + const unsigned nMinViewsFuse(MINF(OPTDENSE::nMinViewsFuse, arrDepthData.size())); const float normalError(COS(FD2R(OPTDENSE::fNormalDiffThreshold))); const IIndex numDMapsReserveFusion(10); CLISTDEF0(Depth*) invalidDepths(0, 32); size_t nDepths(0); typedef TImage DepthIndex; typedef cList DepthIndexArr; - DepthIndexArr arrDepthIdx(scene.images.size()); - const size_t nPointsEstimate(scene.images.size() * 9000); //TODO: better estimate number of points + DepthIndexArr arrDepthIdx(arrDepthData.size()); + const size_t nPointsEstimate(arrDepthData.size() * 9000); //TODO: better estimate number of points ProjsArr projs(0, nPointsEstimate); pointcloud.points.reserve(nPointsEstimate); pointcloud.pointViews.reserve(nPointsEstimate); @@ -1399,12 +1399,14 @@ void DepthMapsData::FuseDepthMaps(PointCloud& pointcloud, bool bEstimateColor, b pointcloud.normals.reserve(nPointsEstimate); depthDataLoadFlags |= HeaderDepthDataRaw::HAS_NORMAL; } - Util::Progress progress(_T("Fused depth-maps"), scene.images.size()); + Util::Progress progress(_T("Fused depth-maps"), arrDepthData.size()); GET_LOGCONSOLE().Pause(); - std::unordered_set fusedDMaps; + BoolArr fusedDMaps(arrDepthData.size()); + fusedDMaps.Memset(0); DMapCache cacheDMaps(arrDepthData, depthDataLoadFlags, GetAvailableMemory(arrDepthData, fusedDMaps, numDMapsReserveFusion)); unsigned totalNumImageNeighborsInCache = 0, totalNumImagesInCache = 0; - while (fusedDMaps.size() < arrDepthData.size()) { + IIndex numDMapsFused = 0; + for (; numDMapsFused < arrDepthData.size(); ++numDMapsFused) { TD_TIMER_STARTD(); // find the best depth-map to fuse next as the one with the most neighbors already in cache const auto [idxImage, numImageNeighborsInCache, numImagesInCache] = FetchBestNextDMapIndex(arrDepthData, cacheDMaps, fusedDMaps); @@ -1412,7 +1414,6 @@ void DepthMapsData::FuseDepthMaps(PointCloud& pointcloud, bool bEstimateColor, b break; // no more depth-maps to fuse (only invalid depth-maps left) totalNumImageNeighborsInCache += numImageNeighborsInCache; totalNumImagesInCache += numImagesInCache; - fusedDMaps.emplace(idxImage); // fuse depth-map cacheDMaps.UseImage(idxImage); cacheDMaps.SkipMemoryCheckIdxImage(idxImage); @@ -1422,6 +1423,7 @@ void DepthMapsData::FuseDepthMaps(PointCloud& pointcloud, bool bEstimateColor, b if (bEstimateNormal && depthData.normalMap.empty()) EstimateNormalMaps(); ASSERT(!depthData.images.empty() && !depthData.neighbors.empty()); + IIndex numNeighbors(0); #ifdef DENSE_USE_OPENMP #pragma omp parallel for for (int64_t i=0; i<(int64_t)depthData.neighbors.size(); ++i) { @@ -1435,6 +1437,12 @@ void DepthMapsData::FuseDepthMaps(PointCloud& pointcloud, bool bEstimateColor, b cacheDMaps.UseImage(neighbor.ID); if (depthDataB.IsEmpty()) continue; + if (++numNeighbors >= OPTDENSE::nMaxViewsFuse) + #ifdef DENSE_USE_OPENMP + continue; + #else + break; + #endif DepthIndex& depthIdxs = arrDepthIdx[neighbor.ID]; if (!depthIdxs.empty()) continue; @@ -1442,12 +1450,12 @@ void DepthMapsData::FuseDepthMaps(PointCloud& pointcloud, bool bEstimateColor, b depthIdxs.memset((uint8_t)NO_ID); } ASSERT(!depthData.IsEmpty()); - ASSERT(depthData.depthMap.size() == depthData.size); const Image& imageData = *depthData.images.front().pImageData; ASSERT(&imageData-scene.images.data() == idxImage); + ASSERT(depthData.depthMap.size() == depthData.size && imageData.GetSize() == depthData.size); DepthIndex& depthIdxs = arrDepthIdx[idxImage]; if (depthIdxs.empty()) { - depthIdxs.create(imageData.GetSize()); + depthIdxs.create(depthData.size); depthIdxs.memset((uint8_t)NO_ID); } const size_t nNumPointsPrev(pointcloud.points.size()); @@ -1551,13 +1559,14 @@ void DepthMapsData::FuseDepthMaps(PointCloud& pointcloud, bool bEstimateColor, b } } } + fusedDMaps[idxImage] = true; ASSERT(pointcloud.points.size() == pointcloud.pointViews.size() && pointcloud.points.size() == pointcloud.pointWeights.size() && pointcloud.points.size() == projs.size()); DEBUG_ULTIMATE("Depths map for reference image %3u fused using %u depths maps: %u new points, %u/%u cached images (%s)", idxImage, depthData.images.size()-1, pointcloud.points.size()-nNumPointsPrev, numImageNeighborsInCache, numImagesInCache, TD_TIMER_GET_FMT().c_str()); - progress.display(fusedDMaps.size()); + progress.display(numDMapsFused); // ensure enough memory is available for the next depth-maps chunk cacheDMaps.SkipMemoryCheckIdxImage(); - if (fusedDMaps.size() % numDMapsReserveFusion == 0) + if (numDMapsFused % numDMapsReserveFusion == 0) cacheDMaps.SetMaxMemory(GetAvailableMemory(arrDepthData, fusedDMaps, numDMapsReserveFusion, cacheDMaps.GetUsedMemory())); } GET_LOGCONSOLE().Play(); @@ -1566,9 +1575,9 @@ void DepthMapsData::FuseDepthMaps(PointCloud& pointcloud, bool bEstimateColor, b cacheDMaps.ClearCache(); DEBUG_EXTRA("Depth-maps fused and filtered: %u depth-maps, %u depths, %u points (%d%%%%), %.2f hits in %.2f cached (%s)", - fusedDMaps.size(), nDepths, pointcloud.points.size(), ROUND2INT((100.f*pointcloud.points.size())/nDepths), - static_cast(totalNumImageNeighborsInCache) / fusedDMaps.size(), - static_cast(totalNumImagesInCache) / fusedDMaps.size(), TD_TIMER_GET_FMT().c_str()); + numDMapsFused, nDepths, pointcloud.points.size(), ROUND2INT((100.f*pointcloud.points.size())/nDepths), + static_cast(totalNumImageNeighborsInCache) / numDMapsFused, + static_cast(totalNumImagesInCache) / numDMapsFused, TD_TIMER_GET_FMT().c_str()); if (bEstimateNormal && !pointcloud.points.empty() && pointcloud.normals.empty()) { // estimate normal also if requested (quite expensive if normal-maps not available) @@ -1597,6 +1606,252 @@ void DepthMapsData::FuseDepthMaps(PointCloud& pointcloud, bool bEstimateColor, b DEBUG_EXTRA("Normals estimated for the dense point-cloud: %u normals (%s)", pointcloud.GetSize(), TD_TIMER_GET_FMT().c_str()); } } // FuseDepthMaps + + +// fuse all valid depth-maps in the same 3D point cloud; +// join points very likely to represent the same 3D point and +// filter out points blocking the view +void DepthMapsData::DenseFuseDepthMaps(PointCloud& pointcloud, bool bEstimateColor, bool _bEstimateNormal) +{ + TD_TIMER_STARTD(); + + typedef SEACAVE::BitMatrix UseMask; + typedef CLISTDEFIDX(UseMask,IIndex) UseMaskArr; + + // fuse all depth-maps, processing the best connected images first + const unsigned nMinViewsFuse(MINF(OPTDENSE::nMinViewsFuse, arrDepthData.size())); + const float normalError(COS(FD2R(OPTDENSE::fNormalDiffThreshold))); + const float minConfidence(1.f - OPTDENSE::fNCCThresholdKeep); + const float maxReprojErrorSq(SQUARE(OPTDENSE::fDepthReprojectionErrorThreshold)); + const IIndex numDMapsReserveFusion(10); + const bool bEstimateNormal(true); // always estimate normals as they are needed for the fusion + size_t nDepths(0); + UseMaskArr arrUseMask(arrDepthData.size()); + const size_t nPointsEstimate(arrDepthData.size() * 9000); //TODO: better estimate number of points + pointcloud.points.reserve(nPointsEstimate); + pointcloud.pointViews.reserve(nPointsEstimate); + pointcloud.pointWeights.reserve(nPointsEstimate); + unsigned depthDataLoadFlags(HeaderDepthDataRaw::HAS_DEPTH | HeaderDepthDataRaw::HAS_CONF); + if (bEstimateColor) + pointcloud.colors.reserve(nPointsEstimate); + if (bEstimateNormal) { + pointcloud.normals.reserve(nPointsEstimate); + depthDataLoadFlags |= HeaderDepthDataRaw::HAS_NORMAL; + } + Util::Progress progress(_T("Dense fused depth-maps"), arrDepthData.size()); + GET_LOGCONSOLE().Pause(); + BoolArr fusedDMaps(arrDepthData.size()); + fusedDMaps.Memset(0); + DMapCache cacheDMaps(arrDepthData, depthDataLoadFlags, GetAvailableMemory(arrDepthData, fusedDMaps, numDMapsReserveFusion)); + unsigned totalNumImageNeighborsInCache = 0, totalNumImagesInCache = 0; + BoolArr neighbors(arrDepthData.size()); + PointCloud::Point refPoint; + PointCloud::Normal refNormal; + CLISTDEF0IDX(float, unsigned) fusedPoints[3]; + PointCloud::ViewArr fusedViews; + FloatArr fusedWeights; + Point3d fusedNormal; + Pixel32F fusedColor; + const auto FusePoint = [&](IIndex ID, const ImageRef& x, unsigned fuseDepth) -> void { + const auto lambda = [&](IIndex ID, const ImageRef& x, unsigned fuseDepth, const auto& FusePointImpl) -> void { + const DepthData& depthData = arrDepthData[ID]; + if (!Image8U::isInside(x, depthData.size)) + return; + // ignore pixel if not estimated + ASSERT(depthData.depthMap.size() == depthData.size); + const Depth depth = depthData.depthMap(x); + if (depth <= Depth(0)) + return; + ASSERT(ISINSIDE(depth, depthData.dMin * 0.95f, depthData.dMax * 1.05f)); + // ignore pixel if already fused + UseMask& useMask = arrUseMask[ID]; + if (useMask(x)) + return; + // ignore pixel if not confident + const float conf(depthData.confMap.empty() ? 1.f : depthData.confMap(x)); + if (conf < minConfidence) + return; + const DepthData::ViewData& image = depthData.GetView(); + // if the fusion depth is greater than zero, the initial reference pixel + // has already been added and we need to check for consistency + PointCloud::Normal normal; + if (fuseDepth > 0) { + // project reference point into current view + const Point3f pt(image.camera.ProjectPointP3(refPoint)); + // check if depth agrees with current depth + ASSERT(pt.z > Depth(0) || !IsDepthSimilar(depth, pt.z, OPTDENSE::fDepthDiffThreshold)); + if (!IsDepthSimilar(depth, pt.z, OPTDENSE::fDepthDiffThreshold)) + return; + // check reprojection error of the reference point in the current view + const Point2f diff(pt.x / pt.z - float(x.x), pt.y / pt.z - float(x.y)); + if (normSq(diff) > maxReprojErrorSq) + return; + // check if normals agree + normal = image.camera.R.t() * Cast(depthData.normalMap(x)); + ASSERT(ISEQUAL(norm(normal), 1.f)); + if (refNormal.dot(normal) < normalError) + return; + } else { + normal = image.camera.R.t() * Cast(depthData.normalMap(x)); + ASSERT(ISEQUAL(norm(normal), 1.f)); + } + // set the current pixel as visited + useMask.set(x); + // compute 3D location of the current depth + const PointCloud::Point X(image.camera.TransformPointI2W(Point3(REAL(x.x), REAL(x.y), REAL(depth)))); + // accumulate statistics for fused point + { + fusedPoints[0].push_back(X(0)); + fusedPoints[1].push_back(X(1)); + fusedPoints[2].push_back(X(2)); + const float weight(Conf2Weight(conf, depth)); + const auto it(fusedViews.InsertSortUnique(ID)); + if (it.second) + fusedWeights[it.first] += weight; + else + fusedWeights.InsertAt(it.first, weight); + if (bEstimateNormal) + fusedNormal += Cast(normal); + if (bEstimateColor) + fusedColor += Cast(image.pImageData->image(x)); + } + // remember the first pixel as the reference. + if (fuseDepth == 0) { + refPoint = X; + refNormal = normal; + } + // do not traverse the graph infinitely in one branch and + // limit the maximum number of pixels fused in one point + // to avoid stack overflow + if (++fuseDepth >= OPTDENSE::nMaxFuseDepth || fusedPoints[0].size() >= OPTDENSE::nMaxPointsFuse) + return; + // traverse the neighbors graph by projecting the point into other views + for (const ViewScore& neighbor : image.pImageData->neighbors) { + const IIndex nextID(neighbor.ID); + ASSERT(nextID != ID); + if (!neighbors[nextID]) + continue; + const DepthData& nextDepthData = arrDepthData[nextID]; + const ImageRef nextx(ROUND2INT(nextDepthData.GetCamera().ProjectPointP(X))); + FusePointImpl(nextID, nextx, fuseDepth, FusePointImpl); + } + }; + lambda(ID, x, fuseDepth, lambda); + }; + // loop over each depth-map + IIndex numDMapsFused = 0; + for (; numDMapsFused < arrDepthData.size(); ++numDMapsFused) { + TD_TIMER_STARTD(); + // find the best depth-map to fuse next as the one with the most neighbors already in cache + const auto [idxImage, numImageNeighborsInCache, numImagesInCache] = FetchBestNextDMapIndex(arrDepthData, cacheDMaps, fusedDMaps); + if (idxImage == NO_ID) + break; // no more depth-maps to fuse (only invalid depth-maps left) + totalNumImageNeighborsInCache += numImageNeighborsInCache; + totalNumImagesInCache += numImagesInCache; + // fuse depth-map + cacheDMaps.UseImage(idxImage); + cacheDMaps.SkipMemoryCheckIdxImage(idxImage); + const DepthData& depthData(arrDepthData[idxImage]); + ASSERT(depthData.GetView().GetLocalID(scene.images) == idxImage); + ASSERT(!depthData.IsEmpty()); + if (bEstimateNormal && depthData.normalMap.empty()) + EstimateNormalMaps(); + // make sure all neighbors are cached + neighbors.Memset(0); + neighbors[idxImage] = true; + IIndex numNeighbors(0); + ASSERT(!depthData.images.empty() && !depthData.neighbors.empty()); + #ifdef DENSE_USE_OPENMP + #pragma omp parallel for + for (int64_t i=0; i<(int64_t)depthData.neighbors.size(); ++i) { + const ViewScore& neighbor = depthData.neighbors[(IIndex)i]; + #else + for (const ViewScore& neighbor: depthData.neighbors) { + #endif + const DepthData& depthDataB(arrDepthData[neighbor.ID]); + if (!depthDataB.IsValid()) + continue; + cacheDMaps.UseImage(neighbor.ID); + if (depthDataB.IsEmpty()) + continue; + neighbors[neighbor.ID] = true; + if (++numNeighbors >= OPTDENSE::nMaxViewsFuse) + #ifdef DENSE_USE_OPENMP + continue; + #else + break; + #endif + UseMask& useMask = arrUseMask[neighbor.ID]; + if (!useMask.empty()) + continue; + useMask.create(depthDataB.depthMap.size()); + useMask.memset(0); + } + ASSERT(!depthData.IsEmpty()); + const Image& imageData = *depthData.images.front().pImageData; + ASSERT(&imageData-scene.images.data() == idxImage); + ASSERT(depthData.depthMap.size() == depthData.size && imageData.GetSize() == depthData.size); + UseMask& useMask = arrUseMask[idxImage]; + if (useMask.empty()) { + useMask.create(depthData.size); + useMask.memset(0); + } + // try to fuse each depth estimate + const size_t nNumPointsPrev(pointcloud.points.size()); + for (int i=0; i= OPTDENSE::nMinPixelsFuse && fusedViews.size() >= nMinViewsFuse) { + // create the corresponding 3D point + pointcloud.points.emplace_back( + fusedPoints[0].GetMedian(), + fusedPoints[1].GetMedian(), + fusedPoints[2].GetMedian() + ); + ASSERT(fusedViews.size() == fusedWeights.size()); + PointCloud::WeightArr& weights = pointcloud.pointWeights.AddEmpty(); + for (float weight: fusedWeights) + weights.push_back(weight); + pointcloud.pointViews.emplace_back(fusedViews); + if (bEstimateNormal) + pointcloud.normals.emplace_back(normalized(fusedNormal)); + if (bEstimateColor) + pointcloud.colors.emplace_back((fusedColor/(float)fusedPoints[0].size()).cast()); + } + if (!fusedViews.empty()) { + nDepths += fusedViews.size(); + fusedPoints[0].clear(); + fusedPoints[1].clear(); + fusedPoints[2].clear(); + fusedViews.clear(); + fusedWeights.clear(); + fusedNormal = Point3d::ZERO; + fusedColor = Pixel32F::BLACK; + } + } + } + fusedDMaps[idxImage] = true; + ASSERT(pointcloud.points.size() == pointcloud.pointViews.size() && pointcloud.points.size() == pointcloud.pointWeights.size()); + DEBUG_ULTIMATE("Depths map for reference image %3u fused using %u depths maps: %u new points, %u/%u cached images (%s)", + idxImage, depthData.images.size() - 1, pointcloud.points.size() - nNumPointsPrev, numImageNeighborsInCache, numImagesInCache, TD_TIMER_GET_FMT().c_str()); + progress.display(numDMapsFused); + // ensure enough memory is available for the next depth-maps chunk + cacheDMaps.SkipMemoryCheckIdxImage(); + if (numDMapsFused % numDMapsReserveFusion == 0) + cacheDMaps.SetMaxMemory(GetAvailableMemory(arrDepthData, fusedDMaps, numDMapsReserveFusion, cacheDMaps.GetUsedMemory())); + } + GET_LOGCONSOLE().Play(); + progress.close(); + arrUseMask.Release(); + cacheDMaps.ClearCache(); + if (!_bEstimateNormal) + pointcloud.normals.Release(); + + DEBUG_EXTRA("Depth-maps dense fused and filtered: %u depth-maps, %u depths, %u points (%d%%%%), %.2f hits in %.2f cached (%s)", + numDMapsFused, nDepths, pointcloud.points.size(), ROUND2INT((100.f*pointcloud.points.size())/nDepths), + static_cast(totalNumImageNeighborsInCache) / numDMapsFused, + static_cast(totalNumImagesInCache) / numDMapsFused, TD_TIMER_GET_FMT().c_str()); +} // DenseFuseDepthMaps /*----------------------------------------------------------------*/ @@ -1645,12 +1900,18 @@ bool Scene::DenseReconstruction(int nFusionMode, bool bCrop2ROI, float fBorderRO // fuse all depth-maps pointcloud.Release(); - if (OPTDENSE::nMinViewsFuse < 2) { + switch (OPTDENSE::nFuseFilter) { + case OPTDENSE::FUSE_NOFILTER: // merge depth-maps data.depthMaps.MergeDepthMaps(pointcloud, OPTDENSE::nEstimateColors == 2, OPTDENSE::nEstimateNormals == 2); - } else { + break; + case OPTDENSE::FUSE_FILTER: // fuse depth-maps data.depthMaps.FuseDepthMaps(pointcloud, OPTDENSE::nEstimateColors == 2, OPTDENSE::nEstimateNormals == 2); + break; + case OPTDENSE::FUSE_DENSEFILTER: + // dense fuse depth-maps + data.depthMaps.DenseFuseDepthMaps(pointcloud, OPTDENSE::nEstimateColors == 2, OPTDENSE::nEstimateNormals == 2); } #if TD_VERBOSE != TD_VERBOSE_OFF if (g_nVerbosityLevel > 2) { diff --git a/libs/MVS/SceneDensify.h b/libs/MVS/SceneDensify.h index 74f29764a..7934ed8e9 100644 --- a/libs/MVS/SceneDensify.h +++ b/libs/MVS/SceneDensify.h @@ -68,6 +68,7 @@ class MVS_API DepthMapsData bool FilterDepthMap(DepthData& depthData, const IIndexArr& idxNeighbors, bool bAdjust=true); void MergeDepthMaps(PointCloud& pointcloud, bool bEstimateColor, bool bEstimateNormal); void FuseDepthMaps(PointCloud& pointcloud, bool bEstimateColor, bool bEstimateNormal); + void DenseFuseDepthMaps(PointCloud& pointcloud, bool bEstimateColor, bool bEstimateNormal); static DepthData ScaleDepthData(const DepthData& inputDeptData, float scale);