diff --git a/gtsam/linear/NoiseModel.cpp b/gtsam/linear/NoiseModel.cpp index 80210c3b7..e539d9579 100644 --- a/gtsam/linear/NoiseModel.cpp +++ b/gtsam/linear/NoiseModel.cpp @@ -807,7 +807,10 @@ GemanMcClure::GemanMcClure(double c, const ReweightScheme reweight) } double GemanMcClure::weight(double error) const { - return c_/(c_ + error*error); + const double c2 = c_*c_; + const double c4 = c2*c2; + const double c2error = c2 + error*error; + return c4/(c2error*c2error); } void GemanMcClure::print(const std::string &s="") const { diff --git a/gtsam/linear/NoiseModel.h b/gtsam/linear/NoiseModel.h index 25b34ce18..db01152f6 100644 --- a/gtsam/linear/NoiseModel.h +++ b/gtsam/linear/NoiseModel.h @@ -827,7 +827,8 @@ namespace gtsam { /// (Zhang97ivc). /// /// Note that Geman-McClure weight function uses the parameter c == 1.0, - /// but here it's allowed to use different values. + /// but here it's allowed to use different values, so we actually have + /// the generalized Geman-McClure from (Agarwal15phd). class GTSAM_EXPORT GemanMcClure : public Base { public: typedef boost::shared_ptr shared_ptr; diff --git a/gtsam/linear/tests/testNoiseModel.cpp b/gtsam/linear/tests/testNoiseModel.cpp index c1e9d95f6..01f653d73 100644 --- a/gtsam/linear/tests/testNoiseModel.cpp +++ b/gtsam/linear/tests/testNoiseModel.cpp @@ -338,8 +338,8 @@ TEST(NoiseModel, robustFunctionGemanMcClure) const mEstimator::GemanMcClure::shared_ptr gmc = mEstimator::GemanMcClure::Create(k); const double weight1 = gmc->weight(error1), weight2 = gmc->weight(error2); - DOUBLES_EQUAL(0.5 , weight1, 1e-8); - DOUBLES_EQUAL(0.00990099, weight2, 1e-8); + DOUBLES_EQUAL(0.25 , weight1, 1e-8); + DOUBLES_EQUAL(9.80296e-5, weight2, 1e-8); } TEST(NoiseModel, robustFunctionDCS) @@ -385,8 +385,12 @@ TEST(NoiseModel, robustNoiseGemanMcClure) robust->WhitenSystem(A, b); - const double sqrt_weight_error1 = sqrt(0.5); - const double sqrt_weight_error2 = sqrt(k/(k + error2*error2)); + const double k2 = k*k; + const double k4 = k2*k2; + const double k2error = k2 + error2*error2; + + const double sqrt_weight_error1 = sqrt(0.25); + const double sqrt_weight_error2 = sqrt(k4/(k2error*k2error)); DOUBLES_EQUAL(sqrt_weight_error1*error1, b(0), 1e-8); DOUBLES_EQUAL(sqrt_weight_error2*error2, b(1), 1e-8);