colmap / pycolmap

Python bindings for COLMAP
BSD 3-Clause "New" or "Revised" License
858 stars 125 forks source link

estimate_triangulation result not equal to incremental_mapper result #269

Closed oneleggedredcow closed 3 months ago

oneleggedredcow commented 3 months ago

I'm trying to use estimate_triangulation to triangulate a set of points. To test that I understood how to use it correctly, I took the results from incremental_mapping and used that as input to estimate_triangulation. Unfortunately, I'm unable to obtain the same result from estimate_triangulation as I get from incremental_mapping.

...
result = pycolmap.incremental_mapping(...)

cams = list(result[0].cameras.values())
imgs = [next(img for img in result[0].images.values() if img.camera_id == cam.camera_id) for cam in cams]
random_pt3d = [x for x in result[0].points3D.values() if x.track.length() == 4][0]
pts_2d = {
    result[0].images[track.image_id].camera_id: result[0].images[track.image_id].points2D[track.point2D_idx]
    for track in random_pt3d.track.elements
}
tri_pts = [pycolmap.PointData(pts_2d[cam.camera_id].xy, cam.cam_from_img(pts_2d[cam.camera_id].xy)) for cam in cams]
tri_pt3d = pycolmap.estimate_triangulation(tri_pts, imgs, cams)

print(f"3D Point: incremental_mapping = {random_pt3d.xyz}, estimate_triangulation = {tri_pt3d['xyz']} -> diff = {random_pt3d.xyz - tri_pt3d['xyz']}")
# 3D Point: incremental_mapping = [ 8.03032441 -0.35841519 38.59303657], estimate_triangulation = [ 8.06614087 -0.37813197 38.73056718] -> diff = [-0.03581646  0.01971678 -0.13753062]

reproj_pts = [cam.img_from_cam(img.cam_from_world * tri_pt3d["xyz"]) for cam, img in zip(cams, imgs)]
reproj_error = np.mean([np.linalg.norm(x - y) for x, y in zip([pts_2d[cam.camera_id].xy for cam in cams], reproj_pts)])

print(f"Reprojection Error: incremental_mapping = {random_pt3d.error}, estimate_triangulation {reproj_error} -> diff = {random_pt3d.error - reproj_error}")
# Reprojection Error: incremental_mapping = 1.730548064280013, estimate_triangulation 1.8416187791990386 -> diff = -0.11107071491902554

The reprojection error from the estimate_triangulation point is slightly higher than the error from incremental_mapping point. Why do they not produce the same result?

I thought it might have to do with the camera_model so I switched all of the camera models to SIMPLE_PINHOLE and I get the same result.

I tried triangulating the points using a different approach (https://gist.github.com/davegreenwood/e1d2227d08e24cc4e353d95d0c18c914) and that gives me a very similar answer to estimate_triangulation. Is there something extra that I need to do in order to get the more accurate result that incremental_mapping provides?

sarlinpe commented 3 months ago

Yes: incremental_mapping subsequently refines the points and camera poses with bundle adjustment: https://github.com/colmap/colmap/blob/06a805386afb3db9dc6704c7ba2a80861bb64f7e/src/colmap/sfm/incremental_mapper.cc#L572-L577

You can accomplish this in Python using pyceres, see the BA example (section Optimize 3D points). Alternatively we could bind colmap/estimators/bundle_adjustment.h.

oneleggedredcow commented 3 months ago

Awesome, thank you!