method null revised and moved to new file kernel.h

jmeidow · jmeidow · commit 4e40205bc2bc · 2026-03-11T08:59:11.000+01:00
diff --git a/src/adjustment.cpp b/src/adjustment.cpp
@@ -21,6 +21,7 @@
 #include "constraints.h"
 #include "geometry/minrot.h"
 #include "global.h"
+#include "kernel.h"
 #include "matfun.h"
 #include "matrix.h"
 
@@ -54,6 +55,7 @@ using Eigen::Vector3d;
 using Eigen::MatrixXd;
 using Eigen::Matrix3d;
 using Eigen::SparseMatrix;
+using Eigen::Matrix;
 
 using Geometry::Rot_ab;
 
@@ -92,7 +94,8 @@ void AdjustmentFramework::update( const VectorXd &x)
         // m.segment(idx3,3).normalize();
 
         // (2) via retraction ......................................................
-        const Vector3d v = null( l0_.segment(idx3, 3) ) * x.segment(2 * s, 2);
+        const Vector3d l0 = l0_.segment(idx3,3);
+        const Vector3d v = null( l0 ) * x.segment(2 * s, 2);
         const double nv = v.norm();
         if ( nv > T_zero ) {
             const Vector3d p = l0_.segment(idx3, 3);
@@ -235,7 +238,8 @@ void AdjustmentFramework::reduce (
         const Index offset3 = 3 * s;
         const Index offset2 = 2 * s;
 
-        const Eigen::Matrix<double, 3, 2> NN = null( l0_.segment(offset3, 3) );
+        const Vector3d l0 = l0_.segment(offset3,3);
+        const Matrix<double, 3, 2> NN = null( l0 );
 
         // (i) reduced coordinates of observations
         lr.segment(offset2,2) = NN.adjoint() * l_.segment(offset3,3);
diff --git a/src/constraints.cpp b/src/constraints.cpp
@@ -26,8 +26,10 @@
 
 #include "geometry/minrot.h"
 #include "geometry/skew.h"
+#include "kernel.h"
 #include "matfun.h"
 
+
 using Geometry::skew;
 using Geometry::Rot_ab;
 
diff --git a/src/greasepad.pro b/src/greasepad.pro
@@ -58,6 +58,7 @@ HEADERS += \
     geometry/minrot.h \
     geometry/skew.h \
     global.h \
+    kernel.h \
     mainwindow.h \
     mainscene.h \
     mainview.h \
diff --git a/src/kernel.h b/src/kernel.h
@@ -0,0 +1,62 @@
+/*
+ * This file is part of the GreasePad distribution (https://github.com/FraunhoferIOSB/GreasePad).
+ * Copyright (c) 2022-2026 Jochen Meidow, Fraunhofer IOSB
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program. If not, see <https://www.gnu.org/licenses/>.
+ */
+
+
+#ifndef KERNEL_H
+#define KERNEL_H
+
+#include <Eigen/Core>
+
+#include <cassert>
+
+namespace Matfun {
+
+using Eigen::Matrix;
+using Eigen::Vector;
+
+
+//! Basis vectors for the nullspace of a row vector
+//!
+//! Förstner & Wrobel (2016) Photogrammetric Computer Vision, eq. (A.120)
+template<int N>
+[[nodiscard]] static Matrix<double,N,N-1>
+null( const Vector<double,N> & xs )
+{
+    constexpr double T_zero = 1e-6;
+    assert( std::fabs( xs.norm()-1. ) < T_zero && "vector norm not 1");
+
+    Vector<double,N-1> x0 = xs.head(N-1);
+    double xN = xs(N-1);
+    if ( xN < 0 ) {
+        x0 = -x0;
+        xN = -xN;
+    }
+
+    const Matrix<double,N,N-1> JJ = (Matrix<double,N,N-1>() <<
+        Matrix<double,N-1,N-1>::Identity(N-1,N-1) -x0*x0.transpose()/(1.+xN),
+        -x0.transpose() ).finished();
+
+    const Vector<double,N-1> check = JJ.adjoint()*xs;
+    assert( check.norm() < T_zero && "not a zero vector");
+
+    return JJ;
+}
+
+} // namespace Matfun
+
+#endif // KERNEL_H
diff --git a/src/matfun.h b/src/matfun.h
@@ -27,7 +27,6 @@
 
 #include <algorithm>
 #include <cassert>
-#include <cmath>
 
 
 namespace Matfun {
@@ -39,35 +38,9 @@ using Eigen::SparseMatrix;
 using Eigen::SparseVector;
 using Eigen::Matrix3d;
 using Eigen::Index;
+using Eigen::Matrix;
 using Eigen::Vector;
 using Eigen::Dynamic;
-using Eigen::placeholders::last;
-
-
-//! Nullspace of row vector
-[[nodiscard,maybe_unused]] static MatrixXd null( const VectorXd & xs )
-{
-    // cf. PCV, eq. (A.120)
-    constexpr double T_zero = 1e-6;
-    assert( std::fabs( xs.norm()-1. ) < T_zero && "vector norm not 1");
-
-    const Index N = xs.size();
-    VectorXd x0 = xs.head(N-1);
-    double   xN = xs(last);
-    if ( xN < 0.0 ) {
-        x0 = -x0;
-        xN = -xN;
-    }
-
-    const MatrixXd JJ = ( MatrixXd(N,N-1)
-         << MatrixXd::Identity(N-1,N-1) -x0*x0.adjoint()/(1.+xN),
-            -x0.adjoint() ).finished();
-
-    const VectorXd check = JJ.adjoint()*xs;
-    assert( check.norm() < T_zero && "not a zero vector");
-
-    return JJ;
-}
 
 
 //! check if the matrix AA is rank-deficient
diff --git a/src/uncertain/uncertain.cpp b/src/uncertain/uncertain.cpp
@@ -16,7 +16,7 @@
  * along with this program. If not, see <https://www.gnu.org/licenses/>.
  */
 
-
+#include "kernel.h"
 #include "matfun.h"
 #include "uncertain.h"
 
