# Copyright (c), ETH Zurich and UNC Chapel Hill. # All rights reserved. # # Redistribution and use in source and binary forms, with or without # modification, are permitted provided that the following conditions are met: # # * Redistributions of source code must retain the above copyright # notice, this list of conditions and the following disclaimer. # # * Redistributions in binary form must reproduce the above copyright # notice, this list of conditions and the following disclaimer in the # documentation and/or other materials provided with the distribution. # # * Neither the name of ETH Zurich and UNC Chapel Hill nor the names of # its contributors may be used to endorse or promote products derived # from this software without specific prior written permission. # # THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" # AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE # IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE # ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDERS OR CONTRIBUTORS BE # LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR # CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF # SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS # INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN # CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) # ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE # POSSIBILITY OF SUCH DAMAGE. import argparse import collections import copy import dataclasses import datetime import functools import multiprocessing import shutil import subprocess from abc import ABC, abstractmethod from pathlib import Path import numpy as np import pycolmap @dataclasses.dataclass(kw_only=True) class SceneInfo: # Dataset name. dataset: str # Category name. category: str # Scene name. scene: str # Path to the workspace directory in the run directory. workspace_path: Path # Path to the input images. image_path: Path # Path to the ground-truth sparse reconstruction. sparse_gt_path: Path # Whether to update camera priors from the ground-truth reconstruction. camera_priors_from_sparse_gt: bool # Additional arguments for the COLMAP reconstruction command. colmap_extra_args: list[str] @dataclasses.dataclass(kw_only=True) class SceneResult: # Scene information for which the result was computed. scene_info: SceneInfo # Flat list of errors. errors: np.ndarray # Number of images in the scene. num_images: int # Number of registered images in the scene (over all components). num_reg_images: int # Number of components in the scene. num_components: int # Number of images in the largest component. largest_component: int @dataclasses.dataclass(kw_only=True) class SceneMetrics: # Area under the curve (AUC) scores at specified error thresholds. aucs: np.ndarray error_thresholds: np.ndarray error_type: str # Number of images in the scene. num_images: int # Number of registered images in the scene (over all components). num_reg_images: int # Number of components in the scene. num_components: int # Number of images in the largest component. largest_component: int class Dataset(ABC): @property @abstractmethod def position_accuracy_gt(self) -> float: """Ground-truth position accuracy in meters.""" pass @abstractmethod def list_scenes(self) -> list[SceneInfo]: """List all scenes to evaluate.""" pass @abstractmethod def prepare_scene(self, scene_info: SceneInfo) -> None: """Prepare the scene for reconstruction.""" pass def parse_args() -> argparse.Namespace: datetime_str = datetime.datetime.now().strftime("%Y-%m-%d_%H-%M-%S") parser = argparse.ArgumentParser() parser.add_argument( "--data_path", default=Path(__file__).parent / "data", type=Path ) parser.add_argument( "--datasets", nargs="+", default=["eth3d", "blended-mvs", "imc2023", "imc2024"], ) parser.add_argument( "--categories", nargs="+", default=[], help="Categories to evaluate, if empty all categories are evaluated.", ) parser.add_argument( "--scenes", nargs="+", default=[], help="Scenes to evaluate, if empty all scenes are evaluated.", ) parser.add_argument( "--run_path", default=Path(__file__).parent / "runs", type=Path ) parser.add_argument("--run_name", default=datetime_str) parser.add_argument("--report_name", default=f"report-{datetime_str}") parser.add_argument( "--overwrite_database", default=False, action="store_true" ) parser.add_argument( "--overwrite_matches", default=False, action="store_true" ) parser.add_argument( "--overwrite_reconstruction", default=False, action="store_true" ) parser.add_argument( "--overwrite_alignment", default=False, action="store_true" ) parser.add_argument("--colmap_path", required=True) parser.add_argument("--use_gpu", default=True, action="store_true") parser.add_argument("--use_cpu", dest="use_gpu", action="store_false") parser.add_argument( "--parallelism", type=int, default=multiprocessing.cpu_count(), help="Number of processes for parallel reconstruction.", ) parser.add_argument("--quality", default="high") parser.add_argument( "--error_type", default="relative", choices=["relative", "absolute"], help="Whether to evaluate relative pairwise pose errors in angular " "distance or absolute pose errors through GT alignment.", ) parser.add_argument( "--rel_error_thresholds", type=float, nargs="+", default=[0.5, 1, 5, 10], help="Evaluation thresholds in degrees.", ) parser.add_argument( "--abs_error_thresholds", type=float, nargs="+", default=[0.02, 0.05, 0.2, 0.5], help="Evaluation thresholds in meters.", ) args = parser.parse_args() args.colmap_path = Path(args.colmap_path).resolve() if args.overwrite_database: pycolmap.logging.info( "Overwriting database also overwrites reconstruction" ) args.overwrite_reconstruction = True if args.overwrite_matches: pycolmap.logging.info( "Overwriting matches also overwrites reconstruction" ) args.overwrite_reconstruction = True if args.overwrite_reconstruction: pycolmap.logging.info( "Overwriting reconstruction also overwrites alignment" ) args.overwrite_alignment = True return args def update_camera_priors_from_sparse_gt( database_path: Path, camera_priors_sparse_gt: pycolmap.Reconstruction ) -> None: pycolmap.logging.info("Setting prior cameras from GT") database = pycolmap.Database() database.open(database_path) camera_id_gt_to_camera_id = {} for camera_id_gt, camera_gt in camera_priors_sparse_gt.cameras.items(): camera_gt.has_prior_focal_length = True camera_id = database.write_camera(camera_gt) camera_id_gt_to_camera_id[camera_id_gt] = camera_id images_gt_by_name = {} for image_gt in camera_priors_sparse_gt.images.values(): images_gt_by_name[image_gt.name] = image_gt for image in database.read_all_images(): if image.name not in images_gt_by_name: pycolmap.logging.warning( f"Not setting prior camera for image {image.name}, " "because it does not exist in GT" ) continue image_gt = images_gt_by_name[image.name] camera_id = camera_id_gt_to_camera_id[image_gt.camera_id] image.camera_id = camera_id database.update_image(image) database.close() def colmap_reconstruction( args: argparse.Namespace, workspace_path: Path, image_path: Path, camera_priors_sparse_gt: pycolmap.Reconstruction = None, colmap_extra_args: list = None, num_threads: int = 1, ) -> None: workspace_path.mkdir(parents=True, exist_ok=True) database_path = workspace_path / "database.db" if args.overwrite_database and database_path.exists(): database_path.unlink() sparse_path = workspace_path / "sparse" if args.overwrite_reconstruction and sparse_path.exists(): shutil.rmtree(sparse_path) if sparse_path.exists(): pycolmap.logging.info("Skipping reconstruction, as it already exists") return if args.overwrite_matches: subprocess.check_call( [ args.colmap_path, "database_cleaner", "--database_path", database_path, "--type", "matches", ], cwd=workspace_path, ) # TODO: Expose automatic reconstruction through pycolmap bindings instead # of using the command line interface. One blocker for this is that we # currently do not produce CUDA enabled pycolmap packages. colmap_args = [ args.colmap_path, "automatic_reconstructor", "--image_path", image_path, "--workspace_path", workspace_path, "--vocab_tree_path", args.data_path / "vocab_tree_flickr100K_words256K.bin", "--use_gpu", "1" if args.use_gpu else "0", "--num_threads", str(num_threads), "--quality", args.quality, ] subprocess.check_call( colmap_args + (colmap_extra_args or []) + [ "--extraction", "1", "--matching", "0", "--sparse", "0", "--dense", "0", ], cwd=workspace_path, ) if camera_priors_sparse_gt is not None: update_camera_priors_from_sparse_gt( database_path, camera_priors_sparse_gt ) subprocess.check_call( colmap_args + (colmap_extra_args or []) + [ "--extraction", "0", "--matching", "1", "--sparse", "0", "--dense", "0", ], cwd=workspace_path, ) # Decouple matching from sparse reconstruction, because matching will # initialize an OpenGL context and Mac on Apple silicon tends to assign GUI # applications to the low efficiency cores but we want to use the # performance cores. subprocess.check_call( colmap_args + (colmap_extra_args or []) + [ "--extraction", "0", "--matching", "0", "--sparse", "1", "--dense", "0", ], cwd=workspace_path, ) def colmap_alignment( args: argparse.Namespace, sparse_path: Path, sparse_gt_path: Path, sparse_aligned_path: Path, max_ref_model_error: float, ) -> None: if args.overwrite_alignment and sparse_aligned_path.exists(): shutil.rmtree(sparse_aligned_path) if sparse_aligned_path.exists(): pycolmap.logging.info("Skipping alignment, as it already exists") return if sparse_path.exists(): sparse_aligned_path.mkdir(parents=True, exist_ok=True) subprocess.call( [ args.colmap_path, "model_aligner", "--input_path", sparse_path, "--ref_model_path", sparse_gt_path, "--output_path", sparse_aligned_path, "--alignment_max_error", str(max_ref_model_error), ] ) def process_scene( args: argparse.Namespace, scene_info: SceneInfo, prepare_scene: callable, position_accuracy_gt: float, num_threads: int, ) -> SceneResult: pycolmap.logging.info( f"Processing dataset={scene_info.dataset}, " f"category={scene_info.category}, " f"scene={scene_info.scene}" ) prepare_scene(scene_info) sparse_gt = pycolmap.Reconstruction(scene_info.sparse_gt_path) colmap_reconstruction( args=args, workspace_path=scene_info.workspace_path, image_path=scene_info.image_path, camera_priors_sparse_gt=( sparse_gt if scene_info.camera_priors_from_sparse_gt else None ), num_threads=num_threads, colmap_extra_args=scene_info.colmap_extra_args, ) # Merge all sub-models into a single reconstruction. Each sub-model will be # "randomly" aligned to the other sub-models. We then compute the overall # error over the merged reconstruction. With this simple appraoch, there is # a small chance that the randomly aligned images in one sub-model are # correctly aligned with other sub-models and the error is therefore # underestimated. However, this is very unlikely to happen. sparse_merged = pycolmap.Reconstruction() num_components = 0 largest_component = 0 for sparse_path in (scene_info.workspace_path / "sparse").iterdir(): if not sparse_path.is_dir(): continue num_components += 1 sparse = None if args.error_type == "relative": sparse = pycolmap.Reconstruction(sparse_path) elif args.error_type == "absolute": sparse_aligned_path = scene_info.workspace_path / "sparse_aligned" colmap_alignment( args=args, sparse_path=sparse_path, sparse_gt_path=scene_info.sparse_gt_path, sparse_aligned_path=sparse_aligned_path, max_ref_model_error=position_accuracy_gt, ) if (sparse_aligned_path / "images.bin").exists(): sparse = pycolmap.Reconstruction(sparse_aligned_path) else: raise ValueError(f"Invalid error type: {args.error_type}") if sparse is not None: largest_component = max(largest_component, sparse.num_images()) for image in sparse.images.values(): if image.image_id in sparse_merged.images: continue if image.camera_id not in sparse_merged.cameras: sparse_merged.add_camera(image.camera) if image.frame_id not in sparse_merged.frames: if image.frame.rig_id not in sparse_merged.rigs: sparse_merged.add_rig(image.frame.rig) image.frame.reset_rig_ptr() sparse_merged.add_frame(image.frame) image.reset_camera_ptr() image.reset_frame_ptr() sparse_merged.add_image(image) if args.error_type == "relative": dts, dRs = compute_rel_errors( sparse_gt=sparse_gt, sparse=sparse_merged, min_proj_center_dist=position_accuracy_gt, ) errors = np.maximum(dts, dRs) elif args.error_type == "absolute": dts, dRs = compute_abs_errors( sparse_gt=sparse_gt, sparse=sparse_merged, ) errors = dts else: raise ValueError(f"Invalid error type: {args.error_type}") return SceneResult( scene_info=scene_info, errors=errors, num_images=sparse_gt.num_images(), num_reg_images=sparse_merged.num_images(), num_components=num_components, largest_component=largest_component, ) def process_scenes( args: argparse.Namespace, scene_infos: list[SceneInfo], prepare_scene: callable, position_accuracy_gt: float, ) -> dict[str, dict[str, SceneMetrics]]: error_thresholds = get_error_thresholds(args) num_threads = min( args.parallelism, 2 * max(1, int(args.parallelism / len(scene_infos))) ) with multiprocessing.Pool(processes=args.parallelism) as p: results = p.map( functools.partial( process_scene, args, prepare_scene=prepare_scene, position_accuracy_gt=position_accuracy_gt, num_threads=num_threads, ), scene_infos, ) metrics = collections.defaultdict(dict) errors_by_category = collections.defaultdict(list) total_num_images = 0 total_num_reg_images = 0 total_num_components = 0 total_largest_components = 0 num_scenes = len(results) for result in results: errors_by_category[result.scene_info.category].extend(result.errors) total_num_images += result.num_images total_num_reg_images += result.num_reg_images total_num_components += result.num_components total_largest_components += result.largest_component metrics[result.scene_info.category][result.scene_info.scene] = ( SceneMetrics( aucs=compute_auc( result.errors, error_thresholds, min_error=position_accuracy_gt, ), error_thresholds=error_thresholds, error_type=args.error_type, num_images=result.num_images, num_reg_images=result.num_reg_images, num_components=result.num_components, largest_component=result.largest_component, ) ) for category, errors in errors_by_category.items(): metrics[category]["__all__"] = SceneMetrics( aucs=compute_auc( errors, error_thresholds, min_error=position_accuracy_gt, ), error_thresholds=error_thresholds, error_type=args.error_type, num_images=total_num_images, num_reg_images=total_num_reg_images, num_components=total_num_components, largest_component=total_largest_components, ) metrics[category]["__avg__"] = SceneMetrics( aucs=compute_avg_auc(metrics[category]), error_thresholds=error_thresholds, error_type=args.error_type, num_images=int(round(total_num_images / num_scenes)), num_reg_images=int(round(total_num_reg_images / num_scenes)), num_components=int(round(total_num_components / num_scenes)), largest_component=int(round(total_largest_components / num_scenes)), ) return metrics def normalize_vec(vec: np.ndarray, eps: float = 1e-10) -> np.ndarray: return vec / max(eps, np.linalg.norm(vec)) def vec_angular_dist_deg(vec1: np.ndarray, vec2: np.ndarray) -> float: cos_dist = np.clip(np.dot(normalize_vec(vec1), normalize_vec(vec2)), -1, 1) return np.rad2deg(np.acos(cos_dist)) def get_error_thresholds(args: argparse.Namespace) -> list[float]: if args.error_type == "relative": return np.array(args.rel_error_thresholds) elif args.error_type == "absolute": return np.array(args.abs_error_thresholds) else: raise ValueError(f"Invalid error type: {args.error_type}") def compute_rel_errors( sparse_gt: pycolmap.Reconstruction, sparse: pycolmap.Reconstruction, min_proj_center_dist: float, ) -> tuple[np.ndarray, np.ndarray]: """Computes angular relative pose errors across all image pairs. Notice that this approach leads to a super-linear decrease in the AUC scores when multiple images fail to register. Consider that we have N images in total in a dataset and M images are registered in the evaluated reconstruction. In this case, we can compute "finite" errors for (N-M)^2 pairs while the dataset has a total of N^2 pairs. In case of many unregistered images, the AUC score will drop much more than the (intuitively) expected (N-M) / N ratio. One could appropriately normalize by computing a single score per image through a suitable normalization of all pairwise errors per image. However, this becomes difficult when multiple sub-components are incorrectly stitched together in the same reconstruction (e.g., in the case of symmetry issues). """ if sparse is None: pycolmap.logging.error("Reconstruction failed") return len(sparse_gt.images) * [np.inf], len(sparse_gt.images) * [180] images = {} for image in sparse.images.values(): images[image.name] = image dts = [] dRs = [] for this_image_gt in sparse_gt.images.values(): if this_image_gt.name not in images: for _ in range(sparse_gt.num_images() - 1): dts.append(np.inf) dRs.append(180) continue this_image = images[this_image_gt.name] for other_image_gt in sparse_gt.images.values(): if this_image_gt.image_id == other_image_gt.image_id: continue if other_image_gt.name not in images: dts.append(np.inf) dRs.append(180) continue other_image = images[other_image_gt.name] other_from_this = ( other_image.cam_from_world() * this_image.cam_from_world().inverse() ) other_from_this_gt = ( other_image_gt.cam_from_world() * this_image_gt.cam_from_world().inverse() ) estimated_from_gt = other_from_this.inverse() * other_from_this_gt if ( np.linalg.norm(other_from_this_gt.translation) < min_proj_center_dist ): # If the cameras almost coincide, then the angular direction # distance is unstable, because a small position change can # cause a large rotational error. In this case, we only measure # rotational relative pose error. dt = 0 else: dt = vec_angular_dist_deg( other_from_this.translation, other_from_this_gt.translation ) dR = np.rad2deg(estimated_from_gt.rotation.angle()) dts.append(dt) dRs.append(dR) return np.array(dts), np.array(dRs) def compute_abs_errors( sparse_gt: pycolmap.Reconstruction, sparse: pycolmap.Reconstruction ) -> tuple[np.ndarray, np.ndarray]: """Computes rotational and translational absolute pose errors. Assumes that the input reconstructions are aligned in the same coordinate system. Computes one error per ground-truth image. """ dts = np.full(len(sparse_gt.images), fill_value=np.inf, dtype=np.float64) dRs = np.full(len(sparse_gt.images), fill_value=180, dtype=np.float64) if sparse is None: pycolmap.logging.error("Reconstruction or alignment failed") return dts, dRs images = {} for image in sparse.images.values(): images[image.name] = image dts = np.full(len(sparse_gt.images), fill_value=np.inf, dtype=np.float64) dRs = np.full(len(sparse_gt.images), fill_value=180, dtype=np.float64) for i, image_gt in enumerate(sparse_gt.images.values()): if image_gt.name not in images: continue image = images[image_gt.name] estimated_from_gt = ( image.cam_from_world() * image_gt.cam_from_world().inverse() ) dts[i] = np.linalg.norm(estimated_from_gt.translation) dRs[i] = np.rad2deg(estimated_from_gt.rotation.angle()) return dts, dRs def compute_recall(errors: np.ndarray) -> tuple[np.ndarray, np.ndarray]: num_elements = len(errors) errors = np.sort(errors) recall = (np.arange(num_elements) + 1) / num_elements return errors, recall def compute_auc( errors: np.ndarray, thresholds: list[float], min_error: float = 0 ) -> list[float]: if len(errors) == 0: raise ValueError("No errors to evaluate") errors, recall = compute_recall(errors) if min_error > 0: min_index = np.searchsorted(errors, min_error, side="right") min_score = min_index / len(errors) recall = np.r_[min_score, min_score, recall[min_index:]] errors = np.r_[0, min_error, errors[min_index:]] else: recall = np.r_[0, recall] errors = np.r_[0, errors] aucs = np.zeros(len(thresholds), dtype=np.float64) for i, t in enumerate(thresholds): last_index = np.searchsorted(errors, t, side="right") r = np.r_[recall[:last_index], recall[last_index - 1]] e = np.r_[errors[:last_index], t] auc = np.trapezoid(r, x=e) / t aucs[i] = auc * 100 return aucs / 1.1 def compute_avg_auc(scene_metrics: dict[str, SceneMetrics]) -> list[float]: auc_sum = None num_scenes = 0 for scene, metrics in scene_metrics.items(): if scene.startswith("__") and scene.endswith("__"): continue num_scenes += 1 if auc_sum is None: auc_sum = copy.copy(metrics.aucs) else: for i in range(len(auc_sum)): auc_sum[i] += metrics.aucs[i] return np.array(auc_sum) / num_scenes def diff_metrics( metrics_a: dict[str, dict[str, SceneMetrics]], metrics_b: dict[str, dict[str, SceneMetrics]], ): """Computes difference between two sets of metrics. Raises exception if the metrics are inconsistent. """ metrics_diff = copy.deepcopy(metrics_a) for dataset, category_metrics_a in metrics_a.items(): if dataset not in metrics_b: raise ValueError(f"Dataset {dataset} not found in metrics_b") category_metrics_b = metrics_b[dataset] for category, scene_metrics_a in category_metrics_a.items(): if category not in category_metrics_b: raise ValueError(f"Category {category} not found in metrics_b") scene_metrics_b = category_metrics_b[category] for scene, metrics_a in scene_metrics_a.items(): if scene not in scene_metrics_b: raise ValueError(f"Scene {scene} not found in metrics_b") metrics_b = scene_metrics_b[scene] if metrics_a.error_type != metrics_b.error_type or not np.all( metrics_a.error_thresholds == metrics_b.error_thresholds ): raise ValueError("Inconsistent error thresholds or types") metrics_diff[dataset][category][scene] = SceneMetrics( aucs=metrics_a.aucs - metrics_b.aucs, error_thresholds=metrics_a.error_thresholds, error_type=metrics_a.error_type, num_images=metrics_a.num_images - metrics_b.num_images, num_reg_images=metrics_a.num_reg_images - metrics_b.num_reg_images, num_components=metrics_a.num_components - metrics_b.num_components, largest_component=metrics_a.largest_component - metrics_b.largest_component, ) return metrics_diff def create_result_table( dataset_metrics: dict[str, dict[str, SceneMetrics]], ) -> str: first_metrics = next( iter(next(iter(next(iter(dataset_metrics.values())).values())).values()) ) if first_metrics.error_type == "relative": label = "AUC @ X deg (%)" thresholds = first_metrics.error_thresholds elif first_metrics.error_type == "absolute": label = "AUC @ X cm (%)" thresholds = 100 * first_metrics.error_thresholds else: raise ValueError(f"Invalid error type: {first_metrics.error_type}") column = "scenes" size_scenes = max( len(column) + 2, max( len(s) for d in dataset_metrics.values() for c in d.values() for s in c ), ) size_aucs = max(len(label) + 2, len(thresholds) * 7 - 1) size_imgs = 12 size_comps = 12 size_sep = size_scenes + size_aucs + size_imgs + size_comps + 3 header = ( f"{column:=^{size_scenes}} {label:=^{size_aucs}} " f"{'images':=^{size_imgs}} {'components':=^{size_comps}}" ) header += "\n" + " " * (size_scenes + 1) header += " ".join(f"{str(t).rstrip('.'):^6}" for t in thresholds) header += " reg all num largest" text = [header] for dataset, category_metrics in dataset_metrics.items(): for category, scene_metrics in category_metrics.items(): text.append(f"\n{dataset + '=' + category:=^{size_sep}}") for scene, metrics in scene_metrics.items(): assert len(metrics.aucs) == len(thresholds) row = "" if scene == "__avg__": scene = "average" row += "-" * size_sep + "\n" if scene == "__all__": scene = "overall" row += "-" * size_sep + "\n" row += f"{scene:<{size_scenes}} " row += " ".join(f"{auc:>6.2f}" for auc in metrics.aucs) row += f" {metrics.num_reg_images:6d}" row += f"{metrics.num_images:6d}" row += f" {metrics.num_components:4d}" row += f"{metrics.largest_component:8d}" text.append(row) return "\n".join(text)