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

.gitmodules

@@ -3,4 +3,4 @@
 	url = https://github.com/pybind/pybind11.git
 [submodule "graph-cut-ransac"]
 	path = graph-cut-ransac
-	url = https://github.com/danini/graph-cut-ransac
+	url = https://github.com/f-dy/graph-cut-ransac

CMakeLists.txt

@@ -75,6 +75,16 @@ include_directories (
 	${PROJECT_SOURCE_DIR}/graph-cut-ransac/src/pygcransac/include
 )
 
+# Check if the graph-cut-ransac submodule provides gamma_values_dof1.cpp (DoF=1 table).
+# If not, use the local copy in src/pymagsac/include/.
+set(GCRANSAC_DOF1_TABLE "${PROJECT_SOURCE_DIR}/graph-cut-ransac/src/pygcransac/include/gamma_values_dof1.cpp")
+if(EXISTS "${GCRANSAC_DOF1_TABLE}")
+	message(STATUS "Using gamma_values_dof1.cpp from graph-cut-ransac submodule")
+else()
+	message(STATUS "gamma_values_dof1.cpp not found in graph-cut-ransac; using local copy")
+	include_directories(BEFORE ${PROJECT_SOURCE_DIR}/src/pymagsac/include)
+endif()
+
 # ==============================================================================
 # Structure: Graph-Cut RANSAC Library
 # ==============================================================================
@@ -84,9 +94,11 @@ file(GLOB_RECURSE HDRS_GraphCutRANSAC
 ) 
 
 # Set source files to be added to the library
+# Note: gamma_values*.cpp are data tables #included by headers, not compiled directly
 file(GLOB_RECURSE SRCS_GraphCutRANSAC 
 	"graph-cut-ransac/src/pygcransac/include/*.cpp"
 )
