1 // Ceres Solver - A fast non-linear least squares minimizer 2 // Copyright 2010, 2011, 2012 Google Inc. All rights reserved. 3 // http://code.google.com/p/ceres-solver/ 4 // 5 // Redistribution and use in source and binary forms, with or without 6 // modification, are permitted provided that the following conditions are met: 7 // 8 // * Redistributions of source code must retain the above copyright notice, 9 // this list of conditions and the following disclaimer. 10 // * Redistributions in binary form must reproduce the above copyright notice, 11 // this list of conditions and the following disclaimer in the documentation 12 // and/or other materials provided with the distribution. 13 // * Neither the name of Google Inc. nor the names of its contributors may be 14 // used to endorse or promote products derived from this software without 15 // specific prior written permission. 16 // 17 // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" 18 // AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE 19 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE 20 // ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE 21 // LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR 22 // CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF 23 // SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS 24 // INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN 25 // CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) 26 // ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE 27 // POSSIBILITY OF SUCH DAMAGE. 28 // 29 // Author: sameeragarwal (at) google.com (Sameer Agarwal) 30 // 31 // Purpose: See .h file. 32 33 #include "ceres/loss_function.h" 34 35 #include <cmath> 36 #include <cstddef> 37 38 namespace ceres { 39 40 void TrivialLoss::Evaluate(double s, double rho[3]) const { 41 rho[0] = s; 42 rho[1] = 1; 43 rho[2] = 0; 44 } 45 46 void HuberLoss::Evaluate(double s, double rho[3]) const { 47 if (s > b_) { 48 // Outlier region. 49 // 'r' is always positive. 50 const double r = sqrt(s); 51 rho[0] = 2 * a_ * r - b_; 52 rho[1] = a_ / r; 53 rho[2] = - rho[1] / (2 * s); 54 } else { 55 // Inlier region. 56 rho[0] = s; 57 rho[1] = 1; 58 rho[2] = 0; 59 } 60 } 61 62 void SoftLOneLoss::Evaluate(double s, double rho[3]) const { 63 const double sum = 1 + s * c_; 64 const double tmp = sqrt(sum); 65 // 'sum' and 'tmp' are always positive, assuming that 's' is. 66 rho[0] = 2 * b_ * (tmp - 1); 67 rho[1] = 1 / tmp; 68 rho[2] = - (c_ * rho[1]) / (2 * sum); 69 } 70 71 void CauchyLoss::Evaluate(double s, double rho[3]) const { 72 const double sum = 1 + s * c_; 73 const double inv = 1 / sum; 74 // 'sum' and 'inv' are always positive, assuming that 's' is. 75 rho[0] = b_ * log(sum); 76 rho[1] = inv; 77 rho[2] = - c_ * (inv * inv); 78 } 79 80 void ArctanLoss::Evaluate(double s, double rho[3]) const { 81 const double sum = 1 + s * s * b_; 82 const double inv = 1 / sum; 83 // 'sum' and 'inv' are always positive. 84 rho[0] = a_ * atan2(s, a_); 85 rho[1] = inv; 86 rho[2] = -2 * s * b_ * (inv * inv); 87 } 88 89 TolerantLoss::TolerantLoss(double a, double b) 90 : a_(a), 91 b_(b), 92 c_(b * log(1.0 + exp(-a / b))) { 93 CHECK_GE(a, 0.0); 94 CHECK_GT(b, 0.0); 95 } 96 97 void TolerantLoss::Evaluate(double s, double rho[3]) const { 98 const double x = (s - a_) / b_; 99 // The basic equation is rho[0] = b ln(1 + e^x). However, if e^x is too 100 // large, it will overflow. Since numerically 1 + e^x == e^x when the 101 // x is greater than about ln(2^53) for doubles, beyond this threshold 102 // we substitute x for ln(1 + e^x) as a numerically equivalent approximation. 103 static const double kLog2Pow53 = 36.7; // ln(MathLimits<double>::kEpsilon). 104 if (x > kLog2Pow53) { 105 rho[0] = s - a_ - c_; 106 rho[1] = 1.0; 107 rho[2] = 0.0; 108 } else { 109 const double e_x = exp(x); 110 rho[0] = b_ * log(1.0 + e_x) - c_; 111 rho[1] = e_x / (1.0 + e_x); 112 rho[2] = 0.5 / (b_ * (1.0 + cosh(x))); 113 } 114 } 115 116 ComposedLoss::ComposedLoss(const LossFunction* f, Ownership ownership_f, 117 const LossFunction* g, Ownership ownership_g) 118 : f_(CHECK_NOTNULL(f)), 119 g_(CHECK_NOTNULL(g)), 120 ownership_f_(ownership_f), 121 ownership_g_(ownership_g) { 122 } 123 124 ComposedLoss::~ComposedLoss() { 125 if (ownership_f_ == DO_NOT_TAKE_OWNERSHIP) { 126 f_.release(); 127 } 128 if (ownership_g_ == DO_NOT_TAKE_OWNERSHIP) { 129 g_.release(); 130 } 131 } 132 133 void ComposedLoss::Evaluate(double s, double rho[3]) const { 134 double rho_f[3], rho_g[3]; 135 g_->Evaluate(s, rho_g); 136 f_->Evaluate(rho_g[0], rho_f); 137 rho[0] = rho_f[0]; 138 // f'(g(s)) * g'(s). 139 rho[1] = rho_f[1] * rho_g[1]; 140 // f''(g(s)) * g'(s) * g'(s) + f'(g(s)) * g''(s). 141 rho[2] = rho_f[2] * rho_g[1] * rho_g[1] + rho_f[1] * rho_g[2]; 142 } 143 144 void ScaledLoss::Evaluate(double s, double rho[3]) const { 145 if (rho_.get() == NULL) { 146 rho[0] = a_ * s; 147 rho[1] = a_; 148 rho[2] = 0.0; 149 } else { 150 rho_->Evaluate(s, rho); 151 rho[0] *= a_; 152 rho[1] *= a_; 153 rho[2] *= a_; 154 } 155 } 156 157 } // namespace ceres 158
