Add uncertainty computation module

INSTALL.md

@@ -42,13 +42,12 @@ Other optional libraries can enable specific features (check "CMake Options" for
 * Alembic (data I/O)
 * CCTag (feature extraction/matching and localization on CPU or GPU)
 * Cuda >= 11.0 (feature extraction and depth map computation)
-* Magma (required for UncertaintyTE)
+* Magma
 * Mosek >= 6 (linear programming)
 * OpenCV >= 3.4.11 (feature extraction, calibration module, video IO), >= 4.5 for colorchecker (mcc)
 * OpenMP (enable multi-threading)
 * PCL (Point Cloud Library) >= 1.12.1 for the registration module
 * PopSift (feature extraction on GPU)
-* UncertaintyTE (Uncertainty computation)
 * Lemon >= 1.3
 * libe57format (support reading .e57 files)
 * SWIG, Python 3 and NumPy 1.26 (Python binding for AliceVision modules)
@@ -222,11 +221,6 @@ At the end of the cmake process, a report shows for each library which version (
   Enable GPU SIFT implementation.  
   `-DPopSift_DIR:PATH=/path/to/popsift/install/lib/cmake/PopSift` (where PopSiftConfig.cmake can be found)
 
-* `ALICEVISION_USE_UNCERTAINTYTE` (default: `AUTO`)  
-  Enable Uncertainty computation.  
-  `-DUNCERTAINTYTE_DIR:PATH=/path/to/uncertaintyTE/install/` (where the `include` and `lib` folders can be found)
-  `-DMAGMA_ROOT:PATH=/path/to/magma/install/` (where the `include` and `lib` folders can be found)
-
 * `ALICEVISION_USE_OPENCV` (default: `OFF`)  
   Build with OpenCV.  
   `-DOpenCV_DIR:PATH=/path/to/opencv/install/share/OpenCV/` (where OpenCVConfig.cmake can be found)

docs/usd/schema.usda

@@ -94,7 +94,8 @@ class "AvPose" (
     bool av:locked = false
     bool av:rotationOnly = false
     bool av:removable = true
-    float6 av:uncertainty = (0, 0, 0, 0, 0, 0)
+    doc = "Symmetric 6x6 uncertainty matrix (upper-triangular: 21 elements)"
+    double[] av:uncertainty = []
 }
 
 class "AvIntrinsic" (

meshroom/aliceVision/ComputeUncertainty.py

@@ -0,0 +1,38 @@
+__version__ = "1.0"
+
+from meshroom.core import desc
+from meshroom.core.utils import VERBOSE_LEVEL
+
+
+class ComputeUncertainty(desc.AVCommandLineNode):
+
+    commandLine = "aliceVision_computeUncertainty {allParams}"
+    size = desc.DynamicNodeSize("input")
+
+    category = "Utils"
+    documentation = """ Compute solution covariance """
+
+    inputs = [
+        desc.File(
+            name="input",
+            label="SfMData",
+            description="Input SfMData file.",
+            value="",
+        ),
+        desc.ChoiceParam(
+            name="verboseLevel",
+            label="Verbose Level",
+            description="Verbosity level (fatal, error, warning, info, debug, trace).",
+            values=VERBOSE_LEVEL,
+            value="info",
+        ),
+    ]
+
+    outputs = [
+        desc.File(
+            name="output",
+            label="SfMData",
+            description="Path to the output SfM file.",
+            value="{nodeCacheFolder}/sfmData.abc",
+        )
+    ]

src/CMakeLists.txt

@@ -68,8 +68,9 @@ trilean_option(ALICEVISION_USE_CCTAG "Enable CCTAG markers" AUTO)
 trilean_option(ALICEVISION_USE_APRILTAG "Enable AprilTag markers" AUTO)
 trilean_option(ALICEVISION_USE_POPSIFT "Enable GPU SIFT implementation" AUTO)
 trilean_option(ALICEVISION_USE_ALEMBIC "Enable Alembic I/O" AUTO)
-trilean_option(ALICEVISION_USE_UNCERTAINTYTE "Enable Uncertainty computation" AUTO)
 trilean_option(ALICEVISION_USE_ONNX "Enable ONNX Runtime" AUTO)
+trilean_option(ALICEVISION_USE_MAGMA "Enable Magma for Uncertainty computation" AUTO)
+
 option(ALICEVISION_USE_ONNX_GPU "Use CUDA with ONNX Runtime" ON)
 if(APPLE)
     trilean_option(ALICEVISION_USE_CUDA "Enable CUDA" OFF)
@@ -219,6 +220,7 @@ set(CMAKE_MODULE_PATH
     ${CMAKE_MODULE_PATH} ${CMAKE_CURRENT_SOURCE_DIR}/cmake)
 include(OFA/OptimizeForArchitecture)
 OptimizeForArchitecture()
+
 set(ALICEVISION_HAVE_SSE 0 CACHE BOOL "SSE2 is available" FORCE)
 if (USE_SSE2 OR TARGET_ARCHITECTURE STREQUAL "native")
     if (MSVC AND NOT ${CMAKE_CL_64})
@@ -506,6 +508,7 @@ if (ALICEVISION_BUILD_SFM)
     endif()
 endif()
 
+
 # ==============================================================================
 # CoinUtils, Clp, Osi
 # ==============================================================================
@@ -636,31 +639,6 @@ if (NOT ALICEVISION_USE_USD STREQUAL "OFF")
     endif()
 endif()
 
-# ==============================================================================
-# UncertaintyTE
-# ==============================================================================
-# - optional, only external and enabled only if ALICEVISION_USE_UNCERTAINTYTE is ON
-# ==============================================================================
-set(ALICEVISION_HAVE_UNCERTAINTYTE 0)
-
-if (ALICEVISION_BUILD_SFM)
-    if (NOT ALICEVISION_USE_UNCERTAINTYTE STREQUAL "OFF")
-        find_package(UncertaintyTE)
-
-        if (UNCERTAINTYTE_FOUND)
-            set(ALICEVISION_HAVE_UNCERTAINTYTE 1)
-            message(STATUS "UncertaintyTE found.")
-        elseif (ALICEVISION_USE_UNCERTAINTYTE STREQUAL "ON")
-            message(SEND_ERROR "Failed to find UncertaintyTE.")
-        endif()
-    endif()
-
-    if (ALICEVISION_HAVE_UNCERTAINTYTE)
-        include_directories(${UNCERTAINTYTE_INCLUDE_DIR})
-        link_directories(${UNCERTAINTYTE_LIBRARY_DIR})
-    endif()
-endif()
-
 # ==============================================================================
 # ZLIB
 # ==============================================================================
@@ -896,6 +874,26 @@ if (ALICEVISION_BUILD_SFM)
     endif()
 endif()
 
+# ==============================================================================
+# Magma
+# ==============================================================================
+
+if ((NOT ALICEVISION_BUILD_SFM) OR (NOT ALICEVISION_HAVE_CUDA))
+    set(ALICEVISION_USE_MAGMA "OFF")
+endif()
+
+set(ALICEVISION_HAVE_MAGMA 0)
+if (NOT ALICEVISION_USE_MAGMA STREQUAL "OFF")
+
+    find_package(magma CONFIG)
+
+    if (TARGET magma::magma)
+        set(MAGMA_LIBRARIES magma::magma)
+        set(ALICEVISION_HAVE_MAGMA 1)
+    elseif (ALICEVISION_USE_MAGMA STREQUAL "ON")
+        message(SEND_ERROR "Failed to find Magma.")
+    endif()
+endif()
 
 # ==============================================================================
 # Include directories
@@ -915,7 +913,6 @@ set(ALICEVISION_INCLUDE_DIRS
         ${LEMON_INCLUDE_DIRS}
         ${EIGEN3_INCLUDE_DIRS}
         ${CERES_INCLUDE_DIRS}
-        ${UNCERTAINTYTE_INCLUDE_DIRS}
         ${MAGMA_INCLUDE_DIRS}
         ${FLANN_INCLUDE_DIRS}
         ${LP_INCLUDE_DIRS}
@@ -1017,7 +1014,6 @@ message("** Build LIDAR part: " ${ALICEVISION_BUILD_LIDAR})
 message("** Build AliceVision tests: " ${ALICEVISION_BUILD_TESTS})
 message("** Build AliceVision documentation: " ${ALICEVISION_HAVE_DOC})
 message("** Build AliceVision+OpenCV samples programs: " ${ALICEVISION_HAVE_OPENCV})
-message("** Build UncertaintyTE: " ${ALICEVISION_HAVE_UNCERTAINTYTE})
 message("** Build MeshSDFilter: " ${ALICEVISION_HAVE_MESHSDFILTER})
 message("** Build Alembic exporter: " ${ALICEVISION_HAVE_ALEMBIC})
 message("** Build SWIG Python binding: " ${ALICEVISION_BUILD_SWIG_BINDING})

src/aliceVision/CMakeLists.txt

@@ -33,6 +33,9 @@ if (ALICEVISION_BUILD_SFM)
     add_subdirectory(track)
     add_subdirectory(voctree)
     add_subdirectory(calibration)
+    if (ALICEVISION_HAVE_MAGMA)
+        add_subdirectory(uncertainty)
+    endif()
     if (ALICEVISION_HAVE_OPENCV)
         add_subdirectory(imageMasking)
         add_subdirectory(keyframe)

src/aliceVision/sfm/CMakeLists.txt

@@ -52,6 +52,7 @@ set(sfm_files_headers
     utils/poseNoise.hpp
     utils/preprocess.hpp
     utils/syntheticScene.hpp
+    utils/gauge.hpp
     LocalBundleAdjustmentGraph.hpp
     FrustumFilter.hpp
     filters.hpp
@@ -96,6 +97,7 @@ set(sfm_files_sources
     utils/poseNoise.cpp
     utils/preprocess.cpp
     utils/syntheticScene.cpp
+    utils/gauge.cpp
     LocalBundleAdjustmentGraph.cpp
     FrustumFilter.cpp
     generateReport.cpp

src/aliceVision/sfm/bundle/BundleAdjustmentCeres.cpp

@@ -245,9 +245,6 @@ void BundleAdjustmentCeres::addExtrinsicsToProblem(const sfmData::SfMData& sfmDa
         double* poseBlockPtr = poseBlock.data();
         problem.AddParameterBlock(poseBlockPtr, 6);
 
-        // add pose parameter to the all parameters blocks pointers list
-        _allParametersBlocks.push_back(poseBlockPtr);
-
         // keep the camera extrinsics constants
         if (cameraPose.isLocked() || isConstant || (!refineCenter && !refineRotation))
         {
@@ -392,9 +389,6 @@ void BundleAdjustmentCeres::addIntrinsicsToProblem(const sfmData::SfMData& sfmDa
         double* intrinsicBlockPtr = intrinsicBlock.data();
         problem.AddParameterBlock(intrinsicBlockPtr, intrinsicBlock.size());
 
-        // add intrinsic parameter to the all parameters blocks pointers list
-        _allParametersBlocks.push_back(intrinsicBlockPtr);
-
         // keep the camera intrinsic constant
         if (intrinsicPtr->isLocked() || !refineIntrinsics || intrinsicPtr->getState() == EEstimatorParameterState::CONSTANT)
         {
@@ -580,9 +574,6 @@ void BundleAdjustmentCeres::addLandmarksToProblem(const sfmData::SfMData& sfmDat
 
         double* fakeDistortionBlockPtr = &_fakeDistortionBlock;
 
-        // add landmark parameter to the all parameters blocks pointers list
-        _allParametersBlocks.push_back(landmarkBlockPtr);
-
         // iterate over 2D observation associated to the 3D landmark
         for (const auto& [viewId, observation] : landmark.getObservations())
         {
@@ -1083,7 +1074,6 @@ void BundleAdjustmentCeres::resetProblem()
 {
     _statistics = Statistics();
 
-    _allParametersBlocks.clear();
     _posesBlocks.clear();
     _intrinsicsBlocks.clear();
     _landmarksBlocks.clear();
@@ -1221,21 +1211,91 @@ void BundleAdjustmentCeres::updateFromSolution(sfmData::SfMData& sfmData, ERefin
     }
 }
 
-void BundleAdjustmentCeres::createJacobian(const sfmData::SfMData& sfmData, ERefineOptions refineOptions, ceres::CRSMatrix& jacobian)
+void BundleAdjustmentCeres::createJacobian(const sfmData::SfMData& sfmData, 
+                                        ceres::CRSMatrix& jacobian, 
+                                        std::map<IndexT, size_t> & poseToPosition, 
+                                        std::map<IndexT, size_t> & intrinsicsToPosition,
+                                        std::map<IndexT, size_t> & distortionToPosition,
+                                        std::map<IndexT, size_t> & landmarkToPosition)
 {
     std::vector<ceres::ResidualBlockId> landmarksBlockIds;
     std::vector<ceres::ResidualBlockId> temporalConstraintBlockIds;
 
-    // create problem
+    poseToPosition.clear();
+    intrinsicsToPosition.clear();
+    distortionToPosition.clear();
+    landmarkToPosition.clear();
+
+    BundleAdjustment::ERefineOptions refineOptions = BundleAdjustment::REFINE_ALL;
+    //BundleAdjustment::ERefineOptions refineOptions = BundleAdjustment::REFINE_ROTATION | BundleAdjustment::REFINE_CENTER | BundleAdjustment::REFINE_STRUCTURE;
+
+    // create problem so that we get the same constraints that the classical bundle
     ceres::Problem::Options problemOptions;
     problemOptions.loss_function_ownership = ceres::DO_NOT_TAKE_OWNERSHIP;
     ceres::Problem problem(problemOptions);
     createProblem(sfmData, refineOptions, problem, landmarksBlockIds, temporalConstraintBlockIds);
 
+
+    // Control manually the order
+    // Make sure the parameters order is stored such that the jacobian 
+    // can be processed in a controlled way.
+
+    std::vector<double*> allParametersBlocks;
+
+    // First, poses
+    for (auto & [id, ptr] : _posesBlocks)
+    {
+        if (problem.IsParameterBlockConstant(ptr.data()))
+        {
+            continue;
+        }
+
+        poseToPosition[id] = allParametersBlocks.size();
+        allParametersBlocks.push_back(ptr.data());
+    }
+
+    // Second, intrinsics
+    for (auto & [id, ptr] : _intrinsicsBlocks)
+    {
+        if (problem.IsParameterBlockConstant(ptr.data()))
+        {
+            continue;
+        }
+
+        intrinsicsToPosition[id] = allParametersBlocks.size();
+        allParametersBlocks.push_back(ptr.data());
+    }
+
+    // Third, distortion
+    for (auto & [id, ptr] : _distortionsBlocks)
+    {
+        if (problem.IsParameterBlockConstant(ptr.data()))
+        {
+            continue;
+        }
+
+        distortionToPosition[id] = allParametersBlocks.size();
+        allParametersBlocks.push_back(ptr.data());
+    }
+
+
+    // Last, landmarks
+    for (auto & [id, ptr] : _landmarksBlocks)
+    {
+        if (problem.IsParameterBlockConstant(ptr.data()))
+        {
+            continue;
+        }
+
+        landmarkToPosition[id] = allParametersBlocks.size();
+        allParametersBlocks.push_back(ptr.data());
+    }
+
+
     // configure Jacobian engine
     double cost = 0.0;
     ceres::Problem::EvaluateOptions evalOpt;
-    evalOpt.parameter_blocks = _allParametersBlocks;
+    evalOpt.parameter_blocks = allParametersBlocks;
     evalOpt.num_threads = 8;
     evalOpt.apply_loss_function = true;
 

src/aliceVision/sfm/bundle/BundleAdjustmentCeres.hpp

@@ -123,10 +123,18 @@ class BundleAdjustmentCeres : public BundleAdjustment, ceres::EvaluationCallback
     /**
      * @brief Create a jacobian CRSMatrix
      * @param[in] sfmData The input SfMData contains all the information about the reconstruction
-     * @param[in] refineOptions The chosen refine flag
      * @param[out] jacobian The jacobian CSRMatrix
+     * @param[out] poseToPosition a map which gives the parameter block index of a pose (wrt to its pose Id) in the jacobian
+     * @param[out] intrinsicsToPosition a map which gives the parameter block index of a intrinsic (wrt to its intrinsic Id) in the jacobian
+     * @param[out] distortionToPosition a map which gives the parameter block index of a distortion (wrt to its intrinsic Id) in the jacobian
+     * @param[out] landmarkToPosition a map which gives the parameter block index of a landmark (wrt to its landmark Id) in the jacobian
      */
-    void createJacobian(const sfmData::SfMData& sfmData, ERefineOptions refineOptions, ceres::CRSMatrix& jacobian);
+    void createJacobian(const sfmData::SfMData& sfmData, 
+                    ceres::CRSMatrix& jacobian, 
+                    std::map<IndexT, size_t> & poseToPosition, 
+                    std::map<IndexT, size_t> & intrinsicsToPosition,
+                    std::map<IndexT, size_t> & distortionToPosition,
+                    std::map<IndexT, size_t> & landmarkToPosition);
 
     /**
      * @brief Perform a Bundle Adjustment on the SfM scene with refinement of the requested parameters
@@ -258,13 +266,9 @@ class BundleAdjustmentCeres : public BundleAdjustment, ceres::EvaluationCallback
 
     // data wrappers for refinement
 
-    /// all parameters blocks pointers
-    std::vector<double*> _allParametersBlocks;
     /// poses blocks wrapper
     /// block: ceres angleAxis(3) + center(3)
-    std::map<IndexT, std::array<double, 6>> _posesBlocks;  // TODO : maybe we can use boost::flat_map instead of std::map ?
-    /// intrinsics blocks wrapper
-    /// block: intrinsics params
+    std::map<IndexT, std::array<double, 6>> _posesBlocks; 
     std::map<IndexT, std::vector<double>> _intrinsicsBlocks;
     std::map<IndexT, std::vector<double>> _distortionsBlocks;
     std::map<IndexT, std::shared_ptr<camera::IntrinsicBase>> _intrinsicObjects;

src/aliceVision/sfm/sfmStatistics.hpp

@@ -119,5 +119,6 @@ void computeObservationsLengthsPerView(const sfmData::SfMData& sfmData,
                                        std::vector<double>& nbResidualsPerViewFirstQuartile,
                                        std::vector<double>& nbResidualsPerViewThirdQuartile);
 
+
 }  // namespace sfm
 }  // namespace aliceVision