Update ceres to the latest version in g3android-wear-5.0.0_r1android-cts-5.1_r9android-cts-5.1_r8android-cts-5.1_r7android-cts-5.1_r6android-cts-5.1_r5android-cts-5.1_r4android-cts-5.1_r3android-cts-5.1_r28android-cts-5.1_r27android-cts-5.1_r26android-cts-5.1_r25android-cts-5.1_r24android-cts-5.1_r23android-cts-5.1_r22android-cts-5.1_r21android-cts-5.1_r20android-cts-5.1_r2android-cts-5.1_r19android-cts-5.1_r18android-cts-5.1_r17android-cts-5.1_r16android-cts-5.1_r15android-cts-5.1_r14android-cts-5.1_r13android-cts-5.1_r10android-cts-5.1_r1android-cts-5.0_r9android-cts-5.0_r8android-cts-5.0_r7android-cts-5.0_r6android-cts-5.0_r5android-cts-5.0_r4android-cts-5.0_r3android-5.1.1_r9android-5.1.1_r8android-5.1.1_r7android-5.1.1_r6android-5.1.1_r5android-5.1.1_r4android-5.1.1_r38android-5.1.1_r37android-5.1.1_r36android-5.1.1_r35android-5.1.1_r34android-5.1.1_r33android-5.1.1_r30android-5.1.1_r3android-5.1.1_r29android-5.1.1_r28android-5.1.1_r26android-5.1.1_r25android-5.1.1_r24android-5.1.1_r23android-5.1.1_r22android-5.1.1_r20android-5.1.1_r2android-5.1.1_r19android-5.1.1_r18android-5.1.1_r17android-5.1.1_r16android-5.1.1_r15android-5.1.1_r14android-5.1.1_r13android-5.1.1_r12android-5.1.1_r10android-5.1.1_r1android-5.1.0_r5android-5.1.0_r4android-5.1.0_r3android-5.1.0_r1android-5.0.2_r3android-5.0.2_r1android-5.0.1_r1android-5.0.0_r7android-5.0.0_r6android-5.0.0_r5.1android-5.0.0_r5android-5.0.0_r4android-5.0.0_r3android-5.0.0_r2android-5.0.0_r1lollipop-wear-releaselollipop-releaselollipop-mr1-wfc-releaselollipop-mr1-releaselollipop-mr1-fi-releaselollipop-mr1-devlollipop-mr1-cts-releaselollipop-devlollipop-cts-release

Please pay special attention to the changes in Android.mk.
They are the only real changes I had to make.

Bug: 16953678

Change-Id: I44a644358e779aaff99a2ea822387fe49ac26888

15 files changed, 1690 insertions, 482 deletions

diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt
index 9dfc80b..dbbcb81 100644
--- a/examples/CMakeLists.txt
+++ b/examples/CMakeLists.txt
@@ -28,6 +28,13 @@
 #
 # Author: keir@google.com (Keir Mierle)
 
+# Only Ceres itself should be compiled with CERES_BUILDING_SHARED_LIBRARY
+# defined, any users of Ceres will have CERES_USING_SHARED_LIBRARY defined
+# for them in Ceres' config.h if appropriate.
+IF (BUILD_SHARED_LIBS)
+  REMOVE_DEFINITIONS(-DCERES_BUILDING_SHARED_LIBRARY)
+ENDIF()
+
 ADD_EXECUTABLE(helloworld helloworld.cc)
 TARGET_LINK_LIBRARIES(helloworld ceres)
 
@@ -42,7 +49,15 @@ TARGET_LINK_LIBRARIES(curve_fitting ceres)
 
 ADD_EXECUTABLE(curve_fitting_c curve_fitting.c)
 TARGET_LINK_LIBRARIES(curve_fitting_c ceres)
+# As this is a C file #including <math.h> we have to explicitly add the math
+# library (libm). Although some compilers (dependent upon options) will accept
+# the indirect link to libm via Ceres, at least GCC 4.8 on pure Debian won't.
+IF (NOT MSVC)
+  TARGET_LINK_LIBRARIES(curve_fitting_c m)
+ENDIF (NOT MSVC)
 
+ADD_EXECUTABLE(ellipse_approximation ellipse_approximation.cc)
+TARGET_LINK_LIBRARIES(ellipse_approximation ceres)
 
 ADD_EXECUTABLE(robust_curve_fitting robust_curve_fitting.cc)
 TARGET_LINK_LIBRARIES(robust_curve_fitting ceres)
@@ -53,26 +68,37 @@ TARGET_LINK_LIBRARIES(simple_bundle_adjuster ceres)
 
 IF (GFLAGS)
   ADD_EXECUTABLE(powell powell.cc)
-  TARGET_LINK_LIBRARIES(powell ceres)
+  TARGET_LINK_LIBRARIES(powell ceres ${GFLAGS_LIBRARIES})
 
   ADD_EXECUTABLE(nist nist.cc)
-  TARGET_LINK_LIBRARIES(nist ceres)
+  TARGET_LINK_LIBRARIES(nist ceres ${GFLAGS_LIBRARIES})
+
+  ADD_EXECUTABLE(more_garbow_hillstrom more_garbow_hillstrom.cc)
+  TARGET_LINK_LIBRARIES(more_garbow_hillstrom ceres ${GFLAGS_LIBRARIES})
 
   ADD_EXECUTABLE(circle_fit circle_fit.cc)
-  TARGET_LINK_LIBRARIES(circle_fit ceres)
+  TARGET_LINK_LIBRARIES(circle_fit ceres ${GFLAGS_LIBRARIES})
 
   ADD_EXECUTABLE(bundle_adjuster
                  bundle_adjuster.cc
                  bal_problem.cc)
-  TARGET_LINK_LIBRARIES(bundle_adjuster ceres)
+  TARGET_LINK_LIBRARIES(bundle_adjuster ceres ${GFLAGS_LIBRARIES})
 
   ADD_EXECUTABLE(libmv_bundle_adjuster
                  libmv_bundle_adjuster.cc)
-  TARGET_LINK_LIBRARIES(libmv_bundle_adjuster ceres)
+  TARGET_LINK_LIBRARIES(libmv_bundle_adjuster ceres ${GFLAGS_LIBRARIES})
+
+  ADD_EXECUTABLE(libmv_homography
+                 libmv_homography.cc)
+  TARGET_LINK_LIBRARIES(libmv_homography ceres ${GFLAGS_LIBRARIES})
 
   ADD_EXECUTABLE(denoising
                  denoising.cc
                  fields_of_experts.cc)
-  TARGET_LINK_LIBRARIES(denoising ceres)
-ENDIF (GFLAGS)
+  TARGET_LINK_LIBRARIES(denoising ceres ${GFLAGS_LIBRARIES})
+
+  ADD_EXECUTABLE(robot_pose_mle
+                 robot_pose_mle.cc)
+  TARGET_LINK_LIBRARIES(robot_pose_mle ceres ${GFLAGS_LIBRARIES})
 
+ENDIF (GFLAGS)
diff --git a/examples/bal_problem.cc b/examples/bal_problem.cc
index c118f88..40d0fdf 100644
--- a/examples/bal_problem.cc
+++ b/examples/bal_problem.cc
@@ -37,6 +37,7 @@
 #include "Eigen/Core"
 #include "ceres/rotation.h"
 #include "glog/logging.h"
