Update ceres to the latest version in google3.

Change-Id: I0165fffa55f60714f23e0096eac89fa68df75a05

16 files changed, 1879 insertions, 113 deletions

diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt
index 2307a03..9dfc80b 100644
--- a/examples/CMakeLists.txt
+++ b/examples/CMakeLists.txt
@@ -28,39 +28,51 @@
 #
 # Author: keir@google.com (Keir Mierle)
 
-IF (${GFLAGS})
-  ADD_EXECUTABLE(quadratic quadratic.cc)
-  TARGET_LINK_LIBRARIES(quadratic ceres)
+ADD_EXECUTABLE(helloworld helloworld.cc)
+TARGET_LINK_LIBRARIES(helloworld ceres)
 
-  ADD_EXECUTABLE(nist nist.cc)
-  TARGET_LINK_LIBRARIES(nist ceres)
+ADD_EXECUTABLE(helloworld_numeric_diff helloworld_numeric_diff.cc)
+TARGET_LINK_LIBRARIES(helloworld_numeric_diff ceres)
+
+ADD_EXECUTABLE(helloworld_analytic_diff helloworld_analytic_diff.cc)
+TARGET_LINK_LIBRARIES(helloworld_analytic_diff ceres)
+
+ADD_EXECUTABLE(curve_fitting curve_fitting.cc)
+TARGET_LINK_LIBRARIES(curve_fitting ceres)
+
+ADD_EXECUTABLE(curve_fitting_c curve_fitting.c)
+TARGET_LINK_LIBRARIES(curve_fitting_c ceres)
 
-  ADD_EXECUTABLE(quadratic_auto_diff quadratic_auto_diff.cc)
-  TARGET_LINK_LIBRARIES(quadratic_auto_diff ceres)
 
-  ADD_EXECUTABLE(quadratic_numeric_diff quadratic_numeric_diff.cc)
-  TARGET_LINK_LIBRARIES(quadratic_numeric_diff ceres)
+ADD_EXECUTABLE(robust_curve_fitting robust_curve_fitting.cc)
+TARGET_LINK_LIBRARIES(robust_curve_fitting ceres)
+
+ADD_EXECUTABLE(simple_bundle_adjuster
+               simple_bundle_adjuster.cc)
+TARGET_LINK_LIBRARIES(simple_bundle_adjuster ceres)
 
+IF (GFLAGS)
   ADD_EXECUTABLE(powell powell.cc)
   TARGET_LINK_LIBRARIES(powell ceres)
 
+  ADD_EXECUTABLE(nist nist.cc)
+  TARGET_LINK_LIBRARIES(nist ceres)
+
   ADD_EXECUTABLE(circle_fit circle_fit.cc)
   TARGET_LINK_LIBRARIES(circle_fit ceres)
 
-  ADD_EXECUTABLE(data_fitting data_fitting.cc)
-  TARGET_LINK_LIBRARIES(data_fitting ceres)
-
   ADD_EXECUTABLE(bundle_adjuster
                  bundle_adjuster.cc
                  bal_problem.cc)
   TARGET_LINK_LIBRARIES(bundle_adjuster ceres)
 
+  ADD_EXECUTABLE(libmv_bundle_adjuster
+                 libmv_bundle_adjuster.cc)
+  TARGET_LINK_LIBRARIES(libmv_bundle_adjuster ceres)
+
   ADD_EXECUTABLE(denoising
                  denoising.cc
                  fields_of_experts.cc)
   TARGET_LINK_LIBRARIES(denoising ceres)
-ENDIF (${GFLAGS})
+ENDIF (GFLAGS)
 
-ADD_EXECUTABLE(simple_bundle_adjuster
-               simple_bundle_adjuster.cc)
-TARGET_LINK_LIBRARIES(simple_bundle_adjuster ceres)
diff --git a/examples/bal_problem.cc b/examples/bal_problem.cc
index 5733f46..c118f88 100644
--- a/examples/bal_problem.cc
+++ b/examples/bal_problem.cc
@@ -35,7 +35,6 @@
 #include <string>
 #include <vector>
 #include "Eigen/Core"
-#include "ceres/random.h"
 #include "ceres/rotation.h"
 #include "glog/logging.h"
 
@@ -45,6 +44,25 @@ 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 78dbd01..c060aed 100644
--- a/examples/bundle_adjuster.cc
+++ b/examples/bundle_adjuster.cc
@@ -60,7 +60,6 @@
 
 #include "bal_problem.h"
 #include "ceres/ceres.h"
-#include "ceres/random.h"
 #include "gflags/gflags.h"
 #include "glog/logging.h"
 #include "snavely_reprojection_error.h"
@@ -97,9 +96,6 @@ DEFINE_double(eta, 1e-2, "Default value for eta. Eta determines the "
              "accuracy of each linear solve of the truncated newton step. "
              "Changing this parameter can affect solve performance.");
 
-DEFINE_bool(use_block_amd, true, "Use a block oriented fill reducing "
-            "ordering.");
-
 DEFINE_int32(num_threads, 1, "Number of threads.");
 DEFINE_int32(num_iterations, 5, "Number of iterations.");
 DEFINE_double(max_solver_time, 1e32, "Maximum solve time in seconds.");
