feat: add new model estimators (plane, PnP, radial homography)

[f-dy]

New model estimators

Add Python bindings for models available in graph-cut-ransac but not previously exposed in pymagsac:

Function	Model	Min sample	Use case
findPlane3D	3D plane	3 points	Point cloud segmentation
findPnP	Camera pose (P3P)	3 2D-3D pairs	Visual localization, AR
findRadialHomography	Homography + radial distortion	5 correspondences	Wide-angle/fisheye cameras
findHomographyAffine	Homography from affine correspondences	2 affine pairs	Matching with oriented features
findFundamentalMatrixAffine	Fundamental from affine/SIFT	4 affine pairs	Feature-rich matching
findEssentialMatrixPlanar	Essential matrix (planar motion)	2 correspondences	Ground vehicles, indoor robots
findEssentialMatrixGravity	Essential matrix with gravity	3 correspondences	Phone/drone visual odometry

Fix: correct degrees of freedom for all MAGSAC estimators

The MAGSAC scoring hardcoded DoF=4 for all models. This is only correct for Sampson distance (fundamental/essential). All other models had wrong chi² tail shape, leading to suboptimal inlier/outlier discrimination.

Estimator	Residual	Correct DoF	Was
Fundamental/Essential	Sampson distance	4	4 ✓
Homography/PnP/Radial/Affine	2D transfer error	2	4 ✗
Rigid transformation	3D Euclidean distance	3	4 ✗
Line2D/Plane3D	Scalar signed distance	1	4 ✗

Implementation

The upper incomplete gamma Γ((DoF-1)/2, x) used in MAGSAC scoring is computed per-DoF:

• DoF=1: Precomputed E₁(x) table (gamma_values_dof1.cpp)
• DoF=2: √π · erfc(√x) — closed form via std::erfc
• DoF=3: e⁻ˣ — closed form via std::exp
• DoF=4: Existing stored_gamma_values[] table (unchanged)

The lower incomplete gamma γ((DoF+1)/2, x) used in MAGSAC++ scoring:

• DoF=1: 1 − e⁻ˣ — closed form
• DoF=2: (√π/2)·erf(√x) − √x·e⁻ˣ — closed form via std::erf
• DoF=3: 1 − (1+x)·e⁻ˣ — closed form
• DoF=4: Existing stored_lower_incomplete_gamma_values[] (unchanged)

Dispatch is via if constexpr with a static_assert in the else branch — adding a model with unsupported DoF triggers a compile-time error.

Backward compatibility

• DoF=4 path uses the exact same precomputed tables — numerically identical results for fundamental/essential/homography (homography was wrong before but existing users may rely on the behavior).
• Independent of fix: use per-estimator degrees of freedom in MAGSAC scoring graph-cut-ransac#50 (ships its own gamma_values_dof1.cpp; CMake detects if the submodule provides one).

Files changed

• src/pymagsac/include/estimators.h — new estimator wrappers + correct DoF/constants for all models
• src/pymagsac/include/magsac.h — DoF-dispatched gamma lookups in IRLS and MAGSAC++ scoring; empty-vector guard
• src/pymagsac/include/gamma_values_dof1.cpp — precomputed E₁(x) table for DoF=1
• src/pymagsac/src/magsac_python.cpp — new bindings + fix memory leaks (stack-allocate MAGSAC objects)
• src/pymagsac/include/magsac_python.hpp — declarations
• src/pymagsac/src/bindings.cpp — Python bindings
• CMakeLists.txt — detect gamma_values_dof1.cpp from submodule; exclude gamma tables from SRCS globs
• tests/test_plane3d.py, tests/test_pnp.py, etc. — unit tests

Additional bugfixes

