From 118dc6469e15530af5b6d642b4d73c8795a3cbc2 Mon Sep 17 00:00:00 2001 From: Fabien Servant Date: Sat, 20 Sep 2025 12:53:23 +0200 Subject: [PATCH 1/2] remove use of eigen large headers globally --- .../calibration/checkerDetector.cpp | 1 + src/aliceVision/feature/akaze/AKAZE.cpp | 2 + src/aliceVision/fuseCut/Mesher.cpp | 1 + src/aliceVision/geometry/Pose3.hpp | 1 + .../lightingEstimation/ellipseGeometry.cpp | 1 + .../lightingCalibration.cpp | 1 + .../lightingEstimation/lightingEstimation.cpp | 1 + src/aliceVision/matching/svgVisualization.cpp | 2 + .../geometricFilterUtils.cpp | 2 + src/aliceVision/multiview/affineSolver.cpp | 2 + src/aliceVision/multiview/essential.cpp | 2 + .../relativePose/Essential5PSolver.cpp | 4 + .../relativePose/Essential8PSolver.cpp | 1 + .../relativePose/Fundamental10PSolver.cpp | 2 + .../relativePose/Fundamental7PSolver.cpp | 3 + .../relativePose/Fundamental8PSolver.cpp | 1 + .../multiview/resection/EPnPSolver.cpp | 1 + .../multiview/resection/P4PfSolver.cpp | 4 + .../multiview/resection/P5PfrSolver.cpp | 5 + .../multiview/resection/Resection6PSolver.cpp | 3 + .../multiview/rotationAveraging/l2.cpp | 3 + src/aliceVision/numeric/CMakeLists.txt | 3 +- .../{numeric.cpp => NumericFunctions.cpp} | 5 +- src/aliceVision/numeric/NumericFunctions.hpp | 339 ++++++++++++++++++ src/aliceVision/numeric/numeric.hpp | 336 +---------------- src/aliceVision/numeric/projection.cpp | 1 + .../robustEstimation/LineKernel.hpp | 2 + .../ReconstructionEngine_panorama.cpp | 1 + src/aliceVision/sfm/utils/alignment.cpp | 2 + .../MeshSDFilter/MeshNormalFilter.h | 2 + .../pipeline/main_checkerboardCalibration.cpp | 2 + src/software/pipeline/main_panoramaInit.cpp | 1 + 32 files changed, 401 insertions(+), 336 deletions(-) rename src/aliceVision/numeric/{numeric.cpp => NumericFunctions.cpp} (97%) create mode 100644 src/aliceVision/numeric/NumericFunctions.hpp diff --git a/src/aliceVision/calibration/checkerDetector.cpp b/src/aliceVision/calibration/checkerDetector.cpp index ed62160fd8..fd4df54e35 100644 --- a/src/aliceVision/calibration/checkerDetector.cpp +++ b/src/aliceVision/calibration/checkerDetector.cpp @@ -19,6 +19,7 @@ #include #include +#include // TODO: to remove when moving to eigen 3.4 namespace Eigen { diff --git a/src/aliceVision/feature/akaze/AKAZE.cpp b/src/aliceVision/feature/akaze/AKAZE.cpp index 22e347e7e9..84e01ae430 100644 --- a/src/aliceVision/feature/akaze/AKAZE.cpp +++ b/src/aliceVision/feature/akaze/AKAZE.cpp @@ -13,6 +13,8 @@ #include #include +#include + namespace aliceVision { namespace feature { diff --git a/src/aliceVision/fuseCut/Mesher.cpp b/src/aliceVision/fuseCut/Mesher.cpp index 38d5a3745a..b2985eb243 100644 --- a/src/aliceVision/fuseCut/Mesher.cpp +++ b/src/aliceVision/fuseCut/Mesher.cpp @@ -9,6 +9,7 @@ #include #include #include +#include namespace aliceVision { namespace fuseCut { diff --git a/src/aliceVision/geometry/Pose3.hpp b/src/aliceVision/geometry/Pose3.hpp index 57f0e8753b..70d2ab14c0 100644 --- a/src/aliceVision/geometry/Pose3.hpp +++ b/src/aliceVision/geometry/Pose3.hpp @@ -10,6 +10,7 @@ #include #include #include +#include namespace aliceVision { namespace geometry { diff --git a/src/aliceVision/lightingEstimation/ellipseGeometry.cpp b/src/aliceVision/lightingEstimation/ellipseGeometry.cpp index 525294717f..71aaf48e80 100644 --- a/src/aliceVision/lightingEstimation/ellipseGeometry.cpp +++ b/src/aliceVision/lightingEstimation/ellipseGeometry.cpp @@ -13,6 +13,7 @@ #include #include +#include namespace aliceVision { namespace lightingEstimation { diff --git a/src/aliceVision/lightingEstimation/lightingCalibration.cpp b/src/aliceVision/lightingEstimation/lightingCalibration.cpp index 102c914aa0..9a5fc0e3b6 100644 --- a/src/aliceVision/lightingEstimation/lightingCalibration.cpp +++ b/src/aliceVision/lightingEstimation/lightingCalibration.cpp @@ -25,6 +25,7 @@ // Eigen #include #include +#include #include #include diff --git a/src/aliceVision/lightingEstimation/lightingEstimation.cpp b/src/aliceVision/lightingEstimation/lightingEstimation.cpp index 71eefd0e60..de9313d22d 100644 --- a/src/aliceVision/lightingEstimation/lightingEstimation.cpp +++ b/src/aliceVision/lightingEstimation/lightingEstimation.cpp @@ -10,6 +10,7 @@ #include #include +#include #include diff --git a/src/aliceVision/matching/svgVisualization.cpp b/src/aliceVision/matching/svgVisualization.cpp index d4b9a7d053..1e58c08630 100644 --- a/src/aliceVision/matching/svgVisualization.cpp +++ b/src/aliceVision/matching/svgVisualization.cpp @@ -9,6 +9,8 @@ #include #include +#include + #if ALICEVISION_IS_DEFINED(ALICEVISION_HAVE_CCTAG) #include #endif diff --git a/src/aliceVision/matchingImageCollection/geometricFilterUtils.cpp b/src/aliceVision/matchingImageCollection/geometricFilterUtils.cpp index 075bdc0952..bf59d3bd23 100644 --- a/src/aliceVision/matchingImageCollection/geometricFilterUtils.cpp +++ b/src/aliceVision/matchingImageCollection/geometricFilterUtils.cpp @@ -5,7 +5,9 @@ // You can obtain one at https://mozilla.org/MPL/2.0/. #include "geometricFilterUtils.hpp" + #include +#include namespace aliceVision { namespace matchingImageCollection { diff --git a/src/aliceVision/multiview/affineSolver.cpp b/src/aliceVision/multiview/affineSolver.cpp index c3e502f6e7..0a22040578 100644 --- a/src/aliceVision/multiview/affineSolver.cpp +++ b/src/aliceVision/multiview/affineSolver.cpp @@ -8,6 +8,8 @@ #include "affineSolver.hpp" +#include + namespace aliceVision { namespace multiview { diff --git a/src/aliceVision/multiview/essential.cpp b/src/aliceVision/multiview/essential.cpp index d4c11ff3c5..4d2e5cfb37 100644 --- a/src/aliceVision/multiview/essential.cpp +++ b/src/aliceVision/multiview/essential.cpp @@ -11,6 +11,8 @@ #include #include +#include + namespace aliceVision { // HZ 9.6 page 257 (formula 9.12) diff --git a/src/aliceVision/multiview/relativePose/Essential5PSolver.cpp b/src/aliceVision/multiview/relativePose/Essential5PSolver.cpp index f94703fdef..2f2caed32b 100644 --- a/src/aliceVision/multiview/relativePose/Essential5PSolver.cpp +++ b/src/aliceVision/multiview/relativePose/Essential5PSolver.cpp @@ -9,6 +9,10 @@ #include "Essential5PSolver.hpp" #include +#include +#include +#include + namespace aliceVision { namespace multiview { namespace relativePose { diff --git a/src/aliceVision/multiview/relativePose/Essential8PSolver.cpp b/src/aliceVision/multiview/relativePose/Essential8PSolver.cpp index ca06393ab9..97a6053d80 100644 --- a/src/aliceVision/multiview/relativePose/Essential8PSolver.cpp +++ b/src/aliceVision/multiview/relativePose/Essential8PSolver.cpp @@ -9,6 +9,7 @@ #include "Essential8PSolver.hpp" #include #include +#include namespace aliceVision { namespace multiview { diff --git a/src/aliceVision/multiview/relativePose/Fundamental10PSolver.cpp b/src/aliceVision/multiview/relativePose/Fundamental10PSolver.cpp index 2778a6be15..77f1ef307c 100644 --- a/src/aliceVision/multiview/relativePose/Fundamental10PSolver.cpp +++ b/src/aliceVision/multiview/relativePose/Fundamental10PSolver.cpp @@ -10,6 +10,8 @@ #include #include +#include + namespace aliceVision { namespace multiview { namespace relativePose { diff --git a/src/aliceVision/multiview/relativePose/Fundamental7PSolver.cpp b/src/aliceVision/multiview/relativePose/Fundamental7PSolver.cpp index 298c9d2437..68dcebb6e6 100644 --- a/src/aliceVision/multiview/relativePose/Fundamental7PSolver.cpp +++ b/src/aliceVision/multiview/relativePose/Fundamental7PSolver.cpp @@ -10,6 +10,9 @@ #include #include #include + +#include + namespace aliceVision { namespace multiview { namespace relativePose { diff --git a/src/aliceVision/multiview/relativePose/Fundamental8PSolver.cpp b/src/aliceVision/multiview/relativePose/Fundamental8PSolver.cpp index 658ec6a503..ca05196f5b 100644 --- a/src/aliceVision/multiview/relativePose/Fundamental8PSolver.cpp +++ b/src/aliceVision/multiview/relativePose/Fundamental8PSolver.cpp @@ -9,6 +9,7 @@ #include "Fundamental8PSolver.hpp" #include #include +#include namespace aliceVision { namespace multiview { diff --git a/src/aliceVision/multiview/resection/EPnPSolver.cpp b/src/aliceVision/multiview/resection/EPnPSolver.cpp index d97f52563d..b21085690c 100644 --- a/src/aliceVision/multiview/resection/EPnPSolver.cpp +++ b/src/aliceVision/multiview/resection/EPnPSolver.cpp @@ -7,6 +7,7 @@ #include "EPnPSolver.hpp" #include +#include namespace aliceVision { namespace multiview { diff --git a/src/aliceVision/multiview/resection/P4PfSolver.cpp b/src/aliceVision/multiview/resection/P4PfSolver.cpp index ea85cf951f..6c7a770520 100644 --- a/src/aliceVision/multiview/resection/P4PfSolver.cpp +++ b/src/aliceVision/multiview/resection/P4PfSolver.cpp @@ -5,9 +5,13 @@ // You can obtain one at https://mozilla.org/MPL/2.0/. #include "P4PfSolver.hpp" + #include #include +#include +#include + #include #include #include diff --git a/src/aliceVision/multiview/resection/P5PfrSolver.cpp b/src/aliceVision/multiview/resection/P5PfrSolver.cpp index 807b3e8cd8..e8ea6dcb78 100644 --- a/src/aliceVision/multiview/resection/P5PfrSolver.cpp +++ b/src/aliceVision/multiview/resection/P5PfrSolver.cpp @@ -5,11 +5,16 @@ // You can obtain one at https://mozilla.org/MPL/2.0/. #include "P5PfrSolver.hpp" + #include #include #include #include +#include +#include +#include + #include #include diff --git a/src/aliceVision/multiview/resection/Resection6PSolver.cpp b/src/aliceVision/multiview/resection/Resection6PSolver.cpp index 3207b0d432..030a65e8f9 100644 --- a/src/aliceVision/multiview/resection/Resection6PSolver.cpp +++ b/src/aliceVision/multiview/resection/Resection6PSolver.cpp @@ -6,9 +6,12 @@ // You can obtain one at https://mozilla.org/MPL/2.0/. #include "Resection6PSolver.hpp" + #include #include +#include + namespace aliceVision { namespace multiview { namespace resection { diff --git a/src/aliceVision/multiview/rotationAveraging/l2.cpp b/src/aliceVision/multiview/rotationAveraging/l2.cpp index e1d2849d27..f45b7a78e3 100644 --- a/src/aliceVision/multiview/rotationAveraging/l2.cpp +++ b/src/aliceVision/multiview/rotationAveraging/l2.cpp @@ -10,6 +10,9 @@ #include #include +#include +#include + #include #include diff --git a/src/aliceVision/numeric/CMakeLists.txt b/src/aliceVision/numeric/CMakeLists.txt index c9829e8c71..e6a81ddce0 100644 --- a/src/aliceVision/numeric/CMakeLists.txt +++ b/src/aliceVision/numeric/CMakeLists.txt @@ -1,6 +1,7 @@ # Headers set(numeric_files_headers numeric.hpp + NumericFunctions.hpp Container.hpp projection.hpp gps.hpp @@ -10,7 +11,7 @@ set(numeric_files_headers # Sources set(numeric_files_sources - numeric.cpp + NumericFunctions.cpp Container.cpp projection.cpp gps.cpp diff --git a/src/aliceVision/numeric/numeric.cpp b/src/aliceVision/numeric/NumericFunctions.cpp similarity index 97% rename from src/aliceVision/numeric/numeric.cpp rename to src/aliceVision/numeric/NumericFunctions.cpp index 6f7147f727..8e01a98e0e 100644 --- a/src/aliceVision/numeric/numeric.cpp +++ b/src/aliceVision/numeric/NumericFunctions.cpp @@ -1,12 +1,13 @@ // This file is part of the AliceVision project. -// Copyright (c) 2016 AliceVision contributors. +// Copyright (c) 2025 AliceVision contributors. // Copyright (c) 2012 openMVG contributors. // Copyright (c) 2007 libmv 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 "numeric.hpp" +#include +#include "NumericFunctions.hpp" #include #include diff --git a/src/aliceVision/numeric/NumericFunctions.hpp b/src/aliceVision/numeric/NumericFunctions.hpp new file mode 100644 index 0000000000..c087ee10b8 --- /dev/null +++ b/src/aliceVision/numeric/NumericFunctions.hpp @@ -0,0 +1,339 @@ +// This file is part of the AliceVision project. +// Copyright (c) 2025 AliceVision contributors. +// Copyright (c) 2012 openMVG contributors. +// Copyright (c) 2007 libmv 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/. + +#pragma once + +#include + +namespace aliceVision { + +/// Return the square of a number. + +template +inline T Square(T x) +{ + return x * x; +} + +/// Clamp return the number if inside range, else min or max range. + +template +inline T clamp(const T& val, const T& min, const T& max) +{ + return std::max(min, std::min(val, max)); + //(val < min) ? val : ((val>max) ? val : max); +} + +inline bool isSimilar(double a, double b) +{ + const double diff = a - b; + return std::abs(diff) < 1e-8; +} +inline bool isSimilar(float a, float b) +{ + const float diff = a - b; + return std::abs(diff) < 1e-8f; +} + +/** + * @brief Create a minimal skew matrix from a 2d vector. + * @param[in] x A 2d vector whose 3rd coordinate is supposed to be 1. + * @return The minimal ske matrix: [0, -1, x(1); 1, 0, -x(0);] + */ +Mat23 SkewMatMinimal(const Vec2& x); + +/** + * @brief Create a cross product matrix from a 3d vector. + * @param x A 3d vector. + * @return the cross matrix representation of the input vector. + */ +Mat3 CrossProductMatrix(const Vec3& x); + +// Create a rotation matrix around axis X with the provided radian angle +Mat3 RotationAroundX(double angle); + +// Create a rotation matrix around axis Y with the provided radian angle +Mat3 RotationAroundY(double angle); + +// Create a rotation matrix around axis Z with the provided radian angle +Mat3 RotationAroundZ(double angle); + +Mat3 rotationXYZ(double angleX, double angleY, double angleZ); + +// Degree to Radian (suppose input in [0;360]) +template +inline T degreeToRadian(T degree) +{ + static_assert(std::is_floating_point::value, "degreeToRadian: must be floating point."); + return degree * boost::math::constants::pi() / 180.0; +} + +// Radian to degree +template +inline T radianToDegree(T radian) +{ + static_assert(std::is_floating_point::value, "radianToDegree: must be floating point."); + return radian / boost::math::constants::pi() * 180.0; +} + +/// Return in radian the mean rotation amplitude of the given rotation matrix +/// Computed as the mean of matrix column dot products to an Identity matrix +double getRotationMagnitude(const Mat3& R2); + +/** + * @brief Compute the angle between two rotation matrices. + * @param[in] R1 The first rotation matrix. + * @param[in] R2 The second rotation matrix. + * @return The angle between the two rotations as the angle of the rotation + * matrix R1*R2.transpose(). + */ +double rotationDifference(const Mat3& R1, const Mat3& R2); + +inline double SIGN(double x) { return x < 0.0 ? -1.0 : 1.0; } + +// L1 norm = Sum (|x0| + |x1| + |xn|) + +template +inline double NormL1(const TVec& x) +{ + return x.array().abs().sum(); +} + +// L2 norm = Sqrt (Sum (x0^2 + x1^2 + xn^2)) + +template +inline double NormL2(const TVec& x) +{ + return x.norm(); +} + +// LInfinity norm = max (|x0|, |x1|, ..., |xn|) + +template +inline double NormLInfinity(const TVec& x) +{ + return x.array().abs().maxCoeff(); +} + +template +inline double DistanceL1(const TVec& x, const TVec& y) +{ + return (x - y).array().abs().sum(); +} + +template +inline double DistanceL2(const TVec& x, const TVec& y) +{ + return (x - y).norm(); +} + +template +inline double DistanceLInfinity(const TVec& x, const TVec& y) +{ + return NormLInfinity(x - y); +} + +template +inline bool AreVecNearEqual(const TVec& x, const TVec& y, const double epsilon) +{ + assert(x.cols() == y.cols()); + for (typename TVec::Index i = 0; i < x.cols(); ++i) + { + if ((y(i) - epsilon > x(i)) || (x(i) > y(i) + epsilon)) + return false; + } + return true; +} + +template +inline bool AreMatNearEqual(const TMat& X, const TMat& Y, const double epsilon) +{ + assert(X.cols() == Y.cols()); + assert(X.rows() == Y.rows()); + for (typename TMat::Index i = 0; i < X.rows(); ++i) + { + for (typename TMat::Index j = 0; j < X.cols(); ++j) + { + if ((Y(i, j) - epsilon > X(i, j)) || (X(i, j) > Y(i, j) + epsilon)) + return false; + } + } + return true; +} + +// Make a rotation matrix such that center becomes the direction of the +// positive z-axis, and y is oriented close to up by default. +Mat3 LookAt(const Vec3& center, const Vec3& up = Vec3::UnitY()); + +Mat3 LookAt2(const Vec3& eyePosition3D, const Vec3& center3D = Vec3::Zero(), const Vec3& upVector3D = Vec3::UnitY()); + +#define SUM_OR_DYNAMIC(x, y) (x == Eigen::Dynamic || y == Eigen::Dynamic) ? Eigen::Dynamic : (x + y) + +template +struct hstack_return +{ + using Scalar = typename Derived1::Scalar; + + enum + { + RowsAtCompileTime = Derived1::RowsAtCompileTime, + ColsAtCompileTime = SUM_OR_DYNAMIC(Derived1::ColsAtCompileTime, Derived2::ColsAtCompileTime), + Options = Derived1::Flags & Eigen::RowMajorBit ? Eigen::RowMajor : 0, + MaxRowsAtCompileTime = Derived1::MaxRowsAtCompileTime, + MaxColsAtCompileTime = SUM_OR_DYNAMIC(Derived1::MaxColsAtCompileTime, Derived2::MaxColsAtCompileTime) + }; + typedef Eigen::Matrix type; +}; + +template +typename hstack_return::type HStack(const Eigen::MatrixBase& lhs, const Eigen::MatrixBase& rhs) +{ + typename hstack_return::type res; + res.resize(lhs.rows(), lhs.cols() + rhs.cols()); + res << lhs, rhs; + return res; +} + +template +struct vstack_return +{ + using Scalar = typename Derived1::Scalar; + + enum + { + RowsAtCompileTime = SUM_OR_DYNAMIC(Derived1::RowsAtCompileTime, Derived2::RowsAtCompileTime), + ColsAtCompileTime = Derived1::ColsAtCompileTime, + Options = Derived1::Flags & Eigen::RowMajorBit ? Eigen::RowMajor : 0, + MaxRowsAtCompileTime = SUM_OR_DYNAMIC(Derived1::MaxRowsAtCompileTime, Derived2::MaxRowsAtCompileTime), + MaxColsAtCompileTime = Derived1::MaxColsAtCompileTime + }; + typedef Eigen::Matrix type; +}; + +template +typename vstack_return::type VStack(const Eigen::MatrixBase& lhs, const Eigen::MatrixBase& rhs) +{ + typename vstack_return::type res; + res.resize(lhs.rows() + rhs.rows(), lhs.cols()); + res << lhs, rhs; + return res; +} +#undef SUM_OR_DYNAMIC + +template +inline double FrobeniusNorm(const TMat& A) +{ + return A.norm(); +} + +template +inline double FrobeniusDistance(const TMat& A, const TMat& B) +{ + return FrobeniusNorm(A - B); +} + +template +double CosinusBetweenMatrices(const TMat& a, const TMat& b) +{ + return (a.array() * b.array()).sum() / FrobeniusNorm(a) / FrobeniusNorm(b); +} + +/** + * @brief Given a vector of element and a vector containing a selection of its indices, + * it returns a new vector containing only the selected elements of the original vector + * @tparam T The type of data contained in the vector. + * @param[out] result The output vector containing only the selected original elements. + * @param[in] input The input vector. + * @param[in] selection The vector containing the selection of elements of \p input + * through the indices specified in \p selection. + */ +template +void pick(std::vector& result, const std::vector& input, const std::vector::size_type>& selection) +{ + result.reserve(selection.size()); + std::transform( + selection.begin(), selection.end(), std::back_inserter(result), [&input](typename std::vector::size_type idx) { return input.at(idx); }); +} + +void MeanAndVarianceAlongRows(const Mat& A, Vec* mean_pointer, Vec* variance_pointer); + +bool exportMatToTextFile(const Mat& mat, const std::string& filename, const std::string& sPrefix = "A"); + +inline int is_finite(const double val) +{ +#ifdef _MSC_VER + return _finite(val); +#else + return std::isfinite(val); +#endif +} + +/** + ** Split a range [ a ; b [ into a set of n ranges : + [ a ; c1 [ U [ c1 ; c2 [ U ... U [ c(n-1) ; b [ + ** + Output range vector only store [ a , c1 , c2 , ... , b ] + + ** if input range can't be split (range [a;b[ size is less than nb_split, only return [a;b[ range + ** + ** @param range_start Start of range to split + ** @param range_end End of range to split + ** @param nb_split Number of desired split + ** @param d_range Output splitted range + **/ +template +void SplitRange(const T range_start, const T range_end, const int nb_split, std::vector& d_range) +{ + const T range_length = range_end - range_start; + if (range_length < nb_split) + { + d_range.push_back(range_start); + d_range.push_back(range_end); + } + else + { + const T delta_range = range_length / nb_split; + + d_range.push_back(range_start); + for (int i = 1; i < nb_split; ++i) + { + d_range.push_back(range_start + i * delta_range); + } + d_range.push_back(range_end); + } +} + +template +constexpr T divideRoundUp(T x, T y) +{ + static_assert(std::is_integral::value, "divideRoundUp only works with integer arguments"); + assert(y != 0); // Prevents division by zero + const auto xPos = x >= 0; + const auto yPos = y >= 0; + if (xPos == yPos) + { + return x / y + T((x % y) != 0); + } + else + { + // negative result, rounds towards zero anyways + return x / y; + } +} + +/** + * This function initializes the global state of random number generators that e.g. our tests + * depend on. This makes it possible to have exactly reproducible program runtime behavior + * without random changes. In order to introduce variation in execution, + * ALICEVISION_RANDOM_SEED environment variable can be set to an integer value. + * + * This function is especially useful in tests where it allows to prevent random failures. + */ +void makeRandomOperationsReproducible(); + +} // namespace aliceVision diff --git a/src/aliceVision/numeric/numeric.hpp b/src/aliceVision/numeric/numeric.hpp index 2075980c28..711c223d4e 100644 --- a/src/aliceVision/numeric/numeric.hpp +++ b/src/aliceVision/numeric/numeric.hpp @@ -31,14 +31,9 @@ // http://eigen.tuxfamily.org/dox-devel/QuickRefPage.html //-- #include -#include #include -#include -#include #include -#include -#include -#include + #include #include @@ -122,332 +117,7 @@ using Mat9 = Eigen::Matrix; using sMat = Eigen::SparseMatrix; using sRMat = Eigen::SparseMatrix; -//-------------- -//-- Function -- -//-------------- - -/// Return the square of a number. - -template -inline T Square(T x) -{ - return x * x; -} - -/// Clamp return the number if inside range, else min or max range. - -template -inline T clamp(const T& val, const T& min, const T& max) -{ - return std::max(min, std::min(val, max)); - //(val < min) ? val : ((val>max) ? val : max); -} - -inline bool isSimilar(double a, double b) -{ - const double diff = a - b; - return std::abs(diff) < 1e-8; -} -inline bool isSimilar(float a, float b) -{ - const float diff = a - b; - return std::abs(diff) < 1e-8f; -} - -/** - * @brief Create a minimal skew matrix from a 2d vector. - * @param[in] x A 2d vector whose 3rd coordinate is supposed to be 1. - * @return The minimal ske matrix: [0, -1, x(1); 1, 0, -x(0);] - */ -Mat23 SkewMatMinimal(const Vec2& x); - -/** - * @brief Create a cross product matrix from a 3d vector. - * @param x A 3d vector. - * @return the cross matrix representation of the input vector. - */ -Mat3 CrossProductMatrix(const Vec3& x); - -// Create a rotation matrix around axis X with the provided radian angle -Mat3 RotationAroundX(double angle); - -// Create a rotation matrix around axis Y with the provided radian angle -Mat3 RotationAroundY(double angle); - -// Create a rotation matrix around axis Z with the provided radian angle -Mat3 RotationAroundZ(double angle); - -Mat3 rotationXYZ(double angleX, double angleY, double angleZ); - -// Degree to Radian (suppose input in [0;360]) -template -inline T degreeToRadian(T degree) -{ - static_assert(std::is_floating_point::value, "degreeToRadian: must be floating point."); - return degree * boost::math::constants::pi() / 180.0; -} - -// Radian to degree -template -inline T radianToDegree(T radian) -{ - static_assert(std::is_floating_point::value, "radianToDegree: must be floating point."); - return radian / boost::math::constants::pi() * 180.0; -} - -/// Return in radian the mean rotation amplitude of the given rotation matrix -/// Computed as the mean of matrix column dot products to an Identity matrix -double getRotationMagnitude(const Mat3& R2); - -/** - * @brief Compute the angle between two rotation matrices. - * @param[in] R1 The first rotation matrix. - * @param[in] R2 The second rotation matrix. - * @return The angle between the two rotations as the angle of the rotation - * matrix R1*R2.transpose(). - */ -double rotationDifference(const Mat3& R1, const Mat3& R2); - -inline double SIGN(double x) { return x < 0.0 ? -1.0 : 1.0; } - -// L1 norm = Sum (|x0| + |x1| + |xn|) - -template -inline double NormL1(const TVec& x) -{ - return x.array().abs().sum(); -} - -// L2 norm = Sqrt (Sum (x0^2 + x1^2 + xn^2)) - -template -inline double NormL2(const TVec& x) -{ - return x.norm(); -} - -// LInfinity norm = max (|x0|, |x1|, ..., |xn|) - -template -inline double NormLInfinity(const TVec& x) -{ - return x.array().abs().maxCoeff(); -} - -template -inline double DistanceL1(const TVec& x, const TVec& y) -{ - return (x - y).array().abs().sum(); -} - -template -inline double DistanceL2(const TVec& x, const TVec& y) -{ - return (x - y).norm(); -} - -template -inline double DistanceLInfinity(const TVec& x, const TVec& y) -{ - return NormLInfinity(x - y); -} - -template -inline bool AreVecNearEqual(const TVec& x, const TVec& y, const double epsilon) -{ - assert(x.cols() == y.cols()); - for (typename TVec::Index i = 0; i < x.cols(); ++i) - { - if ((y(i) - epsilon > x(i)) || (x(i) > y(i) + epsilon)) - return false; - } - return true; -} - -template -inline bool AreMatNearEqual(const TMat& X, const TMat& Y, const double epsilon) -{ - assert(X.cols() == Y.cols()); - assert(X.rows() == Y.rows()); - for (typename TMat::Index i = 0; i < X.rows(); ++i) - { - for (typename TMat::Index j = 0; j < X.cols(); ++j) - { - if ((Y(i, j) - epsilon > X(i, j)) || (X(i, j) > Y(i, j) + epsilon)) - return false; - } - } - return true; -} - -// Make a rotation matrix such that center becomes the direction of the -// positive z-axis, and y is oriented close to up by default. -Mat3 LookAt(const Vec3& center, const Vec3& up = Vec3::UnitY()); - -Mat3 LookAt2(const Vec3& eyePosition3D, const Vec3& center3D = Vec3::Zero(), const Vec3& upVector3D = Vec3::UnitY()); - -#define SUM_OR_DYNAMIC(x, y) (x == Eigen::Dynamic || y == Eigen::Dynamic) ? Eigen::Dynamic : (x + y) - -template -struct hstack_return -{ - using Scalar = typename Derived1::Scalar; - - enum - { - RowsAtCompileTime = Derived1::RowsAtCompileTime, - ColsAtCompileTime = SUM_OR_DYNAMIC(Derived1::ColsAtCompileTime, Derived2::ColsAtCompileTime), - Options = Derived1::Flags & Eigen::RowMajorBit ? Eigen::RowMajor : 0, - MaxRowsAtCompileTime = Derived1::MaxRowsAtCompileTime, - MaxColsAtCompileTime = SUM_OR_DYNAMIC(Derived1::MaxColsAtCompileTime, Derived2::MaxColsAtCompileTime) - }; - typedef Eigen::Matrix type; -}; - -template -typename hstack_return::type HStack(const Eigen::MatrixBase& lhs, const Eigen::MatrixBase& rhs) -{ - typename hstack_return::type res; - res.resize(lhs.rows(), lhs.cols() + rhs.cols()); - res << lhs, rhs; - return res; -} - -template -struct vstack_return -{ - using Scalar = typename Derived1::Scalar; - - enum - { - RowsAtCompileTime = SUM_OR_DYNAMIC(Derived1::RowsAtCompileTime, Derived2::RowsAtCompileTime), - ColsAtCompileTime = Derived1::ColsAtCompileTime, - Options = Derived1::Flags & Eigen::RowMajorBit ? Eigen::RowMajor : 0, - MaxRowsAtCompileTime = SUM_OR_DYNAMIC(Derived1::MaxRowsAtCompileTime, Derived2::MaxRowsAtCompileTime), - MaxColsAtCompileTime = Derived1::MaxColsAtCompileTime - }; - typedef Eigen::Matrix type; -}; - -template -typename vstack_return::type VStack(const Eigen::MatrixBase& lhs, const Eigen::MatrixBase& rhs) -{ - typename vstack_return::type res; - res.resize(lhs.rows() + rhs.rows(), lhs.cols()); - res << lhs, rhs; - return res; -} -#undef SUM_OR_DYNAMIC - -template -inline double FrobeniusNorm(const TMat& A) -{ - return A.norm(); -} - -template -inline double FrobeniusDistance(const TMat& A, const TMat& B) -{ - return FrobeniusNorm(A - B); -} - -template -double CosinusBetweenMatrices(const TMat& a, const TMat& b) -{ - return (a.array() * b.array()).sum() / FrobeniusNorm(a) / FrobeniusNorm(b); -} - -/** - * @brief Given a vector of element and a vector containing a selection of its indices, - * it returns a new vector containing only the selected elements of the original vector - * @tparam T The type of data contained in the vector. - * @param[out] result The output vector containing only the selected original elements. - * @param[in] input The input vector. - * @param[in] selection The vector containing the selection of elements of \p input - * through the indices specified in \p selection. - */ -template -void pick(std::vector& result, const std::vector& input, const std::vector::size_type>& selection) -{ - result.reserve(selection.size()); - std::transform( - selection.begin(), selection.end(), std::back_inserter(result), [&input](typename std::vector::size_type idx) { return input.at(idx); }); -} - -void MeanAndVarianceAlongRows(const Mat& A, Vec* mean_pointer, Vec* variance_pointer); - -bool exportMatToTextFile(const Mat& mat, const std::string& filename, const std::string& sPrefix = "A"); - -inline int is_finite(const double val) -{ -#ifdef _MSC_VER - return _finite(val); -#else - return std::isfinite(val); -#endif -} - -/** - ** Split a range [ a ; b [ into a set of n ranges : - [ a ; c1 [ U [ c1 ; c2 [ U ... U [ c(n-1) ; b [ - ** - Output range vector only store [ a , c1 , c2 , ... , b ] - - ** if input range can't be split (range [a;b[ size is less than nb_split, only return [a;b[ range - ** - ** @param range_start Start of range to split - ** @param range_end End of range to split - ** @param nb_split Number of desired split - ** @param d_range Output splitted range - **/ -template -void SplitRange(const T range_start, const T range_end, const int nb_split, std::vector& d_range) -{ - const T range_length = range_end - range_start; - if (range_length < nb_split) - { - d_range.push_back(range_start); - d_range.push_back(range_end); - } - else - { - const T delta_range = range_length / nb_split; - - d_range.push_back(range_start); - for (int i = 1; i < nb_split; ++i) - { - d_range.push_back(range_start + i * delta_range); - } - d_range.push_back(range_end); - } -} - -template -constexpr T divideRoundUp(T x, T y) -{ - static_assert(std::is_integral::value, "divideRoundUp only works with integer arguments"); - assert(y != 0); // Prevents division by zero - const auto xPos = x >= 0; - const auto yPos = y >= 0; - if (xPos == yPos) - { - return x / y + T((x % y) != 0); - } - else - { - // negative result, rounds towards zero anyways - return x / y; - } -} - -/** - * This function initializes the global state of random number generators that e.g. our tests - * depend on. This makes it possible to have exactly reproducible program runtime behavior - * without random changes. In order to introduce variation in execution, - * ALICEVISION_RANDOM_SEED environment variable can be set to an integer value. - * - * This function is especially useful in tests where it allows to prevent random failures. - */ -void makeRandomOperationsReproducible(); } // namespace aliceVision + +#include \ No newline at end of file diff --git a/src/aliceVision/numeric/projection.cpp b/src/aliceVision/numeric/projection.cpp index fa854377ea..c75be7d630 100644 --- a/src/aliceVision/numeric/projection.cpp +++ b/src/aliceVision/numeric/projection.cpp @@ -8,6 +8,7 @@ #include "projection.hpp" +#include #include namespace aliceVision { diff --git a/src/aliceVision/robustEstimation/LineKernel.hpp b/src/aliceVision/robustEstimation/LineKernel.hpp index 85cbb3f894..305e1b2cab 100644 --- a/src/aliceVision/robustEstimation/LineKernel.hpp +++ b/src/aliceVision/robustEstimation/LineKernel.hpp @@ -12,6 +12,8 @@ #include #include +#include + namespace aliceVision { namespace robustEstimation { diff --git a/src/aliceVision/sfm/pipeline/panorama/ReconstructionEngine_panorama.cpp b/src/aliceVision/sfm/pipeline/panorama/ReconstructionEngine_panorama.cpp index fa17643c38..6abcfdc268 100644 --- a/src/aliceVision/sfm/pipeline/panorama/ReconstructionEngine_panorama.cpp +++ b/src/aliceVision/sfm/pipeline/panorama/ReconstructionEngine_panorama.cpp @@ -37,6 +37,7 @@ #include #include +#include #ifdef _MSC_VER #pragma warning(once : 4267) // warning C4267: 'argument' : conversion from 'size_t' to 'const int', possible loss of data diff --git a/src/aliceVision/sfm/utils/alignment.cpp b/src/aliceVision/sfm/utils/alignment.cpp index 4d47b99fdc..c83de37385 100644 --- a/src/aliceVision/sfm/utils/alignment.cpp +++ b/src/aliceVision/sfm/utils/alignment.cpp @@ -25,6 +25,8 @@ #include +#include + namespace bacc = boost::accumulators; namespace aliceVision { diff --git a/src/dependencies/MeshSDFilter/MeshNormalFilter.h b/src/dependencies/MeshSDFilter/MeshNormalFilter.h index b4f658905d..1e9a48adda 100644 --- a/src/dependencies/MeshSDFilter/MeshNormalFilter.h +++ b/src/dependencies/MeshSDFilter/MeshNormalFilter.h @@ -36,6 +36,8 @@ #include "SDFilter.h" #include #include +#include +#include namespace SDFilter { diff --git a/src/software/pipeline/main_checkerboardCalibration.cpp b/src/software/pipeline/main_checkerboardCalibration.cpp index f9a90f2ab1..8aa8fa03e6 100644 --- a/src/software/pipeline/main_checkerboardCalibration.cpp +++ b/src/software/pipeline/main_checkerboardCalibration.cpp @@ -36,6 +36,8 @@ #include +#include + #include #include #include diff --git a/src/software/pipeline/main_panoramaInit.cpp b/src/software/pipeline/main_panoramaInit.cpp index 62b82b2f32..9e9bc840c6 100644 --- a/src/software/pipeline/main_panoramaInit.cpp +++ b/src/software/pipeline/main_panoramaInit.cpp @@ -24,6 +24,7 @@ #include #include #include +#include // These constants define the current software version. // They must be updated when the command line is changed. From a751fbe5fcdadc4661e147a051564a0081cd2639 Mon Sep 17 00:00:00 2001 From: Fabien Servant Date: Sat, 20 Sep 2025 18:43:01 +0200 Subject: [PATCH 2/2] Speed up alignment --- src/aliceVision/numeric/TopN.hpp | 62 ++++++++++++++++++ src/aliceVision/sfm/utils/alignment.cpp | 83 ++++++++++++++++--------- 2 files changed, 115 insertions(+), 30 deletions(-) create mode 100644 src/aliceVision/numeric/TopN.hpp diff --git a/src/aliceVision/numeric/TopN.hpp b/src/aliceVision/numeric/TopN.hpp new file mode 100644 index 0000000000..9a99e0577d --- /dev/null +++ b/src/aliceVision/numeric/TopN.hpp @@ -0,0 +1,62 @@ +// This file is part of the AliceVision project. +// Copyright (c) 2026 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/. + +#pragma once + +#include +#include + +class TopN +{ +public: + /** + @brief Constructur + @param cacheSize the maximal number of element kept + */ + TopN(size_t cacheSize) + : _cacheSize(cacheSize) + { + + } + + void add(const double & value) + { + if (_minHeap.size() < _cacheSize) + { + _minHeap.push(value); + return; + } + + /*Only add if value is greater than minimal value*/ + if (value > _minHeap.top()) + { + _minHeap.pop(); + _minHeap.push(value); + } + } + + std::vector getAndReset() + { + std::vector result; + + //Result will be in ascending order + while (!_minHeap.empty()) + { + result.push_back(_minHeap.top()); + _minHeap.pop(); + } + + return result; + } + +private: + /** + Top value will be the smallest value + */ + std::priority_queue, std::greater> _minHeap; + + size_t _cacheSize; +}; \ No newline at end of file diff --git a/src/aliceVision/sfm/utils/alignment.cpp b/src/aliceVision/sfm/utils/alignment.cpp index c83de37385..93faded272 100644 --- a/src/aliceVision/sfm/utils/alignment.cpp +++ b/src/aliceVision/sfm/utils/alignment.cpp @@ -10,14 +10,7 @@ #include #include #include - -#include -#include -#include - -#include -#include -#include +#include #include #include @@ -27,29 +20,60 @@ #include -namespace bacc = boost::accumulators; - namespace aliceVision { namespace sfm { +std::vector split(const std::string & input, const std::string & delimiters) +{ + std::vector result; + size_t start = 0; + size_t end = input.find_first_of(delimiters); + + while (end != std::string::npos) + { + //Only if content between 2 delimiters + if (end != start) + { + result.push_back(input.substr(start, end - start)); + } + + //Next + start = end + 1; + end = input.find_first_of(delimiters, start); + } + + if (start < input.length()) + { + // Add the last part + result.push_back(input.substr(start)); + } + + return result; +} + std::istream& operator>>(std::istream& in, MarkerWithCoord& marker) { std::string token(std::istreambuf_iterator(in), {}); - std::vector markerCoord; - boost::split(markerCoord, token, boost::algorithm::is_any_of(":=")); + std::vector markerCoord = split(token, ":="); + if (markerCoord.size() != 2) + { throw std::invalid_argument("Failed to parse MarkerWithCoord from: " + token); - marker.id = boost::lexical_cast(markerCoord.front()); + } + + marker.id = std::stoi(markerCoord.front()); - std::vector coord; - boost::split(coord, markerCoord.back(), boost::algorithm::is_any_of(",;_")); + std::vector coord = split(markerCoord.back(), ",;_"); if (coord.size() != 3) + { throw std::invalid_argument("Failed to parse Marker coordinates from: " + markerCoord.back()); + } for (int i = 0; i < 3; ++i) { - marker.coord(i) = boost::lexical_cast(coord[i]); + marker.coord(i) = std::stod(coord[i]); } + return in; } @@ -939,7 +963,7 @@ IndexT getViewIdFromExpression(const sfmData::SfMData& sfmData, const std::strin try { - viewId = boost::lexical_cast(camName); + viewId = std::stoul(camName); if (!sfmData.getViews().count(viewId)) { bool found = false; @@ -968,7 +992,7 @@ IndexT getViewIdFromExpression(const sfmData::SfMData& sfmData, const std::strin } } } - catch (const boost::bad_lexical_cast&) + catch (const std::exception&) { viewId = -1; } @@ -996,17 +1020,16 @@ IndexT getViewIdFromExpression(const sfmData::SfMData& sfmData, const std::strin IndexT getCenterCameraView(const sfmData::SfMData& sfmData) { - using namespace boost::accumulators; - accumulator_set> accX, accY, accZ; + Vec3 sum = Vec3::Zero(); + size_t count = 0; for (auto& [_, pose] : sfmData.getPoses().valueRange()) { - const auto& c = pose.getTransform().center(); - accX(c(0)); - accY(c(1)); - accZ(c(2)); + sum += pose.getTransform().center(); + count++; } - const Vec3 camerasCenter(extract::mean(accX), extract::mean(accY), extract::mean(accZ)); + + const Vec3 camerasCenter = sum / double(count); double minDist = std::numeric_limits::max(); IndexT centerViewId = UndefinedIndexT; @@ -1069,15 +1092,13 @@ void computeNewCoordinateSystemFromLandmarks(const sfmData::SfMData& sfmData, const std::size_t cacheSize = 10000; const double percentile = 0.99; - using namespace boost::accumulators; - using AccumulatorMax = accumulator_set>>; - AccumulatorMax accDist(tag::tail::cache_size = cacheSize); + TopN topN(cacheSize); // Center the point cloud in [0;0;0] for (int i = 0; i < landmarksCount; ++i) { vX.col(i) -= meanPoints; - accDist(vX.col(i).norm()); + topN.add(vX.col(i).norm()); } // Perform an svd over vX*vXT (var-covar) @@ -1092,7 +1113,9 @@ void computeNewCoordinateSystemFromLandmarks(const sfmData::SfMData& sfmData, U.col(2) = -U.col(2); } - const double distMax = quantile(accDist, quantile_probability = percentile); + std::vector vTopN = topN.getAndReset(); + size_t pos = size_t(std::ceil(percentile * double(vTopN.size()))); + const double distMax = vTopN[pos]; out_S = (distMax > 0.00001 ? 1.0 / distMax : 1.0); out_R = U.transpose();