@@ -116,6 +112,8 @@ 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.");
 
 namespace ceres {
 namespace examples {
@@ -140,8 +138,6 @@ void SetOrdering(BALProblem* bal_problem, Solver::Options* options) {
   const int camera_block_size = bal_problem->camera_block_size();
   double* cameras = bal_problem->mutable_cameras();
 
-  options->use_block_amd = FLAGS_use_block_amd;
-
   if (options->use_inner_iterations) {
     if (FLAGS_blocks_for_inner_iterations == "cameras") {
       LOG(INFO) << "Camera blocks for inner iterations";
@@ -225,6 +221,10 @@ void SetMinimizerOptions(Solver::Options* options) {
   options->eta = FLAGS_eta;
   options->max_solver_time_in_seconds = FLAGS_max_solver_time;
   options->use_nonmonotonic_steps = FLAGS_nonmonotonic_steps;
+  if (FLAGS_line_search) {
+    options->minimizer_type = ceres::LINE_SEARCH;
+  }
+
   CHECK(StringToTrustRegionStrategyType(FLAGS_trust_region_strategy,
                                         &options->trust_region_strategy_type));
   CHECK(StringToDoglegType(FLAGS_dogleg, &options->dogleg_type));
@@ -305,7 +305,7 @@ void SolveProblem(const char* filename) {
   BALProblem bal_problem(filename, FLAGS_use_quaternions);
   Problem problem;
 
-  SetRandomState(FLAGS_random_seed);
+  srand(FLAGS_random_seed);
   bal_problem.Normalize();
   bal_problem.Perturb(FLAGS_rotation_sigma,
                       FLAGS_translation_sigma,
diff --git a/examples/curve_fitting.c b/examples/curve_fitting.c
new file mode 100644
index 0000000..3ea2ef1
--- /dev/null
+++ b/examples/curve_fitting.c
@@ -0,0 +1,187 @@
+/* Ceres Solver - A fast non-linear least squares minimizer
+ * Copyright 2013 Google Inc. All rights reserved.
+ * http://code.google.com/p/ceres-solver/
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions are met:
+ *
+ * - Redistributions of source code must retain the above copyright notice,
+ *   this list of conditions and the following disclaimer.
+ * - Redistributions in binary form must reproduce the above copyright notice,
+ *   this list of conditions and the following disclaimer in the documentation
+ *   and/or other materials provided with the distribution.
+ * - Neither the name of Google Inc. nor the names of its contributors may be
+ *   used to endorse or promote products derived from this software without
+ *   specific prior written permission.
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+ * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+ * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+ * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
+ * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
+ * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
+ * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
+ * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
+ * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
+ * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
+ * POSSIBILITY OF SUCH DAMAGE.
+ *
+ * Author: mierle@gmail.com (Keir Mierle)
+ *
+ * This is a port of curve_fitting.cc to the minimal C API for Ceres.
+ */
+
+#include <math.h>
+#include <stdio.h>
+#include <string.h>  // For NULL
+#include "ceres/c_api.h"
+
+/* 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'];
+ *
+ */
+
+int num_observations = 67;
+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,
+};
+
+/* This is the equivalent of a use-defined CostFunction in the C++ Ceres API.
+ * This is passed as a callback to the Ceres C API, which internally converts
+ * the callback into a CostFunction. */
+int exponential_residual(void* user_data,
+                         double** parameters,
+                         double* residuals,
+                         double** jacobians) {
+  double* measurement = (double*) user_data;
+  double x = measurement[0];
+  double y = measurement[1];
+  double m = parameters[0][0];
+  double c = parameters[1][0];
+
+  residuals[0] = y - exp(m * x + c);
+  if (jacobians == NULL) {
+    return 1;
+  }
+  if (jacobians[0] != NULL) {
+    jacobians[0][0] = - x * exp(m * x + c);  /* dr/dm */
+  }
+  if (jacobians[1] != NULL) {
+    jacobians[1][0] =     - exp(m * x + c);  /* dr/dc */
+  }
+  return 1;
+}
+
+int main(int argc, char** argv) {
+  /* Note: Typically it is better to compact m and c into one block,
+   * but in this case use separate blocks to illustrate the use of
+   * multiple parameter blocks. */
+  double m = 0.0;
+  double c = 0.0;
+
+  double *parameter_pointers[] = { &m, &c };
+  int parameter_sizes[] = { 1, 1 };
+  int i;
+
+  ceres_problem_t* problem;
+
+  /* Ceres has some internal stuff that needs to get initialized. */
+  ceres_init();
+
+  problem = ceres_create_problem();
+
+  /* Add all the residuals. */
+  for (i = 0; i < num_observations; ++i) {
+    ceres_problem_add_residual_block(
+        problem,
+        exponential_residual,  /* Cost function */
+        &data[2 * i],          /* Points to the (x,y) measurement */
+        NULL,                  /* No loss function */
+        NULL,                  /* No loss function user data */
+        1,                     /* Number of residuals */
+        2,                     /* Number of parameter blocks */
+        parameter_sizes,
+        parameter_pointers);
+  }
+
+  ceres_solve(problem);
+  ceres_free_problem(problem);
+
+  printf("Initial m: 0.0, c: 0.0\n");
+  printf("Final m: %g, c: %g\n", m, c);
+  return 0;
+}
diff --git a/examples/curve_fitting.cc b/examples/curve_fitting.cc
new file mode 100644
index 0000000..3e3adcd
--- /dev/null
+++ b/examples/curve_fitting.cc
@@ -0,0 +1,163 @@
+// 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 "glog/logging.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,
+};
+
+struct ExponentialResidual {
+  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::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/denoising.cc b/examples/denoising.cc
index 086be00..bfd8118 100644
--- a/examples/denoising.cc
+++ b/examples/denoising.cc
@@ -57,6 +57,8 @@ DEFINE_string(foe_file, "", "FoE file to use");
 DEFINE_string(output, "", "File to which the output image should be written");
 DEFINE_double(sigma, 20.0, "Standard deviation of noise");
 DEFINE_bool(verbose, false, "Prints information about the solver progress.");
+DEFINE_bool(line_search, false, "Use a line search instead of trust region "
+            "algorithm.");
 
 namespace ceres {
 namespace examples {
@@ -143,7 +145,11 @@ void SolveProblem(Problem* problem, PGMImage<double>* solution) {
   if (FLAGS_verbose) {
     options.minimizer_progress_to_stdout = true;
   }
-  options.trust_region_strategy_type = ceres::LEVENBERG_MARQUARDT;
+
+  if (FLAGS_line_search) {
+    options.minimizer_type = ceres::LINE_SEARCH;
+  }
+
   options.linear_solver_type = ceres::SPARSE_NORMAL_CHOLESKY;
   options.function_tolerance = 1e-3;  // Enough for denoising.
 
diff --git a/examples/helloworld.cc b/examples/helloworld.cc
new file mode 100644
index 0000000..6341918
--- /dev/null
+++ b/examples/helloworld.cc
@@ -0,0 +1,83 @@
+// 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 "ceres/ceres.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.
+struct CostFunctor {
+  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::InitGoogleLogging(argv[0]);
+
+  // The variable to solve for with its initial value. It will be
+  // mutated in place by the solver.
+  double x = 0.5;
+  const double initial_x = 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).
+  CostFunction* cost_function =
+      new AutoDiffCostFunction<CostFunctor, 1, 1>(new CostFunctor);
+  problem.AddResidualBlock(cost_function, NULL, &x);
+
+  // Run the solver!
+  Solver::Options options;
+  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/helloworld_analytic_diff.cc b/examples/helloworld_analytic_diff.cc
new file mode 100644
index 0000000..bff4804
--- /dev/null
+++ b/examples/helloworld_analytic_diff.cc
@@ -0,0 +1,107 @@
+// 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 "glog/logging.h"
+
+using ceres::CostFunction;
+using ceres::SizedCostFunction;
+using ceres::Problem;
+using ceres::Solver;
+using ceres::Solve;
+
+// A CostFunction implementing analytically derivatives for the
+// function f(x) = 10 - x.
+class QuadraticCostFunction
+  : public SizedCostFunction<1 /* number of residuals */,
+                             1 /* size of first parameter */> {
+ public:
+  virtual ~QuadraticCostFunction() {}
+
+  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.
+    //
+    // Since the Evaluate function can be called with the jacobians
+    // pointer equal to NULL, the Evaluate function must check to see
+    // if jacobians need to be computed.
+    //
+    // For this simple problem it is overkill to check if jacobians[0]
+    // is NULL, but in general when writing more complex
+    // CostFunctions, it is possible that Ceres may only demand the
+    // derivatives w.r.t. a subset of the parameter blocks.
+    if (jacobians != NULL && jacobians[0] != NULL) {
+      jacobians[0][0] = -1;
+    }
+
+    return true;
+  }
+};
+
+int main(int argc, char** argv) {
+  google::InitGoogleLogging(argv[0]);
+
+  // The variable to solve for with its initial value. It will be
+  // mutated in place by the solver.
+  double x = 0.5;
+  const double initial_x = x;
+
+  // Build the problem.
+  Problem problem;
+
+  // Set up the only cost function (also known as residual).
+  CostFunction* cost_function = new QuadraticCostFunction;
+  problem.AddResidualBlock(cost_function, NULL, &x);
+
+  // Run the solver!
+  Solver::Options options;
+  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/helloworld_numeric_diff.cc b/examples/helloworld_numeric_diff.cc
new file mode 100644
index 0000000..026b155
--- /dev/null
+++ b/examples/helloworld_numeric_diff.cc
@@ -0,0 +1,79 @@
+// 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 "ceres/ceres.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.
+struct CostFunctor {
+  bool operator()(const double* const x, double* residual) const {
+    residual[0] = 10.0 - x[0];
+    return true;
+  }
+};
+
+int main(int argc, char** argv) {
+  google::InitGoogleLogging(argv[0]);
+
+  // The variable to solve for with its initial value. It will be
+  // mutated in place by the solver.
+  double x = 0.5;
+  const double initial_x = x;
+
+  // Build the problem.
+  Problem problem;
+
+  // Set up the only cost function (also known as residual). This uses
+  // numeric differentiation to obtain the derivative (jacobian).
+  CostFunction* cost_function =
+      new NumericDiffCostFunction<CostFunctor, CENTRAL, 1, 1> (new CostFunctor);
+  problem.AddResidualBlock(cost_function, NULL, &x);
+
+  // Run the solver!
+  Solver::Options options;
+  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/libmv_bundle_adjuster.cc b/examples/libmv_bundle_adjuster.cc
new file mode 100644
index 0000000..4b01202
--- /dev/null
+++ b/examples/libmv_bundle_adjuster.cc
@@ -0,0 +1,842 @@
+// Copyright (c) 2013 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: mierle@gmail.com (Keir Mierle)
+//         sergey.vfx@gmail.com (Sergey Sharybin)
+//
+// This is an example application which contains bundle adjustment code used
+// in the Libmv library and Blender. It reads problems from files passed via
+// the command line and runs the bundle adjuster on the problem.
+//
+// File with problem a binary file, for which it is crucial to know in which
+// order bytes of float values are stored in. This information is provided
+// by a single character in the beginning of the file. There're two possible
+// values of this byte:
+//  - V, which means values in the file are stored with big endian type
+//  - v, which means values in the file are stored with little endian type
+//
+// The rest of the file contains data in the following order:
+//   - Space in which markers' coordinates are stored in
+//   - Camera intrinsics
+//   - Number of cameras
+//   - Cameras
+//   - Number of 3D points
+//   - 3D points
+//   - Number of markers
+//   - Markers
+//
+// Markers' space could either be normalized or image (pixels). This is defined
+// by the single character in the file. P means markers in the file is in image
+// space, and N means markers are in normalized space.
+//
+// Camera intrinsics are 8 described by 8 float 8.
+// This values goes in the following order:
+//
+//   - Focal length, principal point X, principal point Y, k1, k2, k3, p1, p2
+//
+// Every camera is described by:
+//
+//   - Image for which camera belongs to (single 4 bytes integer value).
+//   - Column-major camera rotation matrix, 9 float values.
+//   - Camera translation, 3-component vector of float values.
+//
+// Image number shall be greater or equal to zero. Order of cameras does not
+// matter and gaps are possible.
+//
+// Every 3D point is decribed by:
+//
+//  - Track number point belongs to (single 4 bytes integer value).
+//  - 3D position vector, 3-component vector of float values.
+//
+// Track number shall be greater or equal to zero. Order of tracks does not
+// matter and gaps are possible.
+//
+// Finally every marker is described by:
+//
+//  - Image marker belongs to single 4 bytes integer value).
+//  - Track marker belongs to single 4 bytes integer value).
+//  - 2D marker position vector, (two float values).
+//
+// Marker's space is used a default value for refine_intrinsics command line
+// flag. This means if there's no refine_intrinsics flag passed via command
+// line, camera intrinsics will be refined if markers in the problem are
+// stored in image space and camera intrinsics will not be refined if markers
+// are in normalized space.
+//
+// Passing refine_intrinsics command line flag defines explicitly whether
+// refinement of intrinsics will happen. Currently, only none and all
+// intrinsics refinement is supported.
+//
+// There're existing problem files dumped from blender stored in folder
+// ../data/libmv-ba-problems.
+
+#include <cstdio>
+#include <fcntl.h>
+#include <sstream>
+#include <string>
+#include <vector>
+
+#ifdef _MSC_VER
+#  include <io.h>
+#  define open _open
+#  define close _close
+typedef unsigned __int32 uint32_t;
+#else
+# include <stdint.h>
+
+// O_BINARY is not defined on unix like platforms, as there is no
+// difference between binary and text files.
+#define O_BINARY 0
+
+#endif
+
+#include "ceres/ceres.h"
+#include "ceres/rotation.h"
+#include "gflags/gflags.h"
+#include "glog/logging.h"
+
+typedef Eigen::Matrix<double, 3, 3> Mat3;
+typedef Eigen::Matrix<double, 6, 1> Vec6;
+typedef Eigen::Vector3d Vec3;
+typedef Eigen::Vector4d Vec4;
+
+using std::vector;
+
+DEFINE_string(input, "", "Input File name");
+DEFINE_string(refine_intrinsics, "", "Camera intrinsics to be refined. "
+              "Options are: none, radial.");
+
+namespace {
+
+// A EuclideanCamera is the location and rotation of the camera
+// viewing an image.
+//
+// image identifies which image this camera represents.
+// R is a 3x3 matrix representing the rotation of the camera.
+// t is a translation vector representing its positions.
+struct EuclideanCamera {
+  EuclideanCamera() : image(-1) {}
+  EuclideanCamera(const EuclideanCamera &c) : image(c.image), R(c.R), t(c.t) {}
+
+  int image;
+  Mat3 R;
+  Vec3 t;
+};
+
+// A Point is the 3D location of a track.
+//
+// track identifies which track this point corresponds to.
+// X represents the 3D position of the track.
+struct EuclideanPoint {
+  EuclideanPoint() : track(-1) {}
+  EuclideanPoint(const EuclideanPoint &p) : track(p.track), X(p.X) {}
+  int track;
+  Vec3 X;
+};
+
+// A Marker is the 2D location of a tracked point in an image.
+//
+// x and y is the position of the marker in pixels from the top left corner
+// in the image identified by an image. All markers for to the same target
+// form a track identified by a common track number.
+struct Marker {
+  int image;
+  int track;
+  double x, y;
+};
+
+// Cameras intrinsics to be bundled.
+//
+// BUNDLE_RADIAL actually implies bundling of k1 and k2 coefficients only,
+// no bundling of k3 is possible at this moment.
+enum BundleIntrinsics {
+  BUNDLE_NO_INTRINSICS = 0,
+  BUNDLE_FOCAL_LENGTH = 1,
+  BUNDLE_PRINCIPAL_POINT = 2,
+  BUNDLE_RADIAL_K1 = 4,
+  BUNDLE_RADIAL_K2 = 8,
+  BUNDLE_RADIAL = 12,
+  BUNDLE_TANGENTIAL_P1 = 16,
+  BUNDLE_TANGENTIAL_P2 = 32,
+  BUNDLE_TANGENTIAL = 48,
+};
+
+// Denotes which blocks to keep constant during bundling.
+// For example it is useful to keep camera translations constant
+// when bundling tripod motions.
+enum BundleConstraints {
+  BUNDLE_NO_CONSTRAINTS = 0,
+  BUNDLE_NO_TRANSLATION = 1,
+};
+
+// The intrinsics need to get combined into a single parameter block; use these
+// enums to index instead of numeric constants.
+enum {
+  OFFSET_FOCAL_LENGTH,
+  OFFSET_PRINCIPAL_POINT_X,
+  OFFSET_PRINCIPAL_POINT_Y,
+  OFFSET_K1,
+  OFFSET_K2,
+  OFFSET_K3,
+  OFFSET_P1,
+  OFFSET_P2,
+};
+
+// Returns a pointer to the camera corresponding to a image.
+EuclideanCamera *CameraForImage(vector<EuclideanCamera> *all_cameras,
+                                const int image) {
+  if (image < 0 || image >= all_cameras->size()) {
+    return NULL;
+  }
+  EuclideanCamera *camera = &(*all_cameras)[image];
+  if (camera->image == -1) {
+    return NULL;
+  }
+  return camera;
+}
+
+const EuclideanCamera *CameraForImage(
+    const vector<EuclideanCamera> &all_cameras,
+    const int image) {
+  if (image < 0 || image >= all_cameras.size()) {
+    return NULL;
+  }
+  const EuclideanCamera *camera = &all_cameras[image];
+  if (camera->image == -1) {
+    return NULL;
+  }
+  return camera;
+}
+
+// Returns maximal image number at which marker exists.
+int MaxImage(const vector<Marker> &all_markers) {
+  if (all_markers.size() == 0) {
+    return -1;
+  }
+
+  int max_image = all_markers[0].image;
+  for (int i = 1; i < all_markers.size(); i++) {
+    max_image = std::max(max_image, all_markers[i].image);
+  }
+  return max_image;
+}
+
+// Returns a pointer to the point corresponding to a track.
+EuclideanPoint *PointForTrack(vector<EuclideanPoint> *all_points,
+                              const int track) {
+  if (track < 0 || track >= all_points->size()) {
+    return NULL;
+  }
+  EuclideanPoint *point = &(*all_points)[track];
+  if (point->track == -1) {
+    return NULL;
+  }
+  return point;
+}
+
+// Reader of binary file which makes sure possibly needed endian
+// conversion happens when loading values like floats and integers.
+//
+// File's endian type is reading from a first character of file, which
+// could either be V for big endian or v for little endian.  This
+// means you need to design file format assuming first character
+// denotes file endianness in this way.
+class EndianAwareFileReader {
+ public:
+  EndianAwareFileReader(void) : file_descriptor_(-1) {
+    // Get an endian type of the host machine.
+    union {
+      unsigned char bytes[4];
+      uint32_t value;
+    } endian_test = { { 0, 1, 2, 3 } };
+    host_endian_type_ = endian_test.value;
+    file_endian_type_ = host_endian_type_;
+  }
+
+  ~EndianAwareFileReader(void) {
+    if (file_descriptor_ > 0) {
+      close(file_descriptor_);
+    }
+  }
+
+  bool OpenFile(const std::string &file_name) {
+    file_descriptor_ = open(file_name.c_str(), O_RDONLY | O_BINARY);
+    if (file_descriptor_ < 0) {
+      return false;
+    }
+    // Get an endian tpye of data in the file.
+    unsigned char file_endian_type_flag = Read<unsigned char>();
+    if (file_endian_type_flag == 'V') {
+      file_endian_type_ = kBigEndian;
+    } else if (file_endian_type_flag == 'v') {
+      file_endian_type_ = kLittleEndian;
+    } else {
+      LOG(FATAL) << "Problem file is stored in unknown endian type.";
+    }
+    return true;
+  }
+
+  // Read value from the file, will switch endian if needed.
+  template <typename T>
+  T Read(void) const {
+    T value;
+    CHECK_GT(read(file_descriptor_, &value, sizeof(value)), 0);
+    // Switch endian type if file contains data in different type
+    // that current machine.
+    if (file_endian_type_ != host_endian_type_) {
+      value = SwitchEndian<T>(value);
+    }
+    return value;
+  }
+ private:
+  static const long int kLittleEndian = 0x03020100ul;
+  static const long int kBigEndian = 0x00010203ul;
+
+  // Switch endian type between big to little.
+  template <typename T>
+  T SwitchEndian(const T value) const {
+    if (sizeof(T) == 4) {
+      unsigned int temp_value = static_cast<unsigned int>(value);
+      return ((temp_value >> 24)) |
+             ((temp_value << 8) & 0x00ff0000) |
+             ((temp_value >> 8) & 0x0000ff00) |
+             ((temp_value << 24));
+    } else if (sizeof(T) == 1) {
+      return value;
+    } else {
+      LOG(FATAL) << "Entered non-implemented part of endian switching function.";
+    }
+  }
+
+  int host_endian_type_;
+  int file_endian_type_;
+  int file_descriptor_;
+};
+
+// Read 3x3 column-major matrix from the file
+void ReadMatrix3x3(const EndianAwareFileReader &file_reader,
+                   Mat3 *matrix) {
+  for (int i = 0; i < 9; i++) {
+    (*matrix)(i % 3, i / 3) = file_reader.Read<float>();
+  }
+}
+
+// Read 3-vector from file
+void ReadVector3(const EndianAwareFileReader &file_reader,
+                 Vec3 *vector) {
+  for (int i = 0; i < 3; i++) {
+    (*vector)(i) = file_reader.Read<float>();
+  }
+}
+
+// Reads a bundle adjustment problem from the file.
+//
+// file_name denotes from which file to read the problem.
+// camera_intrinsics will contain initial camera intrinsics values.
+//
+// all_cameras is a vector of all reconstructed cameras to be optimized,
+// vector element with number i will contain camera for image i.
+//
+// all_points is a vector of all reconstructed 3D points to be optimized,
+// vector element with number i will contain point for track i.
+//
+// all_markers is a vector of all tracked markers existing in
+// the problem. Only used for reprojection error calculation, stay
+// unchanged during optimization.
+//
+// Returns false if any kind of error happened during
+// reading.
+bool ReadProblemFromFile(const std::string &file_name,
+                         double camera_intrinsics[8],
+                         vector<EuclideanCamera> *all_cameras,
+                         vector<EuclideanPoint> *all_points,
+                         bool *is_image_space,
+                         vector<Marker> *all_markers) {
+  EndianAwareFileReader file_reader;
+  if (!file_reader.OpenFile(file_name)) {
+    return false;
+  }
+
+  // Read markers' space flag.
+  unsigned char is_image_space_flag = file_reader.Read<unsigned char>();
+  if (is_image_space_flag == 'P') {
+    *is_image_space = true;
+  } else if (is_image_space_flag == 'N') {
+    *is_image_space = false;
+  } else {
+    LOG(FATAL) << "Problem file contains markers stored in unknown space.";
+  }
+
+  // Read camera intrinsics.
+  for (int i = 0; i < 8; i++) {
+    camera_intrinsics[i] = file_reader.Read<float>();
+  }
+
+  // Read all cameras.
+  int number_of_cameras = file_reader.Read<int>();
+  for (int i = 0; i < number_of_cameras; i++) {
+    EuclideanCamera camera;
+
+    camera.image = file_reader.Read<int>();
+    ReadMatrix3x3(file_reader, &camera.R);
+    ReadVector3(file_reader, &camera.t);
+
+    if (camera.image >= all_cameras->size()) {
+      all_cameras->resize(camera.image + 1);
+    }
+
+    (*all_cameras)[camera.image].image = camera.image;
+    (*all_cameras)[camera.image].R = camera.R;
+    (*all_cameras)[camera.image].t = camera.t;
+  }
+
+  LOG(INFO) << "Read " << number_of_cameras << " cameras.";
+
+  // Read all reconstructed 3D points.
+  int number_of_points = file_reader.Read<int>();
+  for (int i = 0; i < number_of_points; i++) {
+    EuclideanPoint point;
+
+    point.track = file_reader.Read<int>();
+    ReadVector3(file_reader, &point.X);
+
+    if (point.track >= all_points->size()) {
+      all_points->resize(point.track + 1);
+    }
+
+    (*all_points)[point.track].track = point.track;
+    (*all_points)[point.track].X = point.X;
+  }
+
+  LOG(INFO) << "Read " << number_of_points << " points.";
+
+  // And finally read all markers.
+  int number_of_markers = file_reader.Read<int>();
+  for (int i = 0; i < number_of_markers; i++) {
+    Marker marker;
+
+    marker.image = file_reader.Read<int>();
+    marker.track = file_reader.Read<int>();
+    marker.x = file_reader.Read<float>();
+    marker.y = file_reader.Read<float>();
+
+    all_markers->push_back(marker);
+  }
+
+  LOG(INFO) << "Read " << number_of_markers << " markers.";
+
+  return true;
+}
+
+// Apply camera intrinsics to the normalized point to get image coordinates.
+// This applies the radial lens distortion to a point which is in normalized
+// camera coordinates (i.e. the principal point is at (0, 0)) to get image
+// coordinates in pixels. Templated for use with autodifferentiation.
+template <typename T>
+inline void ApplyRadialDistortionCameraIntrinsics(const T &focal_length_x,
+                                                  const T &focal_length_y,
+                                                  const T &principal_point_x,
+                                                  const T &principal_point_y,
+                                                  const T &k1,
+                                                  const T &k2,
+                                                  const T &k3,
+                                                  const T &p1,
+                                                  const T &p2,
+                                                  const T &normalized_x,
+                                                  const T &normalized_y,
+                                                  T *image_x,
+                                                  T *image_y) {
+  T x = normalized_x;
+  T y = normalized_y;
+
+  // Apply distortion to the normalized points to get (xd, yd).
+  T r2 = x*x + y*y;
+  T r4 = r2 * r2;
+  T r6 = r4 * r2;
+  T r_coeff = (T(1) + k1*r2 + k2*r4 + k3*r6);
+  T xd = x * r_coeff + T(2)*p1*x*y + p2*(r2 + T(2)*x*x);
+  T yd = y * r_coeff + T(2)*p2*x*y + p1*(r2 + T(2)*y*y);
+
+  // Apply focal length and principal point to get the final image coordinates.
+  *image_x = focal_length_x * xd + principal_point_x;
+  *image_y = focal_length_y * yd + principal_point_y;
+}
+
+// Cost functor which computes reprojection error of 3D point X
+// on camera defined by angle-axis rotation and it's translation
+// (which are in the same block due to optimization reasons).
+//
+// This functor uses a radial distortion model.
+struct OpenCVReprojectionError {
+  OpenCVReprojectionError(const double observed_x, const double observed_y)
+      : observed_x(observed_x), observed_y(observed_y) {}
+
+  template <typename T>
+  bool operator()(const T* const intrinsics,
+                  const T* const R_t,  // Rotation denoted by angle axis
+                                       // followed with translation
+                  const T* const X,    // Point coordinates 3x1.
+                  T* residuals) const {
+    // Unpack the intrinsics.
+    const T& focal_length      = intrinsics[OFFSET_FOCAL_LENGTH];
+    const T& principal_point_x = intrinsics[OFFSET_PRINCIPAL_POINT_X];
+    const T& principal_point_y = intrinsics[OFFSET_PRINCIPAL_POINT_Y];
+    const T& k1                = intrinsics[OFFSET_K1];
+    const T& k2                = intrinsics[OFFSET_K2];
+    const T& k3                = intrinsics[OFFSET_K3];
+    const T& p1                = intrinsics[OFFSET_P1];
+    const T& p2                = intrinsics[OFFSET_P2];
+
+    // Compute projective coordinates: x = RX + t.
+    T x[3];
+
+    ceres::AngleAxisRotatePoint(R_t, X, x);
+    x[0] += R_t[3];
+    x[1] += R_t[4];
+    x[2] += R_t[5];
+
+    // Compute normalized coordinates: x /= x[2].
+    T xn = x[0] / x[2];
+    T yn = x[1] / x[2];
+
+    T predicted_x, predicted_y;
+
+    // Apply distortion to the normalized points to get (xd, yd).
+    // TODO(keir): Do early bailouts for zero distortion; these are expensive
+    // jet operations.
+    ApplyRadialDistortionCameraIntrinsics(focal_length,
+                                          focal_length,
+                                          principal_point_x,
+                                          principal_point_y,
+                                          k1, k2, k3,
+                                          p1, p2,
+                                          xn, yn,
+                                          &predicted_x,
+                                          &predicted_y);
+
+    // The error is the difference between the predicted and observed position.
+    residuals[0] = predicted_x - T(observed_x);
+    residuals[1] = predicted_y - T(observed_y);
+
+    return true;
+  }
+
+  const double observed_x;
+  const double observed_y;
+};
+
+// Print a message to the log which camera intrinsics are gonna to be optimized.
+void BundleIntrinsicsLogMessage(const int bundle_intrinsics) {
+  if (bundle_intrinsics == BUNDLE_NO_INTRINSICS) {
+    LOG(INFO) << "Bundling only camera positions.";
+  } else {
+    std::string bundling_message = "";
+
+#define APPEND_BUNDLING_INTRINSICS(name, flag) \
+    if (bundle_intrinsics & flag) { \
+      if (!bundling_message.empty()) { \
+        bundling_message += ", "; \
+      } \
+      bundling_message += name; \
+    } (void)0
+
+    APPEND_BUNDLING_INTRINSICS("f",      BUNDLE_FOCAL_LENGTH);
+    APPEND_BUNDLING_INTRINSICS("px, py", BUNDLE_PRINCIPAL_POINT);
+    APPEND_BUNDLING_INTRINSICS("k1",     BUNDLE_RADIAL_K1);
+    APPEND_BUNDLING_INTRINSICS("k2",     BUNDLE_RADIAL_K2);
+    APPEND_BUNDLING_INTRINSICS("p1",     BUNDLE_TANGENTIAL_P1);
+    APPEND_BUNDLING_INTRINSICS("p2",     BUNDLE_TANGENTIAL_P2);
+
+    LOG(INFO) << "Bundling " << bundling_message << ".";
+  }
+}
+
+// Print a message to the log containing all the camera intriniscs values.
+void PrintCameraIntrinsics(const char *text, const double *camera_intrinsics) {
+  std::ostringstream intrinsics_output;
+
+  intrinsics_output << "f=" << camera_intrinsics[OFFSET_FOCAL_LENGTH];
+
+  intrinsics_output <<
+    " cx=" << camera_intrinsics[OFFSET_PRINCIPAL_POINT_X] <<
+    " cy=" << camera_intrinsics[OFFSET_PRINCIPAL_POINT_Y];
+
+#define APPEND_DISTORTION_COEFFICIENT(name, offset) \
+  { \
+    if (camera_intrinsics[offset] != 0.0) { \
+      intrinsics_output << " " name "=" << camera_intrinsics[offset];  \
+    } \
+  } (void)0
+
+  APPEND_DISTORTION_COEFFICIENT("k1", OFFSET_K1);
+  APPEND_DISTORTION_COEFFICIENT("k2", OFFSET_K2);
+  APPEND_DISTORTION_COEFFICIENT("k3", OFFSET_K3);
+  APPEND_DISTORTION_COEFFICIENT("p1", OFFSET_P1);
+  APPEND_DISTORTION_COEFFICIENT("p2", OFFSET_P2);
+
+#undef APPEND_DISTORTION_COEFFICIENT
+
+  LOG(INFO) << text << intrinsics_output.str();
+}
+
+// Get a vector of camera's rotations denoted by angle axis
+// conjuncted with translations into single block
+//
+// Element with index i matches to a rotation+translation for
+// camera at image i.
+vector<Vec6> PackCamerasRotationAndTranslation(
+    const vector<Marker> &all_markers,
+    const vector<EuclideanCamera> &all_cameras) {
+  vector<Vec6> all_cameras_R_t;
+  int max_image = MaxImage(all_markers);
+
+  all_cameras_R_t.resize(max_image + 1);
+
+  for (int i = 0; i <= max_image; i++) {
+    const EuclideanCamera *camera = CameraForImage(all_cameras, i);
+
+    if (!camera) {
+      continue;
+    }
+
+    ceres::RotationMatrixToAngleAxis(&camera->R(0, 0),
+                                     &all_cameras_R_t[i](0));
+    all_cameras_R_t[i].tail<3>() = camera->t;
+  }
+
+  return all_cameras_R_t;
+}
+
+// Convert cameras rotations fro mangle axis back to rotation matrix.
+void UnpackCamerasRotationAndTranslation(
+    const vector<Marker> &all_markers,
+    const vector<Vec6> &all_cameras_R_t,
+    vector<EuclideanCamera> *all_cameras) {
+  int max_image = MaxImage(all_markers);
+
+  for (int i = 0; i <= max_image; i++) {
+    EuclideanCamera *camera = CameraForImage(all_cameras, i);
+
+    if (!camera) {
+      continue;
+    }
+
+    ceres::AngleAxisToRotationMatrix(&all_cameras_R_t[i](0),
+                                     &camera->R(0, 0));
+    camera->t = all_cameras_R_t[i].tail<3>();
+  }
+}
+
+void EuclideanBundleCommonIntrinsics(const vector<Marker> &all_markers,
+                                     const int bundle_intrinsics,
+                                     const int bundle_constraints,
+                                     double *camera_intrinsics,
+                                     vector<EuclideanCamera> *all_cameras,
+                                     vector<EuclideanPoint> *all_points) {
+  PrintCameraIntrinsics("Original intrinsics: ", camera_intrinsics);
+
+  ceres::Problem::Options problem_options;
+  ceres::Problem problem(problem_options);
+
+  // Convert cameras rotations to angle axis and merge with translation
+  // into single parameter block for maximal minimization speed
+  //
+  // Block for minimization has got the following structure:
+  //   <3 elements for angle-axis> <3 elements for translation>
+  vector<Vec6> all_cameras_R_t =
+    PackCamerasRotationAndTranslation(all_markers, *all_cameras);
+
+  // Parameterization used to restrict camera motion for modal solvers.
+  ceres::SubsetParameterization *constant_transform_parameterization = NULL;
+  if (bundle_constraints & BUNDLE_NO_TRANSLATION) {
+      std::vector<int> constant_translation;
+
+      // First three elements are rotation, last three are translation.
+      constant_translation.push_back(3);
+      constant_translation.push_back(4);
+      constant_translation.push_back(5);
+
+      constant_transform_parameterization =
+        new ceres::SubsetParameterization(6, constant_translation);
+  }
+
+  int num_residuals = 0;
+  bool have_locked_camera = false;
+  for (int i = 0; i < all_markers.size(); ++i) {
+    const Marker &marker = all_markers[i];
+    EuclideanCamera *camera = CameraForImage(all_cameras, marker.image);
+    EuclideanPoint *point = PointForTrack(all_points, marker.track);
+    if (camera == NULL || point == NULL) {
+      continue;
+    }
+
+    // Rotation of camera denoted in angle axis followed with
+    // camera translaiton.
+    double *current_camera_R_t = &all_cameras_R_t[camera->image](0);
+
+    problem.AddResidualBlock(new ceres::AutoDiffCostFunction<
+        OpenCVReprojectionError, 2, 8, 6, 3>(
+            new OpenCVReprojectionError(
+                marker.x,
+                marker.y)),
+        NULL,
+        camera_intrinsics,
+        current_camera_R_t,
+        &point->X(0));
+
+    // We lock the first camera to better deal with scene orientation ambiguity.
+    if (!have_locked_camera) {
+      problem.SetParameterBlockConstant(current_camera_R_t);
+      have_locked_camera = true;
+    }
+
+    if (bundle_constraints & BUNDLE_NO_TRANSLATION) {
+      problem.SetParameterization(current_camera_R_t,
+                                  constant_transform_parameterization);
+    }
+
+    num_residuals++;
+  }
+  LOG(INFO) << "Number of residuals: " << num_residuals;
+
+  if (!num_residuals) {
+    LOG(INFO) << "Skipping running minimizer with zero residuals";
+    return;
+  }
+
+  BundleIntrinsicsLogMessage(bundle_intrinsics);
+
+  if (bundle_intrinsics == BUNDLE_NO_INTRINSICS) {
+    // No camera intrinsics are being refined,
+    // set the whole parameter block as constant for best performance.
+    problem.SetParameterBlockConstant(camera_intrinsics);
+  } else {
+    // Set the camera intrinsics that are not to be bundled as
+    // constant using some macro trickery.
+
+    std::vector<int> constant_intrinsics;
+#define MAYBE_SET_CONSTANT(bundle_enum, offset) \
+    if (!(bundle_intrinsics & bundle_enum)) { \
+      constant_intrinsics.push_back(offset); \
+    }
+    MAYBE_SET_CONSTANT(BUNDLE_FOCAL_LENGTH,    OFFSET_FOCAL_LENGTH);
+    MAYBE_SET_CONSTANT(BUNDLE_PRINCIPAL_POINT, OFFSET_PRINCIPAL_POINT_X);
+    MAYBE_SET_CONSTANT(BUNDLE_PRINCIPAL_POINT, OFFSET_PRINCIPAL_POINT_Y);
+    MAYBE_SET_CONSTANT(BUNDLE_RADIAL_K1,       OFFSET_K1);
+    MAYBE_SET_CONSTANT(BUNDLE_RADIAL_K2,       OFFSET_K2);
+    MAYBE_SET_CONSTANT(BUNDLE_TANGENTIAL_P1,   OFFSET_P1);
+    MAYBE_SET_CONSTANT(BUNDLE_TANGENTIAL_P2,   OFFSET_P2);
+#undef MAYBE_SET_CONSTANT
+
+    // Always set K3 constant, it's not used at the moment.
+    constant_intrinsics.push_back(OFFSET_K3);
+
+    ceres::SubsetParameterization *subset_parameterization =
+      new ceres::SubsetParameterization(8, constant_intrinsics);
+
+    problem.SetParameterization(camera_intrinsics, subset_parameterization);
+  }
+
+  // Configure the solver.
+  ceres::Solver::Options options;
+  options.use_nonmonotonic_steps = true;
+  options.preconditioner_type = ceres::SCHUR_JACOBI;
+  options.linear_solver_type = ceres::ITERATIVE_SCHUR;
+  options.use_inner_iterations = true;
+  options.max_num_iterations = 100;
+  options.minimizer_progress_to_stdout = true;
+
+  // Solve!
+  ceres::Solver::Summary summary;
+  ceres::Solve(options, &problem, &summary);
+
+  std::cout << "Final report:\n" << summary.FullReport();
+
+  // Copy rotations and translations back.
+  UnpackCamerasRotationAndTranslation(all_markers,
+                                      all_cameras_R_t,
+                                      all_cameras);
+
+  PrintCameraIntrinsics("Final intrinsics: ", camera_intrinsics);
+}
+}  // namespace
+
+int main(int argc, char **argv) {
+  google::ParseCommandLineFlags(&argc, &argv, true);
+  google::InitGoogleLogging(argv[0]);
+
+  if (FLAGS_input.empty()) {
+    LOG(ERROR) << "Usage: libmv_bundle_adjuster --input=blender_problem";
+    return EXIT_FAILURE;
+  }
+
+  double camera_intrinsics[8];
+  vector<EuclideanCamera> all_cameras;
+  vector<EuclideanPoint> all_points;
+  bool is_image_space;
+  vector<Marker> all_markers;
+
+  if (!ReadProblemFromFile(FLAGS_input,
+                           camera_intrinsics,
+                           &all_cameras,
+                           &all_points,
+                           &is_image_space,
+                           &all_markers)) {
+    LOG(ERROR) << "Error reading problem file";
+    return EXIT_FAILURE;
+  }
+
+  // If there's no refine_intrinsics passed via command line
+  // (in this case FLAGS_refine_intrinsics will be an empty string)
+  // we use problem's settings to detect whether intrinsics
+  // shall be refined or not.
+  //
+  // Namely, if problem has got markers stored in image (pixel)
+  // space, we do full intrinsics refinement. If markers are
+  // stored in normalized space, and refine_intrinsics is not
+  // set, no refining will happen.
+  //
+  // Using command line argument refine_intrinsics will explicitly
+  // declare which intrinsics need to be refined and in this case
+  // refining flags does not depend on problem at all.
+  int bundle_intrinsics = BUNDLE_NO_INTRINSICS;
+  if (FLAGS_refine_intrinsics.empty()) {
+    if (is_image_space) {
+      bundle_intrinsics = BUNDLE_FOCAL_LENGTH | BUNDLE_RADIAL;
+    }
+  } else {
+    if (FLAGS_refine_intrinsics == "radial") {
+      bundle_intrinsics = BUNDLE_FOCAL_LENGTH | BUNDLE_RADIAL;
+    } else if (FLAGS_refine_intrinsics != "none") {
+      LOG(ERROR) << "Unsupported value for refine-intrinsics";
+      return EXIT_FAILURE;
+    }
+  }
+
+  // Run the bundler.
+  EuclideanBundleCommonIntrinsics(all_markers,
+                                  bundle_intrinsics,
+                                  BUNDLE_NO_CONSTRAINTS,
+                                  camera_intrinsics,
+                                  &all_cameras,
+                                  &all_points);
+
+  return EXIT_SUCCESS;
+}
diff --git a/examples/nist.cc b/examples/nist.cc
index 440ab5c..1773a0f 100644
--- a/examples/nist.cc
+++ b/examples/nist.cc
@@ -28,29 +28,61 @@
 //
 // Author: sameeragarwal@google.com (Sameer Agarwal)
 //
-// NIST non-linear regression problems solved using Ceres.
+// The National Institute of Standards and Technology has released a
+// set of problems to test non-linear least squares solvers.
 //
-// The data was obtained from
-// http://www.itl.nist.gov/div898/strd/nls/nls_main.shtml, where more
-// background on these problems can also be found.
+// More information about the background on these problems and
+// suggested evaluation methodology can be found at:
 //
-// Currently not all problems are solved successfully. Some of the
-// failures are due to convergence to a local minimum, and some fail
-// because of numerical issues.
+//   http://www.itl.nist.gov/div898/strd/nls/nls_info.shtml
 //
-// TODO(sameeragarwal): Fix numerical issues so that all the problems
-// converge and then look at convergence to the wrong solution issues.
+// The problem data themselves can be found at
+//
+//   http://www.itl.nist.gov/div898/strd/nls/nls_main.shtml
+//
+// The problems are divided into three levels of difficulty, Easy,
+// Medium and Hard. For each problem there are two starting guesses,
+// the first one far away from the global minimum and the second
+// closer to it.
+//
+// A problem is considered successfully solved, if every components of
+// the solution matches the globally optimal solution in at least 4
+// digits or more.
+//
+// This dataset was used for an evaluation of Non-linear least squares
+// solvers:
+//
+// P. F. Mondragon & B. Borchers, A Comparison of Nonlinear Regression
+// Codes, Journal of Modern Applied Statistical Methods, 4(1):343-351,
+// 2005.
+//
+// The results from Mondragon & Borchers can be summarized as
+//               Excel  Gnuplot  GaussFit  HBN  MinPack
+// Average LRE     2.3      4.3       4.0  6.8      4.4
+//      Winner       1        5        12   29       12
+//
+// Where the row Winner counts, the number of problems for which the
+// solver had the highest LRE.
+
+// In this file, we implement the same evaluation methodology using
+// Ceres. Currently using Levenberg-Marquard with DENSE_QR, we get
+//
+//               Excel  Gnuplot  GaussFit  HBN  MinPack  Ceres
+// Average LRE     2.3      4.3       4.0  6.8      4.4    9.4
+//      Winner       0        0         5   11        2     41
 
 #include <iostream>
+#include <iterator>
 #include <fstream>
 #include "ceres/ceres.h"
-#include "ceres/split.h"
 #include "gflags/gflags.h"
 #include "glog/logging.h"
 #include "Eigen/Core"
 
 DEFINE_string(nist_data_dir, "", "Directory containing the NIST non-linear"
               "regression examples");
+DEFINE_string(minimizer, "trust_region",
+              "Minimizer type to use, choices are: line_search & trust_region");
 DEFINE_string(trust_region_strategy, "levenberg_marquardt",
               "Options are: levenberg_marquardt, dogleg");
 DEFINE_string(dogleg, "traditional_dogleg",
@@ -60,21 +92,63 @@ DEFINE_string(linear_solver, "dense_qr", "Options are: "
               "cgnr");
 DEFINE_string(preconditioner, "jacobi", "Options are: "
               "identity, jacobi");
+DEFINE_string(line_search, "armijo",
+              "Line search algorithm to use, choices are: armijo and wolfe.");
+DEFINE_string(line_search_direction, "lbfgs",
+              "Line search direction algorithm to use, choices: lbfgs, bfgs");
+DEFINE_int32(max_line_search_iterations, 20,
+             "Maximum number of iterations for each line search.");
+DEFINE_int32(max_line_search_restarts, 10,
+             "Maximum number of restarts of line search direction algorithm.");
+DEFINE_string(line_search_interpolation, "cubic",
+              "Degree of polynomial aproximation in line search, "
+              "choices are: bisection, quadratic & cubic.");
+DEFINE_int32(lbfgs_rank, 20,
+             "Rank of L-BFGS inverse Hessian approximation in line search.");
+DEFINE_bool(approximate_eigenvalue_bfgs_scaling, false,
+            "Use approximate eigenvalue scaling in (L)BFGS line search.");
+DEFINE_double(sufficient_decrease, 1.0e-4,
+              "Line search Armijo sufficient (function) decrease factor.");
+DEFINE_double(sufficient_curvature_decrease, 0.9,
+              "Line search Wolfe sufficient curvature decrease factor.");
 DEFINE_int32(num_iterations, 10000, "Number of iterations");
 DEFINE_bool(nonmonotonic_steps, false, "Trust region algorithm can use"
             " nonmonotic steps");
 DEFINE_double(initial_trust_region_radius, 1e4, "Initial trust region radius");
 
+namespace ceres {
+namespace examples {
+
 using Eigen::Dynamic;
 using Eigen::RowMajor;
 typedef Eigen::Matrix<double, Dynamic, 1> Vector;
 typedef Eigen::Matrix<double, Dynamic, Dynamic, RowMajor> Matrix;
 
+void SplitStringUsingChar(const string& full,
+                          const char delim,
+                          vector<string>* result) {
+  back_insert_iterator< vector<string> > it(*result);
+
+  const char* p = full.data();
+  const char* end = p + full.size();
+  while (p != end) {
+    if (*p == delim) {
+      ++p;
+    } else {
+      const char* start = p;
+      while (++p != end && *p != delim) {
+        // Skip to the next occurence of the delimiter.
+      }
+      *it++ = string(start, p - start);
+    }
+  }
+}
+
 bool GetAndSplitLine(std::ifstream& ifs, std::vector<std::string>* pieces) {
   pieces->clear();
   char buf[256];
   ifs.getline(buf, 256);
-  ceres::SplitStringUsing(std::string(buf), " ", pieces);
+  SplitStringUsingChar(std::string(buf), ' ', pieces);
   return true;
 }
 
@@ -339,7 +413,7 @@ struct Nelson {
 
 template <typename Model, int num_residuals, int num_parameters>
 int RegressionDriver(const std::string& filename,
-                      const ceres::Solver::Options& options) {
+                     const ceres::Solver::Options& options) {
   NISTProblem nist_problem(FLAGS_nist_data_dir + filename);
   CHECK_EQ(num_residuals, nist_problem.response_size());
   CHECK_EQ(num_parameters, nist_problem.num_parameters());
@@ -347,11 +421,12 @@ int RegressionDriver(const std::string& filename,
   Matrix predictor = nist_problem.predictor();
   Matrix response = nist_problem.response();
   Matrix final_parameters = nist_problem.final_parameters();
-  std::vector<ceres::Solver::Summary> summaries(nist_problem.num_starts() + 1);
-  std::cerr << filename << std::endl;
+
+  printf("%s\n", filename.c_str());
 
   // Each NIST problem comes with multiple starting points, so we
   // construct the problem from scratch for each case and solve it.
+  int num_success = 0;
   for (int start = 0; start < nist_problem.num_starts(); ++start) {
     Matrix initial_parameters = nist_problem.initial_parameters(start);
 
@@ -365,43 +440,49 @@ int RegressionDriver(const std::string& filename,
           initial_parameters.data());
     }
 
-    Solve(options, &problem, &summaries[start]);
-  }
-
-  const double certified_cost = nist_problem.certified_cost();
-
-  int num_success = 0;
-  const int kMinNumMatchingDigits = 4;
-  for (int start = 0; start < nist_problem.num_starts(); ++start) {
-    const ceres::Solver::Summary& summary = summaries[start];
-
-    int num_matching_digits = 0;
-    if (IsSuccessfulTermination(summary.termination_type)
-        && summary.final_cost < certified_cost) {
-      num_matching_digits = kMinNumMatchingDigits + 1;
-    } else {
-      num_matching_digits =
-          -std::log10(fabs(summary.final_cost - certified_cost) / certified_cost);
+    ceres::Solver::Summary summary;
+    Solve(options, &problem, &summary);
+
+    // Compute the LRE by comparing each component of the solution
+    // with the ground truth, and taking the minimum.
+    Matrix final_parameters = nist_problem.final_parameters();
+    const double kMaxNumSignificantDigits = 11;
+    double log_relative_error = kMaxNumSignificantDigits + 1;
+    for (int i = 0; i < num_parameters; ++i) {
+      const double tmp_lre =
+          -std::log10(std::fabs(final_parameters(i) - initial_parameters(i)) /
+                      std::fabs(final_parameters(i)));
+      // The maximum LRE is capped at 11 - the precision at which the
+      // ground truth is known.
+      //
+      // The minimum LRE is capped at 0 - no digits match between the
+      // computed solution and the ground truth.
+      log_relative_error =
+          std::min(log_relative_error,
+                   std::max(0.0, std::min(kMaxNumSignificantDigits, tmp_lre)));
     }
 
-    std::cerr << "start " << start + 1 << " " ;
-    if (num_matching_digits <= kMinNumMatchingDigits) {
-      std::cerr <<  "FAILURE";
-    } else {
-      std::cerr <<  "SUCCESS";
+    const int kMinNumMatchingDigits = 4;
+    if (log_relative_error >= kMinNumMatchingDigits) {
       ++num_success;
     }
-    std::cerr << " summary: "
-              << summary.BriefReport()
-              << " Certified cost: " << certified_cost
-              << std::endl;
 
+    printf("start: %d status: %s lre: %4.1f initial cost: %e final cost:%e "
+           "certified cost: %e total iterations: %d\n",
+           start + 1,
+           log_relative_error < kMinNumMatchingDigits ? "FAILURE" : "SUCCESS",
+           log_relative_error,
+           summary.initial_cost,
+           summary.final_cost,
+           nist_problem.certified_cost(),
+           (summary.num_successful_steps + summary.num_unsuccessful_steps));
   }
-
   return num_success;
 }
 
 void SetMinimizerOptions(ceres::Solver::Options* options) {
+  CHECK(ceres::StringToMinimizerType(FLAGS_minimizer,
+                                     &options->minimizer_type));
   CHECK(ceres::StringToLinearSolverType(FLAGS_linear_solver,
                                         &options->linear_solver_type));
   CHECK(ceres::StringToPreconditionerType(FLAGS_preconditioner,
@@ -410,10 +491,28 @@ void SetMinimizerOptions(ceres::Solver::Options* options) {
             FLAGS_trust_region_strategy,
             &options->trust_region_strategy_type));
   CHECK(ceres::StringToDoglegType(FLAGS_dogleg, &options->dogleg_type));
+  CHECK(ceres::StringToLineSearchDirectionType(
+      FLAGS_line_search_direction,
+      &options->line_search_direction_type));
+  CHECK(ceres::StringToLineSearchType(FLAGS_line_search,
+                                      &options->line_search_type));
+  CHECK(ceres::StringToLineSearchInterpolationType(
+      FLAGS_line_search_interpolation,
+      &options->line_search_interpolation_type));
 
   options->max_num_iterations = FLAGS_num_iterations;
   options->use_nonmonotonic_steps = FLAGS_nonmonotonic_steps;
   options->initial_trust_region_radius = FLAGS_initial_trust_region_radius;
+  options->max_lbfgs_rank = FLAGS_lbfgs_rank;
+  options->line_search_sufficient_function_decrease = FLAGS_sufficient_decrease;
+  options->line_search_sufficient_curvature_decrease =
+      FLAGS_sufficient_curvature_decrease;
+  options->max_num_line_search_step_size_iterations =
+      FLAGS_max_line_search_iterations;
+  options->max_num_line_search_direction_restarts =
+      FLAGS_max_line_search_restarts;
+  options->use_approximate_eigenvalue_bfgs_scaling =
+      FLAGS_approximate_eigenvalue_bfgs_scaling;
   options->function_tolerance = 1e-18;
   options->gradient_tolerance = 1e-18;
   options->parameter_tolerance = 1e-18;
@@ -427,7 +526,7 @@ void SolveNISTProblems() {
   ceres::Solver::Options options;
   SetMinimizerOptions(&options);
 
-  std::cerr << "Lower Difficulty\n";
+  std::cout << "Lower Difficulty\n";
   int easy_success = 0;
   easy_success += RegressionDriver<Misra1a,  1, 2>("Misra1a.dat",  options);
   easy_success += RegressionDriver<Chwirut,  1, 3>("Chwirut1.dat", options);
@@ -438,7 +537,7 @@ void SolveNISTProblems() {
   easy_success += RegressionDriver<DanWood,  1, 2>("DanWood.dat",  options);
   easy_success += RegressionDriver<Misra1b,  1, 2>("Misra1b.dat",  options);
 
-  std::cerr << "\nMedium Difficulty\n";
+  std::cout << "\nMedium Difficulty\n";
   int medium_success = 0;
   medium_success += RegressionDriver<Kirby2,   1, 5>("Kirby2.dat",   options);
   medium_success += RegressionDriver<Hahn1,    1, 7>("Hahn1.dat",    options);
@@ -452,7 +551,7 @@ void SolveNISTProblems() {
   medium_success += RegressionDriver<Roszman1, 1, 4>("Roszman1.dat", options);
   medium_success += RegressionDriver<ENSO,     1, 9>("ENSO.dat",     options);
 
-  std::cerr << "\nHigher Difficulty\n";
+  std::cout << "\nHigher Difficulty\n";
   int hard_success = 0;
   hard_success += RegressionDriver<MGH09,    1, 4>("MGH09.dat",    options);
   hard_success += RegressionDriver<Thurber,  1, 7>("Thurber.dat",  options);
@@ -464,16 +563,19 @@ void SolveNISTProblems() {
   hard_success += RegressionDriver<Rat43,    1, 4>("Rat43.dat",    options);
   hard_success += RegressionDriver<Bennet5,  1, 3>("Bennett5.dat", options);
 
-  std::cerr << "\n";
-  std::cerr << "Easy    : " << easy_success << "/16\n";
-  std::cerr << "Medium  : " << medium_success << "/22\n";
-  std::cerr << "Hard    : " << hard_success << "/16\n";
-  std::cerr << "Total   : " << easy_success + medium_success + hard_success << "/54\n";
+  std::cout << "\n";
+  std::cout << "Easy    : " << easy_success << "/16\n";
+  std::cout << "Medium  : " << medium_success << "/22\n";
+  std::cout << "Hard    : " << hard_success << "/16\n";
+  std::cout << "Total   : " << easy_success + medium_success + hard_success << "/54\n";
 }
 
+}  // namespace examples
+}  // namespace ceres
+
 int main(int argc, char** argv) {
   google::ParseCommandLineFlags(&argc, &argv, true);
   google::InitGoogleLogging(argv[0]);
-  SolveNISTProblems();
+  ceres::examples::SolveNISTProblems();
   return 0;
 };
diff --git a/examples/powell.cc b/examples/powell.cc
index 6cd3611..c0cba02 100644
--- a/examples/powell.cc
+++ b/examples/powell.cc
@@ -55,8 +55,7 @@ using ceres::Problem;
 using ceres::Solver;
 using ceres::Solve;
 
-class F1 {
- public:
+struct F1 {
   template <typename T> bool operator()(const T* const x1,
                                         const T* const x2,
                                         T* residual) const {
@@ -66,8 +65,7 @@ class F1 {
   }
 };
 
-class F2 {
- public:
+struct F2 {
   template <typename T> bool operator()(const T* const x3,
                                         const T* const x4,
                                         T* residual) const {
@@ -77,8 +75,7 @@ class F2 {
   }
 };
 
-class F3 {
- public:
+struct F3 {
   template <typename T> bool operator()(const T* const x2,
                                         const T* const x4,
                                         T* residual) const {
@@ -88,8 +85,7 @@ class F3 {
   }
 };
 
-class F4 {
- public:
+struct F4 {
   template <typename T> bool operator()(const T* const x1,
                                         const T* const x4,
                                         T* residual) const {
@@ -99,6 +95,9 @@ class F4 {
   }
 };
 
+DEFINE_string(minimizer, "trust_region",
+              "Minimizer type to use, choices are: line_search & trust_region");
+
 int main(int argc, char** argv) {
   google::ParseCommandLineFlags(&argc, &argv, true);
   google::InitGoogleLogging(argv[0]);
@@ -125,23 +124,27 @@ int main(int argc, char** argv) {
                            NULL,
                            &x1, &x4);
 
-  // Run the solver!
   Solver::Options options;
-  options.max_num_iterations = 30;
+  LOG_IF(FATAL, !ceres::StringToMinimizerType(FLAGS_minimizer,
+                                              &options.minimizer_type))
+      << "Invalid minimizer: " << FLAGS_minimizer
+      << ", valid options are: trust_region and line_search.";
+
+  options.max_num_iterations = 100;
   options.linear_solver_type = ceres::DENSE_QR;
   options.minimizer_progress_to_stdout = true;
 
-  Solver::Summary summary;
-
   std::cout << "Initial x1 = " << x1
             << ", x2 = " << x2
             << ", x3 = " << x3
             << ", x4 = " << x4
             << "\n";
 
+  // Run the solver!
+  Solver::Summary summary;
   Solve(options, &problem, &summary);
 
-  std::cout << summary.BriefReport() << "\n";
+  std::cout << summary.FullReport() << "\n";
   std::cout << "Final x1 = " << x1
             << ", x2 = " << x2
             << ", x3 = " << x3
diff --git a/examples/quadratic_auto_diff.cc b/examples/quadratic_auto_diff.cc
index ea7fae9..1e2f3ef 100644
--- a/examples/quadratic_auto_diff.cc
+++ b/examples/quadratic_auto_diff.cc
@@ -44,10 +44,11 @@ using ceres::Problem;
 using ceres::Solver;
 using ceres::Solve;
 
-// A templated cost function that implements the residual r = 10 - x. The method
-// Map is templated so that we can then use an automatic differentiation wrapper
-// around it to generate its derivatives.
-class QuadraticCostFunction {
+// 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];
@@ -69,8 +70,8 @@ int main(int argc, char** argv) {
   // Set up the only cost function (also known as residual). This uses
   // auto-differentiation to obtain the derivative (jacobian).
   problem.AddResidualBlock(
-      new AutoDiffCostFunction<QuadraticCostFunction, 1, 1>(
-          new QuadraticCostFunction),
+      new AutoDiffCostFunction<QuadraticCostFunctor, 1, 1>(
+          new QuadraticCostFunctor),
       NULL,
       &x);
 
diff --git a/examples/quadratic_numeric_diff.cc b/examples/quadratic_numeric_diff.cc
index 8ec88ef..1082616 100644
--- a/examples/quadratic_numeric_diff.cc
+++ b/examples/quadratic_numeric_diff.cc
@@ -38,24 +38,16 @@
 
 using ceres::NumericDiffCostFunction;
 using ceres::CENTRAL;
-using ceres::SizedCostFunction;
 using ceres::CostFunction;
 using ceres::Problem;
 using ceres::Solver;
 using ceres::Solve;
 
-class ResidualWithNoDerivative
-  : public SizedCostFunction<1 /* number of residuals */,
-                             1 /* size of first parameter */> {
+// A cost functor that implements the residual r = 10 - x.
+class QuadraticCostFunctor {
  public:
-  virtual ~ResidualWithNoDerivative() {}
-  virtual bool Evaluate(double const* const* parameters,
-                        double* residuals,
-                        double** jacobians) const {
-    (void) jacobians;  // Ignored; filled in by numeric differentiation.
-
-    // f(x) = 10 - x.
-    residuals[0] = 10 - parameters[0][0];
+  bool operator()(const double* const x, double* residual) const {
+    residual[0] = 10.0 - x[0];
     return true;
   }
 };
@@ -71,8 +63,8 @@ int main(int argc, char** argv) {
   // Set up the only cost function (also known as residual). This uses
   // numeric differentiation to obtain the derivative (jacobian).
   CostFunction* cost =
-      new NumericDiffCostFunction<ResidualWithNoDerivative, CENTRAL, 1, 1> (
-          new ResidualWithNoDerivative, ceres::TAKE_OWNERSHIP);
+      new NumericDiffCostFunction<QuadraticCostFunctor, CENTRAL, 1, 1> (
+          new QuadraticCostFunctor);
 
   // Build the problem.
   Problem problem;
diff --git a/examples/robust_curve_fitting.cc b/examples/robust_curve_fitting.cc
new file mode 100644
index 0000000..ad71a5d
--- /dev/null
+++ b/examples/robust_curve_fitting.cc
@@ -0,0 +1,163 @@
+// 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 "glog/logging.h"
+
+// 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;
+//   outlier_noise = rand(size(x)) < 0.05;
+//   y_observed = y + noise + outlier_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, 5.543078e+00, // Outlier point
+1.500000e+00, 5.664015e+00, // Outlier point
+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
+};
+
+using ceres::AutoDiffCostFunction;
+using ceres::CostFunction;
+using ceres::CauchyLoss;
+using ceres::Problem;
+using ceres::Solve;
+using ceres::Solver;
+
+struct ExponentialResidual {
+  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::InitGoogleLogging(argv[0]);
+
+  double m = 0.0;
+  double c = 0.0;
+
+  Problem problem;
+  for (int i = 0; i < kNumObservations; ++i) {
+    CostFunction* cost_function =
+        new AutoDiffCostFunction<ExponentialResidual, 1, 1, 1>(
+            new ExponentialResidual(data[2 * i], data[2 * i + 1]));
+    problem.AddResidualBlock(cost_function, NULL, &m, &c);
+  }
+
+  Solver::Options options;
+  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/simple_bundle_adjuster.cc b/examples/simple_bundle_adjuster.cc
index cc6f04a..736d4a2 100644
--- a/examples/simple_bundle_adjuster.cc
+++ b/examples/simple_bundle_adjuster.cc
@@ -160,6 +160,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;
 };
@@ -177,6 +185,8 @@ int main(int argc, char** argv) {
     return 1;
   }
 
+  const double* observations = bal_problem.observations();
+
   // Create residuals for each observation in the bundle adjustment problem. The
   // parameters for cameras and points are added automatically.
   ceres::Problem problem;
@@ -184,12 +194,10 @@ int main(int argc, char** argv) {
     // Each Residual block takes a point and a camera as input and outputs a 2
     // dimensional residual. Internally, the cost function stores the observed
     // image location and compares the reprojection against the observation.
-    ceres::CostFunction* cost_function =
-        new ceres::AutoDiffCostFunction<SnavelyReprojectionError, 2, 9, 3>(
-            new SnavelyReprojectionError(
-                bal_problem.observations()[2 * i + 0],
-                bal_problem.observations()[2 * i + 1]));
 
+    ceres::CostFunction* cost_function =
+        SnavelyReprojectionError::Create(observations[2 * i + 0],
+                                         observations[2 * i + 1]);
     problem.AddResidualBlock(cost_function,
                              NULL /* squared loss */,
                              bal_problem.mutable_camera_for_observation(i),