On-the-fly resection candidates

Instead of recomputing resection candidates from scratch, which has O(N2) for the entire processing (and takes up to 1-2 seconds on large projects), we instead keep track of candidates per images. This makes everything constant time instead.

On helenenschacht (https://github.com/OpenDroneMap/odm_data_helenenschacht/tree/main) :

BEFORE : around 120 ms
2025-08-15 13:01:55,007 INFO: Reconstruction now has 3 shots.
2025-08-15 13:01:55,138 INFO: -------------------------------------------------------

2025-08-15 13:01:56,832 INFO: Reconstruction now has 4 shots.
2025-08-15 13:01:56,965 INFO: -------------------------------------------------------

2025-08-15 13:01:57,359 INFO: Reconstruction now has 5 shots.
2025-08-15 13:01:57,492 INFO: -------------------------------------------------------

AFTER : constant < 1 ms
2025-08-15 13:10:39,368 INFO: Reconstruction now has 3 shots.
2025-08-15 13:10:39,368 INFO: -------------------------------------------------------

2025-08-15 13:10:41,195 INFO: Reconstruction now has 4 shots.
2025-08-15 13:10:41,195 INFO: -------------------------------------------------------

2025-08-15 13:10:41,706 INFO: Reconstruction now has 5 shots.
2025-08-15 13:10:41,706 INFO: -------------------------------------------------------

On 7K images datasets (https://github.com/zivillian/odm_ziegeleipark/tree/master), this takes up to 1-2 seconds after a few thousands images

Author: Yann N.

opensfm/reconstruction.py

@@ -705,7 +705,7 @@ def resect(
     bs = np.array(bs)
     Xs = np.array(Xs)
     if len(bs) < 5:
-        return False, set(), {"num_common_points": len(bs)}
+        return False, set(), set(), {"num_common_points": len(bs)}
 
     T = multiview.absolute_pose_ransac(bs, Xs, threshold, 1000, 0.999)
 
@@ -736,15 +736,17 @@ def resect(
             triangulate_shot_features(
                 tracks_manager, reconstruction, new_shots, data.config
             )
+        tracks = set()
         for i, succeed in enumerate(inliers):
             if succeed:
                 add_observation_to_reconstruction(
                     tracks_manager, reconstruction, shot_id, ids[i]
                 )
+                tracks.add(ids[i])
         report["shots"] = list(new_shots)
-        return True, new_shots, report
+        return True, new_shots, tracks, report
     else:
-        return False, set(), report
+        return False, set(), set(), report
 
 
 def corresponding_tracks(
@@ -911,7 +913,7 @@ def triangulate_robust(
         min_ray_angle_degrees: float,
         min_depth: float,
         iterations: int,
-    ) -> None:
+    ) -> bool:
         """Triangulate track in a RANSAC way and add point to reconstruction."""
         os, bs, ids = [], [], []
         for shot_id, obs in self.tracks_handler.get_observations(track).items():
@@ -923,7 +925,7 @@ def triangulate_robust(
             ids.append(shot_id)
 
         if len(ids) < 2:
-            return
+            return False
 
         os = np.array(os)
         bs = np.array(bs)
@@ -1003,6 +1005,7 @@ def triangulate_robust(
             self.tracks_handler.store_track_coordinates(track, best_point)
             for i in best_inliers:
                 self.tracks_handler.store_inliers_observation(track, ids[i])
+        return len(best_inliers) > 1
 
     def triangulate(
         self,
@@ -1011,7 +1014,7 @@ def triangulate(
         min_ray_angle_degrees: float,
         min_depth: float,
         iterations: int,
-    ) -> None:
+    ) -> bool:
         """Triangulate track and add point to reconstruction."""
         os, bs, ids = [], [], []
         for shot_id, obs in self.tracks_handler.get_observations(track).items():
@@ -1039,6 +1042,9 @@ def triangulate(
                 self.tracks_handler.store_track_coordinates(track, X.tolist())
                 for shot_id in ids:
                     self.tracks_handler.store_inliers_observation(track, shot_id)
+                return True
+
+        return False            
 
     def triangulate_dlt(
         self,
@@ -1047,7 +1053,7 @@ def triangulate_dlt(
         min_ray_angle_degrees: float,
         min_depth: float,
         iterations: int,
-    ) -> None:
+    ) -> bool:
         """Triangulate track using DLT and add point to reconstruction."""
         Rts, bs, os, ids = [], [], [], []
         for shot_id, obs in self.tracks_handler.get_observations(track).items():
@@ -1075,6 +1081,8 @@ def triangulate_dlt(
                 self.tracks_handler.store_track_coordinates(track, X.tolist())
                 for shot_id in ids:
                     self.tracks_handler.store_inliers_observation(track, shot_id)
+                return True
+        return False
 
     def _shot_origin(self, shot: pymap.Shot) -> NDArray:
         if shot.id in self.origins:
@@ -1106,7 +1114,7 @@ def triangulate_shot_features(
     reconstruction: types.Reconstruction,
     shot_ids: Set[str],
     config: Dict[str, Any],
-) -> None:
+) -> List[str]:
     """Reconstruct as many tracks seen in shot_id as possible."""
     reproj_threshold = config["triangulation_threshold"]
     min_ray_angle = config["triangulation_min_ray_angle"]
@@ -1124,24 +1132,29 @@ def triangulate_shot_features(
         if s in all_shots_ids
         for t in tracks_manager.get_shot_observations(s)
     }
+
+    created_tracks = []
     for track in tracks_ids:
         if track not in reconstruction.points:
             if config["triangulation_type"] == "ROBUST":
-                triangulator.triangulate_robust(
+                if(triangulator.triangulate_robust(
                     track,
                     reproj_threshold,
                     min_ray_angle,
                     min_depth,
                     refinement_iterations,
-                )
+                )):
+                    created_tracks.append(track)
             elif config["triangulation_type"] == "FULL":
-                triangulator.triangulate(
+                if(triangulator.triangulate(
                     track,
                     reproj_threshold,
                     min_ray_angle,
                     min_depth,
                     refinement_iterations,
-                )
+                )):
+                    created_tracks.append(track)
+    return created_tracks
 
 
 def retriangulate(
@@ -1214,7 +1227,7 @@ def remove_outliers(
     reconstruction: types.Reconstruction,
     config: Dict[str, Any],
     points: Optional[Dict[str, pymap.Landmark]] = None,
-) -> int:
+) -> Tuple[List[Tuple[str, str]], Set[str]]:
     """Remove points with large reprojection error.
 
     A list of point ids to be processed can be given in ``points``.
@@ -1236,14 +1249,16 @@ def remove_outliers(
         reconstruction.map.remove_observation(shot_id, track)
         track_ids.add(track)
 
+    removed_tracks = set()
     for track in track_ids:
         if track in reconstruction.points:
             lm = reconstruction.points[track]
             if lm.number_of_observations() < 2:
                 reconstruction.map.remove_landmark(lm)
+                removed_tracks.add(track)
 
     logger.info("Removed outliers: {}".format(len(outliers)))
-    return len(outliers)
+    return outliers, removed_tracks
 
 
 def shot_lla_and_compass(
@@ -1407,6 +1422,75 @@ def should(self) -> bool:
     def done(self) -> None:
         self.num_points_last = len(self.reconstruction.points)
 
+class ResectionCandidates:
+    """Helper to keep track of resection candidates."""
+
+    def __init__(self, tracks_manager: pymap.TracksManager, reconstruction: types.Reconstruction) -> None:
+        self.tracks_manager = tracks_manager
+        self.reconstruction = reconstruction
+        self.candidates: Dict[str, Set[str]] = {}
+        self.tracks = set()
+
+    def update(self, reconstruction: types.Reconstruction) -> None:
+        """Update candidates based on a new reconstruction."""
+        self.reconstruction = reconstruction
+        rec_tracks = set(reconstruction.points.keys())
+
+        # add new tracks
+        for track_id in rec_tracks:
+            if track_id not in self.tracks:
+                self.tracks.add(track_id)
+                for shot_id in self.tracks_manager.get_track_observations(track_id):
+                    if shot_id not in self.candidates:
+                        self.candidates[shot_id] = set()
+                    self.candidates[shot_id].add(track_id)
+
+        # remove tracks that are not in the reconstruction
+        to_remove = set()
+        for track_id in self.tracks:
+            if track_id not in rec_tracks:
+                for shot_id, shot_candidates in self.candidates.items():
+                    if track_id in shot_candidates:
+                        shot_candidates.remove(track_id)
+                to_remove.add(track_id)
+
+        for track_id in to_remove:
+            self.tracks.remove(track_id)
+                
+    def add(self, shots: List[str], tracks: Set[str]) -> None:
+        """Add newly resected tracks and the shots they belong to."""
+        for shot_id in shots:
+            if shot_id in self.candidates:
+                del self.candidates[shot_id]
+
+        for track_id in tracks:
+            self.tracks.add(track_id)
+            for shot_id in self.tracks_manager.get_track_observations(track_id):
+                if shot_id in self.reconstruction.shots:
+                    continue
+                if shot_id not in self.candidates:
+                    self.candidates[shot_id] = set()
+                self.candidates[shot_id].add(track_id)
+
+    def remove(self, outliers: List[Tuple[str, str]], removed_tracks: List[str]) -> None:
+        """Remove outliers from candidates."""
+        for shot_id, track_id in outliers:
+            if shot_id in self.candidates and track_id in self.candidates[shot_id]:
+                self.candidates[shot_id].remove(track_id)
+                if not self.candidates[shot_id]:
+                    del self.candidates[shot_id]
+
+        for track_id in removed_tracks:
+            for shot_id, shot_candidates in self.candidates.items():
+                if track_id in shot_candidates:
+                    shot_candidates.remove(track_id)
+            self.tracks.remove(track_id)
+
+
+    def get_candidates(self, images: Set[str]) -> List[Tuple[str, int]]:
+        """Get sorted candidates for resection."""
+        return [(k, len(v)) for k,v in filter(lambda x: x[0] in images, sorted(self.candidates.items(), key=lambda x: -len(x[1])))]
+
 
 def grow_reconstruction(
     data: DataSetBase,
@@ -1429,8 +1513,15 @@ def grow_reconstruction(
     remove_outliers(reconstruction, config)
     paint_reconstruction(data, tracks_manager, reconstruction)
 
+    resection_candidates = ResectionCandidates(tracks_manager, reconstruction)
     should_bundle = ShouldBundle(data, reconstruction)
     should_retriangulate = ShouldRetriangulate(data, reconstruction)
+
+    resection_candidates.add(
+        list(reconstruction.shots.keys()),
+        list(reconstruction.points.keys()),
+    )
+    
     while True:
         if config["save_partial_reconstructions"]:
             paint_reconstruction(data, tracks_manager, reconstruction)
@@ -1441,17 +1532,15 @@ def grow_reconstruction(
                 ),
             )
 
-        candidates = reconstructed_points_for_images(
-            tracks_manager, reconstruction, images
-        )
+        candidates = resection_candidates.get_candidates(images)
         if not candidates:
             break
 
         logger.info("-------------------------------------------------------")
         threshold = data.config["resection_threshold"]
         min_inliers = data.config["resection_min_inliers"]
         for image, _ in candidates:
-            ok, new_shots, resrep = resect(
+            ok, new_shots, resected_tracks, resrep = resect(
                 data,
                 tracks_manager,
                 reconstruction,
@@ -1470,6 +1559,7 @@ def grow_reconstruction(
                 rig_camera_priors,
                 data.config,
             )
+            resection_candidates.add(new_shots, resected_tracks)
 
             logger.info(f"Adding {' and '.join(new_shots)} to the reconstruction")
             step: Dict[str, Union[List[int], List[str], int, List[int], Any]] = {
@@ -1480,7 +1570,8 @@ def grow_reconstruction(
             report["steps"].append(step)
 
             np_before = len(reconstruction.points)
-            triangulate_shot_features(tracks_manager, reconstruction, new_shots, config)
+            new_tracks = triangulate_shot_features(tracks_manager, reconstruction, new_shots, config)
+            resection_candidates.add([], new_tracks)
             np_after = len(reconstruction.points)
             step["triangulated_points"] = np_after - np_before
 
@@ -1491,10 +1582,11 @@ def grow_reconstruction(
                     reconstruction, camera_priors, rig_camera_priors, None, config
                 )
                 rrep = retriangulate(tracks_manager, reconstruction, config)
+                resection_candidates.update(reconstruction)
                 b2rep = bundle(
                     reconstruction, camera_priors, rig_camera_priors, None, config
                 )
-                remove_outliers(reconstruction, config)
+                resection_candidates.remove(*remove_outliers(reconstruction, config))
                 step["bundle"] = b1rep
                 step["retriangulation"] = rrep
                 step["bundle_after_retriangulation"] = b2rep
@@ -1505,7 +1597,7 @@ def grow_reconstruction(
                 brep = bundle(
                     reconstruction, camera_priors, rig_camera_priors, None, config
                 )
-                remove_outliers(reconstruction, config)
+                resection_candidates.remove(*remove_outliers(reconstruction, config))
                 step["bundle"] = brep
                 should_bundle.done()
             elif config["local_bundle_radius"] > 0:
@@ -1517,7 +1609,7 @@ def grow_reconstruction(
                     image,
                     config,
                 )
-                remove_outliers(reconstruction, config, bundled_points)
+                resection_candidates.remove(*remove_outliers(reconstruction, config, bundled_points))
                 step["local_bundle"] = brep
 
             logger.info(f"Reconstruction now has {len(reconstruction.shots)} shots.")
@@ -1536,7 +1628,7 @@ def grow_reconstruction(
         config = overidden_config
 
     bundle(reconstruction, camera_priors, rig_camera_priors, gcp, config)
-    remove_outliers(reconstruction, config)
+    resection_candidates.remove(*remove_outliers(reconstruction, config))
     paint_reconstruction(data, tracks_manager, reconstruction)
     return reconstruction, report