• Segfault in solver_linear_model.h: when estimateModelNonminimal passes sample_=nullptr (meaning "use all normalized points") and sampleNumber_ == sampleSize(), the solver's fast-path dereferenced the null pointer. Fixed with && sample_ != nullptr guard. (This guard is the standalone gcransac PR fix: guard LinearModelSolver minimal fast-path against a null sample graph-cut-ransac#52)
• Memory leaks in magsac_python.cpp: all 13 new MAGSAC<...> calls had no corresponding delete. Replaced with stack allocation.
• tgamma(0) UB for DoF=1: guarded with if constexpr lambda to guarantee the call is never emitted.
• Empty-vector UB: &(sigma_inliers)[0] on potentially empty vectors replaced with .data() + explicit empty check before solver calls.
• Unit-dependent reference threshold for the new estimators: the MAGSAC reference threshold (interrupting_threshold) defaults to 1.0 — a pixel-unit assumption that drives the adaptive iteration count N = log(1−conf) / log(1 − (inliers/N)^sample_size). 1.0 is unprincipled even for pixels (SIFT at coarse scales has inlier localization error of several pixels) and unit-dependent for metric data. This PR sets it from the data for every new estimator it introduces — setReferenceThreshold(estimator.getSigmaQuantile() * sigma_max) (= the statistical inlier threshold τ = k·σ_max): findPlane3D, findPnP, findPnPAC, findPnPSIFT, findRadialHomography, findHomographyAffine, findFundamentalMatrixAffine, findEssentialMatrixPlanar, findEssentialMatrixGravity. Since these are new bindings there is no prior behavior to preserve. The same fix for the pre-existing estimators (rigid, fundamental, essential, line, homography) is a separate, independently-mergeable PR against master (it changes existing behavior and needs failure-rate benchmarking). Affects only the RANSAC search budget, not the per-hypothesis weight/loss.
• DoF=1 weight singularity (Bug A / Bug C): the E₁(x) table has E₁(0)=∞ at index 0. Small residuals (relative to σ_max) round to index 0 and produced infinite weights → the weighted LSQ received NaN/inf and findPlane3D returned no model (Bug A) or wandered (Bug C). Fixed by capping table[0] = E₁(step) (the first finite entry). Verified: plane fitting now recovers the truth to <0.001° / <0.1 mm, matching a reference numpy MAGSAC++.

Known gamma-weight issues (NOT fixed here — flagged for maintainers)

A review against the MAGSAC++ paper (Barath et al., CVPR 2020, Eq. 2.1–2.2) found two
pre-existing inaccuracies in the IRLS weight path (sigmaConsensusPlusPlus) that
also affect graph-cut-ransac's scoring_function.h. They are not fixed in this PR
because they change results for all models (including n=4) and need accuracy
benchmarks the repo lacks. The MAGSAC++ scoring function (getModelQualityPlusPlus)
is correct and matches the paper exactly.

• B2 — table-step/precision mismatch: the coarse table stores Γ(a, i/1000) (step 1/1000) but is indexed with precision_of_stored_gammas = 1e4, so the IRLS evaluates the gamma at 10× the intended argument.
• B3 — σ_max definition: the IRLS uses σ_max = maximum_threshold directly, whereas the paper (and the scoring function) use σ_max = maximum_threshold / k.
• B2 and B3 partially cancel (net factor 10/k² ≈ 0.75 for n=4), which is why MAGSAC++ works well on real data despite both. Fixing B2 alone (e.g. precision_of_stored_gammas = 1e3) is harmful — it removes the cancellation and makes the gamma argument 1/k² of correct (13× too small for n=4). B2 and B3 must be fixed together (use σ_max = threshold/k + the fine table). See the maintainer notes for the full analysis and numerical verification.

Reference threshold for the pre-existing estimators (separate PR)

The same interrupting_threshold = 1.0 default also affects the pre-existing estimators
(rigid, fundamental, essential, line, homography). Fixing those changes behavior existing
users rely on, so it is a separate, independently-mergeable PR against master
(needs failure-rate benchmarking on Adelaide/EVD/homogr). This PR only sets the reference
threshold for the new estimators it introduces (see the bullet above), where there is no
prior behavior to preserve.