Better handle outlier matches, bugfixes

Author: Phil26AT

gluefactory/geometry/depth.py

@@ -60,7 +60,6 @@ def project(
     camera_i,
     camera_j,
     T_itoj,
-    validi,
     ccth=None,
     sample_depth_fun=sample_depth,
     sample_depth_kwargs=None,
@@ -71,20 +70,22 @@ def project(
     kpi_3d_i = camera_i.image2cam(kpi)
     kpi_3d_i = kpi_3d_i * di[..., None]
     kpi_3d_j = T_itoj.transform(kpi_3d_i)
-    kpi_j, validj = camera_j.cam2image(kpi_3d_j)
+    kpi_j, valid = camera_j.cam2image(kpi_3d_j)
+    invalid = ~valid
     # di_j = kpi_3d_j[..., -1]
-    validi = validi & validj
     if depthj is None or ccth is None:
-        return kpi_j, validi & validj
+        return kpi_j, valid, invalid
     else:
         # circle consistency
         dj, validj = sample_depth_fun(kpi_j, depthj, **sample_depth_kwargs)
         kpi_j_3d_j = camera_j.image2cam(kpi_j) * dj[..., None]
         kpi_j_i, validj_i = camera_i.cam2image(T_itoj.inv().transform(kpi_j_3d_j))
-        consistent = ((kpi - kpi_j_i) ** 2).sum(-1) < ccth
-        visible = validi & consistent & validj_i & validj
+        reproj_error = ((kpi - kpi_j_i) ** 2).sum(-1)
+        consistent = reproj_error < ccth**2
+        visible = valid & consistent & validj_i & validj
+        invalid = invalid | (validj & ((~validj_i) | (~consistent)))
         # visible = validi
-        return kpi_j, visible
+        return kpi_j, visible, invalid
 
 
 def dense_warp_consistency(
@@ -100,7 +101,8 @@ def dense_warp_consistency(
         -2,
     )
     validi = di > 0
-    kpir, validir = project(kpi, di, depthj, camerai, cameraj, T_itoj, validi, **kwargs)
+    kpir, validir, _ = project(kpi, di, depthj, camerai, cameraj, T_itoj, **kwargs)
+    validir = validir & validi
 
     return kpir.unflatten(-2, depthi.shape[-2:]), validir.unflatten(
         -1, (depthi.shape[-2:])
@@ -120,12 +122,10 @@ def symmetric_reprojection_error(
     d0, valid0 = sample_depth(pts0, depth0)
     d1, valid1 = sample_depth(pts1, depth1)
 
-    pts0_1, visible0 = project(
-        pts0, d0, depth1, camera0, camera1, T_0to1, valid0, ccth=None
-    )
-    pts1_0, visible1 = project(
-        pts1, d1, depth0, camera1, camera0, T_1to0, valid1, ccth=None
-    )
+    pts0_1, visible0, _ = project(pts0, d0, depth1, camera0, camera1, T_0to1, ccth=None)
+    visible0 = visible0 & valid0
+    pts1_0, visible1, _ = project(pts1, d1, depth0, camera1, camera0, T_1to0, ccth=None)
+    visible1 = visible1 & valid1
 
     reprojection_errors_px = 0.5 * (
         (pts0_1 - pts1).norm(dim=-1) + (pts1_0 - pts0).norm(dim=-1)

gluefactory/geometry/epipolar.py

@@ -162,3 +162,37 @@ def relative_pose_error(T_0to1, R, t, ignore_gt_t_thr=0.0, eps=1e-10):
     r_err = angle_error_mat(R, R_gt)
 
     return t_err, r_err
+
+
+def check_epipolar_intersection(
+    x_i0: torch.Tensor,
+    i1_F_i0: torch.Tensor,
+    width1: int | torch.Tensor,
+    height1: int | torch.Tensor,
+) -> torch.BoolTensor:
+    x_i0 = tr.to_homogeneous(x_i0)  # (..., 3)
+    L_B = x_i0 @ i1_F_i0.T
+
+    l1, l2, l3 = L_B.split(1, dim=1)
+
+    eps = 1e-6
+
+    # Vertical boundary checks (x=0, x=W_B)
+    mask_v = torch.abs(l2) > eps
+    l2_inv = torch.where(mask_v, 1.0 / l2, torch.tensor(0.0, device=x_i0.device))
+
+    y0 = -l3 * l2_inv
+    yW = -(l1 * width1 + l3) * l2_inv
+
+    check_v = mask_v & (((y0 >= 0) & (y0 <= height1)) | ((yW >= 0) & (yW <= height1)))
+
+    # Horizontal boundary checks (y=0, y=H_B)
+    mask_h = torch.abs(l1) > eps
+    l1_inv = torch.where(mask_h, 1.0 / l1, torch.tensor(0.0, device=x_i0.device))
+
+    x0 = -l3 * l1_inv
+    xH = -(l2 * height1 + l3) * l1_inv
+
+    check_h = mask_h & (((x0 >= 0) & (x0 <= width1)) | ((xH >= 0) & (xH <= width1)))
+
+    return (check_v | check_h).squeeze(-1)

gluefactory/geometry/gt_generation.py

@@ -16,6 +16,7 @@ def gt_matches_from_pose_depth(
     epi_th=None,
     cc_th=None,
     min_overlap: float | None = None,
+    add_epi_outliers=True,
     **kw,
 ):
     if kp0.shape[1] == 0 or kp1.shape[1] == 0:
@@ -41,12 +42,32 @@ def gt_matches_from_pose_depth(
         d0, valid0 = depth.sample_depth(kp0, depth0)
         d1, valid1 = depth.sample_depth(kp1, depth1)
 
-    kp0_1, visible0 = depth.project(
-        kp0, d0, depth1, camera0, camera1, T_0to1, valid0, ccth=cc_th
+    kp0_1, visible0, unmatchable0 = depth.project(
+        kp0, d0, depth1, camera0, camera1, T_0to1, ccth=cc_th
     )
-    kp1_0, visible1 = depth.project(
-        kp1, d1, depth0, camera1, camera0, T_1to0, valid1, ccth=cc_th
+    visible0 = visible0 & valid0
+    kp1_0, visible1, unmatchable1 = depth.project(
+        kp1, d1, depth0, camera1, camera0, T_1to0, ccth=cc_th
     )
+    visible1 = visible1 & valid1
+
+    unmatchable0 = valid0 & unmatchable0
+    unmatchable1 = valid1 & unmatchable1
+
+    if add_epi_outliers:
+        i1_F_i0 = epipolar.T_to_F(camera0, camera1, T_0to1)
+        image_size0 = data["view0"]["image_size"]
+        image_size1 = data["view1"]["image_size"]
+
+        evalid0 = torch.vmap(
+            epipolar.check_epipolar_intersection,
+        )(kp0, i1_F_i0, image_size1[:, 0], image_size1[:, 1])
+        evalid1 = torch.vmap(
+            epipolar.check_epipolar_intersection,
+        )(kp1, i1_F_i0.transpose(-1, -2), image_size0[:, 0], image_size0[:, 1])
+
+        unmatchable0 = unmatchable0 | (~evalid0)
+        unmatchable1 = unmatchable1 | (~evalid1)
     if min_overlap is not None and "overlap_0to1" in data:
         has_overlap = (
             torch.max(data["overlap_0to1"], data["overlap_1to0"]) > min_overlap
@@ -118,6 +139,8 @@ def gt_matches_from_pose_depth(
         "proj_1to0": kp1_0,
         "visible0": visible0,
         "visible1": visible1,
+        "unmatchable0": unmatchable0,
+        "unmatchable1": unmatchable1,
         "has_overlap": has_overlap,
         "xyz_keypoints0": c0_t_w.inv() @ (camera0.image2cam(kp0) * d0.unsqueeze(-1)),
         "xyz_keypoints1": c1_t_w.inv() @ (camera1.image2cam(kp1) * d1.unsqueeze(-1)),

gluefactory/utils/experiments.py

@@ -72,6 +72,7 @@ def compose_config(
     default_config_dir: str = "configs/",
     overrides: Optional[list[str]] = None,
     sweep_idx: int | None = None,
+    resolve: bool = True,
 ) -> tuple[Path, OmegaConf]:
 
     conf_path = parse_config_path(name_or_path, default_config_dir)
@@ -88,6 +89,8 @@ def compose_config(
         OmegaConf.set_struct(custom_conf, False)
         custom_conf = OmegaConf.merge(custom_conf, sweep_conf)
         del custom_conf["sweep"]
+    if resolve:
+        OmegaConf.resolve(custom_conf)
     return conf_path, custom_conf
 
 

gluefactory/utils/metrics.py

@@ -48,5 +48,7 @@ def ranking_ap(m, gt_m, scores):
         f"{prefix}average_precision": ap,
         "num_matchable": (data[f"gt_{prefix_gt}matches0"] > -1).sum(1),
         "num_unmatchable": (data[f"gt_{prefix_gt}matches0"] == -1).sum(1),
+        "num_matches": (pred[f"{prefix}matches0"] > -1).sum(1),
+        "average_match_score": pred[f"{prefix}matching_scores0"].mean(1),
     }
     return metrics

gluefactory/utils/misc.py

@@ -471,6 +471,13 @@ def resize_image(
     return image_out
 
 
+def l2_normalize(
+    tensor: torch.Tensor, dim: int = -1, eps: float = 1e-10
+) -> torch.Tensor:
+    norm = torch.norm(tensor, p=2, dim=dim, keepdim=True).clamp_min(eps)
+    return tensor / norm
+
+
 def is_image_of_shape(image: torch.Tensor, hw: tuple[int, int]) -> bool:
     h, w = hw
     return h in image.shape and w in image.shape
@@ -806,7 +813,7 @@ def interpolate_patches(
     is_chw: bool = False,
     align_corners: bool = False,
     padding_mode: str = "zeros",
-) -> tuple[torch.Tensor, torch.Tensor]:  # B x N x D x ps x ps, B x N x 2
+) -> tuple[torch.Tensor, torch.Tensor, torch.Tensor]:  # B x N x D x ps x ps, B x N x 2
     if not is_chw:
         features = chw_from_hwc(features)
     if normalize:
@@ -817,8 +824,8 @@ def interpolate_patches(
         (1, 1, ps, ps), device=features.device, dtype=features.dtype
     )
     p_xy = get_image_coords(dummy_patch)
-    corners = torch.round(pts_i - ps / 2 - 0.5)
-    p_xy_i = corners[:, :, None, None, :] + p_xy[:, None]
+    cxy_i = torch.round(pts_i - ps / 2 - 0.5)
+    p_xy_i = cxy_i[:, :, None, None, :] + p_xy[:, None]
     p_xy_n = normalize_coords(p_xy_i, features.shape[-2:])
     patches = torch.vmap(grid_sample, in_dims=(None, 1), out_dims=1)(
         features,
@@ -829,7 +836,32 @@ def interpolate_patches(
     )
     if not is_chw:
         patches = hwc_from_chw(patches)
-    return patches, corners
+    return patches, p_xy_n, cxy_i
+
+
+def patch_interpolate_points(
+    pts: torch.Tensor,  # B x N X 2
+    patches: torch.Tensor,  # B x N x D x ps x ps OR B x N x ps x ps x D
+    **kwargs,
+):
+    return torch.vmap(interpolate_points, in_dims=0, out_dims=0)(
+        pts[:, :, None],
+        patches,
+        **kwargs,
+    )[..., 0, :]
+
+
+def log_softmax(scores: torch.Tensor, dim: int | tuple = -1) -> torch.Tensor:
+    """Numerically stable log softmax."""
+    if isinstance(dim, int):
+        return torch.log_softmax(scores, dim=dim)
+    else:
+        last = tuple(range(-len(dim), 0))
+        scores = scores.moveaxis(dim, last)
+        log_probs = torch.log_softmax(scores.flatten(-len(dim)), dim=-1).reshape(
+            *scores.shape
+        )
+        return log_probs.moveaxis(last, dim)
 
 
 def interpolate_matches(

gluefactory/visualization/viz2d.py

@@ -164,7 +164,7 @@ def plot_image_grid(
         return axs
 
 
-def plot_keypoints(kpts, colors="lime", ps=4, axes=None, a=1.0):
+def plot_keypoints(kpts, colors="lime", ps=4, axes=None, a=1.0, **kwargs):
     """Plot keypoints for existing images.
     Args:
         kpts: list of ndarrays of size (N, 2).
@@ -180,7 +180,7 @@ def plot_keypoints(kpts, colors="lime", ps=4, axes=None, a=1.0):
     for ax, k, c, alpha in zip(axes, kpts, colors, a):
         if isinstance(k, torch.Tensor):
             k = k.detach().cpu().numpy()
-        ax.scatter(k[:, 0], k[:, 1], c=c, s=ps, linewidths=0, alpha=alpha)
+        ax.scatter(k[:, 0], k[:, 1], c=c, s=ps, linewidths=0, alpha=alpha, **kwargs)
 
 
 def plot_matches(kpts0, kpts1, color=None, lw=1.5, ps=4, a=1.0, labels=None, axes=None):