Fix GCP handling

conda.yml

@@ -6,7 +6,9 @@ dependencies:
     - libxcrypt
     - libopencv *=headless*
     - py-opencv
-    - conda-forge::libvorbis
-    - ceres-solver=2.1
+    - anaconda::ceres-solver=2.2
     - conda-forge::llvm-openmp
     - conda-forge::cxx-compiler
+    - anaconda::suitesparse
+    - conda-forge::openblas
+    - anaconda::metis

opensfm/actions/create_tracks.py

@@ -13,9 +13,10 @@ def run_dataset(data: DataSetBase) -> None:
         data, data.images()
     )
     features_end = timer()
-    matches = tracking.load_matches(data, data.images())
+    matches = lambda: tracking.load_matches(data, data.images())
     matches_end = timer()
-    tracks_manager = tracking.create_tracks_manager(
+
+    tracks_manager = tracking.create_tracks_manager_from_matches_iter(
         features,
         colors,
         segmentations,

opensfm/actions/reconstruct.py

@@ -1,4 +1,5 @@
 # pyre-strict
+import gc
 from opensfm import io, reconstruction
 from opensfm.dataset_base import DataSetBase
 
@@ -21,5 +22,8 @@ def run_dataset(
     else:
         raise RuntimeError(f"Unsupported algorithm for reconstruction {algorithm}")
 
-    data.save_reconstruction(reconstructions)
+    del tracks_manager
+    gc.collect()
+
     data.save_report(io.json_dumps(report), "reconstruction.json")
+    data.save_reconstruction(reconstructions)

opensfm/align.py

@@ -127,11 +127,13 @@ def alignment_constraints(
     X, Xp = [], []
     # Get Ground Control Point correspondences
     if gcp and config["bundle_use_gcp"]:
-        triangulated, measured = triangulate_all_gcp(reconstruction, gcp)
+        triangulated, measured = triangulate_all_gcp(
+            reconstruction, gcp, config["gcp_reprojection_error_threshold"]
+        )
         X.extend(triangulated)
         Xp.extend(measured)
     # Get camera center correspondences
-    if use_gps and config["bundle_use_gps"]:
+    elif use_gps and config["bundle_use_gps"]:
         for rig_instance in reconstruction.rig_instances.values():
             gpses = [
                 shot.metadata.gps_position.value
@@ -433,14 +435,17 @@ def get_horizontal_and_vertical_directions(
 
 
 def triangulate_all_gcp(
-    reconstruction: types.Reconstruction, gcp: List[pymap.GroundControlPoint]
+    reconstruction: types.Reconstruction,
+    gcp: List[pymap.GroundControlPoint],
+    threshold: float,
 ) -> Tuple[List[NDArray], List[NDArray]]:
     """Group and triangulate Ground Control Points seen in 2+ images."""
     triangulated, measured = [], []
     for point in gcp:
         x = multiview.triangulate_gcp(
             point,
             reconstruction.shots,
+            threshold,
         )
         if x is not None and len(point.lla):
             point_enu = np.array(

opensfm/config.py

@@ -267,7 +267,7 @@ class OpenSfMConfig:
     # The default vertical standard deviation of the GCPs (in meters)
     gcp_vertical_sd: float = 0.1
     # Global weight for GCPs, expressed a ratio of the sum of (# projections) + (# shots) + (# relative motions)
-    gcp_global_weight: float = 0.01
+    gcp_global_weight: float = 0.1
     # The standard deviation of the rig translation
     rig_translation_sd: float = 0.1
     # The standard deviation of the rig rotation
@@ -283,6 +283,8 @@ class OpenSfMConfig:
     # Maximum optimizer iterations.
     bundle_max_iterations: int = 100
 
+    # Ratio of (resection candidates / total tracks) of a given image so that it is culled at resection and resected later
+    resect_redundancy_threshold: float = 0.7
     # Retriangulate all points from time to time
     retriangulation: bool = True
     # Retriangulate when the number of points grows by this ratio
@@ -299,6 +301,10 @@ class OpenSfMConfig:
     local_bundle_min_common_points: int = 20
     # Max number of shots to optimize during local bundle adjustment
     local_bundle_max_shots: int = 30
+    # Number of grid division for seleccting tracks in local bundle adjustment
+    local_bundle_grid: int = 8
+    # Number of grid division for selecting tracks in final bundle adjustment
+    final_bundle_grid: int = 32
 
     # Remove uncertain and isolated points from the final point cloud
     filter_final_point_cloud: bool = False
@@ -307,7 +313,7 @@ class OpenSfMConfig:
     save_partial_reconstructions: bool = False
 
     ##################################
-    # Params for GPS alignment
+    # Params for GPS/GCP alignment
     ##################################
     # Use or ignore EXIF altitude tag
     use_altitude_tag: bool = True
@@ -321,6 +327,8 @@ class OpenSfMConfig:
     bundle_use_gcp: bool = True
     # Compensate GPS with a per-camera similarity transform
     bundle_compensate_gps_bias: bool = False
+    # Thrershold for the reprojection error of GCPs to be considered outliers
+    gcp_reprojection_error_threshold: float = 0.05
 
     ##################################
     # Params for rigs

opensfm/dense.py

@@ -30,7 +30,7 @@ def compute_depthmaps(
     num_neighbors = config["depthmap_num_neighbors"]
 
     neighbors = {}
-    common_tracks = common_tracks_double_dict(graph)
+    common_tracks = common_tracks_double_dict(graph, processes)
     for shot in reconstruction.shots.values():
         neighbors[shot.id] = find_neighboring_images(
             shot, common_tracks, reconstruction, num_neighbors
@@ -394,13 +394,14 @@ def compute_depth_range(
 
 def common_tracks_double_dict(
     tracks_manager: pymap.TracksManager,
+    processes: int,
 ) -> Dict[str, Dict[str, List[str]]]:
     """List of track ids observed by each image pair.
 
     Return a dict, ``res``, such that ``res[im1][im2]`` is the list of
     common tracks between ``im1`` and ``im2``.
     """
-    common_tracks_per_pair = tracking.all_common_tracks_without_features(tracks_manager)
+    common_tracks_per_pair = tracking.all_common_tracks_without_features(tracks_manager, processes=processes)
     res = {image: {} for image in tracks_manager.get_shot_ids()}
     for (im1, im2), v in common_tracks_per_pair.items():
         res[im1][im2] = v

opensfm/exif.py

@@ -3,6 +3,7 @@
 import datetime
 import logging
 from codecs import decode, encode
+import math
 from typing import Any, BinaryIO, Callable, Dict, List, Optional, Tuple, Union
 
 import exifread
@@ -60,22 +61,52 @@ def get_tag_as_float(tags: Dict[str, Any], key: str, index: int = 0) -> Optional
         return None
 
 
+def focal35_to_focal_ratio(
+    focal35_or_ratio: float, width: int, height: int, inverse=False
+) -> float:
+    """
+    Convert focal length in 35mm film equivalent to focal ratio (and vice versa).
+    We follow https://en.wikipedia.org/wiki/35_mm_equivalent_focal_length
+    """
+    image_ratio = float(max(width, height)) / min(width, height)
+    is_32 = math.fabs(image_ratio - 3.0 / 2.0)
+    is_43 = math.fabs(image_ratio - 4.0 / 3.0)
+    if is_32 < is_43:
+        # 3:2 aspect ratio : use 36mm for 35mm film
+        film_width = 36.0
+        if inverse:
+            return focal35_or_ratio * film_width
+        else:
+            return focal35_or_ratio / film_width
+    else:
+        # 4:3 aspect ratio : use 34mm for 35mm film
+        film_width = 34
+        if inverse:
+            return focal35_or_ratio * film_width
+        else:
+            return focal35_or_ratio / film_width
+
+
 def compute_focal(
+    pixel_width: int,
+    pixel_height: int,
     focal_35: Optional[float],
     focal: Optional[float],
     sensor_width: Optional[float],
     sensor_string: Optional[str],
 ) -> Tuple[float, float]:
     if focal_35 is not None and focal_35 > 0:
-        focal_ratio = focal_35 / 36.0  # 35mm film produces 36x24mm pictures.
+        focal_ratio = focal35_to_focal_ratio(focal_35, pixel_width, pixel_height)
     else:
         if not sensor_width:
             sensor_width = (
                 sensor_data().get(sensor_string, None) if sensor_string else None
             )
         if sensor_width and focal:
             focal_ratio = focal / sensor_width
-            focal_35 = 36.0 * focal_ratio
+            focal_35 = focal35_to_focal_ratio(
+                focal, pixel_width, pixel_height, inverse=True
+            )
         else:
             focal_35 = 0.0
             focal_ratio = 0.0
@@ -248,7 +279,10 @@ def extract_projection_type(self) -> str:
 
     def extract_focal(self) -> Tuple[float, float]:
         make, model = self.extract_make(), self.extract_model()
+        width, height = self.extract_image_size()
         focal_35, focal_ratio = compute_focal(
+            width,
+            height,
             get_tag_as_float(self.tags, "EXIF FocalLengthIn35mmFilm"),
             get_tag_as_float(self.tags, "EXIF FocalLength"),
             self.extract_sensor_width(),
@@ -636,7 +670,13 @@ def extract_exif(self) -> Dict[str, Any]:
 
 def hard_coded_calibration(exif: Dict[str, Any]) -> Optional[Dict[str, Any]]:
     focal = exif["focal_ratio"]
-    fmm35 = int(round(focal * 36.0))
+    fmm35 = int(
+        round(
+            focal35_to_focal_ratio(
+                focal, int(exif["width"]), int(exif["height"]), inverse=True
+            )
+        )
+    )
     make = exif["make"].strip().lower()
     model = exif["model"].strip().lower()
     raw_calibrations = camera_calibration()[0]

opensfm/multiview.py

@@ -533,6 +533,7 @@ def triangulate_gcp(
     reproj_threshold: float = 0.02,
     min_ray_angle_degrees: float = 1.0,
     min_depth: float = 0.001,
+    iterations: int = 10,
 ) -> Optional[NDArray]:
     """Compute the reconstructed position of a GCP from observations."""
 
@@ -550,13 +551,17 @@ def triangulate_gcp(
 
     if len(os) >= 2:
         thresholds = len(os) * [reproj_threshold]
+        os = np.asarray(os)
+        bs = np.asarray(bs)
         valid_triangulation, X = pygeometry.triangulate_bearings_midpoint(
-            np.asarray(os),
-            np.asarray(bs),
+            os,
+            bs,
             thresholds,
             np.radians(min_ray_angle_degrees),
             min_depth,
         )
+
         if valid_triangulation:
+            X = pygeometry.point_refinement(os, bs, X, iterations)
             return X
     return None