diff --git a/docs/source/solving.rst b/docs/source/solving.rst index f9d4fb36f..8202d1664 100644 --- a/docs/source/solving.rst +++ b/docs/source/solving.rst @@ -1147,6 +1147,28 @@ elimination group [LiSaad]_. ``CLUSTER_JACOBI`` and ``CLUSTER_TRIDIAGONAL``. See :ref:`section-preconditioner` for more details. +.. member:: VisibilityClusteringType Solver::Options::visibility_clustering_type + + Default: ``CANONICAL_VIEWS`` + + Type of clustering algorithm to use when constructing a visibility + based preconditioner. The original visibility based preconditioning + paper and implementation only used the canonical views algorithm. + + This algorithm gives high quality results but for large dense + graphs can be particularly expensive. As its worst case complexity + is cubic in size of the graph. + + Another option is to use ``SINGLE_LINKAGE`` which is a simple + thresholded single linkage clustering algorithm that only pays + attention to tightly coupled blocks in the Schur complement. This + is a fast algorithm that works well. + + The optimal choice of the clustering algorithm depends on the + sparsity structure of the problem, but generally speaking we + recommend that you try ``CANONICAL_VIEWS`` first and if it is too + expensive try ``SINGLE_LINKAGE``. + .. member:: DenseLinearAlgebraLibrary Solver::Options::dense_linear_algebra_library_type Default:``EIGEN`` @@ -2020,6 +2042,13 @@ The three arrays will be: Type of preconditioner used for solving the trust region step. Only meaningful when an iterative linear solver is used. +.. member:: VisibilityClusteringType Solver::Summary::visibility_clustering_type + + Type of clustering algorithm used for visibility based + preconditioning. Only meaningful when the + :member:`Solver::Summary::preconditioner_type` is + ``CLUSTER_JACOBI`` or ``CLUSTER_TRIDIAGONAL``. + .. member:: TrustRegionStrategyType Solver::Summary::trust_region_strategy_type Type of trust region strategy. diff --git a/examples/bundle_adjuster.cc b/examples/bundle_adjuster.cc index 224ad748e..fb827194b 100644 --- a/examples/bundle_adjuster.cc +++ b/examples/bundle_adjuster.cc @@ -82,6 +82,9 @@ DEFINE_string(linear_solver, "sparse_schur", "Options are: " DEFINE_string(preconditioner, "jacobi", "Options are: " "identity, jacobi, schur_jacobi, cluster_jacobi, " "cluster_tridiagonal."); +DEFINE_string(visibility_clustering, "canonical_views", + "single_linkage, canonical_views"); + DEFINE_string(sparse_linear_algebra_library, "suite_sparse", "Options are: suite_sparse and cx_sparse."); DEFINE_string(dense_linear_algebra_library, "eigen", @@ -125,6 +128,8 @@ void SetLinearSolver(Solver::Options* options) { &options->linear_solver_type)); CHECK(StringToPreconditionerType(FLAGS_preconditioner, &options->preconditioner_type)); + CHECK(StringToVisibilityClusteringType(FLAGS_visibility_clustering, + &options->visibility_clustering_type)); CHECK(StringToSparseLinearAlgebraLibraryType( FLAGS_sparse_linear_algebra_library, &options->sparse_linear_algebra_library_type)); diff --git a/include/ceres/solver.h b/include/ceres/solver.h index 0b3a18e5b..7d8e4e28b 100644 --- a/include/ceres/solver.h +++ b/include/ceres/solver.h @@ -98,7 +98,7 @@ class Solver { #endif preconditioner_type = JACOBI; - + visibility_clustering_type = CANONICAL_VIEWS; dense_linear_algebra_library_type = EIGEN; sparse_linear_algebra_library_type = SUITE_SPARSE; #if defined(CERES_NO_SUITESPARSE) && !defined(CERES_NO_CXSPARSE) @@ -385,6 +385,11 @@ class Solver { // Type of preconditioner to use with the iterative linear solvers. PreconditionerType preconditioner_type; + // Type of clustering algorithm to use visibility based + // preconditioning. This option is used only when the + // preconditioner_type is CLUSTER_JACOBI or CLUSTER_TRIDIAGONAL. + VisibilityClusteringType visibility_clustering_type; + // Ceres supports using multiple dense linear algebra libraries // for dense matrix factorizations. Currently EIGEN and LAPACK are // the valid choices. EIGEN is always available, LAPACK refers to @@ -886,6 +891,11 @@ class Solver { // step. Only meaningful when an iterative linear solver is used. PreconditionerType preconditioner_type; + // Type of clustering algorithm used for visibility based + // preconditioning. Only meaningful when the preconditioner_type + // is CLUSTER_JACOBI or CLUSTER_TRIDIAGONAL. + VisibilityClusteringType visibility_clustering_type; + // Type of trust region strategy. TrustRegionStrategyType trust_region_strategy_type; diff --git a/include/ceres/types.h b/include/ceres/types.h index ffa743a2d..c8c9cfed9 100644 --- a/include/ceres/types.h +++ b/include/ceres/types.h @@ -102,21 +102,49 @@ enum PreconditionerType { // Block diagonal of the Gauss-Newton Hessian. JACOBI, + // Note: The following three preconditioners can only be used with + // the ITERATIVE_SCHUR solver. They are well suited for Structure + // from Motion problems. + // Block diagonal of the Schur complement. This preconditioner may // only be used with the ITERATIVE_SCHUR solver. SCHUR_JACOBI, // Visibility clustering based preconditioners. // - // These preconditioners are well suited for Structure from Motion - // problems, particularly problems arising from community photo - // collections. These preconditioners use the visibility structure - // of the scene to determine the sparsity structure of the - // preconditioner. Requires SuiteSparse/CHOLMOD. + // The following two preconditioners use the visibility structure of + // the scene to determine the sparsity structure of the + // preconditioner. This is done using a clustering algorithm. The + // available visibility clustering algorithms are described below. + // + // Note: Requires SuiteSparse. CLUSTER_JACOBI, CLUSTER_TRIDIAGONAL }; +enum VisibilityClusteringType { + // Canonical views algorithm as described in + // + // "Scene Summarization for Online Image Collections", Ian Simon, Noah + // Snavely, Steven M. Seitz, ICCV 2007. + // + // This clustering algorithm can be quite slow, but gives high + // quality clusters. The original visibility based clustering paper + // used this algorithm. + CANONICAL_VIEWS, + + // The classic single linkage algorithm. It is extremely fast as + // compared to CANONICAL_VIEWS, but can give slightly poor + // results. For problems with large number of cameras though, this + // is generally a pretty good option. + // + // If you are using SCHUR_JACOBI preconditioner and have SuiteSparse + // available, CLUSTER_JACOBI and CLUSTER_TRIDIAGONAL in combination + // with the SINGLE_LINKAGE algorithm will generally give better + // results. + SINGLE_LINKAGE +}; + enum SparseLinearAlgebraLibraryType { // High performance sparse Cholesky factorization and approximate // minimum degree ordering. @@ -400,6 +428,10 @@ bool StringToLinearSolverType(string value, LinearSolverType* type); const char* PreconditionerTypeToString(PreconditionerType type); bool StringToPreconditionerType(string value, PreconditionerType* type); +const char* VisibilityClusteringTypeToString(VisibilityClusteringType type); +bool StringToVisibilityClusteringType(string value, + VisibilityClusteringType* type); + const char* SparseLinearAlgebraLibraryTypeToString( SparseLinearAlgebraLibraryType type); bool StringToSparseLinearAlgebraLibraryType( diff --git a/internal/ceres/CMakeLists.txt b/internal/ceres/CMakeLists.txt index d963fd64a..f4a457819 100644 --- a/internal/ceres/CMakeLists.txt +++ b/internal/ceres/CMakeLists.txt @@ -91,6 +91,7 @@ SET(CERES_INTERNAL_SRC schur_eliminator.cc schur_jacobi_preconditioner.cc scratch_evaluate_preparer.cc + single_linkage_clustering.cc solver.cc solver_impl.cc sparse_matrix.cc @@ -247,6 +248,7 @@ IF (BUILD_TESTING AND GFLAGS) CERES_TEST(rotation) CERES_TEST(schur_complement_solver) CERES_TEST(schur_eliminator) + CERES_TEST(single_linkage_clustering) CERES_TEST(small_blas) CERES_TEST(solver_impl) diff --git a/internal/ceres/canonical_views_clustering.cc b/internal/ceres/canonical_views_clustering.cc index 653194571..044d4381b 100644 --- a/internal/ceres/canonical_views_clustering.cc +++ b/internal/ceres/canonical_views_clustering.cc @@ -57,8 +57,8 @@ class CanonicalViewsClustering { // 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 Graph& graph, - const CanonicalViewsClusteringOptions& options, + void ComputeClustering(const CanonicalViewsClusteringOptions& options, + const Graph& graph, vector* centers, IntMap* membership); @@ -81,21 +81,21 @@ class CanonicalViewsClustering { }; void ComputeCanonicalViewsClustering( - const Graph& graph, const CanonicalViewsClusteringOptions& options, + const Graph& graph, vector* centers, IntMap* membership) { time_t start_time = time(NULL); CanonicalViewsClustering cv; - cv.ComputeClustering(graph, options, centers, membership); + cv.ComputeClustering(options, graph, centers, membership); VLOG(2) << "Canonical views clustering time (secs): " << time(NULL) - start_time; } // Implementation of CanonicalViewsClustering void CanonicalViewsClustering::ComputeClustering( - const Graph& graph, const CanonicalViewsClusteringOptions& options, + const Graph& graph, vector* centers, IntMap* membership) { options_ = options; diff --git a/internal/ceres/canonical_views_clustering.h b/internal/ceres/canonical_views_clustering.h index 48d1ed210..06d80c89e 100644 --- a/internal/ceres/canonical_views_clustering.h +++ b/internal/ceres/canonical_views_clustering.h @@ -47,9 +47,6 @@ #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 { @@ -100,8 +97,8 @@ struct CanonicalViewsClusteringOptions; // 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 Graph& graph, const CanonicalViewsClusteringOptions& options, + const Graph& graph, vector* centers, HashMap* membership); diff --git a/internal/ceres/canonical_views_clustering_test.cc b/internal/ceres/canonical_views_clustering_test.cc index 78d563563..f4f9e2deb 100644 --- a/internal/ceres/canonical_views_clustering_test.cc +++ b/internal/ceres/canonical_views_clustering_test.cc @@ -69,7 +69,7 @@ class CanonicalViewsTest : public ::testing::Test { } void ComputeClustering() { - ComputeCanonicalViewsClustering(graph_, options_, ¢ers_, &membership_); + ComputeCanonicalViewsClustering(options_, graph_, ¢ers_, &membership_); } Graph graph_; diff --git a/internal/ceres/iterative_schur_complement_solver.cc b/internal/ceres/iterative_schur_complement_solver.cc index 90013ff16..e9c54c86c 100644 --- a/internal/ceres/iterative_schur_complement_solver.cc +++ b/internal/ceres/iterative_schur_complement_solver.cc @@ -110,6 +110,8 @@ LinearSolver::Summary IterativeSchurComplementSolver::SolveImpl( Preconditioner::Options preconditioner_options; preconditioner_options.type = options_.preconditioner_type; + preconditioner_options.visibility_clustering_type = + options_.visibility_clustering_type; preconditioner_options.sparse_linear_algebra_library_type = options_.sparse_linear_algebra_library_type; preconditioner_options.num_threads = options_.num_threads; diff --git a/internal/ceres/linear_solver.h b/internal/ceres/linear_solver.h index 22691b33e..b0ab80dd7 100644 --- a/internal/ceres/linear_solver.h +++ b/internal/ceres/linear_solver.h @@ -74,6 +74,7 @@ class LinearSolver { Options() : type(SPARSE_NORMAL_CHOLESKY), preconditioner_type(JACOBI), + visibility_clustering_type(CANONICAL_VIEWS), dense_linear_algebra_library_type(EIGEN), sparse_linear_algebra_library_type(SUITE_SPARSE), use_postordering(false), @@ -89,6 +90,7 @@ class LinearSolver { LinearSolverType type; PreconditionerType preconditioner_type; + VisibilityClusteringType visibility_clustering_type; DenseLinearAlgebraLibraryType dense_linear_algebra_library_type; SparseLinearAlgebraLibraryType sparse_linear_algebra_library_type; diff --git a/internal/ceres/preconditioner.h b/internal/ceres/preconditioner.h index af64e3c9a..21cbc00b5 100644 --- a/internal/ceres/preconditioner.h +++ b/internal/ceres/preconditioner.h @@ -48,6 +48,7 @@ class Preconditioner : public LinearOperator { struct Options { Options() : type(JACOBI), + visibility_clustering_type(CANONICAL_VIEWS), sparse_linear_algebra_library_type(SUITE_SPARSE), num_threads(1), row_block_size(Eigen::Dynamic), @@ -56,7 +57,7 @@ class Preconditioner : public LinearOperator { } PreconditionerType type; - + VisibilityClusteringType visibility_clustering_type; SparseLinearAlgebraLibraryType sparse_linear_algebra_library_type; // If possible, how many threads the preconditioner can use. diff --git a/internal/ceres/single_linkage_clustering.cc b/internal/ceres/single_linkage_clustering.cc new file mode 100644 index 000000000..54b4379d5 --- /dev/null +++ b/internal/ceres/single_linkage_clustering.cc @@ -0,0 +1,107 @@ +// Ceres Solver - A fast non-linear least squares minimizer +// Copyright 2013 Google Inc. All rights reserved. +// http://code.google.com/p/ceres-solver/ +// +// 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_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 Graph& graph, + HashMap* membership) { + CHECK_NOTNULL(membership)->clear(); + + // Initially each vertex is in its own cluster. + const HashSet& vertices = graph.vertices(); + for (HashSet::const_iterator it = vertices.begin(); + it != vertices.end(); + ++it) { + (*membership)[*it] = *it; + } + + for (HashSet::const_iterator it1 = vertices.begin(); + it1 != vertices.end(); + ++it1) { + const int vertex1 = *it1; + const HashSet& neighbors = graph.Neighbors(vertex1); + for (HashSet::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::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/internal/ceres/single_linkage_clustering.h b/internal/ceres/single_linkage_clustering.h new file mode 100644 index 000000000..a0fe8f96c --- /dev/null +++ b/internal/ceres/single_linkage_clustering.h @@ -0,0 +1,71 @@ +// Ceres Solver - A fast non-linear least squares minimizer +// Copyright 2013 Google Inc. All rights reserved. +// http://code.google.com/p/ceres-solver/ +// +// 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_ + +#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 that 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 Graph& graph, + HashMap* membership); + +} // namespace internal +} // namespace ceres + +#endif // CERES_NO_SUITESPARSE +#endif // CERES_INTERNAL_SINGLE_LINKAGE_CLUSTERING_H_ diff --git a/internal/ceres/single_linkage_clustering_test.cc b/internal/ceres/single_linkage_clustering_test.cc new file mode 100644 index 000000000..8fbf84be2 --- /dev/null +++ b/internal/ceres/single_linkage_clustering_test.cc @@ -0,0 +1,125 @@ +// Ceres Solver - A fast non-linear least squares minimizer +// Copyright 2013 Google Inc. All rights reserved. +// http://code.google.com/p/ceres-solver/ +// +// 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: Sameer Agarwal (sameeragarwal@google.com) + +#ifndef CERES_NO_SUITESPARSE + +#include "ceres/single_linkage_clustering.h" + +#include "ceres/collections_port.h" +#include "ceres/graph.h" +#include "gtest/gtest.h" + +namespace ceres { +namespace internal { + +TEST(SingleLinkageClustering, GraphHasTwoComponents) { + Graph graph; + const int kNumVertices = 6; + for (int i = 0; i < kNumVertices; ++i) { + graph.AddVertex(i); + } + // Graph structure: + // + // 0-1-2-3 4-5 + graph.AddEdge(0, 1, 1.0); + graph.AddEdge(1, 2, 1.0); + graph.AddEdge(2, 3, 1.0); + graph.AddEdge(4, 5, 1.0); + + SingleLinkageClusteringOptions options; + HashMap membership; + ComputeSingleLinkageClustering(options, graph, &membership); + EXPECT_EQ(membership.size(), kNumVertices); + + EXPECT_EQ(membership[1], membership[0]); + EXPECT_EQ(membership[2], membership[0]); + EXPECT_EQ(membership[3], membership[0]); + EXPECT_EQ(membership[4], membership[5]); +} + +TEST(SingleLinkageClustering, ComponentWithWeakLink) { + Graph graph; + const int kNumVertices = 6; + for (int i = 0; i < kNumVertices; ++i) { + graph.AddVertex(i); + } + // Graph structure: + // + // 0-1-2-3 4-5 + graph.AddEdge(0, 1, 1.0); + graph.AddEdge(1, 2, 1.0); + graph.AddEdge(2, 3, 1.0); + + // This component should break up into two. + graph.AddEdge(4, 5, 0.5); + + SingleLinkageClusteringOptions options; + HashMap membership; + ComputeSingleLinkageClustering(options, graph, &membership); + EXPECT_EQ(membership.size(), kNumVertices); + + EXPECT_EQ(membership[1], membership[0]); + EXPECT_EQ(membership[2], membership[0]); + EXPECT_EQ(membership[3], membership[0]); + EXPECT_NE(membership[4], membership[5]); +} + +TEST(SingleLinkageClustering, ComponentWithWeakLinkAndStrongLink) { + Graph graph; + const int kNumVertices = 6; + for (int i = 0; i < kNumVertices; ++i) { + graph.AddVertex(i); + } + // Graph structure: + // + // 0-1-2-3 4-5 + graph.AddEdge(0, 1, 1.0); + graph.AddEdge(1, 2, 1.0); + graph.AddEdge(2, 3, 0.5); // Weak link + graph.AddEdge(0, 3, 1.0); + + // This component should break up into two. + graph.AddEdge(4, 5, 1.0); + + SingleLinkageClusteringOptions options; + HashMap membership; + ComputeSingleLinkageClustering(options, graph, &membership); + EXPECT_EQ(membership.size(), kNumVertices); + + EXPECT_EQ(membership[1], membership[0]); + EXPECT_EQ(membership[2], membership[0]); + EXPECT_EQ(membership[3], membership[0]); + EXPECT_EQ(membership[4], membership[5]); +} + +} // namespace internal +} // namespace ceres + +#endif // CERES_NO_SUITESPARSE diff --git a/internal/ceres/solver.cc b/internal/ceres/solver.cc index 3b6774604..2750e2300 100644 --- a/internal/ceres/solver.cc +++ b/internal/ceres/solver.cc @@ -119,6 +119,7 @@ Solver::Summary::Summary() inner_iterations_given(false), inner_iterations_used(false), preconditioner_type(IDENTITY), + visibility_clustering_type(CANONICAL_VIEWS), trust_region_strategy_type(LEVENBERG_MARQUARDT), dense_linear_algebra_library_type(EIGEN), sparse_linear_algebra_library_type(SUITE_SPARSE), @@ -237,6 +238,12 @@ string Solver::Summary::FullReport() const { PreconditionerTypeToString(preconditioner_type)); } + if (preconditioner_type == CLUSTER_JACOBI || + preconditioner_type == CLUSTER_TRIDIAGONAL) { + StringAppendF(&report, "Visibility clustering%24s%25s\n", + VisibilityClusteringTypeToString(visibility_clustering_type), + VisibilityClusteringTypeToString(visibility_clustering_type)); + } StringAppendF(&report, "Threads % 25d% 25d\n", num_threads_given, num_threads_used); StringAppendF(&report, "Linear solver threads % 23d% 25d\n", diff --git a/internal/ceres/solver_impl.cc b/internal/ceres/solver_impl.cc index ae7fc3ffc..f544f8f5c 100644 --- a/internal/ceres/solver_impl.cc +++ b/internal/ceres/solver_impl.cc @@ -512,6 +512,7 @@ void SolverImpl::TrustRegionSolve(const Solver::Options& original_options, summary->linear_solver_type_used = options.linear_solver_type; summary->preconditioner_type = options.preconditioner_type; + summary->visibility_clustering_type = options.visibility_clustering_type; summary->num_linear_solver_threads_given = original_options.num_linear_solver_threads; @@ -1244,6 +1245,8 @@ LinearSolver* SolverImpl::CreateLinearSolver(Solver::Options* options, options->max_linear_solver_iterations; linear_solver_options.type = options->linear_solver_type; linear_solver_options.preconditioner_type = options->preconditioner_type; + linear_solver_options.visibility_clustering_type = + options->visibility_clustering_type; linear_solver_options.sparse_linear_algebra_library_type = options->sparse_linear_algebra_library_type; linear_solver_options.dense_linear_algebra_library_type = diff --git a/internal/ceres/types.cc b/internal/ceres/types.cc index a97f1a55e..80f8d39e8 100644 --- a/internal/ceres/types.cc +++ b/internal/ceres/types.cc @@ -278,6 +278,25 @@ bool StringToCovarianceAlgorithmType( return false; } +const char* VisibilityClusteringTypeToString( + VisibilityClusteringType type) { + switch (type) { + CASESTR(CANONICAL_VIEWS); + CASESTR(SINGLE_LINKAGE); + default: + return "UNKNOWN"; + } +} + +bool StringToVisibilityClusteringType( + string value, + VisibilityClusteringType* type) { + UpperCase(&value); + STRENUM(CANONICAL_VIEWS); + STRENUM(SINGLE_LINKAGE); + return false; +} + const char* SolverTerminationTypeToString(SolverTerminationType type) { switch (type) { CASESTR(NO_CONVERGENCE); diff --git a/internal/ceres/visibility_based_preconditioner.cc b/internal/ceres/visibility_based_preconditioner.cc index 6535b4252..c760c5b51 100644 --- a/internal/ceres/visibility_based_preconditioner.cc +++ b/internal/ceres/visibility_based_preconditioner.cc @@ -48,6 +48,7 @@ #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" @@ -60,8 +61,9 @@ namespace internal { // // This will require some more work on the clustering algorithm and // possibly some more refactoring of the code. -static const double kSizePenaltyWeight = 3.0; -static const double kSimilarityPenaltyWeight = 0.0; +static const double kCanonicalViewsSizePenaltyWeight = 3.0; +static const double kCanonicalViewsSimilarityPenaltyWeight = 0.0; +static const double kSingleLinkageMinSimilarity = 0.9; VisibilityBasedPreconditioner::VisibilityBasedPreconditioner( const CompressedRowBlockStructure& bs, @@ -187,17 +189,29 @@ void VisibilityBasedPreconditioner::ClusterCameras( scoped_ptr > schur_complement_graph( CHECK_NOTNULL(CreateSchurComplementGraph(visibility))); - CanonicalViewsClusteringOptions options; - options.size_penalty_weight = kSizePenaltyWeight; - options.similarity_penalty_weight = kSimilarityPenaltyWeight; - - vector centers; HashMap membership; - ComputeCanonicalViewsClustering(*schur_complement_graph, - options, - ¢ers, - &membership); - num_clusters_ = centers.size(); + + if (options_.visibility_clustering_type == CANONICAL_VIEWS) { + vector 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); + } + CHECK_GT(num_clusters_, 0); VLOG(2) << "num_clusters: " << num_clusters_; FlattenMembershipMap(membership, &cluster_membership_); @@ -542,11 +556,17 @@ Graph* VisibilityBasedPreconditioner::CreateClusterGraph( // 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& membership_map, vector* membership_vector) const { CHECK_NOTNULL(membership_vector)->resize(0); membership_vector->resize(num_blocks_, -1); + + HashMap 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. @@ -567,7 +587,17 @@ void VisibilityBasedPreconditioner::FlattenMembershipMap( cluster_id = camera_id % num_clusters_; } - membership_vector->at(camera_id) = cluster_id; + 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_NE(index, cluster_id_to_index.size()); + } + + CHECK_NE(index, num_clusters_); + membership_vector->at(camera_id) = index; } }