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