11d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// Ceres Solver - A fast non-linear least squares minimizer
21d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// Copyright 2013 Google Inc. All rights reserved.
31d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// http://code.google.com/p/ceres-solver/
41d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling//
51d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// Redistribution and use in source and binary forms, with or without
61d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// modification, are permitted provided that the following conditions are met:
71d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling//
81d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// * Redistributions of source code must retain the above copyright notice,
91d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling//   this list of conditions and the following disclaimer.
101d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// * Redistributions in binary form must reproduce the above copyright notice,
111d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling//   this list of conditions and the following disclaimer in the documentation
121d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling//   and/or other materials provided with the distribution.
131d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// * Neither the name of Google Inc. nor the names of its contributors may be
141d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling//   used to endorse or promote products derived from this software without
151d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling//   specific prior written permission.
161d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling//
171d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
181d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
191d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
201d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
211d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
221d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
231d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
241d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
251d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
261d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
271d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// POSSIBILITY OF SUCH DAMAGE.
281d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling//
291d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// Author: mierle@gmail.com (Keir Mierle)
301d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling//
311d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// An incomplete C API for Ceres.
321d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling//
331d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// TODO(keir): Figure out why logging does not seem to work.
341d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling
351d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling#include "ceres/c_api.h"
361d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling
371d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling#include <vector>
381d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling#include <iostream>
391d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling#include <string>
401d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling#include "ceres/cost_function.h"
411d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling#include "ceres/loss_function.h"
421d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling#include "ceres/problem.h"
431d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling#include "ceres/solver.h"
441d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling#include "ceres/types.h"  // for std
451d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling#include "glog/logging.h"
461d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling
471d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberlingusing ceres::Problem;
481d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling
491d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberlingvoid ceres_init() {
501d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling  // This is not ideal, but it's not clear what to do if there is no gflags and
511d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling  // no access to command line arguments.
52399f7d09e0c45af54b77b4ab9508d6f23759b927Scott Ettinger  char message[] = "<unknown>";
53399f7d09e0c45af54b77b4ab9508d6f23759b927Scott Ettinger  google::InitGoogleLogging(message);
541d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling}
551d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling
561d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberlingceres_problem_t* ceres_create_problem() {
571d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling  return reinterpret_cast<ceres_problem_t*>(new Problem);
581d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling}
591d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling
601d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberlingvoid ceres_free_problem(ceres_problem_t* problem) {
611d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling  delete reinterpret_cast<Problem*>(problem);
621d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling}
631d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling
641d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// This cost function wraps a C-level function pointer from the user, to bridge
651d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// between C and C++.
661d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberlingclass CallbackCostFunction : public ceres::CostFunction {
671d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling public:
681d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling  CallbackCostFunction(ceres_cost_function_t cost_function,
691d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling                       void* user_data,
701d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling                       int num_residuals,
711d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling                       int num_parameter_blocks,
721d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling                       int* parameter_block_sizes)
731d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling      : cost_function_(cost_function),
741d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling        user_data_(user_data) {
751d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling    set_num_residuals(num_residuals);
761d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling    for (int i = 0; i < num_parameter_blocks; ++i) {
771d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling      mutable_parameter_block_sizes()->push_back(parameter_block_sizes[i]);
781d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling    }
791d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling  }
801d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling
811d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling  virtual ~CallbackCostFunction() {}
821d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling
831d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling  virtual bool Evaluate(double const* const* parameters,
841d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling                        double* residuals,
851d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling                        double** jacobians) const {
861d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling    return (*cost_function_)(user_data_,
871d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling                             const_cast<double**>(parameters),
881d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling                             residuals,
891d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling                             jacobians);
901d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling  }
911d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling
921d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling private:
931d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling  ceres_cost_function_t cost_function_;
941d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling  void* user_data_;
951d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling};
961d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling
971d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// This loss function wraps a C-level function pointer from the user, to bridge
981d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// between C and C++.
991d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberlingclass CallbackLossFunction : public ceres::LossFunction {
1001d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling public:
1011d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling  explicit CallbackLossFunction(ceres_loss_function_t loss_function,
1021d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling                                void* user_data)
1031d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling    : loss_function_(loss_function), user_data_(user_data) {}
1041d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling  virtual void Evaluate(double sq_norm, double* rho) const {
1051d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling    (*loss_function_)(user_data_, sq_norm, rho);
1061d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling  }
1071d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling
1081d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling private:
1091d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling  ceres_loss_function_t loss_function_;
1101d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling  void* user_data_;
1111d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling};
1121d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling
1131d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// Wrappers for the stock loss functions.
1141d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberlingvoid* ceres_create_huber_loss_function_data(double a) {
1151d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling  return new ceres::HuberLoss(a);
1161d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling}
1171d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberlingvoid* ceres_create_softl1_loss_function_data(double a) {
1181d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling  return new ceres::SoftLOneLoss(a);
1191d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling}
1201d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberlingvoid* ceres_create_cauchy_loss_function_data(double a) {
1211d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling  return new ceres::CauchyLoss(a);
1221d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling}
1231d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberlingvoid* ceres_create_arctan_loss_function_data(double a) {
1241d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling  return new ceres::ArctanLoss(a);
1251d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling}
1261d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberlingvoid* ceres_create_tolerant_loss_function_data(double a, double b) {
1271d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling  return new ceres::TolerantLoss(a, b);
1281d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling}
1291d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling
1301d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberlingvoid ceres_free_stock_loss_function_data(void* loss_function_data) {
1311d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling  delete reinterpret_cast<ceres::LossFunction*>(loss_function_data);
1321d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling}
1331d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling
1341d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberlingvoid ceres_stock_loss_function(void* user_data,
1351d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling                               double squared_norm,
1361d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling                               double out[3]) {
1371d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling  reinterpret_cast<ceres::LossFunction*>(user_data)
1381d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling      ->Evaluate(squared_norm, out);
1391d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling}
1401d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling
1411d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberlingceres_residual_block_id_t* ceres_problem_add_residual_block(
1421d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling    ceres_problem_t* problem,
1431d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling    ceres_cost_function_t cost_function,
1441d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling    void* cost_function_data,
1451d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling    ceres_loss_function_t loss_function,
1461d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling    void* loss_function_data,
1471d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling    int num_residuals,
1481d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling    int num_parameter_blocks,
1491d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling    int* parameter_block_sizes,
1501d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling    double** parameters) {
1511d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling  Problem* ceres_problem = reinterpret_cast<Problem*>(problem);
1521d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling
1531d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling  ceres::CostFunction* callback_cost_function =
1541d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling      new CallbackCostFunction(cost_function,
1551d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling                               cost_function_data,
1561d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling                               num_residuals,
1571d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling                               num_parameter_blocks,
1581d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling                               parameter_block_sizes);
1591d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling
1601d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling  ceres::LossFunction* callback_loss_function = NULL;
1611d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling  if (loss_function != NULL) {
1621d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling    callback_loss_function = new CallbackLossFunction(loss_function,
1631d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling                                                      loss_function_data);
1641d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling  }
1651d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling
1661d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling  std::vector<double*> parameter_blocks(parameters,
1671d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling                                        parameters + num_parameter_blocks);
1681d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling  return reinterpret_cast<ceres_residual_block_id_t*>(
1691d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling      ceres_problem->AddResidualBlock(callback_cost_function,
1701d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling                                      callback_loss_function,
1711d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling                                      parameter_blocks));
1721d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling}
1731d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling
1741d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberlingvoid ceres_solve(ceres_problem_t* c_problem) {
1751d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling  Problem* problem = reinterpret_cast<Problem*>(c_problem);
176399f7d09e0c45af54b77b4ab9508d6f23759b927Scott Ettinger
1771d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling  // TODO(keir): Obviously, this way of setting options won't scale or last.
1781d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling  // Instead, figure out a way to specify some of the options without
1791d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling  // duplicating everything.
1801d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling  ceres::Solver::Options options;
1811d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling  options.max_num_iterations = 100;
1821d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling  options.linear_solver_type = ceres::DENSE_QR;
1831d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling  options.minimizer_progress_to_stdout = true;
1841d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling
1851d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling  ceres::Solver::Summary summary;
1861d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling  ceres::Solve(options, problem, &summary);
1871d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling  std::cout << summary.FullReport() << "\n";
1881d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling}
189