diff options
11 files changed, 5 insertions, 2034 deletions
diff --git a/extern/libmv/third_party/ceres/CMakeLists.txt b/extern/libmv/third_party/ceres/CMakeLists.txt index f9a7d0631c0..cc72faa392e 100644 --- a/extern/libmv/third_party/ceres/CMakeLists.txt +++ b/extern/libmv/third_party/ceres/CMakeLists.txt @@ -52,7 +52,6 @@ set(SRC internal/ceres/block_sparse_matrix.cc internal/ceres/block_structure.cc internal/ceres/callbacks.cc - internal/ceres/canonical_views_clustering.cc internal/ceres/c_api.cc internal/ceres/cgnr_solver.cc internal/ceres/compressed_col_sparse_matrix_utils.cc @@ -64,7 +63,6 @@ set(SRC internal/ceres/corrector.cc internal/ceres/covariance.cc internal/ceres/covariance_impl.cc - internal/ceres/cxsparse.cc internal/ceres/dense_normal_cholesky_solver.cc internal/ceres/dense_qr_solver.cc internal/ceres/dense_sparse_matrix.cc @@ -110,21 +108,17 @@ set(SRC internal/ceres/schur_eliminator.cc internal/ceres/schur_jacobi_preconditioner.cc internal/ceres/scratch_evaluate_preparer.cc - internal/ceres/single_linkage_clustering.cc internal/ceres/solver.cc internal/ceres/solver_utils.cc internal/ceres/sparse_matrix.cc internal/ceres/sparse_normal_cholesky_solver.cc internal/ceres/split.cc internal/ceres/stringprintf.cc - internal/ceres/suitesparse.cc internal/ceres/triplet_sparse_matrix.cc internal/ceres/trust_region_minimizer.cc internal/ceres/trust_region_preprocessor.cc internal/ceres/trust_region_strategy.cc internal/ceres/types.cc - internal/ceres/visibility_based_preconditioner.cc - internal/ceres/visibility.cc internal/ceres/wall_time.cc include/ceres/autodiff_cost_function.h @@ -178,7 +172,6 @@ set(SRC internal/ceres/block_sparse_matrix.h internal/ceres/block_structure.h internal/ceres/callbacks.h - internal/ceres/canonical_views_clustering.h internal/ceres/casts.h internal/ceres/cgnr_linear_operator.h internal/ceres/cgnr_solver.h @@ -242,7 +235,6 @@ set(SRC internal/ceres/schur_eliminator_impl.h internal/ceres/schur_jacobi_preconditioner.h internal/ceres/scratch_evaluate_preparer.h - internal/ceres/single_linkage_clustering.h internal/ceres/small_blas.h internal/ceres/solver_utils.h internal/ceres/sparse_matrix.h @@ -256,7 +248,6 @@ set(SRC internal/ceres/trust_region_preprocessor.h internal/ceres/trust_region_strategy.h internal/ceres/visibility_based_preconditioner.h - internal/ceres/visibility.h internal/ceres/wall_time.h ) diff --git a/extern/libmv/third_party/ceres/files.txt b/extern/libmv/third_party/ceres/files.txt index 976200b6229..f49f1fb0ded 100644 --- a/extern/libmv/third_party/ceres/files.txt +++ b/extern/libmv/third_party/ceres/files.txt @@ -8,6 +8,7 @@ include/ceres/cost_function_to_functor.h include/ceres/covariance.h include/ceres/crs_matrix.h include/ceres/dynamic_autodiff_cost_function.h +include/ceres/dynamic_cost_function_to_functor.h include/ceres/dynamic_numeric_diff_cost_function.h include/ceres/fpclassify.h include/ceres/gradient_checker.h @@ -30,6 +31,7 @@ include/ceres/local_parameterization.h include/ceres/loss_function.h include/ceres/normal_prior.h include/ceres/numeric_diff_cost_function.h +include/ceres/numeric_diff_options.h include/ceres/ordered_groups.h include/ceres/problem.h include/ceres/rotation.h @@ -61,8 +63,6 @@ internal/ceres/block_structure.cc internal/ceres/block_structure.h internal/ceres/callbacks.cc internal/ceres/callbacks.h -internal/ceres/canonical_views_clustering.cc -internal/ceres/canonical_views_clustering.h internal/ceres/c_api.cc internal/ceres/casts.h internal/ceres/cgnr_linear_operator.h @@ -85,7 +85,6 @@ internal/ceres/corrector.h internal/ceres/covariance.cc internal/ceres/covariance_impl.cc internal/ceres/covariance_impl.h -internal/ceres/cxsparse.cc internal/ceres/cxsparse.h internal/ceres/dense_jacobian_writer.h internal/ceres/dense_normal_cholesky_solver.cc @@ -114,6 +113,7 @@ internal/ceres/generated/partitioned_matrix_view_2_2_4.cc internal/ceres/generated/partitioned_matrix_view_2_2_d.cc internal/ceres/generated/partitioned_matrix_view_2_3_3.cc internal/ceres/generated/partitioned_matrix_view_2_3_4.cc +internal/ceres/generated/partitioned_matrix_view_2_3_6.cc internal/ceres/generated/partitioned_matrix_view_2_3_9.cc internal/ceres/generated/partitioned_matrix_view_2_3_d.cc internal/ceres/generated/partitioned_matrix_view_2_4_3.cc @@ -133,6 +133,7 @@ internal/ceres/generated/schur_eliminator_2_2_4.cc internal/ceres/generated/schur_eliminator_2_2_d.cc internal/ceres/generated/schur_eliminator_2_3_3.cc internal/ceres/generated/schur_eliminator_2_3_4.cc +internal/ceres/generated/schur_eliminator_2_3_6.cc internal/ceres/generated/schur_eliminator_2_3_9.cc internal/ceres/generated/schur_eliminator_2_3_d.cc internal/ceres/generated/schur_eliminator_2_4_3.cc @@ -155,6 +156,7 @@ internal/ceres/gradient_problem_evaluator.h internal/ceres/gradient_problem_solver.cc internal/ceres/graph_algorithms.h internal/ceres/graph.h +internal/ceres/householder_vector.h internal/ceres/implicit_schur_complement.cc internal/ceres/implicit_schur_complement.h internal/ceres/integral_types.h @@ -221,8 +223,6 @@ internal/ceres/schur_jacobi_preconditioner.cc internal/ceres/schur_jacobi_preconditioner.h internal/ceres/scratch_evaluate_preparer.cc internal/ceres/scratch_evaluate_preparer.h -internal/ceres/single_linkage_clustering.cc -internal/ceres/single_linkage_clustering.h internal/ceres/small_blas.h internal/ceres/solver.cc internal/ceres/solver_utils.cc @@ -236,7 +236,6 @@ internal/ceres/split.h internal/ceres/stl_util.h internal/ceres/stringprintf.cc internal/ceres/stringprintf.h -internal/ceres/suitesparse.cc internal/ceres/suitesparse.h internal/ceres/triplet_sparse_matrix.cc internal/ceres/triplet_sparse_matrix.h @@ -247,10 +246,7 @@ internal/ceres/trust_region_preprocessor.h internal/ceres/trust_region_strategy.cc internal/ceres/trust_region_strategy.h internal/ceres/types.cc -internal/ceres/visibility_based_preconditioner.cc internal/ceres/visibility_based_preconditioner.h -internal/ceres/visibility.cc -internal/ceres/visibility.h internal/ceres/wall_time.cc internal/ceres/wall_time.h config/ceres/internal/config.h diff --git a/extern/libmv/third_party/ceres/internal/ceres/canonical_views_clustering.cc b/extern/libmv/third_party/ceres/internal/ceres/canonical_views_clustering.cc deleted file mode 100644 index b655b1eccd3..00000000000 --- a/extern/libmv/third_party/ceres/internal/ceres/canonical_views_clustering.cc +++ /dev/null @@ -1,247 +0,0 @@ -// Ceres Solver - A fast non-linear least squares minimizer -// Copyright 2015 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: David Gallup (dgallup@google.com) -// Sameer Agarwal (sameeragarwal@google.com) - -// This include must come before any #ifndef check on Ceres compile options. -#include "ceres/internal/port.h" - -#ifndef CERES_NO_SUITESPARSE - -#include "ceres/canonical_views_clustering.h" - -#include "ceres/collections_port.h" -#include "ceres/graph.h" -#include "ceres/internal/macros.h" -#include "ceres/map_util.h" -#include "glog/logging.h" - -namespace ceres { -namespace internal { - -using std::vector; - -typedef HashMap<int, int> IntMap; -typedef HashSet<int> IntSet; - -class CanonicalViewsClustering { - public: - CanonicalViewsClustering() {} - - // Compute the canonical views clustering of the vertices of the - // graph. centers will contain the vertices that are the identified - // as the canonical views/cluster centers, and membership is a map - // from vertices to cluster_ids. The i^th cluster center corresponds - // to the i^th cluster. It is possible depending on the - // configuration of the clustering algorithm that some of the - // vertices may not be assigned to any cluster. In this case they - // are assigned to a cluster with id = kInvalidClusterId. - void ComputeClustering(const CanonicalViewsClusteringOptions& options, - const WeightedGraph<int>& graph, - vector<int>* centers, - IntMap* membership); - - private: - void FindValidViews(IntSet* valid_views) const; - double ComputeClusteringQualityDifference(const int candidate, - const vector<int>& centers) const; - void UpdateCanonicalViewAssignments(const int canonical_view); - void ComputeClusterMembership(const vector<int>& centers, - IntMap* membership) const; - - CanonicalViewsClusteringOptions options_; - const WeightedGraph<int>* graph_; - // Maps a view to its representative canonical view (its cluster - // center). - IntMap view_to_canonical_view_; - // Maps a view to its similarity to its current cluster center. - HashMap<int, double> view_to_canonical_view_similarity_; - CERES_DISALLOW_COPY_AND_ASSIGN(CanonicalViewsClustering); -}; - -void ComputeCanonicalViewsClustering( - const CanonicalViewsClusteringOptions& options, - const WeightedGraph<int>& graph, - vector<int>* centers, - IntMap* membership) { - time_t start_time = time(NULL); - CanonicalViewsClustering cv; - cv.ComputeClustering(options, graph, centers, membership); - VLOG(2) << "Canonical views clustering time (secs): " - << time(NULL) - start_time; -} - -// Implementation of CanonicalViewsClustering -void CanonicalViewsClustering::ComputeClustering( - const CanonicalViewsClusteringOptions& options, - const WeightedGraph<int>& graph, - vector<int>* centers, - IntMap* membership) { - options_ = options; - CHECK_NOTNULL(centers)->clear(); - CHECK_NOTNULL(membership)->clear(); - graph_ = &graph; - - IntSet valid_views; - FindValidViews(&valid_views); - while (valid_views.size() > 0) { - // Find the next best canonical view. - double best_difference = -std::numeric_limits<double>::max(); - int best_view = 0; - - // TODO(sameeragarwal): Make this loop multi-threaded. - for (IntSet::const_iterator view = valid_views.begin(); - view != valid_views.end(); - ++view) { - const double difference = - ComputeClusteringQualityDifference(*view, *centers); - if (difference > best_difference) { - best_difference = difference; - best_view = *view; - } - } - - CHECK_GT(best_difference, -std::numeric_limits<double>::max()); - - // Add canonical view if quality improves, or if minimum is not - // yet met, otherwise break. - if ((best_difference <= 0) && - (centers->size() >= options_.min_views)) { - break; - } - - centers->push_back(best_view); - valid_views.erase(best_view); - UpdateCanonicalViewAssignments(best_view); - } - - ComputeClusterMembership(*centers, membership); -} - -// Return the set of vertices of the graph which have valid vertex -// weights. -void CanonicalViewsClustering::FindValidViews( - IntSet* valid_views) const { - const IntSet& views = graph_->vertices(); - for (IntSet::const_iterator view = views.begin(); - view != views.end(); - ++view) { - if (graph_->VertexWeight(*view) != WeightedGraph<int>::InvalidWeight()) { - valid_views->insert(*view); - } - } -} - -// Computes the difference in the quality score if 'candidate' were -// added to the set of canonical views. -double CanonicalViewsClustering::ComputeClusteringQualityDifference( - const int candidate, - const vector<int>& centers) const { - // View score. - double difference = - options_.view_score_weight * graph_->VertexWeight(candidate); - - // Compute how much the quality score changes if the candidate view - // was added to the list of canonical views and its nearest - // neighbors became members of its cluster. - const IntSet& neighbors = graph_->Neighbors(candidate); - for (IntSet::const_iterator neighbor = neighbors.begin(); - neighbor != neighbors.end(); - ++neighbor) { - const double old_similarity = - FindWithDefault(view_to_canonical_view_similarity_, *neighbor, 0.0); - const double new_similarity = graph_->EdgeWeight(*neighbor, candidate); - if (new_similarity > old_similarity) { - difference += new_similarity - old_similarity; - } - } - - // Number of views penalty. - difference -= options_.size_penalty_weight; - - // Orthogonality. - for (int i = 0; i < centers.size(); ++i) { - difference -= options_.similarity_penalty_weight * - graph_->EdgeWeight(centers[i], candidate); - } - - return difference; -} - -// Reassign views if they're more similar to the new canonical view. -void CanonicalViewsClustering::UpdateCanonicalViewAssignments( - const int canonical_view) { - const IntSet& neighbors = graph_->Neighbors(canonical_view); - for (IntSet::const_iterator neighbor = neighbors.begin(); - neighbor != neighbors.end(); - ++neighbor) { - const double old_similarity = - FindWithDefault(view_to_canonical_view_similarity_, *neighbor, 0.0); - const double new_similarity = - graph_->EdgeWeight(*neighbor, canonical_view); - if (new_similarity > old_similarity) { - view_to_canonical_view_[*neighbor] = canonical_view; - view_to_canonical_view_similarity_[*neighbor] = new_similarity; - } - } -} - -// Assign a cluster id to each view. -void CanonicalViewsClustering::ComputeClusterMembership( - const vector<int>& centers, - IntMap* membership) const { - CHECK_NOTNULL(membership)->clear(); - - // The i^th cluster has cluster id i. - IntMap center_to_cluster_id; - for (int i = 0; i < centers.size(); ++i) { - center_to_cluster_id[centers[i]] = i; - } - - static const int kInvalidClusterId = -1; - - const IntSet& views = graph_->vertices(); - for (IntSet::const_iterator view = views.begin(); - view != views.end(); - ++view) { - IntMap::const_iterator it = - view_to_canonical_view_.find(*view); - int cluster_id = kInvalidClusterId; - if (it != view_to_canonical_view_.end()) { - cluster_id = FindOrDie(center_to_cluster_id, it->second); - } - - InsertOrDie(membership, *view, cluster_id); - } -} - -} // namespace internal -} // namespace ceres - -#endif // CERES_NO_SUITESPARSE diff --git a/extern/libmv/third_party/ceres/internal/ceres/canonical_views_clustering.h b/extern/libmv/third_party/ceres/internal/ceres/canonical_views_clustering.h deleted file mode 100644 index 6b0b38a854d..00000000000 --- a/extern/libmv/third_party/ceres/internal/ceres/canonical_views_clustering.h +++ /dev/null @@ -1,136 +0,0 @@ -// Ceres Solver - A fast non-linear least squares minimizer -// Copyright 2015 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) -// -// An implementation of the Canonical Views clustering algorithm from -// "Scene Summarization for Online Image Collections", Ian Simon, Noah -// Snavely, Steven M. Seitz, ICCV 2007. -// -// More details can be found at -// http://grail.cs.washington.edu/projects/canonview/ -// -// Ceres uses this algorithm to perform view clustering for -// constructing visibility based preconditioners. - -#ifndef CERES_INTERNAL_CANONICAL_VIEWS_CLUSTERING_H_ -#define CERES_INTERNAL_CANONICAL_VIEWS_CLUSTERING_H_ - -// This include must come before any #ifndef check on Ceres compile options. -#include "ceres/internal/port.h" - -#ifndef CERES_NO_SUITESPARSE - -#include <vector> - -#include "ceres/collections_port.h" -#include "ceres/graph.h" - -namespace ceres { -namespace internal { - -struct CanonicalViewsClusteringOptions; - -// Compute a partitioning of the vertices of the graph using the -// canonical views clustering algorithm. -// -// In the following we will use the terms vertices and views -// interchangably. Given a weighted Graph G(V,E), the canonical views -// of G are the the set of vertices that best "summarize" the content -// of the graph. If w_ij i s the weight connecting the vertex i to -// vertex j, and C is the set of canonical views. Then the objective -// of the canonical views algorithm is -// -// E[C] = sum_[i in V] max_[j in C] w_ij -// - size_penalty_weight * |C| -// - similarity_penalty_weight * sum_[i in C, j in C, j > i] w_ij -// -// alpha is the size penalty that penalizes large number of canonical -// views. -// -// beta is the similarity penalty that penalizes canonical views that -// are too similar to other canonical views. -// -// Thus the canonical views algorithm tries to find a canonical view -// for each vertex in the graph which best explains it, while trying -// to minimize the number of canonical views and the overlap between -// them. -// -// We further augment the above objective function by allowing for per -// vertex weights, higher weights indicating a higher preference for -// being chosen as a canonical view. Thus if w_i is the vertex weight -// for vertex i, the objective function is then -// -// E[C] = sum_[i in V] max_[j in C] w_ij -// - size_penalty_weight * |C| -// - similarity_penalty_weight * sum_[i in C, j in C, j > i] w_ij -// + view_score_weight * sum_[i in C] w_i -// -// centers will contain the vertices that are the identified -// as the canonical views/cluster centers, and membership is a map -// from vertices to cluster_ids. The i^th cluster center corresponds -// to the i^th cluster. -// -// It is possible depending on the configuration of the clustering -// algorithm that some of the vertices may not be assigned to any -// cluster. In this case they are assigned to a cluster with id = -1; -void ComputeCanonicalViewsClustering( - const CanonicalViewsClusteringOptions& options, - const WeightedGraph<int>& graph, - std::vector<int>* centers, - HashMap<int, int>* membership); - -struct CanonicalViewsClusteringOptions { - CanonicalViewsClusteringOptions() - : min_views(3), - size_penalty_weight(5.75), - similarity_penalty_weight(100.0), - view_score_weight(0.0) { - } - // The minimum number of canonical views to compute. - int min_views; - - // Penalty weight for the number of canonical views. A higher - // number will result in fewer canonical views. - double size_penalty_weight; - - // Penalty weight for the diversity (orthogonality) of the - // canonical views. A higher number will encourage less similar - // canonical views. - double similarity_penalty_weight; - - // Weight for per-view scores. Lower weight places less - // confidence in the view scores. - double view_score_weight; -}; - -} // namespace internal -} // namespace ceres - -#endif // CERES_NO_SUITESPARSE -#endif // CERES_INTERNAL_CANONICAL_VIEWS_CLUSTERING_H_ diff --git a/extern/libmv/third_party/ceres/internal/ceres/cxsparse.cc b/extern/libmv/third_party/ceres/internal/ceres/cxsparse.cc deleted file mode 100644 index 60b58089bb5..00000000000 --- a/extern/libmv/third_party/ceres/internal/ceres/cxsparse.cc +++ /dev/null @@ -1,218 +0,0 @@ -// Ceres Solver - A fast non-linear least squares minimizer -// Copyright 2015 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: strandmark@google.com (Petter Strandmark) - -// This include must come before any #ifndef check on Ceres compile options. -#include "ceres/internal/port.h" - -#ifndef CERES_NO_CXSPARSE - -#include "ceres/cxsparse.h" - -#include <vector> -#include "ceres/compressed_col_sparse_matrix_utils.h" -#include "ceres/compressed_row_sparse_matrix.h" -#include "ceres/triplet_sparse_matrix.h" -#include "glog/logging.h" - -namespace ceres { -namespace internal { - -using std::vector; - -CXSparse::CXSparse() : scratch_(NULL), scratch_size_(0) { -} - -CXSparse::~CXSparse() { - if (scratch_size_ > 0) { - cs_di_free(scratch_); - } -} - - -bool CXSparse::SolveCholesky(cs_di* A, - cs_dis* symbolic_factorization, - double* b) { - // Make sure we have enough scratch space available. - if (scratch_size_ < A->n) { - if (scratch_size_ > 0) { - cs_di_free(scratch_); - } - scratch_ = - reinterpret_cast<CS_ENTRY*>(cs_di_malloc(A->n, sizeof(CS_ENTRY))); - scratch_size_ = A->n; - } - - // Solve using Cholesky factorization - csn* numeric_factorization = cs_di_chol(A, symbolic_factorization); - if (numeric_factorization == NULL) { - LOG(WARNING) << "Cholesky factorization failed."; - return false; - } - - // When the Cholesky factorization succeeded, these methods are - // guaranteed to succeeded as well. In the comments below, "x" - // refers to the scratch space. - // - // Set x = P * b. - cs_di_ipvec(symbolic_factorization->pinv, b, scratch_, A->n); - // Set x = L \ x. - cs_di_lsolve(numeric_factorization->L, scratch_); - // Set x = L' \ x. - cs_di_ltsolve(numeric_factorization->L, scratch_); - // Set b = P' * x. - cs_di_pvec(symbolic_factorization->pinv, scratch_, b, A->n); - - // Free Cholesky factorization. - cs_di_nfree(numeric_factorization); - return true; -} - -cs_dis* CXSparse::AnalyzeCholesky(cs_di* A) { - // order = 1 for Cholesky factorization. - return cs_schol(1, A); -} - -cs_dis* CXSparse::AnalyzeCholeskyWithNaturalOrdering(cs_di* A) { - // order = 0 for Natural ordering. - return cs_schol(0, A); -} - -cs_dis* CXSparse::BlockAnalyzeCholesky(cs_di* A, - const vector<int>& row_blocks, - const vector<int>& col_blocks) { - const int num_row_blocks = row_blocks.size(); - const int num_col_blocks = col_blocks.size(); - - vector<int> block_rows; - vector<int> block_cols; - CompressedColumnScalarMatrixToBlockMatrix(A->i, - A->p, - row_blocks, - col_blocks, - &block_rows, - &block_cols); - cs_di block_matrix; - block_matrix.m = num_row_blocks; - block_matrix.n = num_col_blocks; - block_matrix.nz = -1; - block_matrix.nzmax = block_rows.size(); - block_matrix.p = &block_cols[0]; - block_matrix.i = &block_rows[0]; - block_matrix.x = NULL; - - int* ordering = cs_amd(1, &block_matrix); - vector<int> block_ordering(num_row_blocks, -1); - std::copy(ordering, ordering + num_row_blocks, &block_ordering[0]); - cs_free(ordering); - - vector<int> scalar_ordering; - BlockOrderingToScalarOrdering(row_blocks, block_ordering, &scalar_ordering); - - cs_dis* symbolic_factorization = - reinterpret_cast<cs_dis*>(cs_calloc(1, sizeof(cs_dis))); - symbolic_factorization->pinv = cs_pinv(&scalar_ordering[0], A->n); - cs* permuted_A = cs_symperm(A, symbolic_factorization->pinv, 0); - - symbolic_factorization->parent = cs_etree(permuted_A, 0); - int* postordering = cs_post(symbolic_factorization->parent, A->n); - int* column_counts = cs_counts(permuted_A, - symbolic_factorization->parent, - postordering, - 0); - cs_free(postordering); - cs_spfree(permuted_A); - - symbolic_factorization->cp = (int*) cs_malloc(A->n+1, sizeof(int)); - symbolic_factorization->lnz = cs_cumsum(symbolic_factorization->cp, - column_counts, - A->n); - symbolic_factorization->unz = symbolic_factorization->lnz; - - cs_free(column_counts); - - if (symbolic_factorization->lnz < 0) { - cs_sfree(symbolic_factorization); - symbolic_factorization = NULL; - } - - return symbolic_factorization; -} - -cs_di CXSparse::CreateSparseMatrixTransposeView(CompressedRowSparseMatrix* A) { - cs_di At; - At.m = A->num_cols(); - At.n = A->num_rows(); - At.nz = -1; - At.nzmax = A->num_nonzeros(); - At.p = A->mutable_rows(); - At.i = A->mutable_cols(); - At.x = A->mutable_values(); - return At; -} - -cs_di* CXSparse::CreateSparseMatrix(TripletSparseMatrix* tsm) { - cs_di_sparse tsm_wrapper; - tsm_wrapper.nzmax = tsm->num_nonzeros(); - tsm_wrapper.nz = tsm->num_nonzeros(); - tsm_wrapper.m = tsm->num_rows(); - tsm_wrapper.n = tsm->num_cols(); - tsm_wrapper.p = tsm->mutable_cols(); - tsm_wrapper.i = tsm->mutable_rows(); - tsm_wrapper.x = tsm->mutable_values(); - - return cs_compress(&tsm_wrapper); -} - -void CXSparse::ApproximateMinimumDegreeOrdering(cs_di* A, int* ordering) { - int* cs_ordering = cs_amd(1, A); - std::copy(cs_ordering, cs_ordering + A->m, ordering); - cs_free(cs_ordering); -} - -cs_di* CXSparse::TransposeMatrix(cs_di* A) { - return cs_di_transpose(A, 1); -} - -cs_di* CXSparse::MatrixMatrixMultiply(cs_di* A, cs_di* B) { - return cs_di_multiply(A, B); -} - -void CXSparse::Free(cs_di* sparse_matrix) { - cs_di_spfree(sparse_matrix); -} - -void CXSparse::Free(cs_dis* symbolic_factorization) { - cs_di_sfree(symbolic_factorization); -} - -} // namespace internal -} // namespace ceres - -#endif // CERES_NO_CXSPARSE diff --git a/extern/libmv/third_party/ceres/internal/ceres/single_linkage_clustering.cc b/extern/libmv/third_party/ceres/internal/ceres/single_linkage_clustering.cc deleted file mode 100644 index 490fcd87cad..00000000000 --- a/extern/libmv/third_party/ceres/internal/ceres/single_linkage_clustering.cc +++ /dev/null @@ -1,110 +0,0 @@ -// Ceres Solver - A fast non-linear least squares minimizer -// Copyright 2015 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) - -// This include must come before any #ifndef check on Ceres compile options. -#include "ceres/internal/port.h" - -#ifndef CERES_NO_SUITESPARSE - -#include "ceres/single_linkage_clustering.h" - -#include "ceres/graph.h" -#include "ceres/collections_port.h" -#include "ceres/graph_algorithms.h" - -namespace ceres { -namespace internal { - -int ComputeSingleLinkageClustering( - const SingleLinkageClusteringOptions& options, - const WeightedGraph<int>& graph, - HashMap<int, int>* membership) { - CHECK_NOTNULL(membership)->clear(); - - // Initially each vertex is in its own cluster. - const HashSet<int>& vertices = graph.vertices(); - for (HashSet<int>::const_iterator it = vertices.begin(); - it != vertices.end(); - ++it) { - (*membership)[*it] = *it; - } - - for (HashSet<int>::const_iterator it1 = vertices.begin(); - it1 != vertices.end(); - ++it1) { - const int vertex1 = *it1; - const HashSet<int>& neighbors = graph.Neighbors(vertex1); - for (HashSet<int>::const_iterator it2 = neighbors.begin(); - it2 != neighbors.end(); - ++it2) { - const int vertex2 = *it2; - - // Since the graph is undirected, only pay attention to one side - // of the edge and ignore weak edges. - if ((vertex1 > vertex2) || - (graph.EdgeWeight(vertex1, vertex2) < options.min_similarity)) { - continue; - } - - // Use a union-find algorithm to keep track of the clusters. - const int c1 = FindConnectedComponent(vertex1, membership); - const int c2 = FindConnectedComponent(vertex2, membership); - - if (c1 == c2) { - continue; - } - - if (c1 < c2) { - (*membership)[c2] = c1; - } else { - (*membership)[c1] = c2; - } - } - } - - // Make sure that every vertex is connected directly to the vertex - // identifying the cluster. - int num_clusters = 0; - for (HashMap<int, int>::iterator it = membership->begin(); - it != membership->end(); - ++it) { - it->second = FindConnectedComponent(it->first, membership); - if (it->first == it->second) { - ++num_clusters; - } - } - - return num_clusters; -} - -} // namespace internal -} // namespace ceres - -#endif // CERES_NO_SUITESPARSE diff --git a/extern/libmv/third_party/ceres/internal/ceres/single_linkage_clustering.h b/extern/libmv/third_party/ceres/internal/ceres/single_linkage_clustering.h deleted file mode 100644 index fb02f01de45..00000000000 --- a/extern/libmv/third_party/ceres/internal/ceres/single_linkage_clustering.h +++ /dev/null @@ -1,74 +0,0 @@ -// Ceres Solver - A fast non-linear least squares minimizer -// Copyright 2015 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) - -#ifndef CERES_INTERNAL_SINGLE_LINKAGE_CLUSTERING_H_ -#define CERES_INTERNAL_SINGLE_LINKAGE_CLUSTERING_H_ - -// This include must come before any #ifndef check on Ceres compile options. -#include "ceres/internal/port.h" - -#ifndef CERES_NO_SUITESPARSE - -#include "ceres/collections_port.h" -#include "ceres/graph.h" - -namespace ceres { -namespace internal { - -struct SingleLinkageClusteringOptions { - SingleLinkageClusteringOptions() - : min_similarity(0.99) { - } - - // Graph edges with edge weight less than min_similarity are ignored - // during the clustering process. - double min_similarity; -}; - -// Compute a partitioning of the vertices of the graph using the -// single linkage clustering algorithm. Edges with weight less than -// SingleLinkageClusteringOptions::min_similarity will be ignored. -// -// membership upon return will contain a mapping from the vertices of -// the graph to an integer indicating the identity of the cluster that -// it belongs to. -// -// The return value of this function is the number of clusters -// identified by the algorithm. -int ComputeSingleLinkageClustering( - const SingleLinkageClusteringOptions& options, - const WeightedGraph<int>& graph, - HashMap<int, int>* membership); - -} // namespace internal -} // namespace ceres - -#endif // CERES_NO_SUITESPARSE -#endif // CERES_INTERNAL_SINGLE_LINKAGE_CLUSTERING_H_ diff --git a/extern/libmv/third_party/ceres/internal/ceres/suitesparse.cc b/extern/libmv/third_party/ceres/internal/ceres/suitesparse.cc deleted file mode 100644 index 200daa2d8bf..00000000000 --- a/extern/libmv/third_party/ceres/internal/ceres/suitesparse.cc +++ /dev/null @@ -1,350 +0,0 @@ -// Ceres Solver - A fast non-linear least squares minimizer -// Copyright 2015 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) - -// This include must come before any #ifndef check on Ceres compile options. -#include "ceres/internal/port.h" - -#ifndef CERES_NO_SUITESPARSE -#include "ceres/suitesparse.h" - -#include <vector> -#include "cholmod.h" -#include "ceres/compressed_col_sparse_matrix_utils.h" -#include "ceres/compressed_row_sparse_matrix.h" -#include "ceres/linear_solver.h" -#include "ceres/triplet_sparse_matrix.h" - -namespace ceres { -namespace internal { - -using std::string; -using std::vector; - -SuiteSparse::SuiteSparse() { - cholmod_start(&cc_); -} - -SuiteSparse::~SuiteSparse() { - cholmod_finish(&cc_); -} - -cholmod_sparse* SuiteSparse::CreateSparseMatrix(TripletSparseMatrix* A) { - cholmod_triplet triplet; - - triplet.nrow = A->num_rows(); - triplet.ncol = A->num_cols(); - triplet.nzmax = A->max_num_nonzeros(); - triplet.nnz = A->num_nonzeros(); - triplet.i = reinterpret_cast<void*>(A->mutable_rows()); - triplet.j = reinterpret_cast<void*>(A->mutable_cols()); - triplet.x = reinterpret_cast<void*>(A->mutable_values()); - triplet.stype = 0; // Matrix is not symmetric. - triplet.itype = CHOLMOD_INT; - triplet.xtype = CHOLMOD_REAL; - triplet.dtype = CHOLMOD_DOUBLE; - - return cholmod_triplet_to_sparse(&triplet, triplet.nnz, &cc_); -} - - -cholmod_sparse* SuiteSparse::CreateSparseMatrixTranspose( - TripletSparseMatrix* A) { - cholmod_triplet triplet; - - triplet.ncol = A->num_rows(); // swap row and columns - triplet.nrow = A->num_cols(); - triplet.nzmax = A->max_num_nonzeros(); - triplet.nnz = A->num_nonzeros(); - - // swap rows and columns - triplet.j = reinterpret_cast<void*>(A->mutable_rows()); - triplet.i = reinterpret_cast<void*>(A->mutable_cols()); - triplet.x = reinterpret_cast<void*>(A->mutable_values()); - triplet.stype = 0; // Matrix is not symmetric. - triplet.itype = CHOLMOD_INT; - triplet.xtype = CHOLMOD_REAL; - triplet.dtype = CHOLMOD_DOUBLE; - - return cholmod_triplet_to_sparse(&triplet, triplet.nnz, &cc_); -} - -cholmod_sparse SuiteSparse::CreateSparseMatrixTransposeView( - CompressedRowSparseMatrix* A) { - cholmod_sparse m; - m.nrow = A->num_cols(); - m.ncol = A->num_rows(); - m.nzmax = A->num_nonzeros(); - m.nz = NULL; - m.p = reinterpret_cast<void*>(A->mutable_rows()); - m.i = reinterpret_cast<void*>(A->mutable_cols()); - m.x = reinterpret_cast<void*>(A->mutable_values()); - m.z = NULL; - m.stype = 0; // Matrix is not symmetric. - m.itype = CHOLMOD_INT; - m.xtype = CHOLMOD_REAL; - m.dtype = CHOLMOD_DOUBLE; - m.sorted = 1; - m.packed = 1; - - return m; -} - -cholmod_dense* SuiteSparse::CreateDenseVector(const double* x, - int in_size, - int out_size) { - CHECK_LE(in_size, out_size); - cholmod_dense* v = cholmod_zeros(out_size, 1, CHOLMOD_REAL, &cc_); - if (x != NULL) { - memcpy(v->x, x, in_size*sizeof(*x)); - } - return v; -} - -cholmod_factor* SuiteSparse::AnalyzeCholesky(cholmod_sparse* A, - string* message) { - // Cholmod can try multiple re-ordering strategies to find a fill - // reducing ordering. Here we just tell it use AMD with automatic - // matrix dependence choice of supernodal versus simplicial - // factorization. - cc_.nmethods = 1; - cc_.method[0].ordering = CHOLMOD_AMD; - cc_.supernodal = CHOLMOD_AUTO; - - cholmod_factor* factor = cholmod_analyze(A, &cc_); - if (VLOG_IS_ON(2)) { - cholmod_print_common(const_cast<char*>("Symbolic Analysis"), &cc_); - } - - if (cc_.status != CHOLMOD_OK) { - *message = StringPrintf("cholmod_analyze failed. error code: %d", - cc_.status); - return NULL; - } - - return CHECK_NOTNULL(factor); -} - -cholmod_factor* SuiteSparse::BlockAnalyzeCholesky( - cholmod_sparse* A, - const vector<int>& row_blocks, - const vector<int>& col_blocks, - string* message) { - vector<int> ordering; - if (!BlockAMDOrdering(A, row_blocks, col_blocks, &ordering)) { - return NULL; - } - return AnalyzeCholeskyWithUserOrdering(A, ordering, message); -} - -cholmod_factor* SuiteSparse::AnalyzeCholeskyWithUserOrdering( - cholmod_sparse* A, - const vector<int>& ordering, - string* message) { - CHECK_EQ(ordering.size(), A->nrow); - - cc_.nmethods = 1; - cc_.method[0].ordering = CHOLMOD_GIVEN; - - cholmod_factor* factor = - cholmod_analyze_p(A, const_cast<int*>(&ordering[0]), NULL, 0, &cc_); - if (VLOG_IS_ON(2)) { - cholmod_print_common(const_cast<char*>("Symbolic Analysis"), &cc_); - } - if (cc_.status != CHOLMOD_OK) { - *message = StringPrintf("cholmod_analyze failed. error code: %d", - cc_.status); - return NULL; - } - - return CHECK_NOTNULL(factor); -} - -cholmod_factor* SuiteSparse::AnalyzeCholeskyWithNaturalOrdering( - cholmod_sparse* A, - string* message) { - cc_.nmethods = 1; - cc_.method[0].ordering = CHOLMOD_NATURAL; - cc_.postorder = 0; - - cholmod_factor* factor = cholmod_analyze(A, &cc_); - if (VLOG_IS_ON(2)) { - cholmod_print_common(const_cast<char*>("Symbolic Analysis"), &cc_); - } - if (cc_.status != CHOLMOD_OK) { - *message = StringPrintf("cholmod_analyze failed. error code: %d", - cc_.status); - return NULL; - } - - return CHECK_NOTNULL(factor); -} - -bool SuiteSparse::BlockAMDOrdering(const cholmod_sparse* A, - const vector<int>& row_blocks, - const vector<int>& col_blocks, - vector<int>* ordering) { - const int num_row_blocks = row_blocks.size(); - const int num_col_blocks = col_blocks.size(); - - // Arrays storing the compressed column structure of the matrix - // incoding the block sparsity of A. - vector<int> block_cols; - vector<int> block_rows; - - CompressedColumnScalarMatrixToBlockMatrix(reinterpret_cast<const int*>(A->i), - reinterpret_cast<const int*>(A->p), - row_blocks, - col_blocks, - &block_rows, - &block_cols); - - cholmod_sparse_struct block_matrix; - block_matrix.nrow = num_row_blocks; - block_matrix.ncol = num_col_blocks; - block_matrix.nzmax = block_rows.size(); - block_matrix.p = reinterpret_cast<void*>(&block_cols[0]); - block_matrix.i = reinterpret_cast<void*>(&block_rows[0]); - block_matrix.x = NULL; - block_matrix.stype = A->stype; - block_matrix.itype = CHOLMOD_INT; - block_matrix.xtype = CHOLMOD_PATTERN; - block_matrix.dtype = CHOLMOD_DOUBLE; - block_matrix.sorted = 1; - block_matrix.packed = 1; - - vector<int> block_ordering(num_row_blocks); - if (!cholmod_amd(&block_matrix, NULL, 0, &block_ordering[0], &cc_)) { - return false; - } - - BlockOrderingToScalarOrdering(row_blocks, block_ordering, ordering); - return true; -} - -LinearSolverTerminationType SuiteSparse::Cholesky(cholmod_sparse* A, - cholmod_factor* L, - string* message) { - CHECK_NOTNULL(A); - CHECK_NOTNULL(L); - - // Save the current print level and silence CHOLMOD, otherwise - // CHOLMOD is prone to dumping stuff to stderr, which can be - // distracting when the error (matrix is indefinite) is not a fatal - // failure. - const int old_print_level = cc_.print; - cc_.print = 0; - - cc_.quick_return_if_not_posdef = 1; - int cholmod_status = cholmod_factorize(A, L, &cc_); - cc_.print = old_print_level; - - // TODO(sameeragarwal): This switch statement is not consistent. It - // treats all kinds of CHOLMOD failures as warnings. Some of these - // like out of memory are definitely not warnings. The problem is - // that the return value Cholesky is two valued, but the state of - // the linear solver is really three valued. SUCCESS, - // NON_FATAL_FAILURE (e.g., indefinite matrix) and FATAL_FAILURE - // (e.g. out of memory). - switch (cc_.status) { - case CHOLMOD_NOT_INSTALLED: - *message = "CHOLMOD failure: Method not installed."; - return LINEAR_SOLVER_FATAL_ERROR; - case CHOLMOD_OUT_OF_MEMORY: - *message = "CHOLMOD failure: Out of memory."; - return LINEAR_SOLVER_FATAL_ERROR; - case CHOLMOD_TOO_LARGE: - *message = "CHOLMOD failure: Integer overflow occured."; - return LINEAR_SOLVER_FATAL_ERROR; - case CHOLMOD_INVALID: - *message = "CHOLMOD failure: Invalid input."; - return LINEAR_SOLVER_FATAL_ERROR; - case CHOLMOD_NOT_POSDEF: - *message = "CHOLMOD warning: Matrix not positive definite."; - return LINEAR_SOLVER_FAILURE; - case CHOLMOD_DSMALL: - *message = "CHOLMOD warning: D for LDL' or diag(L) or " - "LL' has tiny absolute value."; - return LINEAR_SOLVER_FAILURE; - case CHOLMOD_OK: - if (cholmod_status != 0) { - return LINEAR_SOLVER_SUCCESS; - } - - *message = "CHOLMOD failure: cholmod_factorize returned false " - "but cholmod_common::status is CHOLMOD_OK." - "Please report this to ceres-solver@googlegroups.com."; - return LINEAR_SOLVER_FATAL_ERROR; - default: - *message = - StringPrintf("Unknown cholmod return code: %d. " - "Please report this to ceres-solver@googlegroups.com.", - cc_.status); - return LINEAR_SOLVER_FATAL_ERROR; - } - - return LINEAR_SOLVER_FATAL_ERROR; -} - -cholmod_dense* SuiteSparse::Solve(cholmod_factor* L, - cholmod_dense* b, - string* message) { - if (cc_.status != CHOLMOD_OK) { - *message = "cholmod_solve failed. CHOLMOD status is not CHOLMOD_OK"; - return NULL; - } - - return cholmod_solve(CHOLMOD_A, L, b, &cc_); -} - -bool SuiteSparse::ApproximateMinimumDegreeOrdering(cholmod_sparse* matrix, - int* ordering) { - return cholmod_amd(matrix, NULL, 0, ordering, &cc_); -} - -bool SuiteSparse::ConstrainedApproximateMinimumDegreeOrdering( - cholmod_sparse* matrix, - int* constraints, - int* ordering) { -#ifndef CERES_NO_CAMD - return cholmod_camd(matrix, NULL, 0, constraints, ordering, &cc_); -#else - LOG(FATAL) << "Congratulations you have found a bug in Ceres." - << "Ceres Solver was compiled with SuiteSparse " - << "version 4.1.0 or less. Calling this function " - << "in that case is a bug. Please contact the" - << "the Ceres Solver developers."; - return false; -#endif -} - -} // namespace internal -} // namespace ceres - -#endif // CERES_NO_SUITESPARSE diff --git a/extern/libmv/third_party/ceres/internal/ceres/visibility.cc b/extern/libmv/third_party/ceres/internal/ceres/visibility.cc deleted file mode 100644 index 55b5322aee6..00000000000 --- a/extern/libmv/third_party/ceres/internal/ceres/visibility.cc +++ /dev/null @@ -1,166 +0,0 @@ -// Ceres Solver - A fast non-linear least squares minimizer -// Copyright 2015 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: kushalav@google.com (Avanish Kushal) - -// This include must come before any #ifndef check on Ceres compile options. -#include "ceres/internal/port.h" - -#ifndef CERES_NO_SUITESPARSE - -#include "ceres/visibility.h" - -#include <cmath> -#include <ctime> -#include <algorithm> -#include <set> -#include <vector> -#include <utility> -#include "ceres/block_structure.h" -#include "ceres/collections_port.h" -#include "ceres/graph.h" -#include "glog/logging.h" - -namespace ceres { -namespace internal { - -using std::make_pair; -using std::max; -using std::pair; -using std::set; -using std::vector; - -void ComputeVisibility(const CompressedRowBlockStructure& block_structure, - const int num_eliminate_blocks, - vector< set<int> >* visibility) { - CHECK_NOTNULL(visibility); - - // Clear the visibility vector and resize it to hold a - // vector for each camera. - visibility->resize(0); - visibility->resize(block_structure.cols.size() - num_eliminate_blocks); - - for (int i = 0; i < block_structure.rows.size(); ++i) { - const vector<Cell>& cells = block_structure.rows[i].cells; - int block_id = cells[0].block_id; - // If the first block is not an e_block, then skip this row block. - if (block_id >= num_eliminate_blocks) { - continue; - } - - for (int j = 1; j < cells.size(); ++j) { - int camera_block_id = cells[j].block_id - num_eliminate_blocks; - DCHECK_GE(camera_block_id, 0); - DCHECK_LT(camera_block_id, visibility->size()); - (*visibility)[camera_block_id].insert(block_id); - } - } -} - -WeightedGraph<int>* CreateSchurComplementGraph( - const vector<set<int> >& visibility) { - const time_t start_time = time(NULL); - // Compute the number of e_blocks/point blocks. Since the visibility - // set for each e_block/camera contains the set of e_blocks/points - // visible to it, we find the maximum across all visibility sets. - int num_points = 0; - for (int i = 0; i < visibility.size(); i++) { - if (visibility[i].size() > 0) { - num_points = max(num_points, (*visibility[i].rbegin()) + 1); - } - } - - // Invert the visibility. The input is a camera->point mapping, - // which tells us which points are visible in which - // cameras. However, to compute the sparsity structure of the Schur - // Complement efficiently, its better to have the point->camera - // mapping. - vector<set<int> > inverse_visibility(num_points); - for (int i = 0; i < visibility.size(); i++) { - const set<int>& visibility_set = visibility[i]; - for (set<int>::const_iterator it = visibility_set.begin(); - it != visibility_set.end(); - ++it) { - inverse_visibility[*it].insert(i); - } - } - - // Map from camera pairs to number of points visible to both cameras - // in the pair. - HashMap<pair<int, int>, int > camera_pairs; - - // Count the number of points visible to each camera/f_block pair. - for (vector<set<int> >::const_iterator it = inverse_visibility.begin(); - it != inverse_visibility.end(); - ++it) { - const set<int>& inverse_visibility_set = *it; - for (set<int>::const_iterator camera1 = inverse_visibility_set.begin(); - camera1 != inverse_visibility_set.end(); - ++camera1) { - set<int>::const_iterator camera2 = camera1; - for (++camera2; camera2 != inverse_visibility_set.end(); ++camera2) { - ++(camera_pairs[make_pair(*camera1, *camera2)]); - } - } - } - - WeightedGraph<int>* graph = new WeightedGraph<int>; - - // Add vertices and initialize the pairs for self edges so that self - // edges are guaranteed. This is needed for the Canonical views - // algorithm to work correctly. - static const double kSelfEdgeWeight = 1.0; - for (int i = 0; i < visibility.size(); ++i) { - graph->AddVertex(i); - graph->AddEdge(i, i, kSelfEdgeWeight); - } - - // Add an edge for each camera pair. - for (HashMap<pair<int, int>, int>::const_iterator it = camera_pairs.begin(); - it != camera_pairs.end(); - ++it) { - const int camera1 = it->first.first; - const int camera2 = it->first.second; - CHECK_NE(camera1, camera2); - - const int count = it->second; - // Static cast necessary for Windows. - const double weight = static_cast<double>(count) / - (sqrt(static_cast<double>( - visibility[camera1].size() * visibility[camera2].size()))); - graph->AddEdge(camera1, camera2, weight); - } - - VLOG(2) << "Schur complement graph time: " << (time(NULL) - start_time); - return graph; -} - -} // namespace internal -} // namespace ceres - -#endif // CERES_NO_SUITESPARSE diff --git a/extern/libmv/third_party/ceres/internal/ceres/visibility.h b/extern/libmv/third_party/ceres/internal/ceres/visibility.h deleted file mode 100644 index 142eb2838a1..00000000000 --- a/extern/libmv/third_party/ceres/internal/ceres/visibility.h +++ /dev/null @@ -1,84 +0,0 @@ -// Ceres Solver - A fast non-linear least squares minimizer -// Copyright 2015 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: kushalav@google.com (Avanish Kushal) -// sameeragarwal@google.com (Sameer Agarwal) -// -// Functions to manipulate visibility information from the block -// structure of sparse matrices. - -#ifndef CERES_INTERNAL_VISIBILITY_H_ -#define CERES_INTERNAL_VISIBILITY_H_ - -// This include must come before any #ifndef check on Ceres compile options. -#include "ceres/internal/port.h" - -#ifndef CERES_NO_SUITESPARSE - -#include <set> -#include <vector> -#include "ceres/graph.h" - -namespace ceres { -namespace internal { - -struct CompressedRowBlockStructure; - -// Given a compressed row block structure, computes the set of -// e_blocks "visible" to each f_block. If an e_block co-occurs with an -// f_block in a residual block, it is visible to the f_block. The -// first num_eliminate_blocks columns blocks are e_blocks and the rest -// f_blocks. -// -// In a structure from motion problem, e_blocks correspond to 3D -// points and f_blocks correspond to cameras. -void ComputeVisibility(const CompressedRowBlockStructure& block_structure, - int num_eliminate_blocks, - std::vector<std::set<int> >* visibility); - -// Given f_block visibility as computed by the ComputeVisibility -// function above, construct and return a graph whose vertices are -// f_blocks and an edge connects two vertices if they have atleast one -// e_block in common. The weight of this edge is normalized dot -// product between the visibility vectors of the two -// vertices/f_blocks. -// -// This graph reflects the sparsity structure of reduced camera -// matrix/Schur complement matrix obtained by eliminating the e_blocks -// from the normal equations. -// -// Caller acquires ownership of the returned WeightedGraph pointer -// (heap-allocated). -WeightedGraph<int>* CreateSchurComplementGraph( - const std::vector<std::set<int> >& visibility); - -} // namespace internal -} // namespace ceres - -#endif // CERES_NO_SUITESPARSE -#endif // CERES_INTERNAL_VISIBILITY_H_ diff --git a/extern/libmv/third_party/ceres/internal/ceres/visibility_based_preconditioner.cc b/extern/libmv/third_party/ceres/internal/ceres/visibility_based_preconditioner.cc deleted file mode 100644 index b0000cdb1fd..00000000000 --- a/extern/libmv/third_party/ceres/internal/ceres/visibility_based_preconditioner.cc +++ /dev/null @@ -1,631 +0,0 @@ -// Ceres Solver - A fast non-linear least squares minimizer -// Copyright 2015 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) - -// This include must come before any #ifndef check on Ceres compile options. -#include "ceres/internal/port.h" - -#ifndef CERES_NO_SUITESPARSE - -#include "ceres/visibility_based_preconditioner.h" - -#include <algorithm> -#include <functional> -#include <iterator> -#include <set> -#include <utility> -#include <vector> -#include "Eigen/Dense" -#include "ceres/block_random_access_sparse_matrix.h" -#include "ceres/block_sparse_matrix.h" -#include "ceres/canonical_views_clustering.h" -#include "ceres/collections_port.h" -#include "ceres/graph.h" -#include "ceres/graph_algorithms.h" -#include "ceres/internal/scoped_ptr.h" -#include "ceres/linear_solver.h" -#include "ceres/schur_eliminator.h" -#include "ceres/single_linkage_clustering.h" -#include "ceres/visibility.h" -#include "glog/logging.h" - -namespace ceres { -namespace internal { - -using std::make_pair; -using std::pair; -using std::set; -using std::swap; -using std::vector; - -// TODO(sameeragarwal): Currently these are magic weights for the -// preconditioner construction. Move these higher up into the Options -// struct and provide some guidelines for choosing them. -// -// This will require some more work on the clustering algorithm and -// possibly some more refactoring of the code. -static const double kCanonicalViewsSizePenaltyWeight = 3.0; -static const double kCanonicalViewsSimilarityPenaltyWeight = 0.0; -static const double kSingleLinkageMinSimilarity = 0.9; - -VisibilityBasedPreconditioner::VisibilityBasedPreconditioner( - const CompressedRowBlockStructure& bs, - const Preconditioner::Options& options) - : options_(options), - num_blocks_(0), - num_clusters_(0), - factor_(NULL) { - CHECK_GT(options_.elimination_groups.size(), 1); - CHECK_GT(options_.elimination_groups[0], 0); - CHECK(options_.type == CLUSTER_JACOBI || - options_.type == CLUSTER_TRIDIAGONAL) - << "Unknown preconditioner type: " << options_.type; - num_blocks_ = bs.cols.size() - options_.elimination_groups[0]; - CHECK_GT(num_blocks_, 0) - << "Jacobian should have atleast 1 f_block for " - << "visibility based preconditioning."; - - // Vector of camera block sizes - block_size_.resize(num_blocks_); - for (int i = 0; i < num_blocks_; ++i) { - block_size_[i] = bs.cols[i + options_.elimination_groups[0]].size; - } - - const time_t start_time = time(NULL); - switch (options_.type) { - case CLUSTER_JACOBI: - ComputeClusterJacobiSparsity(bs); - break; - case CLUSTER_TRIDIAGONAL: - ComputeClusterTridiagonalSparsity(bs); - break; - default: - LOG(FATAL) << "Unknown preconditioner type"; - } - const time_t structure_time = time(NULL); - InitStorage(bs); - const time_t storage_time = time(NULL); - InitEliminator(bs); - const time_t eliminator_time = time(NULL); - - // Allocate temporary storage for a vector used during - // RightMultiply. - tmp_rhs_ = CHECK_NOTNULL(ss_.CreateDenseVector(NULL, - m_->num_rows(), - m_->num_rows())); - const time_t init_time = time(NULL); - VLOG(2) << "init time: " - << init_time - start_time - << " structure time: " << structure_time - start_time - << " storage time:" << storage_time - structure_time - << " eliminator time: " << eliminator_time - storage_time; -} - -VisibilityBasedPreconditioner::~VisibilityBasedPreconditioner() { - if (factor_ != NULL) { - ss_.Free(factor_); - factor_ = NULL; - } - if (tmp_rhs_ != NULL) { - ss_.Free(tmp_rhs_); - tmp_rhs_ = NULL; - } -} - -// Determine the sparsity structure of the CLUSTER_JACOBI -// preconditioner. It clusters cameras using their scene -// visibility. The clusters form the diagonal blocks of the -// preconditioner matrix. -void VisibilityBasedPreconditioner::ComputeClusterJacobiSparsity( - const CompressedRowBlockStructure& bs) { - vector<set<int> > visibility; - ComputeVisibility(bs, options_.elimination_groups[0], &visibility); - CHECK_EQ(num_blocks_, visibility.size()); - ClusterCameras(visibility); - cluster_pairs_.clear(); - for (int i = 0; i < num_clusters_; ++i) { - cluster_pairs_.insert(make_pair(i, i)); - } -} - -// Determine the sparsity structure of the CLUSTER_TRIDIAGONAL -// preconditioner. It clusters cameras using using the scene -// visibility and then finds the strongly interacting pairs of -// clusters by constructing another graph with the clusters as -// vertices and approximating it with a degree-2 maximum spanning -// forest. The set of edges in this forest are the cluster pairs. -void VisibilityBasedPreconditioner::ComputeClusterTridiagonalSparsity( - const CompressedRowBlockStructure& bs) { - vector<set<int> > visibility; - ComputeVisibility(bs, options_.elimination_groups[0], &visibility); - CHECK_EQ(num_blocks_, visibility.size()); - ClusterCameras(visibility); - - // Construct a weighted graph on the set of clusters, where the - // edges are the number of 3D points/e_blocks visible in both the - // clusters at the ends of the edge. Return an approximate degree-2 - // maximum spanning forest of this graph. - vector<set<int> > cluster_visibility; - ComputeClusterVisibility(visibility, &cluster_visibility); - scoped_ptr<WeightedGraph<int> > cluster_graph( - CHECK_NOTNULL(CreateClusterGraph(cluster_visibility))); - scoped_ptr<WeightedGraph<int> > forest( - CHECK_NOTNULL(Degree2MaximumSpanningForest(*cluster_graph))); - ForestToClusterPairs(*forest, &cluster_pairs_); -} - -// Allocate storage for the preconditioner matrix. -void VisibilityBasedPreconditioner::InitStorage( - const CompressedRowBlockStructure& bs) { - ComputeBlockPairsInPreconditioner(bs); - m_.reset(new BlockRandomAccessSparseMatrix(block_size_, block_pairs_)); -} - -// Call the canonical views algorithm and cluster the cameras based on -// their visibility sets. The visibility set of a camera is the set of -// e_blocks/3D points in the scene that are seen by it. -// -// The cluster_membership_ vector is updated to indicate cluster -// memberships for each camera block. -void VisibilityBasedPreconditioner::ClusterCameras( - const vector<set<int> >& visibility) { - scoped_ptr<WeightedGraph<int> > schur_complement_graph( - CHECK_NOTNULL(CreateSchurComplementGraph(visibility))); - - HashMap<int, int> membership; - - if (options_.visibility_clustering_type == CANONICAL_VIEWS) { - vector<int> centers; - CanonicalViewsClusteringOptions clustering_options; - clustering_options.size_penalty_weight = - kCanonicalViewsSizePenaltyWeight; - clustering_options.similarity_penalty_weight = - kCanonicalViewsSimilarityPenaltyWeight; - ComputeCanonicalViewsClustering(clustering_options, - *schur_complement_graph, - ¢ers, - &membership); - num_clusters_ = centers.size(); - } else if (options_.visibility_clustering_type == SINGLE_LINKAGE) { - SingleLinkageClusteringOptions clustering_options; - clustering_options.min_similarity = - kSingleLinkageMinSimilarity; - num_clusters_ = ComputeSingleLinkageClustering(clustering_options, - *schur_complement_graph, - &membership); - } else { - LOG(FATAL) << "Unknown visibility clustering algorithm."; - } - - CHECK_GT(num_clusters_, 0); - VLOG(2) << "num_clusters: " << num_clusters_; - FlattenMembershipMap(membership, &cluster_membership_); -} - -// Compute the block sparsity structure of the Schur complement -// matrix. For each pair of cameras contributing a non-zero cell to -// the schur complement, determine if that cell is present in the -// preconditioner or not. -// -// A pair of cameras contribute a cell to the preconditioner if they -// are part of the same cluster or if the the two clusters that they -// belong have an edge connecting them in the degree-2 maximum -// spanning forest. -// -// For example, a camera pair (i,j) where i belonges to cluster1 and -// j belongs to cluster2 (assume that cluster1 < cluster2). -// -// The cell corresponding to (i,j) is present in the preconditioner -// if cluster1 == cluster2 or the pair (cluster1, cluster2) were -// connected by an edge in the degree-2 maximum spanning forest. -// -// Since we have already expanded the forest into a set of camera -// pairs/edges, including self edges, the check can be reduced to -// checking membership of (cluster1, cluster2) in cluster_pairs_. -void VisibilityBasedPreconditioner::ComputeBlockPairsInPreconditioner( - const CompressedRowBlockStructure& bs) { - block_pairs_.clear(); - for (int i = 0; i < num_blocks_; ++i) { - block_pairs_.insert(make_pair(i, i)); - } - - int r = 0; - const int num_row_blocks = bs.rows.size(); - const int num_eliminate_blocks = options_.elimination_groups[0]; - - // Iterate over each row of the matrix. The block structure of the - // matrix is assumed to be sorted in order of the e_blocks/point - // blocks. Thus all row blocks containing an e_block/point occur - // contiguously. Further, if present, an e_block is always the first - // parameter block in each row block. These structural assumptions - // are common to all Schur complement based solvers in Ceres. - // - // For each e_block/point block we identify the set of cameras - // seeing it. The cross product of this set with itself is the set - // of non-zero cells contibuted by this e_block. - // - // The time complexity of this is O(nm^2) where, n is the number of - // 3d points and m is the maximum number of cameras seeing any - // point, which for most scenes is a fairly small number. - while (r < num_row_blocks) { - int e_block_id = bs.rows[r].cells.front().block_id; - if (e_block_id >= num_eliminate_blocks) { - // Skip the rows whose first block is an f_block. - break; - } - - set<int> f_blocks; - for (; r < num_row_blocks; ++r) { - const CompressedRow& row = bs.rows[r]; - if (row.cells.front().block_id != e_block_id) { - break; - } - - // Iterate over the blocks in the row, ignoring the first block - // since it is the one to be eliminated and adding the rest to - // the list of f_blocks associated with this e_block. - for (int c = 1; c < row.cells.size(); ++c) { - const Cell& cell = row.cells[c]; - const int f_block_id = cell.block_id - num_eliminate_blocks; - CHECK_GE(f_block_id, 0); - f_blocks.insert(f_block_id); - } - } - - for (set<int>::const_iterator block1 = f_blocks.begin(); - block1 != f_blocks.end(); - ++block1) { - set<int>::const_iterator block2 = block1; - ++block2; - for (; block2 != f_blocks.end(); ++block2) { - if (IsBlockPairInPreconditioner(*block1, *block2)) { - block_pairs_.insert(make_pair(*block1, *block2)); - } - } - } - } - - // The remaining rows which do not contain any e_blocks. - for (; r < num_row_blocks; ++r) { - const CompressedRow& row = bs.rows[r]; - CHECK_GE(row.cells.front().block_id, num_eliminate_blocks); - for (int i = 0; i < row.cells.size(); ++i) { - const int block1 = row.cells[i].block_id - num_eliminate_blocks; - for (int j = 0; j < row.cells.size(); ++j) { - const int block2 = row.cells[j].block_id - num_eliminate_blocks; - if (block1 <= block2) { - if (IsBlockPairInPreconditioner(block1, block2)) { - block_pairs_.insert(make_pair(block1, block2)); - } - } - } - } - } - - VLOG(1) << "Block pair stats: " << block_pairs_.size(); -} - -// Initialize the SchurEliminator. -void VisibilityBasedPreconditioner::InitEliminator( - const CompressedRowBlockStructure& bs) { - LinearSolver::Options eliminator_options; - eliminator_options.elimination_groups = options_.elimination_groups; - eliminator_options.num_threads = options_.num_threads; - eliminator_options.e_block_size = options_.e_block_size; - eliminator_options.f_block_size = options_.f_block_size; - eliminator_options.row_block_size = options_.row_block_size; - eliminator_.reset(SchurEliminatorBase::Create(eliminator_options)); - eliminator_->Init(eliminator_options.elimination_groups[0], &bs); -} - -// Update the values of the preconditioner matrix and factorize it. -bool VisibilityBasedPreconditioner::UpdateImpl(const BlockSparseMatrix& A, - const double* D) { - const time_t start_time = time(NULL); - const int num_rows = m_->num_rows(); - CHECK_GT(num_rows, 0); - - // We need a dummy rhs vector and a dummy b vector since the Schur - // eliminator combines the computation of the reduced camera matrix - // with the computation of the right hand side of that linear - // system. - // - // TODO(sameeragarwal): Perhaps its worth refactoring the - // SchurEliminator::Eliminate function to allow NULL for the rhs. As - // of now it does not seem to be worth the effort. - Vector rhs = Vector::Zero(m_->num_rows()); - Vector b = Vector::Zero(A.num_rows()); - - // Compute a subset of the entries of the Schur complement. - eliminator_->Eliminate(&A, b.data(), D, m_.get(), rhs.data()); - - // Try factorizing the matrix. For CLUSTER_JACOBI, this should - // always succeed modulo some numerical/conditioning problems. For - // CLUSTER_TRIDIAGONAL, in general the preconditioner matrix as - // constructed is not positive definite. However, we will go ahead - // and try factorizing it. If it works, great, otherwise we scale - // all the cells in the preconditioner corresponding to the edges in - // the degree-2 forest and that guarantees positive - // definiteness. The proof of this fact can be found in Lemma 1 in - // "Visibility Based Preconditioning for Bundle Adjustment". - // - // Doing the factorization like this saves us matrix mass when - // scaling is not needed, which is quite often in our experience. - LinearSolverTerminationType status = Factorize(); - - if (status == LINEAR_SOLVER_FATAL_ERROR) { - return false; - } - - // The scaling only affects the tri-diagonal case, since - // ScaleOffDiagonalBlocks only pays attenion to the cells that - // belong to the edges of the degree-2 forest. In the CLUSTER_JACOBI - // case, the preconditioner is guaranteed to be positive - // semidefinite. - if (status == LINEAR_SOLVER_FAILURE && options_.type == CLUSTER_TRIDIAGONAL) { - VLOG(1) << "Unscaled factorization failed. Retrying with off-diagonal " - << "scaling"; - ScaleOffDiagonalCells(); - status = Factorize(); - } - - VLOG(2) << "Compute time: " << time(NULL) - start_time; - return (status == LINEAR_SOLVER_SUCCESS); -} - -// Consider the preconditioner matrix as meta-block matrix, whose -// blocks correspond to the clusters. Then cluster pairs corresponding -// to edges in the degree-2 forest are off diagonal entries of this -// matrix. Scaling these off-diagonal entries by 1/2 forces this -// matrix to be positive definite. -void VisibilityBasedPreconditioner::ScaleOffDiagonalCells() { - for (set<pair<int, int> >::const_iterator it = block_pairs_.begin(); - it != block_pairs_.end(); - ++it) { - const int block1 = it->first; - const int block2 = it->second; - if (!IsBlockPairOffDiagonal(block1, block2)) { - continue; - } - - int r, c, row_stride, col_stride; - CellInfo* cell_info = m_->GetCell(block1, block2, - &r, &c, - &row_stride, &col_stride); - CHECK(cell_info != NULL) - << "Cell missing for block pair (" << block1 << "," << block2 << ")" - << " cluster pair (" << cluster_membership_[block1] - << " " << cluster_membership_[block2] << ")"; - - // Ah the magic of tri-diagonal matrices and diagonal - // dominance. See Lemma 1 in "Visibility Based Preconditioning - // For Bundle Adjustment". - MatrixRef m(cell_info->values, row_stride, col_stride); - m.block(r, c, block_size_[block1], block_size_[block2]) *= 0.5; - } -} - -// Compute the sparse Cholesky factorization of the preconditioner -// matrix. -LinearSolverTerminationType VisibilityBasedPreconditioner::Factorize() { - // Extract the TripletSparseMatrix that is used for actually storing - // S and convert it into a cholmod_sparse object. - cholmod_sparse* lhs = ss_.CreateSparseMatrix( - down_cast<BlockRandomAccessSparseMatrix*>( - m_.get())->mutable_matrix()); - - // The matrix is symmetric, and the upper triangular part of the - // matrix contains the values. - lhs->stype = 1; - - // TODO(sameeragarwal): Refactor to pipe this up and out. - std::string status; - - // Symbolic factorization is computed if we don't already have one handy. - if (factor_ == NULL) { - factor_ = ss_.BlockAnalyzeCholesky(lhs, block_size_, block_size_, &status); - } - - const LinearSolverTerminationType termination_type = - (factor_ != NULL) - ? ss_.Cholesky(lhs, factor_, &status) - : LINEAR_SOLVER_FATAL_ERROR; - - ss_.Free(lhs); - return termination_type; -} - -void VisibilityBasedPreconditioner::RightMultiply(const double* x, - double* y) const { - CHECK_NOTNULL(x); - CHECK_NOTNULL(y); - SuiteSparse* ss = const_cast<SuiteSparse*>(&ss_); - - const int num_rows = m_->num_rows(); - memcpy(CHECK_NOTNULL(tmp_rhs_)->x, x, m_->num_rows() * sizeof(*x)); - // TODO(sameeragarwal): Better error handling. - std::string status; - cholmod_dense* solution = - CHECK_NOTNULL(ss->Solve(factor_, tmp_rhs_, &status)); - memcpy(y, solution->x, sizeof(*y) * num_rows); - ss->Free(solution); -} - -int VisibilityBasedPreconditioner::num_rows() const { - return m_->num_rows(); -} - -// Classify camera/f_block pairs as in and out of the preconditioner, -// based on whether the cluster pair that they belong to is in the -// preconditioner or not. -bool VisibilityBasedPreconditioner::IsBlockPairInPreconditioner( - const int block1, - const int block2) const { - int cluster1 = cluster_membership_[block1]; - int cluster2 = cluster_membership_[block2]; - if (cluster1 > cluster2) { - swap(cluster1, cluster2); - } - return (cluster_pairs_.count(make_pair(cluster1, cluster2)) > 0); -} - -bool VisibilityBasedPreconditioner::IsBlockPairOffDiagonal( - const int block1, - const int block2) const { - return (cluster_membership_[block1] != cluster_membership_[block2]); -} - -// Convert a graph into a list of edges that includes self edges for -// each vertex. -void VisibilityBasedPreconditioner::ForestToClusterPairs( - const WeightedGraph<int>& forest, - HashSet<pair<int, int> >* cluster_pairs) const { - CHECK_NOTNULL(cluster_pairs)->clear(); - const HashSet<int>& vertices = forest.vertices(); - CHECK_EQ(vertices.size(), num_clusters_); - - // Add all the cluster pairs corresponding to the edges in the - // forest. - for (HashSet<int>::const_iterator it1 = vertices.begin(); - it1 != vertices.end(); - ++it1) { - const int cluster1 = *it1; - cluster_pairs->insert(make_pair(cluster1, cluster1)); - const HashSet<int>& neighbors = forest.Neighbors(cluster1); - for (HashSet<int>::const_iterator it2 = neighbors.begin(); - it2 != neighbors.end(); - ++it2) { - const int cluster2 = *it2; - if (cluster1 < cluster2) { - cluster_pairs->insert(make_pair(cluster1, cluster2)); - } - } - } -} - -// The visibilty set of a cluster is the union of the visibilty sets -// of all its cameras. In other words, the set of points visible to -// any camera in the cluster. -void VisibilityBasedPreconditioner::ComputeClusterVisibility( - const vector<set<int> >& visibility, - vector<set<int> >* cluster_visibility) const { - CHECK_NOTNULL(cluster_visibility)->resize(0); - cluster_visibility->resize(num_clusters_); - for (int i = 0; i < num_blocks_; ++i) { - const int cluster_id = cluster_membership_[i]; - (*cluster_visibility)[cluster_id].insert(visibility[i].begin(), - visibility[i].end()); - } -} - -// Construct a graph whose vertices are the clusters, and the edge -// weights are the number of 3D points visible to cameras in both the -// vertices. -WeightedGraph<int>* VisibilityBasedPreconditioner::CreateClusterGraph( - const vector<set<int> >& cluster_visibility) const { - WeightedGraph<int>* cluster_graph = new WeightedGraph<int>; - - for (int i = 0; i < num_clusters_; ++i) { - cluster_graph->AddVertex(i); - } - - for (int i = 0; i < num_clusters_; ++i) { - const set<int>& cluster_i = cluster_visibility[i]; - for (int j = i+1; j < num_clusters_; ++j) { - vector<int> intersection; - const set<int>& cluster_j = cluster_visibility[j]; - set_intersection(cluster_i.begin(), cluster_i.end(), - cluster_j.begin(), cluster_j.end(), - back_inserter(intersection)); - - if (intersection.size() > 0) { - // Clusters interact strongly when they share a large number - // of 3D points. The degree-2 maximum spanning forest - // alorithm, iterates on the edges in decreasing order of - // their weight, which is the number of points shared by the - // two cameras that it connects. - cluster_graph->AddEdge(i, j, intersection.size()); - } - } - } - return cluster_graph; -} - -// Canonical views clustering returns a HashMap from vertices to -// cluster ids. Convert this into a flat array for quick lookup. It is -// possible that some of the vertices may not be associated with any -// cluster. In that case, randomly assign them to one of the clusters. -// -// The cluster ids can be non-contiguous integers. So as we flatten -// the membership_map, we also map the cluster ids to a contiguous set -// of integers so that the cluster ids are in [0, num_clusters_). -void VisibilityBasedPreconditioner::FlattenMembershipMap( - const HashMap<int, int>& membership_map, - vector<int>* membership_vector) const { - CHECK_NOTNULL(membership_vector)->resize(0); - membership_vector->resize(num_blocks_, -1); - - HashMap<int, int> cluster_id_to_index; - // Iterate over the cluster membership map and update the - // cluster_membership_ vector assigning arbitrary cluster ids to - // the few cameras that have not been clustered. - for (HashMap<int, int>::const_iterator it = membership_map.begin(); - it != membership_map.end(); - ++it) { - const int camera_id = it->first; - int cluster_id = it->second; - - // If the view was not clustered, randomly assign it to one of the - // clusters. This preserves the mathematical correctness of the - // preconditioner. If there are too many views which are not - // clustered, it may lead to some quality degradation though. - // - // TODO(sameeragarwal): Check if a large number of views have not - // been clustered and deal with it? - if (cluster_id == -1) { - cluster_id = camera_id % num_clusters_; - } - - const int index = FindWithDefault(cluster_id_to_index, - cluster_id, - cluster_id_to_index.size()); - - if (index == cluster_id_to_index.size()) { - cluster_id_to_index[cluster_id] = index; - } - - CHECK_LT(index, num_clusters_); - membership_vector->at(camera_id) = index; - } -} - -} // namespace internal -} // namespace ceres - -#endif // CERES_NO_SUITESPARSE |