feat: angular residual errors

Summary: This Diff add angular errors to OpenSfM quality report.

Reviewed By: paulinus

Differential Revision: D28475730

fbshipit-source-id: 7f22e6b7fd7c461739a4363c699ebc8ecabbccbe

Author: YanNoun

opensfm/pymap.pyi

@@ -162,6 +162,20 @@ class Map:
     def update_rig_instance(self, arg0) -> RigInstance: ...
     def update_shot(self, arg0) -> Shot: ...
 
+class ErrorType:
+    Pixel: Any = ...
+    Normalized: Any = ...
+    Angular: Any = ...
+    def __init__(self, arg0: int) -> None: ...
+    def __eq__(self, arg0: ErrorType) -> bool: ...
+    def __getstate__(self) -> tuple: ...
+    def __hash__(self) -> int: ...
+    def __int__(self) -> int: ...
+    def __ne__(self, arg0: ErrorType) -> bool: ...
+    def __setstate__(self, arg0: tuple) -> None: ...
+    @property
+    def __members__(self) -> dict: ...
+
 class PanoShotView:
     def __init__(self, arg0: Map) -> None: ...
     def get(self, arg0: str) -> Shot: ...

opensfm/report.py

@@ -319,10 +319,11 @@ def make_reconstruction_details(self):
 
         rows = [
             [
-                "Average Reprojection Error (normalized / pixels)",
+                "Average Reprojection Error (normalized / pixels / angular)",
                 (
                     f"{self.stats['reconstruction_statistics']['reprojection_error_normalized']:.2f} / "
-                    f"{self.stats['reconstruction_statistics']['reprojection_error_pixels']:.2f}"
+                    f"{self.stats['reconstruction_statistics']['reprojection_error_pixels']:.2f} / "
+                    f"{self.stats['reconstruction_statistics']['reprojection_error_angular']:.5f}"
                 ),
             ],
             [

opensfm/src/map/map.h

@@ -169,9 +169,10 @@ class Map {
   TracksManager ToTracksManager() const;
 
   // Tracks manager x Reconstruction intersection functions
+  enum ErrorType { Pixel = 0x0, Normalized = 0x1, Angular = 0x2 };
   std::unordered_map<ShotId, std::unordered_map<LandmarkId, Vec2d> >
   ComputeReprojectionErrors(const TracksManager& tracks_manager,
-                            bool scaled) const;
+                            const ErrorType& error_type) const;
   std::unordered_map<ShotId, std::unordered_map<LandmarkId, Observation> >
   GetValidObservations(const TracksManager& tracks_manager) const;
 

opensfm/src/map/python/pybind.cc

@@ -45,6 +45,12 @@ void DeclareShotMeasurement(py::module &m, const std::string &type_name) {
           }));
 }
 PYBIND11_MODULE(pymap, m) {
+  py::enum_<map::Map::ErrorType>(m, "ErrorType")
+      .value("Pixel", map::Map::Pixel)
+      .value("Normalized", map::Map::Normalized)
+      .value("Angular", map::Map::Angular)
+      .export_values();
+
   py::class_<map::Map>(m, "Map")
       .def(py::init())
       // Camera

opensfm/src/map/src/map.cc

@@ -4,6 +4,7 @@
 #include <map/rig.h>
 #include <map/shot.h>
 
+#include <cmath>
 #include <unordered_set>
 
 namespace map {
@@ -364,7 +365,7 @@ RigInstance& Map::GetRigInstance(const RigInstanceId& instance_id) {
   return it->second;
 }
 
-const RigInstance& Map::GetRigInstance(const RigInstanceId& instance_id) const{
+const RigInstance& Map::GetRigInstance(const RigInstanceId& instance_id) const {
   const auto& it = rig_instances_.find(instance_id);
   if (it == rig_instances_.end()) {
     throw std::runtime_error("Accessing invalid RigInstance index");
@@ -378,7 +379,7 @@ bool Map::HasRigInstance(const RigInstanceId& instance_id) const {
 
 std::unordered_map<ShotId, std::unordered_map<LandmarkId, Vec2d> >
 Map::ComputeReprojectionErrors(const TracksManager& tracks_manager,
-                               bool scaled) const {
+                               const Map::ErrorType& error_type) const {
   std::unordered_map<ShotId, std::unordered_map<LandmarkId, Vec2d> > errors;
   for (const auto& shot_id : tracks_manager.GetShotIds()) {
     const auto find_shot = shots_.find(shot_id);
@@ -394,11 +395,27 @@ Map::ComputeReprojectionErrors(const TracksManager& tracks_manager,
         continue;
       }
 
-      const Vec2d error_2d =
-          (track_n_obs.second.point -
-           shot.Project(find_landmark->second.GetGlobalPos()));
-      const auto scale = scaled ? track_n_obs.second.scale : 1.0;
-      per_shot[track_n_obs.first] = error_2d / scale;
+      if (error_type == Map::ErrorType::Pixel) {
+        const Vec2d error_2d =
+            (track_n_obs.second.point -
+             shot.Project(find_landmark->second.GetGlobalPos()));
+        per_shot[track_n_obs.first] = error_2d;
+      }
+      if (error_type == Map::ErrorType::Normalized) {
+        const Vec2d error_2d =
+            (track_n_obs.second.point -
+             shot.Project(find_landmark->second.GetGlobalPos()));
+        per_shot[track_n_obs.first] = error_2d / track_n_obs.second.scale;
+      }
+      if (error_type == Map::ErrorType::Angular) {
+        const Vec3d point =
+            (find_landmark->second.GetGlobalPos() - shot.GetPose()->GetOrigin())
+                .normalized();
+        const Vec3d bearing =
+            shot.Bearing(track_n_obs.second.point).normalized();
+        const double angle = std::acos(point.dot(bearing));
+        per_shot[track_n_obs.first] = Vec2d::Constant(angle);
+      }
     }
   }
   return errors;

opensfm/src/map/test/map_test.cc

@@ -245,7 +245,7 @@ TEST_F(EmptyMapFixture, ConstructSmallProblem) {
   ASSERT_EQ(map.GetLandmarks().size(), n_points);
 }
 
-TEST_F(OneCameraMapFixture, ComputeReprojectionError) {
+TEST_F(OneCameraMapFixture, ComputeReprojectionErrorNormalized) {
   const auto& shot = map.CreateShot("0", "0");
   Eigen::Vector3d pos = Eigen::Vector3d::Random();
   map.CreateLandmark("1", pos);
@@ -257,7 +257,8 @@ TEST_F(OneCameraMapFixture, ComputeReprojectionError) {
   const Observation o(0., 0., scale, 1, 1, 1, 1, 1, 1);
   manager.AddObservation("0", "1", o);
 
-  auto errors = map.ComputeReprojectionErrors(manager, true);
+  auto errors =
+      map.ComputeReprojectionErrors(manager, map::Map::ErrorType::Normalized);
   const auto computed = errors["0"]["1"];
   ASSERT_NEAR(expected[0] / scale, computed[0], 1e-8);
   ASSERT_NEAR(expected[1] / scale, computed[1], 1e-8);

opensfm/stats.py

@@ -11,7 +11,7 @@
 import matplotlib.colors as colors
 import matplotlib.pyplot as plt
 import numpy as np
-from opensfm import io, multiview, feature_loader
+from opensfm import io, multiview, feature_loader, pymap
 from opensfm.dataset import DataSet, DataSetBase
 
 RESIDUAL_PIXEL_CUTOFF = 4
@@ -98,9 +98,10 @@ def gcp_errors(data: DataSetBase, reconstructions):
 
 def _compute_errors(reconstructions, tracks_manager):
     @lru_cache(10)
-    def _compute_errors_cached(index, normalized):
+    def _compute_errors_cached(index, error_type):
         return reconstructions[index].map.compute_reprojection_errors(
-            tracks_manager, normalized
+            tracks_manager,
+            error_type,
         )
 
     return _compute_errors_cached
@@ -115,38 +116,52 @@ def _get_valid_observations_cached(index):
 
 
 def _projection_error(tracks_manager, reconstructions):
-    all_errors_normalized, all_errors_pixels = [], []
-    average_error_normalized, average_error_pixels = 0, 0
+    all_errors_normalized, all_errors_pixels, all_errors_angular = [], [], []
+    average_error_normalized, average_error_pixels, average_error_angular = 0, 0, 0
     for i in range(len(reconstructions)):
-        errors_normalized = _compute_errors(reconstructions, tracks_manager)(i, True)
-        errors_unnormalized = _compute_errors(reconstructions, tracks_manager)(i, False)
+        errors_normalized = _compute_errors(reconstructions, tracks_manager)(
+            i, pymap.ErrorType.Normalized
+        )
+        errors_unnormalized = _compute_errors(reconstructions, tracks_manager)(
+            i, pymap.ErrorType.Pixel
+        )
+        errors_angular = _compute_errors(reconstructions, tracks_manager)(
+            i, pymap.ErrorType.Angular
+        )
 
         for shot_id, shot_errors_normalized in errors_normalized.items():
             shot = reconstructions[i].get_shot(shot_id)
             normalizer = max(shot.camera.width, shot.camera.height)
 
-            for error_normalized, error_unnormalized in zip(
-                shot_errors_normalized.values(), errors_unnormalized[shot_id].values()
+            for error_normalized, error_unnormalized, error_angular in zip(
+                shot_errors_normalized.values(),
+                errors_unnormalized[shot_id].values(),
+                errors_angular[shot_id].values(),
             ):
                 norm_pixels = _norm2d(error_unnormalized * normalizer)
                 norm_normalized = _norm2d(error_normalized)
+                norm_angle = error_angular[0]
                 if norm_pixels > RESIDUAL_PIXEL_CUTOFF:
                     continue
                 average_error_normalized += norm_normalized
                 average_error_pixels += norm_pixels
+                average_error_angular += norm_angle
                 all_errors_normalized.append(norm_normalized)
                 all_errors_pixels.append(norm_pixels)
+                all_errors_angular.append(norm_angle)
 
     error_count = len(all_errors_normalized)
     if error_count == 0:
-        return (-1.0, -1.0, ([], []), ([], []))
+        return (-1.0, -1.0, -1.0, ([], []), ([], []), ([], []))
 
     bins = 30
     return (
         average_error_normalized / error_count,
         average_error_pixels / error_count,
+        average_error_angular / error_count,
         np.histogram(all_errors_normalized, bins),
         np.histogram(all_errors_pixels, bins),
+        np.histogram(all_errors_angular, bins),
     )
 
 
@@ -203,11 +218,14 @@ def reconstruction_statistics(data: DataSetBase, tracks_manager, reconstructions
     (
         avg_normalized,
         avg_pixels,
+        avg_angular,
         (hist_normalized, bins_normalized),
         (hist_pixels, bins_pixels),
+        (hist_angular, bins_angular),
     ) = _projection_error(tracks_manager, reconstructions)
     stats["reprojection_error_normalized"] = avg_normalized
     stats["reprojection_error_pixels"] = avg_pixels
+    stats["reprojection_error_angular"] = avg_angular
     stats["reprojection_histogram_normalized"] = (
         list(map(float, hist_normalized)),
         list(map(float, bins_normalized)),
@@ -216,6 +234,10 @@ def reconstruction_statistics(data: DataSetBase, tracks_manager, reconstructions
         list(map(float, hist_pixels)),
         list(map(float, bins_pixels)),
     )
+    stats["reprojection_histogram_angular"] = (
+        list(map(float, hist_angular)),
+        list(map(float, bins_angular)),
+    )
 
     return stats
 
@@ -475,7 +497,7 @@ def save_residual_histogram(
     io_handler,
 ):
     backup = dict(mpl.rcParams)
-    fig, axs = plt.subplots(1, 2, tight_layout=True, figsize=(15, 3))
+    fig, axs = plt.subplots(1, 3, tight_layout=True, figsize=(15, 3))
 
     h_norm, b_norm = stats["reconstruction_statistics"][
         "reprojection_histogram_normalized"
@@ -493,8 +515,19 @@ def save_residual_histogram(
     for i in range(len(p_pixel)):
         p_pixel[i].set_facecolor(plt.cm.viridis(n[i] / max(n)))
 
+    h_angular, b_angular = stats["reconstruction_statistics"][
+        "reprojection_histogram_angular"
+    ]
+    n, _, p_angular, = axs[
+        2
+    ].hist(b_angular[:-1], b_angular, weights=h_angular)
+    n = n.astype("int")
+    for i in range(len(p_angular)):
+        p_angular[i].set_facecolor(plt.cm.viridis(n[i] / max(n)))
+
     axs[0].set_title("Normalized Residual")
     axs[1].set_title("Pixel Residual")
+    axs[2].set_title("Angular Residual")
 
     with io_handler.open(
         os.path.join(output_path, "residual_histogram.png"), "wb"
@@ -768,8 +801,8 @@ def save_residual_grids(
 
     for i in range(len(reconstructions)):
         valid_observations = _get_valid_observations(reconstructions, tracks_manager)(i)
-        errors_scaled = _compute_errors(reconstructions, tracks_manager)(i, True)
-        errors_unscaled = _compute_errors(reconstructions, tracks_manager)(i, False)
+        errors_scaled = _compute_errors(reconstructions, tracks_manager)(i, pymap.ErrorType.Normalized)
+        errors_unscaled = _compute_errors(reconstructions, tracks_manager)(i, pymap.ErrorType.Pixel)
 
         for shot_id, shot_errors in errors_scaled.items():
             shot = reconstructions[i].get_shot(shot_id)