|  | // Ceres Solver - A fast non-linear least squares minimizer | 
|  | // Copyright 2022 Google Inc. All rights reserved. | 
|  | // http://ceres-solver.org/ | 
|  | // | 
|  | // Redistribution and use in source and binary forms, with or without | 
|  | // modification, are permitted provided that the following conditions are met: | 
|  | // | 
|  | // * Redistributions of source code must retain the above copyright notice, | 
|  | //   this list of conditions and the following disclaimer. | 
|  | // * Redistributions in binary form must reproduce the above copyright notice, | 
|  | //   this list of conditions and the following disclaimer in the documentation | 
|  | //   and/or other materials provided with the distribution. | 
|  | // * Neither the name of Google Inc. nor the names of its contributors may be | 
|  | //   used to endorse or promote products derived from this software without | 
|  | //   specific prior written permission. | 
|  | // | 
|  | // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" | 
|  | // AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE | 
|  | // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE | 
|  | // ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE | 
|  | // LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR | 
|  | // CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF | 
|  | // SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS | 
|  | // INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN | 
|  | // CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) | 
|  | // ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE | 
|  | // POSSIBILITY OF SUCH DAMAGE. | 
|  | // | 
|  | // Author: sameeragarwal@google.com (Sameer Agarwal) | 
|  |  | 
|  | #ifndef CERES_PUBLIC_MANIFOLD_H_ | 
|  | #define CERES_PUBLIC_MANIFOLD_H_ | 
|  |  | 
|  | #include <Eigen/Core> | 
|  | #include <algorithm> | 
|  | #include <array> | 
|  | #include <memory> | 
|  | #include <utility> | 
|  | #include <vector> | 
|  |  | 
|  | #include "ceres/internal/disable_warnings.h" | 
|  | #include "ceres/internal/export.h" | 
|  | #include "ceres/types.h" | 
|  | #include "glog/logging.h" | 
|  |  | 
|  | namespace ceres { | 
|  |  | 
|  | // In sensor fusion problems, often we have to model quantities that live in | 
|  | // spaces known as Manifolds, for example the rotation/orientation of a sensor | 
|  | // that is represented by a quaternion. | 
|  | // | 
|  | // Manifolds are spaces which locally look like Euclidean spaces. More | 
|  | // precisely, at each point on the manifold there is a linear space that is | 
|  | // tangent to the manifold. It has dimension equal to the intrinsic dimension of | 
|  | // the manifold itself, which is less than or equal to the ambient space in | 
|  | // which the manifold is embedded. | 
|  | // | 
|  | // For example, the tangent space to a point on a sphere in three dimensions is | 
|  | // the two dimensional plane that is tangent to the sphere at that point. There | 
|  | // are two reasons tangent spaces are interesting: | 
|  | // | 
|  | // 1. They are Eucliean spaces so the usual vector space operations apply there, | 
|  | //    which makes numerical operations easy. | 
|  | // 2. Movement in the tangent space translate into movements along the manifold. | 
|  | //    Movements perpendicular to the tangent space do not translate into | 
|  | //    movements on the manifold. | 
|  | // | 
|  | // Returning to our sphere example, moving in the 2 dimensional plane | 
|  | // tangent to the sphere and projecting back onto the sphere will move you away | 
|  | // from the point you started from but moving along the normal at the same point | 
|  | // and the projecting back onto the sphere brings you back to the point. | 
|  | // | 
|  | // The Manifold interface defines two operations (and their derivatives) | 
|  | // involving the tangent space, allowing filtering and optimization to be | 
|  | // performed on said manifold: | 
|  | // | 
|  | // 1. x_plus_delta = Plus(x, delta) | 
|  | // 2. delta = Minus(x_plus_delta, x) | 
|  | // | 
|  | // "Plus" computes the result of moving along delta in the tangent space at x, | 
|  | // and then projecting back onto the manifold that x belongs to. In Differential | 
|  | // Geometry this is known as a "Retraction". It is a generalization of vector | 
|  | // addition in Euclidean spaces. | 
|  | // | 
|  | // Given two points on the manifold, "Minus" computes the change delta to x in | 
|  | // the tangent space at x, that will take it to x_plus_delta. | 
|  | // | 
|  | // Let us now consider two examples. | 
|  | // | 
|  | // The Euclidean space R^n is the simplest example of a manifold. It has | 
|  | // dimension n (and so does its tangent space) and Plus and Minus are the | 
|  | // familiar vector sum and difference operations. | 
|  | // | 
|  | //  Plus(x, delta) = x + delta = y, | 
|  | //  Minus(y, x) = y - x = delta. | 
|  | // | 
|  | // A more interesting case is SO(3), the special orthogonal group in three | 
|  | // dimensions - the space of 3x3 rotation matrices. SO(3) is a three dimensional | 
|  | // manifold embedded in R^9 or R^(3x3). So points on SO(3) are represented using | 
|  | // 9 dimensional vectors or 3x3 matrices, and points in its tangent spaces are | 
|  | // represented by 3 dimensional vectors. | 
|  | // | 
|  | // Defining Plus and Minus are defined in terms of the matrix Exp and Log | 
|  | // operations as follows: | 
|  | // | 
|  | // Let Exp(p, q, r) = [cos(theta) + cp^2, -sr + cpq        ,  sq + cpr        ] | 
|  | //                    [sr + cpq         , cos(theta) + cq^2, -sp + cqr        ] | 
|  | //                    [-sq + cpr        , sp + cqr         , cos(theta) + cr^2] | 
|  | // | 
|  | // where: theta = sqrt(p^2 + q^2 + r^2) | 
|  | //            s = sinc(theta) | 
|  | //            c = (1 - cos(theta))/theta^2 | 
|  | // | 
|  | // and Log(x) = 1/(2 sinc(theta))[x_32 - x_23, x_13 - x_31, x_21 - x_12] | 
|  | // | 
|  | // where: theta = acos((Trace(x) - 1)/2) | 
|  | // | 
|  | // Then, | 
|  | // | 
|  | // Plus(x, delta) = x Exp(delta) | 
|  | // Minus(y, x) = Log(x^T y) | 
|  | // | 
|  | // For Plus and Minus to be mathematically consistent, the following identities | 
|  | // must be satisfied at all points x on the manifold: | 
|  | // | 
|  | // 1.  Plus(x, 0) = x. | 
|  | // 2.  For all y, Plus(x, Minus(y, x)) = y. | 
|  | // 3.  For all delta, Minus(Plus(x, delta), x) = delta. | 
|  | // 4.  For all delta_1, delta_2 | 
|  | //    |Minus(Plus(x, delta_1), Plus(x, delta_2)) <= |delta_1 - delta_2| | 
|  | // | 
|  | // Briefly: | 
|  | // (1) Ensures that the tangent space is "centered" at x, and the zero vector is | 
|  | //     the identity element. | 
|  | // (2) Ensures that any y can be reached from x. | 
|  | // (3) Ensures that Plus is an injective (one-to-one) map. | 
|  | // (4) Allows us to define a metric on the manifold. | 
|  | // | 
|  | // Additionally we require that Plus and Minus be sufficiently smooth. In | 
|  | // particular they need to be differentiable everywhere on the manifold. | 
|  | // | 
|  | // For more details, please see | 
|  | // | 
|  | // "Integrating Generic Sensor Fusion Algorithms with Sound State | 
|  | // Representations through Encapsulation of Manifolds" | 
|  | // By C. Hertzberg, R. Wagner, U. Frese and L. Schroder | 
|  | // https://arxiv.org/pdf/1107.1119.pdf | 
|  | class CERES_EXPORT Manifold { | 
|  | public: | 
|  | virtual ~Manifold(); | 
|  |  | 
|  | // Dimension of the ambient space in which the manifold is embedded. | 
|  | virtual int AmbientSize() const = 0; | 
|  |  | 
|  | // Dimension of the manifold/tangent space. | 
|  | virtual int TangentSize() const = 0; | 
|  |  | 
|  | //   x_plus_delta = Plus(x, delta), | 
|  | // | 
|  | // A generalization of vector addition in Euclidean space, Plus computes the | 
|  | // result of moving along delta in the tangent space at x, and then projecting | 
|  | // back onto the manifold that x belongs to. | 
|  | // | 
|  | // x and x_plus_delta are AmbientSize() vectors. | 
|  | // delta is a TangentSize() vector. | 
|  | // | 
|  | // Return value indicates if the operation was successful or not. | 
|  | virtual bool Plus(const double* x, | 
|  | const double* delta, | 
|  | double* x_plus_delta) const = 0; | 
|  |  | 
|  | // Compute the derivative of Plus(x, delta) w.r.t delta at delta = 0, i.e. | 
|  | // | 
|  | // (D_2 Plus)(x, 0) | 
|  | // | 
|  | // jacobian is a row-major AmbientSize() x TangentSize() matrix. | 
|  | // | 
|  | // Return value indicates whether the operation was successful or not. | 
|  | virtual bool PlusJacobian(const double* x, double* jacobian) const = 0; | 
|  |  | 
|  | // tangent_matrix = ambient_matrix * (D_2 Plus)(x, 0) | 
|  | // | 
|  | // ambient_matrix is a row-major num_rows x AmbientSize() matrix. | 
|  | // tangent_matrix is a row-major num_rows x TangentSize() matrix. | 
|  | // | 
|  | // Return value indicates whether the operation was successful or not. | 
|  | // | 
|  | // This function is only used by the GradientProblemSolver, where the | 
|  | // dimension of the parameter block can be large and it may be more efficient | 
|  | // to compute this product directly rather than first evaluating the Jacobian | 
|  | // into a matrix and then doing a matrix vector product. | 
|  | // | 
|  | // Because this is not an often used function, we provide a default | 
|  | // implementation for convenience. If performance becomes an issue then the | 
|  | // user should consider implementing a specialization. | 
|  | virtual bool RightMultiplyByPlusJacobian(const double* x, | 
|  | const int num_rows, | 
|  | const double* ambient_matrix, | 
|  | double* tangent_matrix) const; | 
|  |  | 
|  | // y_minus_x = Minus(y, x) | 
|  | // | 
|  | // Given two points on the manifold, Minus computes the change to x in the | 
|  | // tangent space at x, that will take it to y. | 
|  | // | 
|  | // x and y are AmbientSize() vectors. | 
|  | // y_minus_x is a TangentSize() vector. | 
|  | // | 
|  | // Return value indicates if the operation was successful or not. | 
|  | virtual bool Minus(const double* y, | 
|  | const double* x, | 
|  | double* y_minus_x) const = 0; | 
|  |  | 
|  | // Compute the derivative of Minus(y, x) w.r.t y at y = x, i.e | 
|  | // | 
|  | //   (D_1 Minus) (x, x) | 
|  | // | 
|  | // Jacobian is a row-major TangentSize() x AmbientSize() matrix. | 
|  | // | 
|  | // Return value indicates whether the operation was successful or not. | 
|  | virtual bool MinusJacobian(const double* x, double* jacobian) const = 0; | 
|  | }; | 
|  |  | 
|  | // The Euclidean manifold is another name for the ordinary vector space R^size, | 
|  | // where the plus and minus operations are the usual vector addition and | 
|  | // subtraction: | 
|  | //   Plus(x, delta) = x + delta | 
|  | //   Minus(y, x) = y - x. | 
|  | // | 
|  | // The class works with dynamic and static ambient space dimensions. If the | 
|  | // ambient space dimensions is know at compile time use | 
|  | // | 
|  | //    EuclideanManifold<3> manifold; | 
|  | // | 
|  | // If the ambient space dimensions is not known at compile time the template | 
|  | // parameter needs to be set to ceres::DYNAMIC and the actual dimension needs | 
|  | // to be provided as a constructor argument: | 
|  | // | 
|  | //    EuclideanManifold<ceres::DYNAMIC> manifold(ambient_dim); | 
|  | template <int Size> | 
|  | class EuclideanManifold final : public Manifold { | 
|  | public: | 
|  | static_assert(Size == ceres::DYNAMIC || Size >= 0, | 
|  | "The size of the manifold needs to be non-negative."); | 
|  | static_assert(ceres::DYNAMIC == Eigen::Dynamic, | 
|  | "ceres::DYNAMIC needs to be the same as Eigen::Dynamic."); | 
|  |  | 
|  | EuclideanManifold() : size_{Size} { | 
|  | static_assert( | 
|  | Size != ceres::DYNAMIC, | 
|  | "The size is set to dynamic. Please call the constructor with a size."); | 
|  | } | 
|  |  | 
|  | explicit EuclideanManifold(int size) : size_(size) { | 
|  | if (Size != ceres::DYNAMIC) { | 
|  | CHECK_EQ(Size, size) | 
|  | << "Specified size by template parameter differs from the supplied " | 
|  | "one."; | 
|  | } else { | 
|  | CHECK_GE(size_, 0) | 
|  | << "The size of the manifold needs to be non-negative."; | 
|  | } | 
|  | } | 
|  |  | 
|  | int AmbientSize() const override { return size_; } | 
|  | int TangentSize() const override { return size_; } | 
|  |  | 
|  | bool Plus(const double* x_ptr, | 
|  | const double* delta_ptr, | 
|  | double* x_plus_delta_ptr) const override { | 
|  | Eigen::Map<const AmbientVector> x(x_ptr, size_); | 
|  | Eigen::Map<const AmbientVector> delta(delta_ptr, size_); | 
|  | Eigen::Map<AmbientVector> x_plus_delta(x_plus_delta_ptr, size_); | 
|  | x_plus_delta = x + delta; | 
|  | return true; | 
|  | } | 
|  |  | 
|  | bool PlusJacobian(const double* x_ptr, double* jacobian_ptr) const override { | 
|  | Eigen::Map<MatrixJacobian> jacobian(jacobian_ptr, size_, size_); | 
|  | jacobian.setIdentity(); | 
|  | return true; | 
|  | } | 
|  |  | 
|  | bool RightMultiplyByPlusJacobian(const double* x, | 
|  | const int num_rows, | 
|  | const double* ambient_matrix, | 
|  | double* tangent_matrix) const override { | 
|  | std::copy_n(ambient_matrix, num_rows * size_, tangent_matrix); | 
|  | return true; | 
|  | } | 
|  |  | 
|  | bool Minus(const double* y_ptr, | 
|  | const double* x_ptr, | 
|  | double* y_minus_x_ptr) const override { | 
|  | Eigen::Map<const AmbientVector> x(x_ptr, size_); | 
|  | Eigen::Map<const AmbientVector> y(y_ptr, size_); | 
|  | Eigen::Map<AmbientVector> y_minus_x(y_minus_x_ptr, size_); | 
|  | y_minus_x = y - x; | 
|  | return true; | 
|  | } | 
|  |  | 
|  | bool MinusJacobian(const double* x_ptr, double* jacobian_ptr) const override { | 
|  | Eigen::Map<MatrixJacobian> jacobian(jacobian_ptr, size_, size_); | 
|  | jacobian.setIdentity(); | 
|  | return true; | 
|  | } | 
|  |  | 
|  | private: | 
|  | static constexpr bool IsDynamic = (Size == ceres::DYNAMIC); | 
|  | using AmbientVector = Eigen::Matrix<double, Size, 1>; | 
|  | using MatrixJacobian = Eigen::Matrix<double, Size, Size, Eigen::RowMajor>; | 
|  |  | 
|  | int size_{}; | 
|  | }; | 
|  |  | 
|  | // Hold a subset of the parameters inside a parameter block constant. | 
|  | class CERES_EXPORT SubsetManifold final : public Manifold { | 
|  | public: | 
|  | SubsetManifold(int size, const std::vector<int>& constant_parameters); | 
|  | int AmbientSize() const override; | 
|  | int TangentSize() const override; | 
|  |  | 
|  | bool Plus(const double* x, | 
|  | const double* delta, | 
|  | double* x_plus_delta) const override; | 
|  | bool PlusJacobian(const double* x, double* jacobian) const override; | 
|  | bool RightMultiplyByPlusJacobian(const double* x, | 
|  | const int num_rows, | 
|  | const double* ambient_matrix, | 
|  | double* tangent_matrix) const override; | 
|  | bool Minus(const double* y, | 
|  | const double* x, | 
|  | double* y_minus_x) const override; | 
|  | bool MinusJacobian(const double* x, double* jacobian) const override; | 
|  |  | 
|  | private: | 
|  | const int tangent_size_ = 0; | 
|  | std::vector<bool> constancy_mask_; | 
|  | }; | 
|  |  | 
|  | // Implements the manifold for a Hamilton quaternion as defined in | 
|  | // https://en.wikipedia.org/wiki/Quaternion. Quaternions are represented as | 
|  | // unit norm 4-vectors, i.e. | 
|  | // | 
|  | // q = [q0; q1; q2; q3], |q| = 1 | 
|  | // | 
|  | // is the ambient space representation. | 
|  | // | 
|  | //   q0  scalar part. | 
|  | //   q1  coefficient of i. | 
|  | //   q2  coefficient of j. | 
|  | //   q3  coefficient of k. | 
|  | // | 
|  | // where: i*i = j*j = k*k = -1 and i*j = k, j*k = i, k*i = j. | 
|  | // | 
|  | // The tangent space is R^3, which relates to the ambient space through the | 
|  | // Plus and Minus operations defined as: | 
|  | // | 
|  | // Plus(x, delta) = [cos(|delta|); sin(|delta|) * delta / |delta|] * x | 
|  | //    Minus(y, x) = to_delta(y * x^{-1}) | 
|  | // | 
|  | // where "*" is the quaternion product and because q is a unit quaternion | 
|  | // (|q|=1), q^-1 = [q0; -q1; -q2; -q3] | 
|  | // | 
|  | // and to_delta( [q0; u_{3x1}] ) = u / |u| * atan2(|u|, q0) | 
|  | class CERES_EXPORT QuaternionManifold final : public Manifold { | 
|  | public: | 
|  | int AmbientSize() const override { return 4; } | 
|  | int TangentSize() const override { return 3; } | 
|  |  | 
|  | bool Plus(const double* x, | 
|  | const double* delta, | 
|  | double* x_plus_delta) const override; | 
|  | bool PlusJacobian(const double* x, double* jacobian) const override; | 
|  | bool Minus(const double* y, | 
|  | const double* x, | 
|  | double* y_minus_x) const override; | 
|  | bool MinusJacobian(const double* x, double* jacobian) const override; | 
|  | }; | 
|  |  | 
|  | // Implements the quaternion manifold for Eigen's representation of the | 
|  | // Hamilton quaternion. Geometrically it is exactly the same as the | 
|  | // QuaternionManifold defined above. However, Eigen uses a different internal | 
|  | // memory layout for the elements of the quaternion than what is commonly | 
|  | // used. It stores the quaternion in memory as [q1, q2, q3, q0] or | 
|  | // [x, y, z, w] where the real (scalar) part is last. | 
|  | // | 
|  | // Since Ceres operates on parameter blocks which are raw double pointers this | 
|  | // difference is important and requires a different manifold. | 
|  | class CERES_EXPORT EigenQuaternionManifold final : public Manifold { | 
|  | public: | 
|  | int AmbientSize() const override { return 4; } | 
|  | int TangentSize() const override { return 3; } | 
|  |  | 
|  | bool Plus(const double* x, | 
|  | const double* delta, | 
|  | double* x_plus_delta) const override; | 
|  | bool PlusJacobian(const double* x, double* jacobian) const override; | 
|  | bool Minus(const double* y, | 
|  | const double* x, | 
|  | double* y_minus_x) const override; | 
|  | bool MinusJacobian(const double* x, double* jacobian) const override; | 
|  | }; | 
|  |  | 
|  | }  // namespace ceres | 
|  |  | 
|  | // clang-format off | 
|  | #include "ceres/internal/reenable_warnings.h" | 
|  | // clang-format on | 
|  |  | 
|  | #endif  // CERES_PUBLIC_MANIFOLD_H_ |