+list(FILTER SRCS_GraphCutRANSAC EXCLUDE REGEX "gamma_values")
 
 add_library(GraphCutRANSAC STATIC
 	${HDRS_GraphCutRANSAC}
@@ -114,6 +126,7 @@ file(GLOB_RECURSE HDRS_MAGSAC
 file(GLOB_RECURSE SRCS_MAGSAC 
 	"${SOURCE_DIR}/include/*.cpp"
 )
+list(FILTER SRCS_MAGSAC EXCLUDE REGEX "gamma_values")
 
 # ==============================================================================
 # Structure: Python binding

src/pymagsac/include/estimators.h

@@ -5,6 +5,23 @@
 #include "estimators/fundamental_estimator.h"
 #include "estimators/homography_estimator.h"
 #include "estimators/linear_model_estimator.h"
+#include "estimators/perspective_n_point_estimator.h"
+#include "estimators/radial_homography_estimator.h"
+#include "estimators/homography_affine_correspondence_estimator.h"
+#include "estimators/solver_p3p.h"
+#include "estimators/solver_epnp_lm.h"
+#include "estimators/solver_pnp_bundle_adjustment.h"
+#include "estimators/solver_acp1p_cayley.h"
+#include "estimators/solver_acp1p_query.h"
+#include "estimators/solver_siftp1p.h"
+#include "estimators/solver_siftp2p.h"
+#include "estimators/solver_radial_homography_5pc.h"
+#include "estimators/solver_radial_homography_6pc.h"
+#include "estimators/solver_homography_two_affine.h"
+#include "estimators/solver_fundamental_matrix_four_sift.h"
+#include "estimators/solver_fundamental_matrix_three_affine.h"
+#include "estimators/solver_essential_matrix_two_points_planar.h"
+#include "estimators/solver_essential_matrix_three_points_gravity.h"
 #include "model.h"
 
 namespace magsac
@@ -35,34 +52,34 @@ namespace magsac
 
 			static constexpr double getSigmaQuantile()
 			{
-				return 3.64;
+				return 3.0348542587702925;
 			}
 
 			static constexpr size_t getDegreesOfFreedom()
 			{
-				return 4;
+				return 2;
 			}
 
 			static constexpr double getC()
 			{
-				return 0.25;
+				return 0.14104739588693907;
 			}
 
 			// Calculating the upper incomplete gamma value of (DoF - 1) / 2 with k^2 / 2.
 			static constexpr double getGammaFunction()
 			{
-				return 1.0;
+				return 1.7724538509055159;
 			}
 
 			static constexpr double getUpperIncompleteGammaOfK()
 			{
-				return 0.0036572608340910764;
+				return 0.0042654446820694645;
 			}
 
 			// Calculating the lower incomplete gamma value of (DoF + 1) / 2 with k^2 / 2.
 			static constexpr double getLowerIncompleteGammaOfK()
 			{
-				return 1.3012265540498875;
+				return 0.86263454284882968;
 			}
 
 			static constexpr double getChiSquareParamCp()
@@ -101,12 +118,12 @@ namespace magsac
 
 			static constexpr double getSigmaQuantile()
 			{
-				return 3.64;
+				return 3.3682141752187271;
 			}
 
 			static constexpr size_t getDegreesOfFreedom()
 			{
-				return 4;
+				return 3;
 			}
 
 			static constexpr double getGammaFunction()
@@ -122,13 +139,13 @@ namespace magsac
 			// Calculating the upper incomplete gamma value of (DoF - 1) / 2 with k^2 / 2.
 			static constexpr double getUpperIncompleteGammaOfK()
 			{
-				return 0.0036572608340910764;
+				return 0.0034394855607548570;
 			}
 
 			// Calculating the lower incomplete gamma value of (DoF + 1) / 2 with k^2 / 2.
 			static constexpr double getLowerIncompleteGammaOfK()
 			{
-				return 1.3012265540498875;
+				return 0.97705026178573529;
 			}
 
 			static constexpr double getChiSquareParamCp()
@@ -613,12 +630,12 @@ namespace magsac
 
 			static constexpr double getSigmaQuantile()
 			{
-				return 3.64;
+				return 2.5758293035489008;
 			}
 
 			static constexpr size_t getDegreesOfFreedom()
 			{
-				return 4;
+				return 1;
 			}
 
 			static constexpr double getC()
@@ -628,7 +645,7 @@ namespace magsac
 
 			static constexpr double getGammaFunction()
 			{
-				return 1.0;
+				return 1.0; // Not used for DoF=1 (Γ(0) diverges); kept for interface compatibility
 			}
 
 			static constexpr bool doesNormalizationForNonMinimalFitting()
@@ -639,20 +656,158 @@ namespace magsac
 			// Calculating the upper incomplete gamma value of (DoF - 1) / 2 with k^2 / 2.
 			static constexpr double getUpperIncompleteGammaOfK()
 			{
-				return 0.0036572608340910764;
+				return 0.0087462271161219563;
 			}
 
 			// Calculating the lower incomplete gamma value of (DoF + 1) / 2 with k^2 / 2.
 			static constexpr double getLowerIncompleteGammaOfK()
 			{
-				return 1.3012265540498875;
+				return 0.96375479928483043;
+			}
+
+			static constexpr double getChiSquareParamCp()
+			{
+				return 1.0 / (4.0 * getGammaFunction());
+			}
+		};
+
+		// This is the estimator class for estimating a 3D plane from a point cloud.
+		template<class _MinimalSolverEngine,  // The solver used for estimating the model from a minimal sample
+			class _NonMinimalSolverEngine> // The solver used for estimating the model from a non-minimal sample
+			class Plane3DEstimator :
+			public gcransac::estimator::LinearModelEstimator<_MinimalSolverEngine, _NonMinimalSolverEngine, 3>
+		{
+		public:
+			using gcransac::estimator::LinearModelEstimator<_MinimalSolverEngine, _NonMinimalSolverEngine, 3>::residual;
+
+			Plane3DEstimator() :
+				gcransac::estimator::LinearModelEstimator<_MinimalSolverEngine, _NonMinimalSolverEngine, 3>()
+			{}
+
+			inline double residualForScoring(const cv::Mat& point_,
+                const gcransac::Model& model_) const
+			{
+				return residual(point_, model_.descriptor);
+			}
+
+			static constexpr double getSigmaQuantile()
+			{
+				return 2.5758293035489008;
+			}
+
+			static constexpr size_t getDegreesOfFreedom()
+			{
+				return 1;
+			}
+
+			static constexpr double getC()
+			{
+				return 0.25;
+			}
+
+			static constexpr double getGammaFunction()
+			{
+				return 1.0; // Not used for DoF=1 (Γ(0) diverges); kept for interface compatibility
+			}
+
+			static constexpr bool doesNormalizationForNonMinimalFitting()
+			{
+				return false;
+			}
+
+			static constexpr double getUpperIncompleteGammaOfK()
+			{
+				return 0.0087462271161219563;
+			}
+
+			static constexpr double getLowerIncompleteGammaOfK()
+			{
+				return 0.96375479928483043;
 			}
 
 			static constexpr double getChiSquareParamCp()
 			{
 				return 1.0 / (4.0 * getGammaFunction());
 			}
 		};
+
+		// MAGSAC wrapper for PnP estimation
+		template<class _MinimalSolverEngine, class _NonMinimalSolverEngine>
+			class PerspectiveNPointEstimator :
+			public gcransac::estimator::PerspectiveNPointEstimator<_MinimalSolverEngine, _NonMinimalSolverEngine>
+		{
+		public:
+			using gcransac::estimator::PerspectiveNPointEstimator<_MinimalSolverEngine, _NonMinimalSolverEngine>::residual;
+
+			PerspectiveNPointEstimator() :
+				gcransac::estimator::PerspectiveNPointEstimator<_MinimalSolverEngine, _NonMinimalSolverEngine>()
+			{}
+
+			inline double residualForScoring(const cv::Mat& point_,
+                const gcransac::Model& model_) const
+			{ return residual(point_, model_.descriptor); }
+
+			static constexpr double getSigmaQuantile() { return 3.0348542587702925; }
+			static constexpr size_t getDegreesOfFreedom() { return 2; }
+			static constexpr double getC() { return 0.14104739588693907; }
+			static constexpr double getGammaFunction() { return 1.7724538509055159; }
+			static constexpr bool doesNormalizationForNonMinimalFitting() { return false; }
+			static constexpr double getUpperIncompleteGammaOfK() { return 0.0042654446820694645; }
+			static constexpr double getLowerIncompleteGammaOfK() { return 0.86263454284882968; }
+			static constexpr double getChiSquareParamCp() { return 1.0 / (4.0 * getGammaFunction()); }
+		};
+
+		// MAGSAC wrapper for radial homography estimation
+		template<class _MinimalSolverEngine, class _NonMinimalSolverEngine>
+			class RadialHomographyEstimator :
+			public gcransac::estimator::RadialHomographyEstimator<_MinimalSolverEngine, _NonMinimalSolverEngine>
+		{
+		public:
+			using gcransac::estimator::RadialHomographyEstimator<_MinimalSolverEngine, _NonMinimalSolverEngine>::residual;
+
+			RadialHomographyEstimator() :
+				gcransac::estimator::RadialHomographyEstimator<_MinimalSolverEngine, _NonMinimalSolverEngine>()
+			{}
+
+			inline double residualForScoring(const cv::Mat& point_,
+                const gcransac::Model& model_) const
+			{ return residual(point_, model_.descriptor); }
+
+			static constexpr double getSigmaQuantile() { return 3.0348542587702925; }
+			static constexpr size_t getDegreesOfFreedom() { return 2; }
+			static constexpr double getC() { return 0.14104739588693907; }
+			static constexpr double getGammaFunction() { return 1.7724538509055159; }
+			static constexpr bool doesNormalizationForNonMinimalFitting() { return true; }
+			static constexpr double getUpperIncompleteGammaOfK() { return 0.0042654446820694645; }
+			static constexpr double getLowerIncompleteGammaOfK() { return 0.86263454284882968; }
+			static constexpr double getChiSquareParamCp() { return 1.0 / (4.0 * getGammaFunction()); }
+		};
+
+		// MAGSAC wrapper for homography from affine correspondences
+		template<class _MinimalSolverEngine, class _NonMinimalSolverEngine>
+			class HomographyAffineEstimator :
+			public gcransac::estimator::RobustHomographyAffineCorrespondenceEstimator<_MinimalSolverEngine, _NonMinimalSolverEngine>
+		{
+		public:
+			using gcransac::estimator::RobustHomographyAffineCorrespondenceEstimator<_MinimalSolverEngine, _NonMinimalSolverEngine>::residual;
+
+			HomographyAffineEstimator() :
+				gcransac::estimator::RobustHomographyAffineCorrespondenceEstimator<_MinimalSolverEngine, _NonMinimalSolverEngine>()
+			{}
+
+			inline double residualForScoring(const cv::Mat& point_,
+                const gcransac::Model& model_) const
+			{ return residual(point_, model_.descriptor); }
+
+			static constexpr double getSigmaQuantile() { return 3.0348542587702925; }
+			static constexpr size_t getDegreesOfFreedom() { return 2; }
+			static constexpr double getC() { return 0.14104739588693907; }
+			static constexpr double getGammaFunction() { return 1.7724538509055159; }
+			static constexpr bool doesNormalizationForNonMinimalFitting() { return true; }
+			static constexpr double getUpperIncompleteGammaOfK() { return 0.0042654446820694645; }
+			static constexpr double getLowerIncompleteGammaOfK() { return 0.86263454284882968; }
+			static constexpr double getChiSquareParamCp() { return 1.0 / (4.0 * getGammaFunction()); }
+		};
 	}
 
 	namespace utils
