diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt
new file mode 100644
index 0000000..b07a3dd
--- /dev/null
+++ b/examples/CMakeLists.txt
@@ -0,0 +1,53 @@
+# 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)
+
+ADD_EXECUTABLE(quadratic quadratic.cc)
+TARGET_LINK_LIBRARIES(quadratic 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(powell powell.cc)
+TARGET_LINK_LIBRARIES(powell ceres)
+
+ADD_EXECUTABLE(circle_fit circle_fit.cc)
+TARGET_LINK_LIBRARIES(circle_fit ceres)
+
+ADD_EXECUTABLE(bundle_adjuster
+               bundle_adjuster.cc
+               bal_problem.cc)
+TARGET_LINK_LIBRARIES(bundle_adjuster ceres)
+
+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
new file mode 100644
index 0000000..2ab35e4
--- /dev/null
+++ b/examples/bal_problem.cc
@@ -0,0 +1,119 @@
+// 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 "bal_problem.h"
+
+#include <cstdio>
+#include <cstdlib>
+#include <string>
+#include <glog/logging.h>
+#include "ceres/rotation.h"
+
+namespace ceres {
+namespace examples {
+
+template<typename T>
+void FscanfOrDie(FILE *fptr, const char *format, T *value) {
+  int num_scanned = fscanf(fptr, format, value);
+  if (num_scanned != 1) {
+    LOG(FATAL) << "Invalid UW data file.";
+  }
+}
+
+BALProblem::BALProblem(const std::string filename, bool use_quaternions) {
+  FILE* fptr = fopen(filename.c_str(), "r");
+
+  if (!fptr) {
+    LOG(FATAL) << "Error: unable to open file " << filename;
+    return;
+  };
+
+  // This wil die horribly on invalid files. Them's the breaks.
+  FscanfOrDie(fptr, "%d", &num_cameras_);
+  FscanfOrDie(fptr, "%d", &num_points_);
+  FscanfOrDie(fptr, "%d", &num_observations_);
+
+  VLOG(1) << "Header: " << num_cameras_
+          << " " << num_points_
+          << " " << num_observations_;
+
+  point_index_ = new int[num_observations_];
+  camera_index_ = new int[num_observations_];
+  observations_ = new double[2 * num_observations_];
+
+  num_parameters_ = 9 * num_cameras_ + 3 * num_points_;
+  parameters_ = new double[num_parameters_];
+
+  for (int i = 0; i < num_observations_; ++i) {
+    FscanfOrDie(fptr, "%d", camera_index_ + i);
+    FscanfOrDie(fptr, "%d", point_index_ + i);
+    for (int j = 0; j < 2; ++j) {
+      FscanfOrDie(fptr, "%lf", observations_ + 2*i + j);
+    }
+  }
+
+  for (int i = 0; i < num_parameters_; ++i) {
+    FscanfOrDie(fptr, "%lf", parameters_ + i);
+  }
+
+  use_quaternions_ = use_quaternions;
+  if (use_quaternions) {
+    // Switch the angle-axis rotations to quaternions.
+    num_parameters_ = 10 * num_cameras_ + 3 * num_points_;
+    double* quaternion_parameters = new double[num_parameters_];
+    double* original_cursor = parameters_;
+    double* quaternion_cursor = quaternion_parameters;
+    for (int i = 0; i < num_cameras_; ++i) {
+      AngleAxisToQuaternion(original_cursor, quaternion_cursor);
+      quaternion_cursor += 4;
+      original_cursor += 3;
+      for (int j = 4; j < 10; ++j) {
+       *quaternion_cursor++ = *original_cursor++;
+      }
+    }
+    // Copy the rest of the points.
+    for (int i = 0; i < 3 * num_points_; ++i) {
+      *quaternion_cursor++ = *original_cursor++;
+    }
+    // Swap in the quaternion parameters.
+    delete []parameters_;
+    parameters_ = quaternion_parameters;
+  }
+}
+
+BALProblem::~BALProblem() {
+  delete []point_index_;
+  delete []camera_index_;
+  delete []observations_;
+  delete []parameters_;
+}
+
+}  // namespace examples
+}  // namespace ceres
diff --git a/examples/bal_problem.h b/examples/bal_problem.h
new file mode 100644
index 0000000..ff411b8
--- /dev/null
+++ b/examples/bal_problem.h
@@ -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: sameeragarwal@google.com (Sameer Agarwal)
+//
+// Class for loading and holding in memory bundle adjustment problems
+// from the BAL (Bundle Adjustment in the Large) dataset from the
+// University of Washington.
+//
+// For more details see http://grail.cs.washington.edu/projects/bal/
+
+#ifndef CERES_EXAMPLES_BAL_PROBLEM_H_
+#define CERES_EXAMPLES_BAL_PROBLEM_H_
+
+#include <string>
+
+namespace ceres {
+namespace examples {
+
+class BALProblem {
+ public:
+  explicit BALProblem(const std::string filename, bool use_quaternions);
+  ~BALProblem();
+
+  int camera_block_size()      const { return use_quaternions_ ? 10 : 9; }
+  int point_block_size()       const { return 3;                         }
+  int num_cameras()            const { return num_cameras_;              }
+  int num_points()             const { return num_points_;               }
+  int num_observations()       const { return num_observations_;         }
+  int num_parameters()         const { return num_parameters_;           }
+  const int* point_index()     const { return point_index_;              }
+  const int* camera_index()    const { return camera_index_;             }
+  const double* observations() const { return observations_;             }
+  const double* parameters()   const { return parameters_;               }
+  double* mutable_cameras()          { return parameters_;               }
+  double* mutable_points() {
+    return parameters_  + camera_block_size() * num_cameras_;
+  }
+
+ private:
+  int num_cameras_;
+  int num_points_;
+  int num_observations_;
+  int num_parameters_;
+  bool use_quaternions_;
+
+  int* point_index_;
+  int* camera_index_;
+  double* observations_;
+  // The parameter vector is laid out as follows
+  // [camera_1, ..., camera_n, point_1, ..., point_m]
+  double* parameters_;
+};
+
+}  // namespace examples
+}  // namespace ceres
+
+#endif  // CERES_EXAMPLES_BAL_PROBLEM_H_
diff --git a/examples/bundle_adjuster.cc b/examples/bundle_adjuster.cc
new file mode 100644
index 0000000..bafdf25
--- /dev/null
+++ b/examples/bundle_adjuster.cc
@@ -0,0 +1,281 @@
+// 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)
+//
+// An example of solving a dynamically sized problem with various
+// solvers and loss functions.
+//
+// For a simpler bare bones example of doing bundle adjustment with
+// Ceres, please see simple_bundle_adjuster.cc.
+//
+// NOTE: This example will not compile without gflags and SuiteSparse.
+//
+// The problem being solved here is known as a Bundle Adjustment
+// problem in computer vision. Given a set of 3d points X_1, ..., X_n,
+// a set of cameras P_1, ..., P_m. If the point X_i is visible in
+// image j, then there is a 2D observation u_ij that is the expected
+// projection of X_i using P_j. The aim of this optimization is to
+// find values of X_i and P_j such that the reprojection error
+//
+//    E(X,P) =  sum_ij  |u_ij - P_j X_i|^2
+//
+// is minimized.
+//
+// The problem used here comes from a collection of bundle adjustment
+// problems published at University of Washington.
+// http://grail.cs.washington.edu/projects/bal
+
+#include <algorithm>
+#include <cmath>
+#include <cstdio>
+#include <string>
+#include <vector>
+
+#include <gflags/gflags.h>
+#include <glog/logging.h>
+#include "bal_problem.h"
+#include "snavely_reprojection_error.h"
+#include "ceres/ceres.h"
+
+DEFINE_string(input, "", "Input File name");
+
+DEFINE_string(solver_type, "sparse_schur", "Options are: "
+              "sparse_schur, dense_schur, iterative_schur, cholesky "
+              "and dense_qr");
+
+DEFINE_string(preconditioner_type, "jacobi", "Options are: "
+              "identity, jacobi, schur_jacobi, cluster_jacobi, "
+              "cluster_tridiagonal");
+
+DEFINE_int32(num_iterations, 5, "Number of iterations");
+DEFINE_int32(num_threads, 1, "Number of threads");
+DEFINE_bool(use_schur_ordering, false, "Use automatic Schur ordering.");
+DEFINE_bool(use_quaternions, false, "If true, uses quaternions to represent "
+            "rotations. If false, angle axis is used");
+DEFINE_bool(use_local_parameterization, false, "For quaternions, use a local "
+            "parameterization.");
+DEFINE_bool(robustify, false, "Use a robust loss function");
+
+namespace ceres {
+namespace examples {
+
+void SetLinearSolver(Solver::Options* options) {
+  if (FLAGS_solver_type == "sparse_schur") {
+    options->linear_solver_type = ceres::SPARSE_SCHUR;
+  } else if (FLAGS_solver_type == "dense_schur") {
+    options->linear_solver_type = ceres::DENSE_SCHUR;
+  } else if (FLAGS_solver_type == "iterative_schur") {
+    options->linear_solver_type = ceres::ITERATIVE_SCHUR;
+  } else if (FLAGS_solver_type == "cholesky") {
+    options->linear_solver_type = ceres::SPARSE_NORMAL_CHOLESKY;
+  } else if (FLAGS_solver_type == "dense_qr") {
+    // DENSE_QR is included here for completeness, but actually using
+    // this opttion is a bad idea due to the amount of memory needed
+    // to store even the smallest of the bundle adjustment jacobian
+    // arrays
+    options->linear_solver_type = ceres::DENSE_QR;
+  } else {
+    LOG(FATAL) << "Unknown ceres solver type: "
+               << FLAGS_solver_type;
+  }
+
+  if (options->linear_solver_type == ceres::ITERATIVE_SCHUR) {
+    options->linear_solver_min_num_iterations = 5;
+    if (FLAGS_preconditioner_type == "identity") {
+      options->preconditioner_type = ceres::IDENTITY;
+    } else if (FLAGS_preconditioner_type == "jacobi") {
+      options->preconditioner_type = ceres::JACOBI;
+    } else if (FLAGS_preconditioner_type == "schur_jacobi") {
+      options->preconditioner_type = ceres::SCHUR_JACOBI;
+    } else if (FLAGS_preconditioner_type == "cluster_jacobi") {
+      options->preconditioner_type = ceres::CLUSTER_JACOBI;
+    } else if (FLAGS_preconditioner_type == "cluster_tridiagonal") {
+      options->preconditioner_type = ceres::CLUSTER_TRIDIAGONAL;
+    } else {
+      LOG(FATAL) << "Unknown ceres preconditioner type: "
+               << FLAGS_preconditioner_type;
+    }
+  }
+
+  options->num_linear_solver_threads = FLAGS_num_threads;
+}
+
+void SetOrdering(BALProblem* bal_problem, Solver::Options* options) {
+  // Bundle adjustment problems have a sparsity structure that makes
+  // them amenable to more specialized and much more efficient
+  // solution strategies. The SPARSE_SCHUR, DENSE_SCHUR and
+  // ITERATIVE_SCHUR solvers make use of this specialized
+  // structure. Using them however requires that the ParameterBlocks
+  // are in a particular order (points before cameras) and
+  // Solver::Options::num_eliminate_blocks is set to the number of
+  // points.
+  //
+  // This can either be done by specifying Options::ordering_type =
+  // ceres::SCHUR, in which case Ceres will automatically determine
+  // the right ParameterBlock ordering, or by manually specifying a
+  // suitable ordering vector and defining
+  // Options::num_eliminate_blocks.
+  if (FLAGS_use_schur_ordering) {
+    options->ordering_type = ceres::SCHUR;
+    return;
+  }
+
+  options->ordering_type = ceres::USER;
+  const int num_points = bal_problem->num_points();
+  const int point_block_size = bal_problem->point_block_size();
+  double* points = bal_problem->mutable_points();
+  const int num_cameras = bal_problem->num_cameras();
+  const int camera_block_size = bal_problem->camera_block_size();
+  double* cameras = bal_problem->mutable_cameras();
+
+  // The points come before the cameras.
+  for (int i = 0; i < num_points; ++i) {
+    options->ordering.push_back(points + point_block_size * i);
+  }
+
+  for (int i = 0; i < num_cameras; ++i) {
+    // When using axis-angle, there is a single parameter block for
+    // the entire camera.
+    options->ordering.push_back(cameras + camera_block_size * i);
+
+    // If quaternions are used, there are two blocks, so add the
+    // second block to the ordering.
+    if (FLAGS_use_quaternions) {
+      options->ordering.push_back(cameras + camera_block_size * i + 4);
+    }
+  }
+
+  options->num_eliminate_blocks = num_points;
+}
+
+void SetMinimizerOptions(Solver::Options* options) {
+  options->max_num_iterations = FLAGS_num_iterations;
+  options->minimizer_progress_to_stdout = true;
+  options->num_threads = FLAGS_num_threads;
+}
+
+void SetSolverOptionsFromFlags(BALProblem* bal_problem,
+                               Solver::Options* options) {
+  SetLinearSolver(options);
+  SetOrdering(bal_problem, options);
+  SetMinimizerOptions(options);
+}
+
+void BuildProblem(BALProblem* bal_problem, Problem* problem) {
+  const int point_block_size = bal_problem->point_block_size();
+  const int camera_block_size = bal_problem->camera_block_size();
+  double* points = bal_problem->mutable_points();
+  double* cameras = bal_problem->mutable_cameras();
+
+  // Observations is 2*num_observations long array observations =
+  // [u_1, u_2, ... , u_n], where each u_i is two dimensional, the x
+  // and y positions of the observation.
+  const double* observations = bal_problem->observations();
+
+  for (int i = 0; i < bal_problem->num_observations(); ++i) {
+    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<
+          SnavelyReprojectionErrorWitQuaternions, 2, 4, 6, 3>(
+              new SnavelyReprojectionErrorWitQuaternions(
+                  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]));
+    }
+
+    // If enabled use Huber's loss function.
+    LossFunction* loss_function = FLAGS_robustify ? new HuberLoss(1.0) : NULL;
+
+    // Each observation correponds to a pair of a camera and a point
+    // which are identified by camera_index()[i] and point_index()[i]
+    // respectively.
+    double* camera =
+        cameras + camera_block_size * bal_problem->camera_index()[i];
+    double* point = points + point_block_size * bal_problem->point_index()[i];
+
+    if (FLAGS_use_quaternions) {
+      // When using quaternions, we split the camera into two
+      // parameter blocks. One of size 4 for the quaternion and the
+      // other of size 6 containing the translation, focal length and
+      // the radial distortion parameters.
+      problem->AddResidualBlock(cost_function,
+                                loss_function,
+                                camera,
+                                camera + 4,
+                                point);
+    } else {
+      problem->AddResidualBlock(cost_function, loss_function, camera, point);
+    }
+  }
+
+  if (FLAGS_use_quaternions && FLAGS_use_local_parameterization) {
+    LocalParameterization* quaternion_parameterization =
+         new QuaternionParameterization;
+    for (int i = 0; i < bal_problem->num_cameras(); ++i) {
+      problem->SetParameterization(cameras + camera_block_size * i,
+                                   quaternion_parameterization);
+    }
+  }
+}
+
+void SolveProblem(const char* filename) {
+  BALProblem bal_problem(filename, FLAGS_use_quaternions);
+  Problem problem;
+  BuildProblem(&bal_problem, &problem);
+  Solver::Options options;
+  SetSolverOptionsFromFlags(&bal_problem, &options);
+
+  Solver::Summary summary;
+  Solve(options, &problem, &summary);
+  std::cout << summary.FullReport() << "\n";
+}
+
+}  // namespace examples
+}  // namespace ceres
+
+int main(int argc, char** argv) {
+  google::ParseCommandLineFlags(&argc, &argv, true);
+  google::InitGoogleLogging(argv[0]);
+  if (FLAGS_input.empty()) {
+    LOG(ERROR) << "Usage: bundle_adjustment_example --input=bal_problem";
+    return 1;
+  }
+
+  CHECK(FLAGS_use_quaternions || !FLAGS_use_local_parameterization)
+      << "--use_local_parameterization can only be used with "
+      << "--use_quaternions.";
+  ceres::examples::SolveProblem(FLAGS_input.c_str());
+  return 0;
+}
diff --git a/examples/circle_fit.cc b/examples/circle_fit.cc
new file mode 100644
index 0000000..a044134
--- /dev/null
+++ b/examples/circle_fit.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: keir@google.com (Keir Mierle)
+//
+// This fits circles to a collection of points, where the error is related to
+// the distance of a point from the circle. This uses auto-differentiation to
+// take the derivatives.
+//
+// The input format is simple text. Feed on standard in:
+//
+//   x_initial y_initial r_initial
+//   x1 y1
+//   x2 y2
+//   y3 y3
+//   ...
+//
+// And the result after solving will be printed to stdout:
+//
+//   x y r
+//
+// There are closed form solutions [1] to this problem which you may want to
+// consider instead of using this one. If you already have a decent guess, Ceres
+// can squeeze down the last bit of error.
+//
+//   [1] http://www.mathworks.com/matlabcentral/fileexchange/5557-circle-fit/content/circfit.m
+
+#include <cstdio>
+#include <vector>
+
+#include <gflags/gflags.h>
+#include "ceres/ceres.h"
+
+using ceres::AutoDiffCostFunction;
+using ceres::CauchyLoss;
+using ceres::CostFunction;
+using ceres::LossFunction;
+using ceres::Problem;
+using ceres::Solve;
+using ceres::Solver;
+
+DEFINE_double(robust_threshold, -1.0, "Robust loss parameter. Set to -1 for "
+              "normal squared error (no robustification).");
+
+// The cost for a single sample. The returned residual is related to the
+// distance of the point from the circle (passed in as x, y, m parameters).
+//
+// Note that the radius is parameterized as r = m^2 to constrain the radius to
+// positive values.
+class DistanceFromCircleCost {
+ public:
+  DistanceFromCircleCost(double xx, double yy) : xx_(xx), yy_(yy) {}
+  template <typename T> bool operator()(const T* const x,
+                                        const T* const y,
+                                        const T* const m,  // r = m^2
+                                        T* residual) const {
+    // Since the radius is parameterized as m^2, unpack m to get r.
+    T r = *m * *m;
+
+    // Get the position of the sample in the circle's coordinate system.
+    T xp = xx_ - *x;
+    T yp = yy_ - *y;
+
+    // It is tempting to use the following cost:
+    //
+    //   residual[0] = *r - sqrt(xp*xp + yp*yp);
+    //
+    // which is the distance of the sample from the circle. This works
+    // reasonably well, but the sqrt() adds strong nonlinearities to the cost
+    // function. Instead, a different cost is used, which while not strictly a
+    // distance in the metric sense (it has units distance^2) it produces more
+    // robust fits when there are outliers. This is because the cost surface is
+    // more convex.
+    residual[0] = r*r - xp*xp - yp*yp;
+    return true;
+  }
+
+ private:
+  // The measured x,y coordinate that should be on the circle.
+  double xx_, yy_;
+};
+
+int main(int argc, char** argv) {
+  google::ParseCommandLineFlags(&argc, &argv, true);
+  google::InitGoogleLogging(argv[0]);
+
+  double x, y, r;
+  if (scanf("%lg %lg %lg", &x, &y, &r) != 3) {
+    fprintf(stderr, "Couldn't read first line.\n");
+    return 1;
+  }
+  fprintf(stderr, "Got x, y, r %lg, %lg, %lg\n", x, y, r);
+
+  // Save initial values for comparison.
+  double initial_x = x;
+  double initial_y = y;
+  double initial_r = r;
+
+  // Parameterize r as m^2 so that it can't be negative.
+  double m = sqrt(r);
+
+  Problem problem;
+
+  // Configure the loss function.
+  LossFunction* loss = NULL;
+  if (FLAGS_robust_threshold) {
+    loss = new CauchyLoss(FLAGS_robust_threshold);
+  }
+
+  // Add the residuals.
+  double xx, yy;
+  int num_points = 0;
+  while (scanf("%lf %lf\n", &xx, &yy) == 2) {
+    CostFunction *cost =
+        new AutoDiffCostFunction<DistanceFromCircleCost, 1, 1, 1, 1>(
+            new DistanceFromCircleCost(xx, yy));
+    problem.AddResidualBlock(cost, loss, &x, &y, &m);
+    num_points++;
+  }
+
+  std::cout << "Got " << num_points << " points.\n";
+
+  // Build and solve the problem.
+  Solver::Options options;
+  options.max_num_iterations = 500;
+  options.linear_solver_type = ceres::DENSE_QR;
+  Solver::Summary summary;
+  Solve(options, &problem, &summary);
+
+  // Recover r from m.
+  r = m * m;
+
+  std::cout << summary.BriefReport() << "\n";
+  std::cout << "x : " << initial_x << " -> " << x << "\n";
+  std::cout << "y : " << initial_y << " -> " << y << "\n";
+  std::cout << "r : " << initial_r << " -> " << r << "\n";
+  return 0;
+}
diff --git a/examples/data_fitting.cc b/examples/data_fitting.cc
new file mode 100644
index 0000000..f097172
--- /dev/null
+++ b/examples/data_fitting.cc
@@ -0,0 +1,164 @@
+// 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"
+
+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/powell.cc b/examples/powell.cc
new file mode 100644
index 0000000..1c7cc68
--- /dev/null
+++ b/examples/powell.cc
@@ -0,0 +1,150 @@
+// 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)
+//
+// An example program that minimizes Powell's singular function.
+//
+//   F = 1/2 (f1^2 + f2^2 + f3^2 + f4^2)
+//
+//   f1 = x1 + 10*x2;
+//   f2 = sqrt(5) * (x3 - x4)
+//   f3 = (x2 - 2*x3)^2
+//   f4 = sqrt(10) * (x1 - x4)^2
+//
+// The starting values are x1 = 3, x2 = -1, x3 = 0, x4 = 1.
+// The minimum is 0 at (x1, x2, x3, x4) = 0.
+//
+// From: Testing Unconstrained Optimization Software by Jorge J. More, Burton S.
+// Garbow and Kenneth E. Hillstrom in ACM Transactions on Mathematical Software,
+// Vol 7(1), March 1981.
+
+#include <vector>
+
+#include "ceres/ceres.h"
+
+using ceres::AutoDiffCostFunction;
+using ceres::CostFunction;
+using ceres::Problem;
+using ceres::Solver;
+using ceres::Solve;
+
+class F1 {
+ public:
+  template <typename T> bool operator()(const T* const x1,
+                                        const T* const x2,
+                                        T* residual) const {
+    // f1 = x1 + 10 * x2;
+    residual[0] = x1[0] + T(10.0) * x2[0];
+    return true;
+  }
+};
+
+class F2 {
+ public:
+  template <typename T> bool operator()(const T* const x3,
+                                        const T* const x4,
+                                        T* residual) const {
+    // f2 = sqrt(5) (x3 - x4)
+    residual[0] = T(sqrt(5.0)) * (x3[0] - x4[0]);
+    return true;
+  }
+};
+
+class F3 {
+ public:
+  template <typename T> bool operator()(const T* const x2,
+                                        const T* const x4,
+                                        T* residual) const {
+    // f3 = (x2 - 2 x3)^2
+    residual[0] = (x2[0] - T(2.0) * x4[0]) * (x2[0] - T(2.0) * x4[0]);
+    return true;
+  }
+};
+
+class F4 {
+ public:
+  template <typename T> bool operator()(const T* const x1,
+                                        const T* const x4,
+                                        T* residual) const {
+    // f4 = sqrt(10) (x1 - x4)^2
+    residual[0] = T(sqrt(10.0)) * (x1[0] - x4[0]) * (x1[0] - x4[0]);
+    return true;
+  }
+};
+
+int main(int argc, char** argv) {
+  google::ParseCommandLineFlags(&argc, &argv, true);
+  google::InitGoogleLogging(argv[0]);
+
+  double x1 =  3.0;
+  double x2 = -1.0;
+  double x3 =  0.0;
+  double x4 =  1.0;
+
+  Problem problem;
+  // Add residual terms to the problem using the using the autodiff
+  // wrapper to get the derivatives automatically. The parameters, x1 through
+  // x4, are modified in place.
+  problem.AddResidualBlock(new AutoDiffCostFunction<F1, 1, 1, 1>(new F1),
+                           NULL,
+                           &x1, &x2);
+  problem.AddResidualBlock(new AutoDiffCostFunction<F2, 1, 1, 1>(new F2),
+                           NULL,
+                           &x3, &x4);
+  problem.AddResidualBlock(new AutoDiffCostFunction<F3, 1, 1, 1>(new F3),
+                           NULL,
+                           &x2, &x3);
+  problem.AddResidualBlock(new AutoDiffCostFunction<F4, 1, 1, 1>(new F4),
+                           NULL,
+                           &x1, &x4);
+
+  // Run the solver!
+  Solver::Options options;
+  options.max_num_iterations = 30;
+  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";
+
+  Solve(options, &problem, &summary);
+
+  std::cout << summary.BriefReport() << "\n";
+  std::cout << "Final x1 = " << x1
+            << ", x2 = " << x2
+            << ", x3 = " << x3
+            << ", x4 = " << x4
+            << "\n";
+  return 0;
+}
diff --git a/examples/quadratic.cc b/examples/quadratic.cc
new file mode 100644
index 0000000..81ba4f9
--- /dev/null
+++ b/examples/quadratic.cc
@@ -0,0 +1,89 @@
+// 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"
+
+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
new file mode 100644
index 0000000..71b216b
--- /dev/null
+++ b/examples/quadratic_auto_diff.cc
@@ -0,0 +1,87 @@
+// 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"
+
+using ceres::AutoDiffCostFunction;
+using ceres::CostFunction;
+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 {
+ 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<QuadraticCostFunction, 1, 1>(
+          new QuadraticCostFunction),
+      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
new file mode 100644
index 0000000..933dbc7
--- /dev/null
+++ b/examples/quadratic_numeric_diff.cc
@@ -0,0 +1,91 @@
+// 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"
+
+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 */> {
+ 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];
+    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<ResidualWithNoDerivative, CENTRAL, 1, 1> (
+          new ResidualWithNoDerivative, ceres::TAKE_OWNERSHIP);
+
+  // 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/simple_bundle_adjuster.cc b/examples/simple_bundle_adjuster.cc
new file mode 100644
index 0000000..cb1143c
--- /dev/null
+++ b/examples/simple_bundle_adjuster.cc
@@ -0,0 +1,210 @@
+// 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 minimal, self-contained bundle adjuster using Ceres, that reads
+// files from University of Washington' Bundle Adjustment in the Large dataset:
+// http://grail.cs.washington.edu/projects/bal
+//
+// This does not use the best configuration for solving; see the more involved
+// bundle_adjuster.cc file for details.
+
+#include <cmath>
+#include <cstdio>
+#include <iostream>
+
+#include "ceres/ceres.h"
+#include "ceres/rotation.h"
+
+// Read a Bundle Adjustment in the Large dataset.
+class BALProblem {
+ public:
+  ~BALProblem() {
+    delete[] point_index_;
+    delete[] camera_index_;
+    delete[] observations_;
+    delete[] parameters_;
+  }
+
+  int num_observations()       const { return num_observations_;               }
+  const double* observations() const { return observations_;                   }
+  double* mutable_cameras()          { return parameters_;                     }
+  double* mutable_points()           { return parameters_  + 9 * num_cameras_; }
+
+  double* mutable_camera_for_observation(int i) {
+    return mutable_cameras() + camera_index_[i] * 9;
+  }
+  double* mutable_point_for_observation(int i) {
+    return mutable_points() + point_index_[i] * 3;
+  }
+
+  bool LoadFile(const char* filename) {
+    FILE* fptr = fopen(filename, "r");
+    if (fptr == NULL) {
+      return false;
+    };
+
+    FscanfOrDie(fptr, "%d", &num_cameras_);
+    FscanfOrDie(fptr, "%d", &num_points_);
+    FscanfOrDie(fptr, "%d", &num_observations_);
+
+    point_index_ = new int[num_observations_];
+    camera_index_ = new int[num_observations_];
+    observations_ = new double[2 * num_observations_];
+
+    num_parameters_ = 9 * num_cameras_ + 3 * num_points_;
+    parameters_ = new double[num_parameters_];
+
+    for (int i = 0; i < num_observations_; ++i) {
+      FscanfOrDie(fptr, "%d", camera_index_ + i);
+      FscanfOrDie(fptr, "%d", point_index_ + i);
+      for (int j = 0; j < 2; ++j) {
+        FscanfOrDie(fptr, "%lf", observations_ + 2*i + j);
+      }
+    }
+
+    for (int i = 0; i < num_parameters_; ++i) {
+      FscanfOrDie(fptr, "%lf", parameters_ + i);
+    }
+    return true;
+  }
+
+ private:
+  template<typename T>
+  void FscanfOrDie(FILE *fptr, const char *format, T *value) {
+    int num_scanned = fscanf(fptr, format, value);
+    if (num_scanned != 1) {
+      LOG(FATAL) << "Invalid UW data file.";
+    }
+  }
+
+  int num_cameras_;
+  int num_points_;
+  int num_observations_;
+  int num_parameters_;
+
+  int* point_index_;
+  int* camera_index_;
+  double* observations_;
+  double* parameters_;
+};
+
+// Templated pinhole camera model for used with Ceres.  The camera is
+// parameterized using 9 parameters: 3 for rotation, 3 for translation, 1 for
+// focal length and 2 for radial distortion. The principal point is not modeled
+// (i.e. it is assumed be located at the image center).
+struct SnavelyReprojectionError {
+  SnavelyReprojectionError(double observed_x, double observed_y)
+      : observed_x(observed_x), observed_y(observed_y) {}
+
+  template <typename T>
+  bool operator()(const T* const camera,
+                  const T* const point,
+                  T* residuals) const {
+    // camera[0,1,2] are the angle-axis rotation.
+    T p[3];
+    ceres::AngleAxisRotatePoint(camera, point, p);
+
+    // camera[3,4,5] are the translation.
+    p[0] += camera[3];
+    p[1] += camera[4];
+    p[2] += camera[5];
+
+    // Compute the center of distortion. The sign change comes from
+    // the camera model that Noah Snavely's Bundler assumes, whereby
+    // the camera coordinate system has a negative z axis.
+    const T& focal = camera[6];
+    T xp = - focal * p[0] / p[2];
+    T yp = - focal * p[1] / p[2];
+
+    // Apply second and fourth order radial distortion.
+    const T& l1 = camera[7];
+    const T& l2 = camera[8];
+    T r2 = xp*xp + yp*yp;
+    T distortion = T(1.0) + r2  * (l1 + l2  * r2);
+
+    // Compute final projected point position.
+    T predicted_x = distortion * xp;
+    T predicted_y = distortion * yp;
+
+    // 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;
+  }
+
+  double observed_x;
+  double observed_y;
+};
+
+int main(int argc, char** argv) {
+  if (argc != 2) {
+    std::cerr << "usage: simple_bundle_adjuster <bal_problem>\n";
+    return 1;
+  }
+
+  BALProblem bal_problem;
+  if (!bal_problem.LoadFile(argv[1])) {
+    std::cerr << "ERROR: unable to open file " << argv[1] << "\n";
+    return 1;
+  }
+
+  // Create residuals for each observation in the bundle adjustment problem. The
+  // parameters for cameras and points are added automatically.
+  ceres::Problem problem;
+  for (int i = 0; i < bal_problem.num_observations(); ++i) {
+    // 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]));
+
+    problem.AddResidualBlock(cost_function,
+                             NULL /* squared loss */,
+                             bal_problem.mutable_camera_for_observation(i),
+                             bal_problem.mutable_point_for_observation(i));
+  }
+
+  // Make Ceres automatically detect the bundle structure. Note that the
+  // standard solver, SPARSE_NORMAL_CHOLESKY, also works fine but it is slower
+  // for standard bundle adjustment problems.
+  ceres::Solver::Options options;
+  options.linear_solver_type = ceres::DENSE_SCHUR;
+  options.ordering_type = ceres::SCHUR;
+  options.minimizer_progress_to_stdout = true;
+
+  ceres::Solver::Summary summary;
+  ceres::Solve(options, &problem, &summary);
+  std::cout << summary.FullReport() << "\n";
+  return 0;
+}
diff --git a/examples/snavely_reprojection_error.h b/examples/snavely_reprojection_error.h
new file mode 100644
index 0000000..e43538f
--- /dev/null
+++ b/examples/snavely_reprojection_error.h
@@ -0,0 +1,156 @@
+// 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)
+//
+// Templated struct implementing the camera model and residual
+// computation for bundle adjustment used by Noah Snavely's Bundler
+// SfM system. This is also the camera model/residual for the bundle
+// adjustment problems in the BAL dataset. It is templated so that we
+// can use Ceres's automatic differentiation to compute analytic
+// jacobians.
+//
+// For details see: http://phototour.cs.washington.edu/bundler/
+// and http://grail.cs.washington.edu/projects/bal/
+
+#ifndef CERES_EXAMPLES_SNAVELY_REPROJECTION_ERROR_H_
+#define CERES_EXAMPLES_SNAVELY_REPROJECTION_ERROR_H_
+
+#include "ceres/rotation.h"
+
+namespace ceres {
+namespace examples {
+
+// Templated pinhole camera model for used with Ceres.  The camera is
+// parameterized using 9 parameters: 3 for rotation, 3 for translation, 1 for
+// focal length and 2 for radial distortion. The principal point is not modeled
+// (i.e. it is assumed be located at the image center).
+struct SnavelyReprojectionError {
+  SnavelyReprojectionError(double observed_x, double observed_y)
+      : observed_x(observed_x), observed_y(observed_y) {}
+
+  template <typename T>
+  bool operator()(const T* const camera,
+                  const T* const point,
+                  T* residuals) const {
+    // camera[0,1,2] are the angle-axis rotation.
+    T p[3];
+    ceres::AngleAxisRotatePoint(camera, point, p);
+
+    // camera[3,4,5] are the translation.
+    p[0] += camera[3];
+    p[1] += camera[4];
+    p[2] += camera[5];
+
+    // Compute the center of distortion. The sign change comes from
+    // the camera model that Noah Snavely's Bundler assumes, whereby
+    // the camera coordinate system has a negative z axis.
+    const T& focal = camera[6];
+    T xp = - focal * p[0] / p[2];
+    T yp = - focal * p[1] / p[2];
+
+    // Apply second and fourth order radial distortion.
+    const T& l1 = camera[7];
+    const T& l2 = camera[8];
+    T r2 = xp*xp + yp*yp;
+    T distortion = T(1.0) + r2  * (l1 + l2  * r2);
+
+    // Compute final projected point position.
+    T predicted_x = distortion * xp;
+    T predicted_y = distortion * yp;
+
+    // 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;
+  }
+
+  double observed_x;
+  double observed_y;
+};
+
+// Templated pinhole camera model for used with Ceres.  The camera is
+// parameterized using 10 parameters. 4 for rotation, 3 for
+// translation, 1 for focal length and 2 for radial distortion. The
+// principal point is not modeled (i.e. it is assumed be located at
+// the image center).
+struct SnavelyReprojectionErrorWitQuaternions {
+  // (u, v): the position of the observation with respect to the image
+  // center point.
+  SnavelyReprojectionErrorWitQuaternions(double observed_x, double observed_y)
+      : observed_x(observed_x), observed_y(observed_y) {}
+
+  template <typename T>
+  bool operator()(const T* const camera_rotation,
+                  const T* const camera_translation_and_intrinsics,
+                  const T* const point,
+                  T* residuals) const {
+    const T& focal = camera_translation_and_intrinsics[3];
+    const T& l1 = camera_translation_and_intrinsics[4];
+    const T& l2 = camera_translation_and_intrinsics[5];
+
+    // Use a quaternion rotation that doesn't assume the quaternion is
+    // normalized, since one of the ways to run the bundler is to let Ceres
+    // optimize all 4 quaternion parameters unconstrained.
+    T p[3];
+    QuaternionRotatePoint(camera_rotation, point, p);
+
+    p[0] += camera_translation_and_intrinsics[0];
+    p[1] += camera_translation_and_intrinsics[1];
+    p[2] += camera_translation_and_intrinsics[2];
+
+    // Compute the center of distortion. The sign change comes from
+    // the camera model that Noah Snavely's Bundler assumes, whereby
+    // the camera coordinate system has a negative z axis.
+    T xp = - focal * p[0] / p[2];
+    T yp = - focal * p[1] / p[2];
+
+    // Apply second and fourth order radial distortion.
+    T r2 = xp*xp + yp*yp;
+    T distortion = T(1.0) + r2  * (l1 + l2  * r2);
+
+    // Compute final projected point position.
+    T predicted_x = distortion * xp;
+    T predicted_y = distortion * yp;
+
+    // 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;
+  }
+
+  double observed_x;
+  double observed_y;
+};
+
+}  // namespace examples
+}  // namespace ceres
+
+#endif  // CERES_EXAMPLES_SNAVELY_REPROJECTION_ERROR_H_