ICPR 2026 solvers

HomLib/CMakeLists.txt

@@ -1,10 +1,12 @@
 set(SRCS
+    helpers/affine2sift.cpp
     helpers/generate_problem_instance.cpp
     helpers/gj.cpp
     helpers/normalize2dpts.cpp
     helpers/radial.cpp
     helpers/roots.cpp
     robust/refinement.cpp
+    solvers/barath_visapp_2016/get_barath_visapp_2016_affine.cpp
     solvers/fitzgibbon_cvpr_2001/get_fitzgibbon_cvpr_2001.cpp
     solvers/fitzgibbon_cvpr_2001/get_fitzgibbon_cvpr_2001_single.cpp
     solvers/kukelova_cvpr_2015/get_kukelova_cvpr_2015.cpp
@@ -13,6 +15,8 @@ set(SRCS
     solvers/nakano_icpr_2025/get_nakano_icpr_2025.cpp
     solvers/valtonenornhag_icpr_2020/get_valtonenornhag_icpr_2020_fHf.cpp
     solvers/valtonenornhag_icpr_2020/solver_valtonenornhag_icpr_2020_fHf.cpp
+    solvers/valtonenornhag_icpr_2026/get_valtonenornhag_icpr_2026_affine.cpp
+    solvers/valtonenornhag_icpr_2026/get_valtonenornhag_icpr_2026_ori.cpp
     solvers/valtonenornhag_wacv_2021/get_valtonenornhag_wacv_2021_fHf.cpp
     solvers/valtonenornhag_wacv_2021/get_valtonenornhag_wacv_2021_frHfr.cpp
     solvers/valtonenornhag_wacv_2021/solver_valtonenornhag_wacv_2021_fHf.cpp

HomLib/helpers/affine2sift.cpp

@@ -0,0 +1,154 @@
+#include <vector>
+#include <numeric>
+#include <Eigen/Dense>
+
+
+
+
+
+using namespace Eigen;
+
+
+
+void affine2sift_solver(Eigen::MatrixXd const& A, double q1, std::vector<double>* r1, std::vector<double>* r2)
+{
+    // Compute coefficients
+double const* A_data = A.data();
+double const A11 = A_data[0];
+double const A21 = A_data[1];
+double const A12 = A_data[2];
+double const A22 = A_data[3];
+    VectorXd coeffs(14);
+  double _t2_ = A12*2.0;
+  double _t3_ = A22*2.0;
+  double _t4_ = q1*2.0;
+  double _t5_ = -A11;
+  double _t6_ = -A21;
+  double _t7_ = -q1;
+  double _t8_ = -_t4_;
+  coeffs[0] = q1+_t5_;
+  coeffs[1] = _t5_+_t7_;
+  coeffs[2] = _t2_;
+  coeffs[3] = _t2_;
+  coeffs[4] = A11+q1;
+  coeffs[5] = A11+_t7_;
+  coeffs[6] = _t6_;
+  coeffs[7] = _t8_;
+  coeffs[8] = _t6_;
+  coeffs[9] = _t3_;
+  coeffs[10] = _t3_;
+  coeffs[11] = A21;
+  coeffs[12] = _t8_;
+  coeffs[13] = A21;
+
+
+    // Setup elimination template
+    static const int coeffs0_ind[] = { 0,6,0,6,7,0,2,6,9,1,7,8,1,8,1,3,7,8,10,0,2,6,7,9,2,4,9,11 };
+    static const int coeffs1_ind[] = { 5,13,3,5,10,13,1,3,8,10,3,5,10,12,13,2,4,9,11,12,5,12,13,4,11,12,4,11 };
+
+
+static const int C0_ind[] = {3,7,10,14,15,17,19,21,23,27,30,31,34,38,41,43,44,45,47,48,50,52,53,54,57,59,61,63};
+
+static const int C1_ind[] = {0,4,8,10,12,14,16,18,20,22,25,27,29,30,31,32,34,36,38,39,41,44,45,48,52,53,57,61};
+
+MatrixXd C0 = MatrixXd::Zero(8,8);
+MatrixXd C1 = MatrixXd::Zero(8,8);
+for (int i = 0; i < 28; i++) {
+    C0(C0_ind[i]) = coeffs(coeffs0_ind[i]);
+}
+
+for (int i = 0; i < 28; i++) {
+    C1(C1_ind[i]) = coeffs(coeffs1_ind[i]);
+}
+
+/*
+Eigen::Matrix<double,8,4> b = Eigen::Matrix<double,8,4>::Zero();
+b(4,0) = -1;
+b(5,1) = -1;
+b(6,2) = -1;
+b(7,3) = -1;
+Eigen::Matrix<double,8,4> alpha = C0.transpose().fullPivLu().solve(b);
+Eigen::Matrix<double,12,4> RR;
+RR << alpha.transpose()*C1, Eigen::Matrix<double,8,8>::Identity();
+//AM_ind = [6,7,1,2,3,8,9,4];
+//AM = RR(AM_ind,:);
+Eigen::Matrix<double,8,8> AM;
+AM << RR.col(5), RR.col(6), RR.col(0), RR.col(1), RR.col(2), RR.col(7), RR.col(8), RR.col(3);
+Eigen::EigenSolver< Eigen::Matrix<double,8,8> > AMsolver(AM);
+Eigen::MatrixXcd V = AMsolver.eigenvectors();
+V = V.array() * (Eigen::Matrix<double,8,1>::Ones()*V.col(0)).array();
+Eigen::VectorXcd r1c = AMsolver.eigenvalues();
+Eigen::VectorXcd r2c = V.row(5);
+*/
+
+//[V,D] = eig(AM);
+//V = V ./ (ones(size(V,1),1)*V(1,:));
+//sols(1,:) = diag(D).';
+//sols(2,:) = V(6,:);
+
+MatrixXd C12 = C0.fullPivLu().solve(C1);
+
+
+
+    // Setup action matrix
+    Matrix<double,12, 8> RR;
+    RR << -C12.bottomRows(4), Matrix<double,8,8>::Identity(8, 8);
+
+    static const int AM_ind[] = { 5,6,0,1,2,7,8,3 };
+    Matrix<double, 8, 8> AM;
+    for (int i = 0; i < 8; i++) {
+        AM.row(i) = RR.row(AM_ind[i]);
+    }
+
+    MatrixXcd sols(2, 8);
+    sols.setZero();
+
+    // Solve eigenvalue problem
+    EigenSolver<Matrix<double, 8, 8> > es(AM);
+    ArrayXcd D = es.eigenvalues();
+    ArrayXXcd V = es.eigenvectors();
+
+    V = (V / V.row(0).array().replicate(8, 1)).eval();
+
+
+        sols.row(0) = D.transpose().array();
+    sols.row(1) = V.row(5).array();
+
+
+
+
+    Eigen::VectorXcd r1c = sols.row(0);
+    Eigen::VectorXcd r2c = sols.row(1);
+    int nsols = r1c.size();
+    for (int isol = 0; isol < nsols; ++isol) {
+        if ( r1c(isol).imag() == 0 && r2c(isol).imag() == 0 )
+        {
+            r1->push_back(r1c(isol).real());
+            r2->push_back(r2c(isol).real());
+        }
+    }
+}
+
+// Action =
+// Quotient ring basis (V) =  r1^2, r1*r2, r1*r2^2, r2, r2^2, r2^3
+// Available monomials (RR*V) = r1^2*r2, r1^2*r2^2, r1*r2^3, 1, r1, r1^2, r1*r2, r1*r2^2, r2, r2^2, r2^3
+
+void affine2sift(const Eigen::Matrix2d &A, double &s1, double &c1, double &s2, double &c2, double &q )
+{
+    q = sqrt(A.determinant());
+    std::vector<double> r1solns, r2solns;
+    affine2sift_solver(A, q, &r1solns, &r2solns);
+
+    double r1 = r1solns[0];
+    double r2 = r2solns[0];
+    c1 = (1-r1*r1)/(1+r1*r1);
+    s1 = (2*r1)/(1+r1*r1);
+    c2 = (1-r2*r2)/(1+r2*r2);
+    s2 = (2*r2)/(1+r2*r2);
+
+    // check residuals
+    //double res1 = c1*s2*A(0,0) + s1*s2*A(0,1) - c1*c2*A(1,0) - c2*s1*A(1,1);
+    //double res2 = A(0,1)*A(1,0)-A(0,0)*A(1,1)+q*q;
+    //double res3 = A(0,0)*c1 + A(0,1)*s1 - c2*q;
+    //double res4 = A(1,0)*c1 + A(1,1)*s1 - s2*q;
+}

HomLib/helpers/affine2sift.hpp

@@ -0,0 +1,18 @@
+#ifndef SRC_HELPERS_AFFINE2SIFT_HPP_
+#define SRC_HELPERS_AFFINE2SIFT_HPP_
+
+#include <vector>
+#include <numeric>
+#include <Eigen/Dense>
+
+using namespace Eigen;
+
+// Converts an affine correspondence to a scale-and-orientation correspondence.
+// There are eight possible solutions, but this arbitrarily returns the first one.
+void affine2sift(const Eigen::Matrix2d &A, // input affine transformation matrix
+                 double &s_ref, double &c_ref, // sine and cosine of feature orientation in reference image
+                 double &s_query, double &c_query, // sine and cosine of feature orientation in query image
+                 double &q // ratio of feature scales (scale in query image / scale in reference image)
+);
+
+#endif  // SRC_HELPERS_AFFINE2SIFT_HPP_

HomLib/helpers/generate_problem_instance.cpp

@@ -21,15 +21,40 @@
 
 #include <Eigen/Dense>
 #include <random>
+#include <cmath>
 #include <chrono>
 
 #include "problem_instance.hpp"
 #include "radial.hpp"
 #include "generate_problem_instance.hpp"
+#include "affine2sift.hpp"
 
 namespace HomLib {
 
     static const double kPI = 3.14159265358979323846;
+
+    static Eigen::Matrix2d affineFromHomography( const Eigen::Matrix3d &H, const Eigen::Vector2d &x, const Eigen::Vector2d &y , double k)
+    {
+        // x are the distorted coeffs
+        double h1 = H(0,0), h2 = H(0,1), h3 = H(0,2),
+               h4 = H(1,0), h5 = H(1,1), h6 = H(1,2),
+               h7 = H(2,0), h8 = H(2,1), h9 = H(2,2);
+        double u1 = x(0), v1 = x(1),
+               u2 = y(0), v2 = y(1);
+        double dist_fact = k*(u1*u1 + v1*v1) + 1;
+        double s = h7*u1 + h8*v1 + h9*(dist_fact);
+         // Note that
+        // u2 = (h1*u1 + h2*v1 + h3*(dist_fact))/(s)
+        // v2 = (h4*u1 + h5*v1 + h6*(dist_fact))/(s)
+        Eigen::Matrix2d A;
+        A << (h1 + 2*h3*k*u1)/(s) - ((h7 + 2*h9*k*u1)*(h1*u1 + h2*v1 + h3*(dist_fact)))/(s*s),
+             (h2 + 2*h3*k*v1)/(s) - ((h8 + 2*h9*k*v1)*(h1*u1 + h2*v1 + h3*(dist_fact)))/(s*s),
+             (h4 + 2*h6*k*u1)/(s) - ((h7 + 2*h9*k*u1)*(h4*u1 + h5*v1 + h6*(dist_fact)))/(s*s),
+             (h5 + 2*h6*k*v1)/(s) - ((h8 + 2*h9*k*v1)*(h4*u1 + h5*v1 + h6*(dist_fact)))/(s*s);
+
+
+        return A;
+    }
 
     HomLib::ProblemInstance generate_problem_instance(const ProblemConfig &config) {
 
@@ -42,7 +67,7 @@ namespace HomLib {
 		random_engine.seed(std::chrono::system_clock::now().time_since_epoch().count());
 		std::uniform_real_distribution<double> depth_gen(config.min_depth_, config.max_depth_);
 		std::uniform_real_distribution<double> coord_gen(-fov_scale, fov_scale);
-		// std::uniform_real_distribution<double> focal_gen(config.min_focal_, config.max_focal_);
+		std::uniform_real_distribution<double> focal_gen(config.min_focal_, config.max_focal_);
 		std::normal_distribution<double> direction_gen(0.0, 1.0);
 		std::uniform_real_distribution<double> dist_gen(config.min_dist_, config.max_dist_);
 
@@ -54,13 +79,17 @@ namespace HomLib {
 		    t.normalize();
 		    Eigen::Matrix3d R = Eigen::Quaternion<double>::UnitRandom().toRotationMatrix();
 
-		    // double focal_gt = focal_gen(random_engine);
+		    double focal_gt = focal_gen(random_engine);
 
 	        // Point to point correspondences
 	        instance.x1.clear();
 	        instance.x2.clear();
+	        instance.A.clear();
+	        instance.ori.clear();
 	        instance.x1.reserve(config.number_points);
 	        instance.x2.reserve(config.number_points);
+	        instance.A.reserve(config.number_points);
+	        instance.ori.reserve(config.number_points);
 
 		    // Generate plane
 		    Eigen::Vector3d n;
@@ -74,6 +103,30 @@ namespace HomLib {
 
 		    // ground truth homography
 		    instance.posedata.homography = alpha * R + t * n.transpose();
+
+		    // Distort
+		    switch (config.distortion) {
+        	    case HomLib::DistortionCase::NO_DISTORTION:
+                    instance.posedata.distortion_parameter = 0.0;
+		            instance.posedata.distortion_parameter2 = 0.0;
+		            break;
+        	    case HomLib::DistortionCase::ONE_SIDED_LEFT:
+                    instance.posedata.distortion_parameter = 0.0;
+		            instance.posedata.distortion_parameter2 = dist_gen(random_engine);
+		            break;
+        	    case HomLib::DistortionCase::ONE_SIDED_RIGHT:
+                    instance.posedata.distortion_parameter = dist_gen(random_engine);
+		            instance.posedata.distortion_parameter2 = 0.0;
+		            break;
+        	    case HomLib::DistortionCase::TWO_SIDED_EQUAL:
+                    instance.posedata.distortion_parameter = dist_gen(random_engine);
+		            instance.posedata.distortion_parameter2 = instance.posedata.distortion_parameter;
+		            break;
+        	    case HomLib::DistortionCase::TWO_SIDED: 
+        	        instance.posedata.distortion_parameter = dist_gen(random_engine);
+		            instance.posedata.distortion_parameter2 = dist_gen(random_engine);
+		            break;
+    	    }
 
 		    bool failed_instance = false;
 		    for (int j = 0; j < config.number_points; ++j) {
@@ -103,11 +156,51 @@ namespace HomLib {
 		                // try to generate another point
 		                continue;
 		            }
+			        Eigen::Vector2d x1h = x1.hnormalized();
+
+			        // Distort
+			        Eigen::Vector2d x1hd, x2hd;
+			        x1hd = HomLib::radialdistort(x1h, instance.posedata.distortion_parameter);
+		            x2hd = HomLib::radialdistort(x2h, instance.posedata.distortion_parameter2);
+		            // calculate affine from homography
+		            // This assumes DistortionCase.NO_DISTORTION or DistortionCase.ONE_SIDED_RIGHT
+		            Eigen::Matrix2d A = affineFromHomography(instance.posedata.homography, x1hd, x2hd, instance.posedata.distortion_parameter);
+
+		            // check if determinant is positive
+		            if ( A.determinant() < 0 ) continue;
 
 					//
-					Eigen::Vector2d x1h = x1.hnormalized();
-		            instance.x1.push_back(x1h);
-		            instance.x2.push_back(x2h);
+					double s_1, c_1, s_2, c_2, q;
+					affine2sift(A, s_1, c_1, s_2, c_2, q);
+                    /*  ONLY WORKS FOR RIGHT-SIDED AND NO_DISTORTION
+                    double h_1 = instance.posedata.homography(0,0),
+                        h_2 = instance.posedata.homography(0,1),
+                        h_3 = instance.posedata.homography(0,2),
+                        h_4 = instance.posedata.homography(1,0),
+                        h_5 = instance.posedata.homography(1,1),
+                        h_6 = instance.posedata.homography(1,2),
+                        h_7 = instance.posedata.homography(2,0),
+                        h_8 = instance.posedata.homography(2,1),
+                        h_9 = instance.posedata.homography(2,2);
+                    double u_1 = x1hd[0],
+                        v_1 = x1hd[1],
+                        u_2 = x2hd[0],
+                        v_2 = x2hd[1];
+                    lambda = instance.posedata.distortion_parameter;
+
+					double res = -h_1*s_2*c_1 - h_2*s_1*s_2 + h_4*c_1*c_2 + h_5*s_1*c_2 + h_7*u_2*s_2*c_1 - h_7*v_2*c_1*c_2 + h_8*u_2*s_1*s_2
+                        -h_8*v_2*s_1*c_2 -2*h_3*u_1*s_2*c_1*lambda - 2*h_3*v_1*s_1*s_2*lambda + 2*h_6*u_1*c_1*c_2*lambda + 2*h_6*v_1*s_1*c_2*lambda
+                        + 2*h_9*u_1*u_2*s_2*c_1*lambda - 2*h_9*u_1*v_2*c_1*c_2*lambda + 2*h_9*v_1*u_2*s_1*s_2*lambda - 2*h_9*v_1*v_2*s_1*c_2*lambda;
+                    */
+
+					Eigen::Vector2d ori;
+					ori[0] = std::atan2(s_1, c_1);
+					ori[1] = std::atan2(s_2, c_2);
+
+		            instance.x1.push_back(x1hd);
+		            instance.x2.push_back(x2hd);
+		            instance.A.push_back(A);
+		            instance.ori.push_back(ori);
 
 		            point_okay = true;
 		            break;
@@ -120,28 +213,6 @@ namespace HomLib {
 		    if (failed_instance) {
 		        continue;
 		    }
-
-		    // Distort
-		    if (config.no_distortion) {
-		        instance.posedata.distortion_parameter = 0.0;
-		        instance.posedata.distortion_parameter2 = 0.0;
-		    } else {
-		        instance.posedata.distortion_parameter2 = dist_gen(random_engine);
-		        if (config.one_sided) {
-			    instance.posedata.distortion_parameter = 0.0;
-		        } else {
-			    if (config.equal) {
-			        instance.posedata.distortion_parameter = instance.posedata.distortion_parameter2;
-			    } else {
-			        instance.posedata.distortion_parameter = dist_gen(random_engine);
-			    }
-		        }
-
-		        if (!config.one_sided) { 
-			    HomLib::radialdistort(instance.x1, &instance.x1, instance.posedata.distortion_parameter);
-		        }
-		        HomLib::radialdistort(instance.x2, &instance.x2, instance.posedata.distortion_parameter2);
-		    }
 
 		    // Focal length		
 			//instance.x1 *= focal_gt;