fix: correct MAGSAC++ IRLS gamma weight (B2 + B3) in getScore

scripts/generate_gamma_values.py

@@ -0,0 +1,69 @@
+#!/usr/bin/env python3
+"""Generate src/pygcransac/include/gamma_values.cpp — the DoF=4 MAGSAC++ gamma
+lookup tables, trimmed to the useful argument interval.
+
+How it is generated: entry x holds the incomplete gamma at argument
+t = x / PRECISION (PRECISION = 1e4, i.e. step 1e-4), via scipy.special:
+  stored_complete_gamma_values[x]         = Gamma(1.5, t)   (upper incomplete; weight)
+  stored_lower_incomplete_gamma_values[x] = gamma(2.5, t)   (lower incomplete; loss)
+for DoF=4 (a_upper=(n-1)/2=1.5, a_lower=(n+1)/2=2.5).
+
+Why this length: MAGSAC++ truncates the weight at the inlier cutoff r <= k*sigma_max,
+so the argument t = r^2/(2 sigma_max^2) never exceeds k^2/2, with
+k = sqrt(chi2_0.99(4)). The tables therefore cover exactly [0, k^2/2]; the lookup
+clamps x to stored_incomplete_gamma_number, so any residual at/above the cutoff maps
+to the last entry and nothing beyond k^2/2 is ever read. (Was [0, 10] = 100001 entries.)
+
+Usage:
+  python3 generate_gamma_values.py > ../src/pygcransac/include/gamma_values.cpp
+"""
+import math
+from scipy.special import gammaincc, gammainc
+from scipy.stats import chi2
+
+PRECISION = 10000          # x = round(PRECISION * t); step = 1 / PRECISION = 1e-4
+DOF = 4
+A_UPPER = (DOF - 1) / 2.0  # 1.5
+A_LOWER = (DOF + 1) / 2.0  # 2.5
+PER_LINE = 100
+
+k = math.sqrt(chi2.ppf(0.99, DOF))
+max_t = k * k / 2.0
+MAX_INDEX = round(PRECISION * max_t)   # last reachable index given r <= k*sigma_max
+N = MAX_INDEX + 1
+
+g_upper = math.gamma(A_UPPER)
+g_lower = math.gamma(A_LOWER)
+
+
+def emit_array(name, fn):
+    print(f"constexpr double {name}[] = {{")
+    rows = []
+    for start in range(0, N, PER_LINE):
+        rows.append("\t" + ", ".join(repr(float(fn(i / PRECISION))) for i in range(start, min(start + PER_LINE, N))))
+    print(",\n".join(rows))
+    print("};\n")
+
+
+print("// AUTO-GENERATED by scripts/generate_gamma_values.py — do not edit by hand.")
+print(f"// DoF={DOF} MAGSAC++ gamma lookup tables, indexed by x = round({PRECISION} * t),")
+print(f"// i.e. argument t = x / {PRECISION} (step {1.0/PRECISION:g}). Computed with scipy.special:")
+print(f"//   stored_complete_gamma_values[x]         = Gamma({A_UPPER}, t)   (upper incomplete; IRLS/scoring weight)")
+print(f"//   stored_lower_incomplete_gamma_values[x] = gamma({A_LOWER}, t)   (lower incomplete; MAGSAC++ loss)")
+print("//")
+print("// Length rationale: MAGSAC++ truncates the weight at the inlier cutoff")
+print(f"//   r <= k*sigma_max,  so t = r^2/(2*sigma_max^2) never exceeds k^2/2.")
+print(f"// The table is sized for the PRECISE k = sqrt(chi2_0.99({DOF})) = {k!r}")
+print(f"// (k^2/2 = {max_t!r}), which is the largest plausible cutoff. It therefore")
+print("// also covers consumers that use a coarser/imprecise k (e.g. the literal 3.64")
+print("// currently in scoring_function.h: k^2/2 = 6.6248 < the size below). Lookups")
+print("// clamp x to stored_incomplete_gamma_number, so the exact-cutoff index and any")
+print("// residual at/above the cutoff map to the last entry; nothing beyond is read.")
+print(f"// Indices 0..{MAX_INDEX} ({N} entries); was [0, 10] = 100001 entries.")
+print("#pragma once")
+print()
+print(f"constexpr unsigned int stored_incomplete_gamma_number = {MAX_INDEX};")
+print("constexpr double precision_of_stored_incomplete_gammas = 1e4;")
+print()
+emit_array("stored_complete_gamma_values", lambda t: g_upper * gammaincc(A_UPPER, t))
+emit_array("stored_lower_incomplete_gamma_values", lambda t: g_lower * gammainc(A_LOWER, t))