1 // Ceres Solver - A fast non-linear least squares minimizer 2 // Copyright 2010, 2011, 2012 Google Inc. All rights reserved. 3 // http://code.google.com/p/ceres-solver/ 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 (at) 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 48 #ifndef CERES_INTERNAL_VISIBILITY_BASED_PRECONDITIONER_H_ 49 #define CERES_INTERNAL_VISIBILITY_BASED_PRECONDITIONER_H_ 50 51 #include <set> 52 #include <vector> 53 #include <utility> 54 #include "ceres/collections_port.h" 55 #include "ceres/graph.h" 56 #include "ceres/internal/macros.h" 57 #include "ceres/internal/scoped_ptr.h" 58 #include "ceres/linear_solver.h" 59 #include "ceres/preconditioner.h" 60 #include "ceres/suitesparse.h" 61 62 namespace ceres { 63 namespace internal { 64 65 class BlockRandomAccessSparseMatrix; 66 class BlockSparseMatrix; 67 struct CompressedRowBlockStructure; 68 class SchurEliminatorBase; 69 70 // This class implements visibility based preconditioners for 71 // Structure from Motion/Bundle Adjustment problems. The name 72 // VisibilityBasedPreconditioner comes from the fact that the sparsity 73 // structure of the preconditioner matrix is determined by analyzing 74 // the visibility structure of the scene, i.e. which cameras see which 75 // points. 76 // 77 // The key idea of visibility based preconditioning is to identify 78 // cameras that we expect have strong interactions, and then using the 79 // entries in the Schur complement matrix corresponding to these 80 // camera pairs as an approximation to the full Schur complement. 81 // 82 // CLUSTER_JACOBI identifies these camera pairs by clustering cameras, 83 // and considering all non-zero camera pairs within each cluster. The 84 // clustering in the current implementation is done using the 85 // Canonical Views algorithm of Simon et al. (see 86 // canonical_views_clustering.h). For the purposes of clustering, the 87 // similarity or the degree of interaction between a pair of cameras 88 // is measured by counting the number of points visible in both the 89 // cameras. Thus the name VisibilityBasedPreconditioner. Further, if we 90 // were to permute the parameter blocks such that all the cameras in 91 // the same cluster occur contiguously, the preconditioner matrix will 92 // be a block diagonal matrix with blocks corresponding to the 93 // clusters. Thus in analogy with the Jacobi preconditioner we refer 94 // to this as the CLUSTER_JACOBI preconditioner. 95 // 96 // CLUSTER_TRIDIAGONAL adds more mass to the CLUSTER_JACOBI 97 // preconditioner by considering the interaction between clusters and 98 // identifying strong interactions between cluster pairs. This is done 99 // by constructing a weighted graph on the clusters, with the weight 100 // on the edges connecting two clusters proportional to the number of 101 // 3D points visible to cameras in both the clusters. A degree-2 102 // maximum spanning forest is identified in this graph and the camera 103 // pairs contained in the edges of this forest are added to the 104 // preconditioner. The detailed reasoning for this construction is 105 // explained in the paper mentioned above. 106 // 107 // Degree-2 spanning trees and forests have the property that they 108 // correspond to tri-diagonal matrices. Thus there exist a permutation 109 // of the camera blocks under which the CLUSTER_TRIDIAGONAL 110 // preconditioner matrix is a block tridiagonal matrix, and thus the 111 // name for the preconditioner. 112 // 113 // Thread Safety: This class is NOT thread safe. 114 // 115 // Example usage: 116 // 117 // LinearSolver::Options options; 118 // options.preconditioner_type = CLUSTER_JACOBI; 119 // options.elimination_groups.push_back(num_points); 120 // options.elimination_groups.push_back(num_cameras); 121 // VisibilityBasedPreconditioner preconditioner( 122 // *A.block_structure(), options); 123 // preconditioner.Update(A, NULL); 124 // preconditioner.RightMultiply(x, y); 125 // 126 #ifndef CERES_NO_SUITESPARSE 127 class VisibilityBasedPreconditioner : public BlockSparseMatrixPreconditioner { 128 public: 129 // Initialize the symbolic structure of the preconditioner. bs is 130 // the block structure of the linear system to be solved. It is used 131 // to determine the sparsity structure of the preconditioner matrix. 132 // 133 // It has the same structural requirement as other Schur complement 134 // based solvers. Please see schur_eliminator.h for more details. 135 VisibilityBasedPreconditioner(const CompressedRowBlockStructure& bs, 136 const Preconditioner::Options& options); 137 virtual ~VisibilityBasedPreconditioner(); 138 139 // Preconditioner interface 140 virtual void RightMultiply(const double* x, double* y) const; 141 virtual int num_rows() const; 142 143 friend class VisibilityBasedPreconditionerTest; 144 145 private: 146 virtual bool UpdateImpl(const BlockSparseMatrix& A, const double* D); 147 void ComputeClusterJacobiSparsity(const CompressedRowBlockStructure& bs); 148 void ComputeClusterTridiagonalSparsity(const CompressedRowBlockStructure& bs); 149 void InitStorage(const CompressedRowBlockStructure& bs); 150 void InitEliminator(const CompressedRowBlockStructure& bs); 151 LinearSolverTerminationType Factorize(); 152 void ScaleOffDiagonalCells(); 153 154 void ClusterCameras(const vector< set<int> >& visibility); 155 void FlattenMembershipMap(const HashMap<int, int>& membership_map, 156 vector<int>* membership_vector) const; 157 void ComputeClusterVisibility(const vector<set<int> >& visibility, 158 vector<set<int> >* cluster_visibility) const; 159 Graph<int>* CreateClusterGraph(const vector<set<int> >& visibility) const; 160 void ForestToClusterPairs(const Graph<int>& forest, 161 HashSet<pair<int, int> >* cluster_pairs) const; 162 void ComputeBlockPairsInPreconditioner(const CompressedRowBlockStructure& bs); 163 bool IsBlockPairInPreconditioner(int block1, int block2) const; 164 bool IsBlockPairOffDiagonal(int block1, int block2) const; 165 166 Preconditioner::Options options_; 167 168 // Number of parameter blocks in the schur complement. 169 int num_blocks_; 170 int num_clusters_; 171 172 // Sizes of the blocks in the schur complement. 173 vector<int> block_size_; 174 175 // Mapping from cameras to clusters. 176 vector<int> cluster_membership_; 177 178 // Non-zero camera pairs from the schur complement matrix that are 179 // present in the preconditioner, sorted by row (first element of 180 // each pair), then column (second). 181 set<pair<int, int> > block_pairs_; 182 183 // Set of cluster pairs (including self pairs (i,i)) in the 184 // preconditioner. 185 HashSet<pair<int, int> > cluster_pairs_; 186 scoped_ptr<SchurEliminatorBase> eliminator_; 187 188 // Preconditioner matrix. 189 scoped_ptr<BlockRandomAccessSparseMatrix> m_; 190 191 // RightMultiply is a const method for LinearOperators. It is 192 // implemented using CHOLMOD's sparse triangular matrix solve 193 // function. This however requires non-const access to the 194 // SuiteSparse context object, even though it does not result in any 195 // of the state of the preconditioner being modified. 196 SuiteSparse ss_; 197 198 // Symbolic and numeric factorization of the preconditioner. 199 cholmod_factor* factor_; 200 201 // Temporary vector used by RightMultiply. 202 cholmod_dense* tmp_rhs_; 203 CERES_DISALLOW_COPY_AND_ASSIGN(VisibilityBasedPreconditioner); 204 }; 205 #else // SuiteSparse 206 // If SuiteSparse is not compiled in, the preconditioner is not 207 // available. 208 class VisibilityBasedPreconditioner : public BlockSparseMatrixPreconditioner { 209 public: 210 VisibilityBasedPreconditioner(const CompressedRowBlockStructure& bs, 211 const Preconditioner::Options& options) { 212 LOG(FATAL) << "Visibility based preconditioning is not available. Please " 213 "build Ceres with SuiteSparse."; 214 } 215 virtual ~VisibilityBasedPreconditioner() {} 216 virtual void RightMultiply(const double* x, double* y) const {} 217 virtual void LeftMultiply(const double* x, double* y) const {} 218 virtual int num_rows() const { return -1; } 219 virtual int num_cols() const { return -1; } 220 221 private: 222 bool UpdateImpl(const BlockSparseMatrix& A, const double* D) { 223 return false; 224 } 225 }; 226 #endif // CERES_NO_SUITESPARSE 227 228 } // namespace internal 229 } // namespace ceres 230 231 #endif // CERES_INTERNAL_VISIBILITY_BASED_PRECONDITIONER_H_ 232
