123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319 |
- // Ceres Solver - A fast non-linear least squares minimizer
- // Copyright 2023 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)
- // sergiu.deitsch@gmail.com (Sergiu Deitsch)
- //
- #ifndef CERES_PUBLIC_PRODUCT_MANIFOLD_H_
- #define CERES_PUBLIC_PRODUCT_MANIFOLD_H_
- #include <algorithm>
- #include <array>
- #include <cassert>
- #include <cstddef>
- #include <numeric>
- #include <tuple>
- #include <type_traits>
- #include <utility>
- #include "ceres/internal/eigen.h"
- #include "ceres/internal/fixed_array.h"
- #include "ceres/internal/port.h"
- #include "ceres/manifold.h"
- namespace ceres {
- // Construct a manifold by taking the Cartesian product of a number of other
- // manifolds. This is useful, when a parameter block is the Cartesian product
- // of two or more manifolds. For example the parameters of a camera consist of
- // a rotation and a translation, i.e., SO(3) x R^3.
- //
- // Example usage:
- //
- // ProductManifold<QuaternionManifold, EuclideanManifold<3>> se3;
- //
- // is the manifold for a rigid transformation, where the rotation is
- // represented using a quaternion.
- //
- // Manifolds can be copied and moved to ProductManifold:
- //
- // SubsetManifold manifold1(5, {2});
- // SubsetManifold manifold2(3, {0, 1});
- // ProductManifold<SubsetManifold, SubsetManifold> manifold(manifold1,
- // manifold2);
- //
- // In advanced use cases, manifolds can be dynamically allocated and passed as
- // (smart) pointers:
- //
- // ProductManifold<std::unique_ptr<QuaternionManifold>, EuclideanManifold<3>>
- // se3{std::make_unique<QuaternionManifold>(), EuclideanManifold<3>{}};
- //
- // In C++17, the template parameters can be left out as they are automatically
- // deduced making the initialization much simpler:
- //
- // ProductManifold se3{QuaternionManifold{}, EuclideanManifold<3>{}};
- //
- // The manifold implementations must be either default constructible, copyable
- // or moveable to be usable in a ProductManifold.
- template <typename Manifold0, typename Manifold1, typename... ManifoldN>
- class ProductManifold final : public Manifold {
- public:
- // ProductManifold constructor perfect forwards arguments to store manifolds.
- //
- // Either use default construction or if you need to copy or move-construct a
- // manifold instance, you need to pass an instance as an argument for all
- // types given as class template parameters.
- template <typename... Args,
- std::enable_if_t<std::is_constructible<
- std::tuple<Manifold0, Manifold1, ManifoldN...>,
- Args...>::value>* = nullptr>
- explicit ProductManifold(Args&&... manifolds)
- : ProductManifold{std::make_index_sequence<kNumManifolds>{},
- std::forward<Args>(manifolds)...} {}
- int AmbientSize() const override { return ambient_size_; }
- int TangentSize() const override { return tangent_size_; }
- bool Plus(const double* x,
- const double* delta,
- double* x_plus_delta) const override {
- return PlusImpl(
- x, delta, x_plus_delta, std::make_index_sequence<kNumManifolds>{});
- }
- bool Minus(const double* y,
- const double* x,
- double* y_minus_x) const override {
- return MinusImpl(
- y, x, y_minus_x, std::make_index_sequence<kNumManifolds>{});
- }
- bool PlusJacobian(const double* x, double* jacobian_ptr) const override {
- MatrixRef jacobian(jacobian_ptr, AmbientSize(), TangentSize());
- jacobian.setZero();
- internal::FixedArray<double> buffer(buffer_size_);
- return PlusJacobianImpl(
- x, jacobian, buffer, std::make_index_sequence<kNumManifolds>{});
- }
- bool MinusJacobian(const double* x, double* jacobian_ptr) const override {
- MatrixRef jacobian(jacobian_ptr, TangentSize(), AmbientSize());
- jacobian.setZero();
- internal::FixedArray<double> buffer(buffer_size_);
- return MinusJacobianImpl(
- x, jacobian, buffer, std::make_index_sequence<kNumManifolds>{});
- }
- private:
- static constexpr std::size_t kNumManifolds = 2 + sizeof...(ManifoldN);
- template <std::size_t... Indices, typename... Args>
- explicit ProductManifold(std::index_sequence<Indices...>, Args&&... manifolds)
- : manifolds_{std::forward<Args>(manifolds)...},
- buffer_size_{(std::max)(
- {(Dereference(std::get<Indices>(manifolds_)).TangentSize() *
- Dereference(std::get<Indices>(manifolds_)).AmbientSize())...})},
- ambient_sizes_{
- Dereference(std::get<Indices>(manifolds_)).AmbientSize()...},
- tangent_sizes_{
- Dereference(std::get<Indices>(manifolds_)).TangentSize()...},
- ambient_offsets_{ExclusiveScan(ambient_sizes_)},
- tangent_offsets_{ExclusiveScan(tangent_sizes_)},
- ambient_size_{
- std::accumulate(ambient_sizes_.begin(), ambient_sizes_.end(), 0)},
- tangent_size_{
- std::accumulate(tangent_sizes_.begin(), tangent_sizes_.end(), 0)} {}
- template <std::size_t Index0, std::size_t... Indices>
- bool PlusImpl(const double* x,
- const double* delta,
- double* x_plus_delta,
- std::index_sequence<Index0, Indices...>) const {
- if (!Dereference(std::get<Index0>(manifolds_))
- .Plus(x + ambient_offsets_[Index0],
- delta + tangent_offsets_[Index0],
- x_plus_delta + ambient_offsets_[Index0])) {
- return false;
- }
- return PlusImpl(x, delta, x_plus_delta, std::index_sequence<Indices...>{});
- }
- static constexpr bool PlusImpl(const double* /*x*/,
- const double* /*delta*/,
- double* /*x_plus_delta*/,
- std::index_sequence<>) noexcept {
- return true;
- }
- template <std::size_t Index0, std::size_t... Indices>
- bool MinusImpl(const double* y,
- const double* x,
- double* y_minus_x,
- std::index_sequence<Index0, Indices...>) const {
- if (!Dereference(std::get<Index0>(manifolds_))
- .Minus(y + ambient_offsets_[Index0],
- x + ambient_offsets_[Index0],
- y_minus_x + tangent_offsets_[Index0])) {
- return false;
- }
- return MinusImpl(y, x, y_minus_x, std::index_sequence<Indices...>{});
- }
- static constexpr bool MinusImpl(const double* /*y*/,
- const double* /*x*/,
- double* /*y_minus_x*/,
- std::index_sequence<>) noexcept {
- return true;
- }
- template <std::size_t Index0, std::size_t... Indices>
- bool PlusJacobianImpl(const double* x,
- MatrixRef& jacobian,
- internal::FixedArray<double>& buffer,
- std::index_sequence<Index0, Indices...>) const {
- if (!Dereference(std::get<Index0>(manifolds_))
- .PlusJacobian(x + ambient_offsets_[Index0], buffer.data())) {
- return false;
- }
- jacobian.block(ambient_offsets_[Index0],
- tangent_offsets_[Index0],
- ambient_sizes_[Index0],
- tangent_sizes_[Index0]) =
- MatrixRef(
- buffer.data(), ambient_sizes_[Index0], tangent_sizes_[Index0]);
- return PlusJacobianImpl(
- x, jacobian, buffer, std::index_sequence<Indices...>{});
- }
- static constexpr bool PlusJacobianImpl(
- const double* /*x*/,
- MatrixRef& /*jacobian*/,
- internal::FixedArray<double>& /*buffer*/,
- std::index_sequence<>) noexcept {
- return true;
- }
- template <std::size_t Index0, std::size_t... Indices>
- bool MinusJacobianImpl(const double* x,
- MatrixRef& jacobian,
- internal::FixedArray<double>& buffer,
- std::index_sequence<Index0, Indices...>) const {
- if (!Dereference(std::get<Index0>(manifolds_))
- .MinusJacobian(x + ambient_offsets_[Index0], buffer.data())) {
- return false;
- }
- jacobian.block(tangent_offsets_[Index0],
- ambient_offsets_[Index0],
- tangent_sizes_[Index0],
- ambient_sizes_[Index0]) =
- MatrixRef(
- buffer.data(), tangent_sizes_[Index0], ambient_sizes_[Index0]);
- return MinusJacobianImpl(
- x, jacobian, buffer, std::index_sequence<Indices...>{});
- }
- static constexpr bool MinusJacobianImpl(
- const double* /*x*/,
- MatrixRef& /*jacobian*/,
- internal::FixedArray<double>& /*buffer*/,
- std::index_sequence<>) noexcept {
- return true;
- }
- template <typename T, std::size_t N>
- static std::array<T, N> ExclusiveScan(const std::array<T, N>& values) {
- std::array<T, N> result;
- // TODO Replace with std::exclusive_scan once all platforms have full C++17
- // STL support.
- T init = 0;
- for (std::size_t i = 0; i != N; ++i) {
- result[i] = init;
- init += values[i];
- }
- return result;
- }
- template <typename T, typename E = void>
- struct IsDereferenceable : std::false_type {};
- template <typename T>
- struct IsDereferenceable<T, std::void_t<decltype(*std::declval<T>())>>
- : std::true_type {};
- template <typename T,
- std::enable_if_t<!IsDereferenceable<T>::value>* = nullptr>
- static constexpr decltype(auto) Dereference(T& value) {
- return value;
- }
- // Support dereferenceable types such as std::unique_ptr, std::shared_ptr, raw
- // pointers etc.
- template <typename T,
- std::enable_if_t<IsDereferenceable<T>::value>* = nullptr>
- static constexpr decltype(auto) Dereference(T& value) {
- return *value;
- }
- template <typename T>
- static constexpr decltype(auto) Dereference(T* p) {
- assert(p != nullptr);
- return *p;
- }
- std::tuple<Manifold0, Manifold1, ManifoldN...> manifolds_;
- int buffer_size_;
- std::array<int, kNumManifolds> ambient_sizes_;
- std::array<int, kNumManifolds> tangent_sizes_;
- std::array<int, kNumManifolds> ambient_offsets_;
- std::array<int, kNumManifolds> tangent_offsets_;
- int ambient_size_;
- int tangent_size_;
- };
- // C++17 deduction guide that allows the user to avoid explicitly specifying
- // the template parameters of ProductManifold. The class can instead be
- // instantiated as follows:
- //
- // ProductManifold manifold{QuaternionManifold{}, EuclideanManifold<3>{}};
- //
- template <typename Manifold0, typename Manifold1, typename... Manifolds>
- ProductManifold(Manifold0&&, Manifold1&&, Manifolds&&...)
- -> ProductManifold<Manifold0, Manifold1, Manifolds...>;
- } // namespace ceres
- #endif // CERES_PUBLIC_PRODUCT_MANIFOLD_H_
|