@@ -681,5 +836,55 @@ namespace magsac
 		typedef estimator::Line2DEstimator<gcransac::estimator::solver::LinearModelSolver<2>, // The solver used for fitting a model to a minimal sample
 			gcransac::estimator::solver::LinearModelSolver<2>>  // The solver used for fitting a model to a non-minimal sample
 			Default2DLineEstimator;
+
+		// The default estimator for 3D plane fitting
+		typedef estimator::Plane3DEstimator<gcransac::estimator::solver::LinearModelSolver<3>, // The solver used for fitting a model to a minimal sample
+			gcransac::estimator::solver::LinearModelSolver<3>>  // The solver used for fitting a model to a non-minimal sample
+			Default3DPlaneEstimator;
+
+		// The default estimator for PnP (P3P minimal, PnP BA non-minimal) — matches upstream
+		typedef estimator::PerspectiveNPointEstimator<gcransac::estimator::solver::P3PSolver,
+			gcransac::estimator::solver::PnPBundleAdjustment>
+			DefaultPnPEstimator;
+
+		// PnP from 1 affine correspondence (ACP1P)
+		typedef estimator::PerspectiveNPointEstimator<gcransac::estimator::solver::P1ACQuerySolver,
+			gcransac::estimator::solver::PnPBundleAdjustment>
+			DefaultACP1PEstimator;
+
+		// PnP from 1 SIFT correspondence (SIFTP1P)
+		typedef estimator::PerspectiveNPointEstimator<gcransac::estimator::solver::SIFTP1PSolver,
+			gcransac::estimator::solver::PnPBundleAdjustment>
+			DefaultSIFTP1PEstimator;
+
+		// PnP from 2 SIFT correspondences (SIFTP2P)
+		typedef estimator::PerspectiveNPointEstimator<gcransac::estimator::solver::SIFTP2PQuerySolver,
+			gcransac::estimator::solver::PnPBundleAdjustment>
+			DefaultSIFTP2PEstimator;
+
+		// The default estimator for radial homography — matches upstream
+		typedef estimator::RadialHomographyEstimator<gcransac::estimator::solver::RadialHomography5PC,
+			gcransac::estimator::solver::RadialHomography6PC>
+			DefaultRadialHomographyEstimator;
+
+		// The default estimator for homography from affine correspondences — matches upstream
+		typedef estimator::HomographyAffineEstimator<gcransac::estimator::solver::HomographyTwoAffineSolver,
+			gcransac::estimator::solver::HomographyFourPointSolver>
+			DefaultHomographyAffineEstimator;
+
+		// The default estimator for fundamental matrix from SIFT correspondences — matches upstream
+		typedef estimator::FundamentalMatrixEstimator<gcransac::estimator::solver::FundamentalMatrixFourSIFTSolver,
+			gcransac::estimator::solver::FundamentalMatrixEightPointSolver>
+			DefaultFundamentalMatrixAffineSIFTEstimator;
+
+		// The default estimator for essential matrix (planar motion, 2 points) — matches upstream
+		typedef estimator::EssentialMatrixEstimator<gcransac::estimator::solver::EssentialMatrixTwoPointsPlanar,
+			gcransac::estimator::solver::EssentialMatrixBundleAdjustmentSolver>
+			DefaultPlanarEssentialMatrixEstimator;
+
+		// The default estimator for essential matrix with gravity (3 points) — matches upstream
+		typedef estimator::EssentialMatrixEstimator<gcransac::estimator::solver::EssentialMatrixThreePointsGravity,
+			gcransac::estimator::solver::EssentialMatrixBundleAdjustmentSolver>
+			DefaultGravityEssentialMatrixEstimator;
 	}
 }