Add the rotation averager to GLOMAP

README.md

@@ -20,6 +20,7 @@ If you use this project for your research, please cite
     year={2024},
 }
 ```
+To use the seperate rotation averaging module, refer to [this README](docs/rotation_averager.md).
 
 ## Getting Started
 

docs/rotation_averager.md

@@ -0,0 +1,62 @@
+# Gravity-aligned Rotation Averaging with Circular Regression 
+
+[Project page](https://lpanaf.github.io/eccv24_ra1dof/) | [Paper](https://www.ecva.net/papers/eccv_2024/papers_ECCV/papers/05651.pdf) | [Supp.](https://lpanaf.github.io/assets/pdf/eccv24_ra1dof_sm.pdf)
+---
+
+## About
+
+This project aims at solving the rotation averaging problem with gravity prior. 
+To achieve this, circular regression is leveraged.
+
+If you use this project for your research, please cite
+```
+@inproceedings{pan2024ra1dof,
+    author={Pan, Linfei and Pollefeys, Marc and Barath, Daniel},
+    title={{Gravity-aligned Rotation Averaging with Circular Regression}},
+    booktitle={European Conference on Computer Vision (ECCV)},
+    year={2024},
+}
+```
+
+## Getting Started
+Install GLOMAP as instrcucted in [README](../README.md).
+Then, call the rotation averager (3 degree-of-freedom) via 
+```
+glomap rotation_averager --relpose_path RELPOSE_PATH --output_path OUTPUT_PATH
+```
+
+If gravity directions are available, call the rotation averager (1 degree-of-freedom) via
+```
+glomap rotation_averager \
+    --relpose_path RELPOSE_PATH \
+    --output_path OUTPUT_PATH \
+    --gravity_path GRAVTIY PATH 
+```
+It is recommended to set `--use_stratified=1` if only a subset of images have gravity direction. 
+If gravity measurements are subject to i.i.d. noise, they can be refined by setting `--refine_gravity=1`.
+
+
+## File Formats
+### Relative Pose
+The relative pose file is expected to be of the following format
+```
+IMAGE_NAME_1 IMAGE_NAME_2 QW QX QY QZ TX TY TZ
+```
+Only images contained in at least one relative pose will be included in the following procedure.
+The relative pose should be <sub>2</sub>R<sub>1</sub> x<sub>1</sub> + <sub>2</sub>t<sub>1</sub> = x<sub>2</sub>.
+
+### Gravity Direction
+The gravity direction file is expected to be of the following format
+```
+IMAGE_NAME GX GY GZ
+```
+The gravity direction $g$ should $[0, 1, 0]$ if the image is parallel to the ground plane, and the estimated rotation would have the property that $R_i \cdot [0, 1, 0]^\top = g$.
+If is acceptable if only a subset of all images have gravity direciton.
+If the specified image name does not match any known image name from relative pose, it is ignored.
+
+### Output
+The estimated global rotation will be in the following format
+```
+IMAGE_NAME QW QX QY QZ
+```
+Any images that are not within the largest connected component of the view-graph formed by the relative pose will be ignored.

glomap/CMakeLists.txt

@@ -1,6 +1,7 @@
 set(SOURCES
     controllers/global_mapper.cc
     controllers/option_manager.cc
+    controllers/rotation_averager.cc
     controllers/track_establishment.cc
     controllers/track_retriangulation.cc
     estimators/bundle_adjustment.cc
@@ -11,7 +12,7 @@ set(SOURCES
     estimators/view_graph_calibration.cc
     io/colmap_converter.cc
     io/colmap_io.cc
-    io/gravity_io.cc
+    io/pose_io.cc
     math/gravity.cc
     math/rigid3d.cc
     math/tree.cc
@@ -28,6 +29,7 @@ set(SOURCES
 set(HEADERS
     controllers/global_mapper.h
     controllers/option_manager.h
+    controllers/rotation_averager.h
     controllers/track_establishment.h
     controllers/track_retriangulation.h
     estimators/bundle_adjustment.h
@@ -40,7 +42,7 @@ set(HEADERS
     estimators/view_graph_calibration.h
     io/colmap_converter.h
     io/colmap_io.h
-    io/gravity_io.h
+    io/pose_io.h
     math/gravity.h
     math/l1_solver.h
     math/rigid3d.h
@@ -99,7 +101,10 @@ endif()
 add_executable(glomap_main
     glomap.cc
     exe/global_mapper.h
-    exe/global_mapper.cc)
+    exe/global_mapper.cc
+    exe/rotation_averager.h
+    exe/rotation_averager.cc
+)
 target_link_libraries(glomap_main glomap)
 
 set_target_properties(glomap_main PROPERTIES OUTPUT_NAME glomap)

glomap/controllers/option_manager.cc

@@ -1,6 +1,7 @@
 #include "option_manager.h"
 
 #include "glomap/controllers/global_mapper.h"
+#include "glomap/estimators/gravity_refinement.h"
 
 #include <boost/filesystem/operations.hpp>
 #include <boost/property_tree/ini_parser.hpp>
@@ -14,6 +15,7 @@ OptionManager::OptionManager(bool add_project_options) {
   image_path = std::make_shared<std::string>();
 
   mapper = std::make_shared<GlobalMapperOptions>();
+  gravity_refiner = std::make_shared<GravityRefinerOptions>();
   Reset();
 
   desc_->add_options()("help,h", "");
@@ -248,6 +250,18 @@ void OptionManager::AddInlierThresholdOptions() {
                               &mapper->inlier_thresholds.max_rotation_error);
 }
 
+void OptionManager::AddGravityRefinerOptions() {
+  if (added_gravity_refiner_options_) {
+    return;
+  }
+  added_gravity_refiner_options_ = true;
+  AddAndRegisterDefaultOption("GravityRefiner.max_outlier_ratio",
+                              &gravity_refiner->max_outlier_ratio);
+  AddAndRegisterDefaultOption("GravityRefiner.max_gravity_error",
+                              &gravity_refiner->max_gravity_error);
+  AddAndRegisterDefaultOption("GravityRefiner.min_num_neighbors",
+                              &gravity_refiner->min_num_neighbors);
+}
 void OptionManager::Reset() {
   const bool kResetPaths = true;
   ResetOptions(kResetPaths);
@@ -276,6 +290,7 @@ void OptionManager::ResetOptions(const bool reset_paths) {
     *image_path = "";
   }
   *mapper = GlobalMapperOptions();
+  *gravity_refiner = GravityRefinerOptions();
 }
 
 void OptionManager::Parse(const int argc, char** argv) {

glomap/controllers/option_manager.h

@@ -18,6 +18,7 @@ struct GlobalPositionerOptions;
 struct BundleAdjusterOptions;
 struct TriangulatorOptions;
 struct InlierThresholdOptions;
+struct GravityRefinerOptions;
 
 class OptionManager {
  public:
@@ -37,6 +38,7 @@ class OptionManager {
   void AddBundleAdjusterOptions();
   void AddTriangulatorOptions();
   void AddInlierThresholdOptions();
+  void AddGravityRefinerOptions();
 
   template <typename T>
   void AddRequiredOption(const std::string& name,
@@ -56,6 +58,7 @@ class OptionManager {
   std::shared_ptr<std::string> image_path;
 
   std::shared_ptr<GlobalMapperOptions> mapper;
+  std::shared_ptr<GravityRefinerOptions> gravity_refiner;
 
  private:
   template <typename T>
@@ -88,6 +91,7 @@ class OptionManager {
   bool added_bundle_adjustment_options_ = false;
   bool added_triangulation_options_ = false;
   bool added_inliers_options_ = false;
+  bool added_gravity_refiner_options_ = false;
 };
 
 template <typename T>

glomap/controllers/rotation_averager.cc

@@ -0,0 +1,64 @@
+#include "rotation_averager.h"
+
+namespace glomap {
+
+bool RotationAverager::Solve(ViewGraph& view_graph,
+                             std::unordered_map<image_t, Image>& images) {
+  view_graph.KeepLargestConnectedComponents(images);
+
+  bool solve_1dof_system = false;
+  solve_1dof_system = solve_1dof_system || options_.use_gravity;
+  solve_1dof_system = solve_1dof_system && options_.use_stratified;
+
+  ViewGraph view_graph_grav;
+  image_pair_t total_pairs = 0;
+  image_pair_t grav_pairs = 0;
+  if (solve_1dof_system) {
+    // Prepare two sets: ones all with gravity, and one does not have gravity.
+    // Solve them separately first, then solve them in a single system
+    for (const auto& [pair_id, image_pair] : view_graph.image_pairs) {
+      if (!image_pair.is_valid) continue;
+
+      image_t image_id1 = image_pair.image_id1;
+      image_t image_id2 = image_pair.image_id2;
+
+      Image& image1 = images[image_id1];
+      Image& image2 = images[image_id2];
+
+      if (!image1.is_registered || !image2.is_registered) continue;
+
+      total_pairs++;
+
+      if (image1.gravity_info.has_gravity && image2.gravity_info.has_gravity) {
+        view_graph_grav.image_pairs.emplace(
+            pair_id,
+            ImagePair(image_id1, image_id2, image_pair.cam2_from_cam1));
+        grav_pairs++;
+      }
+    }
+  }
+
+  // If there is no image pairs with gravity or most image pairs are with
+  // gravity, then just run the 3dof version
+  bool status = false;
+  status = status || grav_pairs == 0;
+  status = status || grav_pairs > total_pairs * 0.95;
+  solve_1dof_system = solve_1dof_system && (!status);
+
+  if (solve_1dof_system) {
+    // Run the 1dof optimization
+    LOG(INFO) << "Solving subset 1DoF rotation averaging problem in the mixed "
+                 "prior system";
+    int num_img_grv = view_graph_grav.KeepLargestConnectedComponents(images);
+    RotationEstimator rotation_estimator_grav(options_);
+    if (!rotation_estimator_grav.EstimateRotations(view_graph_grav, images)) {
+      return false;
+    }
+    view_graph.KeepLargestConnectedComponents(images);
+  }
+
+  RotationEstimator rotation_estimator(options_);
+  return rotation_estimator.EstimateRotations(view_graph, images);
+}
+
+}  // namespace glomap

glomap/controllers/rotation_averager.h

@@ -0,0 +1,22 @@
+#pragma once
+
+#include "glomap/estimators/global_rotation_averaging.h"
+
+namespace glomap {
+
+struct RotationAveragerOptions : public RotationEstimatorOptions {
+  bool use_stratified = true;
+};
+
+class RotationAverager {
+ public:
+  RotationAverager(const RotationAveragerOptions& options)
+      : options_(options){};
+
+  bool Solve(ViewGraph& view_graph, std::unordered_map<image_t, Image>& images);
+
+ private:
+  const RotationAveragerOptions options_;
+};
+
+}  // namespace glomap

glomap/estimators/global_rotation_averaging.cc

@@ -56,6 +56,9 @@ bool RotationEstimator::EstimateRotations(
       image.cam_from_world.rotation = Eigen::Quaterniond(AngleAxisToRotation(
           rotation_estimated_.segment(image_id_to_idx_[image_id], 3)));
     }
+    // Restore the prior position (t = -Rc = R * R_ori * t_ori = R * t_ori)
+    image.cam_from_world.translation =
+        (image.cam_from_world.rotation * image.cam_from_world.translation);
   }
 
   return true;

glomap/estimators/gravity_refinement.cc

@@ -12,6 +12,10 @@ void GravityRefiner::RefineGravity(const ViewGraph& view_graph,
       view_graph.image_pairs;
   const std::unordered_map<image_t, std::unordered_set<image_t>>&
       adjacency_list = view_graph.GetAdjacencyList();
+  if (adjacency_list.empty()) {
+    LOG(INFO) << "Adjacency list not established" << std::endl;
+    return;
+  }
 
   // Identify the images that are error prone
   int counter_rect = 0;
@@ -74,6 +78,7 @@ void GravityRefiner::RefineGravity(const ViewGraph& view_graph,
     if (gravities.size() < options_.min_num_neighbors) continue;
 
     // Then, run refinment
+    gravity = AverageGravity(gravities);
     colmap::SetSphereManifold<3>(&problem, gravity.data());
     ceres::Solver::Summary summary_solver;
     ceres::Solve(options_.solver_options, &problem, &summary_solver);

glomap/exe/rotation_averager.cc

@@ -0,0 +1,87 @@
+
+#include "glomap/controllers/rotation_averager.h"
+
+#include "glomap/controllers/option_manager.h"
+#include "glomap/estimators/gravity_refinement.h"
+#include "glomap/io/colmap_io.h"
+#include "glomap/io/pose_io.h"
+#include "glomap/types.h"
+
+#include <colmap/util/misc.h>
+#include <colmap/util/timer.h>
+
+namespace glomap {
+// -------------------------------------
+// Running Global Rotation Averager
+// -------------------------------------
+int RunRotationAverager(int argc, char** argv) {
+  std::string relpose_path;
+  std::string output_path;
+  std::string gravity_path = "";
+
+  bool use_stratified = true;
+  bool refine_gravity = false;
+
+  OptionManager options;
+  options.AddRequiredOption("relpose_path", &relpose_path);
+  options.AddRequiredOption("output_path", &output_path);
+  options.AddDefaultOption("gravity_path", &gravity_path);
+  options.AddDefaultOption("use_stratified", &use_stratified);
+  options.AddDefaultOption("refine_gravity", &refine_gravity);
+  options.AddGravityRefinerOptions();
+  options.Parse(argc, argv);
+
+  if (!colmap::ExistsFile(relpose_path)) {
+    LOG(ERROR) << "`relpose_path` is not a file";
+    return EXIT_FAILURE;
+  }
+
+  if (gravity_path != "" && !colmap::ExistsFile(gravity_path)) {
+    LOG(ERROR) << "`gravity_path` is not a file";
+    return EXIT_FAILURE;
+  }
+
+  RotationAveragerOptions rotation_averager_options;
+  rotation_averager_options.skip_initialization = true;
+  rotation_averager_options.use_gravity = true;
+
+  rotation_averager_options.use_stratified = use_stratified;
+
+  // Load the database
+  ViewGraph view_graph;
+  std::unordered_map<image_t, Image> images;
+
+  ReadRelPose(relpose_path, images, view_graph);
+
+  if (gravity_path != "") {
+    ReadGravity(gravity_path, images);
+  }
+
+  int num_img = view_graph.KeepLargestConnectedComponents(images);
+  LOG(INFO) << num_img << " / " << images.size()
+            << " are within the largest connected component";
+
+  if (refine_gravity && gravity_path != "") {
+    GravityRefiner grav_refiner(*options.gravity_refiner);
+    grav_refiner.RefineGravity(view_graph, images);
+  }
+
+  RotationAverager rotation_averager(rotation_averager_options);
+  colmap::Timer run_timer;
+  run_timer.Start();
+  if (!rotation_averager.Solve(view_graph, images)) {
+    LOG(ERROR) << "Failed to solve global rotation averaging";
+    return EXIT_FAILURE;
+  }
+  run_timer.Pause();
+  LOG(INFO) << "Global rotation averaging done in "
+            << run_timer.ElapsedSeconds() << " seconds";
+
+  // Write out the estimated rotation
+  WriteGlobalRotation(output_path, images);
+  LOG(INFO) << "Global rotation averaging done" << std::endl;
+
+  return EXIT_SUCCESS;
+}
+
+}  // namespace glomap