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@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 62namespace ceres { 63namespace internal { 64 65class BlockRandomAccessSparseMatrix; 66class BlockSparseMatrix; 67struct CompressedRowBlockStructure; 68class 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 127class 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. 208class 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