Refactor system_test 1. Move common test infrastructure into test_util. 2. system_test now only contains powells function. 3. Add bundle_adjustment_test. Instead of a single function which computes everything, there is now a test for each solver configuration which uses the reference solution computed by the fixture. Change-Id: I16a9a9a83a845a7aaf28762bcecf1a8ff5aee805
diff --git a/internal/ceres/CMakeLists.txt b/internal/ceres/CMakeLists.txt index 55643de..6c71c29 100644 --- a/internal/ceres/CMakeLists.txt +++ b/internal/ceres/CMakeLists.txt
@@ -285,6 +285,7 @@ ceres_test(block_random_access_diagonal_matrix) ceres_test(block_random_access_sparse_matrix) ceres_test(block_sparse_matrix) + ceres_test(bundle_adjustment) ceres_test(c_api) ceres_test(canonical_views_clustering) ceres_test(compressed_row_sparse_matrix)
diff --git a/internal/ceres/bundle_adjustment_test.cc b/internal/ceres/bundle_adjustment_test.cc new file mode 100644 index 0000000..c2f8843 --- /dev/null +++ b/internal/ceres/bundle_adjustment_test.cc
@@ -0,0 +1,560 @@ +// 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: keir@google.com (Keir Mierle) +// sameeragarwal@google.com (Sameer Agarwal) +// +// End-to-end bundle adjustment tests for Ceres. It uses a bundle +// adjustment problem with 16 cameras and two thousand points. + +#include <cmath> +#include <cstdio> +#include <cstdlib> +#include <string> + +#include "ceres/internal/port.h" + +#include "ceres/autodiff_cost_function.h" +#include "ceres/ordered_groups.h" +#include "ceres/problem.h" +#include "ceres/rotation.h" +#include "ceres/solver.h" +#include "ceres/stringprintf.h" +#include "ceres/test_util.h" +#include "ceres/types.h" +#include "gflags/gflags.h" +#include "glog/logging.h" +#include "gtest/gtest.h" + +namespace ceres { +namespace internal { + +using std::string; +using std::vector; + +const bool kAutomaticOrdering = true; +const bool kUserOrdering = false; + +// This class implements the SystemTestProblem interface and provides +// access to a bundle adjustment problem. It is based on +// examples/bundle_adjustment_example.cc. Currently a small 16 camera +// problem is hard coded in the constructor. +class BundleAdjustmentProblem { + public: + BundleAdjustmentProblem() { + const string input_file = TestFileAbsolutePath("problem-16-22106-pre.txt"); + ReadData(input_file); + BuildProblem(); + } + + ~BundleAdjustmentProblem() { + delete []point_index_; + delete []camera_index_; + delete []observations_; + delete []parameters_; + } + + Problem* mutable_problem() { return &problem_; } + Solver::Options* mutable_solver_options() { return &options_; } + + int num_cameras() const { return num_cameras_; } + int num_points() const { return num_points_; } + int num_observations() const { return num_observations_; } + const int* point_index() const { return point_index_; } + const int* camera_index() const { return camera_index_; } + const double* observations() const { return observations_; } + double* mutable_cameras() { return parameters_; } + double* mutable_points() { return parameters_ + 9 * num_cameras_; } + + static double kResidualTolerance; + + private: + void ReadData(const string& filename) { + FILE * fptr = fopen(filename.c_str(), "r"); + + if (!fptr) { + LOG(FATAL) << "File Error: unable to open file " << filename; + } + + // This will die horribly on invalid files. Them's the breaks. + FscanfOrDie(fptr, "%d", &num_cameras_); + FscanfOrDie(fptr, "%d", &num_points_); + FscanfOrDie(fptr, "%d", &num_observations_); + + VLOG(1) << "Header: " << num_cameras_ + << " " << num_points_ + << " " << num_observations_; + + point_index_ = new int[num_observations_]; + camera_index_ = new int[num_observations_]; + observations_ = new double[2 * num_observations_]; + + num_parameters_ = 9 * num_cameras_ + 3 * num_points_; + parameters_ = new double[num_parameters_]; + + for (int i = 0; i < num_observations_; ++i) { + FscanfOrDie(fptr, "%d", camera_index_ + i); + FscanfOrDie(fptr, "%d", point_index_ + i); + for (int j = 0; j < 2; ++j) { + FscanfOrDie(fptr, "%lf", observations_ + 2*i + j); + } + } + + for (int i = 0; i < num_parameters_; ++i) { + FscanfOrDie(fptr, "%lf", parameters_ + i); + } + } + + void BuildProblem() { + double* points = mutable_points(); + double* cameras = mutable_cameras(); + + for (int i = 0; i < num_observations(); ++i) { + // Each Residual block takes a point and a camera as input and + // outputs a 2 dimensional residual. + CostFunction* cost_function = + new AutoDiffCostFunction<BundlerResidual, 2, 9, 3>( + new BundlerResidual(observations_[2*i + 0], + observations_[2*i + 1])); + + // Each observation correponds to a pair of a camera and a point + // which are identified by camera_index()[i] and + // point_index()[i] respectively. + double* camera = cameras + 9 * camera_index_[i]; + double* point = points + 3 * point_index()[i]; + problem_.AddResidualBlock(cost_function, NULL, camera, point); + } + + options_.linear_solver_ordering.reset(new ParameterBlockOrdering); + + // The points come before the cameras. + for (int i = 0; i < num_points_; ++i) { + options_.linear_solver_ordering->AddElementToGroup(points + 3 * i, 0); + } + + for (int i = 0; i < num_cameras_; ++i) { + options_.linear_solver_ordering->AddElementToGroup(cameras + 9 * i, 1); + } + + options_.linear_solver_type = DENSE_SCHUR; + options_.max_num_iterations = 25; + options_.function_tolerance = 1e-10; + options_.gradient_tolerance = 1e-10; + options_.parameter_tolerance = 1e-10; + } + + template<typename T> + void FscanfOrDie(FILE *fptr, const char *format, T *value) { + int num_scanned = fscanf(fptr, format, value); + if (num_scanned != 1) { + LOG(FATAL) << "Invalid UW data file."; + } + } + + // Templated pinhole camera model. The camera is parameterized + // using 9 parameters. 3 for rotation, 3 for translation, 1 for + // focal length and 2 for radial distortion. The principal point is + // not modeled (i.e. it is assumed be located at the image center). + struct BundlerResidual { + // (u, v): the position of the observation with respect to the image + // center point. + BundlerResidual(double u, double v): u(u), v(v) {} + + template <typename T> + bool operator()(const T* const camera, + const T* const point, + T* residuals) const { + T p[3]; + AngleAxisRotatePoint(camera, point, p); + + // Add the translation vector + p[0] += camera[3]; + p[1] += camera[4]; + p[2] += camera[5]; + + const T& focal = camera[6]; + const T& l1 = camera[7]; + const T& l2 = camera[8]; + + // Compute the center of distortion. The sign change comes from + // the camera model that Noah Snavely's Bundler assumes, whereby + // the camera coordinate system has a negative z axis. + T xp = - focal * p[0] / p[2]; + T yp = - focal * p[1] / p[2]; + + // Apply second and fourth order radial distortion. + T r2 = xp*xp + yp*yp; + T distortion = T(1.0) + r2 * (l1 + l2 * r2); + + residuals[0] = distortion * xp - T(u); + residuals[1] = distortion * yp - T(v); + + return true; + } + + double u; + double v; + }; + + Problem problem_; + Solver::Options options_; + + int num_cameras_; + int num_points_; + int num_observations_; + int num_parameters_; + + int* point_index_; + int* camera_index_; + double* observations_; + // The parameter vector is laid out as follows + // [camera_1, ..., camera_n, point_1, ..., point_m] + double* parameters_; +}; + +double BundleAdjustmentProblem::kResidualTolerance = 1e-4; +typedef SystemTest<BundleAdjustmentProblem> BundleAdjustmentTest; + +TEST_F(BundleAdjustmentTest, DenseSchurWithAutomaticOrdering) { + RunSolverForConfigAndExpectResidualsMatch( + SolverConfig(DENSE_SCHUR, NO_SPARSE, kAutomaticOrdering)); +} + +TEST_F(BundleAdjustmentTest, DenseSchurWithUserOrdering) { + RunSolverForConfigAndExpectResidualsMatch( + SolverConfig(DENSE_SCHUR, NO_SPARSE, kUserOrdering)); +} + +TEST_F(BundleAdjustmentTest, IterativeSchurWithJacobiAndAutomaticOrdering) { + RunSolverForConfigAndExpectResidualsMatch( + SolverConfig(ITERATIVE_SCHUR, NO_SPARSE, kAutomaticOrdering, JACOBI)); +} + +TEST_F(BundleAdjustmentTest, IterativeSchurWithJacobiAndUserOrdering) { + RunSolverForConfigAndExpectResidualsMatch( + SolverConfig(ITERATIVE_SCHUR, NO_SPARSE, kUserOrdering, JACOBI)); +} + +TEST_F(BundleAdjustmentTest, + IterativeSchurWithSchurJacobiAndAutomaticOrdering) { + RunSolverForConfigAndExpectResidualsMatch( + SolverConfig(ITERATIVE_SCHUR, + NO_SPARSE, + kAutomaticOrdering, + SCHUR_JACOBI)); +} + +TEST_F(BundleAdjustmentTest, IterativeSchurWithSchurJacobiAndUserOrdering) { + RunSolverForConfigAndExpectResidualsMatch( + SolverConfig(ITERATIVE_SCHUR, NO_SPARSE, kUserOrdering, SCHUR_JACOBI)); +} + +#ifndef CERES_NO_SUITESPARSE +TEST_F(BundleAdjustmentTest, + SparseNormalCholeskyWithAutomaticOrderingUsingSuiteSparse) { + RunSolverForConfigAndExpectResidualsMatch( + SolverConfig(SPARSE_NORMAL_CHOLESKY, SUITE_SPARSE, kAutomaticOrdering)); +} + +TEST_F(BundleAdjustmentTest, + SparseNormalCholeskyWithUserOrderingUsingSuiteSparse) { + RunSolverForConfigAndExpectResidualsMatch( + SolverConfig(SPARSE_NORMAL_CHOLESKY, SUITE_SPARSE, kUserOrdering)); +} + +TEST_F(BundleAdjustmentTest, + SparseSchurWithAutomaticOrderingUsingSuiteSparse) { + RunSolverForConfigAndExpectResidualsMatch( + SolverConfig(SPARSE_SCHUR, SUITE_SPARSE, kAutomaticOrdering)); +} + +TEST_F(BundleAdjustmentTest, SparseSchurWithUserOrderingUsingSuiteSparse) { + RunSolverForConfigAndExpectResidualsMatch( + SolverConfig(SPARSE_SCHUR, SUITE_SPARSE, kUserOrdering)); +} + +TEST_F(BundleAdjustmentTest, + IterativeSchurWithClusterJacobiAndAutomaticOrderingUsingSuiteSparse) { + RunSolverForConfigAndExpectResidualsMatch( + SolverConfig(ITERATIVE_SCHUR, + SUITE_SPARSE, + kAutomaticOrdering, + CLUSTER_JACOBI)); +} + +TEST_F(BundleAdjustmentTest, + IterativeSchurWithClusterJacobiAndUserOrderingUsingSuiteSparse) { + RunSolverForConfigAndExpectResidualsMatch( + SolverConfig(ITERATIVE_SCHUR, + SUITE_SPARSE, + kUserOrdering, + CLUSTER_JACOBI)); +} + +TEST_F(BundleAdjustmentTest, + IterativeSchurWithClusterTridiagonalAndAutomaticOrderingUsingSuiteSparse) { + RunSolverForConfigAndExpectResidualsMatch( + SolverConfig(ITERATIVE_SCHUR, + SUITE_SPARSE, + kAutomaticOrdering, + CLUSTER_TRIDIAGONAL)); +} + +TEST_F(BundleAdjustmentTest, + IterativeSchurWithClusterTridiagonalAndUserOrderingUsingSuiteSparse) { + RunSolverForConfigAndExpectResidualsMatch( + SolverConfig(ITERATIVE_SCHUR, + SUITE_SPARSE, + kUserOrdering, + CLUSTER_TRIDIAGONAL)); +} +#endif // CERES_NO_SUITESPARSE + +#ifndef CERES_NO_CXSPARSE +TEST_F(BundleAdjustmentTest, + SparseNormalCholeskyWithAutomaticOrderingUsingCXSparse) { + RunSolverForConfigAndExpectResidualsMatch( + SolverConfig(SPARSE_NORMAL_CHOLESKY, CX_SPARSE, kAutomaticOrdering)); +} + +TEST_F(BundleAdjustmentTest, + SparseNormalCholeskyWithUserOrderingUsingCXSparse) { + RunSolverForConfigAndExpectResidualsMatch( + SolverConfig(SPARSE_NORMAL_CHOLESKY, CX_SPARSE, kUserOrdering)); +} + +TEST_F(BundleAdjustmentTest, SparseSchurWithAutomaticOrderingUsingCXSparse) { + RunSolverForConfigAndExpectResidualsMatch( + SolverConfig(SPARSE_SCHUR, CX_SPARSE, kAutomaticOrdering)); +} + +TEST_F(BundleAdjustmentTest, SparseSchurWithUserOrderingUsingCXSparse) { + RunSolverForConfigAndExpectResidualsMatch( + SolverConfig(SPARSE_SCHUR, CX_SPARSE, kUserOrdering)); +} +#endif // CERES_NO_CXSPARSE + +#ifdef CERES_USE_EIGEN_SPARSE +TEST_F(BundleAdjustmentTest, + SparseNormalCholeskyWithAutomaticOrderingUsingEigenSparse) { + RunSolverForConfigAndExpectResidualsMatch( + SolverConfig(SPARSE_NORMAL_CHOLESKY, EIGEN_SPARSE, kAutomaticOrdering)); +} + +TEST_F(BundleAdjustmentTest, + SparseNormalCholeskyWithUserOrderingUsingEigenSparse) { + RunSolverForConfigAndExpectResidualsMatch( + SolverConfig(SPARSE_NORMAL_CHOLESKY, EIGEN_SPARSE, kUserOrdering)); +} + +TEST_F(BundleAdjustmentTest, + SparseSchurWithAutomaticOrderingUsingEigenSparse) { + RunSolverForConfigAndExpectResidualsMatch( + SolverConfig(SPARSE_SCHUR, EIGEN_SPARSE, kAutomaticOrdering)); +} + +TEST_F(BundleAdjustmentTest, SparseSchurWithUserOrderingUsingEigenSparse) { + RunSolverForConfigAndExpectResidualsMatch( + SolverConfig(SPARSE_SCHUR, EIGEN_SPARSE, kUserOrdering)); +} +#endif // CERES_USE_EIGEN_SPARSE + +#ifdef CERES_USE_OPENMP + +TEST_F(BundleAdjustmentTest, MultiThreadedDenseSchurWithAutomaticOrdering) { + RunSolverForConfigAndExpectResidualsMatch( + ThreadedSolverConfig(DENSE_SCHUR, NO_SPARSE, kAutomaticOrdering)); +} + +TEST_F(BundleAdjustmentTest, MultiThreadedDenseSchurWithUserOrdering) { + RunSolverForConfigAndExpectResidualsMatch( + ThreadedSolverConfig(DENSE_SCHUR, NO_SPARSE, kUserOrdering)); +} + +TEST_F(BundleAdjustmentTest, + MultiThreadedIterativeSchurWithJacobiAndAutomaticOrdering) { + RunSolverForConfigAndExpectResidualsMatch( + ThreadedSolverConfig(ITERATIVE_SCHUR, + NO_SPARSE, + kAutomaticOrdering, + JACOBI)); +} + +TEST_F(BundleAdjustmentTest, + MultiThreadedIterativeSchurWithJacobiAndUserOrdering) { + RunSolverForConfigAndExpectResidualsMatch( + ThreadedSolverConfig(ITERATIVE_SCHUR, NO_SPARSE, kUserOrdering, JACOBI)); +} + +TEST_F(BundleAdjustmentTest, + MultiThreadedIterativeSchurWithSchurJacobiAndAutomaticOrdering) { + RunSolverForConfigAndExpectResidualsMatch( + ThreadedSolverConfig(ITERATIVE_SCHUR, + NO_SPARSE, + kAutomaticOrdering, + SCHUR_JACOBI)); +} + +TEST_F(BundleAdjustmentTest, + MultiThreadedIterativeSchurWithSchurJacobiAndUserOrdering) { + RunSolverForConfigAndExpectResidualsMatch( + ThreadedSolverConfig(ITERATIVE_SCHUR, + NO_SPARSE, + kUserOrdering, + SCHUR_JACOBI)); +} + +#ifndef CERES_NO_SUITESPARSE +TEST_F(BundleAdjustmentTest, + MultiThreadedSparseNormalCholeskyWithAutomaticOrderingUsingSuiteSparse) { + RunSolverForConfigAndExpectResidualsMatch( + ThreadedSolverConfig(SPARSE_NORMAL_CHOLESKY, + SUITE_SPARSE, + kAutomaticOrdering)); +} + +TEST_F(BundleAdjustmentTest, + MultiThreadedSparseNormalCholeskyWithUserOrderingUsingSuiteSparse) { + RunSolverForConfigAndExpectResidualsMatch( + ThreadedSolverConfig(SPARSE_NORMAL_CHOLESKY, + SUITE_SPARSE, + kUserOrdering)); +} + +TEST_F(BundleAdjustmentTest, + MultiThreadedSparseSchurWithAutomaticOrderingUsingSuiteSparse) { + RunSolverForConfigAndExpectResidualsMatch( + ThreadedSolverConfig(SPARSE_SCHUR, + SUITE_SPARSE, + kAutomaticOrdering)); +} + +TEST_F(BundleAdjustmentTest, + MultiThreadedSparseSchurWithUserOrderingUsingSuiteSparse) { + RunSolverForConfigAndExpectResidualsMatch( + ThreadedSolverConfig(SPARSE_SCHUR, SUITE_SPARSE, kUserOrdering)); +} + +TEST_F(BundleAdjustmentTest, + MultiThreadedIterativeSchurWithClusterJacobiAndAutomaticOrderingUsingSuiteSparse) { // NOLINT + RunSolverForConfigAndExpectResidualsMatch( + ThreadedSolverConfig(ITERATIVE_SCHUR, + SUITE_SPARSE, + kAutomaticOrdering, + CLUSTER_JACOBI)); +} + +TEST_F(BundleAdjustmentTest, + MultiThreadedIterativeSchurWithClusterJacobiAndUserOrderingUsingSuiteSparse) { // NOLINT + RunSolverForConfigAndExpectResidualsMatch( + ThreadedSolverConfig(ITERATIVE_SCHUR, + SUITE_SPARSE, + kUserOrdering, + CLUSTER_JACOBI)); +} + +TEST_F(BundleAdjustmentTest, + MultiThreadedIterativeSchurWithClusterTridiagonalAndAutomaticOrderingUsingSuiteSparse) { // NOLINT + RunSolverForConfigAndExpectResidualsMatch( + ThreadedSolverConfig(ITERATIVE_SCHUR, + SUITE_SPARSE, + kAutomaticOrdering, + CLUSTER_TRIDIAGONAL)); +} + +TEST_F(BundleAdjustmentTest, + MultiThreadedIterativeSchurWithClusterTridiagonalAndUserOrderingUsingSuiteSparse) { // NOTLINT + RunSolverForConfigAndExpectResidualsMatch( + ThreadedSolverConfig(ITERATIVE_SCHUR, + SUITE_SPARSE, + kUserOrdering, + CLUSTER_TRIDIAGONAL)); +} +#endif // CERES_NO_SUITESPARSE + +#ifndef CERES_NO_CXSPARSE +TEST_F(BundleAdjustmentTest, + MultiThreadedSparseNormalCholeskyWithAutomaticOrderingUsingCXSparse) { + RunSolverForConfigAndExpectResidualsMatch( + ThreadedSolverConfig(SPARSE_NORMAL_CHOLESKY, + CX_SPARSE, + kAutomaticOrdering)); +} + +TEST_F(BundleAdjustmentTest, + MultiThreadedSparseNormalCholeskyWithUserOrderingUsingCXSparse) { + RunSolverForConfigAndExpectResidualsMatch( + ThreadedSolverConfig(SPARSE_NORMAL_CHOLESKY, CX_SPARSE, kUserOrdering)); +} + +TEST_F(BundleAdjustmentTest, + MultiThreadedSparseSchurWithAutomaticOrderingUsingCXSparse) { + RunSolverForConfigAndExpectResidualsMatch( + ThreadedSolverConfig(SPARSE_SCHUR, CX_SPARSE, kAutomaticOrdering)); +} + +TEST_F(BundleAdjustmentTest, + MultiThreadedSparseSchurWithUserOrderingUsingCXSparse) { + RunSolverForConfigAndExpectResidualsMatch( + ThreadedSolverConfig(SPARSE_SCHUR, CX_SPARSE, kUserOrdering)); +} +#endif // CERES_NO_CXSPARSE + +#ifdef CERES_USE_EIGEN_SPARSE +TEST_F(BundleAdjustmentTest, + MultiThreadedSparseNormalCholeskyWithAutomaticOrderingUsingEigenSparse) { + RunSolverForConfigAndExpectResidualsMatch( + ThreadedSolverConfig(SPARSE_NORMAL_CHOLESKY, + EIGEN_SPARSE, + kAutomaticOrdering)); +} + +TEST_F(BundleAdjustmentTest, + MultiThreadedSparseNormalCholeskyWithUserOrderingUsingEigenSparse) { + RunSolverForConfigAndExpectResidualsMatch( + ThreadedSolverConfig(SPARSE_NORMAL_CHOLESKY, + EIGEN_SPARSE, + kUserOrdering)); +} + +TEST_F(BundleAdjustmentTest, + MultiThreadedSparseSchurWithAutomaticOrderingUsingEigenSparse) { + RunSolverForConfigAndExpectResidualsMatch( + ThreadedSolverConfig(SPARSE_SCHUR, EIGEN_SPARSE, kAutomaticOrdering)); +} + +TEST_F(BundleAdjustmentTest, + MultiThreadedSparseSchurWithUserOrderingUsingEigenSparse) { + RunSolverForConfigAndExpectResidualsMatch( + ThreadedSolverConfig(SPARSE_SCHUR, EIGEN_SPARSE, kUserOrdering)); +} +#endif // CERES_USE_EIGEN_SPARSE +#endif // CERES_USE_OPENMP + +} // namespace internal +} // namespace ceres
diff --git a/internal/ceres/system_test.cc b/internal/ceres/system_test.cc index 51812cc..5fb608e 100644 --- a/internal/ceres/system_test.cc +++ b/internal/ceres/system_test.cc
@@ -29,173 +29,22 @@ // Author: keir@google.com (Keir Mierle) // sameeragarwal@google.com (Sameer Agarwal) // -// System level tests for Ceres. The current suite of two tests. The -// first test is a small test based on Powell's Function. It is a -// scalar problem with 4 variables. The second problem is a bundle -// adjustment problem with 16 cameras and two thousand cameras. The -// first problem is to test the sanity test the factorization based -// solvers. The second problem is used to test the various -// combinations of solvers, orderings, preconditioners and -// multithreading. +// End-to-end tests for Ceres using Powell's function. #include <cmath> -#include <cstdio> #include <cstdlib> -#include <string> - -#include "ceres/internal/port.h" #include "ceres/autodiff_cost_function.h" -#include "ceres/ordered_groups.h" #include "ceres/problem.h" -#include "ceres/rotation.h" #include "ceres/solver.h" -#include "ceres/stringprintf.h" #include "ceres/test_util.h" #include "ceres/types.h" -#include "gflags/gflags.h" #include "glog/logging.h" #include "gtest/gtest.h" namespace ceres { namespace internal { -using std::string; -using std::vector; - -const bool kAutomaticOrdering = true; -const bool kUserOrdering = false; - -// Struct used for configuring the solver. -struct SolverConfig { - SolverConfig( - LinearSolverType linear_solver_type, - SparseLinearAlgebraLibraryType sparse_linear_algebra_library_type, - bool use_automatic_ordering) - : linear_solver_type(linear_solver_type), - sparse_linear_algebra_library_type(sparse_linear_algebra_library_type), - use_automatic_ordering(use_automatic_ordering), - preconditioner_type(IDENTITY), - num_threads(1) { - } - - SolverConfig( - LinearSolverType linear_solver_type, - SparseLinearAlgebraLibraryType sparse_linear_algebra_library_type, - bool use_automatic_ordering, - PreconditionerType preconditioner_type) - : linear_solver_type(linear_solver_type), - sparse_linear_algebra_library_type(sparse_linear_algebra_library_type), - use_automatic_ordering(use_automatic_ordering), - preconditioner_type(preconditioner_type), - num_threads(1) { - } - - string ToString() const { - return StringPrintf( - "(%s, %s, %s, %s, %d)", - LinearSolverTypeToString(linear_solver_type), - SparseLinearAlgebraLibraryTypeToString( - sparse_linear_algebra_library_type), - use_automatic_ordering ? "AUTOMATIC" : "USER", - PreconditionerTypeToString(preconditioner_type), - num_threads); - } - - LinearSolverType linear_solver_type; - SparseLinearAlgebraLibraryType sparse_linear_algebra_library_type; - bool use_automatic_ordering; - PreconditionerType preconditioner_type; - int num_threads; -}; - -// Templated function that given a set of solver configurations, -// instantiates a new copy of SystemTestProblem for each configuration -// and solves it. The solutions are expected to have residuals with -// coordinate-wise maximum absolute difference less than or equal to -// max_abs_difference. -// -// The template parameter SystemTestProblem is expected to implement -// the following interface. -// -// class SystemTestProblem { -// public: -// SystemTestProblem(); -// Problem* mutable_problem(); -// Solver::Options* mutable_solver_options(); -// }; -template <typename SystemTestProblem> -void RunSolversAndCheckTheyMatch( - const vector<SolverConfig>& configurations, - const double max_abs_difference) { - int num_configurations = configurations.size(); - vector<SystemTestProblem*> problems; - vector<vector<double> > final_residuals(num_configurations); - - for (int i = 0; i < num_configurations; ++i) { - SystemTestProblem* system_test_problem = new SystemTestProblem(); - - const SolverConfig& config = configurations[i]; - - Solver::Options& options = *(system_test_problem->mutable_solver_options()); - options.linear_solver_type = config.linear_solver_type; - options.sparse_linear_algebra_library_type = - config.sparse_linear_algebra_library_type; - options.preconditioner_type = config.preconditioner_type; - options.num_threads = config.num_threads; - options.num_linear_solver_threads = config.num_threads; - - if (config.use_automatic_ordering) { - options.linear_solver_ordering.reset(); - } - - LOG(INFO) << "Running solver configuration: " - << config.ToString(); - - Solver::Summary summary; - Solve(options, - system_test_problem->mutable_problem(), - &summary); - - system_test_problem - ->mutable_problem() - ->Evaluate(Problem::EvaluateOptions(), - NULL, - &final_residuals[i], - NULL, - NULL); - - CHECK_NE(summary.termination_type, ceres::FAILURE) - << "Solver configuration " << i << " failed."; - problems.push_back(system_test_problem); - - // Compare the resulting solutions to each other. Arbitrarily take - // SPARSE_NORMAL_CHOLESKY as the golden solve. We compare - // solutions by comparing their residual vectors. We do not - // compare parameter vectors because it is much more brittle and - // error prone to do so, since the same problem can have nearly - // the same residuals at two completely different positions in - // parameter space. - if (i > 0) { - const vector<double>& reference_residuals = final_residuals[0]; - const vector<double>& current_residuals = final_residuals[i]; - - for (int j = 0; j < reference_residuals.size(); ++j) { - EXPECT_NEAR(current_residuals[j], - reference_residuals[j], - max_abs_difference) - << "Not close enough residual:" << j - << " reference " << reference_residuals[j] - << " current " << current_residuals[j]; - } - } - } - - for (int i = 0; i < num_configurations; ++i) { - delete problems[i]; - } -} - // This class implements the SystemTestProblem interface and provides // access to an implementation of Powell's singular function. // @@ -229,12 +78,17 @@ problem_.AddResidualBlock( new AutoDiffCostFunction<F4, 1, 1, 1>(new F4), NULL, &x_[0], &x_[3]); + // Settings for the reference solution. + options_.linear_solver_type = ceres::DENSE_QR; options_.max_num_iterations = 10; + options_.num_threads = 1; } Problem* mutable_problem() { return &problem_; } Solver::Options* mutable_solver_options() { return &options_; } + static double kResidualTolerance; + private: // Templated functions used for automatically differentiated cost // functions. @@ -287,274 +141,51 @@ Solver::Options options_; }; -TEST(SystemTest, PowellsFunction) { - vector<SolverConfig> configs; -#define CONFIGURE(linear_solver, sparse_linear_algebra_library_type, ordering) \ - configs.push_back(SolverConfig(linear_solver, \ - sparse_linear_algebra_library_type, \ - ordering)) +double PowellsFunction::kResidualTolerance = 1e-8; - CONFIGURE(DENSE_QR, SUITE_SPARSE, kAutomaticOrdering); - CONFIGURE(DENSE_NORMAL_CHOLESKY, SUITE_SPARSE, kAutomaticOrdering); - CONFIGURE(DENSE_SCHUR, SUITE_SPARSE, kAutomaticOrdering); +typedef SystemTest<PowellsFunction> PowellTest; +const bool kAutomaticOrdering = true; -#ifndef CERES_NO_SUITESPARSE - CONFIGURE(SPARSE_NORMAL_CHOLESKY, SUITE_SPARSE, kAutomaticOrdering); -#endif // CERES_NO_SUITESPARSE - -#ifndef CERES_NO_CXSPARSE - CONFIGURE(SPARSE_NORMAL_CHOLESKY, CX_SPARSE, kAutomaticOrdering); -#endif // CERES_NO_CXSPARSE - - CONFIGURE(ITERATIVE_SCHUR, SUITE_SPARSE, kAutomaticOrdering); - -#undef CONFIGURE - - const double kMaxAbsoluteDifference = 1e-8; - RunSolversAndCheckTheyMatch<PowellsFunction>(configs, kMaxAbsoluteDifference); +TEST_F(PowellTest, DenseQR) { + RunSolverForConfigAndExpectResidualsMatch( + SolverConfig(DENSE_QR, NO_SPARSE)); } -// This class implements the SystemTestProblem interface and provides -// access to a bundle adjustment problem. It is based on -// examples/bundle_adjustment_example.cc. Currently a small 16 camera -// problem is hard coded in the constructor. Going forward we may -// extend this to a larger number of problems. -class BundleAdjustmentProblem { - public: - BundleAdjustmentProblem() { - const string input_file = TestFileAbsolutePath("problem-16-22106-pre.txt"); - ReadData(input_file); - BuildProblem(); - } +TEST_F(PowellTest, DenseNormalCholesky) { + RunSolverForConfigAndExpectResidualsMatch( + SolverConfig(DENSE_NORMAL_CHOLESKY)); +} - ~BundleAdjustmentProblem() { - delete []point_index_; - delete []camera_index_; - delete []observations_; - delete []parameters_; - } +TEST_F(PowellTest, DenseSchur) { + RunSolverForConfigAndExpectResidualsMatch( + SolverConfig(DENSE_SCHUR)); +} - Problem* mutable_problem() { return &problem_; } - Solver::Options* mutable_solver_options() { return &options_; } - - int num_cameras() const { return num_cameras_; } - int num_points() const { return num_points_; } - int num_observations() const { return num_observations_; } - const int* point_index() const { return point_index_; } - const int* camera_index() const { return camera_index_; } - const double* observations() const { return observations_; } - double* mutable_cameras() { return parameters_; } - double* mutable_points() { return parameters_ + 9 * num_cameras_; } - - private: - void ReadData(const string& filename) { - FILE * fptr = fopen(filename.c_str(), "r"); - - if (!fptr) { - LOG(FATAL) << "File Error: unable to open file " << filename; - } - - // This will die horribly on invalid files. Them's the breaks. - FscanfOrDie(fptr, "%d", &num_cameras_); - FscanfOrDie(fptr, "%d", &num_points_); - FscanfOrDie(fptr, "%d", &num_observations_); - - VLOG(1) << "Header: " << num_cameras_ - << " " << num_points_ - << " " << num_observations_; - - point_index_ = new int[num_observations_]; - camera_index_ = new int[num_observations_]; - observations_ = new double[2 * num_observations_]; - - num_parameters_ = 9 * num_cameras_ + 3 * num_points_; - parameters_ = new double[num_parameters_]; - - for (int i = 0; i < num_observations_; ++i) { - FscanfOrDie(fptr, "%d", camera_index_ + i); - FscanfOrDie(fptr, "%d", point_index_ + i); - for (int j = 0; j < 2; ++j) { - FscanfOrDie(fptr, "%lf", observations_ + 2*i + j); - } - } - - for (int i = 0; i < num_parameters_; ++i) { - FscanfOrDie(fptr, "%lf", parameters_ + i); - } - } - - void BuildProblem() { - double* points = mutable_points(); - double* cameras = mutable_cameras(); - - for (int i = 0; i < num_observations(); ++i) { - // Each Residual block takes a point and a camera as input and - // outputs a 2 dimensional residual. - CostFunction* cost_function = - new AutoDiffCostFunction<BundlerResidual, 2, 9, 3>( - new BundlerResidual(observations_[2*i + 0], - observations_[2*i + 1])); - - // Each observation correponds to a pair of a camera and a point - // which are identified by camera_index()[i] and - // point_index()[i] respectively. - double* camera = cameras + 9 * camera_index_[i]; - double* point = points + 3 * point_index()[i]; - problem_.AddResidualBlock(cost_function, NULL, camera, point); - } - - options_.linear_solver_ordering.reset(new ParameterBlockOrdering); - - // The points come before the cameras. - for (int i = 0; i < num_points_; ++i) { - options_.linear_solver_ordering->AddElementToGroup(points + 3 * i, 0); - } - - for (int i = 0; i < num_cameras_; ++i) { - options_.linear_solver_ordering->AddElementToGroup(cameras + 9 * i, 1); - } - - options_.max_num_iterations = 25; - options_.function_tolerance = 1e-10; - options_.gradient_tolerance = 1e-10; - options_.parameter_tolerance = 1e-10; - } - - template<typename T> - void FscanfOrDie(FILE *fptr, const char *format, T *value) { - int num_scanned = fscanf(fptr, format, value); - if (num_scanned != 1) { - LOG(FATAL) << "Invalid UW data file."; - } - } - - // Templated pinhole camera model. The camera is parameterized - // using 9 parameters. 3 for rotation, 3 for translation, 1 for - // focal length and 2 for radial distortion. The principal point is - // not modeled (i.e. it is assumed be located at the image center). - struct BundlerResidual { - // (u, v): the position of the observation with respect to the image - // center point. - BundlerResidual(double u, double v): u(u), v(v) {} - - template <typename T> - bool operator()(const T* const camera, - const T* const point, - T* residuals) const { - T p[3]; - AngleAxisRotatePoint(camera, point, p); - - // Add the translation vector - p[0] += camera[3]; - p[1] += camera[4]; - p[2] += camera[5]; - - const T& focal = camera[6]; - const T& l1 = camera[7]; - const T& l2 = camera[8]; - - // Compute the center of distortion. The sign change comes from - // the camera model that Noah Snavely's Bundler assumes, whereby - // the camera coordinate system has a negative z axis. - T xp = - focal * p[0] / p[2]; - T yp = - focal * p[1] / p[2]; - - // Apply second and fourth order radial distortion. - T r2 = xp*xp + yp*yp; - T distortion = T(1.0) + r2 * (l1 + l2 * r2); - - residuals[0] = distortion * xp - T(u); - residuals[1] = distortion * yp - T(v); - - return true; - } - - double u; - double v; - }; - - - Problem problem_; - Solver::Options options_; - - int num_cameras_; - int num_points_; - int num_observations_; - int num_parameters_; - - int* point_index_; - int* camera_index_; - double* observations_; - // The parameter vector is laid out as follows - // [camera_1, ..., camera_n, point_1, ..., point_m] - double* parameters_; -}; - -TEST(SystemTest, BundleAdjustmentProblem) { - vector<SolverConfig> configs; - -#define CONFIGURE(linear_solver, sparse_linear_algebra_library_type, ordering, preconditioner) \ - configs.push_back(SolverConfig(linear_solver, \ - sparse_linear_algebra_library_type, \ - ordering, \ - preconditioner)) - - CONFIGURE(DENSE_SCHUR, SUITE_SPARSE, kAutomaticOrdering, IDENTITY); - CONFIGURE(DENSE_SCHUR, SUITE_SPARSE, kUserOrdering, IDENTITY); - - CONFIGURE(CGNR, SUITE_SPARSE, kAutomaticOrdering, JACOBI); - - CONFIGURE(ITERATIVE_SCHUR, SUITE_SPARSE, kUserOrdering, JACOBI); - CONFIGURE(ITERATIVE_SCHUR, SUITE_SPARSE, kAutomaticOrdering, JACOBI); - - CONFIGURE(ITERATIVE_SCHUR, SUITE_SPARSE, kUserOrdering, SCHUR_JACOBI); - CONFIGURE(ITERATIVE_SCHUR, SUITE_SPARSE, kAutomaticOrdering, SCHUR_JACOBI); +TEST_F(PowellTest, IterativeSchurWithJacobi) { + RunSolverForConfigAndExpectResidualsMatch( + SolverConfig(ITERATIVE_SCHUR, NO_SPARSE, kAutomaticOrdering, JACOBI)); +} #ifndef CERES_NO_SUITESPARSE - CONFIGURE(SPARSE_NORMAL_CHOLESKY, SUITE_SPARSE, kAutomaticOrdering, IDENTITY); - CONFIGURE(SPARSE_NORMAL_CHOLESKY, SUITE_SPARSE, kUserOrdering, IDENTITY); - - CONFIGURE(SPARSE_SCHUR, SUITE_SPARSE, kAutomaticOrdering, IDENTITY); - CONFIGURE(SPARSE_SCHUR, SUITE_SPARSE, kUserOrdering, IDENTITY); - - CONFIGURE(ITERATIVE_SCHUR, SUITE_SPARSE, kAutomaticOrdering, CLUSTER_JACOBI); - CONFIGURE(ITERATIVE_SCHUR, SUITE_SPARSE, kUserOrdering, CLUSTER_JACOBI); - - CONFIGURE(ITERATIVE_SCHUR, SUITE_SPARSE, kAutomaticOrdering, CLUSTER_TRIDIAGONAL); - CONFIGURE(ITERATIVE_SCHUR, SUITE_SPARSE, kUserOrdering, CLUSTER_TRIDIAGONAL); +TEST_F(PowellTest, SparseNormalCholeskyUsingSuiteSparse) { + RunSolverForConfigAndExpectResidualsMatch( + SolverConfig(SPARSE_NORMAL_CHOLESKY, SUITE_SPARSE, kAutomaticOrdering)); +} #endif // CERES_NO_SUITESPARSE #ifndef CERES_NO_CXSPARSE - CONFIGURE(SPARSE_NORMAL_CHOLESKY, CX_SPARSE, kAutomaticOrdering, IDENTITY); - CONFIGURE(SPARSE_NORMAL_CHOLESKY, CX_SPARSE, kUserOrdering, IDENTITY); - - CONFIGURE(SPARSE_SCHUR, CX_SPARSE, kAutomaticOrdering, IDENTITY); - CONFIGURE(SPARSE_SCHUR, CX_SPARSE, kUserOrdering, IDENTITY); +TEST_F(PowellTest, SparseNormalCholeskyUsingCXSparse) { + RunSolverForConfigAndExpectResidualsMatch( + SolverConfig(SPARSE_NORMAL_CHOLESKY, CX_SPARSE, kAutomaticOrdering)); +} #endif // CERES_NO_CXSPARSE #ifdef CERES_USE_EIGEN_SPARSE - CONFIGURE(SPARSE_SCHUR, EIGEN_SPARSE, kAutomaticOrdering, IDENTITY); - CONFIGURE(SPARSE_SCHUR, EIGEN_SPARSE, kUserOrdering, IDENTITY); - CONFIGURE(SPARSE_NORMAL_CHOLESKY, EIGEN_SPARSE, kAutomaticOrdering, IDENTITY); - CONFIGURE(SPARSE_NORMAL_CHOLESKY, EIGEN_SPARSE, kUserOrdering, IDENTITY); -#endif // CERES_USE_EIGEN_SPARSE - -#undef CONFIGURE - - // Single threaded evaluators and linear solvers. - const double kMaxAbsoluteDifference = 1e-4; - RunSolversAndCheckTheyMatch<BundleAdjustmentProblem>(configs, - kMaxAbsoluteDifference); - -#ifdef CERES_USE_OPENMP - // Multithreaded evaluators and linear solvers. - for (int i = 0; i < configs.size(); ++i) { - configs[i].num_threads = 2; - } - RunSolversAndCheckTheyMatch<BundleAdjustmentProblem>(configs, - kMaxAbsoluteDifference); -#endif // CERES_USE_OPENMP +TEST_F(PowellTest, SparseNormalCholeskyUsingEigenSparse) { + RunSolverForConfigAndExpectResidualsMatch( + SolverConfig(SPARSE_NORMAL_CHOLESKY, EIGEN_SPARSE, kAutomaticOrdering)); } +#endif // CERES_USE_EIGEN_SPARSE } // namespace internal } // namespace ceres
diff --git a/internal/ceres/test_util.cc b/internal/ceres/test_util.cc index 695d4d8..5edd4fc 100644 --- a/internal/ceres/test_util.cc +++ b/internal/ceres/test_util.cc
@@ -36,6 +36,7 @@ #include <cmath> #include "ceres/file.h" #include "ceres/stringprintf.h" +#include "ceres/types.h" #include "gflags/gflags.h" #include "glog/logging.h" #include "gtest/gtest.h" @@ -125,6 +126,18 @@ filename); } +SolverConfig ThreadedSolverConfig( + LinearSolverType linear_solver_type, + SparseLinearAlgebraLibraryType sparse_linear_algebra_library_type, + bool use_automatic_ordering, + PreconditionerType preconditioner_type) { + const int kNumThreads = 4; + return SolverConfig(linear_solver_type, + sparse_linear_algebra_library_type, + use_automatic_ordering, + preconditioner_type, + kNumThreads); +} } // namespace internal } // namespace ceres
diff --git a/internal/ceres/test_util.h b/internal/ceres/test_util.h index 65cc7dc..6aff541 100644 --- a/internal/ceres/test_util.h +++ b/internal/ceres/test_util.h
@@ -30,6 +30,11 @@ #include <string> #include "ceres/internal/port.h" +#include "ceres/problem.h" +#include "ceres/solver.h" +#include "ceres/stringprintf.h" +#include "gtest/gtest.h" + #ifndef CERES_INTERNAL_TEST_UTIL_H_ #define CERES_INTERNAL_TEST_UTIL_H_ @@ -65,6 +70,117 @@ // local build/testing environment. std::string TestFileAbsolutePath(const std::string& filename); +// Struct used for configuring the solver. Used by end-to-end tests. +struct SolverConfig { + SolverConfig( + LinearSolverType linear_solver_type, + SparseLinearAlgebraLibraryType sparse_linear_algebra_library_type = NO_SPARSE, + bool use_automatic_ordering = true, + PreconditionerType preconditioner_type = IDENTITY, + int num_threads = 1) + : linear_solver_type(linear_solver_type), + sparse_linear_algebra_library_type(sparse_linear_algebra_library_type), + use_automatic_ordering(use_automatic_ordering), + preconditioner_type(preconditioner_type), + num_threads(num_threads) { + } + + std::string ToString() const { + return StringPrintf( + "(%s, %s, %s, %s, %d)", + LinearSolverTypeToString(linear_solver_type), + SparseLinearAlgebraLibraryTypeToString( + sparse_linear_algebra_library_type), + use_automatic_ordering ? "AUTOMATIC" : "USER", + PreconditionerTypeToString(preconditioner_type), + num_threads); + } + + void UpdateOptions(Solver::Options* options) const { + options->linear_solver_type = linear_solver_type; + options->sparse_linear_algebra_library_type = + sparse_linear_algebra_library_type; + options->preconditioner_type = preconditioner_type; + options->num_threads = num_threads; + options->num_linear_solver_threads = num_threads; + + if (use_automatic_ordering) { + options->linear_solver_ordering.reset(); + } + } + + LinearSolverType linear_solver_type; + SparseLinearAlgebraLibraryType sparse_linear_algebra_library_type; + bool use_automatic_ordering; + PreconditionerType preconditioner_type; + int num_threads; +}; + +SolverConfig ThreadedSolverConfig( + LinearSolverType linear_solver_type, + SparseLinearAlgebraLibraryType sparse_linear_algebra_library_type = NO_SPARSE, + bool use_automatic_ordering = true, + PreconditionerType preconditioner_type = IDENTITY); + +// A templated test fixture, that is used for testing Ceres end to end +// by computing a solution to the problem for a given solver +// configuration and comparing it to a reference solver configuration. +// +// It is assumed that the SystemTestProblem has an Solver::Options +// struct that contains the reference Solver configuration. +template <class SystemTestProblem> +class SystemTest : public::testing::Test { + protected: + virtual void SetUp() { + SystemTestProblem system_test_problem; + SolveAndEvaluateFinalResiduals( + *system_test_problem.mutable_solver_options(), + system_test_problem.mutable_problem(), + &expected_final_residuals_); + } + + void RunSolverForConfigAndExpectResidualsMatch(const SolverConfig& config) { + LOG(INFO) << "Running solver configuration: " + << config.ToString(); + + SystemTestProblem system_test_problem; + config.UpdateOptions(system_test_problem.mutable_solver_options()); + std::vector<double> final_residuals; + SolveAndEvaluateFinalResiduals( + *system_test_problem.mutable_solver_options(), + system_test_problem.mutable_problem(), + &final_residuals); + + // We compare solutions by comparing their residual vectors. We do + // not compare parameter vectors because it is much more brittle + // and error prone to do so, since the same problem can have + // nearly the same residuals at two completely different positions + // in parameter space. + CHECK_EQ(expected_final_residuals_.size(), final_residuals.size()); + for (int i = 0; i < final_residuals.size(); ++i) { + EXPECT_NEAR(final_residuals[i], + expected_final_residuals_[i], + SystemTestProblem::kResidualTolerance) + << "Not close enough residual:" << i; + } + } + + void SolveAndEvaluateFinalResiduals(const Solver::Options& options, + Problem* problem, + std::vector<double>* final_residuals) { + Solver::Summary summary; + Solve(options, problem, &summary); + CHECK_NE(summary.termination_type, ceres::FAILURE); + problem->Evaluate(Problem::EvaluateOptions(), + NULL, + final_residuals, + NULL, + NULL); + } + + std::vector<double> expected_final_residuals_; +}; + } // namespace internal } // namespace ceres