From d99edd50d5871c67b1565e3435bf8fc87361e55d Mon Sep 17 00:00:00 2001 From: Abdelrahman AL MAROUK Date: Wed, 2 Aug 2023 11:47:37 +0200 Subject: [PATCH 01/37] [prepareDenseScene] add landmarks-based mask for depth maps --- .../pipeline/main_prepareDenseScene.cpp | 74 ++++++++++++++----- 1 file changed, 54 insertions(+), 20 deletions(-) diff --git a/src/software/pipeline/main_prepareDenseScene.cpp b/src/software/pipeline/main_prepareDenseScene.cpp index bd15188639..f88c0c572a 100644 --- a/src/software/pipeline/main_prepareDenseScene.cpp +++ b/src/software/pipeline/main_prepareDenseScene.cpp @@ -92,7 +92,8 @@ bool prepareDenseScene(const SfMData& sfmData, image::EImageFileType outputFileType, bool saveMetadata, bool saveMatricesFiles, - bool evCorrection) + bool evCorrection, + float landmarksMaskScale) { // defined view Ids std::set viewIds; @@ -128,6 +129,8 @@ bool prepareDenseScene(const SfMData& sfmData, const double medianCameraExposure = sfmData.getMedianCameraExposureSetting().getExposure(); ALICEVISION_LOG_INFO("Median Camera Exposure: " << medianCameraExposure << ", Median EV: " << std::log2(1.0 / medianCameraExposure)); + const LandmarksPerView landmarksPerView = getLandmarksPerViews(sfmData); + #pragma omp parallel for num_threads(3) for (int i = 0; i < viewIds.size(); ++i) { @@ -251,29 +254,55 @@ bool prepareDenseScene(const SfMData& sfmData, << " Ev compensation: " + std::to_string(exposureCompensation)); } - image::Image mask; - if (tryLoadMask(&mask, masksFolders, viewId, srcImage, maskExtension)) + image::Image maskLandmarks; + if(landmarksMaskScale > 0.f) { - process>( - dstColorImage, cam, metadata, srcImage, evCorrection, exposureCompensation, [&mask](Image& image) { - if (image.width() * image.height() != mask.width() * mask.height()) - { - ALICEVISION_LOG_WARNING("Invalid image mask size: mask is ignored."); - return; - } - - for (int pix = 0; pix < image.width() * image.height(); ++pix) - { - const bool masked = (mask(pix) == 0); - image(pix).a() = masked ? 0.f : 1.f; - } - }); + const LandmarkIdSet& landmarksSet = landmarksPerView.at(viewId); + image::Image image; + readImage(srcImage, image, image::EImageColorSpace::LINEAR); + maskLandmarks = image::Image(image.Width(), image.Height(), true, 0); + int r = (int)(landmarksMaskScale * 0.5f * (image.Width() + image.Height())); + for(const auto& landmarkId : landmarksSet) + { + const sfmData::Landmark& landmark = sfmData.getLandmarks().at(landmarkId); + Observation obs = landmark.observations.at(viewId); + // maskLandmarks(obs.x.y(), obs.x.x()) = 1; + for(int y = std::max(obs.x.y() - r, 0.); + y <= std::min(obs.x.y() + r, (double)maskLandmarks.Height() - 1); y++) + { + for(int x = std::max(obs.x.x() - r, 0.); + x <= std::min(obs.x.x() + r, (double)maskLandmarks.Width() - 1); x++) + { + maskLandmarks(y, x) = 1; + } + } + } } + + image::Image mask; + if(tryLoadMask(&mask, masksFolders, viewId, srcImage, maskExtension)) + mask = mask * maskLandmarks; else + mask = maskLandmarks; + process>(dstColorImage, cam, metadata, srcImage, evCorrection, exposureCompensation, [&mask] (Image & image) + { + if(image.Width() * image.Height() != mask.Width() * mask.Height()) + { + ALICEVISION_LOG_WARNING("Invalid image mask size: mask is ignored."); + return; + } + + for(int pix = 0; pix < image.Width() * image.Height(); ++pix) + { + const bool masked = (mask(pix) == 0); + image(pix).a() = masked ? 0.f : 1.f; + } + }); + /* else { const auto noMaskingFunc = [](Image& image) {}; process>(dstColorImage, cam, metadata, srcImage, evCorrection, exposureCompensation, noMaskingFunc); - } + }*/ } ++progressDisplay; @@ -298,6 +327,7 @@ int aliceVision_main(int argc, char* argv[]) bool saveMetadata = true; bool saveMatricesTxtFiles = false; bool evCorrection = false; + float landmarksMaskScale = 0.f; // clang-format off po::options_description requiredParams("Required parameters"); @@ -328,7 +358,10 @@ int aliceVision_main(int argc, char* argv[]) ("rangeSize", po::value(&rangeSize)->default_value(rangeSize), "Range size.") ("evCorrection", po::value(&evCorrection)->default_value(evCorrection), - "Correct exposure value."); + "Correct exposure value.") + ("landmarksMaskScale", po::value(&landmarksMaskScale)->default_value(landmarksMaskScale), + "Scale (relative to image size) of the projection of landmarks to mask images for depth computation.\n" + "If 0, masking using landmarks will not be used."); // clang-format on CmdLine cmdline("AliceVision prepareDenseScene"); @@ -392,7 +425,8 @@ int aliceVision_main(int argc, char* argv[]) outputFileType, saveMetadata, saveMatricesTxtFiles, - evCorrection)) + evCorrection, + landmarksMaskScale)) return EXIT_SUCCESS; return EXIT_FAILURE; From feaa520dd73415758d8c03e8c31d47e6621286de Mon Sep 17 00:00:00 2001 From: Abdelrahman AL MAROUK Date: Mon, 31 Jul 2023 09:46:07 +0100 Subject: [PATCH 02/37] [filterSfM] add main for filtering Sfm data --- src/software/pipeline/CMakeLists.txt | 12 ++ src/software/pipeline/main_filterSfM.cpp | 144 +++++++++++++++++++++++ 2 files changed, 156 insertions(+) create mode 100644 src/software/pipeline/main_filterSfM.cpp diff --git a/src/software/pipeline/CMakeLists.txt b/src/software/pipeline/CMakeLists.txt index 507162cfaf..7daa29a57e 100644 --- a/src/software/pipeline/CMakeLists.txt +++ b/src/software/pipeline/CMakeLists.txt @@ -330,6 +330,18 @@ if(ALICEVISION_BUILD_SFM) Boost::boost Boost::timer ) + + # Filter SfM data + alicevision_add_software(aliceVision_filterSfM + SOURCE main_filterSfM.cpp + FOLDER ${FOLDER_SOFTWARE_PIPELINE} + LINKS aliceVision_system + aliceVision_cmdline + aliceVision_sfmData + aliceVision_sfmDataIO + Boost::program_options + Boost::filesystem + ) endif() # if(ALICEVISION_BUILD_SFM) diff --git a/src/software/pipeline/main_filterSfM.cpp b/src/software/pipeline/main_filterSfM.cpp new file mode 100644 index 0000000000..b77a65c8ed --- /dev/null +++ b/src/software/pipeline/main_filterSfM.cpp @@ -0,0 +1,144 @@ +// This file is part of the AliceVision project. +// Copyright (c) 2017 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 +#include +#include + +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include + +// These constants define the current software version. +// They must be updated when the command line is changed. +#define ALICEVISION_SOFTWARE_VERSION_MAJOR 1 +#define ALICEVISION_SOFTWARE_VERSION_MINOR 0 + +using namespace aliceVision; +using namespace aliceVision::camera; +using namespace aliceVision::geometry; +using namespace aliceVision::image; +using namespace aliceVision::sfmData; +using namespace aliceVision::sfmDataIO; + +namespace po = boost::program_options; +namespace fs = boost::filesystem; + +bool filterSfMData(SfMData& sfmData, int maxNbObservationsPerLandmark) +{ + #pragma omp parallel for + for(auto i = 0; i < sfmData.getLandmarks().size(); i++) + { + auto landmarkPair = sfmData.getLandmarks().begin(); + std::advance(landmarkPair, i); + sfmData::Landmark& landmark = landmarkPair->second; + + // check number of observations + if(landmark.observations.size() <= maxNbObservationsPerLandmark) + continue; + + // // check angle between observations + // if(!checkLandmarkMinObservationAngle(sfmData, landmark, minObservationAngle)) + + // compute observation scores + + std::vector > observationScores; + observationScores.reserve(landmark.observations.size()); + + for(const auto& observationPair : landmark.observations) + { + const IndexT viewId = observationPair.first; + const sfmData::View& view = *(sfmData.getViews().at(viewId)); + const geometry::Pose3 pose = sfmData.getPose(view).getTransform(); + + observationScores.push_back(std::pair( + (pose.center() - landmark.X).squaredNorm(), + viewId + )); + } + + // sort observations by ascending score order + std::stable_sort(observationScores.begin(), observationScores.end()); + // take only best observations + observationScores.resize(maxNbObservationsPerLandmark); + + // replace the observations + Observations filteredObservations; + for(auto observationScorePair : observationScores) + { + filteredObservations[observationScorePair.second] = landmark.observations[observationScorePair.second]; + } + landmark.observations = filteredObservations; + } + return true; +} + +int aliceVision_main(int argc, char *argv[]) +{ + // command-line parameters + + // std::string verboseLevel = system::EVerboseLevel_enumToString(system::Logger::getDefaultVerboseLevel()); + std::string inputSfmFilename; + std::string outputSfmFilename; + int maxNbObservationsPerLandmark = 5; + + po::options_description requiredParams("Required parameters"); + requiredParams.add_options() + ("input,i", po::value(&inputSfmFilename)->required(), + "Input SfMData file.") + ("output,o", po::value(&outputSfmFilename)->required(), + "Output SfMData file."); + + po::options_description optionalParams("Optional parameters"); + optionalParams.add_options() + ("maxNbObservationsPerLandmark", po::value(&maxNbObservationsPerLandmark)->default_value(maxNbObservationsPerLandmark), + "Maximum number of allowed observations per landmark."); + + CmdLine cmdline("AliceVision SfM filtering."); + cmdline.add(requiredParams); + cmdline.add(optionalParams); + if (!cmdline.execute(argc, argv)) + { + return EXIT_FAILURE; + } + + // Create output dir + { + auto outDir = fs::path(outputSfmFilename).parent_path().string(); + if (!fs::exists(outDir)) + fs::create_directory(outDir); + } + + // Read the input SfM scene + SfMData sfmData; + if(!sfmDataIO::Load(sfmData, inputSfmFilename, sfmDataIO::ESfMData::ALL)) + { + ALICEVISION_LOG_ERROR("The input SfMData file '" << inputSfmFilename << "' cannot be read."); + return EXIT_FAILURE; + } + + // filter SfM data + if(filterSfMData(sfmData, maxNbObservationsPerLandmark)) + { + sfmDataIO::Save(sfmData, outputSfmFilename, sfmDataIO::ESfMData::ALL); + return EXIT_SUCCESS; + } + + return EXIT_FAILURE; +} From 6aa4ef3cc84715e595adecc5c2b17af194cf4d1c Mon Sep 17 00:00:00 2001 From: Abdelrahman AL MAROUK Date: Thu, 10 Aug 2023 10:12:24 +0200 Subject: [PATCH 03/37] [filterSfM] add support for filtering landmarks --- src/software/pipeline/CMakeLists.txt | 24 +-- src/software/pipeline/main_filterSfM.cpp | 183 ++++++++++++++++++++++- 2 files changed, 193 insertions(+), 14 deletions(-) diff --git a/src/software/pipeline/CMakeLists.txt b/src/software/pipeline/CMakeLists.txt index 7daa29a57e..040eec66ad 100644 --- a/src/software/pipeline/CMakeLists.txt +++ b/src/software/pipeline/CMakeLists.txt @@ -331,17 +331,19 @@ if(ALICEVISION_BUILD_SFM) Boost::timer ) - # Filter SfM data - alicevision_add_software(aliceVision_filterSfM - SOURCE main_filterSfM.cpp - FOLDER ${FOLDER_SOFTWARE_PIPELINE} - LINKS aliceVision_system - aliceVision_cmdline - aliceVision_sfmData - aliceVision_sfmDataIO - Boost::program_options - Boost::filesystem - ) + # Filter SfM data + alicevision_add_software(aliceVision_filterSfM + SOURCE main_filterSfM.cpp + FOLDER ${FOLDER_SOFTWARE_PIPELINE} + LINKS aliceVision_system + aliceVision_cmdline + aliceVision_sfmData + aliceVision_mvsUtils + aliceVision_sfmDataIO + Boost::program_options + Boost::filesystem + nanoflann + ) endif() # if(ALICEVISION_BUILD_SFM) diff --git a/src/software/pipeline/main_filterSfM.cpp b/src/software/pipeline/main_filterSfM.cpp index b77a65c8ed..be76d68396 100644 --- a/src/software/pipeline/main_filterSfM.cpp +++ b/src/software/pipeline/main_filterSfM.cpp @@ -11,6 +11,9 @@ #include #include #include +#include + +#include "nanoflann.hpp" #include #include @@ -40,7 +43,171 @@ using namespace aliceVision::sfmDataIO; namespace po = boost::program_options; namespace fs = boost::filesystem; -bool filterSfMData(SfMData& sfmData, int maxNbObservationsPerLandmark) +static const std::size_t MAX_LEAF_ELEMENTS = 10; + +struct LandmarksAdaptator +{ + using Derived = LandmarksAdaptator; //!< In this case the dataset class is myself. + using T = double; + + /// CRTP helper method + inline const Derived& derived() const { return *static_cast(this); } + /// CRTP helper method + inline Derived& derived() { return *static_cast(this); } + + const sfmData::Landmarks& _data; + LandmarksAdaptator(const sfmData::Landmarks& data) + : _data(data) + { + } + + // Must return the number of data points + inline size_t kdtree_get_point_count() const { return _data.size(); } + + // Returns the dim'th component of the idx'th point in the class: + inline T kdtree_get_pt(const size_t idx, int dim) const + { + auto it = _data.begin(); + std::advance(it, idx); + return it->second.X(dim); + } + + // Optional bounding-box computation: return false to default to a standard bbox computation loop. + // Return true if the BBOX was already computed by the class and returned in "bb" so it can be avoided to redo it + // again. Look at bb.size() to find out the expected dimensionality (e.g. 2 or 3 for point clouds) + template + bool kdtree_get_bbox(BBOX& bb) const + { + return false; + } +}; + +using KdTree = nanoflann::KDTreeSingleIndexAdaptor< + nanoflann::L2_Simple_Adaptor, + LandmarksAdaptator, + 3, /* dim */ + IndexT +>; + +/** + * A result-set class used when performing a radius based search. + */ +class PixSizeSearch +{ +public: + const double _radius_sq; + const std::vector& _pixSize; + const double _pixSize_i; + bool found = false; + + inline PixSizeSearch(double radius, const std::vector& pixSize, int i) + : _radius_sq(radius * radius) + , _pixSize(pixSize) + , _pixSize_i(pixSize[i]) + { + } + + inline bool full() const { return found; } + + inline bool addPoint(double dist, IndexT index) + { + if(dist < _radius_sq && _pixSize[index] < _pixSize_i) + { + found = true; + return false; + } + return true; + } + + inline double worstDist() const { return _radius_sq; } +}; + +bool filterLandmarks(SfMData& sfmData, double radiusScale) +{ + mvsUtils::MultiViewParams mp(sfmData, "", "", "", false); + std::vector landmarksPixSize(sfmData.getLandmarks().size()); + + ALICEVISION_LOG_INFO("Computing pixel size: started."); + #pragma omp parallel for + for(auto i = 0; i < sfmData.getLandmarks().size(); i++) + { + auto landmarkPair = sfmData.getLandmarks().begin(); + std::advance(landmarkPair, i); + const sfmData::Landmark& landmark = landmarkPair->second; + + // compute landmark pixSize + double pixSize = 0.; + int n = 0; + for(const auto& observationPair : landmark.observations) + { + const IndexT viewId = observationPair.first; + pixSize += mp.getCamPixelSize(Point3d(landmark.X.x(), landmark.X.y(), landmark.X.z()), + mp.getIndexFromViewId(viewId), observationPair.second.scale); + n++; + } + pixSize /= n; + + landmarksPixSize[i] = pixSize; + } + ALICEVISION_LOG_INFO("Computing pixel size: done."); + + //// sort landmarks by descending pixSize order + //std::stable_sort(landmarksPixSize.begin(), landmarksPixSize.end(), std::greater<>{}); + + ALICEVISION_LOG_INFO("Build nanoflann KdTree index."); + LandmarksAdaptator data(sfmData.getLandmarks()); + KdTree tree(3, data, nanoflann::KDTreeSingleIndexAdaptorParams(10)); + tree.buildIndex(); + ALICEVISION_LOG_INFO("KdTree created for " << sfmData.getLandmarks().size() << " points."); + std::vector newIdx(sfmData.getLandmarks().size()); + IndexT currentIdx = 0; + + ALICEVISION_LOG_INFO("Identifying landmarks to remove: started."); + #pragma omp parallel for + for (auto i = 0; i < sfmData.getLandmarks().size(); i++) + { + PixSizeSearch search(landmarksPixSize[i] * radiusScale, landmarksPixSize, i); + bool found = tree.findNeighbors(search, sfmData.getLandmarks().at(i).X.data(), nanoflann::SearchParams()); + if (!found) + { + newIdx[i] = -1; + } + else + { + newIdx[i] = currentIdx++; + } + } + ALICEVISION_LOG_INFO("Identifying landmarks to remove: done."); + + ALICEVISION_LOG_INFO("Removing landmarks: started."); + Landmarks filteredLandmarks; + #pragma omp parallel for + for (auto i = 0; i < sfmData.getLandmarks().size(); i++) + { + if(newIdx[i] != -1) + { + filteredLandmarks[newIdx[i]] = sfmData.getLandmarks().at(i); + } + } + sfmData.getLandmarks() = filteredLandmarks; + ALICEVISION_LOG_INFO("Removing landmarks: done."); + + //// take only best observations + //observationScores.resize(maxNbObservationsPerLandmark); + + + + //// replace the observations + //Observations filteredObservations; + //for(auto observationScorePair : observationScores) + //{ + // filteredObservations[observationScorePair.second] = landmark.observations[observationScorePair.second]; + //} + //landmark.observations = filteredObservations; + return true; +} + +bool filterObservations(SfMData& sfmData, int maxNbObservationsPerLandmark) { #pragma omp parallel for for(auto i = 0; i < sfmData.getLandmarks().size(); i++) @@ -97,6 +264,7 @@ int aliceVision_main(int argc, char *argv[]) std::string inputSfmFilename; std::string outputSfmFilename; int maxNbObservationsPerLandmark = 5; + double radiusScale = 2; po::options_description requiredParams("Required parameters"); requiredParams.add_options() @@ -108,7 +276,9 @@ int aliceVision_main(int argc, char *argv[]) po::options_description optionalParams("Optional parameters"); optionalParams.add_options() ("maxNbObservationsPerLandmark", po::value(&maxNbObservationsPerLandmark)->default_value(maxNbObservationsPerLandmark), - "Maximum number of allowed observations per landmark."); + "Maximum number of allowed observations per landmark.") + ("radiusScale", po::value(&radiusScale)->default_value(radiusScale), + "Scale factor applied to pixel size based radius filter applied to landmarks."); CmdLine cmdline("AliceVision SfM filtering."); cmdline.add(requiredParams); @@ -134,7 +304,14 @@ int aliceVision_main(int argc, char *argv[]) } // filter SfM data - if(filterSfMData(sfmData, maxNbObservationsPerLandmark)) + ALICEVISION_LOG_INFO("Filtering landmarks: started."); + bool success2 = filterLandmarks(sfmData, radiusScale); + ALICEVISION_LOG_INFO("Filtering landmarks: done."); + ALICEVISION_LOG_INFO("Filtering observations: started."); + bool success1 = filterObservations(sfmData, maxNbObservationsPerLandmark); + ALICEVISION_LOG_INFO("Filtering observations: done."); + + if(success1) { sfmDataIO::Save(sfmData, outputSfmFilename, sfmDataIO::ESfMData::ALL); return EXIT_SUCCESS; From 5ab39e0f466b07eb1d396c1de56d4858fc77ea57 Mon Sep 17 00:00:00 2001 From: Abdelrahman AL MAROUK Date: Fri, 11 Aug 2023 10:20:32 +0200 Subject: [PATCH 04/37] [filterSfM] fix bug - OpenMP loop race problem - inverse if condition --- src/software/pipeline/main_filterSfM.cpp | 19 +++++++++++++------ 1 file changed, 13 insertions(+), 6 deletions(-) diff --git a/src/software/pipeline/main_filterSfM.cpp b/src/software/pipeline/main_filterSfM.cpp index be76d68396..f1fcaf7744 100644 --- a/src/software/pipeline/main_filterSfM.cpp +++ b/src/software/pipeline/main_filterSfM.cpp @@ -43,7 +43,7 @@ using namespace aliceVision::sfmDataIO; namespace po = boost::program_options; namespace fs = boost::filesystem; -static const std::size_t MAX_LEAF_ELEMENTS = 10; +static const std::size_t MAX_LEAF_ELEMENTS = 64; struct LandmarksAdaptator { @@ -156,7 +156,7 @@ bool filterLandmarks(SfMData& sfmData, double radiusScale) ALICEVISION_LOG_INFO("Build nanoflann KdTree index."); LandmarksAdaptator data(sfmData.getLandmarks()); - KdTree tree(3, data, nanoflann::KDTreeSingleIndexAdaptorParams(10)); + KdTree tree(3, data, nanoflann::KDTreeSingleIndexAdaptorParams(MAX_LEAF_ELEMENTS)); tree.buildIndex(); ALICEVISION_LOG_INFO("KdTree created for " << sfmData.getLandmarks().size() << " points."); std::vector newIdx(sfmData.getLandmarks().size()); @@ -168,7 +168,7 @@ bool filterLandmarks(SfMData& sfmData, double radiusScale) { PixSizeSearch search(landmarksPixSize[i] * radiusScale, landmarksPixSize, i); bool found = tree.findNeighbors(search, sfmData.getLandmarks().at(i).X.data(), nanoflann::SearchParams()); - if (!found) + if (found) { newIdx[i] = -1; } @@ -177,19 +177,26 @@ bool filterLandmarks(SfMData& sfmData, double radiusScale) newIdx[i] = currentIdx++; } } + + ALICEVISION_LOG_INFO( + "Identified " << (sfmData.getLandmarks().size() - currentIdx) << + " landmarks to remove out of " << (sfmData.getLandmarks().size()) << + ", i.e. " << ((sfmData.getLandmarks().size() - currentIdx) * 100.f / sfmData.getLandmarks().size()) << + " % " + ); ALICEVISION_LOG_INFO("Identifying landmarks to remove: done."); ALICEVISION_LOG_INFO("Removing landmarks: started."); - Landmarks filteredLandmarks; + std::vector> filteredLandmarks(currentIdx); #pragma omp parallel for for (auto i = 0; i < sfmData.getLandmarks().size(); i++) { if(newIdx[i] != -1) { - filteredLandmarks[newIdx[i]] = sfmData.getLandmarks().at(i); + filteredLandmarks[newIdx[i]] = std::make_pair(newIdx[i], sfmData.getLandmarks().at(i)); } } - sfmData.getLandmarks() = filteredLandmarks; + sfmData.getLandmarks() = Landmarks(filteredLandmarks.begin(), filteredLandmarks.end()); ALICEVISION_LOG_INFO("Removing landmarks: done."); //// take only best observations From 31d10b5ce34b96a34c7ba954a766ee3de353cb38 Mon Sep 17 00:00:00 2001 From: Abdelrahman AL MAROUK Date: Fri, 11 Aug 2023 17:45:07 +0200 Subject: [PATCH 05/37] [filterSfM] optimize code for performance --- src/software/pipeline/main_filterSfM.cpp | 136 ++++++++++++----------- 1 file changed, 74 insertions(+), 62 deletions(-) diff --git a/src/software/pipeline/main_filterSfM.cpp b/src/software/pipeline/main_filterSfM.cpp index f1fcaf7744..98e39cc469 100644 --- a/src/software/pipeline/main_filterSfM.cpp +++ b/src/software/pipeline/main_filterSfM.cpp @@ -55,8 +55,8 @@ struct LandmarksAdaptator /// CRTP helper method inline Derived& derived() { return *static_cast(this); } - const sfmData::Landmarks& _data; - LandmarksAdaptator(const sfmData::Landmarks& data) + const std::vector _data; + LandmarksAdaptator(const std::vector& data) : _data(data) { } @@ -67,9 +67,7 @@ struct LandmarksAdaptator // Returns the dim'th component of the idx'th point in the class: inline T kdtree_get_pt(const size_t idx, int dim) const { - auto it = _data.begin(); - std::advance(it, idx); - return it->second.X(dim); + return _data[idx].X(dim); } // Optional bounding-box computation: return false to default to a standard bbox computation loop. @@ -100,7 +98,7 @@ class PixSizeSearch const double _pixSize_i; bool found = false; - inline PixSizeSearch(double radius, const std::vector& pixSize, int i) + inline PixSizeSearch(double radius, const std::vector& pixSize, size_t i) : _radius_sq(radius * radius) , _pixSize(pixSize) , _pixSize_i(pixSize[i]) @@ -124,16 +122,23 @@ class PixSizeSearch bool filterLandmarks(SfMData& sfmData, double radiusScale) { + std::vector landmarksData(sfmData.getLandmarks().size()); + { + size_t i = 0; + for(const auto& it : sfmData.getLandmarks()) + { + landmarksData[i++] = it.second; + } + } + mvsUtils::MultiViewParams mp(sfmData, "", "", "", false); - std::vector landmarksPixSize(sfmData.getLandmarks().size()); + std::vector landmarksPixSize(landmarksData.size()); ALICEVISION_LOG_INFO("Computing pixel size: started."); #pragma omp parallel for - for(auto i = 0; i < sfmData.getLandmarks().size(); i++) + for(auto i = 0; i < landmarksData.size(); i++) { - auto landmarkPair = sfmData.getLandmarks().begin(); - std::advance(landmarkPair, i); - const sfmData::Landmark& landmark = landmarkPair->second; + const Landmark& landmark = landmarksData[i]; // compute landmark pixSize double pixSize = 0.; @@ -155,73 +160,66 @@ bool filterLandmarks(SfMData& sfmData, double radiusScale) //std::stable_sort(landmarksPixSize.begin(), landmarksPixSize.end(), std::greater<>{}); ALICEVISION_LOG_INFO("Build nanoflann KdTree index."); - LandmarksAdaptator data(sfmData.getLandmarks()); + LandmarksAdaptator data(landmarksData); KdTree tree(3, data, nanoflann::KDTreeSingleIndexAdaptorParams(MAX_LEAF_ELEMENTS)); tree.buildIndex(); - ALICEVISION_LOG_INFO("KdTree created for " << sfmData.getLandmarks().size() << " points."); - std::vector newIdx(sfmData.getLandmarks().size()); - IndexT currentIdx = 0; + ALICEVISION_LOG_INFO("KdTree created for " << landmarksData.size() << " points."); + std::vector toRemove(landmarksData.size(), false); ALICEVISION_LOG_INFO("Identifying landmarks to remove: started."); - #pragma omp parallel for - for (auto i = 0; i < sfmData.getLandmarks().size(); i++) + + size_t nbToRemove = 0; + #pragma omp parallel for reduction(+:nbToRemove) + for(auto i = 0; i < landmarksData.size(); i++) { PixSizeSearch search(landmarksPixSize[i] * radiusScale, landmarksPixSize, i); - bool found = tree.findNeighbors(search, sfmData.getLandmarks().at(i).X.data(), nanoflann::SearchParams()); + bool found = tree.findNeighbors(search, landmarksData[i].X.data(), nanoflann::SearchParams()); if (found) { - newIdx[i] = -1; - } - else - { - newIdx[i] = currentIdx++; + toRemove[i] = true; + nbToRemove++; } } ALICEVISION_LOG_INFO( - "Identified " << (sfmData.getLandmarks().size() - currentIdx) << - " landmarks to remove out of " << (sfmData.getLandmarks().size()) << - ", i.e. " << ((sfmData.getLandmarks().size() - currentIdx) * 100.f / sfmData.getLandmarks().size()) << + "Identified " << (landmarksData.size() - nbToRemove) << + " landmarks to remove out of " << (landmarksData.size()) << + ", i.e. " << ((landmarksData.size() - nbToRemove) * 100.f / landmarksData.size()) << " % " ); ALICEVISION_LOG_INFO("Identifying landmarks to remove: done."); ALICEVISION_LOG_INFO("Removing landmarks: started."); - std::vector> filteredLandmarks(currentIdx); - #pragma omp parallel for - for (auto i = 0; i < sfmData.getLandmarks().size(); i++) + std::vector> filteredLandmarks(landmarksData.size() - nbToRemove); + IndexT newIdx = 0; + for (size_t i = 0; i < landmarksData.size(); i++) { - if(newIdx[i] != -1) + if(!toRemove[i]) { - filteredLandmarks[newIdx[i]] = std::make_pair(newIdx[i], sfmData.getLandmarks().at(i)); + filteredLandmarks[newIdx++] = std::make_pair(newIdx, landmarksData[i]); } } sfmData.getLandmarks() = Landmarks(filteredLandmarks.begin(), filteredLandmarks.end()); ALICEVISION_LOG_INFO("Removing landmarks: done."); - //// take only best observations - //observationScores.resize(maxNbObservationsPerLandmark); - - - - //// replace the observations - //Observations filteredObservations; - //for(auto observationScorePair : observationScores) - //{ - // filteredObservations[observationScorePair.second] = landmark.observations[observationScorePair.second]; - //} - //landmark.observations = filteredObservations; return true; } bool filterObservations(SfMData& sfmData, int maxNbObservationsPerLandmark) { + std::vector landmarksData(sfmData.getLandmarks().size()); + { + size_t i = 0; + for(const auto& it : sfmData.getLandmarks()) + { + landmarksData[i++] = it.second; + } + } + #pragma omp parallel for - for(auto i = 0; i < sfmData.getLandmarks().size(); i++) + for(auto i = 0; i < landmarksData.size(); i++) { - auto landmarkPair = sfmData.getLandmarks().begin(); - std::advance(landmarkPair, i); - sfmData::Landmark& landmark = landmarkPair->second; + sfmData::Landmark landmark = landmarksData[i]; // check number of observations if(landmark.observations.size() <= maxNbObservationsPerLandmark) @@ -232,8 +230,7 @@ bool filterObservations(SfMData& sfmData, int maxNbObservationsPerLandmark) // compute observation scores - std::vector > observationScores; - observationScores.reserve(landmark.observations.size()); + std::vector> observationScores(landmark.observations.size()); for(const auto& observationPair : landmark.observations) { @@ -254,7 +251,7 @@ bool filterObservations(SfMData& sfmData, int maxNbObservationsPerLandmark) // replace the observations Observations filteredObservations; - for(auto observationScorePair : observationScores) + for(const auto& observationScorePair : observationScores) { filteredObservations[observationScorePair.second] = landmark.observations[observationScorePair.second]; } @@ -273,6 +270,11 @@ int aliceVision_main(int argc, char *argv[]) int maxNbObservationsPerLandmark = 5; double radiusScale = 2; + // user optional parameters + std::vector featuresFolders; + std::vector matchesFolders; + std::string describerTypesName = feature::EImageDescriberType_enumToString(feature::EImageDescriberType::SIFT); + po::options_description requiredParams("Required parameters"); requiredParams.add_options() ("input,i", po::value(&inputSfmFilename)->required(), @@ -285,7 +287,13 @@ int aliceVision_main(int argc, char *argv[]) ("maxNbObservationsPerLandmark", po::value(&maxNbObservationsPerLandmark)->default_value(maxNbObservationsPerLandmark), "Maximum number of allowed observations per landmark.") ("radiusScale", po::value(&radiusScale)->default_value(radiusScale), - "Scale factor applied to pixel size based radius filter applied to landmarks."); + "Scale factor applied to pixel size based radius filter applied to landmarks.") + ("featuresFolders,f", po::value>(&featuresFolders)->multitoken(), + "Path to folder(s) containing the extracted features.") + ("matchesFolders,m", po::value>(&matchesFolders)->multitoken(), + "Path to folder(s) in which computed matches are stored.") + ("describerTypes,d", po::value(&describerTypesName)->default_value(describerTypesName), + feature::EImageDescriberType_informations().c_str()); CmdLine cmdline("AliceVision SfM filtering."); cmdline.add(requiredParams); @@ -311,18 +319,22 @@ int aliceVision_main(int argc, char *argv[]) } // filter SfM data - ALICEVISION_LOG_INFO("Filtering landmarks: started."); - bool success2 = filterLandmarks(sfmData, radiusScale); - ALICEVISION_LOG_INFO("Filtering landmarks: done."); - ALICEVISION_LOG_INFO("Filtering observations: started."); - bool success1 = filterObservations(sfmData, maxNbObservationsPerLandmark); - ALICEVISION_LOG_INFO("Filtering observations: done."); - - if(success1) + if(radiusScale > 0) + { + ALICEVISION_LOG_INFO("Filtering landmarks: started."); + filterLandmarks(sfmData, radiusScale); + ALICEVISION_LOG_INFO("Filtering landmarks: done."); + } + + if(maxNbObservationsPerLandmark > 0) { - sfmDataIO::Save(sfmData, outputSfmFilename, sfmDataIO::ESfMData::ALL); - return EXIT_SUCCESS; + ALICEVISION_LOG_INFO("Filtering observations: started."); + filterObservations(sfmData, maxNbObservationsPerLandmark); + ALICEVISION_LOG_INFO("Filtering observations: done."); } - return EXIT_FAILURE; + + sfmDataIO::Save(sfmData, outputSfmFilename, sfmDataIO::ESfMData::ALL); + return EXIT_SUCCESS; + } From ea9a1537969956fee7e932ec207c357e76844f9b Mon Sep 17 00:00:00 2001 From: Abdelrahman AL MAROUK Date: Mon, 14 Aug 2023 07:41:46 +0200 Subject: [PATCH 06/37] [prepareDenseScene] fix landmarks-based mask for depth maps - put masked pixels to 0.5 instead of 0 in order to be taken into consideration for Tc cameras - fix index out of range thrown in case there are no observations for a certain view --- .../pipeline/main_prepareDenseScene.cpp | 61 ++++++++++--------- 1 file changed, 33 insertions(+), 28 deletions(-) diff --git a/src/software/pipeline/main_prepareDenseScene.cpp b/src/software/pipeline/main_prepareDenseScene.cpp index f88c0c572a..d3ef386f06 100644 --- a/src/software/pipeline/main_prepareDenseScene.cpp +++ b/src/software/pipeline/main_prepareDenseScene.cpp @@ -129,7 +129,7 @@ bool prepareDenseScene(const SfMData& sfmData, const double medianCameraExposure = sfmData.getMedianCameraExposureSetting().getExposure(); ALICEVISION_LOG_INFO("Median Camera Exposure: " << medianCameraExposure << ", Median EV: " << std::log2(1.0 / medianCameraExposure)); - const LandmarksPerView landmarksPerView = getLandmarksPerViews(sfmData); + const LandmarksPerView& landmarksPerView = getLandmarksPerViews(sfmData); #pragma omp parallel for num_threads(3) for (int i = 0; i < viewIds.size(); ++i) @@ -255,49 +255,54 @@ bool prepareDenseScene(const SfMData& sfmData, } image::Image maskLandmarks; - if(landmarksMaskScale > 0.f) + bool doMaskLandmarks = landmarksMaskScale > 0.f; + if(doMaskLandmarks) { - const LandmarkIdSet& landmarksSet = landmarksPerView.at(viewId); image::Image image; readImage(srcImage, image, image::EImageColorSpace::LINEAR); - maskLandmarks = image::Image(image.Width(), image.Height(), true, 0); + // for the T camera, image alpha should be at least 0.4f * 255 (masking) + maskLandmarks = image::Image(image.Width(), image.Height(), true, 127); int r = (int)(landmarksMaskScale * 0.5f * (image.Width() + image.Height())); - for(const auto& landmarkId : landmarksSet) + const auto& landmarksSetIt = landmarksPerView.find(viewId); + if(landmarksSetIt != landmarksPerView.end()) { - const sfmData::Landmark& landmark = sfmData.getLandmarks().at(landmarkId); - Observation obs = landmark.observations.at(viewId); - // maskLandmarks(obs.x.y(), obs.x.x()) = 1; - for(int y = std::max(obs.x.y() - r, 0.); - y <= std::min(obs.x.y() + r, (double)maskLandmarks.Height() - 1); y++) + const LandmarkIdSet& landmarksSet = landmarksSetIt->second; + for(const auto& landmarkId : landmarksSet) { - for(int x = std::max(obs.x.x() - r, 0.); - x <= std::min(obs.x.x() + r, (double)maskLandmarks.Width() - 1); x++) + const sfmData::Landmark& landmark = sfmData.getLandmarks().at(landmarkId); + Observation obs = landmark.observations.at(viewId); + for(int y = std::max(obs.x.y() - r, 0.); + y <= std::min(obs.x.y() + r, (double)maskLandmarks.Height() - 1); y++) { - maskLandmarks(y, x) = 1; + for(int x = std::max(obs.x.x() - r, 0.); + x <= std::min(obs.x.x() + r, (double)maskLandmarks.Width() - 1); x++) + { + maskLandmarks(y, x) = std::numeric_limits::max(); + } } } } } image::Image mask; + bool maskLoaded = false; if(tryLoadMask(&mask, masksFolders, viewId, srcImage, maskExtension)) - mask = mask * maskLandmarks; - else - mask = maskLandmarks; - process>(dstColorImage, cam, metadata, srcImage, evCorrection, exposureCompensation, [&mask] (Image & image) - { - if(image.Width() * image.Height() != mask.Width() * mask.Height()) + maskLoaded = true; + process>( + dstColorImage, cam, metadata, srcImage, evCorrection, exposureCompensation, + [&maskLoaded, &mask, &maskLandmarks, &doMaskLandmarks](Image& image) { - ALICEVISION_LOG_WARNING("Invalid image mask size: mask is ignored."); - return; - } + if(maskLoaded && (image.Width() * image.Height() != mask.Width() * mask.Height())) + { + ALICEVISION_LOG_WARNING("Invalid image mask size: mask is ignored."); + return; + } - for(int pix = 0; pix < image.Width() * image.Height(); ++pix) - { - const bool masked = (mask(pix) == 0); - image(pix).a() = masked ? 0.f : 1.f; - } - }); + for(int pix = 0; pix < image.Width() * image.Height(); ++pix) + { + image(pix).a() = (maskLoaded && mask(pix) == 0) ? 0.f : (doMaskLandmarks && maskLandmarks(pix) == 127) ? .5f : 1.f; + } + }); /* else { const auto noMaskingFunc = [](Image& image) {}; From 160c7bb7a402f61f57b987135b625ebd179f443d Mon Sep 17 00:00:00 2001 From: Abdelrahman AL MAROUK Date: Mon, 14 Aug 2023 12:58:35 +0200 Subject: [PATCH 07/37] [filterSfM] add option for not using feature scale when computing pixel size --- src/software/pipeline/main_filterSfM.cpp | 15 ++++++++++----- 1 file changed, 10 insertions(+), 5 deletions(-) diff --git a/src/software/pipeline/main_filterSfM.cpp b/src/software/pipeline/main_filterSfM.cpp index 98e39cc469..f5b85f993f 100644 --- a/src/software/pipeline/main_filterSfM.cpp +++ b/src/software/pipeline/main_filterSfM.cpp @@ -120,7 +120,7 @@ class PixSizeSearch inline double worstDist() const { return _radius_sq; } }; -bool filterLandmarks(SfMData& sfmData, double radiusScale) +bool filterLandmarks(SfMData& sfmData, double radiusScale, bool useFeatureScale) { std::vector landmarksData(sfmData.getLandmarks().size()); { @@ -146,8 +146,10 @@ bool filterLandmarks(SfMData& sfmData, double radiusScale) for(const auto& observationPair : landmark.observations) { const IndexT viewId = observationPair.first; - pixSize += mp.getCamPixelSize(Point3d(landmark.X.x(), landmark.X.y(), landmark.X.z()), - mp.getIndexFromViewId(viewId), observationPair.second.scale); + pixSize += mp.getCamPixelSize( + Point3d(landmark.X.x(), landmark.X.y(), landmark.X.z()), + mp.getIndexFromViewId(viewId), + useFeatureScale ? observationPair.second.scale : 1); n++; } pixSize /= n; @@ -269,6 +271,7 @@ int aliceVision_main(int argc, char *argv[]) std::string outputSfmFilename; int maxNbObservationsPerLandmark = 5; double radiusScale = 2; + bool useFeatureScale = true; // user optional parameters std::vector featuresFolders; @@ -287,7 +290,9 @@ int aliceVision_main(int argc, char *argv[]) ("maxNbObservationsPerLandmark", po::value(&maxNbObservationsPerLandmark)->default_value(maxNbObservationsPerLandmark), "Maximum number of allowed observations per landmark.") ("radiusScale", po::value(&radiusScale)->default_value(radiusScale), - "Scale factor applied to pixel size based radius filter applied to landmarks.") + "Scale factor applied to pixel size based radius filter applied to landmarks.")( + "useFeatureScale", po::value(&useFeatureScale)->default_value(useFeatureScale), + "If true, use feature scale for computing pixel size. Otherwise, use a scale of 1 pixel.") ("featuresFolders,f", po::value>(&featuresFolders)->multitoken(), "Path to folder(s) containing the extracted features.") ("matchesFolders,m", po::value>(&matchesFolders)->multitoken(), @@ -322,7 +327,7 @@ int aliceVision_main(int argc, char *argv[]) if(radiusScale > 0) { ALICEVISION_LOG_INFO("Filtering landmarks: started."); - filterLandmarks(sfmData, radiusScale); + filterLandmarks(sfmData, radiusScale, useFeatureScale); ALICEVISION_LOG_INFO("Filtering landmarks: done."); } From 309e771f1fd8f7b386d0367e6880b9611d4ad348 Mon Sep 17 00:00:00 2001 From: Abdelrahman AL MAROUK Date: Wed, 16 Aug 2023 14:57:29 +0200 Subject: [PATCH 08/37] [filterSfM] fix bug - fix reference and constantness problem - fix pre-allocation of vector --- src/software/pipeline/main_filterSfM.cpp | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/src/software/pipeline/main_filterSfM.cpp b/src/software/pipeline/main_filterSfM.cpp index f5b85f993f..1c11b0df04 100644 --- a/src/software/pipeline/main_filterSfM.cpp +++ b/src/software/pipeline/main_filterSfM.cpp @@ -209,19 +209,19 @@ bool filterLandmarks(SfMData& sfmData, double radiusScale, bool useFeatureScale) bool filterObservations(SfMData& sfmData, int maxNbObservationsPerLandmark) { - std::vector landmarksData(sfmData.getLandmarks().size()); + std::vector landmarksData(sfmData.getLandmarks().size()); { size_t i = 0; - for(const auto& it : sfmData.getLandmarks()) + for(auto& it : sfmData.getLandmarks()) { - landmarksData[i++] = it.second; + landmarksData[i++] = &it.second; } } #pragma omp parallel for for(auto i = 0; i < landmarksData.size(); i++) { - sfmData::Landmark landmark = landmarksData[i]; + sfmData::Landmark& landmark = *landmarksData[i]; // check number of observations if(landmark.observations.size() <= maxNbObservationsPerLandmark) @@ -232,7 +232,8 @@ bool filterObservations(SfMData& sfmData, int maxNbObservationsPerLandmark) // compute observation scores - std::vector> observationScores(landmark.observations.size()); + std::vector> observationScores; + observationScores.reserve(landmark.observations.size()); for(const auto& observationPair : landmark.observations) { From 73c6222500cc740e9f14f4e46b923ec0bcfff65d Mon Sep 17 00:00:00 2001 From: Abdelrahman AL MAROUK Date: Wed, 16 Aug 2023 17:33:57 +0200 Subject: [PATCH 09/37] [filterSfM] add new landmark filtering option filter unstable landmarks that have an insufficient number of observations (specified as a parameter) --- src/software/pipeline/main_filterSfM.cpp | 67 ++++++++++++++++++------ 1 file changed, 52 insertions(+), 15 deletions(-) diff --git a/src/software/pipeline/main_filterSfM.cpp b/src/software/pipeline/main_filterSfM.cpp index 1c11b0df04..438189c6f9 100644 --- a/src/software/pipeline/main_filterSfM.cpp +++ b/src/software/pipeline/main_filterSfM.cpp @@ -120,17 +120,50 @@ class PixSizeSearch inline double worstDist() const { return _radius_sq; } }; -bool filterLandmarks(SfMData& sfmData, double radiusScale, bool useFeatureScale) +bool filterLandmarks(SfMData& sfmData, double radiusScale, bool useFeatureScale, int minNbObservationsPerLandmark) { - std::vector landmarksData(sfmData.getLandmarks().size()); + const auto initialNbLandmarks = sfmData.getLandmarks().size(); + std::vector landmarksData(initialNbLandmarks); { size_t i = 0; - for(const auto& it : sfmData.getLandmarks()) + if (minNbObservationsPerLandmark > 0) { - landmarksData[i++] = it.second; + ALICEVISION_LOG_INFO("Removing landmarks having an insufficient number of observations: started."); + for(auto& it : sfmData.getLandmarks()) + { + if(it.second.observations.size() < minNbObservationsPerLandmark) + continue; + landmarksData[i++] = it.second; + } + landmarksData.resize(i); + ALICEVISION_LOG_INFO( + "Removed " << (initialNbLandmarks - i) << + " landmarks out of " << initialNbLandmarks << + ", i.e. " << ((initialNbLandmarks - i) * 100.f / initialNbLandmarks) << + " % " + ); + ALICEVISION_LOG_INFO("Removing landmarks having an insufficient number of observations: done."); + } + else + { + for(auto& it : sfmData.getLandmarks()) + { + landmarksData[i++] = it.second; + } } } + if (radiusScale <= 0) + { + std::vector> filteredLandmarks(landmarksData.size()); + for(IndexT newIdx = 0; newIdx < landmarksData.size(); newIdx++) + { + filteredLandmarks[newIdx] = std::make_pair(newIdx, landmarksData[newIdx]); + } + sfmData.getLandmarks() = Landmarks(filteredLandmarks.begin(), filteredLandmarks.end()); + return true; + } + mvsUtils::MultiViewParams mp(sfmData, "", "", "", false); std::vector landmarksPixSize(landmarksData.size()); @@ -168,7 +201,7 @@ bool filterLandmarks(SfMData& sfmData, double radiusScale, bool useFeatureScale) ALICEVISION_LOG_INFO("KdTree created for " << landmarksData.size() << " points."); std::vector toRemove(landmarksData.size(), false); - ALICEVISION_LOG_INFO("Identifying landmarks to remove: started."); + ALICEVISION_LOG_INFO("Identifying landmarks to remove based on pixel size: started."); size_t nbToRemove = 0; #pragma omp parallel for reduction(+:nbToRemove) @@ -184,14 +217,14 @@ bool filterLandmarks(SfMData& sfmData, double radiusScale, bool useFeatureScale) } ALICEVISION_LOG_INFO( - "Identified " << (landmarksData.size() - nbToRemove) << - " landmarks to remove out of " << (landmarksData.size()) << - ", i.e. " << ((landmarksData.size() - nbToRemove) * 100.f / landmarksData.size()) << + "Identified " << nbToRemove << + " landmarks to remove out of " << initialNbLandmarks << + ", i.e. " << (nbToRemove * 100.f / initialNbLandmarks) << " % " ); ALICEVISION_LOG_INFO("Identifying landmarks to remove: done."); - ALICEVISION_LOG_INFO("Removing landmarks: started."); + ALICEVISION_LOG_INFO("Removing landmarks based on pixel size: started."); std::vector> filteredLandmarks(landmarksData.size() - nbToRemove); IndexT newIdx = 0; for (size_t i = 0; i < landmarksData.size(); i++) @@ -202,7 +235,7 @@ bool filterLandmarks(SfMData& sfmData, double radiusScale, bool useFeatureScale) } } sfmData.getLandmarks() = Landmarks(filteredLandmarks.begin(), filteredLandmarks.end()); - ALICEVISION_LOG_INFO("Removing landmarks: done."); + ALICEVISION_LOG_INFO("Removing landmarks based on pixel size: done."); return true; } @@ -270,7 +303,8 @@ int aliceVision_main(int argc, char *argv[]) // std::string verboseLevel = system::EVerboseLevel_enumToString(system::Logger::getDefaultVerboseLevel()); std::string inputSfmFilename; std::string outputSfmFilename; - int maxNbObservationsPerLandmark = 5; + int maxNbObservationsPerLandmark = 2; + int minNbObservationsPerLandmark = 5; double radiusScale = 2; bool useFeatureScale = true; @@ -290,9 +324,11 @@ int aliceVision_main(int argc, char *argv[]) optionalParams.add_options() ("maxNbObservationsPerLandmark", po::value(&maxNbObservationsPerLandmark)->default_value(maxNbObservationsPerLandmark), "Maximum number of allowed observations per landmark.") + ("minNbObservationsPerLandmark", po::value(&minNbObservationsPerLandmark)->default_value(minNbObservationsPerLandmark), + "Minimum number of observations required to keep a landmark.") ("radiusScale", po::value(&radiusScale)->default_value(radiusScale), - "Scale factor applied to pixel size based radius filter applied to landmarks.")( - "useFeatureScale", po::value(&useFeatureScale)->default_value(useFeatureScale), + "Scale factor applied to pixel size based radius filter applied to landmarks.") + ("useFeatureScale", po::value(&useFeatureScale)->default_value(useFeatureScale), "If true, use feature scale for computing pixel size. Otherwise, use a scale of 1 pixel.") ("featuresFolders,f", po::value>(&featuresFolders)->multitoken(), "Path to folder(s) containing the extracted features.") @@ -325,10 +361,11 @@ int aliceVision_main(int argc, char *argv[]) } // filter SfM data - if(radiusScale > 0) + + if(radiusScale > 0 || minNbObservationsPerLandmark > 0) { ALICEVISION_LOG_INFO("Filtering landmarks: started."); - filterLandmarks(sfmData, radiusScale, useFeatureScale); + filterLandmarks(sfmData, radiusScale, useFeatureScale, minNbObservationsPerLandmark); ALICEVISION_LOG_INFO("Filtering landmarks: done."); } From 36abb243bacba4392b3a3abf86662fa1b43b98ec Mon Sep 17 00:00:00 2001 From: Abdelrahman AL MAROUK Date: Tue, 29 Aug 2023 14:34:40 +0200 Subject: [PATCH 10/37] [filterSfM] small change for nanoflann update --- src/software/pipeline/main_filterSfM.cpp | 2 +- src/software/pipeline/main_prepareDenseScene.cpp | 5 ----- 2 files changed, 1 insertion(+), 6 deletions(-) diff --git a/src/software/pipeline/main_filterSfM.cpp b/src/software/pipeline/main_filterSfM.cpp index 438189c6f9..2afb2b5960 100644 --- a/src/software/pipeline/main_filterSfM.cpp +++ b/src/software/pipeline/main_filterSfM.cpp @@ -208,7 +208,7 @@ bool filterLandmarks(SfMData& sfmData, double radiusScale, bool useFeatureScale, for(auto i = 0; i < landmarksData.size(); i++) { PixSizeSearch search(landmarksPixSize[i] * radiusScale, landmarksPixSize, i); - bool found = tree.findNeighbors(search, landmarksData[i].X.data(), nanoflann::SearchParams()); + bool found = tree.findNeighbors(search, landmarksData[i].X.data(), nanoflann::SearchParameters()); if (found) { toRemove[i] = true; diff --git a/src/software/pipeline/main_prepareDenseScene.cpp b/src/software/pipeline/main_prepareDenseScene.cpp index d3ef386f06..dc4b0307ea 100644 --- a/src/software/pipeline/main_prepareDenseScene.cpp +++ b/src/software/pipeline/main_prepareDenseScene.cpp @@ -303,11 +303,6 @@ bool prepareDenseScene(const SfMData& sfmData, image(pix).a() = (maskLoaded && mask(pix) == 0) ? 0.f : (doMaskLandmarks && maskLandmarks(pix) == 127) ? .5f : 1.f; } }); - /* else - { - const auto noMaskingFunc = [](Image& image) {}; - process>(dstColorImage, cam, metadata, srcImage, evCorrection, exposureCompensation, noMaskingFunc); - }*/ } ++progressDisplay; From a5e8959857a0628ad8787f0e1609f75a4517c1db Mon Sep 17 00:00:00 2001 From: Abdelrahman AL MAROUK Date: Thu, 7 Sep 2023 15:42:06 +0200 Subject: [PATCH 11/37] [prepareDenseScene] add 2D observation filtering for depth map mask If a 2D projected observation in a view doesn't have enough other neighboring observations, it will not be considered for creating the depth map mask --- src/software/pipeline/CMakeLists.txt | 2 + .../pipeline/main_prepareDenseScene.cpp | 128 ++++++++++++++++-- 2 files changed, 118 insertions(+), 12 deletions(-) diff --git a/src/software/pipeline/CMakeLists.txt b/src/software/pipeline/CMakeLists.txt index 040eec66ad..7011a54703 100644 --- a/src/software/pipeline/CMakeLists.txt +++ b/src/software/pipeline/CMakeLists.txt @@ -327,8 +327,10 @@ if(ALICEVISION_BUILD_SFM) aliceVision_sfmData aliceVision_sfmDataIO Boost::program_options + Boost::filesystem Boost::boost Boost::timer + nanoflann ) # Filter SfM data diff --git a/src/software/pipeline/main_prepareDenseScene.cpp b/src/software/pipeline/main_prepareDenseScene.cpp index dc4b0307ea..46a1c97edf 100644 --- a/src/software/pipeline/main_prepareDenseScene.cpp +++ b/src/software/pipeline/main_prepareDenseScene.cpp @@ -14,6 +14,8 @@ #include #include +#include "nanoflann.hpp" + #include #include @@ -41,6 +43,90 @@ using namespace aliceVision::sfmDataIO; namespace po = boost::program_options; namespace fs = std::filesystem; +static const std::size_t MAX_LEAF_ELEMENTS = 64; + +struct ObservationsAdaptator +{ + using Derived = ObservationsAdaptator; //!< In this case the dataset class is myself. + using T = double; + + /// CRTP helper method + inline const Derived& derived() const { return *static_cast(this); } + /// CRTP helper method + inline Derived& derived() { return *static_cast(this); } + + const std::vector _data; + ObservationsAdaptator(const std::vector& data) + : _data(data) + { + } + + // Must return the number of data points + inline size_t kdtree_get_point_count() const { return _data.size(); } + + // Returns the dim'th component of the idx'th point in the class: + inline T kdtree_get_pt(const size_t idx, int dim) const { return _data[idx]->x(dim); } + + // Optional bounding-box computation: return false to default to a standard bbox computation loop. + // Return true if the BBOX was already computed by the class and returned in "bb" so it can be avoided to redo it + // again. Look at bb.size() to find out the expected dimensionality (e.g. 2 or 3 for point clouds) + template + bool kdtree_get_bbox(BBOX& bb) const + { + return false; + } +}; + +using KdTree = nanoflann::KDTreeSingleIndexAdaptor, + ObservationsAdaptator, 2, /* dim */ + size_t>; + +/** + * A result-set class used when performing a radius based search. + */ +class RadiusKnnSearch +{ +public: + const double _radius_sq; + const int _nb_neighbors; + int nb_found = 0; + + inline RadiusKnnSearch(double radius, int k) + : _radius_sq(radius * radius) + , _nb_neighbors(k) + { + } + + inline bool full() const { return nb_found == _nb_neighbors; } + + inline bool addPoint(double dist, IndexT index) + { + if(dist < _radius_sq) + { + nb_found++; + return nb_found < _nb_neighbors; + } + return true; + } + + inline double worstDist() const { return _radius_sq; } +}; + +const auto getObservationsPerViews(const SfMData& sfmData) +{ + stl::flat_map> observationsPerView; + for(const auto& landIt : sfmData.getLandmarks()) + { + for(const auto& obsIt : landIt.second.observations) + { + IndexT viewId = obsIt.first; + auto& landmarksSet = observationsPerView[viewId]; + landmarksSet.push_back(&obsIt.second); + } + } + return observationsPerView; +} + template void process(const std::string& dstColorImage, const IntrinsicBase* cam, @@ -129,9 +215,9 @@ bool prepareDenseScene(const SfMData& sfmData, const double medianCameraExposure = sfmData.getMedianCameraExposureSetting().getExposure(); ALICEVISION_LOG_INFO("Median Camera Exposure: " << medianCameraExposure << ", Median EV: " << std::log2(1.0 / medianCameraExposure)); - const LandmarksPerView& landmarksPerView = getLandmarksPerViews(sfmData); + const auto& observationsPerView = getObservationsPerViews(sfmData); -#pragma omp parallel for num_threads(3) +// #pragma omp parallel for num_threads(3) for (int i = 0; i < viewIds.size(); ++i) { auto itView = viewIds.begin(); @@ -258,19 +344,37 @@ bool prepareDenseScene(const SfMData& sfmData, bool doMaskLandmarks = landmarksMaskScale > 0.f; if(doMaskLandmarks) { - image::Image image; - readImage(srcImage, image, image::EImageColorSpace::LINEAR); // for the T camera, image alpha should be at least 0.4f * 255 (masking) - maskLandmarks = image::Image(image.Width(), image.Height(), true, 127); - int r = (int)(landmarksMaskScale * 0.5f * (image.Width() + image.Height())); - const auto& landmarksSetIt = landmarksPerView.find(viewId); - if(landmarksSetIt != landmarksPerView.end()) + maskLandmarks = image::Image(view->getWidth(), view->getHeight(), true, 127); + int r = (int)(landmarksMaskScale * 0.5f * (view->getWidth() + view->getHeight())); + const auto& observationsIt = observationsPerView.find(viewId); + if(observationsIt != observationsPerView.end()) { - const LandmarkIdSet& landmarksSet = landmarksSetIt->second; - for(const auto& landmarkId : landmarksSet) + const auto& observations = observationsIt->second; + + ALICEVISION_LOG_INFO("Build nanoflann KdTree index."); + ObservationsAdaptator data(observations); + KdTree tree(2, data, nanoflann::KDTreeSingleIndexAdaptorParams(MAX_LEAF_ELEMENTS)); + tree.buildIndex(); + ALICEVISION_LOG_INFO("KdTree created for " << observations.size() << " points."); + + /*index = std::numeric_limits::max(); + sq_dist = std::numeric_limits::max(); + nanoflann::KNNResultSet resultSet(1); + resultSet.init(&index, &sq_dist); + if(!_tree->findNeighbors(resultSet, p.m, nanoflann::SearchParameters())) + { + return false; + } + return true;*/ + + for(const auto& observation : observations) { - const sfmData::Landmark& landmark = sfmData.getLandmarks().at(landmarkId); - Observation obs = landmark.observations.at(viewId); + const auto& obs = *observation; + RadiusKnnSearch search(r, 5); + bool found = tree.findNeighbors(search, obs.x.data(), nanoflann::SearchParameters()); + if(!found) + continue; for(int y = std::max(obs.x.y() - r, 0.); y <= std::min(obs.x.y() + r, (double)maskLandmarks.Height() - 1); y++) { From 6340c92c4289509c02edc6e033329f5b04220f37 Mon Sep 17 00:00:00 2001 From: Abdelrahman AL MAROUK Date: Fri, 8 Sep 2023 14:05:25 +0200 Subject: [PATCH 12/37] [prepareDenseScene] improve 2D observation filtering - estimate the mask radius dynamically for each view based on average distance to neighbors - filter 2D observations based on quantile on distribution of average distance to neighbors --- src/software/pipeline/CMakeLists.txt | 1 - .../pipeline/main_prepareDenseScene.cpp | 73 +++++++++++++------ 2 files changed, 50 insertions(+), 24 deletions(-) diff --git a/src/software/pipeline/CMakeLists.txt b/src/software/pipeline/CMakeLists.txt index 7011a54703..4ff9f2ab39 100644 --- a/src/software/pipeline/CMakeLists.txt +++ b/src/software/pipeline/CMakeLists.txt @@ -327,7 +327,6 @@ if(ALICEVISION_BUILD_SFM) aliceVision_sfmData aliceVision_sfmDataIO Boost::program_options - Boost::filesystem Boost::boost Boost::timer nanoflann diff --git a/src/software/pipeline/main_prepareDenseScene.cpp b/src/software/pipeline/main_prepareDenseScene.cpp index 46a1c97edf..2e4c1e65ce 100644 --- a/src/software/pipeline/main_prepareDenseScene.cpp +++ b/src/software/pipeline/main_prepareDenseScene.cpp @@ -17,6 +17,9 @@ #include "nanoflann.hpp" #include +#include +#include +#include #include #include @@ -42,6 +45,7 @@ using namespace aliceVision::sfmDataIO; namespace po = boost::program_options; namespace fs = std::filesystem; +using namespace boost::accumulators; static const std::size_t MAX_LEAF_ELEMENTS = 64; @@ -179,7 +183,9 @@ bool prepareDenseScene(const SfMData& sfmData, bool saveMetadata, bool saveMatricesFiles, bool evCorrection, - float landmarksMaskScale) + float landmarksMaskScale, + int nbNeighborObservations, + float percentile) { // defined view Ids std::set viewIds; @@ -217,7 +223,7 @@ bool prepareDenseScene(const SfMData& sfmData, const auto& observationsPerView = getObservationsPerViews(sfmData); -// #pragma omp parallel for num_threads(3) +#pragma omp parallel for for (int i = 0; i < viewIds.size(); ++i) { auto itView = viewIds.begin(); @@ -341,7 +347,7 @@ bool prepareDenseScene(const SfMData& sfmData, } image::Image maskLandmarks; - bool doMaskLandmarks = landmarksMaskScale > 0.f; + bool doMaskLandmarks = nbNeighborObservations > 0 || percentile < 1.f; if(doMaskLandmarks) { // for the T camera, image alpha should be at least 0.4f * 255 (masking) @@ -358,32 +364,45 @@ bool prepareDenseScene(const SfMData& sfmData, tree.buildIndex(); ALICEVISION_LOG_INFO("KdTree created for " << observations.size() << " points."); - /*index = std::numeric_limits::max(); - sq_dist = std::numeric_limits::max(); - nanoflann::KNNResultSet resultSet(1); - resultSet.init(&index, &sq_dist); - if(!_tree->findNeighbors(resultSet, p.m, nanoflann::SearchParameters())) + int n = std::min(nbNeighborObservations, static_cast(observations.size())); + std::vector means(observations.size()); + const std::size_t cacheSize = 1000; + accumulator_set, tag::mean>> acc( + tag::tail::cache_size = cacheSize); + int j = 0; + for (const auto& observation : observations) { - return false; + const auto& obs = *observation; + std::vector indices_(n); + std::vector distances_(n); + tree.knnSearch(obs.x.data(), n, &indices_[0], &distances_[0]); + + std::transform(distances_.begin(), distances_.end(), distances_.begin(), + static_cast(std::sqrt)); + const auto& mean = std::accumulate(distances_.begin(), distances_.end(), 0.0) / n; + means[j++] = mean; + acc(mean); } - return true;*/ + double mean_max = std::numeric_limits::max(); + if(percentile != 1.f) + mean_max = quantile(acc, quantile_probability = percentile); + r = mean(acc); + j = 0; for(const auto& observation : observations) { const auto& obs = *observation; - RadiusKnnSearch search(r, 5); - bool found = tree.findNeighbors(search, obs.x.data(), nanoflann::SearchParameters()); - if(!found) - continue; - for(int y = std::max(obs.x.y() - r, 0.); - y <= std::min(obs.x.y() + r, (double)maskLandmarks.Height() - 1); y++) - { - for(int x = std::max(obs.x.x() - r, 0.); - x <= std::min(obs.x.x() + r, (double)maskLandmarks.Width() - 1); x++) + if(means[j] < mean_max) + for(int y = std::max(obs.x.y() - r, 0.); + y <= std::min(obs.x.y() + r, (double)maskLandmarks.Height() - 1); y++) { - maskLandmarks(y, x) = std::numeric_limits::max(); + for(int x = std::max(obs.x.x() - r, 0.); + x <= std::min(obs.x.x() + r, (double)maskLandmarks.Width() - 1); x++) + { + maskLandmarks(y, x) = std::numeric_limits::max(); + } } - } + j++; } } } @@ -432,6 +451,8 @@ int aliceVision_main(int argc, char* argv[]) bool saveMatricesTxtFiles = false; bool evCorrection = false; float landmarksMaskScale = 0.f; + int nbNeighborObservations = 5; + float percentile = 0.95; // clang-format off po::options_description requiredParams("Required parameters"); @@ -465,7 +486,11 @@ int aliceVision_main(int argc, char* argv[]) "Correct exposure value.") ("landmarksMaskScale", po::value(&landmarksMaskScale)->default_value(landmarksMaskScale), "Scale (relative to image size) of the projection of landmarks to mask images for depth computation.\n" - "If 0, masking using landmarks will not be used."); + "If 0, masking using landmarks will not be used.") + ("nbNeighborObservations", po::value(&nbNeighborObservations)->default_value(nbNeighborObservations), + "Number of neighbor observations to be considered for the landmarks-based masking.") + ("percentile", po::value(&percentile)->default_value(percentile), + "TODO."); // clang-format on CmdLine cmdline("AliceVision prepareDenseScene"); @@ -530,7 +555,9 @@ int aliceVision_main(int argc, char* argv[]) saveMetadata, saveMatricesTxtFiles, evCorrection, - landmarksMaskScale)) + landmarksMaskScale, + nbNeighborObservations, + percentile)) return EXIT_SUCCESS; return EXIT_FAILURE; From c3865a794f2d68e2a7d5c77e5daf61a3123cef2a Mon Sep 17 00:00:00 2001 From: Abdelrahman AL MAROUK Date: Tue, 12 Sep 2023 15:39:18 +0200 Subject: [PATCH 13/37] [filterSfM] propagate neighbors scores when filtering observations --- src/software/pipeline/main_filterSfM.cpp | 181 ++++++++++++++++++++++- 1 file changed, 177 insertions(+), 4 deletions(-) diff --git a/src/software/pipeline/main_filterSfM.cpp b/src/software/pipeline/main_filterSfM.cpp index 2afb2b5960..58cc4a8774 100644 --- a/src/software/pipeline/main_filterSfM.cpp +++ b/src/software/pipeline/main_filterSfM.cpp @@ -84,7 +84,7 @@ using KdTree = nanoflann::KDTreeSingleIndexAdaptor< nanoflann::L2_Simple_Adaptor, LandmarksAdaptator, 3, /* dim */ - IndexT + size_t >; /** @@ -107,7 +107,7 @@ class PixSizeSearch inline bool full() const { return found; } - inline bool addPoint(double dist, IndexT index) + inline bool addPoint(double dist, size_t index) { if(dist < _radius_sq && _pixSize[index] < _pixSize_i) { @@ -291,11 +291,184 @@ bool filterObservations(SfMData& sfmData, int maxNbObservationsPerLandmark) { filteredObservations[observationScorePair.second] = landmark.observations[observationScorePair.second]; } - landmark.observations = filteredObservations; + landmark.observations = std::move(filteredObservations); } return true; } +bool filterObservations2(SfMData& sfmData, int maxNbObservationsPerLandmark, int nbNeighbors = 10, int nbIterations = 5, + double fraction = 0.5) +{ + std::vector landmarksData(sfmData.getLandmarks().size()); + { + size_t i = 0; + for(auto& it : sfmData.getLandmarks()) + { + landmarksData[i++] = it.second; + } + } + + std::vector, std::vector>> viewScoresData(landmarksData.size()); + + ALICEVISION_LOG_INFO("Computing initial observation scores based on distance to observing view: started"); +#pragma omp parallel for + for(auto i = 0; i < landmarksData.size(); i++) + { + const sfmData::Landmark& landmark = landmarksData[i]; + + // compute observation scores + + const auto& nbObservations = landmark.observations.size(); + auto& [viewIds, viewScores] = viewScoresData[i]; + viewIds.reserve(nbObservations); + viewScores.reserve(nbObservations); + double total = 0.; + for(const auto& observationPair : landmark.observations) + { + const IndexT viewId = observationPair.first; + const sfmData::View& view = *(sfmData.getViews().at(viewId)); + const geometry::Pose3 pose = sfmData.getPose(view).getTransform(); + + viewIds.push_back(viewId); + const auto& v = 1. / (pose.center() - landmark.X).squaredNorm(); + total += v; + viewScores.push_back(v); + } + + // normalize view scores + for(auto j = 0; j < nbObservations; j++) + { + viewScores[j] /= total; + } + + // sort by ascending view id order + // for consequent faster access + std::vector idx(nbObservations); + std::iota(idx.begin(), idx.end(), 0); + std::stable_sort(idx.begin(), idx.end(), [&v = viewIds](size_t i1, size_t i2) { return v[i1] < v[i2]; }); + auto ids_temp = viewIds; + auto scores_temp = viewScores; + for(auto j = 0; j < nbObservations; j++) + { + viewIds[j] = ids_temp[idx[j]]; + viewScores[j] = scores_temp[idx[j]]; + } + } + ALICEVISION_LOG_INFO("Computing initial observation scores based on distance to observing view: done"); + + ALICEVISION_LOG_INFO("Computing landmark neighbors and distance-based weights: started"); + ALICEVISION_LOG_INFO("Build nanoflann KdTree index."); + LandmarksAdaptator dataAdaptor(landmarksData); + KdTree tree(3, dataAdaptor, nanoflann::KDTreeSingleIndexAdaptorParams(MAX_LEAF_ELEMENTS)); + tree.buildIndex(); + ALICEVISION_LOG_INFO("KdTree created for " << landmarksData.size() << " points."); + std::vector, std::vector>> neighborsData(landmarksData.size()); +#pragma omp parallel for + for(auto i = 0; i < landmarksData.size(); i++) + { + const sfmData::Landmark& landmark = landmarksData[i]; + auto& [indices_, weights_] = neighborsData[i]; + indices_.resize(nbNeighbors); + weights_.resize(nbNeighbors); + tree.knnSearch(landmark.X.data(), nbNeighbors, &indices_[0], &weights_[0]); + double total = 0.; + for(auto& w : weights_) + { + w = 1. / std::sqrt(w); + total += w; + } + for(auto& w : weights_) + { + w /= total; + } + } + ALICEVISION_LOG_INFO("Computing landmark neighbors and distance-based weights: done"); + + ALICEVISION_LOG_INFO("Propagating neighbors observation scores: started"); + std::vector> viewScoresData_t(landmarksData.size()); + for(auto i = 0; i < nbIterations; i++) + { + ALICEVISION_LOG_INFO("Iteration " << i << "..."); +#pragma omp parallel for + for(auto id = 0; id < landmarksData.size(); id++) + { + const auto& [viewIds, viewScores] = viewScoresData[id]; + auto& viewScores_acc = viewScoresData_t[id]; + viewScores_acc.assign(viewScores.size(), 0.); + double viewScores_total = 0.; + auto& [indices_, weights_] = neighborsData[id]; + for(auto j = 0; j < nbNeighbors; j++) + { + const auto& neighborId = indices_[j]; + const auto& neighborWeight = weights_[j]; + const auto& [viewIds_neighbor, viewScores_neighbor] = viewScoresData[neighborId]; + auto viewIds_it = viewIds.begin(); + auto viewIds_neighbor_it = viewIds_neighbor.begin(); + auto viewScores_neighbor_it = viewScores_neighbor.begin(); + auto viewScores_acc_it = viewScores_acc.begin(); + while(viewIds_it != viewIds.end() && viewIds_neighbor_it != viewIds_neighbor.end()) + { + if(*viewIds_it < *viewIds_neighbor_it) + { + ++viewIds_it; + ++viewScores_acc_it; + } + else + { + if(!(*viewIds_neighbor_it < *viewIds_it)) + { + const auto& v = *viewScores_neighbor_it * neighborWeight; + (*viewScores_acc_it) += v; + viewScores_total += v; + ++viewIds_it; + ++viewScores_acc_it; + } + ++viewIds_neighbor_it; + ++viewScores_neighbor_it; + } + } + } + for(auto j = 0; j < viewScores_acc.size(); j++) + { + viewScores_acc[j] *= fraction / viewScores_total; + viewScores_acc[j] += (1 - fraction) * viewScores[j]; + } + } +#pragma omp parallel for + for(auto id = 0; id < landmarksData.size(); id++) + { + viewScoresData[id].second = std::move(viewScoresData_t[id]); + } + } + ALICEVISION_LOG_INFO("Propagating neighbors observation scores: done"); + + ALICEVISION_LOG_INFO("Selecting observations with best scores: started"); +#pragma omp parallel for + for(auto i = 0; i < landmarksData.size(); i++) + { + sfmData::Landmark& landmark = landmarksData[i]; + const auto& nbObservations = landmark.observations.size(); + auto& [viewIds, viewScores] = viewScoresData[i]; + + // sort by descending view score order + std::vector idx(nbObservations); + std::iota(idx.begin(), idx.end(), 0); + std::stable_sort(idx.begin(), idx.end(), [&v = viewScores](size_t i1, size_t i2) { return v[i1] > v[i2]; }); + + // replace the observations + Observations filteredObservations; + for(auto j = 0; j < maxNbObservationsPerLandmark; j++) + { + const auto& viewId = viewIds[idx[j]]; + filteredObservations[viewId] = landmark.observations[viewId]; + } + landmark.observations = std::move(filteredObservations); + } + ALICEVISION_LOG_INFO("Selecting observations with best scores: done"); + + return true; +} + int aliceVision_main(int argc, char *argv[]) { // command-line parameters @@ -372,7 +545,7 @@ int aliceVision_main(int argc, char *argv[]) if(maxNbObservationsPerLandmark > 0) { ALICEVISION_LOG_INFO("Filtering observations: started."); - filterObservations(sfmData, maxNbObservationsPerLandmark); + filterObservations2(sfmData, maxNbObservationsPerLandmark); ALICEVISION_LOG_INFO("Filtering observations: done."); } From 583342c52e41b9913e65bd69db27f546fdb332f6 Mon Sep 17 00:00:00 2001 From: Abdelrahman AL MAROUK Date: Tue, 12 Sep 2023 17:07:22 +0200 Subject: [PATCH 14/37] [filterSfM] fix KdTree search a point is neighbor to itself if it participated in creating the tree --- src/software/pipeline/main_filterSfM.cpp | 9 ++++++--- src/software/pipeline/main_prepareDenseScene.cpp | 5 +++-- 2 files changed, 9 insertions(+), 5 deletions(-) diff --git a/src/software/pipeline/main_filterSfM.cpp b/src/software/pipeline/main_filterSfM.cpp index 58cc4a8774..5fe1fe5f80 100644 --- a/src/software/pipeline/main_filterSfM.cpp +++ b/src/software/pipeline/main_filterSfM.cpp @@ -368,9 +368,12 @@ bool filterObservations2(SfMData& sfmData, int maxNbObservationsPerLandmark, int { const sfmData::Landmark& landmark = landmarksData[i]; auto& [indices_, weights_] = neighborsData[i]; - indices_.resize(nbNeighbors); - weights_.resize(nbNeighbors); - tree.knnSearch(landmark.X.data(), nbNeighbors, &indices_[0], &weights_[0]); + // a landmark is a neighbor to itself with zero distance + indices_.resize(nbNeighbors + 1); + weights_.resize(nbNeighbors + 1); + tree.knnSearch(landmark.X.data(), nbNeighbors + 1, &indices_[0], &weights_[0]); + indices_.erase(indices_.begin()); + weights_.erase(weights_.begin()); double total = 0.; for(auto& w : weights_) { diff --git a/src/software/pipeline/main_prepareDenseScene.cpp b/src/software/pipeline/main_prepareDenseScene.cpp index 2e4c1e65ce..c32922625a 100644 --- a/src/software/pipeline/main_prepareDenseScene.cpp +++ b/src/software/pipeline/main_prepareDenseScene.cpp @@ -364,7 +364,8 @@ bool prepareDenseScene(const SfMData& sfmData, tree.buildIndex(); ALICEVISION_LOG_INFO("KdTree created for " << observations.size() << " points."); - int n = std::min(nbNeighborObservations, static_cast(observations.size())); + int n = std::min(nbNeighborObservations, static_cast(observations.size() - 1)) + 1; + // note that the observation is a neighbor to itself with zero distance, hence the +/- 1 std::vector means(observations.size()); const std::size_t cacheSize = 1000; accumulator_set, tag::mean>> acc( @@ -379,7 +380,7 @@ bool prepareDenseScene(const SfMData& sfmData, std::transform(distances_.begin(), distances_.end(), distances_.begin(), static_cast(std::sqrt)); - const auto& mean = std::accumulate(distances_.begin(), distances_.end(), 0.0) / n; + const auto& mean = std::accumulate(distances_.begin(), distances_.end(), 0.0) / (n - 1); means[j++] = mean; acc(mean); } From a8c6086d85aa85612619761a2b0bdb8277e036c6 Mon Sep 17 00:00:00 2001 From: Abdelrahman AL MAROUK Date: Tue, 12 Sep 2023 18:04:41 +0200 Subject: [PATCH 15/37] [filterSfM] fix references --- src/software/pipeline/main_filterSfM.cpp | 23 ++++++++++++++++------- 1 file changed, 16 insertions(+), 7 deletions(-) diff --git a/src/software/pipeline/main_filterSfM.cpp b/src/software/pipeline/main_filterSfM.cpp index 5fe1fe5f80..82db8fb0a4 100644 --- a/src/software/pipeline/main_filterSfM.cpp +++ b/src/software/pipeline/main_filterSfM.cpp @@ -56,9 +56,18 @@ struct LandmarksAdaptator inline Derived& derived() { return *static_cast(this); } const std::vector _data; + const std::vector _data_ptr; + bool usePtr; LandmarksAdaptator(const std::vector& data) : _data(data) { + usePtr = false; + } + + LandmarksAdaptator(const std::vector& data) + : _data_ptr(data) + { + usePtr = true; } // Must return the number of data points @@ -67,7 +76,7 @@ struct LandmarksAdaptator // Returns the dim'th component of the idx'th point in the class: inline T kdtree_get_pt(const size_t idx, int dim) const { - return _data[idx].X(dim); + return usePtr ? _data_ptr[idx]->X(dim) : _data[idx].X(dim); } // Optional bounding-box computation: return false to default to a standard bbox computation loop. @@ -234,7 +243,7 @@ bool filterLandmarks(SfMData& sfmData, double radiusScale, bool useFeatureScale, filteredLandmarks[newIdx++] = std::make_pair(newIdx, landmarksData[i]); } } - sfmData.getLandmarks() = Landmarks(filteredLandmarks.begin(), filteredLandmarks.end()); + sfmData.getLandmarks() = std::move(Landmarks(filteredLandmarks.begin(), filteredLandmarks.end())); ALICEVISION_LOG_INFO("Removing landmarks based on pixel size: done."); return true; @@ -299,12 +308,12 @@ bool filterObservations(SfMData& sfmData, int maxNbObservationsPerLandmark) bool filterObservations2(SfMData& sfmData, int maxNbObservationsPerLandmark, int nbNeighbors = 10, int nbIterations = 5, double fraction = 0.5) { - std::vector landmarksData(sfmData.getLandmarks().size()); + std::vector landmarksData(sfmData.getLandmarks().size()); { size_t i = 0; for(auto& it : sfmData.getLandmarks()) { - landmarksData[i++] = it.second; + landmarksData[i++] = &it.second; } } @@ -314,7 +323,7 @@ bool filterObservations2(SfMData& sfmData, int maxNbObservationsPerLandmark, int #pragma omp parallel for for(auto i = 0; i < landmarksData.size(); i++) { - const sfmData::Landmark& landmark = landmarksData[i]; + const sfmData::Landmark& landmark = *landmarksData[i]; // compute observation scores @@ -366,7 +375,7 @@ bool filterObservations2(SfMData& sfmData, int maxNbObservationsPerLandmark, int #pragma omp parallel for for(auto i = 0; i < landmarksData.size(); i++) { - const sfmData::Landmark& landmark = landmarksData[i]; + const sfmData::Landmark& landmark = *landmarksData[i]; auto& [indices_, weights_] = neighborsData[i]; // a landmark is a neighbor to itself with zero distance indices_.resize(nbNeighbors + 1); @@ -449,7 +458,7 @@ bool filterObservations2(SfMData& sfmData, int maxNbObservationsPerLandmark, int #pragma omp parallel for for(auto i = 0; i < landmarksData.size(); i++) { - sfmData::Landmark& landmark = landmarksData[i]; + sfmData::Landmark& landmark = *landmarksData[i]; const auto& nbObservations = landmark.observations.size(); auto& [viewIds, viewScores] = viewScoresData[i]; From 278b5926d783e4f11c1af619182e667519f0e7a0 Mon Sep 17 00:00:00 2001 From: Abdelrahman AL MAROUK Date: Wed, 13 Sep 2023 09:29:00 +0200 Subject: [PATCH 16/37] [filterSfM] add params + remove deprecated function --- src/software/pipeline/main_filterSfM.cpp | 89 ++++++------------------ 1 file changed, 22 insertions(+), 67 deletions(-) diff --git a/src/software/pipeline/main_filterSfM.cpp b/src/software/pipeline/main_filterSfM.cpp index 82db8fb0a4..dd0e9bb819 100644 --- a/src/software/pipeline/main_filterSfM.cpp +++ b/src/software/pipeline/main_filterSfM.cpp @@ -249,64 +249,8 @@ bool filterLandmarks(SfMData& sfmData, double radiusScale, bool useFeatureScale, return true; } -bool filterObservations(SfMData& sfmData, int maxNbObservationsPerLandmark) -{ - std::vector landmarksData(sfmData.getLandmarks().size()); - { - size_t i = 0; - for(auto& it : sfmData.getLandmarks()) - { - landmarksData[i++] = &it.second; - } - } - - #pragma omp parallel for - for(auto i = 0; i < landmarksData.size(); i++) - { - sfmData::Landmark& landmark = *landmarksData[i]; - - // check number of observations - if(landmark.observations.size() <= maxNbObservationsPerLandmark) - continue; - - // // check angle between observations - // if(!checkLandmarkMinObservationAngle(sfmData, landmark, minObservationAngle)) - - // compute observation scores - - std::vector> observationScores; - observationScores.reserve(landmark.observations.size()); - - for(const auto& observationPair : landmark.observations) - { - const IndexT viewId = observationPair.first; - const sfmData::View& view = *(sfmData.getViews().at(viewId)); - const geometry::Pose3 pose = sfmData.getPose(view).getTransform(); - - observationScores.push_back(std::pair( - (pose.center() - landmark.X).squaredNorm(), - viewId - )); - } - - // sort observations by ascending score order - std::stable_sort(observationScores.begin(), observationScores.end()); - // take only best observations - observationScores.resize(maxNbObservationsPerLandmark); - - // replace the observations - Observations filteredObservations; - for(const auto& observationScorePair : observationScores) - { - filteredObservations[observationScorePair.second] = landmark.observations[observationScorePair.second]; - } - landmark.observations = std::move(filteredObservations); - } - return true; -} - -bool filterObservations2(SfMData& sfmData, int maxNbObservationsPerLandmark, int nbNeighbors = 10, int nbIterations = 5, - double fraction = 0.5) +bool filterObservations(SfMData& sfmData, int maxNbObservationsPerLandmark, int nbNeighbors, double neighborsInfluence, + int nbIterations) { std::vector landmarksData(sfmData.getLandmarks().size()); { @@ -371,16 +315,18 @@ bool filterObservations2(SfMData& sfmData, int maxNbObservationsPerLandmark, int KdTree tree(3, dataAdaptor, nanoflann::KDTreeSingleIndexAdaptorParams(MAX_LEAF_ELEMENTS)); tree.buildIndex(); ALICEVISION_LOG_INFO("KdTree created for " << landmarksData.size() << " points."); + // note that the landmark is a neighbor to itself with zero distance, hence the +/- 1 + int nbNeighbors_ = std::min(nbNeighbors, static_cast(landmarksData.size() - 1)) + 1; std::vector, std::vector>> neighborsData(landmarksData.size()); #pragma omp parallel for for(auto i = 0; i < landmarksData.size(); i++) { const sfmData::Landmark& landmark = *landmarksData[i]; auto& [indices_, weights_] = neighborsData[i]; - // a landmark is a neighbor to itself with zero distance - indices_.resize(nbNeighbors + 1); - weights_.resize(nbNeighbors + 1); - tree.knnSearch(landmark.X.data(), nbNeighbors + 1, &indices_[0], &weights_[0]); + indices_.resize(nbNeighbors_); + weights_.resize(nbNeighbors_); + tree.knnSearch(landmark.X.data(), nbNeighbors_, &indices_[0], &weights_[0]); + // a landmark is a neighbor to itself with zero distance, remove it indices_.erase(indices_.begin()); weights_.erase(weights_.begin()); double total = 0.; @@ -442,8 +388,8 @@ bool filterObservations2(SfMData& sfmData, int maxNbObservationsPerLandmark, int } for(auto j = 0; j < viewScores_acc.size(); j++) { - viewScores_acc[j] *= fraction / viewScores_total; - viewScores_acc[j] += (1 - fraction) * viewScores[j]; + viewScores_acc[j] *= neighborsInfluence / viewScores_total; + viewScores_acc[j] += (1 - neighborsInfluence) * viewScores[j]; } } #pragma omp parallel for @@ -485,13 +431,15 @@ int aliceVision_main(int argc, char *argv[]) { // command-line parameters - // std::string verboseLevel = system::EVerboseLevel_enumToString(system::Logger::getDefaultVerboseLevel()); std::string inputSfmFilename; std::string outputSfmFilename; int maxNbObservationsPerLandmark = 2; int minNbObservationsPerLandmark = 5; double radiusScale = 2; bool useFeatureScale = true; + int nbNeighbors = 10; + double neighborsInfluence = 0.5; + int nbIterations = 5; // user optional parameters std::vector featuresFolders; @@ -515,6 +463,13 @@ int aliceVision_main(int argc, char *argv[]) "Scale factor applied to pixel size based radius filter applied to landmarks.") ("useFeatureScale", po::value(&useFeatureScale)->default_value(useFeatureScale), "If true, use feature scale for computing pixel size. Otherwise, use a scale of 1 pixel.") + ("nbNeighbors", po::value(&nbNeighbors)->default_value(nbNeighbors), + "Number of neighbor landmarks used in making the decision for best observations.") + ("neighborsInfluence", po::value(&neighborsInfluence)->default_value(neighborsInfluence), + "Specifies how much influential the neighbors are in selecting the best observations." + "Between 0. and 1., the closer to 1., the more influencial the neighborhood is.") + ("nbIterations", po::value(&nbIterations)->default_value(nbIterations), + "Number of iterations to propagate neighbors information.") ("featuresFolders,f", po::value>(&featuresFolders)->multitoken(), "Path to folder(s) containing the extracted features.") ("matchesFolders,m", po::value>(&matchesFolders)->multitoken(), @@ -522,7 +477,7 @@ int aliceVision_main(int argc, char *argv[]) ("describerTypes,d", po::value(&describerTypesName)->default_value(describerTypesName), feature::EImageDescriberType_informations().c_str()); - CmdLine cmdline("AliceVision SfM filtering."); + CmdLine cmdline("AliceVision SfM filtering."); // TODO add description cmdline.add(requiredParams); cmdline.add(optionalParams); if (!cmdline.execute(argc, argv)) @@ -557,7 +512,7 @@ int aliceVision_main(int argc, char *argv[]) if(maxNbObservationsPerLandmark > 0) { ALICEVISION_LOG_INFO("Filtering observations: started."); - filterObservations2(sfmData, maxNbObservationsPerLandmark); + filterObservations(sfmData, maxNbObservationsPerLandmark, nbNeighbors, neighborsInfluence, nbIterations); ALICEVISION_LOG_INFO("Filtering observations: done."); } From 44cee478cfa3d8d38a203ae60f91b302d12421af Mon Sep 17 00:00:00 2001 From: Abdelrahman AL MAROUK Date: Wed, 13 Sep 2023 11:13:35 +0200 Subject: [PATCH 17/37] [filterSfM] add check for number of observations --- src/software/pipeline/main_filterSfM.cpp | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/src/software/pipeline/main_filterSfM.cpp b/src/software/pipeline/main_filterSfM.cpp index dd0e9bb819..1d2771fae5 100644 --- a/src/software/pipeline/main_filterSfM.cpp +++ b/src/software/pipeline/main_filterSfM.cpp @@ -408,6 +408,10 @@ bool filterObservations(SfMData& sfmData, int maxNbObservationsPerLandmark, int const auto& nbObservations = landmark.observations.size(); auto& [viewIds, viewScores] = viewScoresData[i]; + // check number of observations + if(landmark.observations.size() <= maxNbObservationsPerLandmark) + continue; + // sort by descending view score order std::vector idx(nbObservations); std::iota(idx.begin(), idx.end(), 0); From 52ff5ccf77e78df38e6ffa5bbc2c627581051dd8 Mon Sep 17 00:00:00 2001 From: Abdelrahman AL MAROUK Date: Wed, 13 Sep 2023 11:13:49 +0200 Subject: [PATCH 18/37] [filterSfM] add comments --- src/software/pipeline/main_filterSfM.cpp | 22 ++++++++++++++++++++-- 1 file changed, 20 insertions(+), 2 deletions(-) diff --git a/src/software/pipeline/main_filterSfM.cpp b/src/software/pipeline/main_filterSfM.cpp index 1d2771fae5..3c139f2fff 100644 --- a/src/software/pipeline/main_filterSfM.cpp +++ b/src/software/pipeline/main_filterSfM.cpp @@ -252,6 +252,7 @@ bool filterLandmarks(SfMData& sfmData, double radiusScale, bool useFeatureScale, bool filterObservations(SfMData& sfmData, int maxNbObservationsPerLandmark, int nbNeighbors, double neighborsInfluence, int nbIterations) { + // store in vector for faster access std::vector landmarksData(sfmData.getLandmarks().size()); { size_t i = 0; @@ -261,6 +262,7 @@ bool filterObservations(SfMData& sfmData, int maxNbObservationsPerLandmark, int } } + // contains the observing view ids for each landmark with their corresponding scores std::vector, std::vector>> viewScoresData(landmarksData.size()); ALICEVISION_LOG_INFO("Computing initial observation scores based on distance to observing view: started"); @@ -270,11 +272,11 @@ bool filterObservations(SfMData& sfmData, int maxNbObservationsPerLandmark, int const sfmData::Landmark& landmark = *landmarksData[i]; // compute observation scores - const auto& nbObservations = landmark.observations.size(); auto& [viewIds, viewScores] = viewScoresData[i]; viewIds.reserve(nbObservations); viewScores.reserve(nbObservations); + // accumulator for normalizing the scores double total = 0.; for(const auto& observationPair : landmark.observations) { @@ -283,6 +285,7 @@ bool filterObservations(SfMData& sfmData, int maxNbObservationsPerLandmark, int const geometry::Pose3 pose = sfmData.getPose(view).getTransform(); viewIds.push_back(viewId); + // score is the inverse of distance to observations const auto& v = 1. / (pose.center() - landmark.X).squaredNorm(); total += v; viewScores.push_back(v); @@ -296,9 +299,12 @@ bool filterObservations(SfMData& sfmData, int maxNbObservationsPerLandmark, int // sort by ascending view id order // for consequent faster access + + // indices that sort the view ids std::vector idx(nbObservations); std::iota(idx.begin(), idx.end(), 0); std::stable_sort(idx.begin(), idx.end(), [&v = viewIds](size_t i1, size_t i2) { return v[i1] < v[i2]; }); + // apply sorting to both view ids and view scores for correspondance auto ids_temp = viewIds; auto scores_temp = viewScores; for(auto j = 0; j < nbObservations; j++) @@ -317,6 +323,7 @@ bool filterObservations(SfMData& sfmData, int maxNbObservationsPerLandmark, int ALICEVISION_LOG_INFO("KdTree created for " << landmarksData.size() << " points."); // note that the landmark is a neighbor to itself with zero distance, hence the +/- 1 int nbNeighbors_ = std::min(nbNeighbors, static_cast(landmarksData.size() - 1)) + 1; + // contains the neighbor landmarks ids with their corresponding weights std::vector, std::vector>> neighborsData(landmarksData.size()); #pragma omp parallel for for(auto i = 0; i < landmarksData.size(); i++) @@ -329,12 +336,15 @@ bool filterObservations(SfMData& sfmData, int maxNbObservationsPerLandmark, int // a landmark is a neighbor to itself with zero distance, remove it indices_.erase(indices_.begin()); weights_.erase(weights_.begin()); + // accumulator used for normalisation double total = 0.; for(auto& w : weights_) { + // weight is the inverse of distance between a landmark and its neighbor w = 1. / std::sqrt(w); total += w; } + // normalize weights for(auto& w : weights_) { w /= total; @@ -343,6 +353,7 @@ bool filterObservations(SfMData& sfmData, int maxNbObservationsPerLandmark, int ALICEVISION_LOG_INFO("Computing landmark neighbors and distance-based weights: done"); ALICEVISION_LOG_INFO("Propagating neighbors observation scores: started"); + // new view scores at iteration t std::vector> viewScoresData_t(landmarksData.size()); for(auto i = 0; i < nbIterations; i++) { @@ -352,7 +363,9 @@ bool filterObservations(SfMData& sfmData, int maxNbObservationsPerLandmark, int { const auto& [viewIds, viewScores] = viewScoresData[id]; auto& viewScores_acc = viewScoresData_t[id]; + // initialize to zero, will first contain the weighted average scores of neighbors viewScores_acc.assign(viewScores.size(), 0.); + // accumulator for normalisation double viewScores_total = 0.; auto& [indices_, weights_] = neighborsData[id]; for(auto j = 0; j < nbNeighbors; j++) @@ -360,6 +373,7 @@ bool filterObservations(SfMData& sfmData, int maxNbObservationsPerLandmark, int const auto& neighborId = indices_[j]; const auto& neighborWeight = weights_[j]; const auto& [viewIds_neighbor, viewScores_neighbor] = viewScoresData[neighborId]; + // loop over common views auto viewIds_it = viewIds.begin(); auto viewIds_neighbor_it = viewIds_neighbor.begin(); auto viewScores_neighbor_it = viewScores_neighbor.begin(); @@ -373,6 +387,7 @@ bool filterObservations(SfMData& sfmData, int maxNbObservationsPerLandmark, int } else { + // if same view, accumulate weighted scores if(!(*viewIds_neighbor_it < *viewIds_it)) { const auto& v = *viewScores_neighbor_it * neighborWeight; @@ -388,10 +403,13 @@ bool filterObservations(SfMData& sfmData, int maxNbObservationsPerLandmark, int } for(auto j = 0; j < viewScores_acc.size(); j++) { + // normalize score and apply influence factor viewScores_acc[j] *= neighborsInfluence / viewScores_total; + // combine weghted neighbor scores and the landmark's own scores viewScores_acc[j] += (1 - neighborsInfluence) * viewScores[j]; } } + // update scores at end of iteration #pragma omp parallel for for(auto id = 0; id < landmarksData.size(); id++) { @@ -417,7 +435,7 @@ bool filterObservations(SfMData& sfmData, int maxNbObservationsPerLandmark, int std::iota(idx.begin(), idx.end(), 0); std::stable_sort(idx.begin(), idx.end(), [&v = viewScores](size_t i1, size_t i2) { return v[i1] > v[i2]; }); - // replace the observations + // keep only observations with best scores Observations filteredObservations; for(auto j = 0; j < maxNbObservationsPerLandmark; j++) { From 09c2aba7773b7d98fbd742fd6ec71bd87b0b1cea Mon Sep 17 00:00:00 2001 From: Abdelrahman AL MAROUK Date: Wed, 13 Sep 2023 11:38:54 +0200 Subject: [PATCH 19/37] [filterSfM] add MSE of scores logging between iterations --- src/software/pipeline/main_filterSfM.cpp | 13 ++++++++++++- 1 file changed, 12 insertions(+), 1 deletion(-) diff --git a/src/software/pipeline/main_filterSfM.cpp b/src/software/pipeline/main_filterSfM.cpp index 3c139f2fff..d2d3ae0f3f 100644 --- a/src/software/pipeline/main_filterSfM.cpp +++ b/src/software/pipeline/main_filterSfM.cpp @@ -410,11 +410,22 @@ bool filterObservations(SfMData& sfmData, int maxNbObservationsPerLandmark, int } } // update scores at end of iteration -#pragma omp parallel for + double error = 0.; +#pragma omp parallel for reduction(+:error) for(auto id = 0; id < landmarksData.size(); id++) { + double error_j = 0.; + for(auto j = 0; j < viewScoresData_t[id].size(); j++) + { + const auto& v = viewScoresData_t[id][j] - viewScoresData[id].second[j]; + error_j += v * v; + } + error_j /= viewScoresData_t[id].size(); + error += error_j; viewScoresData[id].second = std::move(viewScoresData_t[id]); } + error /= landmarksData.size(); + ALICEVISION_LOG_INFO("MSE at iteration " << i << ": " << error); } ALICEVISION_LOG_INFO("Propagating neighbors observation scores: done"); From 28299cb1fae6b6f9d17fc15ec4de83ac511cc5be Mon Sep 17 00:00:00 2001 From: Abdelrahman AL MAROUK Date: Wed, 13 Sep 2023 16:43:11 +0200 Subject: [PATCH 20/37] [filterSfM] fix NaN --- src/software/pipeline/main_filterSfM.cpp | 9 ++++++++- 1 file changed, 8 insertions(+), 1 deletion(-) diff --git a/src/software/pipeline/main_filterSfM.cpp b/src/software/pipeline/main_filterSfM.cpp index d2d3ae0f3f..b950b221c4 100644 --- a/src/software/pipeline/main_filterSfM.cpp +++ b/src/software/pipeline/main_filterSfM.cpp @@ -71,7 +71,7 @@ struct LandmarksAdaptator } // Must return the number of data points - inline size_t kdtree_get_point_count() const { return _data.size(); } + inline size_t kdtree_get_point_count() const { return usePtr ? _data_ptr.size() : _data.size(); } // Returns the dim'th component of the idx'th point in the class: inline T kdtree_get_pt(const size_t idx, int dim) const @@ -342,8 +342,12 @@ bool filterObservations(SfMData& sfmData, int maxNbObservationsPerLandmark, int { // weight is the inverse of distance between a landmark and its neighbor w = 1. / std::sqrt(w); + if(std::isinf(w)) + w = std::numeric_limits::max(); total += w; } + if(std::isinf(total)) + total = std::numeric_limits::max(); // normalize weights for(auto& w : weights_) { @@ -401,6 +405,9 @@ bool filterObservations(SfMData& sfmData, int maxNbObservationsPerLandmark, int } } } + // if no common views with neighbor landmarks + if(viewScores_total == 0.) + continue; for(auto j = 0; j < viewScores_acc.size(); j++) { // normalize score and apply influence factor From a7ac8c7f50864c10f4424b749e80d6633ee5cbea Mon Sep 17 00:00:00 2001 From: Abdelrahman AL MAROUK Date: Mon, 18 Sep 2023 11:41:22 +0200 Subject: [PATCH 21/37] [prepareDenseScene] update accessing View from SfmData --- src/software/pipeline/main_prepareDenseScene.cpp | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/software/pipeline/main_prepareDenseScene.cpp b/src/software/pipeline/main_prepareDenseScene.cpp index c32922625a..b416778d3c 100644 --- a/src/software/pipeline/main_prepareDenseScene.cpp +++ b/src/software/pipeline/main_prepareDenseScene.cpp @@ -351,8 +351,9 @@ bool prepareDenseScene(const SfMData& sfmData, if(doMaskLandmarks) { // for the T camera, image alpha should be at least 0.4f * 255 (masking) - maskLandmarks = image::Image(view->getWidth(), view->getHeight(), true, 127); - int r = (int)(landmarksMaskScale * 0.5f * (view->getWidth() + view->getHeight())); + maskLandmarks = + image::Image(view->getImage().getWidth(), view->getImage().getHeight(), true, 127); + int r = (int)(landmarksMaskScale * 0.5f * (view->getImage().getWidth() + view->getImage().getHeight())); const auto& observationsIt = observationsPerView.find(viewId); if(observationsIt != observationsPerView.end()) { From 7edd783391455b0dcb142ba167b2324b6e497e97 Mon Sep 17 00:00:00 2001 From: Abdelrahman AL MAROUK Date: Fri, 22 Sep 2023 07:37:21 +0200 Subject: [PATCH 22/37] [filterSfM] move 2D observations filtering from prepareDenseScene to filterSfM --- src/aliceVision/sfmData/SfMData.cpp | 16 ++ src/aliceVision/sfmData/SfMData.hpp | 10 + src/software/pipeline/main_filterSfM.cpp | 232 ++++++++++++++++-- .../pipeline/main_prepareDenseScene.cpp | 209 +++++----------- 4 files changed, 297 insertions(+), 170 deletions(-) diff --git a/src/aliceVision/sfmData/SfMData.cpp b/src/aliceVision/sfmData/SfMData.cpp index 0cab74643c..b7fc6e3150 100644 --- a/src/aliceVision/sfmData/SfMData.cpp +++ b/src/aliceVision/sfmData/SfMData.cpp @@ -378,5 +378,21 @@ LandmarksPerView getLandmarksPerViews(const SfMData& sfmData) return landmarksPerView; } +ObservationsPerView getObservationsPerViews(SfMData& sfmData) +{ + ObservationsPerView observationsPerView; + for(auto& landIt : sfmData.getLandmarks()) + { + for(const auto& obsIt : landIt.second.observations) + { + IndexT viewId = obsIt.first; + auto& landmarksSet = observationsPerView[viewId]; + landmarksSet.first.push_back(&obsIt.second); + landmarksSet.second.push_back(&landIt.second); + } + } + return observationsPerView; +} + } // namespace sfmData } // namespace aliceVision diff --git a/src/aliceVision/sfmData/SfMData.hpp b/src/aliceVision/sfmData/SfMData.hpp index 8aa5e6c6b1..fd9d7023a5 100644 --- a/src/aliceVision/sfmData/SfMData.hpp +++ b/src/aliceVision/sfmData/SfMData.hpp @@ -599,5 +599,15 @@ using LandmarksPerView = stl::flat_map; LandmarksPerView getLandmarksPerViews(const SfMData& sfmData); +using ObservationsPerView = stl::flat_map, std::vector>>; + +/** + * @brief Get the landmark observations of camera views + * with the corresponding landmarks information. + * @param[in] sfmData: A given SfMData + * @return Observation information per camera view + */ +ObservationsPerView getObservationsPerViews(SfMData& sfmData); + } // namespace sfmData } // namespace aliceVision diff --git a/src/software/pipeline/main_filterSfM.cpp b/src/software/pipeline/main_filterSfM.cpp index b950b221c4..2a76a7ed6f 100644 --- a/src/software/pipeline/main_filterSfM.cpp +++ b/src/software/pipeline/main_filterSfM.cpp @@ -18,6 +18,8 @@ #include #include #include +#include +#include #include #include @@ -42,9 +44,42 @@ using namespace aliceVision::sfmDataIO; namespace po = boost::program_options; namespace fs = boost::filesystem; +using namespace boost::accumulators; static const std::size_t MAX_LEAF_ELEMENTS = 64; +struct ObservationsAdaptator +{ + using Derived = ObservationsAdaptator; //!< In this case the dataset class is myself. + using T = double; + + /// CRTP helper method + inline const Derived& derived() const { return *static_cast(this); } + /// CRTP helper method + inline Derived& derived() { return *static_cast(this); } + + const std::vector _data; + ObservationsAdaptator(const std::vector& data) + : _data(data) + { + } + + // Must return the number of data points + inline size_t kdtree_get_point_count() const { return _data.size(); } + + // Returns the dim'th component of the idx'th point in the class: + inline T kdtree_get_pt(const size_t idx, int dim) const { return _data[idx]->x(dim); } + + // Optional bounding-box computation: return false to default to a standard bbox computation loop. + // Return true if the BBOX was already computed by the class and returned in "bb" so it can be avoided to redo it + // again. Look at bb.size() to find out the expected dimensionality (e.g. 2 or 3 for point clouds) + template + bool kdtree_get_bbox(BBOX& bb) const + { + return false; + } +}; + struct LandmarksAdaptator { using Derived = LandmarksAdaptator; //!< In this case the dataset class is myself. @@ -89,16 +124,42 @@ struct LandmarksAdaptator } }; -using KdTree = nanoflann::KDTreeSingleIndexAdaptor< +using KdTree2D = nanoflann::KDTreeSingleIndexAdaptor< + nanoflann::L2_Simple_Adaptor, + ObservationsAdaptator, 2 /* dim */, size_t>; + +using KdTree3D = nanoflann::KDTreeSingleIndexAdaptor< nanoflann::L2_Simple_Adaptor, - LandmarksAdaptator, - 3, /* dim */ - size_t ->; - -/** - * A result-set class used when performing a radius based search. - */ + LandmarksAdaptator, 3 /* dim */,size_t>; + +class RadiusKnnSearch +{ +public: + const double _radius_sq; + const int _nb_neighbors; + int nb_found = 0; + + inline RadiusKnnSearch(double radius, int k) + : _radius_sq(radius * radius) + , _nb_neighbors(k) + { + } + + inline bool full() const { return nb_found == _nb_neighbors; } + + inline bool addPoint(double dist, IndexT index) + { + if(dist < _radius_sq) + { + nb_found++; + return nb_found < _nb_neighbors; + } + return true; + } + + inline double worstDist() const { return _radius_sq; } +}; + class PixSizeSearch { public: @@ -118,6 +179,7 @@ class PixSizeSearch inline bool addPoint(double dist, size_t index) { + // strict comparison because a point in the tree can be a neighbor of itself if(dist < _radius_sq && _pixSize[index] < _pixSize_i) { found = true; @@ -169,7 +231,7 @@ bool filterLandmarks(SfMData& sfmData, double radiusScale, bool useFeatureScale, { filteredLandmarks[newIdx] = std::make_pair(newIdx, landmarksData[newIdx]); } - sfmData.getLandmarks() = Landmarks(filteredLandmarks.begin(), filteredLandmarks.end()); + sfmData.getLandmarks() = std::move(Landmarks(filteredLandmarks.begin(), filteredLandmarks.end())); return true; } @@ -203,9 +265,9 @@ bool filterLandmarks(SfMData& sfmData, double radiusScale, bool useFeatureScale, //// sort landmarks by descending pixSize order //std::stable_sort(landmarksPixSize.begin(), landmarksPixSize.end(), std::greater<>{}); - ALICEVISION_LOG_INFO("Build nanoflann KdTree index."); + ALICEVISION_LOG_INFO("Build nanoflann KdTree index for landmarks in 3D."); LandmarksAdaptator data(landmarksData); - KdTree tree(3, data, nanoflann::KDTreeSingleIndexAdaptorParams(MAX_LEAF_ELEMENTS)); + KdTree3D tree(3, data, nanoflann::KDTreeSingleIndexAdaptorParams(MAX_LEAF_ELEMENTS)); tree.buildIndex(); ALICEVISION_LOG_INFO("KdTree created for " << landmarksData.size() << " points."); std::vector toRemove(landmarksData.size(), false); @@ -249,7 +311,7 @@ bool filterLandmarks(SfMData& sfmData, double radiusScale, bool useFeatureScale, return true; } -bool filterObservations(SfMData& sfmData, int maxNbObservationsPerLandmark, int nbNeighbors, double neighborsInfluence, +bool filterObservations3D(SfMData& sfmData, int maxNbObservationsPerLandmark, int nbNeighbors3D, double neighborsInfluence, int nbIterations) { // store in vector for faster access @@ -316,13 +378,13 @@ bool filterObservations(SfMData& sfmData, int maxNbObservationsPerLandmark, int ALICEVISION_LOG_INFO("Computing initial observation scores based on distance to observing view: done"); ALICEVISION_LOG_INFO("Computing landmark neighbors and distance-based weights: started"); - ALICEVISION_LOG_INFO("Build nanoflann KdTree index."); + ALICEVISION_LOG_INFO("Build nanoflann KdTree index for landmarks in 3D."); LandmarksAdaptator dataAdaptor(landmarksData); - KdTree tree(3, dataAdaptor, nanoflann::KDTreeSingleIndexAdaptorParams(MAX_LEAF_ELEMENTS)); + KdTree3D tree(3, dataAdaptor, nanoflann::KDTreeSingleIndexAdaptorParams(MAX_LEAF_ELEMENTS)); tree.buildIndex(); ALICEVISION_LOG_INFO("KdTree created for " << landmarksData.size() << " points."); // note that the landmark is a neighbor to itself with zero distance, hence the +/- 1 - int nbNeighbors_ = std::min(nbNeighbors, static_cast(landmarksData.size() - 1)) + 1; + int nbNeighbors_ = std::min(nbNeighbors3D, static_cast(landmarksData.size() - 1)) + 1; // contains the neighbor landmarks ids with their corresponding weights std::vector, std::vector>> neighborsData(landmarksData.size()); #pragma omp parallel for @@ -372,7 +434,7 @@ bool filterObservations(SfMData& sfmData, int maxNbObservationsPerLandmark, int // accumulator for normalisation double viewScores_total = 0.; auto& [indices_, weights_] = neighborsData[id]; - for(auto j = 0; j < nbNeighbors; j++) + for(auto j = 0; j < nbNeighbors3D; j++) { const auto& neighborId = indices_[j]; const auto& neighborWeight = weights_[j]; @@ -467,19 +529,118 @@ bool filterObservations(SfMData& sfmData, int maxNbObservationsPerLandmark, int return true; } +bool filterObservations2D(SfMData& sfmData, int nbNeighbors2D, float percentile, + HashMap& estimatedRadii) +{ + std::set viewIds = sfmData.getValidViews(); + std::vector estimatedRadii_(viewIds.size(), -1.); + auto observationsPerView = getObservationsPerViews(sfmData); + +#pragma omp parallel for + for(int i = 0; i < viewIds.size(); ++i) + { + auto itView = viewIds.begin(); + std::advance(itView, i); + const IndexT viewId = *itView; + + auto& observationsIt = observationsPerView.find(viewId); + if(observationsIt != observationsPerView.end()) + { + auto& observations = observationsIt->second.first; + auto& landmarks = observationsIt->second.second; + + ALICEVISION_LOG_INFO("Build nanoflann KdTree index for projected landmarks in 2D."); + ObservationsAdaptator data(observations); + KdTree2D tree(2, data, nanoflann::KDTreeSingleIndexAdaptorParams(MAX_LEAF_ELEMENTS)); + tree.buildIndex(); + ALICEVISION_LOG_INFO("KdTree created for " << observations.size() << " points."); + + // note that the observation is a neighbor to itself with zero distance, hence the +/- 1 + size_t nbNeighbors_ = std::min(static_cast(nbNeighbors2D), observations.size() - 1) + 1; + std::vector means(observations.size()); + const std::size_t cacheSize = 1000; + accumulator_set, tag::mean>> acc(tag::tail::cache_size = + cacheSize); + for(auto j = 0; j < observations.size(); j++) + { + const auto& obs = *observations[j]; + std::vector indices_(nbNeighbors_); + std::vector distances_(nbNeighbors_); + tree.knnSearch(obs.x.data(), nbNeighbors_, &indices_[0], &distances_[0]); + + std::transform(distances_.begin(), distances_.end(), distances_.begin(), + static_cast(std::sqrt)); + const auto& mean = std::accumulate(distances_.begin(), distances_.end(), 0.0) / (nbNeighbors_ - 1); + means[j] = mean; + acc(mean); + } + double mean_max = std::numeric_limits::max(); + if(percentile != 1.f) + mean_max = quantile(acc, quantile_probability = percentile); + estimatedRadii_[i] = mean(acc); + + std::vector filteredObservations; + std::vector filteredLandmarks; + filteredObservations.reserve(observations.size()); + filteredLandmarks.reserve(landmarks.size()); + for(auto j = 0; j < observations.size(); j++) + if (means[j] < mean_max) + { + filteredObservations.push_back(observations[j]); + filteredLandmarks.push_back(landmarks[j]); + } + filteredObservations.shrink_to_fit(); + filteredLandmarks.shrink_to_fit(); + observations = std::move(filteredObservations); + landmarks = std::move(filteredLandmarks); + } + } + + for(auto& landmark : sfmData.getLandmarks()) + { + landmark.second.observations.clear(); + } + + for(int i = 0; i < viewIds.size(); ++i) + { + auto itView = viewIds.begin(); + std::advance(itView, i); + const IndexT viewId = *itView; + + if(estimatedRadii_[i] != -1.) + estimatedRadii[viewId] = estimatedRadii_[i]; + + auto& observationsIt = observationsPerView.find(viewId); + if(observationsIt != observationsPerView.end()) + { + auto& observations = observationsIt->second.first; + auto& landmarks = observationsIt->second.second; + for(int j = 0; j < observations.size(); j++) + { + landmarks[j]->observations[viewId] = *observations[j]; + } + } + } + + return true; +} + int aliceVision_main(int argc, char *argv[]) { // command-line parameters std::string inputSfmFilename; std::string outputSfmFilename; + std::string outputRadiiFilename; int maxNbObservationsPerLandmark = 2; int minNbObservationsPerLandmark = 5; double radiusScale = 2; bool useFeatureScale = true; - int nbNeighbors = 10; + int nbNeighbors3D = 10; double neighborsInfluence = 0.5; int nbIterations = 5; + int nbNeighbors2D = 5; + float percentile = 0.95; // user optional parameters std::vector featuresFolders; @@ -503,19 +664,25 @@ int aliceVision_main(int argc, char *argv[]) "Scale factor applied to pixel size based radius filter applied to landmarks.") ("useFeatureScale", po::value(&useFeatureScale)->default_value(useFeatureScale), "If true, use feature scale for computing pixel size. Otherwise, use a scale of 1 pixel.") - ("nbNeighbors", po::value(&nbNeighbors)->default_value(nbNeighbors), + ("nbNeighbors3D", po::value(&nbNeighbors3D)->default_value(nbNeighbors3D), "Number of neighbor landmarks used in making the decision for best observations.") ("neighborsInfluence", po::value(&neighborsInfluence)->default_value(neighborsInfluence), "Specifies how much influential the neighbors are in selecting the best observations." "Between 0. and 1., the closer to 1., the more influencial the neighborhood is.") ("nbIterations", po::value(&nbIterations)->default_value(nbIterations), "Number of iterations to propagate neighbors information.") + ("nbNeighbors2D", po::value(&nbNeighbors2D)->default_value(nbNeighbors2D), + "Number of neighbor observations to be considered for the landmarks-based masking.") + ("percentile", po::value(&percentile)->default_value(percentile), + "TODO.") ("featuresFolders,f", po::value>(&featuresFolders)->multitoken(), "Path to folder(s) containing the extracted features.") ("matchesFolders,m", po::value>(&matchesFolders)->multitoken(), "Path to folder(s) in which computed matches are stored.") ("describerTypes,d", po::value(&describerTypesName)->default_value(describerTypesName), - feature::EImageDescriberType_informations().c_str()); + feature::EImageDescriberType_informations().c_str()) + ("outputRadiiFile", po::value(&outputRadiiFilename)->default_value(outputRadiiFilename), + "Output Radii file containing the estimated projection radius of observations per view."); CmdLine cmdline("AliceVision SfM filtering."); // TODO add description cmdline.add(requiredParams); @@ -551,11 +718,30 @@ int aliceVision_main(int argc, char *argv[]) if(maxNbObservationsPerLandmark > 0) { - ALICEVISION_LOG_INFO("Filtering observations: started."); - filterObservations(sfmData, maxNbObservationsPerLandmark, nbNeighbors, neighborsInfluence, nbIterations); - ALICEVISION_LOG_INFO("Filtering observations: done."); + ALICEVISION_LOG_INFO("Filtering observations in 3D: started."); + filterObservations3D(sfmData, maxNbObservationsPerLandmark, nbNeighbors3D, neighborsInfluence, nbIterations); + ALICEVISION_LOG_INFO("Filtering observations in 3D: done."); } + if(nbNeighbors2D > 0 || percentile < 1.f) + { + HashMap estimatedRadii; + ALICEVISION_LOG_INFO("Filtering observations in 2D: started."); + filterObservations2D(sfmData, nbNeighbors2D, percentile, estimatedRadii); + ALICEVISION_LOG_INFO("Filtering observations in 2D: done."); + + if(outputRadiiFilename.empty()) + outputRadiiFilename = (fs::path(outputSfmFilename).parent_path() / "radii.txt").string(); + std::ofstream fs(outputRadiiFilename, std::ios::out); + if(!fs.is_open()) + ALICEVISION_LOG_WARNING("Unable to create the radii file " << outputRadiiFilename); + else + { + for(const auto& radius : estimatedRadii) + fs << radius.first << "\t" << radius.second << std::endl; + fs.close(); + } + } sfmDataIO::Save(sfmData, outputSfmFilename, sfmDataIO::ESfMData::ALL); return EXIT_SUCCESS; diff --git a/src/software/pipeline/main_prepareDenseScene.cpp b/src/software/pipeline/main_prepareDenseScene.cpp index b416778d3c..9170b014fb 100644 --- a/src/software/pipeline/main_prepareDenseScene.cpp +++ b/src/software/pipeline/main_prepareDenseScene.cpp @@ -14,12 +14,8 @@ #include #include -#include "nanoflann.hpp" - #include #include -#include -#include #include #include @@ -45,91 +41,6 @@ using namespace aliceVision::sfmDataIO; namespace po = boost::program_options; namespace fs = std::filesystem; -using namespace boost::accumulators; - -static const std::size_t MAX_LEAF_ELEMENTS = 64; - -struct ObservationsAdaptator -{ - using Derived = ObservationsAdaptator; //!< In this case the dataset class is myself. - using T = double; - - /// CRTP helper method - inline const Derived& derived() const { return *static_cast(this); } - /// CRTP helper method - inline Derived& derived() { return *static_cast(this); } - - const std::vector _data; - ObservationsAdaptator(const std::vector& data) - : _data(data) - { - } - - // Must return the number of data points - inline size_t kdtree_get_point_count() const { return _data.size(); } - - // Returns the dim'th component of the idx'th point in the class: - inline T kdtree_get_pt(const size_t idx, int dim) const { return _data[idx]->x(dim); } - - // Optional bounding-box computation: return false to default to a standard bbox computation loop. - // Return true if the BBOX was already computed by the class and returned in "bb" so it can be avoided to redo it - // again. Look at bb.size() to find out the expected dimensionality (e.g. 2 or 3 for point clouds) - template - bool kdtree_get_bbox(BBOX& bb) const - { - return false; - } -}; - -using KdTree = nanoflann::KDTreeSingleIndexAdaptor, - ObservationsAdaptator, 2, /* dim */ - size_t>; - -/** - * A result-set class used when performing a radius based search. - */ -class RadiusKnnSearch -{ -public: - const double _radius_sq; - const int _nb_neighbors; - int nb_found = 0; - - inline RadiusKnnSearch(double radius, int k) - : _radius_sq(radius * radius) - , _nb_neighbors(k) - { - } - - inline bool full() const { return nb_found == _nb_neighbors; } - - inline bool addPoint(double dist, IndexT index) - { - if(dist < _radius_sq) - { - nb_found++; - return nb_found < _nb_neighbors; - } - return true; - } - - inline double worstDist() const { return _radius_sq; } -}; - -const auto getObservationsPerViews(const SfMData& sfmData) -{ - stl::flat_map> observationsPerView; - for(const auto& landIt : sfmData.getLandmarks()) - { - for(const auto& obsIt : landIt.second.observations) - { - IndexT viewId = obsIt.first; - auto& landmarksSet = observationsPerView[viewId]; - landmarksSet.push_back(&obsIt.second); - } - } - return observationsPerView; -} template void process(const std::string& dstColorImage, @@ -172,7 +83,7 @@ void process(const std::string& dstColorImage, } } -bool prepareDenseScene(const SfMData& sfmData, +bool prepareDenseScene(SfMData& sfmData, const std::vector& imagesFolders, const std::vector& masksFolders, const std::string& maskExtension, @@ -184,8 +95,7 @@ bool prepareDenseScene(const SfMData& sfmData, bool saveMatricesFiles, bool evCorrection, float landmarksMaskScale, - int nbNeighborObservations, - float percentile) + std::string inputRadiiFilename) { // defined view Ids std::set viewIds; @@ -221,7 +131,42 @@ bool prepareDenseScene(const SfMData& sfmData, const double medianCameraExposure = sfmData.getMedianCameraExposureSetting().getExposure(); ALICEVISION_LOG_INFO("Median Camera Exposure: " << medianCameraExposure << ", Median EV: " << std::log2(1.0 / medianCameraExposure)); - const auto& observationsPerView = getObservationsPerViews(sfmData); + bool doMaskLandmarks = landmarksMaskScale > 0.f; + sfmData::ObservationsPerView observationsPerView; + HashMap estimatedRadii; + if (doMaskLandmarks) + { + observationsPerView = std::move(getObservationsPerViews(sfmData)); + + if (!inputRadiiFilename.empty()) + { + std::stringstream stream; + std::string line; + IndexT viewId; + double radius; + + std::fstream fs(inputRadiiFilename, std::ios::in); + if(!fs.is_open()) + { + ALICEVISION_LOG_WARNING("Unable to open the radii file " << inputRadiiFilename + << "\nDefaulting to using image size."); + } + else + { + while(!fs.eof()) + { + std::getline(fs, line); + stream.clear(); + stream.str(line); + stream >> viewId; + stream >> radius; + estimatedRadii[viewId] = radius; + } + fs.close(); + } + } + + } #pragma omp parallel for for (int i = 0; i < viewIds.size(); ++i) @@ -347,63 +292,35 @@ bool prepareDenseScene(const SfMData& sfmData, } image::Image maskLandmarks; - bool doMaskLandmarks = nbNeighborObservations > 0 || percentile < 1.f; if(doMaskLandmarks) { // for the T camera, image alpha should be at least 0.4f * 255 (masking) maskLandmarks = image::Image(view->getImage().getWidth(), view->getImage().getHeight(), true, 127); - int r = (int)(landmarksMaskScale * 0.5f * (view->getImage().getWidth() + view->getImage().getHeight())); + double radius; + const auto& it = estimatedRadii.find(viewId); + if(it != estimatedRadii.end()) + radius = it->second; + else + radius = 0.5 * (view->getImage().getWidth() + view->getImage().getHeight()); + int r = (int)(landmarksMaskScale * radius); const auto& observationsIt = observationsPerView.find(viewId); if(observationsIt != observationsPerView.end()) { - const auto& observations = observationsIt->second; - - ALICEVISION_LOG_INFO("Build nanoflann KdTree index."); - ObservationsAdaptator data(observations); - KdTree tree(2, data, nanoflann::KDTreeSingleIndexAdaptorParams(MAX_LEAF_ELEMENTS)); - tree.buildIndex(); - ALICEVISION_LOG_INFO("KdTree created for " << observations.size() << " points."); - - int n = std::min(nbNeighborObservations, static_cast(observations.size() - 1)) + 1; - // note that the observation is a neighbor to itself with zero distance, hence the +/- 1 - std::vector means(observations.size()); - const std::size_t cacheSize = 1000; - accumulator_set, tag::mean>> acc( - tag::tail::cache_size = cacheSize); + const auto& observations = observationsIt->second.first; int j = 0; - for (const auto& observation : observations) - { - const auto& obs = *observation; - std::vector indices_(n); - std::vector distances_(n); - tree.knnSearch(obs.x.data(), n, &indices_[0], &distances_[0]); - - std::transform(distances_.begin(), distances_.end(), distances_.begin(), - static_cast(std::sqrt)); - const auto& mean = std::accumulate(distances_.begin(), distances_.end(), 0.0) / (n - 1); - means[j++] = mean; - acc(mean); - } - double mean_max = std::numeric_limits::max(); - if(percentile != 1.f) - mean_max = quantile(acc, quantile_probability = percentile); - r = mean(acc); - - j = 0; for(const auto& observation : observations) { const auto& obs = *observation; - if(means[j] < mean_max) - for(int y = std::max(obs.x.y() - r, 0.); - y <= std::min(obs.x.y() + r, (double)maskLandmarks.Height() - 1); y++) + for(int y = std::max(obs.x.y() - r, 0.); + y <= std::min(obs.x.y() + r, (double)maskLandmarks.Height() - 1); y++) + { + for(int x = std::max(obs.x.x() - r, 0.); + x <= std::min(obs.x.x() + r, (double)maskLandmarks.Width() - 1); x++) { - for(int x = std::max(obs.x.x() - r, 0.); - x <= std::min(obs.x.x() + r, (double)maskLandmarks.Width() - 1); x++) - { - maskLandmarks(y, x) = std::numeric_limits::max(); - } + maskLandmarks(y, x) = std::numeric_limits::max(); } + } j++; } } @@ -452,9 +369,8 @@ int aliceVision_main(int argc, char* argv[]) bool saveMetadata = true; bool saveMatricesTxtFiles = false; bool evCorrection = false; + std::string inputRadiiFilename; float landmarksMaskScale = 0.f; - int nbNeighborObservations = 5; - float percentile = 0.95; // clang-format off po::options_description requiredParams("Required parameters"); @@ -487,14 +403,14 @@ int aliceVision_main(int argc, char* argv[]) ("evCorrection", po::value(&evCorrection)->default_value(evCorrection), "Correct exposure value.") ("landmarksMaskScale", po::value(&landmarksMaskScale)->default_value(landmarksMaskScale), - "Scale (relative to image size) of the projection of landmarks to mask images for depth computation.\n" - "If 0, masking using landmarks will not be used.") - ("nbNeighborObservations", po::value(&nbNeighborObservations)->default_value(nbNeighborObservations), - "Number of neighbor observations to be considered for the landmarks-based masking.") - ("percentile", po::value(&percentile)->default_value(percentile), - "TODO."); + "Scale of the projection of landmarks to mask images for depth computation.\n" + "If 0, masking using landmarks will not be used.\n" + "Otherwise, it's used to scale the projection radius \n" + "(either specified by `inputRadiiFile` or by image size if the former is not given).") + ("inputRadiiFile", po::value(&inputRadiiFilename)->default_value(inputRadiiFilename), + "Input Radii file containing the estimated projection radius of landmarks per view. \n" + "If not specified, image size will be used to specify the radius."); // clang-format on - CmdLine cmdline("AliceVision prepareDenseScene"); cmdline.add(requiredParams); cmdline.add(optionalParams); @@ -558,8 +474,7 @@ int aliceVision_main(int argc, char* argv[]) saveMatricesTxtFiles, evCorrection, landmarksMaskScale, - nbNeighborObservations, - percentile)) + inputRadiiFilename)) return EXIT_SUCCESS; return EXIT_FAILURE; From 2d5d4e5326ca727d99589b05b358c0152d10faa2 Mon Sep 17 00:00:00 2001 From: Abdelrahman AL MAROUK Date: Fri, 22 Sep 2023 08:03:55 +0200 Subject: [PATCH 23/37] [filterSfM] add limit check on estimated 2D radius --- src/software/pipeline/main_filterSfM.cpp | 103 ++++++++++++----------- 1 file changed, 56 insertions(+), 47 deletions(-) diff --git a/src/software/pipeline/main_filterSfM.cpp b/src/software/pipeline/main_filterSfM.cpp index 2a76a7ed6f..bc8a522732 100644 --- a/src/software/pipeline/main_filterSfM.cpp +++ b/src/software/pipeline/main_filterSfM.cpp @@ -544,56 +544,65 @@ bool filterObservations2D(SfMData& sfmData, int nbNeighbors2D, float percentile, const IndexT viewId = *itView; auto& observationsIt = observationsPerView.find(viewId); - if(observationsIt != observationsPerView.end()) - { - auto& observations = observationsIt->second.first; - auto& landmarks = observationsIt->second.second; + if(observationsIt == observationsPerView.end()) + continue; - ALICEVISION_LOG_INFO("Build nanoflann KdTree index for projected landmarks in 2D."); - ObservationsAdaptator data(observations); - KdTree2D tree(2, data, nanoflann::KDTreeSingleIndexAdaptorParams(MAX_LEAF_ELEMENTS)); - tree.buildIndex(); - ALICEVISION_LOG_INFO("KdTree created for " << observations.size() << " points."); - - // note that the observation is a neighbor to itself with zero distance, hence the +/- 1 - size_t nbNeighbors_ = std::min(static_cast(nbNeighbors2D), observations.size() - 1) + 1; - std::vector means(observations.size()); - const std::size_t cacheSize = 1000; - accumulator_set, tag::mean>> acc(tag::tail::cache_size = - cacheSize); - for(auto j = 0; j < observations.size(); j++) + auto& observations = observationsIt->second.first; + auto& landmarks = observationsIt->second.second; + + ALICEVISION_LOG_INFO("Build nanoflann KdTree index for projected landmarks in 2D."); + ObservationsAdaptator data(observations); + KdTree2D tree(2, data, nanoflann::KDTreeSingleIndexAdaptorParams(MAX_LEAF_ELEMENTS)); + tree.buildIndex(); + ALICEVISION_LOG_INFO("KdTree created for " << observations.size() << " points."); + + // note that the observation is a neighbor to itself with zero distance, hence the +/- 1 + size_t nbNeighbors_ = std::min(static_cast(nbNeighbors2D), observations.size() - 1) + 1; + std::vector means(observations.size()); + const std::size_t cacheSize = 1000; + accumulator_set, tag::mean>> acc(tag::tail::cache_size = + cacheSize); + for(auto j = 0; j < observations.size(); j++) + { + const auto& obs = *observations[j]; + std::vector indices_(nbNeighbors_); + std::vector distances_(nbNeighbors_); + tree.knnSearch(obs.x.data(), nbNeighbors_, &indices_[0], &distances_[0]); + + std::transform(distances_.begin(), distances_.end(), distances_.begin(), + static_cast(std::sqrt)); + const auto& mean = std::accumulate(distances_.begin(), distances_.end(), 0.0) / (nbNeighbors_ - 1); + means[j] = mean; + acc(mean); + } + double mean_max = std::numeric_limits::max(); + if(percentile != 1.f) + mean_max = quantile(acc, quantile_probability = percentile); + + double radius = mean(acc); + // check if estimated radius is too large + { + const View& view = *(sfmData.getViews().at(viewId)); + double radiusMax = 0.15 * 0.5 * (view.getImage().getWidth() + view.getImage().getHeight()); + if(radius > radiusMax) + radius = radiusMax; + } + estimatedRadii_[i] = radius; + + std::vector filteredObservations; + std::vector filteredLandmarks; + filteredObservations.reserve(observations.size()); + filteredLandmarks.reserve(landmarks.size()); + for(auto j = 0; j < observations.size(); j++) + if(means[j] < mean_max) { - const auto& obs = *observations[j]; - std::vector indices_(nbNeighbors_); - std::vector distances_(nbNeighbors_); - tree.knnSearch(obs.x.data(), nbNeighbors_, &indices_[0], &distances_[0]); - - std::transform(distances_.begin(), distances_.end(), distances_.begin(), - static_cast(std::sqrt)); - const auto& mean = std::accumulate(distances_.begin(), distances_.end(), 0.0) / (nbNeighbors_ - 1); - means[j] = mean; - acc(mean); + filteredObservations.push_back(observations[j]); + filteredLandmarks.push_back(landmarks[j]); } - double mean_max = std::numeric_limits::max(); - if(percentile != 1.f) - mean_max = quantile(acc, quantile_probability = percentile); - estimatedRadii_[i] = mean(acc); - - std::vector filteredObservations; - std::vector filteredLandmarks; - filteredObservations.reserve(observations.size()); - filteredLandmarks.reserve(landmarks.size()); - for(auto j = 0; j < observations.size(); j++) - if (means[j] < mean_max) - { - filteredObservations.push_back(observations[j]); - filteredLandmarks.push_back(landmarks[j]); - } - filteredObservations.shrink_to_fit(); - filteredLandmarks.shrink_to_fit(); - observations = std::move(filteredObservations); - landmarks = std::move(filteredLandmarks); - } + filteredObservations.shrink_to_fit(); + filteredLandmarks.shrink_to_fit(); + observations = std::move(filteredObservations); + landmarks = std::move(filteredLandmarks); } for(auto& landmark : sfmData.getLandmarks()) From 47048d8a56720756219db01a931ee74c442a742d Mon Sep 17 00:00:00 2001 From: Abdelrahman AL MAROUK Date: Fri, 22 Sep 2023 09:45:28 +0200 Subject: [PATCH 24/37] [filterSfM] improve landmarks filtering remove landmarks that don't have enough common observations with their neighbors --- src/software/pipeline/main_filterSfM.cpp | 105 +++++++++++++++++++++++ 1 file changed, 105 insertions(+) diff --git a/src/software/pipeline/main_filterSfM.cpp b/src/software/pipeline/main_filterSfM.cpp index bc8a522732..5b4cf00f84 100644 --- a/src/software/pipeline/main_filterSfM.cpp +++ b/src/software/pipeline/main_filterSfM.cpp @@ -224,6 +224,111 @@ bool filterLandmarks(SfMData& sfmData, double radiusScale, bool useFeatureScale, } } + { + ALICEVISION_LOG_INFO("Build nanoflann KdTree index for landmarks in 3D."); + LandmarksAdaptator data(landmarksData); + KdTree3D tree(3, data, nanoflann::KDTreeSingleIndexAdaptorParams(MAX_LEAF_ELEMENTS)); + tree.buildIndex(); + ALICEVISION_LOG_INFO("KdTree created for " << landmarksData.size() << " points."); + + // TODO as parameter + int nbNeighbors3D = 10; + // note that the landmark is a neighbor to itself with zero distance, hence the +/- 1 + int nbNeighbors_ = std::min(nbNeighbors3D, static_cast(landmarksData.size() - 1)) + 1; + + // contains the observing view ids and neighbors for each landmark + std::vector, std::vector>> viewData(landmarksData.size()); + +#pragma omp parallel for + for(auto i = 0; i < landmarksData.size(); i++) + { + const sfmData::Landmark& landmark = landmarksData[i]; + const auto& nbObservations = landmark.observations.size(); + auto& [viewIds, neighbors] = viewData[i]; + viewIds.reserve(nbObservations); + for(const auto& observationPair : landmark.observations) + { + const IndexT viewId = observationPair.first; + viewIds.push_back(viewId); + } + // sort by ascending view id order for consequent faster access + std::stable_sort(viewIds.begin(), viewIds.end()); + + neighbors.resize(nbNeighbors_); + std::vector weights_(nbNeighbors_); + tree.knnSearch(landmark.X.data(), nbNeighbors_, &neighbors[0], &weights_[0]); + // a landmark is a neighbor to itself with zero distance, remove it + neighbors.erase(neighbors.begin()); + } + + std::vector toRemove(landmarksData.size(), false); + size_t nbToRemove = 0; + // TODO as parameter + int maxNbObservationsPerLandmark = 2; + const auto initialNbLandmarks = sfmData.getLandmarks().size(); +#pragma omp parallel for reduction(+ : nbToRemove) + for(auto i = 0; i < landmarksData.size(); i++) + { + int score = 0; + auto& [viewIds, neighbors] = viewData[i]; + for(auto j = 0; j < neighbors.size(); j++) + { + int scoreNeighbor = 0; + const auto& neighborId = neighbors[j]; + const auto& viewIds_neighbor = viewData[neighborId].first; + // loop over common views + auto viewIds_it = viewIds.begin(); + auto viewIds_neighbor_it = viewIds_neighbor.begin(); + while(viewIds_it != viewIds.end() && viewIds_neighbor_it != viewIds_neighbor.end()) + { + if(*viewIds_it < *viewIds_neighbor_it) + { + ++viewIds_it; + } + else + { + // if same view + if(!(*viewIds_neighbor_it < *viewIds_it)) + { + scoreNeighbor++; + if(scoreNeighbor == maxNbObservationsPerLandmark) + { + score++; + break; + } + ++viewIds_it; + } + ++viewIds_neighbor_it; + } + } + } + if(score < nbNeighbors3D / 2) + { + toRemove[i] = true; + nbToRemove++; + } + } + ALICEVISION_LOG_INFO("Identified " << nbToRemove << " landmarks to remove out of " << initialNbLandmarks + << ", i.e. " << (nbToRemove * 100.f / initialNbLandmarks) << " % "); + ALICEVISION_LOG_INFO("Identifying landmarks to remove: done."); + + ALICEVISION_LOG_INFO("Removing landmarks based on TODO: started."); + std::vector> filteredLandmarks(landmarksData.size() - nbToRemove); + std::vector landmarksData_filtered(landmarksData.size() - nbToRemove); + IndexT newIdx = 0; + for(size_t i = 0; i < landmarksData.size(); i++) + { + if(!toRemove[i]) + { + landmarksData_filtered[newIdx] = landmarksData[i]; + filteredLandmarks[newIdx++] = std::make_pair(newIdx, landmarksData[i]); + } + } + sfmData.getLandmarks() = std::move(Landmarks(filteredLandmarks.begin(), filteredLandmarks.end())); + landmarksData = std::move(landmarksData_filtered); + ALICEVISION_LOG_INFO("Removing landmarks based on TODO: done."); + } + if (radiusScale <= 0) { std::vector> filteredLandmarks(landmarksData.size()); From 95dbe2ca4cd21f290e776edbaf8295f80487cc4a Mon Sep 17 00:00:00 2001 From: Abdelrahman AL MAROUK Date: Fri, 22 Sep 2023 16:08:50 +0200 Subject: [PATCH 25/37] [filterSfM] fix reference issue and refactor code --- src/aliceVision/sfmData/SfMData.cpp | 16 ------- src/aliceVision/sfmData/SfMData.hpp | 10 ----- src/software/pipeline/main_filterSfM.cpp | 42 +++++++++++++++---- .../pipeline/main_prepareDenseScene.cpp | 28 +++++++++++-- 4 files changed, 58 insertions(+), 38 deletions(-) diff --git a/src/aliceVision/sfmData/SfMData.cpp b/src/aliceVision/sfmData/SfMData.cpp index b7fc6e3150..0cab74643c 100644 --- a/src/aliceVision/sfmData/SfMData.cpp +++ b/src/aliceVision/sfmData/SfMData.cpp @@ -378,21 +378,5 @@ LandmarksPerView getLandmarksPerViews(const SfMData& sfmData) return landmarksPerView; } -ObservationsPerView getObservationsPerViews(SfMData& sfmData) -{ - ObservationsPerView observationsPerView; - for(auto& landIt : sfmData.getLandmarks()) - { - for(const auto& obsIt : landIt.second.observations) - { - IndexT viewId = obsIt.first; - auto& landmarksSet = observationsPerView[viewId]; - landmarksSet.first.push_back(&obsIt.second); - landmarksSet.second.push_back(&landIt.second); - } - } - return observationsPerView; -} - } // namespace sfmData } // namespace aliceVision diff --git a/src/aliceVision/sfmData/SfMData.hpp b/src/aliceVision/sfmData/SfMData.hpp index fd9d7023a5..8aa5e6c6b1 100644 --- a/src/aliceVision/sfmData/SfMData.hpp +++ b/src/aliceVision/sfmData/SfMData.hpp @@ -599,15 +599,5 @@ using LandmarksPerView = stl::flat_map; LandmarksPerView getLandmarksPerViews(const SfMData& sfmData); -using ObservationsPerView = stl::flat_map, std::vector>>; - -/** - * @brief Get the landmark observations of camera views - * with the corresponding landmarks information. - * @param[in] sfmData: A given SfMData - * @return Observation information per camera view - */ -ObservationsPerView getObservationsPerViews(SfMData& sfmData); - } // namespace sfmData } // namespace aliceVision diff --git a/src/software/pipeline/main_filterSfM.cpp b/src/software/pipeline/main_filterSfM.cpp index 5b4cf00f84..ea695f6c35 100644 --- a/src/software/pipeline/main_filterSfM.cpp +++ b/src/software/pipeline/main_filterSfM.cpp @@ -58,8 +58,8 @@ struct ObservationsAdaptator /// CRTP helper method inline Derived& derived() { return *static_cast(this); } - const std::vector _data; - ObservationsAdaptator(const std::vector& data) + const std::vector _data; + ObservationsAdaptator(const std::vector& data) : _data(data) { } @@ -68,7 +68,7 @@ struct ObservationsAdaptator inline size_t kdtree_get_point_count() const { return _data.size(); } // Returns the dim'th component of the idx'th point in the class: - inline T kdtree_get_pt(const size_t idx, int dim) const { return _data[idx]->x(dim); } + inline T kdtree_get_pt(const size_t idx, int dim) const { return _data[idx].x(dim); } // Optional bounding-box computation: return false to default to a standard bbox computation loop. // Return true if the BBOX was already computed by the class and returned in "bb" so it can be avoided to redo it @@ -191,6 +191,30 @@ class PixSizeSearch inline double worstDist() const { return _radius_sq; } }; +using ObservationsPerView = stl::flat_map, std::vector>>; + +/** + * @brief Get the landmark observations of camera views + * with the corresponding landmarks information. + * @param[in] sfmData: A given SfMData + * @return Observation information per camera view + */ +ObservationsPerView getObservationsPerViews(SfMData& sfmData) +{ + ObservationsPerView observationsPerView; + for(auto& landIt : sfmData.getLandmarks()) + { + for(const auto& obsIt : landIt.second.observations) + { + IndexT viewId = obsIt.first; + auto& landmarksSet = observationsPerView[viewId]; + landmarksSet.first.push_back(obsIt.second); + landmarksSet.second.push_back(&landIt.second); + } + } + return observationsPerView; +} + bool filterLandmarks(SfMData& sfmData, double radiusScale, bool useFeatureScale, int minNbObservationsPerLandmark) { const auto initialNbLandmarks = sfmData.getLandmarks().size(); @@ -648,7 +672,7 @@ bool filterObservations2D(SfMData& sfmData, int nbNeighbors2D, float percentile, std::advance(itView, i); const IndexT viewId = *itView; - auto& observationsIt = observationsPerView.find(viewId); + auto observationsIt = observationsPerView.find(viewId); if(observationsIt == observationsPerView.end()) continue; @@ -669,7 +693,7 @@ bool filterObservations2D(SfMData& sfmData, int nbNeighbors2D, float percentile, cacheSize); for(auto j = 0; j < observations.size(); j++) { - const auto& obs = *observations[j]; + const auto& obs = observations[j]; std::vector indices_(nbNeighbors_); std::vector distances_(nbNeighbors_); tree.knnSearch(obs.x.data(), nbNeighbors_, &indices_[0], &distances_[0]); @@ -694,7 +718,7 @@ bool filterObservations2D(SfMData& sfmData, int nbNeighbors2D, float percentile, } estimatedRadii_[i] = radius; - std::vector filteredObservations; + std::vector filteredObservations; std::vector filteredLandmarks; filteredObservations.reserve(observations.size()); filteredLandmarks.reserve(landmarks.size()); @@ -724,14 +748,14 @@ bool filterObservations2D(SfMData& sfmData, int nbNeighbors2D, float percentile, if(estimatedRadii_[i] != -1.) estimatedRadii[viewId] = estimatedRadii_[i]; - auto& observationsIt = observationsPerView.find(viewId); + auto observationsIt = observationsPerView.find(viewId); if(observationsIt != observationsPerView.end()) { auto& observations = observationsIt->second.first; auto& landmarks = observationsIt->second.second; for(int j = 0; j < observations.size(); j++) { - landmarks[j]->observations[viewId] = *observations[j]; + landmarks[j]->observations[viewId] = observations[j]; } } } @@ -754,7 +778,7 @@ int aliceVision_main(int argc, char *argv[]) double neighborsInfluence = 0.5; int nbIterations = 5; int nbNeighbors2D = 5; - float percentile = 0.95; + float percentile = 0.95f; // user optional parameters std::vector featuresFolders; diff --git a/src/software/pipeline/main_prepareDenseScene.cpp b/src/software/pipeline/main_prepareDenseScene.cpp index 9170b014fb..2e1daed7e2 100644 --- a/src/software/pipeline/main_prepareDenseScene.cpp +++ b/src/software/pipeline/main_prepareDenseScene.cpp @@ -83,7 +83,29 @@ void process(const std::string& dstColorImage, } } -bool prepareDenseScene(SfMData& sfmData, +using ObservationsPerView = stl::flat_map>; + +/** + * @brief Get the landmark observations of camera views. + * @param[in] sfmData: A given SfMData + * @return Observations per camera view + */ +ObservationsPerView getObservationsPerViews(const SfMData& sfmData) +{ + ObservationsPerView observationsPerView; + for(auto& landIt : sfmData.getLandmarks()) + { + for(const auto& obsIt : landIt.second.observations) + { + IndexT viewId = obsIt.first; + auto& observationsSet = observationsPerView[viewId]; + observationsSet.push_back(&obsIt.second); + } + } + return observationsPerView; +} + +bool prepareDenseScene(const SfMData& sfmData, const std::vector& imagesFolders, const std::vector& masksFolders, const std::string& maskExtension, @@ -132,7 +154,7 @@ bool prepareDenseScene(SfMData& sfmData, ALICEVISION_LOG_INFO("Median Camera Exposure: " << medianCameraExposure << ", Median EV: " << std::log2(1.0 / medianCameraExposure)); bool doMaskLandmarks = landmarksMaskScale > 0.f; - sfmData::ObservationsPerView observationsPerView; + ObservationsPerView observationsPerView; HashMap estimatedRadii; if (doMaskLandmarks) { @@ -307,7 +329,7 @@ bool prepareDenseScene(SfMData& sfmData, const auto& observationsIt = observationsPerView.find(viewId); if(observationsIt != observationsPerView.end()) { - const auto& observations = observationsIt->second.first; + const auto& observations = observationsIt->second; int j = 0; for(const auto& observation : observations) { From 1483c7033a21ccb21a5bb45e4f52ec45269b04e5 Mon Sep 17 00:00:00 2001 From: Abdelrahman AL MAROUK Date: Mon, 25 Sep 2023 11:26:09 +0200 Subject: [PATCH 26/37] [filterSfM] improve observation selection strategy - add a re-distribution of observation scores at each iteration by favouring best N observations that will be selected at the end - fix hard-coded parameters --- src/software/pipeline/main_filterSfM.cpp | 47 ++++++++++++++++-------- 1 file changed, 32 insertions(+), 15 deletions(-) diff --git a/src/software/pipeline/main_filterSfM.cpp b/src/software/pipeline/main_filterSfM.cpp index ea695f6c35..f16f40adaa 100644 --- a/src/software/pipeline/main_filterSfM.cpp +++ b/src/software/pipeline/main_filterSfM.cpp @@ -215,7 +215,8 @@ ObservationsPerView getObservationsPerViews(SfMData& sfmData) return observationsPerView; } -bool filterLandmarks(SfMData& sfmData, double radiusScale, bool useFeatureScale, int minNbObservationsPerLandmark) +bool filterLandmarks(SfMData& sfmData, double radiusScale, bool useFeatureScale, int minNbObservationsPerLandmark, + int nbNeighbors3D, int maxNbObservationsPerLandmark) { const auto initialNbLandmarks = sfmData.getLandmarks().size(); std::vector landmarksData(initialNbLandmarks); @@ -255,11 +256,8 @@ bool filterLandmarks(SfMData& sfmData, double radiusScale, bool useFeatureScale, tree.buildIndex(); ALICEVISION_LOG_INFO("KdTree created for " << landmarksData.size() << " points."); - // TODO as parameter - int nbNeighbors3D = 10; // note that the landmark is a neighbor to itself with zero distance, hence the +/- 1 int nbNeighbors_ = std::min(nbNeighbors3D, static_cast(landmarksData.size() - 1)) + 1; - // contains the observing view ids and neighbors for each landmark std::vector, std::vector>> viewData(landmarksData.size()); @@ -287,8 +285,6 @@ bool filterLandmarks(SfMData& sfmData, double radiusScale, bool useFeatureScale, std::vector toRemove(landmarksData.size(), false); size_t nbToRemove = 0; - // TODO as parameter - int maxNbObservationsPerLandmark = 2; const auto initialNbLandmarks = sfmData.getLandmarks().size(); #pragma omp parallel for reduction(+ : nbToRemove) for(auto i = 0; i < landmarksData.size(); i++) @@ -552,7 +548,6 @@ bool filterObservations3D(SfMData& sfmData, int maxNbObservationsPerLandmark, in std::vector> viewScoresData_t(landmarksData.size()); for(auto i = 0; i < nbIterations; i++) { - ALICEVISION_LOG_INFO("Iteration " << i << "..."); #pragma omp parallel for for(auto id = 0; id < landmarksData.size(); id++) { @@ -603,23 +598,44 @@ bool filterObservations3D(SfMData& sfmData, int maxNbObservationsPerLandmark, in { // normalize score and apply influence factor viewScores_acc[j] *= neighborsInfluence / viewScores_total; - // combine weghted neighbor scores and the landmark's own scores + // combine weighted neighbor scores and the landmark's own scores viewScores_acc[j] += (1 - neighborsInfluence) * viewScores[j]; } + // dampen scores of non-chosen observations + if(viewScores_acc.size() <= maxNbObservationsPerLandmark) + continue; + // sort by descending view score order + std::vector idx(viewScores_acc.size()); + std::iota(idx.begin(), idx.end(), 0); + std::stable_sort(idx.begin(), idx.end(), + [&v = viewScores_acc](size_t i1, size_t i2) { return v[i1] > v[i2]; }); + viewScores_total = 1.; + for(auto j = maxNbObservationsPerLandmark; j < viewScores_acc.size(); j++) + { + const double& v = 0.5 * viewScores_acc[j]; + viewScores_acc[j] = v; + viewScores_total -= v; + } + // re-normalize + for(auto j = 0; j < viewScores_acc.size(); j++) + viewScores_acc[j] /= viewScores_total; } // update scores at end of iteration double error = 0.; #pragma omp parallel for reduction(+:error) for(auto id = 0; id < landmarksData.size(); id++) { - double error_j = 0.; - for(auto j = 0; j < viewScoresData_t[id].size(); j++) + // compute MSE { - const auto& v = viewScoresData_t[id][j] - viewScoresData[id].second[j]; - error_j += v * v; + double error_j = 0.; + for(auto j = 0; j < viewScoresData_t[id].size(); j++) + { + const auto& v = viewScoresData_t[id][j] - viewScoresData[id].second[j]; + error_j += v * v; + } + error_j /= viewScoresData_t[id].size(); + error += error_j; } - error_j /= viewScoresData_t[id].size(); - error += error_j; viewScoresData[id].second = std::move(viewScoresData_t[id]); } error /= landmarksData.size(); @@ -850,7 +866,8 @@ int aliceVision_main(int argc, char *argv[]) if(radiusScale > 0 || minNbObservationsPerLandmark > 0) { ALICEVISION_LOG_INFO("Filtering landmarks: started."); - filterLandmarks(sfmData, radiusScale, useFeatureScale, minNbObservationsPerLandmark); + filterLandmarks(sfmData, radiusScale, useFeatureScale, minNbObservationsPerLandmark, nbNeighbors3D, + maxNbObservationsPerLandmark); ALICEVISION_LOG_INFO("Filtering landmarks: done."); } From 7f344cedbac253651c759058ced8910fcf3dc007 Mon Sep 17 00:00:00 2001 From: Abdelrahman AL MAROUK Date: Tue, 26 Sep 2023 17:10:24 +0200 Subject: [PATCH 27/37] [filterSfM] multiple improvements - update parameters and code structure for a clearer understanding of different steps - fix bug with damping renormalization - add comments --- src/software/pipeline/main_filterSfM.cpp | 635 +++++++++++++---------- 1 file changed, 372 insertions(+), 263 deletions(-) diff --git a/src/software/pipeline/main_filterSfM.cpp b/src/software/pipeline/main_filterSfM.cpp index f16f40adaa..862a7b1759 100644 --- a/src/software/pipeline/main_filterSfM.cpp +++ b/src/software/pipeline/main_filterSfM.cpp @@ -48,6 +48,89 @@ using namespace boost::accumulators; static const std::size_t MAX_LEAF_ELEMENTS = 64; +struct FilterParams +{ + struct FilterLandmarksParams + { + bool enabled = true; + struct FilterLandmarksStep1Params + { + bool enabled = true; + int minNbObservationsPerLandmark = 3; + } step1; + struct FilterLandmarksStep2Params + { + bool enabled = true; + int maxNbObservationsPerLandmark = 2; + int nbNeighbors3D = 10; + } step2; + struct FilterLandmarksStep3Params + { + bool enabled = true; + double radiusScale = 2; + bool useFeatureScale = true; + } step3; + } filterLandmarks; + struct FilterObservations3DParams + { + bool enabled = true; + int maxNbObservationsPerLandmark = 2; + bool propagationEnabled = true; + int nbNeighbors3D = 10; + double neighborsInfluence = 0.5; + int nbIterations = 5; + bool dampingEnabled = true; + double dampingFactor = 0.5; + } filterObservations3D; + struct FilterObservations2DParams + { + bool enabled = true; + int nbNeighbors2D = 5; + float percentile = 0.95f; + double maskRadiusThreshold = 0.1; + } filterObservations2D; +}; + +std::istream& operator>>(std::istream& in, FilterParams& params) +{ + std::string token; + in >> token; + std::vector splitParams; + boost::split(splitParams, token, boost::algorithm::is_any_of(":")); + if(splitParams.size() != 21) + throw std::invalid_argument("Failed to parse FilterParams from: \n" + token); + int i = 0; + + params.filterLandmarks.enabled = boost::to_lower_copy(splitParams[i++]) == "true"; + + params.filterLandmarks.step1.enabled = boost::to_lower_copy(splitParams[i++]) == "true"; + params.filterLandmarks.step1.minNbObservationsPerLandmark = boost::lexical_cast(splitParams[i++]); + + params.filterLandmarks.step2.enabled = boost::to_lower_copy(splitParams[i++]) == "true"; + params.filterLandmarks.step2.maxNbObservationsPerLandmark = boost::lexical_cast(splitParams[i++]); + params.filterLandmarks.step2.nbNeighbors3D = boost::lexical_cast(splitParams[i++]); + + params.filterLandmarks.step3.enabled = boost::to_lower_copy(splitParams[i++]) == "true"; + params.filterLandmarks.step3.radiusScale = boost::lexical_cast(splitParams[i++]); + params.filterLandmarks.step3.useFeatureScale = boost::to_lower_copy(splitParams[i++]) == "true"; + + params.filterObservations3D.enabled = boost::to_lower_copy(splitParams[i++]) == "true"; + params.filterObservations3D.maxNbObservationsPerLandmark = boost::lexical_cast(splitParams[i++]); + params.filterObservations3D.propagationEnabled = boost::to_lower_copy(splitParams[i++]) == "true"; + params.filterObservations3D.nbNeighbors3D = boost::lexical_cast(splitParams[i++]); + params.filterObservations3D.neighborsInfluence = boost::lexical_cast(splitParams[i++]); + params.filterObservations3D.nbIterations = boost::lexical_cast(splitParams[i++]); + params.filterObservations3D.dampingEnabled = boost::to_lower_copy(splitParams[i++]) == "true"; + params.filterObservations3D.dampingFactor = boost::lexical_cast(splitParams[i++]); + + params.filterObservations2D.enabled = boost::to_lower_copy(splitParams[i++]) == "true"; + params.filterObservations2D.nbNeighbors2D = boost::lexical_cast(splitParams[i++]); + params.filterObservations2D.percentile = boost::lexical_cast(splitParams[i++]); + params.filterObservations2D.maskRadiusThreshold = boost::lexical_cast(splitParams[i++]); + + return in; +} + struct ObservationsAdaptator { using Derived = ObservationsAdaptator; //!< In this case the dataset class is myself. @@ -215,19 +298,20 @@ ObservationsPerView getObservationsPerViews(SfMData& sfmData) return observationsPerView; } -bool filterLandmarks(SfMData& sfmData, double radiusScale, bool useFeatureScale, int minNbObservationsPerLandmark, - int nbNeighbors3D, int maxNbObservationsPerLandmark) +bool filterLandmarks(SfMData& sfmData, const FilterParams::FilterLandmarksParams& params) { - const auto initialNbLandmarks = sfmData.getLandmarks().size(); - std::vector landmarksData(initialNbLandmarks); + std::vector landmarksData; + // Step 1 { + const auto& initialNbLandmarks = sfmData.getLandmarks().size(); + landmarksData.resize(initialNbLandmarks); size_t i = 0; - if (minNbObservationsPerLandmark > 0) + if(params.step1.enabled) { ALICEVISION_LOG_INFO("Removing landmarks having an insufficient number of observations: started."); for(auto& it : sfmData.getLandmarks()) { - if(it.second.observations.size() < minNbObservationsPerLandmark) + if(it.second.observations.size() < params.step1.minNbObservationsPerLandmark) continue; landmarksData[i++] = it.second; } @@ -249,18 +333,22 @@ bool filterLandmarks(SfMData& sfmData, double radiusScale, bool useFeatureScale, } } + // Step 2 + if(params.step2.enabled) { + ALICEVISION_LOG_INFO("Removing landmarks with dissimilar observations of 3D landmark neighbors: started."); + ALICEVISION_LOG_INFO("Build nanoflann KdTree index for landmarks in 3D."); LandmarksAdaptator data(landmarksData); KdTree3D tree(3, data, nanoflann::KDTreeSingleIndexAdaptorParams(MAX_LEAF_ELEMENTS)); tree.buildIndex(); ALICEVISION_LOG_INFO("KdTree created for " << landmarksData.size() << " points."); + ALICEVISION_LOG_INFO("Computing landmarks neighbors: started."); // note that the landmark is a neighbor to itself with zero distance, hence the +/- 1 - int nbNeighbors_ = std::min(nbNeighbors3D, static_cast(landmarksData.size() - 1)) + 1; + int nbNeighbors_ = std::min(params.step2.nbNeighbors3D, static_cast(landmarksData.size() - 1)) + 1; // contains the observing view ids and neighbors for each landmark std::vector, std::vector>> viewData(landmarksData.size()); - #pragma omp parallel for for(auto i = 0; i < landmarksData.size(); i++) { @@ -282,10 +370,12 @@ bool filterLandmarks(SfMData& sfmData, double radiusScale, bool useFeatureScale, // a landmark is a neighbor to itself with zero distance, remove it neighbors.erase(neighbors.begin()); } + ALICEVISION_LOG_INFO("Computing landmarks neighbors: done."); + ALICEVISION_LOG_INFO("Identifying landmarks to remove: started."); std::vector toRemove(landmarksData.size(), false); size_t nbToRemove = 0; - const auto initialNbLandmarks = sfmData.getLandmarks().size(); + const auto& initialNbLandmarks = landmarksData.size(); #pragma omp parallel for reduction(+ : nbToRemove) for(auto i = 0; i < landmarksData.size(); i++) { @@ -311,7 +401,7 @@ bool filterLandmarks(SfMData& sfmData, double radiusScale, bool useFeatureScale, if(!(*viewIds_neighbor_it < *viewIds_it)) { scoreNeighbor++; - if(scoreNeighbor == maxNbObservationsPerLandmark) + if(scoreNeighbor == params.step2.maxNbObservationsPerLandmark) { score++; break; @@ -322,7 +412,7 @@ bool filterLandmarks(SfMData& sfmData, double radiusScale, bool useFeatureScale, } } } - if(score < nbNeighbors3D / 2) + if(score < params.step2.nbNeighbors3D / 2) { toRemove[i] = true; nbToRemove++; @@ -332,112 +422,104 @@ bool filterLandmarks(SfMData& sfmData, double radiusScale, bool useFeatureScale, << ", i.e. " << (nbToRemove * 100.f / initialNbLandmarks) << " % "); ALICEVISION_LOG_INFO("Identifying landmarks to remove: done."); - ALICEVISION_LOG_INFO("Removing landmarks based on TODO: started."); - std::vector> filteredLandmarks(landmarksData.size() - nbToRemove); + ALICEVISION_LOG_INFO("Removing identified landmarks: started."); std::vector landmarksData_filtered(landmarksData.size() - nbToRemove); IndexT newIdx = 0; for(size_t i = 0; i < landmarksData.size(); i++) { if(!toRemove[i]) { - landmarksData_filtered[newIdx] = landmarksData[i]; - filteredLandmarks[newIdx++] = std::make_pair(newIdx, landmarksData[i]); + landmarksData_filtered[newIdx++] = landmarksData[i]; } } - sfmData.getLandmarks() = std::move(Landmarks(filteredLandmarks.begin(), filteredLandmarks.end())); landmarksData = std::move(landmarksData_filtered); - ALICEVISION_LOG_INFO("Removing landmarks based on TODO: done."); - } + ALICEVISION_LOG_INFO("Removing identified landmarks: done."); - if (radiusScale <= 0) - { - std::vector> filteredLandmarks(landmarksData.size()); - for(IndexT newIdx = 0; newIdx < landmarksData.size(); newIdx++) - { - filteredLandmarks[newIdx] = std::make_pair(newIdx, landmarksData[newIdx]); - } - sfmData.getLandmarks() = std::move(Landmarks(filteredLandmarks.begin(), filteredLandmarks.end())); - return true; + ALICEVISION_LOG_INFO("Removing landmarks with dissimilar observations of 3D landmark neighbors: done."); } - mvsUtils::MultiViewParams mp(sfmData, "", "", "", false); - std::vector landmarksPixSize(landmarksData.size()); - - ALICEVISION_LOG_INFO("Computing pixel size: started."); - #pragma omp parallel for - for(auto i = 0; i < landmarksData.size(); i++) + // Step 3 + if(params.step3.enabled) { - const Landmark& landmark = landmarksData[i]; + ALICEVISION_LOG_INFO("Removing landmarks with worse resolution than neighbors: started."); - // compute landmark pixSize - double pixSize = 0.; - int n = 0; - for(const auto& observationPair : landmark.observations) + mvsUtils::MultiViewParams mp(sfmData, "", "", "", false); + std::vector landmarksPixSize(landmarksData.size()); + + ALICEVISION_LOG_INFO("Computing pixel size: started."); +#pragma omp parallel for + for(auto i = 0; i < landmarksData.size(); i++) { - const IndexT viewId = observationPair.first; - pixSize += mp.getCamPixelSize( - Point3d(landmark.X.x(), landmark.X.y(), landmark.X.z()), - mp.getIndexFromViewId(viewId), - useFeatureScale ? observationPair.second.scale : 1); - n++; + const Landmark& landmark = landmarksData[i]; + // compute landmark pixSize + double pixSize = 0.; + int n = 0; + for(const auto& observationPair : landmark.observations) + { + const IndexT viewId = observationPair.first; + pixSize += mp.getCamPixelSize(Point3d(landmark.X.x(), landmark.X.y(), landmark.X.z()), + mp.getIndexFromViewId(viewId), + params.step3.useFeatureScale ? observationPair.second.scale : 1); + n++; + } + pixSize /= n; + landmarksPixSize[i] = pixSize; } - pixSize /= n; - - landmarksPixSize[i] = pixSize; - } - ALICEVISION_LOG_INFO("Computing pixel size: done."); + ALICEVISION_LOG_INFO("Computing pixel size: done."); - //// sort landmarks by descending pixSize order - //std::stable_sort(landmarksPixSize.begin(), landmarksPixSize.end(), std::greater<>{}); + ALICEVISION_LOG_INFO("Identifying landmarks to remove based on pixel size: started."); - ALICEVISION_LOG_INFO("Build nanoflann KdTree index for landmarks in 3D."); - LandmarksAdaptator data(landmarksData); - KdTree3D tree(3, data, nanoflann::KDTreeSingleIndexAdaptorParams(MAX_LEAF_ELEMENTS)); - tree.buildIndex(); - ALICEVISION_LOG_INFO("KdTree created for " << landmarksData.size() << " points."); - std::vector toRemove(landmarksData.size(), false); + ALICEVISION_LOG_INFO("Build nanoflann KdTree index for landmarks in 3D."); + LandmarksAdaptator data(landmarksData); + KdTree3D tree(3, data, nanoflann::KDTreeSingleIndexAdaptorParams(MAX_LEAF_ELEMENTS)); + tree.buildIndex(); + ALICEVISION_LOG_INFO("KdTree created for " << landmarksData.size() << " points."); - ALICEVISION_LOG_INFO("Identifying landmarks to remove based on pixel size: started."); + const auto& initialNbLandmarks = landmarksData.size(); + std::vector toRemove(landmarksData.size(), false); + size_t nbToRemove = 0; +#pragma omp parallel for reduction(+ : nbToRemove) + for(auto i = 0; i < landmarksData.size(); i++) + { + PixSizeSearch search(landmarksPixSize[i] * params.step3.radiusScale, landmarksPixSize, i); + bool found = tree.findNeighbors(search, landmarksData[i].X.data(), nanoflann::SearchParameters()); + if(found) + { + toRemove[i] = true; + nbToRemove++; + } + } + ALICEVISION_LOG_INFO("Identified " << nbToRemove << " landmarks to remove out of " << initialNbLandmarks + << ", i.e. " << (nbToRemove * 100.f / initialNbLandmarks) << " % "); + ALICEVISION_LOG_INFO("Identifying landmarks to remove: done."); - size_t nbToRemove = 0; - #pragma omp parallel for reduction(+:nbToRemove) - for(auto i = 0; i < landmarksData.size(); i++) - { - PixSizeSearch search(landmarksPixSize[i] * radiusScale, landmarksPixSize, i); - bool found = tree.findNeighbors(search, landmarksData[i].X.data(), nanoflann::SearchParameters()); - if (found) + ALICEVISION_LOG_INFO("Removing identified landmarks: started."); + std::vector landmarksData_filtered(landmarksData.size() - nbToRemove); + IndexT newIdx = 0; + for(size_t i = 0; i < landmarksData.size(); i++) { - toRemove[i] = true; - nbToRemove++; + if(!toRemove[i]) + { + landmarksData_filtered[newIdx++] = landmarksData[i]; + } } - } + landmarksData = std::move(landmarksData_filtered); + ALICEVISION_LOG_INFO("Removing identified landmarks: done."); - ALICEVISION_LOG_INFO( - "Identified " << nbToRemove << - " landmarks to remove out of " << initialNbLandmarks << - ", i.e. " << (nbToRemove * 100.f / initialNbLandmarks) << - " % " - ); - ALICEVISION_LOG_INFO("Identifying landmarks to remove: done."); + ALICEVISION_LOG_INFO("Removing landmarks with worse resolution than neighbors: done."); + } - ALICEVISION_LOG_INFO("Removing landmarks based on pixel size: started."); - std::vector> filteredLandmarks(landmarksData.size() - nbToRemove); - IndexT newIdx = 0; - for (size_t i = 0; i < landmarksData.size(); i++) + // applying modifications to Sfm data + std::vector> filteredLandmarks(landmarksData.size()); + for(IndexT newIdx = 0; newIdx < landmarksData.size(); newIdx++) { - if(!toRemove[i]) - { - filteredLandmarks[newIdx++] = std::make_pair(newIdx, landmarksData[i]); - } + filteredLandmarks[newIdx] = std::make_pair(newIdx, landmarksData[newIdx]); } sfmData.getLandmarks() = std::move(Landmarks(filteredLandmarks.begin(), filteredLandmarks.end())); - ALICEVISION_LOG_INFO("Removing landmarks based on pixel size: done."); - return true; } -bool filterObservations3D(SfMData& sfmData, int maxNbObservationsPerLandmark, int nbNeighbors3D, double neighborsInfluence, - int nbIterations) +bool filterObservations3D(SfMData& sfmData, const FilterParams::FilterObservations3DParams& params) { // store in vector for faster access std::vector landmarksData(sfmData.getLandmarks().size()); @@ -502,146 +584,152 @@ bool filterObservations3D(SfMData& sfmData, int maxNbObservationsPerLandmark, in } ALICEVISION_LOG_INFO("Computing initial observation scores based on distance to observing view: done"); - ALICEVISION_LOG_INFO("Computing landmark neighbors and distance-based weights: started"); - ALICEVISION_LOG_INFO("Build nanoflann KdTree index for landmarks in 3D."); - LandmarksAdaptator dataAdaptor(landmarksData); - KdTree3D tree(3, dataAdaptor, nanoflann::KDTreeSingleIndexAdaptorParams(MAX_LEAF_ELEMENTS)); - tree.buildIndex(); - ALICEVISION_LOG_INFO("KdTree created for " << landmarksData.size() << " points."); - // note that the landmark is a neighbor to itself with zero distance, hence the +/- 1 - int nbNeighbors_ = std::min(nbNeighbors3D, static_cast(landmarksData.size() - 1)) + 1; - // contains the neighbor landmarks ids with their corresponding weights - std::vector, std::vector>> neighborsData(landmarksData.size()); -#pragma omp parallel for - for(auto i = 0; i < landmarksData.size(); i++) + if(params.propagationEnabled) { - const sfmData::Landmark& landmark = *landmarksData[i]; - auto& [indices_, weights_] = neighborsData[i]; - indices_.resize(nbNeighbors_); - weights_.resize(nbNeighbors_); - tree.knnSearch(landmark.X.data(), nbNeighbors_, &indices_[0], &weights_[0]); - // a landmark is a neighbor to itself with zero distance, remove it - indices_.erase(indices_.begin()); - weights_.erase(weights_.begin()); - // accumulator used for normalisation - double total = 0.; - for(auto& w : weights_) - { - // weight is the inverse of distance between a landmark and its neighbor - w = 1. / std::sqrt(w); - if(std::isinf(w)) - w = std::numeric_limits::max(); - total += w; - } - if(std::isinf(total)) - total = std::numeric_limits::max(); - // normalize weights - for(auto& w : weights_) + ALICEVISION_LOG_INFO("Computing landmark neighbors and distance-based weights: started"); + ALICEVISION_LOG_INFO("Build nanoflann KdTree index for landmarks in 3D."); + LandmarksAdaptator dataAdaptor(landmarksData); + KdTree3D tree(3, dataAdaptor, nanoflann::KDTreeSingleIndexAdaptorParams(MAX_LEAF_ELEMENTS)); + tree.buildIndex(); + ALICEVISION_LOG_INFO("KdTree created for " << landmarksData.size() << " points."); + // note that the landmark is a neighbor to itself with zero distance, hence the +/- 1 + int nbNeighbors_ = std::min(params.nbNeighbors3D, static_cast(landmarksData.size() - 1)) + 1; + // contains the neighbor landmarks ids with their corresponding weights + std::vector, std::vector>> neighborsData(landmarksData.size()); +#pragma omp parallel for + for(auto i = 0; i < landmarksData.size(); i++) { - w /= total; + const sfmData::Landmark& landmark = *landmarksData[i]; + auto& [indices_, weights_] = neighborsData[i]; + indices_.resize(nbNeighbors_); + weights_.resize(nbNeighbors_); + tree.knnSearch(landmark.X.data(), nbNeighbors_, &indices_[0], &weights_[0]); + // a landmark is a neighbor to itself with zero distance, remove it + indices_.erase(indices_.begin()); + weights_.erase(weights_.begin()); + // accumulator used for normalisation + double total = 0.; + for(auto& w : weights_) + { + // weight is the inverse of distance between a landmark and its neighbor + w = 1. / std::sqrt(w); + if(std::isinf(w)) + w = std::numeric_limits::max(); + total += w; + } + if(std::isinf(total)) + total = std::numeric_limits::max(); + // normalize weights + for(auto& w : weights_) + { + w /= total; + } } - } - ALICEVISION_LOG_INFO("Computing landmark neighbors and distance-based weights: done"); + ALICEVISION_LOG_INFO("Computing landmark neighbors and distance-based weights: done"); - ALICEVISION_LOG_INFO("Propagating neighbors observation scores: started"); - // new view scores at iteration t - std::vector> viewScoresData_t(landmarksData.size()); - for(auto i = 0; i < nbIterations; i++) - { -#pragma omp parallel for - for(auto id = 0; id < landmarksData.size(); id++) + ALICEVISION_LOG_INFO("Propagating neighbors observation scores: started"); + // new view scores at iteration t + std::vector> viewScoresData_t(landmarksData.size()); + for(auto i = 0; i < params.nbIterations; i++) { - const auto& [viewIds, viewScores] = viewScoresData[id]; - auto& viewScores_acc = viewScoresData_t[id]; - // initialize to zero, will first contain the weighted average scores of neighbors - viewScores_acc.assign(viewScores.size(), 0.); - // accumulator for normalisation - double viewScores_total = 0.; - auto& [indices_, weights_] = neighborsData[id]; - for(auto j = 0; j < nbNeighbors3D; j++) +#pragma omp parallel for + for(auto id = 0; id < landmarksData.size(); id++) { - const auto& neighborId = indices_[j]; - const auto& neighborWeight = weights_[j]; - const auto& [viewIds_neighbor, viewScores_neighbor] = viewScoresData[neighborId]; - // loop over common views - auto viewIds_it = viewIds.begin(); - auto viewIds_neighbor_it = viewIds_neighbor.begin(); - auto viewScores_neighbor_it = viewScores_neighbor.begin(); - auto viewScores_acc_it = viewScores_acc.begin(); - while(viewIds_it != viewIds.end() && viewIds_neighbor_it != viewIds_neighbor.end()) + const auto& [viewIds, viewScores] = viewScoresData[id]; + auto& viewScores_acc = viewScoresData_t[id]; + // initialize to zero, will first contain the weighted average scores of neighbors + viewScores_acc.assign(viewScores.size(), 0.); + // accumulator for normalisation + double viewScores_total = 0.; + auto& [indices_, weights_] = neighborsData[id]; + for(auto j = 0; j < params.nbNeighbors3D; j++) { - if(*viewIds_it < *viewIds_neighbor_it) - { - ++viewIds_it; - ++viewScores_acc_it; - } - else + const auto& neighborId = indices_[j]; + const auto& neighborWeight = weights_[j]; + const auto& [viewIds_neighbor, viewScores_neighbor] = viewScoresData[neighborId]; + // loop over common views + auto viewIds_it = viewIds.begin(); + auto viewIds_neighbor_it = viewIds_neighbor.begin(); + auto viewScores_neighbor_it = viewScores_neighbor.begin(); + auto viewScores_acc_it = viewScores_acc.begin(); + while(viewIds_it != viewIds.end() && viewIds_neighbor_it != viewIds_neighbor.end()) { - // if same view, accumulate weighted scores - if(!(*viewIds_neighbor_it < *viewIds_it)) + if(*viewIds_it < *viewIds_neighbor_it) { - const auto& v = *viewScores_neighbor_it * neighborWeight; - (*viewScores_acc_it) += v; - viewScores_total += v; ++viewIds_it; ++viewScores_acc_it; } - ++viewIds_neighbor_it; - ++viewScores_neighbor_it; + else + { + // if same view, accumulate weighted scores + if(!(*viewIds_neighbor_it < *viewIds_it)) + { + const auto& v = *viewScores_neighbor_it * neighborWeight; + (*viewScores_acc_it) += v; + viewScores_total += v; + ++viewIds_it; + ++viewScores_acc_it; + } + ++viewIds_neighbor_it; + ++viewScores_neighbor_it; + } } } + // if common views with neighbor landmarks + if(viewScores_total != 0.) + for(auto j = 0; j < viewScores_acc.size(); j++) + { + // normalize score and apply influence factor + viewScores_acc[j] *= params.neighborsInfluence / viewScores_total; + // combine weighted neighbor scores and the landmark's own scores + viewScores_acc[j] += (1 - params.neighborsInfluence) * viewScores[j]; + } + + // dampen scores of non-chosen observations + if(params.dampingEnabled && viewScores_acc.size() <= params.maxNbObservationsPerLandmark) + { + // sort by descending view score order + std::vector idx(viewScores_acc.size()); + std::iota(idx.begin(), idx.end(), 0); + std::stable_sort(idx.begin(), idx.end(), + [&v = viewScores_acc](size_t i1, size_t i2) { return v[i1] > v[i2]; }); + viewScores_total = 1.; + for(auto j = params.maxNbObservationsPerLandmark; j < viewScores_acc.size(); j++) + { + const double& v = params.dampingFactor * viewScores_acc[j]; + viewScores_total += v - viewScores_acc[j]; + viewScores_acc[j] = v; + } + // re-normalize + for(auto j = 0; j < viewScores_acc.size(); j++) + viewScores_acc[j] /= viewScores_total; + } } - // if no common views with neighbor landmarks - if(viewScores_total == 0.) - continue; - for(auto j = 0; j < viewScores_acc.size(); j++) - { - // normalize score and apply influence factor - viewScores_acc[j] *= neighborsInfluence / viewScores_total; - // combine weighted neighbor scores and the landmark's own scores - viewScores_acc[j] += (1 - neighborsInfluence) * viewScores[j]; - } - // dampen scores of non-chosen observations - if(viewScores_acc.size() <= maxNbObservationsPerLandmark) - continue; - // sort by descending view score order - std::vector idx(viewScores_acc.size()); - std::iota(idx.begin(), idx.end(), 0); - std::stable_sort(idx.begin(), idx.end(), - [&v = viewScores_acc](size_t i1, size_t i2) { return v[i1] > v[i2]; }); - viewScores_total = 1.; - for(auto j = maxNbObservationsPerLandmark; j < viewScores_acc.size(); j++) + // Mean Squared Error + double error = 0.; + // update scores at end of iteration +#pragma omp parallel for reduction(+ : error) + for(auto id = 0; id < landmarksData.size(); id++) { - const double& v = 0.5 * viewScores_acc[j]; - viewScores_acc[j] = v; - viewScores_total -= v; - } - // re-normalize - for(auto j = 0; j < viewScores_acc.size(); j++) - viewScores_acc[j] /= viewScores_total; - } - // update scores at end of iteration - double error = 0.; -#pragma omp parallel for reduction(+:error) - for(auto id = 0; id < landmarksData.size(); id++) - { - // compute MSE - { - double error_j = 0.; - for(auto j = 0; j < viewScoresData_t[id].size(); j++) + // compute MSE { - const auto& v = viewScoresData_t[id][j] - viewScoresData[id].second[j]; - error_j += v * v; + double error_j = 0.; + for(auto j = 0; j < viewScoresData_t[id].size(); j++) + { + const auto& v = viewScoresData_t[id][j] - viewScoresData[id].second[j]; + error_j += v * v; + } + error_j /= viewScoresData_t[id].size(); + error += error_j; } - error_j /= viewScoresData_t[id].size(); - error += error_j; + // update scores + viewScoresData[id].second = std::move(viewScoresData_t[id]); } - viewScoresData[id].second = std::move(viewScoresData_t[id]); + error /= landmarksData.size(); + ALICEVISION_LOG_INFO("MSE at iteration " << i << ": " << error); } - error /= landmarksData.size(); - ALICEVISION_LOG_INFO("MSE at iteration " << i << ": " << error); + ALICEVISION_LOG_INFO("Propagating neighbors observation scores: done"); } - ALICEVISION_LOG_INFO("Propagating neighbors observation scores: done"); ALICEVISION_LOG_INFO("Selecting observations with best scores: started"); #pragma omp parallel for @@ -652,7 +740,7 @@ bool filterObservations3D(SfMData& sfmData, int maxNbObservationsPerLandmark, in auto& [viewIds, viewScores] = viewScoresData[i]; // check number of observations - if(landmark.observations.size() <= maxNbObservationsPerLandmark) + if(landmark.observations.size() <= params.maxNbObservationsPerLandmark) continue; // sort by descending view score order @@ -662,7 +750,7 @@ bool filterObservations3D(SfMData& sfmData, int maxNbObservationsPerLandmark, in // keep only observations with best scores Observations filteredObservations; - for(auto j = 0; j < maxNbObservationsPerLandmark; j++) + for(auto j = 0; j < params.maxNbObservationsPerLandmark; j++) { const auto& viewId = viewIds[idx[j]]; filteredObservations[viewId] = landmark.observations[viewId]; @@ -674,7 +762,7 @@ bool filterObservations3D(SfMData& sfmData, int maxNbObservationsPerLandmark, in return true; } -bool filterObservations2D(SfMData& sfmData, int nbNeighbors2D, float percentile, +bool filterObservations2D(SfMData& sfmData, const FilterParams::FilterObservations2DParams& params, HashMap& estimatedRadii) { std::set viewIds = sfmData.getValidViews(); @@ -687,53 +775,57 @@ bool filterObservations2D(SfMData& sfmData, int nbNeighbors2D, float percentile, auto itView = viewIds.begin(); std::advance(itView, i); const IndexT viewId = *itView; - auto observationsIt = observationsPerView.find(viewId); if(observationsIt == observationsPerView.end()) continue; - auto& observations = observationsIt->second.first; auto& landmarks = observationsIt->second.second; - ALICEVISION_LOG_INFO("Build nanoflann KdTree index for projected landmarks in 2D."); + // Build nanoflann KdTree index for projected landmarks in 2D ObservationsAdaptator data(observations); KdTree2D tree(2, data, nanoflann::KDTreeSingleIndexAdaptorParams(MAX_LEAF_ELEMENTS)); tree.buildIndex(); - ALICEVISION_LOG_INFO("KdTree created for " << observations.size() << " points."); // note that the observation is a neighbor to itself with zero distance, hence the +/- 1 - size_t nbNeighbors_ = std::min(static_cast(nbNeighbors2D), observations.size() - 1) + 1; + size_t nbNeighbors_ = std::min(static_cast(params.nbNeighbors2D), observations.size() - 1) + 1; + // average neighbors distance for each observation std::vector means(observations.size()); const std::size_t cacheSize = 1000; - accumulator_set, tag::mean>> acc(tag::tail::cache_size = - cacheSize); + // accumulator for quantile computation + accumulator_set>> acc(tag::tail::cache_size = cacheSize); for(auto j = 0; j < observations.size(); j++) { + // Find neighbors and the corresponding distances const auto& obs = observations[j]; std::vector indices_(nbNeighbors_); std::vector distances_(nbNeighbors_); tree.knnSearch(obs.x.data(), nbNeighbors_, &indices_[0], &distances_[0]); - + // returned distances are L2 -> apply squared root std::transform(distances_.begin(), distances_.end(), distances_.begin(), static_cast(std::sqrt)); + // average distance const auto& mean = std::accumulate(distances_.begin(), distances_.end(), 0.0) / (nbNeighbors_ - 1); means[j] = mean; + // update accumulator acc(mean); } double mean_max = std::numeric_limits::max(); - if(percentile != 1.f) - mean_max = quantile(acc, quantile_probability = percentile); - - double radius = mean(acc); + // check to avoid exception + if(params.percentile != 1.f) + mean_max = quantile(acc, quantile_probability = params.percentile); + // estimated mask radius is the average of distance means + // quantile is used to avoid outlier bias + double radius = quantile(acc, quantile_probability = (1.f - params.percentile) * 0.5f); // check if estimated radius is too large { const View& view = *(sfmData.getViews().at(viewId)); - double radiusMax = 0.15 * 0.5 * (view.getImage().getWidth() + view.getImage().getHeight()); + double radiusMax = params.maskRadiusThreshold * 0.5 * (view.getImage().getWidth() + view.getImage().getHeight()); if(radius > radiusMax) radius = radiusMax; } estimatedRadii_[i] = radius; + // filter outlier observations std::vector filteredObservations; std::vector filteredLandmarks; filteredObservations.reserve(observations.size()); @@ -750,11 +842,11 @@ bool filterObservations2D(SfMData& sfmData, int nbNeighbors2D, float percentile, landmarks = std::move(filteredLandmarks); } + // clear and update landmark observations for(auto& landmark : sfmData.getLandmarks()) { landmark.second.observations.clear(); } - for(int i = 0; i < viewIds.size(); ++i) { auto itView = viewIds.begin(); @@ -764,7 +856,7 @@ bool filterObservations2D(SfMData& sfmData, int nbNeighbors2D, float percentile, if(estimatedRadii_[i] != -1.) estimatedRadii[viewId] = estimatedRadii_[i]; - auto observationsIt = observationsPerView.find(viewId); + const auto& observationsIt = observationsPerView.find(viewId); if(observationsIt != observationsPerView.end()) { auto& observations = observationsIt->second.first; @@ -775,6 +867,14 @@ bool filterObservations2D(SfMData& sfmData, int nbNeighbors2D, float percentile, } } } + // remove landmarks with no remaining observations + for(auto it = sfmData.getLandmarks().begin(); it != sfmData.getLandmarks().end();) + { + if(it->second.observations.size() == 0) + it = sfmData.getLandmarks().erase(it); + else + ++it; + } return true; } @@ -785,18 +885,11 @@ int aliceVision_main(int argc, char *argv[]) std::string inputSfmFilename; std::string outputSfmFilename; - std::string outputRadiiFilename; - int maxNbObservationsPerLandmark = 2; - int minNbObservationsPerLandmark = 5; - double radiusScale = 2; - bool useFeatureScale = true; - int nbNeighbors3D = 10; - double neighborsInfluence = 0.5; - int nbIterations = 5; - int nbNeighbors2D = 5; - float percentile = 0.95f; // user optional parameters + FilterParams params; + std::string outputRadiiFilename; + // required for 2D visualization in meshroom std::vector featuresFolders; std::vector matchesFolders; std::string describerTypesName = feature::EImageDescriberType_enumToString(feature::EImageDescriberType::SIFT); @@ -810,33 +903,49 @@ int aliceVision_main(int argc, char *argv[]) po::options_description optionalParams("Optional parameters"); optionalParams.add_options() - ("maxNbObservationsPerLandmark", po::value(&maxNbObservationsPerLandmark)->default_value(maxNbObservationsPerLandmark), - "Maximum number of allowed observations per landmark.") - ("minNbObservationsPerLandmark", po::value(&minNbObservationsPerLandmark)->default_value(minNbObservationsPerLandmark), - "Minimum number of observations required to keep a landmark.") - ("radiusScale", po::value(&radiusScale)->default_value(radiusScale), - "Scale factor applied to pixel size based radius filter applied to landmarks.") - ("useFeatureScale", po::value(&useFeatureScale)->default_value(useFeatureScale), - "If true, use feature scale for computing pixel size. Otherwise, use a scale of 1 pixel.") - ("nbNeighbors3D", po::value(&nbNeighbors3D)->default_value(nbNeighbors3D), - "Number of neighbor landmarks used in making the decision for best observations.") - ("neighborsInfluence", po::value(&neighborsInfluence)->default_value(neighborsInfluence), - "Specifies how much influential the neighbors are in selecting the best observations." - "Between 0. and 1., the closer to 1., the more influencial the neighborhood is.") - ("nbIterations", po::value(&nbIterations)->default_value(nbIterations), - "Number of iterations to propagate neighbors information.") - ("nbNeighbors2D", po::value(&nbNeighbors2D)->default_value(nbNeighbors2D), - "Number of neighbor observations to be considered for the landmarks-based masking.") - ("percentile", po::value(&percentile)->default_value(percentile), - "TODO.") + ("filterParams", po::value(¶ms)->default_value(params), + "Filter parameters devided into 3 successive stages.\n" + "\n" + "Filter landmarks parameters:\n" + " * Enabled: Filter Landmarks over multiple steps.\n" + " * Step 1 parameters:\n" + " * Enabled: Remove landmarks with insufficient observations.\n" + " * Min Nb of Observations: Minimum number of observations required to keep a landmark.\n" + " * Step 2 parameters:\n" + " * Enabled: Remove landmarks with dissimilar observations of 3D landmark neighbors.\n" + " * Max Nb of Observations: Maximum number of allowed observations per landmark.\n" + " * Nb of Neighbors in 3D: Number of neighbor landmarks used in making the decision for best observations.\n" + " * Step 3 parameters:\n" + " * Enabled: Remove landmarks with worse resolution than neighbors.\n" + " * Radius Scale: Scale factor applied to pixel size based radius filter applied to landmarks.\n" + " * Use Feature Scale: If true, use feature scale for computing pixel size. Otherwise, use a scale of 1 pixel.\n" + "\n" + "Filter Observations 3D parameters:\n" + " * Enabled: Select best observations for observation consistency between 3D neighboring landmarks.\n" + " * Max Nb of Observations: Maximum number of allowed observations per landmark.\n" + " * Enable Neighbors Influence: nable propagating neighbors' scores iteratively.\n" + " * Nb of Neighbors in 3D: Number of neighbor landmarks used in making the decision for best observations.\n" + " * Neighbors Influence: Specifies how much influential the neighbors are in selecting the best observations.\n" + " Between 0. and 1., the closer to 1., the more influencial the neighborhood is.\n" + " * Nb of Iterations: Number of iterations to propagate neighbors information.\n" + " * Enable Damping: Enable additional damping of observations to reject after each iterations.\n" + " * Damping Factor: Multiplicative damping factor.\n" + "\n" + "Filter Observations 2D parameters:\n" + " * Enabled: Select best observations for observation consistency between 2D projected neighboring\n" + " landmarks per view. Eventually remove landmarks with no remaining observations.\n" + " Also estimate depth map mask radius per view based on landmarks.\n" + " * Nb of Neighbors in 2D: Number of neighbor observations to be considered for the landmarks-based masking.\n" + " * Percentile: Used as a quantile probability for filtering relatively outlier observations.\n" + " * Mask Radius Threshold: Percentage of image size to be used as an upper limit for estimated mask radius.") + ("outputRadiiFile", po::value(&outputRadiiFilename)->default_value(outputRadiiFilename), + "Output Radii file containing the estimated projection radius of observations per view.") ("featuresFolders,f", po::value>(&featuresFolders)->multitoken(), "Path to folder(s) containing the extracted features.") ("matchesFolders,m", po::value>(&matchesFolders)->multitoken(), "Path to folder(s) in which computed matches are stored.") ("describerTypes,d", po::value(&describerTypesName)->default_value(describerTypesName), - feature::EImageDescriberType_informations().c_str()) - ("outputRadiiFile", po::value(&outputRadiiFilename)->default_value(outputRadiiFilename), - "Output Radii file containing the estimated projection radius of observations per view."); + feature::EImageDescriberType_informations().c_str()); CmdLine cmdline("AliceVision SfM filtering."); // TODO add description cmdline.add(requiredParams); @@ -863,26 +972,26 @@ int aliceVision_main(int argc, char *argv[]) // filter SfM data - if(radiusScale > 0 || minNbObservationsPerLandmark > 0) + if(params.filterLandmarks.enabled && (params.filterLandmarks.step1.enabled || + params.filterLandmarks.step2.enabled || params.filterLandmarks.step3.enabled)) { ALICEVISION_LOG_INFO("Filtering landmarks: started."); - filterLandmarks(sfmData, radiusScale, useFeatureScale, minNbObservationsPerLandmark, nbNeighbors3D, - maxNbObservationsPerLandmark); + filterLandmarks(sfmData, params.filterLandmarks); ALICEVISION_LOG_INFO("Filtering landmarks: done."); } - if(maxNbObservationsPerLandmark > 0) + if(params.filterObservations3D.enabled) { ALICEVISION_LOG_INFO("Filtering observations in 3D: started."); - filterObservations3D(sfmData, maxNbObservationsPerLandmark, nbNeighbors3D, neighborsInfluence, nbIterations); + filterObservations3D(sfmData, params.filterObservations3D); ALICEVISION_LOG_INFO("Filtering observations in 3D: done."); } - if(nbNeighbors2D > 0 || percentile < 1.f) + if(params.filterObservations2D.enabled) { HashMap estimatedRadii; ALICEVISION_LOG_INFO("Filtering observations in 2D: started."); - filterObservations2D(sfmData, nbNeighbors2D, percentile, estimatedRadii); + filterObservations2D(sfmData, params.filterObservations2D, estimatedRadii); ALICEVISION_LOG_INFO("Filtering observations in 2D: done."); if(outputRadiiFilename.empty()) From 8a0b9aa79afa0115fb40c798a4c80ec322f045eb Mon Sep 17 00:00:00 2001 From: Abdelrahman AL MAROUK Date: Wed, 27 Sep 2023 09:25:08 +0200 Subject: [PATCH 28/37] [filterSfM] fix for linux build --- src/software/pipeline/main_filterSfM.cpp | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/src/software/pipeline/main_filterSfM.cpp b/src/software/pipeline/main_filterSfM.cpp index 862a7b1759..e3df67a73c 100644 --- a/src/software/pipeline/main_filterSfM.cpp +++ b/src/software/pipeline/main_filterSfM.cpp @@ -131,6 +131,12 @@ std::istream& operator>>(std::istream& in, FilterParams& params) return in; } +// required for linux build +inline std::ostream& operator<<(std::ostream& os, const FilterParams& params) +{ + return os; +} + struct ObservationsAdaptator { using Derived = ObservationsAdaptator; //!< In this case the dataset class is myself. From 416e6a861994048c23e296a723a829cb2f0c8c4e Mon Sep 17 00:00:00 2001 From: Abdelrahman AL MAROUK Date: Wed, 27 Sep 2023 17:34:21 +0200 Subject: [PATCH 29/37] [filterSfM] refacto of complex methods --- src/software/pipeline/main_filterSfM.cpp | 787 ++++++++++++----------- 1 file changed, 418 insertions(+), 369 deletions(-) diff --git a/src/software/pipeline/main_filterSfM.cpp b/src/software/pipeline/main_filterSfM.cpp index e3df67a73c..36f658debd 100644 --- a/src/software/pipeline/main_filterSfM.cpp +++ b/src/software/pipeline/main_filterSfM.cpp @@ -304,216 +304,229 @@ ObservationsPerView getObservationsPerViews(SfMData& sfmData) return observationsPerView; } -bool filterLandmarks(SfMData& sfmData, const FilterParams::FilterLandmarksParams& params) +void filterLandmarks_step1(SfMData& sfmData, + const FilterParams::FilterLandmarksParams::FilterLandmarksStep1Params& params, + std::vector& landmarksData) { - std::vector landmarksData; - // Step 1 + const auto& initialNbLandmarks = sfmData.getLandmarks().size(); + landmarksData.resize(initialNbLandmarks); + size_t i = 0; + if(params.enabled) { - const auto& initialNbLandmarks = sfmData.getLandmarks().size(); - landmarksData.resize(initialNbLandmarks); - size_t i = 0; - if(params.step1.enabled) + ALICEVISION_LOG_INFO("Removing landmarks having an insufficient number of observations: started."); + for(auto& it : sfmData.getLandmarks()) { - ALICEVISION_LOG_INFO("Removing landmarks having an insufficient number of observations: started."); - for(auto& it : sfmData.getLandmarks()) - { - if(it.second.observations.size() < params.step1.minNbObservationsPerLandmark) - continue; - landmarksData[i++] = it.second; - } - landmarksData.resize(i); - ALICEVISION_LOG_INFO( - "Removed " << (initialNbLandmarks - i) << - " landmarks out of " << initialNbLandmarks << - ", i.e. " << ((initialNbLandmarks - i) * 100.f / initialNbLandmarks) << - " % " - ); - ALICEVISION_LOG_INFO("Removing landmarks having an insufficient number of observations: done."); + if(it.second.observations.size() < params.minNbObservationsPerLandmark) + continue; + landmarksData[i++] = it.second; } - else + landmarksData.resize(i); + ALICEVISION_LOG_INFO("Removed " << (initialNbLandmarks - i) << " landmarks out of " << initialNbLandmarks + << ", i.e. " << ((initialNbLandmarks - i) * 100.f / initialNbLandmarks) + << " % "); + ALICEVISION_LOG_INFO("Removing landmarks having an insufficient number of observations: done."); + } + else + { + for(auto& it : sfmData.getLandmarks()) { - for(auto& it : sfmData.getLandmarks()) - { - landmarksData[i++] = it.second; - } + landmarksData[i++] = it.second; } } +} - // Step 2 - if(params.step2.enabled) - { - ALICEVISION_LOG_INFO("Removing landmarks with dissimilar observations of 3D landmark neighbors: started."); - - ALICEVISION_LOG_INFO("Build nanoflann KdTree index for landmarks in 3D."); - LandmarksAdaptator data(landmarksData); - KdTree3D tree(3, data, nanoflann::KDTreeSingleIndexAdaptorParams(MAX_LEAF_ELEMENTS)); - tree.buildIndex(); - ALICEVISION_LOG_INFO("KdTree created for " << landmarksData.size() << " points."); - - ALICEVISION_LOG_INFO("Computing landmarks neighbors: started."); - // note that the landmark is a neighbor to itself with zero distance, hence the +/- 1 - int nbNeighbors_ = std::min(params.step2.nbNeighbors3D, static_cast(landmarksData.size() - 1)) + 1; - // contains the observing view ids and neighbors for each landmark - std::vector, std::vector>> viewData(landmarksData.size()); +void filterLandmarks_step2(SfMData& sfmData, + const FilterParams::FilterLandmarksParams::FilterLandmarksStep2Params& params, + std::vector& landmarksData) +{ + ALICEVISION_LOG_INFO("Removing landmarks with dissimilar observations of 3D landmark neighbors: started."); + + ALICEVISION_LOG_INFO("Build nanoflann KdTree index for landmarks in 3D."); + LandmarksAdaptator data(landmarksData); + KdTree3D tree(3, data, nanoflann::KDTreeSingleIndexAdaptorParams(MAX_LEAF_ELEMENTS)); + tree.buildIndex(); + ALICEVISION_LOG_INFO("KdTree created for " << landmarksData.size() << " points."); + + ALICEVISION_LOG_INFO("Computing landmarks neighbors: started."); + // note that the landmark is a neighbor to itself with zero distance, hence the +/- 1 + int nbNeighbors_ = std::min(params.nbNeighbors3D, static_cast(landmarksData.size() - 1)) + 1; + // contains the observing view ids and neighbors for each landmark + std::vector, std::vector>> viewData(landmarksData.size()); #pragma omp parallel for - for(auto i = 0; i < landmarksData.size(); i++) + for(auto i = 0; i < landmarksData.size(); i++) + { + const sfmData::Landmark& landmark = landmarksData[i]; + const auto& nbObservations = landmark.observations.size(); + auto& [viewIds, neighbors] = viewData[i]; + viewIds.reserve(nbObservations); + for(const auto& observationPair : landmark.observations) { - const sfmData::Landmark& landmark = landmarksData[i]; - const auto& nbObservations = landmark.observations.size(); - auto& [viewIds, neighbors] = viewData[i]; - viewIds.reserve(nbObservations); - for(const auto& observationPair : landmark.observations) - { - const IndexT viewId = observationPair.first; - viewIds.push_back(viewId); - } - // sort by ascending view id order for consequent faster access - std::stable_sort(viewIds.begin(), viewIds.end()); - - neighbors.resize(nbNeighbors_); - std::vector weights_(nbNeighbors_); - tree.knnSearch(landmark.X.data(), nbNeighbors_, &neighbors[0], &weights_[0]); - // a landmark is a neighbor to itself with zero distance, remove it - neighbors.erase(neighbors.begin()); + const IndexT viewId = observationPair.first; + viewIds.push_back(viewId); } - ALICEVISION_LOG_INFO("Computing landmarks neighbors: done."); + // sort by ascending view id order for consequent faster access + std::stable_sort(viewIds.begin(), viewIds.end()); + + neighbors.resize(nbNeighbors_); + std::vector weights_(nbNeighbors_); + tree.knnSearch(landmark.X.data(), nbNeighbors_, &neighbors[0], &weights_[0]); + // a landmark is a neighbor to itself with zero distance, remove it + neighbors.erase(neighbors.begin()); + } + ALICEVISION_LOG_INFO("Computing landmarks neighbors: done."); - ALICEVISION_LOG_INFO("Identifying landmarks to remove: started."); - std::vector toRemove(landmarksData.size(), false); - size_t nbToRemove = 0; - const auto& initialNbLandmarks = landmarksData.size(); + ALICEVISION_LOG_INFO("Identifying landmarks to remove: started."); + std::vector toRemove(landmarksData.size(), false); + size_t nbToRemove = 0; + const auto& initialNbLandmarks = landmarksData.size(); #pragma omp parallel for reduction(+ : nbToRemove) - for(auto i = 0; i < landmarksData.size(); i++) + for(auto i = 0; i < landmarksData.size(); i++) + { + int score = 0; + auto& [viewIds, neighbors] = viewData[i]; + for(auto j = 0; j < neighbors.size(); j++) { - int score = 0; - auto& [viewIds, neighbors] = viewData[i]; - for(auto j = 0; j < neighbors.size(); j++) + int scoreNeighbor = 0; + const auto& neighborId = neighbors[j]; + const auto& viewIds_neighbor = viewData[neighborId].first; + // loop over common views + auto viewIds_it = viewIds.begin(); + auto viewIds_neighbor_it = viewIds_neighbor.begin(); + while(viewIds_it != viewIds.end() && viewIds_neighbor_it != viewIds_neighbor.end()) { - int scoreNeighbor = 0; - const auto& neighborId = neighbors[j]; - const auto& viewIds_neighbor = viewData[neighborId].first; - // loop over common views - auto viewIds_it = viewIds.begin(); - auto viewIds_neighbor_it = viewIds_neighbor.begin(); - while(viewIds_it != viewIds.end() && viewIds_neighbor_it != viewIds_neighbor.end()) + if(*viewIds_it < *viewIds_neighbor_it) { - if(*viewIds_it < *viewIds_neighbor_it) - { - ++viewIds_it; - } - else + ++viewIds_it; + } + else + { + // if same view + if(!(*viewIds_neighbor_it < *viewIds_it)) { - // if same view - if(!(*viewIds_neighbor_it < *viewIds_it)) + scoreNeighbor++; + if(scoreNeighbor == params.maxNbObservationsPerLandmark) { - scoreNeighbor++; - if(scoreNeighbor == params.step2.maxNbObservationsPerLandmark) - { - score++; - break; - } - ++viewIds_it; + score++; + break; } - ++viewIds_neighbor_it; + ++viewIds_it; } + ++viewIds_neighbor_it; } } - if(score < params.step2.nbNeighbors3D / 2) - { - toRemove[i] = true; - nbToRemove++; - } } - ALICEVISION_LOG_INFO("Identified " << nbToRemove << " landmarks to remove out of " << initialNbLandmarks - << ", i.e. " << (nbToRemove * 100.f / initialNbLandmarks) << " % "); - ALICEVISION_LOG_INFO("Identifying landmarks to remove: done."); - - ALICEVISION_LOG_INFO("Removing identified landmarks: started."); - std::vector landmarksData_filtered(landmarksData.size() - nbToRemove); - IndexT newIdx = 0; - for(size_t i = 0; i < landmarksData.size(); i++) + if(score < params.nbNeighbors3D / 2) { - if(!toRemove[i]) - { - landmarksData_filtered[newIdx++] = landmarksData[i]; - } + toRemove[i] = true; + nbToRemove++; } - landmarksData = std::move(landmarksData_filtered); - ALICEVISION_LOG_INFO("Removing identified landmarks: done."); - - ALICEVISION_LOG_INFO("Removing landmarks with dissimilar observations of 3D landmark neighbors: done."); } - - // Step 3 - if(params.step3.enabled) + ALICEVISION_LOG_INFO("Identified " << nbToRemove << " landmarks to remove out of " << initialNbLandmarks + << ", i.e. " << (nbToRemove * 100.f / initialNbLandmarks) << " % "); + ALICEVISION_LOG_INFO("Identifying landmarks to remove: done."); + + ALICEVISION_LOG_INFO("Removing identified landmarks: started."); + std::vector landmarksData_filtered(landmarksData.size() - nbToRemove); + IndexT newIdx = 0; + for(size_t i = 0; i < landmarksData.size(); i++) { - ALICEVISION_LOG_INFO("Removing landmarks with worse resolution than neighbors: started."); + if(!toRemove[i]) + { + landmarksData_filtered[newIdx++] = landmarksData[i]; + } + } + landmarksData = std::move(landmarksData_filtered); + ALICEVISION_LOG_INFO("Removing identified landmarks: done."); + + ALICEVISION_LOG_INFO("Removing landmarks with dissimilar observations of 3D landmark neighbors: done."); +} + +void filterLandmarks_step3(SfMData& sfmData, + const FilterParams::FilterLandmarksParams::FilterLandmarksStep3Params& params, + std::vector& landmarksData) +{ + ALICEVISION_LOG_INFO("Removing landmarks with worse resolution than neighbors: started."); - mvsUtils::MultiViewParams mp(sfmData, "", "", "", false); - std::vector landmarksPixSize(landmarksData.size()); + mvsUtils::MultiViewParams mp(sfmData, "", "", "", false); + std::vector landmarksPixSize(landmarksData.size()); - ALICEVISION_LOG_INFO("Computing pixel size: started."); + ALICEVISION_LOG_INFO("Computing pixel size: started."); #pragma omp parallel for - for(auto i = 0; i < landmarksData.size(); i++) + for(auto i = 0; i < landmarksData.size(); i++) + { + const Landmark& landmark = landmarksData[i]; + // compute landmark pixSize + double pixSize = 0.; + int n = 0; + for(const auto& observationPair : landmark.observations) { - const Landmark& landmark = landmarksData[i]; - // compute landmark pixSize - double pixSize = 0.; - int n = 0; - for(const auto& observationPair : landmark.observations) - { - const IndexT viewId = observationPair.first; - pixSize += mp.getCamPixelSize(Point3d(landmark.X.x(), landmark.X.y(), landmark.X.z()), - mp.getIndexFromViewId(viewId), - params.step3.useFeatureScale ? observationPair.second.scale : 1); - n++; - } - pixSize /= n; - landmarksPixSize[i] = pixSize; + const IndexT viewId = observationPair.first; + pixSize += mp.getCamPixelSize(Point3d(landmark.X.x(), landmark.X.y(), landmark.X.z()), + mp.getIndexFromViewId(viewId), + params.useFeatureScale ? observationPair.second.scale : 1); + n++; } - ALICEVISION_LOG_INFO("Computing pixel size: done."); + pixSize /= n; + landmarksPixSize[i] = pixSize; + } + ALICEVISION_LOG_INFO("Computing pixel size: done."); - ALICEVISION_LOG_INFO("Identifying landmarks to remove based on pixel size: started."); + ALICEVISION_LOG_INFO("Identifying landmarks to remove based on pixel size: started."); - ALICEVISION_LOG_INFO("Build nanoflann KdTree index for landmarks in 3D."); - LandmarksAdaptator data(landmarksData); - KdTree3D tree(3, data, nanoflann::KDTreeSingleIndexAdaptorParams(MAX_LEAF_ELEMENTS)); - tree.buildIndex(); - ALICEVISION_LOG_INFO("KdTree created for " << landmarksData.size() << " points."); + ALICEVISION_LOG_INFO("Build nanoflann KdTree index for landmarks in 3D."); + LandmarksAdaptator data(landmarksData); + KdTree3D tree(3, data, nanoflann::KDTreeSingleIndexAdaptorParams(MAX_LEAF_ELEMENTS)); + tree.buildIndex(); + ALICEVISION_LOG_INFO("KdTree created for " << landmarksData.size() << " points."); - const auto& initialNbLandmarks = landmarksData.size(); - std::vector toRemove(landmarksData.size(), false); - size_t nbToRemove = 0; + const auto& initialNbLandmarks = landmarksData.size(); + std::vector toRemove(landmarksData.size(), false); + size_t nbToRemove = 0; #pragma omp parallel for reduction(+ : nbToRemove) - for(auto i = 0; i < landmarksData.size(); i++) + for(auto i = 0; i < landmarksData.size(); i++) + { + PixSizeSearch search(landmarksPixSize[i] * params.radiusScale, landmarksPixSize, i); + bool found = tree.findNeighbors(search, landmarksData[i].X.data(), nanoflann::SearchParameters()); + if(found) { - PixSizeSearch search(landmarksPixSize[i] * params.step3.radiusScale, landmarksPixSize, i); - bool found = tree.findNeighbors(search, landmarksData[i].X.data(), nanoflann::SearchParameters()); - if(found) - { - toRemove[i] = true; - nbToRemove++; - } + toRemove[i] = true; + nbToRemove++; } - ALICEVISION_LOG_INFO("Identified " << nbToRemove << " landmarks to remove out of " << initialNbLandmarks - << ", i.e. " << (nbToRemove * 100.f / initialNbLandmarks) << " % "); - ALICEVISION_LOG_INFO("Identifying landmarks to remove: done."); - - ALICEVISION_LOG_INFO("Removing identified landmarks: started."); - std::vector landmarksData_filtered(landmarksData.size() - nbToRemove); - IndexT newIdx = 0; - for(size_t i = 0; i < landmarksData.size(); i++) + } + ALICEVISION_LOG_INFO("Identified " << nbToRemove << " landmarks to remove out of " << initialNbLandmarks + << ", i.e. " << (nbToRemove * 100.f / initialNbLandmarks) << " % "); + ALICEVISION_LOG_INFO("Identifying landmarks to remove: done."); + + ALICEVISION_LOG_INFO("Removing identified landmarks: started."); + std::vector landmarksData_filtered(landmarksData.size() - nbToRemove); + IndexT newIdx = 0; + for(size_t i = 0; i < landmarksData.size(); i++) + { + if(!toRemove[i]) { - if(!toRemove[i]) - { - landmarksData_filtered[newIdx++] = landmarksData[i]; - } + landmarksData_filtered[newIdx++] = landmarksData[i]; } - landmarksData = std::move(landmarksData_filtered); - ALICEVISION_LOG_INFO("Removing identified landmarks: done."); - - ALICEVISION_LOG_INFO("Removing landmarks with worse resolution than neighbors: done."); } + landmarksData = std::move(landmarksData_filtered); + ALICEVISION_LOG_INFO("Removing identified landmarks: done."); + + ALICEVISION_LOG_INFO("Removing landmarks with worse resolution than neighbors: done."); +} + +bool filterLandmarks(SfMData& sfmData, const FilterParams::FilterLandmarksParams& params) +{ + std::vector landmarksData; + + // Step 1 + filterLandmarks_step1(sfmData, params.step1, landmarksData); + + // Step 2 + if(params.step2.enabled) + filterLandmarks_step2(sfmData, params.step2, landmarksData); + + // Step 3 + if(params.step3.enabled) + filterLandmarks_step3(sfmData, params.step3, landmarksData); // applying modifications to Sfm data std::vector> filteredLandmarks(landmarksData.size()); @@ -525,21 +538,9 @@ bool filterLandmarks(SfMData& sfmData, const FilterParams::FilterLandmarksParams return true; } -bool filterObservations3D(SfMData& sfmData, const FilterParams::FilterObservations3DParams& params) +void computeInitialObservationScores(SfMData& sfmData, std::vector& landmarksData, + std::vector, std::vector>>& viewScoresData) { - // store in vector for faster access - std::vector landmarksData(sfmData.getLandmarks().size()); - { - size_t i = 0; - for(auto& it : sfmData.getLandmarks()) - { - landmarksData[i++] = &it.second; - } - } - - // contains the observing view ids for each landmark with their corresponding scores - std::vector, std::vector>> viewScoresData(landmarksData.size()); - ALICEVISION_LOG_INFO("Computing initial observation scores based on distance to observing view: started"); #pragma omp parallel for for(auto i = 0; i < landmarksData.size(); i++) @@ -589,152 +590,192 @@ bool filterObservations3D(SfMData& sfmData, const FilterParams::FilterObservatio } } ALICEVISION_LOG_INFO("Computing initial observation scores based on distance to observing view: done"); +} - if(params.propagationEnabled) - { - ALICEVISION_LOG_INFO("Computing landmark neighbors and distance-based weights: started"); - ALICEVISION_LOG_INFO("Build nanoflann KdTree index for landmarks in 3D."); - LandmarksAdaptator dataAdaptor(landmarksData); - KdTree3D tree(3, dataAdaptor, nanoflann::KDTreeSingleIndexAdaptorParams(MAX_LEAF_ELEMENTS)); - tree.buildIndex(); - ALICEVISION_LOG_INFO("KdTree created for " << landmarksData.size() << " points."); - // note that the landmark is a neighbor to itself with zero distance, hence the +/- 1 - int nbNeighbors_ = std::min(params.nbNeighbors3D, static_cast(landmarksData.size() - 1)) + 1; - // contains the neighbor landmarks ids with their corresponding weights - std::vector, std::vector>> neighborsData(landmarksData.size()); +void computeNeighborsInfo(std::vector& landmarksData, const FilterParams::FilterObservations3DParams& params, + std::vector, std::vector>>& neighborsData) +{ + ALICEVISION_LOG_INFO("Computing landmark neighbors and distance-based weights: started"); + ALICEVISION_LOG_INFO("Build nanoflann KdTree index for landmarks in 3D."); + LandmarksAdaptator dataAdaptor(landmarksData); + KdTree3D tree(3, dataAdaptor, nanoflann::KDTreeSingleIndexAdaptorParams(MAX_LEAF_ELEMENTS)); + tree.buildIndex(); + ALICEVISION_LOG_INFO("KdTree created for " << landmarksData.size() << " points."); + // note that the landmark is a neighbor to itself with zero distance, hence the +/- 1 + int nbNeighbors_ = std::min(params.nbNeighbors3D, static_cast(landmarksData.size() - 1)) + 1; #pragma omp parallel for - for(auto i = 0; i < landmarksData.size(); i++) + for(auto i = 0; i < landmarksData.size(); i++) + { + const sfmData::Landmark& landmark = *landmarksData[i]; + auto& [indices_, weights_] = neighborsData[i]; + indices_.resize(nbNeighbors_); + weights_.resize(nbNeighbors_); + tree.knnSearch(landmark.X.data(), nbNeighbors_, &indices_[0], &weights_[0]); + // a landmark is a neighbor to itself with zero distance, remove it + indices_.erase(indices_.begin()); + weights_.erase(weights_.begin()); + // accumulator used for normalisation + double total = 0.; + for(auto& w : weights_) { - const sfmData::Landmark& landmark = *landmarksData[i]; - auto& [indices_, weights_] = neighborsData[i]; - indices_.resize(nbNeighbors_); - weights_.resize(nbNeighbors_); - tree.knnSearch(landmark.X.data(), nbNeighbors_, &indices_[0], &weights_[0]); - // a landmark is a neighbor to itself with zero distance, remove it - indices_.erase(indices_.begin()); - weights_.erase(weights_.begin()); - // accumulator used for normalisation - double total = 0.; - for(auto& w : weights_) - { - // weight is the inverse of distance between a landmark and its neighbor - w = 1. / std::sqrt(w); - if(std::isinf(w)) - w = std::numeric_limits::max(); - total += w; - } - if(std::isinf(total)) - total = std::numeric_limits::max(); - // normalize weights - for(auto& w : weights_) - { - w /= total; - } + // weight is the inverse of distance between a landmark and its neighbor + w = 1. / std::sqrt(w); + if(std::isinf(w)) + w = std::numeric_limits::max(); + total += w; } - ALICEVISION_LOG_INFO("Computing landmark neighbors and distance-based weights: done"); - - ALICEVISION_LOG_INFO("Propagating neighbors observation scores: started"); - // new view scores at iteration t - std::vector> viewScoresData_t(landmarksData.size()); - for(auto i = 0; i < params.nbIterations; i++) + if(std::isinf(total)) + total = std::numeric_limits::max(); + // normalize weights + for(auto& w : weights_) { + w /= total; + } + } + ALICEVISION_LOG_INFO("Computing landmark neighbors and distance-based weights: done"); +} + +void computeNewScores(const std::vector& landmarksData, + const FilterParams::FilterObservations3DParams& params, + const std::vector, std::vector>>& neighborsData, + const std::vector, std::vector>>& viewScoresData, + std::vector>& viewScoresData_t) +{ #pragma omp parallel for - for(auto id = 0; id < landmarksData.size(); id++) + for(auto id = 0; id < landmarksData.size(); id++) + { + const auto& [viewIds, viewScores] = viewScoresData[id]; + auto& viewScores_acc = viewScoresData_t[id]; + // initialize to zero, will first contain the weighted average scores of neighbors + viewScores_acc.assign(viewScores.size(), 0.); + // accumulator for normalisation + double viewScores_total = 0.; + auto& [indices_, weights_] = neighborsData[id]; + for(auto j = 0; j < params.nbNeighbors3D; j++) + { + const auto& neighborId = indices_[j]; + const auto& neighborWeight = weights_[j]; + const auto& [viewIds_neighbor, viewScores_neighbor] = viewScoresData[neighborId]; + // loop over common views + auto viewIds_it = viewIds.begin(); + auto viewIds_neighbor_it = viewIds_neighbor.begin(); + auto viewScores_neighbor_it = viewScores_neighbor.begin(); + auto viewScores_acc_it = viewScores_acc.begin(); + while(viewIds_it != viewIds.end() && viewIds_neighbor_it != viewIds_neighbor.end()) { - const auto& [viewIds, viewScores] = viewScoresData[id]; - auto& viewScores_acc = viewScoresData_t[id]; - // initialize to zero, will first contain the weighted average scores of neighbors - viewScores_acc.assign(viewScores.size(), 0.); - // accumulator for normalisation - double viewScores_total = 0.; - auto& [indices_, weights_] = neighborsData[id]; - for(auto j = 0; j < params.nbNeighbors3D; j++) + if(*viewIds_it < *viewIds_neighbor_it) { - const auto& neighborId = indices_[j]; - const auto& neighborWeight = weights_[j]; - const auto& [viewIds_neighbor, viewScores_neighbor] = viewScoresData[neighborId]; - // loop over common views - auto viewIds_it = viewIds.begin(); - auto viewIds_neighbor_it = viewIds_neighbor.begin(); - auto viewScores_neighbor_it = viewScores_neighbor.begin(); - auto viewScores_acc_it = viewScores_acc.begin(); - while(viewIds_it != viewIds.end() && viewIds_neighbor_it != viewIds_neighbor.end()) - { - if(*viewIds_it < *viewIds_neighbor_it) - { - ++viewIds_it; - ++viewScores_acc_it; - } - else - { - // if same view, accumulate weighted scores - if(!(*viewIds_neighbor_it < *viewIds_it)) - { - const auto& v = *viewScores_neighbor_it * neighborWeight; - (*viewScores_acc_it) += v; - viewScores_total += v; - ++viewIds_it; - ++viewScores_acc_it; - } - ++viewIds_neighbor_it; - ++viewScores_neighbor_it; - } - } + ++viewIds_it; + ++viewScores_acc_it; } - // if common views with neighbor landmarks - if(viewScores_total != 0.) - for(auto j = 0; j < viewScores_acc.size(); j++) - { - // normalize score and apply influence factor - viewScores_acc[j] *= params.neighborsInfluence / viewScores_total; - // combine weighted neighbor scores and the landmark's own scores - viewScores_acc[j] += (1 - params.neighborsInfluence) * viewScores[j]; - } - - // dampen scores of non-chosen observations - if(params.dampingEnabled && viewScores_acc.size() <= params.maxNbObservationsPerLandmark) + else { - // sort by descending view score order - std::vector idx(viewScores_acc.size()); - std::iota(idx.begin(), idx.end(), 0); - std::stable_sort(idx.begin(), idx.end(), - [&v = viewScores_acc](size_t i1, size_t i2) { return v[i1] > v[i2]; }); - viewScores_total = 1.; - for(auto j = params.maxNbObservationsPerLandmark; j < viewScores_acc.size(); j++) + // if same view, accumulate weighted scores + if(!(*viewIds_neighbor_it < *viewIds_it)) { - const double& v = params.dampingFactor * viewScores_acc[j]; - viewScores_total += v - viewScores_acc[j]; - viewScores_acc[j] = v; + const auto& v = *viewScores_neighbor_it * neighborWeight; + (*viewScores_acc_it) += v; + viewScores_total += v; + ++viewIds_it; + ++viewScores_acc_it; } - // re-normalize - for(auto j = 0; j < viewScores_acc.size(); j++) - viewScores_acc[j] /= viewScores_total; + ++viewIds_neighbor_it; + ++viewScores_neighbor_it; } } - // Mean Squared Error - double error = 0.; - // update scores at end of iteration + } + // if common views with neighbor landmarks + if(viewScores_total != 0.) + for(auto j = 0; j < viewScores_acc.size(); j++) + { + // normalize score and apply influence factor + viewScores_acc[j] *= params.neighborsInfluence / viewScores_total; + // combine weighted neighbor scores and the landmark's own scores + viewScores_acc[j] += (1 - params.neighborsInfluence) * viewScores[j]; + } + + // dampen scores of non-chosen observations + if(params.dampingEnabled && viewScores_acc.size() <= params.maxNbObservationsPerLandmark) + { + // sort by descending view score order + std::vector idx(viewScores_acc.size()); + std::iota(idx.begin(), idx.end(), 0); + std::stable_sort(idx.begin(), idx.end(), + [&v = viewScores_acc](size_t i1, size_t i2) { return v[i1] > v[i2]; }); + viewScores_total = 1.; + for(auto j = params.maxNbObservationsPerLandmark; j < viewScores_acc.size(); j++) + { + const double& v = params.dampingFactor * viewScores_acc[j]; + viewScores_total += v - viewScores_acc[j]; + viewScores_acc[j] = v; + } + // re-normalize + for(auto j = 0; j < viewScores_acc.size(); j++) + viewScores_acc[j] /= viewScores_total; + } + } +} + +void propagateNeighborsInfo(std::vector& landmarksData, + const FilterParams::FilterObservations3DParams& params, + std::vector, std::vector>>& neighborsData, + std::vector, std::vector>>& viewScoresData) +{ + ALICEVISION_LOG_INFO("Propagating neighbors observation scores: started"); + // new view scores at iteration t + std::vector> viewScoresData_t(landmarksData.size()); + for(auto i = 0; i < params.nbIterations; i++) + { + computeNewScores(landmarksData, params, neighborsData, viewScoresData, viewScoresData_t); + + // Mean Squared Error + double error = 0.; + // update scores at end of iteration #pragma omp parallel for reduction(+ : error) - for(auto id = 0; id < landmarksData.size(); id++) + for(auto id = 0; id < landmarksData.size(); id++) + { + // compute MSE { - // compute MSE + double error_j = 0.; + for(auto j = 0; j < viewScoresData_t[id].size(); j++) { - double error_j = 0.; - for(auto j = 0; j < viewScoresData_t[id].size(); j++) - { - const auto& v = viewScoresData_t[id][j] - viewScoresData[id].second[j]; - error_j += v * v; - } - error_j /= viewScoresData_t[id].size(); - error += error_j; + const auto& v = viewScoresData_t[id][j] - viewScoresData[id].second[j]; + error_j += v * v; } - // update scores - viewScoresData[id].second = std::move(viewScoresData_t[id]); + error_j /= viewScoresData_t[id].size(); + error += error_j; } - error /= landmarksData.size(); - ALICEVISION_LOG_INFO("MSE at iteration " << i << ": " << error); + // update scores + viewScoresData[id].second = std::move(viewScoresData_t[id]); } - ALICEVISION_LOG_INFO("Propagating neighbors observation scores: done"); + error /= landmarksData.size(); + ALICEVISION_LOG_INFO("MSE at iteration " << i << ": " << error); + } + ALICEVISION_LOG_INFO("Propagating neighbors observation scores: done"); +} + +bool filterObservations3D(SfMData& sfmData, const FilterParams::FilterObservations3DParams& params) +{ + // store in vector for faster access + std::vector landmarksData(sfmData.getLandmarks().size()); + { + size_t i = 0; + for(auto& it : sfmData.getLandmarks()) + { + landmarksData[i++] = &it.second; + } + } + + // contains the observing view ids for each landmark with their corresponding scores + std::vector, std::vector>> viewScoresData(landmarksData.size()); + computeInitialObservationScores(sfmData, landmarksData, viewScoresData); + + if(params.propagationEnabled) + { + // contains the neighbor landmarks ids with their corresponding weights + std::vector, std::vector>> neighborsData(landmarksData.size()); + computeNeighborsInfo(landmarksData, params, neighborsData); + propagateNeighborsInfo(landmarksData, params, neighborsData, viewScoresData); } ALICEVISION_LOG_INFO("Selecting observations with best scores: started"); @@ -768,6 +809,72 @@ bool filterObservations3D(SfMData& sfmData, const FilterParams::FilterObservatio return true; } +double filter2DView(SfMData& sfmData, const FilterParams::FilterObservations2DParams& params, const IndexT& viewId, + std::vector& observations, std::vector& landmarks) +{ + // Build nanoflann KdTree index for projected landmarks in 2D + ObservationsAdaptator data(observations); + KdTree2D tree(2, data, nanoflann::KDTreeSingleIndexAdaptorParams(MAX_LEAF_ELEMENTS)); + tree.buildIndex(); + + // note that the observation is a neighbor to itself with zero distance, hence the +/- 1 + size_t nbNeighbors_ = std::min(static_cast(params.nbNeighbors2D), observations.size() - 1) + 1; + // average neighbors distance for each observation + std::vector means(observations.size()); + const std::size_t cacheSize = 1000; + // accumulator for quantile computation + accumulator_set>> acc(tag::tail::cache_size = cacheSize); + for(auto j = 0; j < observations.size(); j++) + { + // Find neighbors and the corresponding distances + const auto& obs = observations[j]; + std::vector indices_(nbNeighbors_); + std::vector distances_(nbNeighbors_); + tree.knnSearch(obs.x.data(), nbNeighbors_, &indices_[0], &distances_[0]); + // returned distances are L2 -> apply squared root + std::transform(distances_.begin(), distances_.end(), distances_.begin(), + static_cast(std::sqrt)); + // average distance + const auto& mean = std::accumulate(distances_.begin(), distances_.end(), 0.0) / (nbNeighbors_ - 1); + means[j] = mean; + // update accumulator + acc(mean); + } + double mean_max = std::numeric_limits::max(); + // check to avoid exception + if(params.percentile != 1.f) + mean_max = quantile(acc, quantile_probability = params.percentile); + // estimated mask radius is the average of distance means + // quantile is used to avoid outlier bias + double radius = quantile(acc, quantile_probability = (1.f - params.percentile) * 0.5f); + // check if estimated radius is too large + { + const View& view = *(sfmData.getViews().at(viewId)); + double radiusMax = + params.maskRadiusThreshold * 0.5 * (view.getImage().getWidth() + view.getImage().getHeight()); + if(radius > radiusMax) + radius = radiusMax; + } + + // filter outlier observations + std::vector filteredObservations; + std::vector filteredLandmarks; + filteredObservations.reserve(observations.size()); + filteredLandmarks.reserve(landmarks.size()); + for(auto j = 0; j < observations.size(); j++) + if(means[j] < mean_max) + { + filteredObservations.push_back(observations[j]); + filteredLandmarks.push_back(landmarks[j]); + } + filteredObservations.shrink_to_fit(); + filteredLandmarks.shrink_to_fit(); + observations = std::move(filteredObservations); + landmarks = std::move(filteredLandmarks); + + return radius; +} + bool filterObservations2D(SfMData& sfmData, const FilterParams::FilterObservations2DParams& params, HashMap& estimatedRadii) { @@ -786,66 +893,7 @@ bool filterObservations2D(SfMData& sfmData, const FilterParams::FilterObservatio continue; auto& observations = observationsIt->second.first; auto& landmarks = observationsIt->second.second; - - // Build nanoflann KdTree index for projected landmarks in 2D - ObservationsAdaptator data(observations); - KdTree2D tree(2, data, nanoflann::KDTreeSingleIndexAdaptorParams(MAX_LEAF_ELEMENTS)); - tree.buildIndex(); - - // note that the observation is a neighbor to itself with zero distance, hence the +/- 1 - size_t nbNeighbors_ = std::min(static_cast(params.nbNeighbors2D), observations.size() - 1) + 1; - // average neighbors distance for each observation - std::vector means(observations.size()); - const std::size_t cacheSize = 1000; - // accumulator for quantile computation - accumulator_set>> acc(tag::tail::cache_size = cacheSize); - for(auto j = 0; j < observations.size(); j++) - { - // Find neighbors and the corresponding distances - const auto& obs = observations[j]; - std::vector indices_(nbNeighbors_); - std::vector distances_(nbNeighbors_); - tree.knnSearch(obs.x.data(), nbNeighbors_, &indices_[0], &distances_[0]); - // returned distances are L2 -> apply squared root - std::transform(distances_.begin(), distances_.end(), distances_.begin(), - static_cast(std::sqrt)); - // average distance - const auto& mean = std::accumulate(distances_.begin(), distances_.end(), 0.0) / (nbNeighbors_ - 1); - means[j] = mean; - // update accumulator - acc(mean); - } - double mean_max = std::numeric_limits::max(); - // check to avoid exception - if(params.percentile != 1.f) - mean_max = quantile(acc, quantile_probability = params.percentile); - // estimated mask radius is the average of distance means - // quantile is used to avoid outlier bias - double radius = quantile(acc, quantile_probability = (1.f - params.percentile) * 0.5f); - // check if estimated radius is too large - { - const View& view = *(sfmData.getViews().at(viewId)); - double radiusMax = params.maskRadiusThreshold * 0.5 * (view.getImage().getWidth() + view.getImage().getHeight()); - if(radius > radiusMax) - radius = radiusMax; - } - estimatedRadii_[i] = radius; - - // filter outlier observations - std::vector filteredObservations; - std::vector filteredLandmarks; - filteredObservations.reserve(observations.size()); - filteredLandmarks.reserve(landmarks.size()); - for(auto j = 0; j < observations.size(); j++) - if(means[j] < mean_max) - { - filteredObservations.push_back(observations[j]); - filteredLandmarks.push_back(landmarks[j]); - } - filteredObservations.shrink_to_fit(); - filteredLandmarks.shrink_to_fit(); - observations = std::move(filteredObservations); - landmarks = std::move(filteredLandmarks); + estimatedRadii_[i] = filter2DView(sfmData, params, viewId, observations, landmarks); } // clear and update landmark observations @@ -887,6 +935,7 @@ bool filterObservations2D(SfMData& sfmData, const FilterParams::FilterObservatio int aliceVision_main(int argc, char *argv[]) { + std::cout << "hello" << std::endl; // command-line parameters std::string inputSfmFilename; From eb36bccdd9c299f818d0829c2370cbe9aff5e093 Mon Sep 17 00:00:00 2001 From: Abdelrahman AL MAROUK Date: Tue, 3 Oct 2023 17:32:14 +0200 Subject: [PATCH 30/37] [WIP][filterSfM] multiple improvements and fixes - resolve duplicate features issue - do not normalize propagated neighbors scores - fix NaN MSE issue - fix radius estimation caused by quantile issue --- src/software/pipeline/main_filterSfM.cpp | 154 ++++++++++++++++------- 1 file changed, 106 insertions(+), 48 deletions(-) diff --git a/src/software/pipeline/main_filterSfM.cpp b/src/software/pipeline/main_filterSfM.cpp index 36f658debd..85aaf2c5da 100644 --- a/src/software/pipeline/main_filterSfM.cpp +++ b/src/software/pipeline/main_filterSfM.cpp @@ -221,27 +221,30 @@ using KdTree3D = nanoflann::KDTreeSingleIndexAdaptor< nanoflann::L2_Simple_Adaptor, LandmarksAdaptator, 3 /* dim */,size_t>; -class RadiusKnnSearch +class PixSizeSearch { public: const double _radius_sq; - const int _nb_neighbors; - int nb_found = 0; + const std::vector& _pixSize; + const double _pixSize_i; + bool found = false; - inline RadiusKnnSearch(double radius, int k) + inline PixSizeSearch(double radius, const std::vector& pixSize, size_t i) : _radius_sq(radius * radius) - , _nb_neighbors(k) + , _pixSize(pixSize) + , _pixSize_i(pixSize[i]) { } - inline bool full() const { return nb_found == _nb_neighbors; } + inline bool full() const { return found; } - inline bool addPoint(double dist, IndexT index) + inline bool addPoint(double dist, size_t index) { - if(dist < _radius_sq) + // strict comparison because a point in the tree can be a neighbor of itself + if(dist < _radius_sq && _pixSize[index] < _pixSize_i) { - nb_found++; - return nb_found < _nb_neighbors; + found = true; + return false; } return true; } @@ -249,35 +252,76 @@ class RadiusKnnSearch inline double worstDist() const { return _radius_sq; } }; -class PixSizeSearch +class KnnNonZeroSearch { + public: - const double _radius_sq; - const std::vector& _pixSize; - const double _pixSize_i; - bool found = false; + size_t* indices; + double* dists; + int capacity; + int count; + double _epsilon = 1e-6; + + inline KnnNonZeroSearch(int capacity_) + : indices(nullptr) + , dists(nullptr) + , capacity(capacity_) + , count(0) + { + } - inline PixSizeSearch(double radius, const std::vector& pixSize, size_t i) - : _radius_sq(radius * radius) - , _pixSize(pixSize) - , _pixSize_i(pixSize[i]) + inline void init(size_t* indices_, double* dists_) { + indices = indices_; + dists = dists_; + count = 0; + if(capacity) + dists[capacity - 1] = (std::numeric_limits::max)(); } - inline bool full() const { return found; } + inline int size() const { return count; } + + inline bool full() const { return count == capacity; } inline bool addPoint(double dist, size_t index) { - // strict comparison because a point in the tree can be a neighbor of itself - if(dist < _radius_sq && _pixSize[index] < _pixSize_i) + if(dist < _epsilon) + return true; + int i; + bool isDuplicate = false; + for(i = count; i > 0; --i) { - found = true; - return false; + if(!(dists[i - 1] > dist)) + { + isDuplicate = (dists[i - 1] == dist); + break; + } } + if(isDuplicate) + { + if(indices[i - 1] > index) + indices[i - 1] = index; + } + else + { + if(i < capacity) + { + if(i > 0) + if(count < capacity) + std::move_backward(&indices[i], &indices[count], &indices[count + 1]); + else + std::move_backward(&indices[i], &indices[capacity - 1], &indices[capacity]); + dists[i] = dist; + indices[i] = index; + } + if(count < capacity) + count++; + } + // tell caller that the search shall continue return true; } - inline double worstDist() const { return _radius_sq; } + inline double worstDist() const { return dists[capacity - 1]; } }; using ObservationsPerView = stl::flat_map, std::vector>>; @@ -349,7 +393,7 @@ void filterLandmarks_step2(SfMData& sfmData, ALICEVISION_LOG_INFO("Computing landmarks neighbors: started."); // note that the landmark is a neighbor to itself with zero distance, hence the +/- 1 - int nbNeighbors_ = std::min(params.nbNeighbors3D, static_cast(landmarksData.size() - 1)) + 1; + int nbNeighbors_ = params.nbNeighbors3D; // contains the observing view ids and neighbors for each landmark std::vector, std::vector>> viewData(landmarksData.size()); #pragma omp parallel for @@ -369,9 +413,12 @@ void filterLandmarks_step2(SfMData& sfmData, neighbors.resize(nbNeighbors_); std::vector weights_(nbNeighbors_); - tree.knnSearch(landmark.X.data(), nbNeighbors_, &neighbors[0], &weights_[0]); - // a landmark is a neighbor to itself with zero distance, remove it - neighbors.erase(neighbors.begin()); + KnnNonZeroSearch resultSet(nbNeighbors_); + resultSet.init(&neighbors[0], &weights_[0]); + tree.findNeighbors(resultSet, landmark.X.data()); + const auto& nbFound = resultSet.size(); + neighbors.resize(nbFound); + weights_.resize(nbFound); } ALICEVISION_LOG_INFO("Computing landmarks neighbors: done."); @@ -601,8 +648,7 @@ void computeNeighborsInfo(std::vector& landmarksData, const FilterPar KdTree3D tree(3, dataAdaptor, nanoflann::KDTreeSingleIndexAdaptorParams(MAX_LEAF_ELEMENTS)); tree.buildIndex(); ALICEVISION_LOG_INFO("KdTree created for " << landmarksData.size() << " points."); - // note that the landmark is a neighbor to itself with zero distance, hence the +/- 1 - int nbNeighbors_ = std::min(params.nbNeighbors3D, static_cast(landmarksData.size() - 1)) + 1; + int nbNeighbors_ = params.nbNeighbors3D; #pragma omp parallel for for(auto i = 0; i < landmarksData.size(); i++) { @@ -610,22 +656,20 @@ void computeNeighborsInfo(std::vector& landmarksData, const FilterPar auto& [indices_, weights_] = neighborsData[i]; indices_.resize(nbNeighbors_); weights_.resize(nbNeighbors_); - tree.knnSearch(landmark.X.data(), nbNeighbors_, &indices_[0], &weights_[0]); - // a landmark is a neighbor to itself with zero distance, remove it - indices_.erase(indices_.begin()); - weights_.erase(weights_.begin()); + KnnNonZeroSearch resultSet(nbNeighbors_); + resultSet.init(&indices_[0], &weights_[0]); + tree.findNeighbors(resultSet, landmark.X.data()); + const auto& nbFound = resultSet.size(); + indices_.resize(nbFound); + weights_.resize(nbFound); // accumulator used for normalisation double total = 0.; for(auto& w : weights_) { // weight is the inverse of distance between a landmark and its neighbor w = 1. / std::sqrt(w); - if(std::isinf(w)) - w = std::numeric_limits::max(); total += w; } - if(std::isinf(total)) - total = std::numeric_limits::max(); // normalize weights for(auto& w : weights_) { @@ -651,7 +695,7 @@ void computeNewScores(const std::vector& landmarksData, // accumulator for normalisation double viewScores_total = 0.; auto& [indices_, weights_] = neighborsData[id]; - for(auto j = 0; j < params.nbNeighbors3D; j++) + for(auto j = 0; j < indices_.size(); j++) { const auto& neighborId = indices_[j]; const auto& neighborWeight = weights_[j]; @@ -689,7 +733,8 @@ void computeNewScores(const std::vector& landmarksData, for(auto j = 0; j < viewScores_acc.size(); j++) { // normalize score and apply influence factor - viewScores_acc[j] *= params.neighborsInfluence / viewScores_total; + //viewScores_acc[j] *= params.neighborsInfluence / viewScores_total; + viewScores_acc[j] *= params.neighborsInfluence; // combine weighted neighbor scores and the landmark's own scores viewScores_acc[j] += (1 - params.neighborsInfluence) * viewScores[j]; } @@ -742,7 +787,8 @@ void propagateNeighborsInfo(std::vector& landmarksData, const auto& v = viewScoresData_t[id][j] - viewScoresData[id].second[j]; error_j += v * v; } - error_j /= viewScoresData_t[id].size(); + if(error_j > 0.) + error_j /= viewScoresData_t[id].size(); error += error_j; } // update scores @@ -799,6 +845,9 @@ bool filterObservations3D(SfMData& sfmData, const FilterParams::FilterObservatio Observations filteredObservations; for(auto j = 0; j < params.maxNbObservationsPerLandmark; j++) { + const auto& viewScore = viewScores[idx[j]]; + if(viewScore < 0.3) + break; const auto& viewId = viewIds[idx[j]]; filteredObservations[viewId] = landmark.observations[viewId]; } @@ -818,24 +867,31 @@ double filter2DView(SfMData& sfmData, const FilterParams::FilterObservations2DPa tree.buildIndex(); // note that the observation is a neighbor to itself with zero distance, hence the +/- 1 - size_t nbNeighbors_ = std::min(static_cast(params.nbNeighbors2D), observations.size() - 1) + 1; + size_t nbNeighbors_ = params.nbNeighbors2D; // average neighbors distance for each observation - std::vector means(observations.size()); + std::vector means(observations.size(), std::numeric_limits::max()); const std::size_t cacheSize = 1000; // accumulator for quantile computation - accumulator_set>> acc(tag::tail::cache_size = cacheSize); + accumulator_set>> acc(tag::tail::cache_size = cacheSize); for(auto j = 0; j < observations.size(); j++) { // Find neighbors and the corresponding distances const auto& obs = observations[j]; std::vector indices_(nbNeighbors_); std::vector distances_(nbNeighbors_); - tree.knnSearch(obs.x.data(), nbNeighbors_, &indices_[0], &distances_[0]); + KnnNonZeroSearch resultSet(nbNeighbors_); + resultSet.init(&indices_[0], &distances_[0]); + tree.findNeighbors(resultSet, obs.x.data()); + const auto& nbFound = resultSet.size(); + if(nbFound == 0) + continue; + indices_.resize(nbFound); + distances_.resize(nbFound); // returned distances are L2 -> apply squared root std::transform(distances_.begin(), distances_.end(), distances_.begin(), static_cast(std::sqrt)); // average distance - const auto& mean = std::accumulate(distances_.begin(), distances_.end(), 0.0) / (nbNeighbors_ - 1); + const auto& mean = std::accumulate(distances_.begin(), distances_.end(), 0.0) / nbFound; means[j] = mean; // update accumulator acc(mean); @@ -846,7 +902,7 @@ double filter2DView(SfMData& sfmData, const FilterParams::FilterObservations2DPa mean_max = quantile(acc, quantile_probability = params.percentile); // estimated mask radius is the average of distance means // quantile is used to avoid outlier bias - double radius = quantile(acc, quantile_probability = (1.f - params.percentile) * 0.5f); + double radius = median(acc); // check if estimated radius is too large { const View& view = *(sfmData.getViews().at(viewId)); @@ -854,6 +910,8 @@ double filter2DView(SfMData& sfmData, const FilterParams::FilterObservations2DPa params.maskRadiusThreshold * 0.5 * (view.getImage().getWidth() + view.getImage().getHeight()); if(radius > radiusMax) radius = radiusMax; + if(mean_max > radiusMax) + mean_max = radiusMax; } // filter outlier observations From 5828fa10f140bc1b73f183de374709e30ac8c609 Mon Sep 17 00:00:00 2001 From: Abdelrahman AL MAROUK Date: Wed, 4 Oct 2023 08:57:02 +0200 Subject: [PATCH 31/37] [WIP][filterSfM] fix bug in neighbor search --- src/software/pipeline/main_filterSfM.cpp | 35 +++++++++++++++--------- 1 file changed, 22 insertions(+), 13 deletions(-) diff --git a/src/software/pipeline/main_filterSfM.cpp b/src/software/pipeline/main_filterSfM.cpp index 85aaf2c5da..f895655f0a 100644 --- a/src/software/pipeline/main_filterSfM.cpp +++ b/src/software/pipeline/main_filterSfM.cpp @@ -306,11 +306,19 @@ class KnnNonZeroSearch { if(i < capacity) { - if(i > 0) - if(count < capacity) + if(i != count) + { + if (count < capacity) + { std::move_backward(&indices[i], &indices[count], &indices[count + 1]); + std::move_backward(&dists[i], &dists[count], &dists[count + 1]); + } else + { std::move_backward(&indices[i], &indices[capacity - 1], &indices[capacity]); + std::move_backward(&dists[i], &dists[capacity - 1], &dists[capacity]); + } + } dists[i] = dist; indices[i] = index; } @@ -729,15 +737,15 @@ void computeNewScores(const std::vector& landmarksData, } } // if common views with neighbor landmarks - if(viewScores_total != 0.) - for(auto j = 0; j < viewScores_acc.size(); j++) - { - // normalize score and apply influence factor - //viewScores_acc[j] *= params.neighborsInfluence / viewScores_total; - viewScores_acc[j] *= params.neighborsInfluence; - // combine weighted neighbor scores and the landmark's own scores - viewScores_acc[j] += (1 - params.neighborsInfluence) * viewScores[j]; - } + //if(viewScores_total != 0.) + for(auto j = 0; j < viewScores_acc.size(); j++) + { + // normalize score and apply influence factor + //viewScores_acc[j] *= params.neighborsInfluence / viewScores_total; + viewScores_acc[j] *= params.neighborsInfluence; + // combine weighted neighbor scores and the landmark's own scores + viewScores_acc[j] += (1 - params.neighborsInfluence) * viewScores[j]; + } // dampen scores of non-chosen observations if(params.dampingEnabled && viewScores_acc.size() <= params.maxNbObservationsPerLandmark) @@ -843,11 +851,12 @@ bool filterObservations3D(SfMData& sfmData, const FilterParams::FilterObservatio // keep only observations with best scores Observations filteredObservations; + //double threshold = 0.1 / params.maxNbObservationsPerLandmark; for(auto j = 0; j < params.maxNbObservationsPerLandmark; j++) { const auto& viewScore = viewScores[idx[j]]; - if(viewScore < 0.3) - break; + /*if(viewScore < threshold) + break;*/ const auto& viewId = viewIds[idx[j]]; filteredObservations[viewId] = landmark.observations[viewId]; } From 5c5b7d4bdda04943f4ef61c1a98db159d0ccd9f8 Mon Sep 17 00:00:00 2001 From: Abdelrahman AL MAROUK Date: Wed, 4 Oct 2023 14:22:22 +0200 Subject: [PATCH 32/37] [WIP][filterSfM] augment observations with neighbors observations --- src/software/pipeline/main_filterSfM.cpp | 105 +++++++++++++++++++---- 1 file changed, 88 insertions(+), 17 deletions(-) diff --git a/src/software/pipeline/main_filterSfM.cpp b/src/software/pipeline/main_filterSfM.cpp index f895655f0a..820f76657e 100644 --- a/src/software/pipeline/main_filterSfM.cpp +++ b/src/software/pipeline/main_filterSfM.cpp @@ -400,7 +400,6 @@ void filterLandmarks_step2(SfMData& sfmData, ALICEVISION_LOG_INFO("KdTree created for " << landmarksData.size() << " points."); ALICEVISION_LOG_INFO("Computing landmarks neighbors: started."); - // note that the landmark is a neighbor to itself with zero distance, hence the +/- 1 int nbNeighbors_ = params.nbNeighbors3D; // contains the observing view ids and neighbors for each landmark std::vector, std::vector>> viewData(landmarksData.size()); @@ -687,6 +686,60 @@ void computeNeighborsInfo(std::vector& landmarksData, const FilterPar ALICEVISION_LOG_INFO("Computing landmark neighbors and distance-based weights: done"); } +void augmentInitialObservations(std::vector, std::vector>>& viewScoresData, + std::vector, std::vector>>& neighborsData) +{ + ALICEVISION_LOG_INFO("Augment initial observations based on neighbors observations: started"); + std::vector> viewScoresData_t(viewScoresData.size()); + // for each landmark, identify the new observations coming from its neighborhood +#pragma omp parallel for + for(int id = 0; id < viewScoresData.size(); id++) + { + const auto& viewIds = viewScoresData[id].first; + auto& viewIds_t = viewScoresData_t[id]; + const auto& neighborIds = neighborsData[id].first; + + for (const auto& neighborId : neighborIds) + { + const auto& viewIds_neighbor = viewScoresData[neighborId].first; + std::vector viewIds_union; + std::set_union(viewIds_neighbor.begin(), viewIds_neighbor.end(), viewIds_t.begin(), viewIds_t.end(), + std::back_inserter(viewIds_union)); + viewIds_t = std::move(viewIds_union); + } + std::vector viewIds_toAdd; + std::set_difference(viewIds_t.begin(), viewIds_t.end(), viewIds.begin(), viewIds.end(), + std::back_inserter(viewIds_toAdd)); + viewIds_t = std::move(viewIds_toAdd); + } + // for each landmark, add the previously identified observations + #pragma omp parallel for + for(int id = 0; id < viewScoresData.size(); id++) + { + auto& [viewIds, viewScores] = viewScoresData[id]; + const auto& viewIds_toAdd = viewScoresData_t[id]; + viewIds.insert(viewIds.end(), viewIds_toAdd.begin(), viewIds_toAdd.end()); + viewScores.insert(viewScores.end(), viewIds_toAdd.size(), 0.); + + // sort by ascending view id order + // for consequent faster access + + // indices that sort the view ids + std::vector idx(viewIds.size()); + std::iota(idx.begin(), idx.end(), 0); + std::stable_sort(idx.begin(), idx.end(), [&v = viewIds](size_t i1, size_t i2) { return v[i1] < v[i2]; }); + // apply sorting to both view ids and view scores for correspondance + auto ids_temp = viewIds; + auto scores_temp = viewScores; + for(auto j = 0; j < viewIds.size(); j++) + { + viewIds[j] = ids_temp[idx[j]]; + viewScores[j] = scores_temp[idx[j]]; + } + } + ALICEVISION_LOG_INFO("Augment initial observations based on neighbors observations: done"); +} + void computeNewScores(const std::vector& landmarksData, const FilterParams::FilterObservations3DParams& params, const std::vector, std::vector>>& neighborsData, @@ -808,6 +861,25 @@ void propagateNeighborsInfo(std::vector& landmarksData, ALICEVISION_LOG_INFO("Propagating neighbors observation scores: done"); } +void removeNonObservedLandmarks(SfMData& sfmData) +{ + // remove landmarks with no remaining observations + ALICEVISION_LOG_INFO("Remove non-observed landmarks: started"); + const auto& initialNbLandmarks = sfmData.getLandmarks().size(); + for(auto it = sfmData.getLandmarks().begin(); it != sfmData.getLandmarks().end();) + { + if(it->second.observations.size() == 0) + it = sfmData.getLandmarks().erase(it); + else + ++it; + } + const auto& modifiedNbLandmarks = sfmData.getLandmarks().size(); + ALICEVISION_LOG_INFO( + "Removed " << (initialNbLandmarks - modifiedNbLandmarks) << " landmarks out of " << initialNbLandmarks + << ", i.e. " << ((initialNbLandmarks - modifiedNbLandmarks) * 100.f / initialNbLandmarks) << " % "); + ALICEVISION_LOG_INFO("Remove non-observed landmarks: done"); +} + bool filterObservations3D(SfMData& sfmData, const FilterParams::FilterObservations3DParams& params) { // store in vector for faster access @@ -829,6 +901,7 @@ bool filterObservations3D(SfMData& sfmData, const FilterParams::FilterObservatio // contains the neighbor landmarks ids with their corresponding weights std::vector, std::vector>> neighborsData(landmarksData.size()); computeNeighborsInfo(landmarksData, params, neighborsData); + augmentInitialObservations(viewScoresData, neighborsData); propagateNeighborsInfo(landmarksData, params, neighborsData, viewScoresData); } @@ -837,12 +910,12 @@ bool filterObservations3D(SfMData& sfmData, const FilterParams::FilterObservatio for(auto i = 0; i < landmarksData.size(); i++) { sfmData::Landmark& landmark = *landmarksData[i]; - const auto& nbObservations = landmark.observations.size(); auto& [viewIds, viewScores] = viewScoresData[i]; + const auto& nbObservations = viewIds.size(); - // check number of observations - if(landmark.observations.size() <= params.maxNbObservationsPerLandmark) - continue; + //// check number of observations + //if(landmark.observations.size() <= params.maxNbObservationsPerLandmark) + // continue; // sort by descending view score order std::vector idx(nbObservations); @@ -854,16 +927,22 @@ bool filterObservations3D(SfMData& sfmData, const FilterParams::FilterObservatio //double threshold = 0.1 / params.maxNbObservationsPerLandmark; for(auto j = 0; j < params.maxNbObservationsPerLandmark; j++) { - const auto& viewScore = viewScores[idx[j]]; + //const auto& viewScore = viewScores[idx[j]]; /*if(viewScore < threshold) break;*/ + + // add observation only if it's an original observation and not augmented const auto& viewId = viewIds[idx[j]]; - filteredObservations[viewId] = landmark.observations[viewId]; + const auto& obsIt = landmark.observations.find(viewId); + if(obsIt != landmark.observations.end()) + filteredObservations[viewId] = landmark.observations[viewId]; } landmark.observations = std::move(filteredObservations); } ALICEVISION_LOG_INFO("Selecting observations with best scores: done"); + removeNonObservedLandmarks(sfmData); + return true; } @@ -875,7 +954,6 @@ double filter2DView(SfMData& sfmData, const FilterParams::FilterObservations2DPa KdTree2D tree(2, data, nanoflann::KDTreeSingleIndexAdaptorParams(MAX_LEAF_ELEMENTS)); tree.buildIndex(); - // note that the observation is a neighbor to itself with zero distance, hence the +/- 1 size_t nbNeighbors_ = params.nbNeighbors2D; // average neighbors distance for each observation std::vector means(observations.size(), std::numeric_limits::max()); @@ -988,21 +1066,14 @@ bool filterObservations2D(SfMData& sfmData, const FilterParams::FilterObservatio } } } - // remove landmarks with no remaining observations - for(auto it = sfmData.getLandmarks().begin(); it != sfmData.getLandmarks().end();) - { - if(it->second.observations.size() == 0) - it = sfmData.getLandmarks().erase(it); - else - ++it; - } + + removeNonObservedLandmarks(sfmData); return true; } int aliceVision_main(int argc, char *argv[]) { - std::cout << "hello" << std::endl; // command-line parameters std::string inputSfmFilename; From 0278a7b964af28486a9c03499aa19d05477e731a Mon Sep 17 00:00:00 2001 From: Abdelrahman AL MAROUK Date: Thu, 12 Oct 2023 10:40:56 +0200 Subject: [PATCH 33/37] [WIP][filterSfM] multiple improvements - add option to propagate neighbor observations - use Jaccard index in step 2 --- src/software/pipeline/main_filterSfM.cpp | 110 +++++++++++------------ 1 file changed, 52 insertions(+), 58 deletions(-) diff --git a/src/software/pipeline/main_filterSfM.cpp b/src/software/pipeline/main_filterSfM.cpp index 820f76657e..0d8a7405f4 100644 --- a/src/software/pipeline/main_filterSfM.cpp +++ b/src/software/pipeline/main_filterSfM.cpp @@ -61,7 +61,7 @@ struct FilterParams struct FilterLandmarksStep2Params { bool enabled = true; - int maxNbObservationsPerLandmark = 2; + float minScore = .3f; int nbNeighbors3D = 10; } step2; struct FilterLandmarksStep3Params @@ -77,6 +77,9 @@ struct FilterParams int maxNbObservationsPerLandmark = 2; bool propagationEnabled = true; int nbNeighbors3D = 10; + bool observationsPropagationEnabled = true; + int observationsPropagationFrequency = 1; + int observationsPropagationCount = 1; double neighborsInfluence = 0.5; int nbIterations = 5; bool dampingEnabled = true; @@ -97,7 +100,7 @@ std::istream& operator>>(std::istream& in, FilterParams& params) in >> token; std::vector splitParams; boost::split(splitParams, token, boost::algorithm::is_any_of(":")); - if(splitParams.size() != 21) + if(splitParams.size() != 24) throw std::invalid_argument("Failed to parse FilterParams from: \n" + token); int i = 0; @@ -107,7 +110,7 @@ std::istream& operator>>(std::istream& in, FilterParams& params) params.filterLandmarks.step1.minNbObservationsPerLandmark = boost::lexical_cast(splitParams[i++]); params.filterLandmarks.step2.enabled = boost::to_lower_copy(splitParams[i++]) == "true"; - params.filterLandmarks.step2.maxNbObservationsPerLandmark = boost::lexical_cast(splitParams[i++]); + params.filterLandmarks.step2.minScore = boost::lexical_cast(splitParams[i++]); params.filterLandmarks.step2.nbNeighbors3D = boost::lexical_cast(splitParams[i++]); params.filterLandmarks.step3.enabled = boost::to_lower_copy(splitParams[i++]) == "true"; @@ -118,6 +121,9 @@ std::istream& operator>>(std::istream& in, FilterParams& params) params.filterObservations3D.maxNbObservationsPerLandmark = boost::lexical_cast(splitParams[i++]); params.filterObservations3D.propagationEnabled = boost::to_lower_copy(splitParams[i++]) == "true"; params.filterObservations3D.nbNeighbors3D = boost::lexical_cast(splitParams[i++]); + params.filterObservations3D.observationsPropagationEnabled = boost::to_lower_copy(splitParams[i++]) == "true"; + params.filterObservations3D.observationsPropagationFrequency = boost::lexical_cast(splitParams[i++]); + params.filterObservations3D.observationsPropagationCount = boost::lexical_cast(splitParams[i++]); params.filterObservations3D.neighborsInfluence = boost::lexical_cast(splitParams[i++]); params.filterObservations3D.nbIterations = boost::lexical_cast(splitParams[i++]); params.filterObservations3D.dampingEnabled = boost::to_lower_copy(splitParams[i++]) == "true"; @@ -436,40 +442,24 @@ void filterLandmarks_step2(SfMData& sfmData, #pragma omp parallel for reduction(+ : nbToRemove) for(auto i = 0; i < landmarksData.size(); i++) { - int score = 0; + float score = 0.f; auto& [viewIds, neighbors] = viewData[i]; for(auto j = 0; j < neighbors.size(); j++) { int scoreNeighbor = 0; const auto& neighborId = neighbors[j]; const auto& viewIds_neighbor = viewData[neighborId].first; - // loop over common views - auto viewIds_it = viewIds.begin(); - auto viewIds_neighbor_it = viewIds_neighbor.begin(); - while(viewIds_it != viewIds.end() && viewIds_neighbor_it != viewIds_neighbor.end()) - { - if(*viewIds_it < *viewIds_neighbor_it) - { - ++viewIds_it; - } - else - { - // if same view - if(!(*viewIds_neighbor_it < *viewIds_it)) - { - scoreNeighbor++; - if(scoreNeighbor == params.maxNbObservationsPerLandmark) - { - score++; - break; - } - ++viewIds_it; - } - ++viewIds_neighbor_it; - } - } + + std::vector viewIds_intersection; + std::set_intersection(viewIds.begin(), viewIds.end(), viewIds_neighbor.begin(), viewIds_neighbor.end(), + std::back_inserter(viewIds_intersection)); + float nbIntersection = static_cast(viewIds_intersection.size()); + // Jaccard index or Intersection over Union (IoU): #overlap / #union + score += nbIntersection / ((viewIds.size() + viewIds_neighbor.size()) - nbIntersection); } - if(score < params.nbNeighbors3D / 2) + //score /= neighbors.size(); + score /= params.nbNeighbors3D; + if(score < params.minScore) { toRemove[i] = true; nbToRemove++; @@ -686,10 +676,10 @@ void computeNeighborsInfo(std::vector& landmarksData, const FilterPar ALICEVISION_LOG_INFO("Computing landmark neighbors and distance-based weights: done"); } -void augmentInitialObservations(std::vector, std::vector>>& viewScoresData, +void propagateNeighborsObservations(std::vector, std::vector>>& viewScoresData, std::vector, std::vector>>& neighborsData) { - ALICEVISION_LOG_INFO("Augment initial observations based on neighbors observations: started"); + ALICEVISION_LOG_INFO("Augment observations with neighbors observations: started"); std::vector> viewScoresData_t(viewScoresData.size()); // for each landmark, identify the new observations coming from its neighborhood #pragma omp parallel for @@ -712,12 +702,14 @@ void augmentInitialObservations(std::vector, std:: std::back_inserter(viewIds_toAdd)); viewIds_t = std::move(viewIds_toAdd); } + double nbAdded = 0; // for each landmark, add the previously identified observations - #pragma omp parallel for +#pragma omp parallel for reduction(+ : nbAdded) for(int id = 0; id < viewScoresData.size(); id++) { auto& [viewIds, viewScores] = viewScoresData[id]; const auto& viewIds_toAdd = viewScoresData_t[id]; + nbAdded += static_cast(viewIds_toAdd.size()) / viewIds.size(); viewIds.insert(viewIds.end(), viewIds_toAdd.begin(), viewIds_toAdd.end()); viewScores.insert(viewScores.end(), viewIds_toAdd.size(), 0.); @@ -737,21 +729,24 @@ void augmentInitialObservations(std::vector, std:: viewScores[j] = scores_temp[idx[j]]; } } - ALICEVISION_LOG_INFO("Augment initial observations based on neighbors observations: done"); + ALICEVISION_LOG_INFO( + "On average, observations per landmark are augmented by: " << nbAdded * 100. / viewScoresData.size() << " %."); + ALICEVISION_LOG_INFO("Augment observations with neighbors observations: done"); } -void computeNewScores(const std::vector& landmarksData, +void propagateNeighborsScores(const std::vector& landmarksData, const FilterParams::FilterObservations3DParams& params, const std::vector, std::vector>>& neighborsData, const std::vector, std::vector>>& viewScoresData, std::vector>& viewScoresData_t) { + ALICEVISION_LOG_INFO("Propagating neighbors observation scores: started"); #pragma omp parallel for for(auto id = 0; id < landmarksData.size(); id++) { const auto& [viewIds, viewScores] = viewScoresData[id]; auto& viewScores_acc = viewScoresData_t[id]; - // initialize to zero, will first contain the weighted average scores of neighbors + // initialize to zero, will contain the weighted average scores of neighbors viewScores_acc.assign(viewScores.size(), 0.); // accumulator for normalisation double viewScores_total = 0.; @@ -789,19 +784,16 @@ void computeNewScores(const std::vector& landmarksData, } } } - // if common views with neighbor landmarks - //if(viewScores_total != 0.) for(auto j = 0; j < viewScores_acc.size(); j++) { - // normalize score and apply influence factor - //viewScores_acc[j] *= params.neighborsInfluence / viewScores_total; + // apply influence factor viewScores_acc[j] *= params.neighborsInfluence; // combine weighted neighbor scores and the landmark's own scores viewScores_acc[j] += (1 - params.neighborsInfluence) * viewScores[j]; } // dampen scores of non-chosen observations - if(params.dampingEnabled && viewScores_acc.size() <= params.maxNbObservationsPerLandmark) + if(params.dampingEnabled && viewScores_acc.size() > params.maxNbObservationsPerLandmark) { // sort by descending view score order std::vector idx(viewScores_acc.size()); @@ -820,6 +812,7 @@ void computeNewScores(const std::vector& landmarksData, viewScores_acc[j] /= viewScores_total; } } + ALICEVISION_LOG_INFO("Propagating neighbors observation scores: done"); } void propagateNeighborsInfo(std::vector& landmarksData, @@ -827,12 +820,19 @@ void propagateNeighborsInfo(std::vector& landmarksData, std::vector, std::vector>>& neighborsData, std::vector, std::vector>>& viewScoresData) { - ALICEVISION_LOG_INFO("Propagating neighbors observation scores: started"); // new view scores at iteration t std::vector> viewScoresData_t(landmarksData.size()); + int observationsPropagationCount = + params.observationsPropagationCount == 0 ? params.nbIterations : params.observationsPropagationCount; for(auto i = 0; i < params.nbIterations; i++) { - computeNewScores(landmarksData, params, neighborsData, viewScoresData, viewScoresData_t); + if(params.observationsPropagationEnabled && (i % params.observationsPropagationFrequency == 0) && + observationsPropagationCount > 0) + { + propagateNeighborsObservations(viewScoresData, neighborsData); + --observationsPropagationCount; + } + propagateNeighborsScores(landmarksData, params, neighborsData, viewScoresData, viewScoresData_t); // Mean Squared Error double error = 0.; @@ -858,7 +858,6 @@ void propagateNeighborsInfo(std::vector& landmarksData, error /= landmarksData.size(); ALICEVISION_LOG_INFO("MSE at iteration " << i << ": " << error); } - ALICEVISION_LOG_INFO("Propagating neighbors observation scores: done"); } void removeNonObservedLandmarks(SfMData& sfmData) @@ -901,7 +900,6 @@ bool filterObservations3D(SfMData& sfmData, const FilterParams::FilterObservatio // contains the neighbor landmarks ids with their corresponding weights std::vector, std::vector>> neighborsData(landmarksData.size()); computeNeighborsInfo(landmarksData, params, neighborsData); - augmentInitialObservations(viewScoresData, neighborsData); propagateNeighborsInfo(landmarksData, params, neighborsData, viewScoresData); } @@ -913,10 +911,6 @@ bool filterObservations3D(SfMData& sfmData, const FilterParams::FilterObservatio auto& [viewIds, viewScores] = viewScoresData[i]; const auto& nbObservations = viewIds.size(); - //// check number of observations - //if(landmark.observations.size() <= params.maxNbObservationsPerLandmark) - // continue; - // sort by descending view score order std::vector idx(nbObservations); std::iota(idx.begin(), idx.end(), 0); @@ -924,13 +918,10 @@ bool filterObservations3D(SfMData& sfmData, const FilterParams::FilterObservatio // keep only observations with best scores Observations filteredObservations; - //double threshold = 0.1 / params.maxNbObservationsPerLandmark; - for(auto j = 0; j < params.maxNbObservationsPerLandmark; j++) + size_t maxNbObservationsPerLandmark = + std::min(static_cast(params.maxNbObservationsPerLandmark), landmark.observations.size()); + for(auto j = 0; j < maxNbObservationsPerLandmark; j++) { - //const auto& viewScore = viewScores[idx[j]]; - /*if(viewScore < threshold) - break;*/ - // add observation only if it's an original observation and not augmented const auto& viewId = viewIds[idx[j]]; const auto& obsIt = landmark.observations.find(viewId); @@ -954,7 +945,7 @@ double filter2DView(SfMData& sfmData, const FilterParams::FilterObservations2DPa KdTree2D tree(2, data, nanoflann::KDTreeSingleIndexAdaptorParams(MAX_LEAF_ELEMENTS)); tree.buildIndex(); - size_t nbNeighbors_ = params.nbNeighbors2D; + int nbNeighbors_ = params.nbNeighbors2D; // average neighbors distance for each observation std::vector means(observations.size(), std::numeric_limits::max()); const std::size_t cacheSize = 1000; @@ -1106,7 +1097,7 @@ int aliceVision_main(int argc, char *argv[]) " * Min Nb of Observations: Minimum number of observations required to keep a landmark.\n" " * Step 2 parameters:\n" " * Enabled: Remove landmarks with dissimilar observations of 3D landmark neighbors.\n" - " * Max Nb of Observations: Maximum number of allowed observations per landmark.\n" + " * Min Score: Minimum similarity score between the observations of the landmark and of its neighbors.\n" " * Nb of Neighbors in 3D: Number of neighbor landmarks used in making the decision for best observations.\n" " * Step 3 parameters:\n" " * Enabled: Remove landmarks with worse resolution than neighbors.\n" @@ -1116,8 +1107,11 @@ int aliceVision_main(int argc, char *argv[]) "Filter Observations 3D parameters:\n" " * Enabled: Select best observations for observation consistency between 3D neighboring landmarks.\n" " * Max Nb of Observations: Maximum number of allowed observations per landmark.\n" - " * Enable Neighbors Influence: nable propagating neighbors' scores iteratively.\n" + " * Enable Neighbors Influence: Enable propagating neighbors' info (scores, observations) iteratively.\n" " * Nb of Neighbors in 3D: Number of neighbor landmarks used in making the decision for best observations.\n" + " * Enable Propagating Observations: Propagate neighbors observations before propagating the scores.\n" + " * Frequency: Specifies every how many iterations should the observations be propagated.\n" + " * Count: Maximum number of times the observations are propagated (0 for no limit).\n" " * Neighbors Influence: Specifies how much influential the neighbors are in selecting the best observations.\n" " Between 0. and 1., the closer to 1., the more influencial the neighborhood is.\n" " * Nb of Iterations: Number of iterations to propagate neighbors information.\n" From 7d394fb5953dea98daba4f88b7374e25002d2205 Mon Sep 17 00:00:00 2001 From: Abdelrahman AL MAROUK Date: Mon, 16 Oct 2023 11:45:17 +0200 Subject: [PATCH 34/37] [filterSfM] fix bug in damping --- src/software/pipeline/main_filterSfM.cpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/software/pipeline/main_filterSfM.cpp b/src/software/pipeline/main_filterSfM.cpp index 0d8a7405f4..b6aa878f5d 100644 --- a/src/software/pipeline/main_filterSfM.cpp +++ b/src/software/pipeline/main_filterSfM.cpp @@ -803,9 +803,9 @@ void propagateNeighborsScores(const std::vector& landmarksData, viewScores_total = 1.; for(auto j = params.maxNbObservationsPerLandmark; j < viewScores_acc.size(); j++) { - const double& v = params.dampingFactor * viewScores_acc[j]; - viewScores_total += v - viewScores_acc[j]; - viewScores_acc[j] = v; + const double& v = params.dampingFactor * viewScores_acc[idx[j]]; + viewScores_total += v - viewScores_acc[idx[j]]; + viewScores_acc[idx[j]] = v; } // re-normalize for(auto j = 0; j < viewScores_acc.size(); j++) From 98cea2f70eea5e0a44e48f0cc99709264e69198d Mon Sep 17 00:00:00 2001 From: Abdelrahman AL MAROUK Date: Mon, 16 Oct 2023 15:42:05 +0200 Subject: [PATCH 35/37] [filterSfM] add option to keep propagated observations at the end of 3D Observations filter --- src/software/pipeline/main_filterSfM.cpp | 14 +++++++++++--- 1 file changed, 11 insertions(+), 3 deletions(-) diff --git a/src/software/pipeline/main_filterSfM.cpp b/src/software/pipeline/main_filterSfM.cpp index b6aa878f5d..8252f3f9d4 100644 --- a/src/software/pipeline/main_filterSfM.cpp +++ b/src/software/pipeline/main_filterSfM.cpp @@ -80,6 +80,7 @@ struct FilterParams bool observationsPropagationEnabled = true; int observationsPropagationFrequency = 1; int observationsPropagationCount = 1; + bool observationsPropagationKeep = false; double neighborsInfluence = 0.5; int nbIterations = 5; bool dampingEnabled = true; @@ -100,7 +101,7 @@ std::istream& operator>>(std::istream& in, FilterParams& params) in >> token; std::vector splitParams; boost::split(splitParams, token, boost::algorithm::is_any_of(":")); - if(splitParams.size() != 24) + if(splitParams.size() != 25) throw std::invalid_argument("Failed to parse FilterParams from: \n" + token); int i = 0; @@ -124,6 +125,7 @@ std::istream& operator>>(std::istream& in, FilterParams& params) params.filterObservations3D.observationsPropagationEnabled = boost::to_lower_copy(splitParams[i++]) == "true"; params.filterObservations3D.observationsPropagationFrequency = boost::lexical_cast(splitParams[i++]); params.filterObservations3D.observationsPropagationCount = boost::lexical_cast(splitParams[i++]); + params.filterObservations3D.observationsPropagationKeep = boost::to_lower_copy(splitParams[i++]) == "true"; params.filterObservations3D.neighborsInfluence = boost::lexical_cast(splitParams[i++]); params.filterObservations3D.nbIterations = boost::lexical_cast(splitParams[i++]); params.filterObservations3D.dampingEnabled = boost::to_lower_copy(splitParams[i++]) == "true"; @@ -924,9 +926,14 @@ bool filterObservations3D(SfMData& sfmData, const FilterParams::FilterObservatio { // add observation only if it's an original observation and not augmented const auto& viewId = viewIds[idx[j]]; - const auto& obsIt = landmark.observations.find(viewId); - if(obsIt != landmark.observations.end()) + if(params.observationsPropagationKeep) filteredObservations[viewId] = landmark.observations[viewId]; + else + { + const auto& obsIt = landmark.observations.find(viewId); + if(obsIt != landmark.observations.end()) + filteredObservations[viewId] = landmark.observations[viewId]; + } } landmark.observations = std::move(filteredObservations); } @@ -1112,6 +1119,7 @@ int aliceVision_main(int argc, char *argv[]) " * Enable Propagating Observations: Propagate neighbors observations before propagating the scores.\n" " * Frequency: Specifies every how many iterations should the observations be propagated.\n" " * Count: Maximum number of times the observations are propagated (0 for no limit).\n" + " * Keep Observations: Specifies if propagated observations are to be kept at the end.\n" " * Neighbors Influence: Specifies how much influential the neighbors are in selecting the best observations.\n" " Between 0. and 1., the closer to 1., the more influencial the neighborhood is.\n" " * Nb of Iterations: Number of iterations to propagate neighbors information.\n" From d76e1660f6f9757b581af6fcb0f56c51bdecb6c9 Mon Sep 17 00:00:00 2001 From: Abdelrahman AL MAROUK Date: Mon, 16 Oct 2023 16:48:57 +0200 Subject: [PATCH 36/37] [WIP][filterSfM] fix option to keep propagated observations - create a fake observation by projecting the landmark in the view --- src/software/pipeline/main_filterSfM.cpp | 16 ++++++++++------ 1 file changed, 10 insertions(+), 6 deletions(-) diff --git a/src/software/pipeline/main_filterSfM.cpp b/src/software/pipeline/main_filterSfM.cpp index 8252f3f9d4..be88b0bc93 100644 --- a/src/software/pipeline/main_filterSfM.cpp +++ b/src/software/pipeline/main_filterSfM.cpp @@ -926,13 +926,17 @@ bool filterObservations3D(SfMData& sfmData, const FilterParams::FilterObservatio { // add observation only if it's an original observation and not augmented const auto& viewId = viewIds[idx[j]]; - if(params.observationsPropagationKeep) - filteredObservations[viewId] = landmark.observations[viewId]; - else + const auto& obsIt = landmark.observations.find(viewId); + if(obsIt != landmark.observations.end()) + filteredObservations[viewId] = obsIt->second; + else if (params.observationsPropagationKeep) { - const auto& obsIt = landmark.observations.find(viewId); - if(obsIt != landmark.observations.end()) - filteredObservations[viewId] = landmark.observations[viewId]; + // project landmark in view to find observation coords + const sfmData::View* view = sfmData.getViews().at(viewId).get(); + const geometry::Pose3 pose = sfmData.getPose(*view).getTransform(); + const camera::IntrinsicBase* intrinsic = sfmData.getIntrinsics().at(view->getIntrinsicId()).get(); + const Vec2& x = intrinsic->project(pose, landmark.X.homogeneous()); + filteredObservations[viewId] = Observation(x, UndefinedIndexT, 0.0); } } landmark.observations = std::move(filteredObservations); From 27a1b6442b51941cfea051cc90882980dd956fd9 Mon Sep 17 00:00:00 2001 From: Abdelrahman AL MAROUK Date: Mon, 18 Mar 2024 15:27:27 +0100 Subject: [PATCH 37/37] [filterSfM] update code to follow latest merges --- src/software/pipeline/main_filterSfM.cpp | 43 +++++++++---------- .../pipeline/main_prepareDenseScene.cpp | 18 ++++---- 2 files changed, 31 insertions(+), 30 deletions(-) diff --git a/src/software/pipeline/main_filterSfM.cpp b/src/software/pipeline/main_filterSfM.cpp index be88b0bc93..230a0809c0 100644 --- a/src/software/pipeline/main_filterSfM.cpp +++ b/src/software/pipeline/main_filterSfM.cpp @@ -165,7 +165,7 @@ struct ObservationsAdaptator inline size_t kdtree_get_point_count() const { return _data.size(); } // Returns the dim'th component of the idx'th point in the class: - inline T kdtree_get_pt(const size_t idx, int dim) const { return _data[idx].x(dim); } + inline T kdtree_get_pt(const size_t idx, int dim) const { return _data[idx].getCoordinates()(dim); } // Optional bounding-box computation: return false to default to a standard bbox computation loop. // Return true if the BBOX was already computed by the class and returned in "bb" so it can be avoided to redo it @@ -353,7 +353,7 @@ ObservationsPerView getObservationsPerViews(SfMData& sfmData) ObservationsPerView observationsPerView; for(auto& landIt : sfmData.getLandmarks()) { - for(const auto& obsIt : landIt.second.observations) + for(const auto& obsIt : landIt.second.getObservations()) { IndexT viewId = obsIt.first; auto& landmarksSet = observationsPerView[viewId]; @@ -376,7 +376,7 @@ void filterLandmarks_step1(SfMData& sfmData, ALICEVISION_LOG_INFO("Removing landmarks having an insufficient number of observations: started."); for(auto& it : sfmData.getLandmarks()) { - if(it.second.observations.size() < params.minNbObservationsPerLandmark) + if (it.second.getObservations().size() < params.minNbObservationsPerLandmark) continue; landmarksData[i++] = it.second; } @@ -415,10 +415,10 @@ void filterLandmarks_step2(SfMData& sfmData, for(auto i = 0; i < landmarksData.size(); i++) { const sfmData::Landmark& landmark = landmarksData[i]; - const auto& nbObservations = landmark.observations.size(); + const auto& nbObservations = landmark.getObservations().size(); auto& [viewIds, neighbors] = viewData[i]; viewIds.reserve(nbObservations); - for(const auto& observationPair : landmark.observations) + for (const auto& observationPair : landmark.getObservations()) { const IndexT viewId = observationPair.first; viewIds.push_back(viewId); @@ -504,12 +504,12 @@ void filterLandmarks_step3(SfMData& sfmData, // compute landmark pixSize double pixSize = 0.; int n = 0; - for(const auto& observationPair : landmark.observations) + for (const auto& observationPair : landmark.getObservations()) { const IndexT viewId = observationPair.first; pixSize += mp.getCamPixelSize(Point3d(landmark.X.x(), landmark.X.y(), landmark.X.z()), mp.getIndexFromViewId(viewId), - params.useFeatureScale ? observationPair.second.scale : 1); + params.useFeatureScale ? observationPair.second.getScale() : 1); n++; } pixSize /= n; @@ -594,13 +594,13 @@ void computeInitialObservationScores(SfMData& sfmData, std::vector& l const sfmData::Landmark& landmark = *landmarksData[i]; // compute observation scores - const auto& nbObservations = landmark.observations.size(); + const auto& nbObservations = landmark.getObservations().size(); auto& [viewIds, viewScores] = viewScoresData[i]; viewIds.reserve(nbObservations); viewScores.reserve(nbObservations); // accumulator for normalizing the scores double total = 0.; - for(const auto& observationPair : landmark.observations) + for (const auto& observationPair : landmark.getObservations()) { const IndexT viewId = observationPair.first; const sfmData::View& view = *(sfmData.getViews().at(viewId)); @@ -869,7 +869,7 @@ void removeNonObservedLandmarks(SfMData& sfmData) const auto& initialNbLandmarks = sfmData.getLandmarks().size(); for(auto it = sfmData.getLandmarks().begin(); it != sfmData.getLandmarks().end();) { - if(it->second.observations.size() == 0) + if (it->second.getObservations().size() == 0) it = sfmData.getLandmarks().erase(it); else ++it; @@ -920,14 +920,13 @@ bool filterObservations3D(SfMData& sfmData, const FilterParams::FilterObservatio // keep only observations with best scores Observations filteredObservations; - size_t maxNbObservationsPerLandmark = - std::min(static_cast(params.maxNbObservationsPerLandmark), landmark.observations.size()); + size_t maxNbObservationsPerLandmark = std::min(static_cast(params.maxNbObservationsPerLandmark), landmark.getObservations().size()); for(auto j = 0; j < maxNbObservationsPerLandmark; j++) { // add observation only if it's an original observation and not augmented const auto& viewId = viewIds[idx[j]]; - const auto& obsIt = landmark.observations.find(viewId); - if(obsIt != landmark.observations.end()) + const auto& obsIt = landmark.getObservations().find(viewId); + if (obsIt != landmark.getObservations().end()) filteredObservations[viewId] = obsIt->second; else if (params.observationsPropagationKeep) { @@ -939,7 +938,7 @@ bool filterObservations3D(SfMData& sfmData, const FilterParams::FilterObservatio filteredObservations[viewId] = Observation(x, UndefinedIndexT, 0.0); } } - landmark.observations = std::move(filteredObservations); + landmark.getObservations() = std::move(filteredObservations); } ALICEVISION_LOG_INFO("Selecting observations with best scores: done"); @@ -970,7 +969,7 @@ double filter2DView(SfMData& sfmData, const FilterParams::FilterObservations2DPa std::vector distances_(nbNeighbors_); KnnNonZeroSearch resultSet(nbNeighbors_); resultSet.init(&indices_[0], &distances_[0]); - tree.findNeighbors(resultSet, obs.x.data()); + tree.findNeighbors(resultSet, obs.getCoordinates().data()); const auto& nbFound = resultSet.size(); if(nbFound == 0) continue; @@ -1023,7 +1022,7 @@ double filter2DView(SfMData& sfmData, const FilterParams::FilterObservations2DPa } bool filterObservations2D(SfMData& sfmData, const FilterParams::FilterObservations2DParams& params, - HashMap& estimatedRadii) + std::map& estimatedRadii) { std::set viewIds = sfmData.getValidViews(); std::vector estimatedRadii_(viewIds.size(), -1.); @@ -1046,7 +1045,7 @@ bool filterObservations2D(SfMData& sfmData, const FilterParams::FilterObservatio // clear and update landmark observations for(auto& landmark : sfmData.getLandmarks()) { - landmark.second.observations.clear(); + landmark.second.getObservations().clear(); } for(int i = 0; i < viewIds.size(); ++i) { @@ -1064,7 +1063,7 @@ bool filterObservations2D(SfMData& sfmData, const FilterParams::FilterObservatio auto& landmarks = observationsIt->second.second; for(int j = 0; j < observations.size(); j++) { - landmarks[j]->observations[viewId] = observations[j]; + landmarks[j]->getObservations()[viewId] = observations[j]; } } } @@ -1163,7 +1162,7 @@ int aliceVision_main(int argc, char *argv[]) // Read the input SfM scene SfMData sfmData; - if(!sfmDataIO::Load(sfmData, inputSfmFilename, sfmDataIO::ESfMData::ALL)) + if(!sfmDataIO::load(sfmData, inputSfmFilename, sfmDataIO::ESfMData::ALL)) { ALICEVISION_LOG_ERROR("The input SfMData file '" << inputSfmFilename << "' cannot be read."); return EXIT_FAILURE; @@ -1188,7 +1187,7 @@ int aliceVision_main(int argc, char *argv[]) if(params.filterObservations2D.enabled) { - HashMap estimatedRadii; + std::map estimatedRadii; ALICEVISION_LOG_INFO("Filtering observations in 2D: started."); filterObservations2D(sfmData, params.filterObservations2D, estimatedRadii); ALICEVISION_LOG_INFO("Filtering observations in 2D: done."); @@ -1206,7 +1205,7 @@ int aliceVision_main(int argc, char *argv[]) } } - sfmDataIO::Save(sfmData, outputSfmFilename, sfmDataIO::ESfMData::ALL); + sfmDataIO::save(sfmData, outputSfmFilename, sfmDataIO::ESfMData::ALL); return EXIT_SUCCESS; } diff --git a/src/software/pipeline/main_prepareDenseScene.cpp b/src/software/pipeline/main_prepareDenseScene.cpp index 2e1daed7e2..8a735468c6 100644 --- a/src/software/pipeline/main_prepareDenseScene.cpp +++ b/src/software/pipeline/main_prepareDenseScene.cpp @@ -95,7 +95,7 @@ ObservationsPerView getObservationsPerViews(const SfMData& sfmData) ObservationsPerView observationsPerView; for(auto& landIt : sfmData.getLandmarks()) { - for(const auto& obsIt : landIt.second.observations) + for (const auto& obsIt : landIt.second.getObservations()) { IndexT viewId = obsIt.first; auto& observationsSet = observationsPerView[viewId]; @@ -155,7 +155,7 @@ bool prepareDenseScene(const SfMData& sfmData, bool doMaskLandmarks = landmarksMaskScale > 0.f; ObservationsPerView observationsPerView; - HashMap estimatedRadii; + std::map estimatedRadii; if (doMaskLandmarks) { observationsPerView = std::move(getObservationsPerViews(sfmData)); @@ -334,11 +334,13 @@ bool prepareDenseScene(const SfMData& sfmData, for(const auto& observation : observations) { const auto& obs = *observation; - for(int y = std::max(obs.x.y() - r, 0.); - y <= std::min(obs.x.y() + r, (double)maskLandmarks.Height() - 1); y++) + for (int y = std::max(obs.getCoordinates().y() - r, 0.); + y <= std::min(obs.getCoordinates().y() + r, (double)maskLandmarks.height() - 1); + y++) { - for(int x = std::max(obs.x.x() - r, 0.); - x <= std::min(obs.x.x() + r, (double)maskLandmarks.Width() - 1); x++) + for (int x = std::max(obs.getCoordinates().x() - r, 0.); + x <= std::min(obs.getCoordinates().x() + r, (double)maskLandmarks.width() - 1); + x++) { maskLandmarks(y, x) = std::numeric_limits::max(); } @@ -356,13 +358,13 @@ bool prepareDenseScene(const SfMData& sfmData, dstColorImage, cam, metadata, srcImage, evCorrection, exposureCompensation, [&maskLoaded, &mask, &maskLandmarks, &doMaskLandmarks](Image& image) { - if(maskLoaded && (image.Width() * image.Height() != mask.Width() * mask.Height())) + if(maskLoaded && (image.width() * image.height() != mask.width() * mask.height())) { ALICEVISION_LOG_WARNING("Invalid image mask size: mask is ignored."); return; } - for(int pix = 0; pix < image.Width() * image.Height(); ++pix) + for(int pix = 0; pix < image.width() * image.height(); ++pix) { image(pix).a() = (maskLoaded && mask(pix) == 0) ? 0.f : (doMaskLandmarks && maskLandmarks(pix) == 127) ? .5f : 1.f; }