Skip to content

Commit ec51503

Browse files
committed
fix: correct MAGSAC++ IRLS gamma weight (B2 + B3) in getScore; trim DoF=4 tables
getScore computed the per-point MAGSAC++ weight at the wrong gamma argument: - B3: sigma_max was increasedThreshold (the inlier threshold) rather than threshold/k. Now sigma_max = increasedThreshold/k drives both the argument and the prefactor; the inlier-collection cutoff stays at increasedThreshold = k*sigma_max. - B2: the index used the 1e4 precision against the coarse stored_gamma_values table (step 1e-3), reading the gamma at 10x the argument. Now it indexes the fine stored_complete_gamma_values table (step 1e-4) whose spacing matches 1e4. Tables (gamma_values.cpp) are now generated by scripts/generate_gamma_values.py (scipy.special) and trimmed to the useful interval [0, k^2/2]: 66385 entries (was [0,10]=100001). The dead coarse table is removed. The generated C++ documents how it was produced and why the length; the table is sized for the precise k = sqrt(chi2_0.99(4)) so it works with either the imprecise k literal kept here (3.64) or the precise per-DoF k introduced by danini#50. k stays at the existing imprecise literal in this PR (focused on B2/B3); precise per-DoF k is danini#50's job. Changes results for existing DoF=4 estimators: validate on Adelaide/EVD/homogr before submission.
1 parent 799a645 commit ec51503

3 files changed

Lines changed: 1436 additions & 2400 deletions

File tree

scripts/generate_gamma_values.py

Lines changed: 69 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,69 @@
1+
#!/usr/bin/env python3
2+
"""Generate src/pygcransac/include/gamma_values.cpp — the DoF=4 MAGSAC++ gamma
3+
lookup tables, trimmed to the useful argument interval.
4+
5+
How it is generated: entry x holds the incomplete gamma at argument
6+
t = x / PRECISION (PRECISION = 1e4, i.e. step 1e-4), via scipy.special:
7+
stored_complete_gamma_values[x] = Gamma(1.5, t) (upper incomplete; weight)
8+
stored_lower_incomplete_gamma_values[x] = gamma(2.5, t) (lower incomplete; loss)
9+
for DoF=4 (a_upper=(n-1)/2=1.5, a_lower=(n+1)/2=2.5).
10+
11+
Why this length: MAGSAC++ truncates the weight at the inlier cutoff r <= k*sigma_max,
12+
so the argument t = r^2/(2 sigma_max^2) never exceeds k^2/2, with
13+
k = sqrt(chi2_0.99(4)). The tables therefore cover exactly [0, k^2/2]; the lookup
14+
clamps x to stored_incomplete_gamma_number, so any residual at/above the cutoff maps
15+
to the last entry and nothing beyond k^2/2 is ever read. (Was [0, 10] = 100001 entries.)
16+
17+
Usage:
18+
python3 generate_gamma_values.py > ../src/pygcransac/include/gamma_values.cpp
19+
"""
20+
import math
21+
from scipy.special import gammaincc, gammainc
22+
from scipy.stats import chi2
23+
24+
PRECISION = 10000 # x = round(PRECISION * t); step = 1 / PRECISION = 1e-4
25+
DOF = 4
26+
A_UPPER = (DOF - 1) / 2.0 # 1.5
27+
A_LOWER = (DOF + 1) / 2.0 # 2.5
28+
PER_LINE = 100
29+
30+
k = math.sqrt(chi2.ppf(0.99, DOF))
31+
max_t = k * k / 2.0
32+
MAX_INDEX = round(PRECISION * max_t) # last reachable index given r <= k*sigma_max
33+
N = MAX_INDEX + 1
34+
35+
g_upper = math.gamma(A_UPPER)
36+
g_lower = math.gamma(A_LOWER)
37+
38+
39+
def emit_array(name, fn):
40+
print(f"constexpr double {name}[] = {{")
41+
rows = []
42+
for start in range(0, N, PER_LINE):
43+
rows.append("\t" + ", ".join(repr(float(fn(i / PRECISION))) for i in range(start, min(start + PER_LINE, N))))
44+
print(",\n".join(rows))
45+
print("};\n")
46+
47+
48+
print("// AUTO-GENERATED by scripts/generate_gamma_values.py — do not edit by hand.")
49+
print(f"// DoF={DOF} MAGSAC++ gamma lookup tables, indexed by x = round({PRECISION} * t),")
50+
print(f"// i.e. argument t = x / {PRECISION} (step {1.0/PRECISION:g}). Computed with scipy.special:")
51+
print(f"// stored_complete_gamma_values[x] = Gamma({A_UPPER}, t) (upper incomplete; IRLS/scoring weight)")
52+
print(f"// stored_lower_incomplete_gamma_values[x] = gamma({A_LOWER}, t) (lower incomplete; MAGSAC++ loss)")
53+
print("//")
54+
print("// Length rationale: MAGSAC++ truncates the weight at the inlier cutoff")
55+
print(f"// r <= k*sigma_max, so t = r^2/(2*sigma_max^2) never exceeds k^2/2.")
56+
print(f"// The table is sized for the PRECISE k = sqrt(chi2_0.99({DOF})) = {k!r}")
57+
print(f"// (k^2/2 = {max_t!r}), which is the largest plausible cutoff. It therefore")
58+
print("// also covers consumers that use a coarser/imprecise k (e.g. the literal 3.64")
59+
print("// currently in scoring_function.h: k^2/2 = 6.6248 < the size below). Lookups")
60+
print("// clamp x to stored_incomplete_gamma_number, so the exact-cutoff index and any")
61+
print("// residual at/above the cutoff map to the last entry; nothing beyond is read.")
62+
print(f"// Indices 0..{MAX_INDEX} ({N} entries); was [0, 10] = 100001 entries.")
63+
print("#pragma once")
64+
print()
65+
print(f"constexpr unsigned int stored_incomplete_gamma_number = {MAX_INDEX};")
66+
print("constexpr double precision_of_stored_incomplete_gammas = 1e4;")
67+
print()
68+
emit_array("stored_complete_gamma_values", lambda t: g_upper * gammaincc(A_UPPER, t))
69+
emit_array("stored_lower_incomplete_gamma_values", lambda t: g_lower * gammainc(A_LOWER, t))

0 commit comments

Comments
 (0)