visibility_based_preconditioner.h 8.8 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203
  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: sameeragarwal@google.com (Sameer Agarwal)
  30. //
  31. // Preconditioners for linear systems that arise in Structure from
  32. // Motion problems. VisibilityBasedPreconditioner implements:
  33. //
  34. // CLUSTER_JACOBI
  35. // CLUSTER_TRIDIAGONAL
  36. //
  37. // Detailed descriptions of these preconditions beyond what is
  38. // documented here can be found in
  39. //
  40. // Visibility Based Preconditioning for Bundle Adjustment
  41. // A. Kushal & S. Agarwal, CVPR 2012.
  42. //
  43. // http://www.cs.washington.edu/homes/sagarwal/vbp.pdf
  44. //
  45. // The two preconditioners share enough code that its most efficient
  46. // to implement them as part of the same code base.
  47. #ifndef CERES_INTERNAL_VISIBILITY_BASED_PRECONDITIONER_H_
  48. #define CERES_INTERNAL_VISIBILITY_BASED_PRECONDITIONER_H_
  49. #include <memory>
  50. #include <set>
  51. #include <unordered_map>
  52. #include <unordered_set>
  53. #include <utility>
  54. #include <vector>
  55. #include "ceres/block_structure.h"
  56. #include "ceres/graph.h"
  57. #include "ceres/linear_solver.h"
  58. #include "ceres/pair_hash.h"
  59. #include "ceres/preconditioner.h"
  60. #include "ceres/sparse_cholesky.h"
  61. namespace ceres::internal {
  62. class BlockRandomAccessSparseMatrix;
  63. class BlockSparseMatrix;
  64. struct CompressedRowBlockStructure;
  65. class SchurEliminatorBase;
  66. // This class implements visibility based preconditioners for
  67. // Structure from Motion/Bundle Adjustment problems. The name
  68. // VisibilityBasedPreconditioner comes from the fact that the sparsity
  69. // structure of the preconditioner matrix is determined by analyzing
  70. // the visibility structure of the scene, i.e. which cameras see which
  71. // points.
  72. //
  73. // The key idea of visibility based preconditioning is to identify
  74. // cameras that we expect have strong interactions, and then using the
  75. // entries in the Schur complement matrix corresponding to these
  76. // camera pairs as an approximation to the full Schur complement.
  77. //
  78. // CLUSTER_JACOBI identifies these camera pairs by clustering cameras,
  79. // and considering all non-zero camera pairs within each cluster. The
  80. // clustering in the current implementation is done using the
  81. // Canonical Views algorithm of Simon et al. (see
  82. // canonical_views_clustering.h). For the purposes of clustering, the
  83. // similarity or the degree of interaction between a pair of cameras
  84. // is measured by counting the number of points visible in both the
  85. // cameras. Thus the name VisibilityBasedPreconditioner. Further, if we
  86. // were to permute the parameter blocks such that all the cameras in
  87. // the same cluster occur contiguously, the preconditioner matrix will
  88. // be a block diagonal matrix with blocks corresponding to the
  89. // clusters. Thus in analogy with the Jacobi preconditioner we refer
  90. // to this as the CLUSTER_JACOBI preconditioner.
  91. //
  92. // CLUSTER_TRIDIAGONAL adds more mass to the CLUSTER_JACOBI
  93. // preconditioner by considering the interaction between clusters and
  94. // identifying strong interactions between cluster pairs. This is done
  95. // by constructing a weighted graph on the clusters, with the weight
  96. // on the edges connecting two clusters proportional to the number of
  97. // 3D points visible to cameras in both the clusters. A degree-2
  98. // maximum spanning forest is identified in this graph and the camera
  99. // pairs contained in the edges of this forest are added to the
  100. // preconditioner. The detailed reasoning for this construction is
  101. // explained in the paper mentioned above.
  102. //
  103. // Degree-2 spanning trees and forests have the property that they
  104. // correspond to tri-diagonal matrices. Thus there exist a permutation
  105. // of the camera blocks under which the CLUSTER_TRIDIAGONAL
  106. // preconditioner matrix is a block tridiagonal matrix, and thus the
  107. // name for the preconditioner.
  108. //
  109. // Thread Safety: This class is NOT thread safe.
  110. //
  111. // Example usage:
  112. //
  113. // LinearSolver::Options options;
  114. // options.preconditioner_type = CLUSTER_JACOBI;
  115. // options.elimination_groups.push_back(num_points);
  116. // options.elimination_groups.push_back(num_cameras);
  117. // VisibilityBasedPreconditioner preconditioner(
  118. // *A.block_structure(), options);
  119. // preconditioner.Update(A, nullptr);
  120. // preconditioner.RightMultiplyAndAccumulate(x, y);
  121. class CERES_NO_EXPORT VisibilityBasedPreconditioner
  122. : public BlockSparseMatrixPreconditioner {
  123. public:
  124. // Initialize the symbolic structure of the preconditioner. bs is
  125. // the block structure of the linear system to be solved. It is used
  126. // to determine the sparsity structure of the preconditioner matrix.
  127. //
  128. // It has the same structural requirement as other Schur complement
  129. // based solvers. Please see schur_eliminator.h for more details.
  130. VisibilityBasedPreconditioner(const CompressedRowBlockStructure& bs,
  131. Preconditioner::Options options);
  132. VisibilityBasedPreconditioner(const VisibilityBasedPreconditioner&) = delete;
  133. void operator=(const VisibilityBasedPreconditioner&) = delete;
  134. ~VisibilityBasedPreconditioner() override;
  135. // Preconditioner interface
  136. void RightMultiplyAndAccumulate(const double* x, double* y) const final;
  137. int num_rows() const final;
  138. friend class VisibilityBasedPreconditionerTest;
  139. private:
  140. bool UpdateImpl(const BlockSparseMatrix& A, const double* D) final;
  141. void ComputeClusterJacobiSparsity(const CompressedRowBlockStructure& bs);
  142. void ComputeClusterTridiagonalSparsity(const CompressedRowBlockStructure& bs);
  143. void InitStorage(const CompressedRowBlockStructure& bs);
  144. void InitEliminator(const CompressedRowBlockStructure& bs);
  145. LinearSolverTerminationType Factorize();
  146. void ScaleOffDiagonalCells();
  147. void ClusterCameras(const std::vector<std::set<int>>& visibility);
  148. void FlattenMembershipMap(const std::unordered_map<int, int>& membership_map,
  149. std::vector<int>* membership_vector) const;
  150. void ComputeClusterVisibility(
  151. const std::vector<std::set<int>>& visibility,
  152. std::vector<std::set<int>>* cluster_visibility) const;
  153. std::unique_ptr<WeightedGraph<int>> CreateClusterGraph(
  154. const std::vector<std::set<int>>& visibility) const;
  155. void ForestToClusterPairs(
  156. const WeightedGraph<int>& forest,
  157. std::unordered_set<std::pair<int, int>, pair_hash>* cluster_pairs) const;
  158. void ComputeBlockPairsInPreconditioner(const CompressedRowBlockStructure& bs);
  159. bool IsBlockPairInPreconditioner(int block1, int block2) const;
  160. bool IsBlockPairOffDiagonal(int block1, int block2) const;
  161. Preconditioner::Options options_;
  162. // Number of parameter blocks in the schur complement.
  163. int num_blocks_;
  164. int num_clusters_;
  165. // Sizes of the blocks in the schur complement.
  166. std::vector<Block> blocks_;
  167. // Mapping from cameras to clusters.
  168. std::vector<int> cluster_membership_;
  169. // Non-zero camera pairs from the schur complement matrix that are
  170. // present in the preconditioner, sorted by row (first element of
  171. // each pair), then column (second).
  172. std::set<std::pair<int, int>> block_pairs_;
  173. // Set of cluster pairs (including self pairs (i,i)) in the
  174. // preconditioner.
  175. std::unordered_set<std::pair<int, int>, pair_hash> cluster_pairs_;
  176. std::unique_ptr<SchurEliminatorBase> eliminator_;
  177. // Preconditioner matrix.
  178. std::unique_ptr<BlockRandomAccessSparseMatrix> m_;
  179. std::unique_ptr<CompressedRowSparseMatrix> m_crs_;
  180. std::unique_ptr<SparseCholesky> sparse_cholesky_;
  181. };
  182. } // namespace ceres::internal
  183. #endif // CERES_INTERNAL_VISIBILITY_BASED_PRECONDITIONER_H_