sphere_manifold.h 9.1 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236
  1. // Ceres Solver - A fast non-linear least squares minimizer
  2. // Copyright 2023 Google Inc. All rights reserved.
  3. // http://ceres-solver.org/
  4. //
  5. // Redistribution and use in source and binary forms, with or without
  6. // modification, are permitted provided that the following conditions are met:
  7. //
  8. // * Redistributions of source code must retain the above copyright notice,
  9. // this list of conditions and the following disclaimer.
  10. // * Redistributions in binary form must reproduce the above copyright notice,
  11. // this list of conditions and the following disclaimer in the documentation
  12. // and/or other materials provided with the distribution.
  13. // * Neither the name of Google Inc. nor the names of its contributors may be
  14. // used to endorse or promote products derived from this software without
  15. // specific prior written permission.
  16. //
  17. // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
  18. // AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
  19. // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
  20. // ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
  21. // LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
  22. // CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
  23. // SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
  24. // INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
  25. // CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
  26. // ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
  27. // POSSIBILITY OF SUCH DAMAGE.
  28. //
  29. // Author: vitus@google.com (Mike Vitus)
  30. // jodebo_beck@gmx.de (Johannes Beck)
  31. #ifndef CERES_PUBLIC_SPHERE_MANIFOLD_H_
  32. #define CERES_PUBLIC_SPHERE_MANIFOLD_H_
  33. #include <Eigen/Core>
  34. #include <algorithm>
  35. #include <array>
  36. #include <memory>
  37. #include <vector>
  38. #include "ceres/internal/disable_warnings.h"
  39. #include "ceres/internal/export.h"
  40. #include "ceres/internal/householder_vector.h"
  41. #include "ceres/internal/sphere_manifold_functions.h"
  42. #include "ceres/manifold.h"
  43. #include "ceres/types.h"
  44. #include "glog/logging.h"
  45. namespace ceres {
  46. // This provides a manifold on a sphere meaning that the norm of the vector
  47. // stays the same. Such cases often arises in Structure for Motion
  48. // problems. One example where they are used is in representing points whose
  49. // triangulation is ill-conditioned. Here it is advantageous to use an
  50. // over-parameterization since homogeneous vectors can represent points at
  51. // infinity.
  52. //
  53. // The plus operator is defined as
  54. // Plus(x, delta) =
  55. // [sin(0.5 * |delta|) * delta / |delta|, cos(0.5 * |delta|)] * x
  56. //
  57. // The minus operator is defined as
  58. // Minus(x, y) = 2 atan2(nhy, y[-1]) / nhy * hy[0 : size_ - 1]
  59. // with nhy = norm(hy[0 : size_ - 1])
  60. //
  61. // with * defined as an operator which applies the update orthogonal to x to
  62. // remain on the sphere. The ambient space dimension is required to be greater
  63. // than 1.
  64. //
  65. // The class works with dynamic and static ambient space dimensions. If the
  66. // ambient space dimensions is known at compile time use
  67. //
  68. // SphereManifold<3> manifold;
  69. //
  70. // If the ambient space dimensions is not known at compile time the template
  71. // parameter needs to be set to ceres::DYNAMIC and the actual dimension needs
  72. // to be provided as a constructor argument:
  73. //
  74. // SphereManifold<ceres::DYNAMIC> manifold(ambient_dim);
  75. //
  76. // See section B.2 (p.25) in "Integrating Generic Sensor Fusion Algorithms
  77. // with Sound State Representations through Encapsulation of Manifolds" by C.
  78. // Hertzberg, R. Wagner, U. Frese and L. Schroder for more details
  79. // (https://arxiv.org/pdf/1107.1119.pdf)
  80. template <int AmbientSpaceDimension>
  81. class SphereManifold final : public Manifold {
  82. public:
  83. static_assert(
  84. AmbientSpaceDimension == ceres::DYNAMIC || AmbientSpaceDimension > 1,
  85. "The size of the homogeneous vector needs to be greater than 1.");
  86. static_assert(ceres::DYNAMIC == Eigen::Dynamic,
  87. "ceres::DYNAMIC needs to be the same as Eigen::Dynamic.");
  88. SphereManifold();
  89. explicit SphereManifold(int size);
  90. int AmbientSize() const override {
  91. return AmbientSpaceDimension == ceres::DYNAMIC ? size_
  92. : AmbientSpaceDimension;
  93. }
  94. int TangentSize() const override { return AmbientSize() - 1; }
  95. bool Plus(const double* x,
  96. const double* delta,
  97. double* x_plus_delta) const override;
  98. bool PlusJacobian(const double* x, double* jacobian) const override;
  99. bool Minus(const double* y,
  100. const double* x,
  101. double* y_minus_x) const override;
  102. bool MinusJacobian(const double* x, double* jacobian) const override;
  103. private:
  104. static constexpr int TangentSpaceDimension =
  105. AmbientSpaceDimension > 0 ? AmbientSpaceDimension - 1 : Eigen::Dynamic;
  106. // NOTE: Eigen does not allow to have a RowMajor column vector.
  107. // In that case, change the storage order
  108. static constexpr int SafeRowMajor =
  109. TangentSpaceDimension == 1 ? Eigen::ColMajor : Eigen::RowMajor;
  110. using AmbientVector = Eigen::Matrix<double, AmbientSpaceDimension, 1>;
  111. using TangentVector = Eigen::Matrix<double, TangentSpaceDimension, 1>;
  112. using MatrixPlusJacobian = Eigen::Matrix<double,
  113. AmbientSpaceDimension,
  114. TangentSpaceDimension,
  115. SafeRowMajor>;
  116. using MatrixMinusJacobian = Eigen::Matrix<double,
  117. TangentSpaceDimension,
  118. AmbientSpaceDimension,
  119. Eigen::RowMajor>;
  120. const int size_{};
  121. };
  122. template <int AmbientSpaceDimension>
  123. SphereManifold<AmbientSpaceDimension>::SphereManifold()
  124. : size_{AmbientSpaceDimension} {
  125. static_assert(
  126. AmbientSpaceDimension != Eigen::Dynamic,
  127. "The size is set to dynamic. Please call the constructor with a size.");
  128. }
  129. template <int AmbientSpaceDimension>
  130. SphereManifold<AmbientSpaceDimension>::SphereManifold(int size) : size_{size} {
  131. if (AmbientSpaceDimension != Eigen::Dynamic) {
  132. CHECK_EQ(AmbientSpaceDimension, size)
  133. << "Specified size by template parameter differs from the supplied "
  134. "one.";
  135. } else {
  136. CHECK_GT(size_, 1)
  137. << "The size of the manifold needs to be greater than 1.";
  138. }
  139. }
  140. template <int AmbientSpaceDimension>
  141. bool SphereManifold<AmbientSpaceDimension>::Plus(
  142. const double* x_ptr,
  143. const double* delta_ptr,
  144. double* x_plus_delta_ptr) const {
  145. Eigen::Map<const AmbientVector> x(x_ptr, size_);
  146. Eigen::Map<const TangentVector> delta(delta_ptr, size_ - 1);
  147. Eigen::Map<AmbientVector> x_plus_delta(x_plus_delta_ptr, size_);
  148. const double norm_delta = delta.norm();
  149. if (norm_delta == 0.0) {
  150. x_plus_delta = x;
  151. return true;
  152. }
  153. AmbientVector v(size_);
  154. double beta;
  155. // NOTE: The explicit template arguments are needed here because
  156. // ComputeHouseholderVector is templated and some versions of MSVC
  157. // have trouble deducing the type of v automatically.
  158. internal::ComputeHouseholderVector<Eigen::Map<const AmbientVector>,
  159. double,
  160. AmbientSpaceDimension>(x, &v, &beta);
  161. internal::ComputeSphereManifoldPlus(
  162. v, beta, x, delta, norm_delta, &x_plus_delta);
  163. return true;
  164. }
  165. template <int AmbientSpaceDimension>
  166. bool SphereManifold<AmbientSpaceDimension>::PlusJacobian(
  167. const double* x_ptr, double* jacobian_ptr) const {
  168. Eigen::Map<const AmbientVector> x(x_ptr, size_);
  169. Eigen::Map<MatrixPlusJacobian> jacobian(jacobian_ptr, size_, size_ - 1);
  170. internal::ComputeSphereManifoldPlusJacobian(x, &jacobian);
  171. return true;
  172. }
  173. template <int AmbientSpaceDimension>
  174. bool SphereManifold<AmbientSpaceDimension>::Minus(const double* y_ptr,
  175. const double* x_ptr,
  176. double* y_minus_x_ptr) const {
  177. AmbientVector y = Eigen::Map<const AmbientVector>(y_ptr, size_);
  178. Eigen::Map<const AmbientVector> x(x_ptr, size_);
  179. Eigen::Map<TangentVector> y_minus_x(y_minus_x_ptr, size_ - 1);
  180. // Apply hoseholder transformation.
  181. AmbientVector v(size_);
  182. double beta;
  183. // NOTE: The explicit template arguments are needed here because
  184. // ComputeHouseholderVector is templated and some versions of MSVC
  185. // have trouble deducing the type of v automatically.
  186. internal::ComputeHouseholderVector<Eigen::Map<const AmbientVector>,
  187. double,
  188. AmbientSpaceDimension>(x, &v, &beta);
  189. internal::ComputeSphereManifoldMinus(v, beta, x, y, &y_minus_x);
  190. return true;
  191. }
  192. template <int AmbientSpaceDimension>
  193. bool SphereManifold<AmbientSpaceDimension>::MinusJacobian(
  194. const double* x_ptr, double* jacobian_ptr) const {
  195. Eigen::Map<const AmbientVector> x(x_ptr, size_);
  196. Eigen::Map<MatrixMinusJacobian> jacobian(jacobian_ptr, size_ - 1, size_);
  197. internal::ComputeSphereManifoldMinusJacobian(x, &jacobian);
  198. return true;
  199. }
  200. } // namespace ceres
  201. // clang-format off
  202. #include "ceres/internal/reenable_warnings.h"
  203. // clang-format on
  204. #endif // CERES_PUBLIC_SPHERE_MANIFOLD_H_