+#include "random.h"
 
 namespace ceres {
 namespace examples {
@@ -44,25 +45,6 @@ namespace {
 typedef Eigen::Map<Eigen::VectorXd> VectorRef;
 typedef Eigen::Map<const Eigen::VectorXd> ConstVectorRef;
 
-inline double RandDouble() {
-  double r = static_cast<double>(rand());
-  return r / RAND_MAX;
-}
-
-// Box-Muller algorithm for normal random number generation.
-// http://en.wikipedia.org/wiki/Box-Muller_transform
-inline double RandNormal() {
-  double x1, x2, w;
-  do {
-    x1 = 2.0 * RandDouble() - 1.0;
-    x2 = 2.0 * RandDouble() - 1.0;
-    w = x1 * x1 + x2 * x2;
-  } while ( w >= 1.0 || w == 0.0 );
-
-  w = sqrt((-2.0 * log(w)) / w);
-  return x1 * w;
-}
-
 template<typename T>
 void FscanfOrDie(FILE* fptr, const char* format, T* value) {
   int num_scanned = fscanf(fptr, format, value);
diff --git a/examples/bundle_adjuster.cc b/examples/bundle_adjuster.cc
index 224ad74..296611f 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",
@@ -113,7 +116,6 @@ DEFINE_double(point_sigma, 0.0, "Standard deviation of the point "
 DEFINE_int32(random_seed, 38401, "Random seed used to set the state "
              "of the pseudo random number generator used to generate "
              "the pertubations.");
-DEFINE_string(solver_log, "", "File to record the solver execution to.");
 DEFINE_bool(line_search, false, "Use a line search instead of trust region "
             "algorithm.");
 
@@ -125,6 +127,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));
@@ -146,19 +150,19 @@ void SetOrdering(BALProblem* bal_problem, Solver::Options* options) {
   if (options->use_inner_iterations) {
     if (FLAGS_blocks_for_inner_iterations == "cameras") {
       LOG(INFO) << "Camera blocks for inner iterations";
-      options->inner_iteration_ordering = new ParameterBlockOrdering;
+      options->inner_iteration_ordering.reset(new ParameterBlockOrdering);
       for (int i = 0; i < num_cameras; ++i) {
         options->inner_iteration_ordering->AddElementToGroup(cameras + camera_block_size * i, 0);
       }
     } else if (FLAGS_blocks_for_inner_iterations == "points") {
       LOG(INFO) << "Point blocks for inner iterations";
-      options->inner_iteration_ordering = new ParameterBlockOrdering;
+      options->inner_iteration_ordering.reset(new ParameterBlockOrdering);
       for (int i = 0; i < num_points; ++i) {
         options->inner_iteration_ordering->AddElementToGroup(points + point_block_size * i, 0);
       }
     } else if (FLAGS_blocks_for_inner_iterations == "cameras,points") {
       LOG(INFO) << "Camera followed by point blocks for inner iterations";
-      options->inner_iteration_ordering = new ParameterBlockOrdering;
+      options->inner_iteration_ordering.reset(new ParameterBlockOrdering);
       for (int i = 0; i < num_cameras; ++i) {
         options->inner_iteration_ordering->AddElementToGroup(cameras + camera_block_size * i, 0);
       }
@@ -167,7 +171,7 @@ void SetOrdering(BALProblem* bal_problem, Solver::Options* options) {
       }
     } else if (FLAGS_blocks_for_inner_iterations == "points,cameras") {
       LOG(INFO) << "Point followed by camera blocks for inner iterations";
-      options->inner_iteration_ordering = new ParameterBlockOrdering;
+      options->inner_iteration_ordering.reset(new ParameterBlockOrdering);
       for (int i = 0; i < num_cameras; ++i) {
         options->inner_iteration_ordering->AddElementToGroup(cameras + camera_block_size * i, 1);
       }
@@ -216,7 +220,7 @@ void SetOrdering(BALProblem* bal_problem, Solver::Options* options) {
     }
   }
 
-  options->linear_solver_ordering = ordering;
+  options->linear_solver_ordering.reset(ordering);
 }
 
 void SetMinimizerOptions(Solver::Options* options) {
@@ -258,18 +262,14 @@ void BuildProblem(BALProblem* bal_problem, Problem* problem) {
     CostFunction* cost_function;
     // Each Residual block takes a point and a camera as input and
     // outputs a 2 dimensional residual.
-    if (FLAGS_use_quaternions) {
-      cost_function = new AutoDiffCostFunction<
-          SnavelyReprojectionErrorWithQuaternions, 2, 4, 6, 3>(
-              new SnavelyReprojectionErrorWithQuaternions(
-                  observations[2 * i + 0],
-                  observations[2 * i + 1]));
-    } else {
-      cost_function =
-          new AutoDiffCostFunction<SnavelyReprojectionError, 2, 9, 3>(
-              new SnavelyReprojectionError(observations[2 * i + 0],
-                                           observations[2 * i + 1]));
-    }
+    cost_function =
+        (FLAGS_use_quaternions)
+        ? SnavelyReprojectionErrorWithQuaternions::Create(
+            observations[2 * i + 0],
+            observations[2 * i + 1])
+        : SnavelyReprojectionError::Create(
+            observations[2 * i + 0],
+            observations[2 * i + 1]);
 
     // If enabled use Huber's loss function.
     LossFunction* loss_function = FLAGS_robustify ? new HuberLoss(1.0) : NULL;
@@ -319,7 +319,6 @@ void SolveProblem(const char* filename) {
   BuildProblem(&bal_problem, &problem);
   Solver::Options options;
   SetSolverOptionsFromFlags(&bal_problem, &options);
-  options.solver_log = FLAGS_solver_log;
   options.gradient_tolerance = 1e-16;
   options.function_tolerance = 1e-16;
   Solver::Summary summary;
diff --git a/examples/data_fitting.cc b/examples/data_fitting.cc
deleted file mode 100644
index 5d54123..0000000
--- a/examples/data_fitting.cc
+++ /dev/null
@@ -1,165 +0,0 @@
-// Ceres Solver - A fast non-linear least squares minimizer
-// Copyright 2010, 2011, 2012 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)
-
-#include "ceres/ceres.h"
-#include "gflags/gflags.h"
-
-using ceres::AutoDiffCostFunction;
-using ceres::CostFunction;
-using ceres::Problem;
-using ceres::Solver;
-using ceres::Solve;
-
-// Data generated using the following octave code.
-//   randn('seed', 23497);
-//   m = 0.3;
-//   c = 0.1;
-//   x=[0:0.075:5];
-//   y = exp(m * x + c);
-//   noise = randn(size(x)) * 0.2;
-//   y_observed = y + noise;
-//   data = [x', y_observed'];
-
-const int kNumObservations = 67;
-const double data[] = {
-  0.000000e+00, 1.133898e+00,
-  7.500000e-02, 1.334902e+00,
-  1.500000e-01, 1.213546e+00,
-  2.250000e-01, 1.252016e+00,
-  3.000000e-01, 1.392265e+00,
-  3.750000e-01, 1.314458e+00,
-  4.500000e-01, 1.472541e+00,
-  5.250000e-01, 1.536218e+00,
-  6.000000e-01, 1.355679e+00,
-  6.750000e-01, 1.463566e+00,
-  7.500000e-01, 1.490201e+00,
-  8.250000e-01, 1.658699e+00,
-  9.000000e-01, 1.067574e+00,
-  9.750000e-01, 1.464629e+00,
-  1.050000e+00, 1.402653e+00,
-  1.125000e+00, 1.713141e+00,
-  1.200000e+00, 1.527021e+00,
-  1.275000e+00, 1.702632e+00,
-  1.350000e+00, 1.423899e+00,
-  1.425000e+00, 1.543078e+00,
-  1.500000e+00, 1.664015e+00,
-  1.575000e+00, 1.732484e+00,
-  1.650000e+00, 1.543296e+00,
-  1.725000e+00, 1.959523e+00,
-  1.800000e+00, 1.685132e+00,
-  1.875000e+00, 1.951791e+00,
-  1.950000e+00, 2.095346e+00,
-  2.025000e+00, 2.361460e+00,
-  2.100000e+00, 2.169119e+00,
-  2.175000e+00, 2.061745e+00,
-  2.250000e+00, 2.178641e+00,
-  2.325000e+00, 2.104346e+00,
-  2.400000e+00, 2.584470e+00,
-  2.475000e+00, 1.914158e+00,
-  2.550000e+00, 2.368375e+00,
-  2.625000e+00, 2.686125e+00,
-  2.700000e+00, 2.712395e+00,
-  2.775000e+00, 2.499511e+00,
-  2.850000e+00, 2.558897e+00,
-  2.925000e+00, 2.309154e+00,
-  3.000000e+00, 2.869503e+00,
-  3.075000e+00, 3.116645e+00,
-  3.150000e+00, 3.094907e+00,
-  3.225000e+00, 2.471759e+00,
-  3.300000e+00, 3.017131e+00,
-  3.375000e+00, 3.232381e+00,
-  3.450000e+00, 2.944596e+00,
-  3.525000e+00, 3.385343e+00,
-  3.600000e+00, 3.199826e+00,
-  3.675000e+00, 3.423039e+00,
-  3.750000e+00, 3.621552e+00,
-  3.825000e+00, 3.559255e+00,
-  3.900000e+00, 3.530713e+00,
-  3.975000e+00, 3.561766e+00,
-  4.050000e+00, 3.544574e+00,
-  4.125000e+00, 3.867945e+00,
-  4.200000e+00, 4.049776e+00,
-  4.275000e+00, 3.885601e+00,
-  4.350000e+00, 4.110505e+00,
-  4.425000e+00, 4.345320e+00,
-  4.500000e+00, 4.161241e+00,
-  4.575000e+00, 4.363407e+00,
-  4.650000e+00, 4.161576e+00,
-  4.725000e+00, 4.619728e+00,
-  4.800000e+00, 4.737410e+00,
-  4.875000e+00, 4.727863e+00,
-  4.950000e+00, 4.669206e+00,
-};
-
-class ExponentialResidual {
- public:
-  ExponentialResidual(double x, double y)
-      : x_(x), y_(y) {}
-
-  template <typename T> bool operator()(const T* const m,
-                                        const T* const c,
-                                        T* residual) const {
-    residual[0] = T(y_) - exp(m[0] * T(x_) + c[0]);
-    return true;
-  }
-
- private:
-  const double x_;
-  const double y_;
-};
-
-int main(int argc, char** argv) {
-  google::ParseCommandLineFlags(&argc, &argv, true);
-  google::InitGoogleLogging(argv[0]);
-
-  double m = 0.0;
-  double c = 0.0;
-
-  Problem problem;
-  for (int i = 0; i < kNumObservations; ++i) {
-    problem.AddResidualBlock(
-        new AutoDiffCostFunction<ExponentialResidual, 1, 1, 1>(
-            new ExponentialResidual(data[2 * i], data[2 * i + 1])),
-        NULL,
-        &m, &c);
-  }
-
-  Solver::Options options;
-  options.max_num_iterations = 25;
-  options.linear_solver_type = ceres::DENSE_QR;
-  options.minimizer_progress_to_stdout = true;
-
-  Solver::Summary summary;
-  Solve(options, &problem, &summary);
-  std::cout << summary.BriefReport() << "\n";
-  std::cout << "Initial m: " << 0.0 << " c: " << 0.0 << "\n";
-  std::cout << "Final   m: " << m << " c: " << c << "\n";
-  return 0;
-}
diff --git a/examples/ellipse_approximation.cc b/examples/ellipse_approximation.cc
new file mode 100644
index 0000000..a5bbe02
--- /dev/null
+++ b/examples/ellipse_approximation.cc
@@ -0,0 +1,451 @@
+// Ceres Solver - A fast non-linear least squares minimizer
+// Copyright 2014 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: richie.stebbing@gmail.com (Richard Stebbing)
+//
+// This fits points randomly distributed on an ellipse with an approximate
+// line segment contour. This is done by jointly optimizing the control points
+// of the line segment contour along with the preimage positions for the data
+// points. The purpose of this example is to show an example use case for
+// dynamic_sparsity, and how it can benefit problems which are numerically
+// dense but dynamically sparse.
+
+#include <cmath>
+#include <vector>
+#include "ceres/ceres.h"
+#include "glog/logging.h"
+
+// Data generated with the following Python code.
+//   import numpy as np
+//   np.random.seed(1337)
+//   t = np.linspace(0.0, 2.0 * np.pi, 212, endpoint=False)
+//   t += 2.0 * np.pi * 0.01 * np.random.randn(t.size)
+//   theta = np.deg2rad(15)
+//   a, b = np.cos(theta), np.sin(theta)
+//   R = np.array([[a, -b],
+//                 [b, a]])
+//   Y = np.dot(np.c_[4.0 * np.cos(t), np.sin(t)], R.T)
+
+const int kYRows = 212;
+const int kYCols = 2;
+const double kYData[kYRows * kYCols] = {
+  +3.871364e+00, +9.916027e-01,
+  +3.864003e+00, +1.034148e+00,
+  +3.850651e+00, +1.072202e+00,
+  +3.868350e+00, +1.014408e+00,
+  +3.796381e+00, +1.153021e+00,
+  +3.857138e+00, +1.056102e+00,
+  +3.787532e+00, +1.162215e+00,
+  +3.704477e+00, +1.227272e+00,
+  +3.564711e+00, +1.294959e+00,
+  +3.754363e+00, +1.191948e+00,
+  +3.482098e+00, +1.322725e+00,
+  +3.602777e+00, +1.279658e+00,
+  +3.585433e+00, +1.286858e+00,
+  +3.347505e+00, +1.356415e+00,
+  +3.220855e+00, +1.378914e+00,
+  +3.558808e+00, +1.297174e+00,
+  +3.403618e+00, +1.343809e+00,
+  +3.179828e+00, +1.384721e+00,
+  +3.054789e+00, +1.398759e+00,
+  +3.294153e+00, +1.366808e+00,
+  +3.247312e+00, +1.374813e+00,
+  +2.988547e+00, +1.404247e+00,
+  +3.114508e+00, +1.392698e+00,
+  +2.899226e+00, +1.409802e+00,
+  +2.533256e+00, +1.414778e+00,
+  +2.654773e+00, +1.415909e+00,
+  +2.565100e+00, +1.415313e+00,
+  +2.976456e+00, +1.405118e+00,
+  +2.484200e+00, +1.413640e+00,
+  +2.324751e+00, +1.407476e+00,
+  +1.930468e+00, +1.378221e+00,
+  +2.329017e+00, +1.407688e+00,
+  +1.760640e+00, +1.360319e+00,
+  +2.147375e+00, +1.396603e+00,
+  +1.741989e+00, +1.358178e+00,
+  +1.743859e+00, +1.358394e+00,
+  +1.557372e+00, +1.335208e+00,
+  +1.280551e+00, +1.295087e+00,
+  +1.429880e+00, +1.317546e+00,
+  +1.213485e+00, +1.284400e+00,
+  +9.168172e-01, +1.232870e+00,
+  +1.311141e+00, +1.299839e+00,
+  +1.231969e+00, +1.287382e+00,
+  +7.453773e-01, +1.200049e+00,
+  +6.151587e-01, +1.173683e+00,
+  +5.935666e-01, +1.169193e+00,
+  +2.538707e-01, +1.094227e+00,
+  +6.806136e-01, +1.187089e+00,
+  +2.805447e-01, +1.100405e+00,
+  +6.184807e-01, +1.174371e+00,
+  +1.170550e-01, +1.061762e+00,
+  +2.890507e-01, +1.102365e+00,
+  +3.834234e-01, +1.123772e+00,
+  +3.980161e-04, +1.033061e+00,
+  -3.651680e-01, +9.370367e-01,
+  -8.386351e-01, +7.987201e-01,
+  -8.105704e-01, +8.073702e-01,
+  -8.735139e-01, +7.878886e-01,
+  -9.913836e-01, +7.506100e-01,
+  -8.784011e-01, +7.863636e-01,
+  -1.181440e+00, +6.882566e-01,
+  -1.229556e+00, +6.720191e-01,
+  -1.035839e+00, +7.362765e-01,
+  -8.031520e-01, +8.096470e-01,
+  -1.539136e+00, +5.629549e-01,
+  -1.755423e+00, +4.817306e-01,
+  -1.337589e+00, +6.348763e-01,
+  -1.836966e+00, +4.499485e-01,
+  -1.913367e+00, +4.195617e-01,
+  -2.126467e+00, +3.314900e-01,
+  -1.927625e+00, +4.138238e-01,
+  -2.339862e+00, +2.379074e-01,
+  -1.881736e+00, +4.322152e-01,
+  -2.116753e+00, +3.356163e-01,
+  -2.255733e+00, +2.754930e-01,
+  -2.555834e+00, +1.368473e-01,
+  -2.770277e+00, +2.895711e-02,
+  -2.563376e+00, +1.331890e-01,
+  -2.826715e+00, -9.000818e-04,
+  -2.978191e+00, -8.457804e-02,
+  -3.115855e+00, -1.658786e-01,
+  -2.982049e+00, -8.678322e-02,
+  -3.307892e+00, -2.902083e-01,
+  -3.038346e+00, -1.194222e-01,
+  -3.190057e+00, -2.122060e-01,
+  -3.279086e+00, -2.705777e-01,
+  -3.322028e+00, -2.999889e-01,
+  -3.122576e+00, -1.699965e-01,
+  -3.551973e+00, -4.768674e-01,
+  -3.581866e+00, -5.032175e-01,
+  -3.497799e+00, -4.315203e-01,
+  -3.565384e+00, -4.885602e-01,
+  -3.699493e+00, -6.199815e-01,
+  -3.585166e+00, -5.061925e-01,
+  -3.758914e+00, -6.918275e-01,
+  -3.741104e+00, -6.689131e-01,
+  -3.688331e+00, -6.077239e-01,
+  -3.810425e+00, -7.689015e-01,
+  -3.791829e+00, -7.386911e-01,
+  -3.789951e+00, -7.358189e-01,
+  -3.823100e+00, -7.918398e-01,
+  -3.857021e+00, -8.727074e-01,
+  -3.858250e+00, -8.767645e-01,
+  -3.872100e+00, -9.563174e-01,
+  -3.864397e+00, -1.032630e+00,
+  -3.846230e+00, -1.081669e+00,
+  -3.834799e+00, -1.102536e+00,
+  -3.866684e+00, -1.022901e+00,
+  -3.808643e+00, -1.139084e+00,
+  -3.868840e+00, -1.011569e+00,
+  -3.791071e+00, -1.158615e+00,
+  -3.797999e+00, -1.151267e+00,
+  -3.696278e+00, -1.232314e+00,
+  -3.779007e+00, -1.170504e+00,
+  -3.622855e+00, -1.270793e+00,
+  -3.647249e+00, -1.259166e+00,
+  -3.655412e+00, -1.255042e+00,
+  -3.573218e+00, -1.291696e+00,
+  -3.638019e+00, -1.263684e+00,
+  -3.498409e+00, -1.317750e+00,
+  -3.304143e+00, -1.364970e+00,
+  -3.183001e+00, -1.384295e+00,
+  -3.202456e+00, -1.381599e+00,
+  -3.244063e+00, -1.375332e+00,
+  -3.233308e+00, -1.377019e+00,
+  -3.060112e+00, -1.398264e+00,
+  -3.078187e+00, -1.396517e+00,
+  -2.689594e+00, -1.415761e+00,
+  -2.947662e+00, -1.407039e+00,
+  -2.854490e+00, -1.411860e+00,
+  -2.660499e+00, -1.415900e+00,
+  -2.875955e+00, -1.410930e+00,
+  -2.675385e+00, -1.415848e+00,
+  -2.813155e+00, -1.413363e+00,
+  -2.417673e+00, -1.411512e+00,
+  -2.725461e+00, -1.415373e+00,
+  -2.148334e+00, -1.396672e+00,
+  -2.108972e+00, -1.393738e+00,
+  -2.029905e+00, -1.387302e+00,
+  -2.046214e+00, -1.388687e+00,
+  -2.057402e+00, -1.389621e+00,
+  -1.650250e+00, -1.347160e+00,
+  -1.806764e+00, -1.365469e+00,
+  -1.206973e+00, -1.283343e+00,
+  -8.029259e-01, -1.211308e+00,
+  -1.229551e+00, -1.286993e+00,
+  -1.101507e+00, -1.265754e+00,
+  -9.110645e-01, -1.231804e+00,
+  -1.110046e+00, -1.267211e+00,
+  -8.465274e-01, -1.219677e+00,
+  -7.594163e-01, -1.202818e+00,
+  -8.023823e-01, -1.211203e+00,
+  -3.732519e-01, -1.121494e+00,
+  -1.918373e-01, -1.079668e+00,
+  -4.671988e-01, -1.142253e+00,
+  -4.033645e-01, -1.128215e+00,
+  -1.920740e-01, -1.079724e+00,
+  -3.022157e-01, -1.105389e+00,
+  -1.652831e-01, -1.073354e+00,
+  +4.671625e-01, -9.085886e-01,
+  +5.940178e-01, -8.721832e-01,
+  +3.147557e-01, -9.508290e-01,
+  +6.383631e-01, -8.591867e-01,
+  +9.888923e-01, -7.514088e-01,
+  +7.076339e-01, -8.386023e-01,
+  +1.326682e+00, -6.386698e-01,
+  +1.149834e+00, -6.988221e-01,
+  +1.257742e+00, -6.624207e-01,
+  +1.492352e+00, -5.799632e-01,
+  +1.595574e+00, -5.421766e-01,
+  +1.240173e+00, -6.684113e-01,
+  +1.706612e+00, -5.004442e-01,
+  +1.873984e+00, -4.353002e-01,
+  +1.985633e+00, -3.902561e-01,
+  +1.722880e+00, -4.942329e-01,
+  +2.095182e+00, -3.447402e-01,
+  +2.018118e+00, -3.768991e-01,
+  +2.422702e+00, -1.999563e-01,
+  +2.370611e+00, -2.239326e-01,
+  +2.152154e+00, -3.205250e-01,
+  +2.525121e+00, -1.516499e-01,
+  +2.422116e+00, -2.002280e-01,
+  +2.842806e+00, +9.536372e-03,
+  +3.030128e+00, +1.146027e-01,
+  +2.888424e+00, +3.433444e-02,
+  +2.991609e+00, +9.226409e-02,
+  +2.924807e+00, +5.445844e-02,
+  +3.007772e+00, +1.015875e-01,
+  +2.781973e+00, -2.282382e-02,
+  +3.164737e+00, +1.961781e-01,
+  +3.237671e+00, +2.430139e-01,
+  +3.046123e+00, +1.240014e-01,
+  +3.414834e+00, +3.669060e-01,
+  +3.436591e+00, +3.833600e-01,
+  +3.626207e+00, +5.444311e-01,
+  +3.223325e+00, +2.336361e-01,
+  +3.511963e+00, +4.431060e-01,
+  +3.698380e+00, +6.187442e-01,
+  +3.670244e+00, +5.884943e-01,
+  +3.558833e+00, +4.828230e-01,
+  +3.661807e+00, +5.797689e-01,
+  +3.767261e+00, +7.030893e-01,
+  +3.801065e+00, +7.532650e-01,
+  +3.828523e+00, +8.024454e-01,
+  +3.840719e+00, +8.287032e-01,
+  +3.848748e+00, +8.485921e-01,
+  +3.865801e+00, +9.066551e-01,
+  +3.870983e+00, +9.404873e-01,
+  +3.870263e+00, +1.001884e+00,
+  +3.864462e+00, +1.032374e+00,
+  +3.870542e+00, +9.996121e-01,
+  +3.865424e+00, +1.028474e+00
+};
+ceres::ConstMatrixRef kY(kYData, kYRows, kYCols);
+
+class PointToLineSegmentContourCostFunction : public ceres::CostFunction {
+ public:
+  PointToLineSegmentContourCostFunction(const int num_segments,
+                                        const Eigen::Vector2d y)
+      : num_segments_(num_segments), y_(y) {
+    // The first parameter is the preimage position.
+    mutable_parameter_block_sizes()->push_back(1);
+    // The next parameters are the control points for the line segment contour.
+    for (int i = 0; i < num_segments_; ++i) {
+      mutable_parameter_block_sizes()->push_back(2);
+    }
+    set_num_residuals(2);
+  }
+
+  virtual bool Evaluate(const double* const* x,
+                        double* residuals,
+                        double** jacobians) const {
+    // Convert the preimage position `t` into a segment index `i0` and the
+    // line segment interpolation parameter `u`. `i1` is the index of the next
+    // control point.
+    const double t = ModuloNumSegments(*x[0]);
+    CHECK_GE(t, 0.0);
+    CHECK_LT(t, num_segments_);
+    const int i0 = floor(t), i1 = (i0 + 1) % num_segments_;
+    const double u = t - i0;
+
+    // Linearly interpolate between control points `i0` and `i1`.
+    residuals[0] = y_[0] - ((1.0 - u) * x[1 + i0][0] + u * x[1 + i1][0]);
+    residuals[1] = y_[1] - ((1.0 - u) * x[1 + i0][1] + u * x[1 + i1][1]);
+
+    if (jacobians == NULL) {
+      return true;
+    }
+
+    if (jacobians[0] != NULL) {
+      jacobians[0][0] = x[1 + i0][0] - x[1 + i1][0];
+      jacobians[0][1] = x[1 + i0][1] - x[1 + i1][1];
+    }
+    for (int i = 0; i < num_segments_; ++i) {
+      if (jacobians[i + 1] != NULL) {
+        ceres::MatrixRef(jacobians[i + 1], 2, 2).setZero();
+        if (i == i0) {
+          jacobians[i + 1][0] = -(1.0 - u);
+          jacobians[i + 1][3] = -(1.0 - u);
+        } else if (i == i1) {
+          jacobians[i + 1][0] = -u;
+          jacobians[i + 1][3] = -u;
+        }
+      }
+    }
+    return true;
+  }
+
+  static ceres::CostFunction* Create(const int num_segments,
+                                     const Eigen::Vector2d y) {
+    return new PointToLineSegmentContourCostFunction(num_segments, y);
+  }
+
+ private:
+  inline double ModuloNumSegments(const double& t) const {
+    return t - num_segments_ * floor(t / num_segments_);
+  }
+
+  const int num_segments_;
+  const Eigen::Vector2d y_;
+};
+
+struct EuclideanDistanceFunctor {
+  EuclideanDistanceFunctor(const double& sqrt_weight)
+      : sqrt_weight_(sqrt_weight) {}
+
+  template <typename T>
+  bool operator()(const T* x0, const T* x1, T* residuals) const {
+    residuals[0] = T(sqrt_weight_) * (x0[0] - x1[0]);
+    residuals[1] = T(sqrt_weight_) * (x0[1] - x1[1]);
+    return true;
+  }
+
+  static ceres::CostFunction* Create(const double& sqrt_weight) {
+    return new ceres::AutoDiffCostFunction<EuclideanDistanceFunctor, 2, 2, 2>(
+        new EuclideanDistanceFunctor(sqrt_weight));
+  }
+
+ private:
+  const double sqrt_weight_;
+};
+
+bool SolveWithFullReport(ceres::Solver::Options options,
+                         ceres::Problem* problem,
+                         bool dynamic_sparsity) {
+  options.dynamic_sparsity = dynamic_sparsity;
+
+  ceres::Solver::Summary summary;
+  ceres::Solve(options, problem, &summary);
+
+  std::cout << "####################" << std::endl;
+  std::cout << "dynamic_sparsity = " << dynamic_sparsity << std::endl;
+  std::cout << "####################" << std::endl;
+  std::cout << summary.FullReport() << std::endl;
+
+  return summary.termination_type == ceres::CONVERGENCE;
+}
+
+int main(int argc, char** argv) {
+  google::InitGoogleLogging(argv[0]);
+
+  // Problem configuration.
+  const int num_segments = 151;
+  const double regularization_weight = 1e-2;
+
+  // Eigen::MatrixXd is column major so we define our own MatrixXd which is
+  // row major. Eigen::VectorXd can be used directly.
+  typedef Eigen::Matrix<double,
+                        Eigen::Dynamic, Eigen::Dynamic,
+                        Eigen::RowMajor> MatrixXd;
+  using Eigen::VectorXd;
+
+  // `X` is the matrix of control points which make up the contour of line
+  // segments. The number of control points is equal to the number of line
+  // segments because the contour is closed.
+  //
+  // Initialize `X` to points on the unit circle.
+  VectorXd w(num_segments + 1);
+  w.setLinSpaced(num_segments + 1, 0.0, 2.0 * M_PI);
+  w.conservativeResize(num_segments);
+  MatrixXd X(num_segments, 2);
+  X.col(0) = w.array().cos();
+  X.col(1) = w.array().sin();
+
+  // Each data point has an associated preimage position on the line segment
+  // contour. For each data point we initialize the preimage positions to
+  // the index of the closest control point.
+  const int num_observations = kY.rows();
+  VectorXd t(num_observations);
+  for (int i = 0; i < num_observations; ++i) {
+    (X.rowwise() - kY.row(i)).rowwise().squaredNorm().minCoeff(&t[i]);
+  }
+
+  ceres::Problem problem;
+
+  // For each data point add a residual which measures its distance to its
+  // corresponding position on the line segment contour.
+  std::vector<double*> parameter_blocks(1 + num_segments);
+  parameter_blocks[0] = NULL;
+  for (int i = 0; i < num_segments; ++i) {
+    parameter_blocks[i + 1] = X.data() + 2 * i;
+  }
+  for (int i = 0; i < num_observations; ++i) {
+    parameter_blocks[0] = &t[i];
+    problem.AddResidualBlock(
+      PointToLineSegmentContourCostFunction::Create(num_segments, kY.row(i)),
+      NULL,
+      parameter_blocks);
+  }
+
+  // Add regularization to minimize the length of the line segment contour.
+  for (int i = 0; i < num_segments; ++i) {
+    problem.AddResidualBlock(
+      EuclideanDistanceFunctor::Create(sqrt(regularization_weight)),
+      NULL,
+      X.data() + 2 * i,
+      X.data() + 2 * ((i + 1) % num_segments));
+  }
+
+  ceres::Solver::Options options;
+  options.max_num_iterations = 100;
+  options.linear_solver_type = ceres::SPARSE_NORMAL_CHOLESKY;
+
+  // First, solve `X` and `t` jointly with dynamic_sparsity = true.
+  MatrixXd X0 = X;
+  VectorXd t0 = t;
+  CHECK(SolveWithFullReport(options, &problem, true));
+
+  // Second, solve with dynamic_sparsity = false.
+  X = X0;
+  t = t0;
+  CHECK(SolveWithFullReport(options, &problem, false));
+
+  return 0;
+}
diff --git a/examples/libmv_homography.cc b/examples/libmv_homography.cc
new file mode 100644
index 0000000..8bc7136
--- /dev/null
+++ b/examples/libmv_homography.cc
@@ -0,0 +1,414 @@
+// Ceres Solver - A fast non-linear least squares minimizer
+// Copyright 2014 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.
+//
+// Copyright (c) 2014 libmv authors.
+//
+// Permission is hereby granted, free of charge, to any person obtaining a copy
+// of this software and associated documentation files (the "Software"), to
+// deal in the Software without restriction, including without limitation the
+// rights to use, copy, modify, merge, publish, distribute, sublicense, and/or
+// sell copies of the Software, and to permit persons to whom the Software is
+// furnished to do so, subject to the following conditions:
+//
+// The above copyright notice and this permission notice shall be included in
+// all copies or substantial portions of the Software.
+//
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS
+// IN THE SOFTWARE.
+//
+// Author: sergey.vfx@gmail.com (Sergey Sharybin)
+//
+// This file demonstrates solving for a homography between two sets of points.
+// A homography describes a transformation between a sets of points on a plane,
+// perspectively projected into two images. The first step is to solve a
+// homogeneous system of equations via singular value decompposition, giving an
+// algebraic solution for the homography, then solving for a final solution by
+// minimizing the symmetric transfer error in image space with Ceres (called the
+// Gold Standard Solution in "Multiple View Geometry"). The routines are based on
+// the routines from the Libmv library.
+//
+// This example demonstrates custom exit criterion by having a callback check
+// for image-space error.
+
+#include "ceres/ceres.h"
+#include "glog/logging.h"
+
+typedef Eigen::NumTraits<double> EigenDouble;
+
+typedef Eigen::MatrixXd Mat;
+typedef Eigen::VectorXd Vec;
+typedef Eigen::Matrix<double, 3, 3> Mat3;
+typedef Eigen::Matrix<double, 2, 1> Vec2;
+typedef Eigen::Matrix<double, Eigen::Dynamic,  8> MatX8;
+typedef Eigen::Vector3d Vec3;
+
+namespace {
+
+// This structure contains options that controls how the homography
+// estimation operates.
+//
+// Defaults should be suitable for a wide range of use cases, but
+// better performance and accuracy might require tweaking.
+struct EstimateHomographyOptions {
+  // Default settings for homography estimation which should be suitable
+  // for a wide range of use cases.
+  EstimateHomographyOptions()
+    :  max_num_iterations(50),
+       expected_average_symmetric_distance(1e-16) {}
+
+  // Maximal number of iterations for the refinement step.
+  int max_num_iterations;
+
+  // Expected average of symmetric geometric distance between
+  // actual destination points and original ones transformed by
+  // estimated homography matrix.
+  //
+  // Refinement will finish as soon as average of symmetric
+  // geometric distance is less or equal to this value.
+  //
+  // This distance is measured in the same units as input points are.
+  double expected_average_symmetric_distance;
+};
+
+// Calculate symmetric geometric cost terms:
+//
+// forward_error = D(H * x1, x2)
+// backward_error = D(H^-1 * x2, x1)
+//
+// Templated to be used with autodifferenciation.
+template <typename T>
+void SymmetricGeometricDistanceTerms(const Eigen::Matrix<T, 3, 3> &H,
+                                     const Eigen::Matrix<T, 2, 1> &x1,
+                                     const Eigen::Matrix<T, 2, 1> &x2,
+                                     T forward_error[2],
+                                     T backward_error[2]) {
+  typedef Eigen::Matrix<T, 3, 1> Vec3;
+  Vec3 x(x1(0), x1(1), T(1.0));
+  Vec3 y(x2(0), x2(1), T(1.0));
+
+  Vec3 H_x = H * x;
+  Vec3 Hinv_y = H.inverse() * y;
+
+  H_x /= H_x(2);
+  Hinv_y /= Hinv_y(2);
+
+  forward_error[0] = H_x(0) - y(0);
+  forward_error[1] = H_x(1) - y(1);
+  backward_error[0] = Hinv_y(0) - x(0);
+  backward_error[1] = Hinv_y(1) - x(1);
+}
+
+// Calculate symmetric geometric cost:
+//
+//   D(H * x1, x2)^2 + D(H^-1 * x2, x1)^2
+//
+double SymmetricGeometricDistance(const Mat3 &H,
+                                  const Vec2 &x1,
+                                  const Vec2 &x2) {
+  Vec2 forward_error, backward_error;
+  SymmetricGeometricDistanceTerms<double>(H,
+                                          x1,
+                                          x2,
+                                          forward_error.data(),
+                                          backward_error.data());
+  return forward_error.squaredNorm() +
+         backward_error.squaredNorm();
+}
+
+// A parameterization of the 2D homography matrix that uses 8 parameters so
+// that the matrix is normalized (H(2,2) == 1).
+// The homography matrix H is built from a list of 8 parameters (a, b,...g, h)
+// as follows
+//
+//         |a b c|
+//     H = |d e f|
+//         |g h 1|
+//
+template<typename T = double>
+class Homography2DNormalizedParameterization {
+ public:
+  typedef Eigen::Matrix<T, 8, 1> Parameters;     // a, b, ... g, h
+  typedef Eigen::Matrix<T, 3, 3> Parameterized;  // H
+
+  // Convert from the 8 parameters to a H matrix.
+  static void To(const Parameters &p, Parameterized *h) {
+    *h << p(0), p(1), p(2),
+          p(3), p(4), p(5),
+          p(6), p(7), 1.0;
+  }
+
+  // Convert from a H matrix to the 8 parameters.
+  static void From(const Parameterized &h, Parameters *p) {
+    *p << h(0, 0), h(0, 1), h(0, 2),
+          h(1, 0), h(1, 1), h(1, 2),
+          h(2, 0), h(2, 1);
+  }
+};
+
+// 2D Homography transformation estimation in the case that points are in
+// euclidean coordinates.
+//
+//   x = H y
+//
+// x and y vector must have the same direction, we could write
+//
+//   crossproduct(|x|, * H * |y| ) = |0|
+//
+//   | 0 -1  x2|   |a b c|   |y1|    |0|
+//   | 1  0 -x1| * |d e f| * |y2| =  |0|
+//   |-x2  x1 0|   |g h 1|   |1 |    |0|
+//
+// That gives:
+//
+//   (-d+x2*g)*y1    + (-e+x2*h)*y2 + -f+x2          |0|
+//   (a-x1*g)*y1     + (b-x1*h)*y2  + c-x1         = |0|
+//   (-x2*a+x1*d)*y1 + (-x2*b+x1*e)*y2 + -x2*c+x1*f  |0|
+//
+bool Homography2DFromCorrespondencesLinearEuc(
+    const Mat &x1,
+    const Mat &x2,
+    Mat3 *H,
+    double expected_precision) {
+  assert(2 == x1.rows());
+  assert(4 <= x1.cols());
+  assert(x1.rows() == x2.rows());
+  assert(x1.cols() == x2.cols());
+
+  int n = x1.cols();
+  MatX8 L = Mat::Zero(n * 3, 8);
+  Mat b = Mat::Zero(n * 3, 1);
+  for (int i = 0; i < n; ++i) {
+    int j = 3 * i;
+    L(j, 0) =  x1(0, i);             // a
+    L(j, 1) =  x1(1, i);             // b
+    L(j, 2) =  1.0;                  // c
+    L(j, 6) = -x2(0, i) * x1(0, i);  // g
+    L(j, 7) = -x2(0, i) * x1(1, i);  // h
+    b(j, 0) =  x2(0, i);             // i
+
+    ++j;
+    L(j, 3) =  x1(0, i);             // d
+    L(j, 4) =  x1(1, i);             // e
+    L(j, 5) =  1.0;                  // f
+    L(j, 6) = -x2(1, i) * x1(0, i);  // g
+    L(j, 7) = -x2(1, i) * x1(1, i);  // h
+    b(j, 0) =  x2(1, i);             // i
+
+    // This ensures better stability
+    // TODO(julien) make a lite version without this 3rd set
+    ++j;
+    L(j, 0) =  x2(1, i) * x1(0, i);  // a
+    L(j, 1) =  x2(1, i) * x1(1, i);  // b
+    L(j, 2) =  x2(1, i);             // c
+    L(j, 3) = -x2(0, i) * x1(0, i);  // d
+    L(j, 4) = -x2(0, i) * x1(1, i);  // e
+    L(j, 5) = -x2(0, i);             // f
+  }
+  // Solve Lx=B
+  const Vec h = L.fullPivLu().solve(b);
+  Homography2DNormalizedParameterization<double>::To(h, H);
+  return (L * h).isApprox(b, expected_precision);
+}
+
+// Cost functor which computes symmetric geometric distance
+// used for homography matrix refinement.
+class HomographySymmetricGeometricCostFunctor {
+ public:
+  HomographySymmetricGeometricCostFunctor(const Vec2 &x,
+                                          const Vec2 &y)
+      : x_(x), y_(y) { }
+
+  template<typename T>
+  bool operator()(const T* homography_parameters, T* residuals) const {
+    typedef Eigen::Matrix<T, 3, 3> Mat3;
+    typedef Eigen::Matrix<T, 2, 1> Vec2;
+
+    Mat3 H(homography_parameters);
+    Vec2 x(T(x_(0)), T(x_(1)));
+    Vec2 y(T(y_(0)), T(y_(1)));
+
+    SymmetricGeometricDistanceTerms<T>(H,
+                                       x,
+                                       y,
+                                       &residuals[0],
+                                       &residuals[2]);
+    return true;
+  }
+
+  const Vec2 x_;
+  const Vec2 y_;
+};
+
+// Termination checking callback. This is needed to finish the
+// optimization when an absolute error threshold is met, as opposed
+// to Ceres's function_tolerance, which provides for finishing when
+// successful steps reduce the cost function by a fractional amount.
+// In this case, the callback checks for the absolute average reprojection
+// error and terminates when it's below a threshold (for example all
+// points < 0.5px error).
+class TerminationCheckingCallback : public ceres::IterationCallback {
+ public:
+  TerminationCheckingCallback(const Mat &x1, const Mat &x2,
+                              const EstimateHomographyOptions &options,
+                              Mat3 *H)
+      : options_(options), x1_(x1), x2_(x2), H_(H) {}
+
+  virtual ceres::CallbackReturnType operator()(
+      const ceres::IterationSummary& summary) {
+    // If the step wasn't successful, there's nothing to do.
+    if (!summary.step_is_successful) {
+      return ceres::SOLVER_CONTINUE;
+    }
+
+    // Calculate average of symmetric geometric distance.
+    double average_distance = 0.0;
+    for (int i = 0; i < x1_.cols(); i++) {
+      average_distance += SymmetricGeometricDistance(*H_,
+                                                     x1_.col(i),
+                                                     x2_.col(i));
+    }
+    average_distance /= x1_.cols();
+
+    if (average_distance <= options_.expected_average_symmetric_distance) {
+      return ceres::SOLVER_TERMINATE_SUCCESSFULLY;
+    }
+
+    return ceres::SOLVER_CONTINUE;
+  }
+
+ private:
+  const EstimateHomographyOptions &options_;
+  const Mat &x1_;
+  const Mat &x2_;
+  Mat3 *H_;
+};
+
+bool EstimateHomography2DFromCorrespondences(
+    const Mat &x1,
+    const Mat &x2,
+    const EstimateHomographyOptions &options,
+    Mat3 *H) {
+  assert(2 == x1.rows());
+  assert(4 <= x1.cols());
+  assert(x1.rows() == x2.rows());
+  assert(x1.cols() == x2.cols());
+
+  // Step 1: Algebraic homography estimation.
+  // Assume algebraic estimation always succeeds.
+  Homography2DFromCorrespondencesLinearEuc(x1,
+                                           x2,
+                                           H,
+                                           EigenDouble::dummy_precision());
+
+  LOG(INFO) << "Estimated matrix after algebraic estimation:\n" << *H;
+
+  // Step 2: Refine matrix using Ceres minimizer.
+  ceres::Problem problem;
+  for (int i = 0; i < x1.cols(); i++) {
+    HomographySymmetricGeometricCostFunctor
+        *homography_symmetric_geometric_cost_function =
+            new HomographySymmetricGeometricCostFunctor(x1.col(i),
+                                                        x2.col(i));
+
+    problem.AddResidualBlock(
+        new ceres::AutoDiffCostFunction<
+            HomographySymmetricGeometricCostFunctor,
+            4,  // num_residuals
+            9>(homography_symmetric_geometric_cost_function),
+        NULL,
+        H->data());
+  }
+
+  // Configure the solve.
+  ceres::Solver::Options solver_options;
+  solver_options.linear_solver_type = ceres::DENSE_QR;
+  solver_options.max_num_iterations = options.max_num_iterations;
+  solver_options.update_state_every_iteration = true;
+
+  // Terminate if the average symmetric distance is good enough.
+  TerminationCheckingCallback callback(x1, x2, options, H);
+  solver_options.callbacks.push_back(&callback);
+
+  // Run the solve.
+  ceres::Solver::Summary summary;
+  ceres::Solve(solver_options, &problem, &summary);
+
+  LOG(INFO) << "Summary:\n" << summary.FullReport();
+  LOG(INFO) << "Final refined matrix:\n" << *H;
+
+  return summary.IsSolutionUsable();
+}
+
+}  // namespace libmv
+
+int main(int argc, char **argv) {
+  google::InitGoogleLogging(argv[0]);
+
+  Mat x1(2, 100);
+  for (int i = 0; i < x1.cols(); ++i) {
+    x1(0, i) = rand() % 1024;
+    x1(1, i) = rand() % 1024;
+  }
+
+  Mat3 homography_matrix;
+  // This matrix has been dumped from a Blender test file of plane tracking.
+  homography_matrix << 1.243715, -0.461057, -111.964454,
+                       0.0,       0.617589, -192.379252,
+                       0.0,      -0.000983,    1.0;
+
+  Mat x2 = x1;
+  for (int i = 0; i < x2.cols(); ++i) {
+    Vec3 homogenous_x1 = Vec3(x1(0, i), x1(1, i), 1.0);
+    Vec3 homogenous_x2 = homography_matrix * homogenous_x1;
+    x2(0, i) = homogenous_x2(0) / homogenous_x2(2);
+    x2(1, i) = homogenous_x2(1) / homogenous_x2(2);
+
+    // Apply some noise so algebraic estimation is not good enough.
+    x2(0, i) += static_cast<double>(rand() % 1000) / 5000.0;
+    x2(1, i) += static_cast<double>(rand() % 1000) / 5000.0;
+  }
+
+  Mat3 estimated_matrix;
+
+  EstimateHomographyOptions options;
+  options.expected_average_symmetric_distance = 0.02;
+  EstimateHomography2DFromCorrespondences(x1, x2, options, &estimated_matrix);
+
+  // Normalize the matrix for easier comparison.
+  estimated_matrix /= estimated_matrix(2 ,2);
+
+  std::cout << "Original matrix:\n" << homography_matrix << "\n";
+  std::cout << "Estimated matrix:\n" << estimated_matrix << "\n";
+
+  return EXIT_SUCCESS;
+}
diff --git a/examples/more_garbow_hillstrom.cc b/examples/more_garbow_hillstrom.cc
new file mode 100644
index 0000000..d98e57c
--- /dev/null
+++ b/examples/more_garbow_hillstrom.cc
@@ -0,0 +1,374 @@
+// Ceres Solver - A fast non-linear least squares minimizer
+// Copyright 2014 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)
+//
+// Bounds constrained test problems from the paper
+//
+// Testing Unconstrained Optimization Software
+// Jorge J. More, Burton S. Garbow and Kenneth E. Hillstrom
+// ACM Transactions on Mathematical Software, 7(1), pp. 17-41, 1981
+//
+// A subset of these problems were augmented with bounds and used for
+// testing bounds constrained optimization algorithms by
+//
+// A Trust Region Approach to Linearly Constrained Optimization
+// David M. Gay
+// Numerical Analysis (Griffiths, D.F., ed.), pp. 72-105
+// Lecture Notes in Mathematics 1066, Springer Verlag, 1984.
+//
+// The latter paper is behind a paywall. We obtained the bounds on the
+// variables and the function values at the global minimums from
+//
+// http://www.mat.univie.ac.at/~neum/glopt/bounds.html
+//
+// A problem is considered solved if of the log relative error of its
+// objective function is at least 5.
+
+
+#include <cmath>
+#include <iostream>  // NOLINT
+#include "ceres/ceres.h"
+#include "gflags/gflags.h"
+#include "glog/logging.h"
+
+namespace ceres {
+namespace examples {
+
+const double kDoubleMax = std::numeric_limits<double>::max();
+
+#define BEGIN_MGH_PROBLEM(name, num_parameters, num_residuals)          \
+  struct name {                                                         \
+    static const int kNumParameters = num_parameters;                   \
+    static const double initial_x[kNumParameters];                      \
+    static const double lower_bounds[kNumParameters];                   \
+    static const double upper_bounds[kNumParameters];                   \
+    static const double constrained_optimal_cost;                       \
+    static const double unconstrained_optimal_cost;                     \
+    static CostFunction* Create() {                                     \
+      return new AutoDiffCostFunction<name,                             \
+                                      num_residuals,                    \
+                                      num_parameters>(new name);        \
+    }                                                                   \
+    template <typename T>                                               \
+    bool operator()(const T* const x, T* residual) const {
+
+#define END_MGH_PROBLEM return true; } };  // NOLINT
+
+// Rosenbrock function.
+BEGIN_MGH_PROBLEM(TestProblem1, 2, 2)
+  const T x1 = x[0];
+  const T x2 = x[1];
+  residual[0] = T(10.0) * (x2 - x1 * x1);
+  residual[1] = T(1.0) - x1;
+END_MGH_PROBLEM;
+
+const double TestProblem1::initial_x[] = {-1.2, 1.0};
+const double TestProblem1::lower_bounds[] = {-kDoubleMax, -kDoubleMax};
+const double TestProblem1::upper_bounds[] = {kDoubleMax, kDoubleMax};
+const double TestProblem1::constrained_optimal_cost =
+    std::numeric_limits<double>::quiet_NaN();
+const double TestProblem1::unconstrained_optimal_cost = 0.0;
+
+// Freudenstein and Roth function.
+BEGIN_MGH_PROBLEM(TestProblem2, 2, 2)
+  const T x1 = x[0];
+  const T x2 = x[1];
+  residual[0] = T(-13.0) + x1 + ((T(5.0) - x2) * x2 - T(2.0)) * x2;
+  residual[1] = T(-29.0) + x1 + ((x2 + T(1.0)) * x2 - T(14.0)) * x2;
+END_MGH_PROBLEM;
+
+const double TestProblem2::initial_x[] = {0.5, -2.0};
+const double TestProblem2::lower_bounds[] = {-kDoubleMax, -kDoubleMax};
+const double TestProblem2::upper_bounds[] = {kDoubleMax, kDoubleMax};
+const double TestProblem2::constrained_optimal_cost =
+    std::numeric_limits<double>::quiet_NaN();
+const double TestProblem2::unconstrained_optimal_cost = 0.0;
+
+// Powell badly scaled function.
+BEGIN_MGH_PROBLEM(TestProblem3, 2, 2)
+  const T x1 = x[0];
+  const T x2 = x[1];
+  residual[0] = T(10000.0) * x1 * x2 - T(1.0);
+  residual[1] = exp(-x1) + exp(-x2) - T(1.0001);
+END_MGH_PROBLEM;
+
+const double TestProblem3::initial_x[] = {0.0, 1.0};
+const double TestProblem3::lower_bounds[] = {0.0, 1.0};
+const double TestProblem3::upper_bounds[] = {1.0, 9.0};
+const double TestProblem3::constrained_optimal_cost = 0.15125900e-9;
+const double TestProblem3::unconstrained_optimal_cost = 0.0;
+
+// Brown badly scaled function.
+BEGIN_MGH_PROBLEM(TestProblem4, 2, 3)
+  const T x1 = x[0];
+  const T x2 = x[1];
+  residual[0] = x1  - T(1000000.0);
+  residual[1] = x2 - T(0.000002);
+  residual[2] = x1 * x2 - T(2.0);
+END_MGH_PROBLEM;
+
+const double TestProblem4::initial_x[] = {1.0, 1.0};
+const double TestProblem4::lower_bounds[] = {0.0, 0.00003};
+const double TestProblem4::upper_bounds[] = {1000000.0, 100.0};
+const double TestProblem4::constrained_optimal_cost = 0.78400000e3;
+const double TestProblem4::unconstrained_optimal_cost = 0.0;
+
+// Beale function.
+BEGIN_MGH_PROBLEM(TestProblem5, 2, 3)
+  const T x1 = x[0];
+  const T x2 = x[1];
+  residual[0] = T(1.5) - x1 * (T(1.0) - x2);
+  residual[1] = T(2.25) - x1 * (T(1.0) - x2 * x2);
+  residual[2] = T(2.625) - x1 * (T(1.0) - x2 * x2 * x2);
+END_MGH_PROBLEM;
+
+const double TestProblem5::initial_x[] = {1.0, 1.0};
+const double TestProblem5::lower_bounds[] = {0.6, 0.5};
+const double TestProblem5::upper_bounds[] = {10.0, 100.0};
+const double TestProblem5::constrained_optimal_cost = 0.0;
+const double TestProblem5::unconstrained_optimal_cost = 0.0;
+
+// Jennrich and Sampson function.
+BEGIN_MGH_PROBLEM(TestProblem6, 2, 10)
+  const T x1 = x[0];
+  const T x2 = x[1];
+  for (int i = 1; i <= 10; ++i) {
+    residual[i - 1] = T(2.0) + T(2.0 * i) -
+        exp(T(static_cast<double>(i)) * x1) -
+        exp(T(static_cast<double>(i) * x2));
+  }
+END_MGH_PROBLEM;
+
+const double TestProblem6::initial_x[] = {1.0, 1.0};
+const double TestProblem6::lower_bounds[] = {-kDoubleMax, -kDoubleMax};
+const double TestProblem6::upper_bounds[] = {kDoubleMax, kDoubleMax};
+const double TestProblem6::constrained_optimal_cost =
+    std::numeric_limits<double>::quiet_NaN();
+const double TestProblem6::unconstrained_optimal_cost = 124.362;
+
+// Helical valley function.
+BEGIN_MGH_PROBLEM(TestProblem7, 3, 3)
+  const T x1 = x[0];
+  const T x2 = x[1];
+  const T x3 = x[2];
+  const T theta = T(0.5 / M_PI)  * atan(x2 / x1) + (x1 > 0.0 ? T(0.0) : T(0.5));
+
+  residual[0] = T(10.0) * (x3 - T(10.0) * theta);
+  residual[1] = T(10.0) * (sqrt(x1 * x1 + x2 * x2) - T(1.0));
+  residual[2] = x3;
+END_MGH_PROBLEM;
+
+const double TestProblem7::initial_x[] = {-1.0, 0.0, 0.0};
+const double TestProblem7::lower_bounds[] = {-100.0, -1.0, -1.0};
+const double TestProblem7::upper_bounds[] = {0.8, 1.0, 1.0};
+const double TestProblem7::constrained_optimal_cost = 0.99042212;
+const double TestProblem7::unconstrained_optimal_cost = 0.0;
+
+// Bard function
+BEGIN_MGH_PROBLEM(TestProblem8, 3, 15)
+  const T x1 = x[0];
+  const T x2 = x[1];
+  const T x3 = x[2];
+
+  double y[] = {0.14, 0.18, 0.22, 0.25,
+                0.29, 0.32, 0.35, 0.39, 0.37, 0.58,
+                0.73, 0.96, 1.34, 2.10, 4.39};
+
+  for (int i = 1; i <=15; ++i) {
+    const T u = T(static_cast<double>(i));
+    const T v = T(static_cast<double>(16 - i));
+    const T w = T(static_cast<double>(std::min(i, 16 - i)));
+    residual[i - 1] = T(y[i - 1]) - x1 + u / (v * x2 + w * x3);
+  }
+END_MGH_PROBLEM;
+
+const double TestProblem8::initial_x[] = {1.0, 1.0, 1.0};
+const double TestProblem8::lower_bounds[] = {
+  -kDoubleMax, -kDoubleMax, -kDoubleMax};
+const double TestProblem8::upper_bounds[] = {
+  kDoubleMax, kDoubleMax, kDoubleMax};
+const double TestProblem8::constrained_optimal_cost =
+    std::numeric_limits<double>::quiet_NaN();
+const double TestProblem8::unconstrained_optimal_cost = 8.21487e-3;
+
+// Gaussian function.
+BEGIN_MGH_PROBLEM(TestProblem9, 3, 15)
+  const T x1 = x[0];
+  const T x2 = x[1];
+  const T x3 = x[2];
+
+  const double y[] = {0.0009, 0.0044, 0.0175, 0.0540, 0.1295, 0.2420, 0.3521,
+                      0.3989,
+                      0.3521, 0.2420, 0.1295, 0.0540, 0.0175, 0.0044, 0.0009};
+  for (int i = 0; i < 15; ++i) {
+    const T t_i = T((8.0 - i - 1.0) / 2.0);
+    const T y_i = T(y[i]);
+    residual[i] = x1 * exp(-x2 * (t_i - x3) * (t_i - x3) / T(2.0)) - y_i;
+  }
+END_MGH_PROBLEM;
+
+const double TestProblem9::initial_x[] = {0.4, 1.0, 0.0};
+const double TestProblem9::lower_bounds[] = {0.398, 1.0, -0.5};
+const double TestProblem9::upper_bounds[] = {4.2, 2.0, 0.1};
+const double TestProblem9::constrained_optimal_cost = 0.11279300e-7;
+const double TestProblem9::unconstrained_optimal_cost = 0.112793e-7;
+
+// Meyer function.
+BEGIN_MGH_PROBLEM(TestProblem10, 3, 16)
+  const T x1 = x[0];
+  const T x2 = x[1];
+  const T x3 = x[2];
+
+  const double y[] = {34780, 28610, 23650, 19630, 16370, 13720, 11540, 9744,
+                      8261, 7030, 6005, 5147, 4427, 3820, 3307, 2872};
+
+  for (int i = 0; i < 16; ++i) {
+    T t = T(45 + 5.0 * (i + 1));
+    residual[i] = x1 * exp(x2 / (t + x3)) - y[i];
+  }
+END_MGH_PROBLEM
+
+
+const double TestProblem10::initial_x[] = {0.02, 4000, 250};
+const double TestProblem10::lower_bounds[] ={
+  -kDoubleMax, -kDoubleMax, -kDoubleMax};
+const double TestProblem10::upper_bounds[] ={
+  kDoubleMax, kDoubleMax, kDoubleMax};
+const double TestProblem10::constrained_optimal_cost =
+    std::numeric_limits<double>::quiet_NaN();
+const double TestProblem10::unconstrained_optimal_cost = 87.9458;
+
+#undef BEGIN_MGH_PROBLEM
+#undef END_MGH_PROBLEM
+
+template<typename TestProblem> string ConstrainedSolve() {
+  double x[TestProblem::kNumParameters];
+  std::copy(TestProblem::initial_x,
+            TestProblem::initial_x + TestProblem::kNumParameters,
+            x);
+
+  Problem problem;
+  problem.AddResidualBlock(TestProblem::Create(), NULL, x);
+  for (int i = 0; i < TestProblem::kNumParameters; ++i) {
+    problem.SetParameterLowerBound(x, i, TestProblem::lower_bounds[i]);
+    problem.SetParameterUpperBound(x, i, TestProblem::upper_bounds[i]);
+  }
+
+  Solver::Options options;
+  options.parameter_tolerance = 1e-18;
+  options.function_tolerance = 1e-18;
+  options.gradient_tolerance = 1e-18;
+  options.max_num_iterations = 1000;
+  options.linear_solver_type = DENSE_QR;
+  Solver::Summary summary;
+  Solve(options, &problem, &summary);
+
+  const double kMinLogRelativeError = 5.0;
+  const double log_relative_error = -std::log10(
+      std::abs(2.0 * summary.final_cost -
+               TestProblem::constrained_optimal_cost) /
+      (TestProblem::constrained_optimal_cost > 0.0
+       ? TestProblem::constrained_optimal_cost
+       : 1.0));
+
+  return (log_relative_error >= kMinLogRelativeError
+          ? "Success\n"
+          : "Failure\n");
+}
+
+template<typename TestProblem> string UnconstrainedSolve() {
+  double x[TestProblem::kNumParameters];
+  std::copy(TestProblem::initial_x,
+            TestProblem::initial_x + TestProblem::kNumParameters,
+            x);
+
+  Problem problem;
+  problem.AddResidualBlock(TestProblem::Create(), NULL, x);
+
+  Solver::Options options;
+  options.parameter_tolerance = 1e-18;
+  options.function_tolerance = 0.0;
+  options.gradient_tolerance = 1e-18;
+  options.max_num_iterations = 1000;
+  options.linear_solver_type = DENSE_QR;
+  Solver::Summary summary;
+  Solve(options, &problem, &summary);
+
+  const double kMinLogRelativeError = 5.0;
+  const double log_relative_error = -std::log10(
+      std::abs(2.0 * summary.final_cost -
+               TestProblem::unconstrained_optimal_cost) /
+      (TestProblem::unconstrained_optimal_cost > 0.0
+       ? TestProblem::unconstrained_optimal_cost
+       : 1.0));
+
+  return (log_relative_error >= kMinLogRelativeError
+          ? "Success\n"
+          : "Failure\n");
+}
+
+}  // namespace examples
+}  // namespace ceres
+
+int main(int argc, char** argv) {
+  google::ParseCommandLineFlags(&argc, &argv, true);
+  google::InitGoogleLogging(argv[0]);
+
+  using ceres::examples::UnconstrainedSolve;
+  using ceres::examples::ConstrainedSolve;
+
+#define UNCONSTRAINED_SOLVE(n)                                          \
+  std::cout << "Problem " << n << " : "                                 \
+            << UnconstrainedSolve<ceres::examples::TestProblem##n>();
+
+#define CONSTRAINED_SOLVE(n)                                            \
+  std::cout << "Problem " << n << " : "                                 \
+            << ConstrainedSolve<ceres::examples::TestProblem##n>();
+
+  std::cout << "Unconstrained problems\n";
+  UNCONSTRAINED_SOLVE(1);
+  UNCONSTRAINED_SOLVE(2);
+  UNCONSTRAINED_SOLVE(3);
+  UNCONSTRAINED_SOLVE(4);
+  UNCONSTRAINED_SOLVE(5);
+  UNCONSTRAINED_SOLVE(6);
+  UNCONSTRAINED_SOLVE(7);
+  UNCONSTRAINED_SOLVE(8);
+  UNCONSTRAINED_SOLVE(9);
+  UNCONSTRAINED_SOLVE(10);
+
+  std::cout << "\nConstrained problems\n";
+  CONSTRAINED_SOLVE(3);
+  CONSTRAINED_SOLVE(4);
+  CONSTRAINED_SOLVE(5);
+  CONSTRAINED_SOLVE(7);
+  CONSTRAINED_SOLVE(9);
+
+  return 0;
+}
diff --git a/examples/nist.cc b/examples/nist.cc
index 1773a0f..b29b285 100644
--- a/examples/nist.cc
+++ b/examples/nist.cc
@@ -159,14 +159,6 @@ void SkipLines(std::ifstream& ifs, int num_lines) {
   }
 }
 
-bool IsSuccessfulTermination(ceres::SolverTerminationType status) {
-  return
-      (status == ceres::FUNCTION_TOLERANCE) ||
-      (status == ceres::GRADIENT_TOLERANCE) ||
-      (status == ceres::PARAMETER_TOLERANCE) ||
-      (status == ceres::USER_SUCCESS);
-}
-
 class NISTProblem {
  public:
   explicit NISTProblem(const std::string& filename) {
diff --git a/examples/pgm_image.h b/examples/pgm_image.h
index 15e99e4..1328d75 100644
--- a/examples/pgm_image.h
+++ b/examples/pgm_image.h
@@ -197,7 +197,7 @@ bool PGMImage<Real>::WriteToFile(std::string filename) const {
     outputfile << static_cast<int>(data_[i] + 0.5) << ' ';
   }
 
-  return outputfile;  // Returns true/false
+  return bool(outputfile);  // Returns true/false
 }
 
 namespace  {
diff --git a/examples/quadratic.cc b/examples/quadratic.cc
deleted file mode 100644
index 8527af3..0000000
--- a/examples/quadratic.cc
+++ /dev/null
@@ -1,90 +0,0 @@
-// Ceres Solver - A fast non-linear least squares minimizer
-// Copyright 2010, 2011, 2012 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: keir@google.com (Keir Mierle)
-//
-// A simple example of using the Ceres minimizer.
-//
-// Minimize 0.5 (10 - x)^2 using analytic jacobian matrix.
-
-#include <vector>
-#include "ceres/ceres.h"
-#include "gflags/gflags.h"
-#include "glog/logging.h"
-
-using ceres::SizedCostFunction;
-using ceres::Problem;
-using ceres::Solver;
-using ceres::Solve;
-
-class SimpleCostFunction
-  : public SizedCostFunction<1 /* number of residuals */,
-                             1 /* size of first parameter */> {
- public:
-  virtual ~SimpleCostFunction() {}
-  virtual bool Evaluate(double const* const* parameters,
-                        double* residuals,
-                        double** jacobians) const {
-    double x = parameters[0][0];
-
-    // f(x) = 10 - x.
-    residuals[0] = 10 - x;
-
-    // f'(x) = -1. Since there's only 1 parameter and that parameter
-    // has 1 dimension, there is only 1 element to fill in the
-    // jacobians.
-    if (jacobians != NULL && jacobians[0] != NULL) {
-      jacobians[0][0] = -1;
-    }
-    return true;
-  }
-};
-
-int main(int argc, char** argv) {
-  google::ParseCommandLineFlags(&argc, &argv, true);
-  google::InitGoogleLogging(argv[0]);
-
-  // The variable with its initial value that we will be solving for.
-  double x = 5.0;
-
-  // Build the problem.
-  Problem problem;
-  // Set up the only cost function (also known as residual).
-  problem.AddResidualBlock(new SimpleCostFunction, NULL, &x);
-
-  // Run the solver!
-  Solver::Options options;
-  options.max_num_iterations = 10;
-  options.linear_solver_type = ceres::DENSE_QR;
-  options.minimizer_progress_to_stdout = true;
-  Solver::Summary summary;
-  Solve(options, &problem, &summary);
-  std::cout << summary.BriefReport() << "\n";
-  std::cout << "x : 5.0 -> " << x << "\n";
-  return 0;
-}
diff --git a/examples/quadratic_auto_diff.cc b/examples/quadratic_auto_diff.cc
deleted file mode 100644
index 1e2f3ef..0000000
--- a/examples/quadratic_auto_diff.cc
+++ /dev/null
@@ -1,89 +0,0 @@
-// Ceres Solver - A fast non-linear least squares minimizer
-// Copyright 2010, 2011, 2012 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: keir@google.com (Keir Mierle)
-//
-// A simple example of using the Ceres minimizer.
-//
-// Minimize 0.5 (10 - x)^2 using jacobian matrix computed using
-// automatic differentiation.
-
-#include <vector>
-#include "ceres/ceres.h"
-#include "gflags/gflags.h"
-#include "glog/logging.h"
-
-using ceres::AutoDiffCostFunction;
-using ceres::CostFunction;
-using ceres::Problem;
-using ceres::Solver;
-using ceres::Solve;
-
-// A templated cost functor that implements the residual r = 10 -
-// x. The method operator() is templated so that we can then use an
-// automatic differentiation wrapper around it to generate its
-// derivatives.
-class QuadraticCostFunctor {
- public:
-  template <typename T> bool operator()(const T* const x, T* residual) const {
-    residual[0] = T(10.0) - x[0];
-    return true;
-  }
-};
-
-int main(int argc, char** argv) {
-  google::ParseCommandLineFlags(&argc, &argv, true);
-  google::InitGoogleLogging(argv[0]);
-
-  // The variable to solve for with its initial value.
-  double initial_x = 5.0;
-  double x = initial_x;
-
-  // Build the problem.
-  Problem problem;
-
-  // Set up the only cost function (also known as residual). This uses
-  // auto-differentiation to obtain the derivative (jacobian).
-  problem.AddResidualBlock(
-      new AutoDiffCostFunction<QuadraticCostFunctor, 1, 1>(
-          new QuadraticCostFunctor),
-      NULL,
-      &x);
-
-  // Run the solver!
-  Solver::Options options;
-  options.max_num_iterations = 10;
-  options.linear_solver_type = ceres::DENSE_QR;
-  options.minimizer_progress_to_stdout = true;
-  Solver::Summary summary;
-  Solve(options, &problem, &summary);
-  std::cout << summary.BriefReport() << "\n";
-  std::cout << "x : " << initial_x
-            << " -> " << x << "\n";
-  return 0;
-}
diff --git a/examples/quadratic_numeric_diff.cc b/examples/quadratic_numeric_diff.cc
deleted file mode 100644
index 1082616..0000000
--- a/examples/quadratic_numeric_diff.cc
+++ /dev/null
@@ -1,84 +0,0 @@
-// Ceres Solver - A fast non-linear least squares minimizer
-// Copyright 2010, 2011, 2012 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: keir@google.com (Keir Mierle)
-//
-// Minimize 0.5 (10 - x)^2 using jacobian matrix computed using
-// numeric differentiation.
-
-#include <vector>
-#include "ceres/ceres.h"
-#include "gflags/gflags.h"
-#include "glog/logging.h"
-
-using ceres::NumericDiffCostFunction;
-using ceres::CENTRAL;
-using ceres::CostFunction;
-using ceres::Problem;
-using ceres::Solver;
-using ceres::Solve;
-
-// A cost functor that implements the residual r = 10 - x.
-class QuadraticCostFunctor {
- public:
-  bool operator()(const double* const x, double* residual) const {
-    residual[0] = 10.0 - x[0];
-    return true;
-  }
-};
-
-int main(int argc, char** argv) {
-  google::ParseCommandLineFlags(&argc, &argv, true);
-  google::InitGoogleLogging(argv[0]);
-
-  // The variable to solve for with its initial value.
-  double initial_x = 5.0;
-  double x = initial_x;
-
-  // Set up the only cost function (also known as residual). This uses
-  // numeric differentiation to obtain the derivative (jacobian).
-  CostFunction* cost =
-      new NumericDiffCostFunction<QuadraticCostFunctor, CENTRAL, 1, 1> (
-          new QuadraticCostFunctor);
-
-  // Build the problem.
-  Problem problem;
-  problem.AddResidualBlock(cost, NULL, &x);
-
-  // Run the solver!
-  Solver::Options options;
-  options.max_num_iterations = 10;
-  options.linear_solver_type = ceres::DENSE_QR;
-  options.minimizer_progress_to_stdout = true;
-  Solver::Summary summary;
-  Solve(options, &problem, &summary);
-  std::cout << summary.BriefReport() << "\n";
-  std::cout << "x : " << initial_x
-            << " -> " << x << "\n";
-  return 0;
-}
diff --git a/examples/random.h b/examples/random.h
new file mode 100644
index 0000000..0d55b49
--- /dev/null
+++ b/examples/random.h
@@ -0,0 +1,64 @@
+// Ceres Solver - A fast non-linear least squares minimizer
+// Copyright 2014 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_EXAMPLES_RANDOM_H_
+#define CERES_EXAMPLES_RANDOM_H_
+
+#include <math.h>
+#include <stdlib.h>
+
+namespace ceres {
+namespace examples {
+
+// Return a random number sampled from a uniform distribution in the range
+// [0,1].
+inline double RandDouble() {
+  double r = static_cast<double>(rand());
+  return r / RAND_MAX;
+}
+
+// Marsaglia Polar method for generation standard normal (pseudo)
+// random numbers http://en.wikipedia.org/wiki/Marsaglia_polar_method
+inline double RandNormal() {
+  double x1, x2, w;
+  do {
+    x1 = 2.0 * RandDouble() - 1.0;
+    x2 = 2.0 * RandDouble() - 1.0;
+    w = x1 * x1 + x2 * x2;
+  } while ( w >= 1.0 || w == 0.0 );
+
+  w = sqrt((-2.0 * log(w)) / w);
+  return x1 * w;
+}
+
+}  // namespace examples
+}  // namespace ceres
+
+#endif  // CERES_EXAMPLES_RANDOM_H_
diff --git a/examples/robot_pose_mle.cc b/examples/robot_pose_mle.cc
new file mode 100644
index 0000000..e1a1dd0
--- /dev/null
+++ b/examples/robot_pose_mle.cc
@@ -0,0 +1,316 @@
+// Ceres Solver - A fast non-linear least squares minimizer
+// Copyright 2014 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: joydeepb@ri.cmu.edu (Joydeep Biswas)
+//
+// This example demonstrates how to use the DynamicAutoDiffCostFunction
+// variant of CostFunction. The DynamicAutoDiffCostFunction is meant to
+// be used in cases where the number of parameter blocks or the sizes are not
+// known at compile time.
+//
+// This example simulates a robot traversing down a 1-dimension hallway with
+// noise odometry readings and noisy range readings of the end of the hallway.
+// By fusing the noisy odometry and sensor readings this example demonstrates
+// how to compute the maximum likelihood estimate (MLE) of the robot's pose at
+// each timestep.
+//
+// The robot starts at the origin, and it is travels to the end of a corridor of
+// fixed length specified by the "--corridor_length" flag. It executes a series
+// of motion commands to move forward a fixed length, specified by the
+// "--pose_separation" flag, at which pose it receives relative odometry
+// measurements as well as a range reading of the distance to the end of the
+// hallway. The odometry readings are drawn with Gaussian noise and standard
+// deviation specified by the "--odometry_stddev" flag, and the range readings
+// similarly with standard deviation specified by the "--range-stddev" flag.
+//
+// There are two types of residuals in this problem:
+// 1) The OdometryConstraint residual, that accounts for the odometry readings
+//    between successive pose estimatess of the robot.
+// 2) The RangeConstraint residual, that accounts for the errors in the observed
+//    range readings from each pose.
+//
+// The OdometryConstraint residual is modeled as an AutoDiffCostFunction with
+// a fixed parameter block size of 1, which is the relative odometry being
+// solved for, between a pair of successive poses of the robot. Differences
+// between observed and computed relative odometry values are penalized weighted
+// by the known standard deviation of the odometry readings.
+//
+// The RangeConstraint residual is modeled as a DynamicAutoDiffCostFunction
+// which sums up the relative odometry estimates to compute the estimated
+// global pose of the robot, and then computes the expected range reading.
+// Differences between the observed and expected range readings are then
+// penalized weighted by the standard deviation of readings of the sensor.
+// Since the number of poses of the robot is not known at compile time, this
+// cost function is implemented as a DynamicAutoDiffCostFunction.
+//
+// The outputs of the example are the initial values of the odometry and range
+// readings, and the range and odometry errors for every pose of the robot.
+// After computing the MLE, the computed poses and corrected odometry values
+// are printed out, along with the corresponding range and odometry errors. Note
+// that as an MLE of a noisy system the errors will not be reduced to zero, but
+// the odometry estimates will be updated to maximize the joint likelihood of
+// all odometry and range readings of the robot.
+//
+// Mathematical Formulation
+// ======================================================
+//
+// Let p_0, .., p_N be (N+1) robot poses, where the robot moves down the
+// corridor starting from p_0 and ending at p_N. We assume that p_0 is the
+// origin of the coordinate system.
+// Odometry u_i is the observed relative odometry between pose p_(i-1) and p_i,
+// and range reading y_i is the range reading of the end of the corridor from
+// pose p_i. Both odometry as well as range readings are noisy, but we wish to
+// compute the maximum likelihood estimate (MLE) of corrected odometry values
+// u*_0 to u*_(N-1), such that the Belief is optimized:
+//
+// Belief(u*_(0:N-1) | u_(0:N-1), y_(0:N-1))                                  1.
+//   =        P(u*_(0:N-1) | u_(0:N-1), y_(0:N-1))                            2.
+//   \propto  P(y_(0:N-1) | u*_(0:N-1), u_(0:N-1)) P(u*_(0:N-1) | u_(0:N-1))  3.
+//   =       \prod_i{ P(y_i | u*_(0:i)) P(u*_i | u_i) }                       4.
+//
+// Here, the subscript "(0:i)" is used as shorthand to indicate entries from all
+// timesteps 0 to i for that variable, both inclusive.
+//
+// Bayes' rule is used to derive eq. 3 from 2, and the independence of
+// odometry observations and range readings is expolited to derive 4 from 3.
+//
+// Thus, the Belief, up to scale, is factored as a product of a number of
+// terms, two for each pose, where for each pose term there is one term for the
+// range reading, P(y_i | u*_(0:i) and one term for the odometry reading,
+// P(u*_i | u_i) . Note that the term for the range reading is dependent on all
+// odometry values u*_(0:i), while the odometry term, P(u*_i | u_i) depends only
+// on a single value, u_i. Both the range reading as well as odoemtry
+// probability terms are modeled as the Normal distribution, and have the form:
+//
+// p(x) \propto \exp{-((x - x_mean) / x_stddev)^2}
+//
+// where x refers to either the MLE odometry u* or range reading y, and x_mean
+// is the corresponding mean value, u for the odometry terms, and y_expected,
+// the expected range reading based on all the previous odometry terms.
+// The MLE is thus found by finding those values x* which minimize:
+//
+// x* = \arg\min{((x - x_mean) / x_stddev)^2}
+//
+// which is in the nonlinear least-square form, suited to being solved by Ceres.
+// The non-linear component arise from the computation of x_mean. The residuals
+// ((x - x_mean) / x_stddev) for the residuals that Ceres will optimize. As
+// mentioned earlier, the odometry term for each pose depends only on one
+// variable, and will be computed by an AutoDiffCostFunction, while the term
+// for the range reading will depend on all previous odometry observations, and
+// will be computed by a DynamicAutoDiffCostFunction since the number of
+// odoemtry observations will only be known at run time.
+
+#include <cstdio>
+#include <math.h>
+#include <vector>
+
+#include "ceres/ceres.h"
+#include "ceres/dynamic_autodiff_cost_function.h"
+#include "gflags/gflags.h"
+#include "glog/logging.h"
+#include "random.h"
+
+using ceres::AutoDiffCostFunction;
+using ceres::DynamicAutoDiffCostFunction;
+using ceres::CauchyLoss;
+using ceres::CostFunction;
+using ceres::LossFunction;
+using ceres::Problem;
+using ceres::Solve;
+using ceres::Solver;
+using ceres::examples::RandNormal;
+using std::min;
+using std::vector;
+
+DEFINE_double(corridor_length, 30.0, "Length of the corridor that the robot is "
+              "travelling down.");
+
+DEFINE_double(pose_separation, 0.5, "The distance that the robot traverses "
+              "between successive odometry updates.");
+
+DEFINE_double(odometry_stddev, 0.1, "The standard deviation of "
+              "odometry error of the robot.");
+
+DEFINE_double(range_stddev, 0.01, "The standard deviation of range readings of "
+              "the robot.");
+
+// The stride length of the dynamic_autodiff_cost_function evaluator.
+static const int kStride = 10;
+
+struct OdometryConstraint {
+  typedef AutoDiffCostFunction<OdometryConstraint, 1, 1> OdometryCostFunction;
+
+  OdometryConstraint(double odometry_mean, double odometry_stddev) :
+      odometry_mean(odometry_mean), odometry_stddev(odometry_stddev) {}
+
+  template <typename T>
+  bool operator()(const T* const odometry, T* residual) const {
+    *residual = (*odometry - T(odometry_mean)) / T(odometry_stddev);
+    return true;
+  }
+
+  static OdometryCostFunction* Create(const double odometry_value) {
+    return new OdometryCostFunction(
+        new OdometryConstraint(odometry_value, FLAGS_odometry_stddev));
+  }
+
+  const double odometry_mean;
+  const double odometry_stddev;
+};
+
+struct RangeConstraint {
+  typedef DynamicAutoDiffCostFunction<RangeConstraint, kStride>
+      RangeCostFunction;
+
+  RangeConstraint(
+      int pose_index,
+      double range_reading,
+      double range_stddev,
+      double corridor_length) :
+      pose_index(pose_index), range_reading(range_reading),
+      range_stddev(range_stddev), corridor_length(corridor_length) {}
+
+  template <typename T>
+  bool operator()(T const* const* relative_poses, T* residuals) const {
+    T global_pose(0);
+    for (int i = 0; i <= pose_index; ++i) {
+      global_pose += relative_poses[i][0];
+    }
+    residuals[0] = (global_pose + T(range_reading) - T(corridor_length)) /
+        T(range_stddev);
+    return true;
+  }
+
+  // Factory method to create a CostFunction from a RangeConstraint to
+  // conveniently add to a ceres problem.
+  static RangeCostFunction* Create(const int pose_index,
+                                   const double range_reading,
+                                   vector<double>* odometry_values,
+                                   vector<double*>* parameter_blocks) {
+    RangeConstraint* constraint = new RangeConstraint(
+        pose_index, range_reading, FLAGS_range_stddev, FLAGS_corridor_length);
+    RangeCostFunction* cost_function = new RangeCostFunction(constraint);
+    // Add all the parameter blocks that affect this constraint.
+    parameter_blocks->clear();
+    for (int i = 0; i <= pose_index; ++i) {
+      parameter_blocks->push_back(&((*odometry_values)[i]));
+      cost_function->AddParameterBlock(1);
+    }
+    cost_function->SetNumResiduals(1);
+    return (cost_function);
+  }
+
+  const int pose_index;
+  const double range_reading;
+  const double range_stddev;
+  const double corridor_length;
+};
+
+void SimulateRobot(vector<double>* odometry_values,
+                   vector<double>* range_readings) {
+  const int num_steps = static_cast<int>(
+      ceil(FLAGS_corridor_length / FLAGS_pose_separation));
+
+  // The robot starts out at the origin.
+  double robot_location = 0.0;
+  for (int i = 0; i < num_steps; ++i) {
+    const double actual_odometry_value = min(
+        FLAGS_pose_separation, FLAGS_corridor_length - robot_location);
+    robot_location += actual_odometry_value;
+    const double actual_range = FLAGS_corridor_length - robot_location;
+    const double observed_odometry =
+        RandNormal() * FLAGS_odometry_stddev + actual_odometry_value;
+    const double observed_range =
+        RandNormal() * FLAGS_range_stddev + actual_range;
+    odometry_values->push_back(observed_odometry);
+    range_readings->push_back(observed_range);
+  }
+}
+
+void PrintState(const vector<double>& odometry_readings,
+                const vector<double>& range_readings) {
+  CHECK_EQ(odometry_readings.size(), range_readings.size());
+  double robot_location = 0.0;
+  printf("pose: location     odom    range  r.error  o.error\n");
+  for (int i = 0; i < odometry_readings.size(); ++i) {
+    robot_location += odometry_readings[i];
+    const double range_error =
+        robot_location + range_readings[i] - FLAGS_corridor_length;
+    const double odometry_error =
+        FLAGS_pose_separation - odometry_readings[i];
+    printf("%4d: %8.3f %8.3f %8.3f %8.3f %8.3f\n",
+           static_cast<int>(i), robot_location, odometry_readings[i],
+           range_readings[i], range_error, odometry_error);
+  }
+}
+
+int main(int argc, char** argv) {
+  google::InitGoogleLogging(argv[0]);
+  google::ParseCommandLineFlags(&argc, &argv, true);
+  // Make sure that the arguments parsed are all positive.
+  CHECK_GT(FLAGS_corridor_length, 0.0);
+  CHECK_GT(FLAGS_pose_separation, 0.0);
+  CHECK_GT(FLAGS_odometry_stddev, 0.0);
+  CHECK_GT(FLAGS_range_stddev, 0.0);
+
+  vector<double> odometry_values;
+  vector<double> range_readings;
+  SimulateRobot(&odometry_values, &range_readings);
+
+  printf("Initial values:\n");
+  PrintState(odometry_values, range_readings);
+  ceres::Problem problem;
+
+  for (int i = 0; i < odometry_values.size(); ++i) {
+    // Create and add a DynamicAutoDiffCostFunction for the RangeConstraint from
+    // pose i.
+    vector<double*> parameter_blocks;
+    RangeConstraint::RangeCostFunction* range_cost_function =
+        RangeConstraint::Create(
+            i, range_readings[i], &odometry_values, &parameter_blocks);
+    problem.AddResidualBlock(range_cost_function, NULL, parameter_blocks);
+
+    // Create and add an AutoDiffCostFunction for the OdometryConstraint for
+    // pose i.
+    problem.AddResidualBlock(OdometryConstraint::Create(odometry_values[i]),
+                             NULL,
+                             &(odometry_values[i]));
+  }
+
+  ceres::Solver::Options solver_options;
+  solver_options.minimizer_progress_to_stdout = true;
+
+  Solver::Summary summary;
+  printf("Solving...\n");
+  Solve(solver_options, &problem, &summary);
+  printf("Done.\n");
+  std::cout << summary.FullReport() << "\n";
+  printf("Final values:\n");
+  PrintState(odometry_values, range_readings);
+  return 0;
+}
diff --git a/examples/snavely_reprojection_error.h b/examples/snavely_reprojection_error.h
index 0704217..d3263f3 100644
--- a/examples/snavely_reprojection_error.h
+++ b/examples/snavely_reprojection_error.h
@@ -91,6 +91,14 @@ struct SnavelyReprojectionError {
     return true;
   }
 
+  // Factory to hide the construction of the CostFunction object from
+  // the client code.
+  static ceres::CostFunction* Create(const double observed_x,
+                                     const double observed_y) {
+    return (new ceres::AutoDiffCostFunction<SnavelyReprojectionError, 2, 9, 3>(
+                new SnavelyReprojectionError(observed_x, observed_y)));
+  }
+
   double observed_x;
   double observed_y;
 };
@@ -146,6 +154,16 @@ struct SnavelyReprojectionErrorWithQuaternions {
     return true;
   }
 
+  // Factory to hide the construction of the CostFunction object from
+  // the client code.
+  static ceres::CostFunction* Create(const double observed_x,
+                                     const double observed_y) {
+    return (new ceres::AutoDiffCostFunction<
+            SnavelyReprojectionErrorWithQuaternions, 2, 4, 6, 3>(
+                new SnavelyReprojectionErrorWithQuaternions(observed_x,
+                                                            observed_y)));
+  }
+
   double observed_x;
   double observed_y;
 };