diff --git a/.cproject b/.cproject
index 62b32938c..412359938 100644
--- a/.cproject
+++ b/.cproject
@@ -2353,6 +2353,38 @@
true
true
+
+ make
+ -j5
+ testSpirit.run
+ true
+ true
+ true
+
+
+ make
+ -j5
+ check.wrap
+ true
+ true
+ true
+
+
+ make
+ -j5
+ testMethod.run
+ true
+ true
+ true
+
+
+ make
+ -j5
+ testClass.run
+ true
+ true
+ true
+
make
-j5
@@ -3222,22 +3254,6 @@
true
true
-
- make
- -j5
- testSpirit.run
- true
- true
- true
-
-
- make
- -j5
- check.wrap
- true
- true
- true
-
make
-j5
diff --git a/CMakeLists.txt b/CMakeLists.txt
index bc7d8aecd..754279e53 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -9,8 +9,8 @@ if(NOT DEFINED CMAKE_MACOSX_RPATH)
endif()
# Set the version number for the library
-set (GTSAM_VERSION_MAJOR 3)
-set (GTSAM_VERSION_MINOR 1)
+set (GTSAM_VERSION_MAJOR 4)
+set (GTSAM_VERSION_MINOR 0)
set (GTSAM_VERSION_PATCH 0)
math (EXPR GTSAM_VERSION_NUMERIC "10000 * ${GTSAM_VERSION_MAJOR} + 100 * ${GTSAM_VERSION_MINOR} + ${GTSAM_VERSION_PATCH}")
set (GTSAM_VERSION_STRING "${GTSAM_VERSION_MAJOR}.${GTSAM_VERSION_MINOR}.${GTSAM_VERSION_PATCH}")
diff --git a/examples/Data/pose3example-rewritten.txt b/examples/Data/pose3example-rewritten.txt
deleted file mode 100644
index 374e6c171..000000000
--- a/examples/Data/pose3example-rewritten.txt
+++ /dev/null
@@ -1,11 +0,0 @@
-VERTEX_SE3:QUAT 0 0 0 0 0 0 0 1
-VERTEX_SE3:QUAT 1 1.00137 0.01539 0.004948 0.190253 0.283162 -0.392318 0.85423
-VERTEX_SE3:QUAT 2 1.9935 0.023275 0.003793 -0.351729 -0.597838 0.584174 0.421446
-VERTEX_SE3:QUAT 3 2.00429 1.02431 0.018047 0.331798 -0.200659 0.919323 0.067024
-VERTEX_SE3:QUAT 4 0.999908 1.05507 0.020212 -0.035697 -0.46249 0.445933 0.765488
-EDGE_SE3:QUAT 0 1 1.00137 0.01539 0.004948 0.190253 0.283162 -0.392318 0.85423 10000 0 0 0 0 0 10000 0 0 0 0 10000 0 0 0 10000 0 0 10000 0 10000
-EDGE_SE3:QUAT 1 2 0.523923 0.776654 0.326659 0.311512 0.656877 -0.678505 0.105373 10000 0 0 0 0 0 10000 0 0 0 0 10000 0 0 0 10000 0 0 10000 0 10000
-EDGE_SE3:QUAT 2 3 0.910927 0.055169 -0.411761 0.595795 -0.561677 0.079353 0.568551 10000 0 0 0 0 0 10000 0 0 0 0 10000 0 0 0 10000 0 0 10000 0 10000
-EDGE_SE3:QUAT 3 4 0.775288 0.228798 -0.596923 -0.592077 0.30338 -0.513226 0.542221 10000 0 0 0 0 0 10000 0 0 0 0 10000 0 0 0 10000 0 0 10000 0 10000
-EDGE_SE3:QUAT 1 4 -0.577841 0.628016 -0.543592 -0.12525 -0.534379 0.769122 0.327419 10000 0 0 0 0 0 10000 0 0 0 0 10000 0 0 0 10000 0 0 10000 0 10000
-EDGE_SE3:QUAT 3 0 -0.623267 0.086928 0.773222 0.104639 0.627755 0.766795 0.083672 10000 0 0 0 0 0 10000 0 0 0 0 10000 0 0 0 10000 0 0 10000 0 10000
diff --git a/gtsam.h b/gtsam.h
index 678e2a3d6..40d204a5f 100644
--- a/gtsam.h
+++ b/gtsam.h
@@ -156,8 +156,14 @@ virtual class Value {
size_t dim() const;
};
+class Vector3 {
+ Vector3(Vector v);
+};
+class Vector6 {
+ Vector6(Vector v);
+};
#include
-virtual class LieScalar : gtsam::Value {
+class LieScalar {
// Standard constructors
LieScalar();
LieScalar(double d);
@@ -186,7 +192,7 @@ virtual class LieScalar : gtsam::Value {
};
#include
-virtual class LieVector : gtsam::Value {
+class LieVector {
// Standard constructors
LieVector();
LieVector(Vector v);
@@ -218,7 +224,7 @@ virtual class LieVector : gtsam::Value {
};
#include
-virtual class LieMatrix : gtsam::Value {
+class LieMatrix {
// Standard constructors
LieMatrix();
LieMatrix(Matrix v);
@@ -253,7 +259,7 @@ virtual class LieMatrix : gtsam::Value {
// geometry
//*************************************************************************
-virtual class Point2 : gtsam::Value {
+class Point2 {
// Standard Constructors
Point2();
Point2(double x, double y);
@@ -290,7 +296,7 @@ virtual class Point2 : gtsam::Value {
void serialize() const;
};
-virtual class StereoPoint2 : gtsam::Value {
+class StereoPoint2 {
// Standard Constructors
StereoPoint2();
StereoPoint2(double uL, double uR, double v);
@@ -325,7 +331,7 @@ virtual class StereoPoint2 : gtsam::Value {
void serialize() const;
};
-virtual class Point3 : gtsam::Value {
+class Point3 {
// Standard Constructors
Point3();
Point3(double x, double y, double z);
@@ -361,7 +367,7 @@ virtual class Point3 : gtsam::Value {
void serialize() const;
};
-virtual class Rot2 : gtsam::Value {
+class Rot2 {
// Standard Constructors and Named Constructors
Rot2();
Rot2(double theta);
@@ -406,7 +412,7 @@ virtual class Rot2 : gtsam::Value {
void serialize() const;
};
-virtual class Rot3 : gtsam::Value {
+class Rot3 {
// Standard Constructors and Named Constructors
Rot3();
Rot3(Matrix R);
@@ -462,7 +468,7 @@ virtual class Rot3 : gtsam::Value {
void serialize() const;
};
-virtual class Pose2 : gtsam::Value {
+class Pose2 {
// Standard Constructor
Pose2();
Pose2(const gtsam::Pose2& pose);
@@ -512,7 +518,7 @@ virtual class Pose2 : gtsam::Value {
void serialize() const;
};
-virtual class Pose3 : gtsam::Value {
+class Pose3 {
// Standard Constructors
Pose3();
Pose3(const gtsam::Pose3& pose);
@@ -564,7 +570,7 @@ virtual class Pose3 : gtsam::Value {
};
#include
-virtual class Unit3 : gtsam::Value {
+class Unit3 {
// Standard Constructors
Unit3();
Unit3(const gtsam::Point3& pose);
@@ -585,7 +591,7 @@ virtual class Unit3 : gtsam::Value {
};
#include
-virtual class EssentialMatrix : gtsam::Value {
+class EssentialMatrix {
// Standard Constructors
EssentialMatrix(const gtsam::Rot3& aRb, const gtsam::Unit3& aTb);
@@ -606,7 +612,7 @@ virtual class EssentialMatrix : gtsam::Value {
double error(Vector vA, Vector vB);
};
-virtual class Cal3_S2 : gtsam::Value {
+class Cal3_S2 {
// Standard Constructors
Cal3_S2();
Cal3_S2(double fx, double fy, double s, double u0, double v0);
@@ -643,7 +649,7 @@ virtual class Cal3_S2 : gtsam::Value {
};
#include
-virtual class Cal3DS2 : gtsam::Value {
+class Cal3DS2 {
// Standard Constructors
Cal3DS2();
Cal3DS2(double fx, double fy, double s, double u0, double v0, double k1, double k2, double k3, double k4);
@@ -699,7 +705,43 @@ class Cal3_S2Stereo {
double baseline() const;
};
-virtual class CalibratedCamera : gtsam::Value {
+#include
+class Cal3Bundler {
+ // Standard Constructors
+ Cal3Bundler();
+ Cal3Bundler(double fx, double k1, double k2, double u0, double v0);
+
+ // Testable
+ void print(string s) const;
+ bool equals(const gtsam::Cal3Bundler& rhs, double tol) const;
+
+ // Manifold
+ static size_t Dim();
+ size_t dim() const;
+ gtsam::Cal3Bundler retract(Vector v) const;
+ Vector localCoordinates(const gtsam::Cal3Bundler& c) const;
+
+ // Action on Point2
+ gtsam::Point2 calibrate(const gtsam::Point2& p, double tol) const;
+ gtsam::Point2 calibrate(const gtsam::Point2& p) const;
+ gtsam::Point2 uncalibrate(const gtsam::Point2& p) const;
+
+ // Standard Interface
+ double fx() const;
+ double fy() const;
+ double k1() const;
+ double k2() const;
+ double u0() const;
+ double v0() const;
+ Vector vector() const;
+ Vector k() const;
+ //Matrix K() const; //FIXME: Uppercase
+
+ // enabling serialization functionality
+ void serialize() const;
+};
+
+class CalibratedCamera {
// Standard Constructors and Named Constructors
CalibratedCamera();
CalibratedCamera(const gtsam::Pose3& pose);
@@ -732,7 +774,7 @@ virtual class CalibratedCamera : gtsam::Value {
void serialize() const;
};
-virtual class SimpleCamera : gtsam::Value {
+class SimpleCamera {
// Standard Constructors and Named Constructors
SimpleCamera();
SimpleCamera(const gtsam::Pose3& pose);
@@ -771,7 +813,7 @@ virtual class SimpleCamera : gtsam::Value {
};
template
-virtual class PinholeCamera : gtsam::Value {
+class PinholeCamera {
// Standard Constructors and Named Constructors
PinholeCamera();
PinholeCamera(const gtsam::Pose3& pose);
@@ -809,7 +851,7 @@ virtual class PinholeCamera : gtsam::Value {
void serialize() const;
};
-virtual class StereoCamera : gtsam::Value {
+class StereoCamera {
// Standard Constructors and Named Constructors
StereoCamera();
StereoCamera(const gtsam::Pose3& pose, const gtsam::Cal3_S2Stereo* K);
@@ -862,7 +904,7 @@ virtual class SymbolicFactor {
};
#include
-class SymbolicFactorGraph {
+virtual class SymbolicFactorGraph {
SymbolicFactorGraph();
SymbolicFactorGraph(const gtsam::SymbolicBayesNet& bayesNet);
SymbolicFactorGraph(const gtsam::SymbolicBayesTree& bayesTree);
@@ -1664,15 +1706,12 @@ class Values {
void print(string s) const;
bool equals(const gtsam::Values& other, double tol) const;
- void insert(size_t j, const gtsam::Value& value);
void insert(const gtsam::Values& values);
- void update(size_t j, const gtsam::Value& val);
void update(const gtsam::Values& values);
void erase(size_t j);
void swap(gtsam::Values& values);
bool exists(size_t j) const;
- gtsam::Value at(size_t j) const;
gtsam::KeyList keys() const;
gtsam::VectorValues zeroVectors() const;
@@ -1682,6 +1721,37 @@ class Values {
// enabling serialization functionality
void serialize() const;
+
+ // New in 4.0, we have to specialize every insert/update/at to generate wrappers
+ // Instead of the old:
+ // void insert(size_t j, const gtsam::Value& value);
+ // void update(size_t j, const gtsam::Value& val);
+ // gtsam::Value at(size_t j) const;
+ void insert(size_t j, const gtsam::Point2& t);
+ void insert(size_t j, const gtsam::Point3& t);
+ void insert(size_t j, const gtsam::Rot2& t);
+ void insert(size_t j, const gtsam::Pose2& t);
+ void insert(size_t j, const gtsam::Rot3& t);
+ void insert(size_t j, const gtsam::Pose3& t);
+ void insert(size_t j, const gtsam::Cal3_S2& t);
+ void insert(size_t j, const gtsam::Cal3DS2& t);
+ void insert(size_t j, const gtsam::Cal3Bundler& t);
+ void insert(size_t j, const gtsam::EssentialMatrix& t);
+
+ void update(size_t j, const gtsam::Point2& t);
+ void update(size_t j, const gtsam::Point3& t);
+ void update(size_t j, const gtsam::Rot2& t);
+ void update(size_t j, const gtsam::Pose2& t);
+ void update(size_t j, const gtsam::Rot3& t);
+ void update(size_t j, const gtsam::Pose3& t);
+ void update(size_t j, const gtsam::Cal3_S2& t);
+ void update(size_t j, const gtsam::Cal3DS2& t);
+ void update(size_t j, const gtsam::Cal3Bundler& t);
+ void update(size_t j, const gtsam::EssentialMatrix& t);
+
+ template
+ T at(size_t j);
};
// Actually a FastList
@@ -2077,7 +2147,7 @@ class NonlinearISAM {
#include
#include
-template
+template
virtual class PriorFactor : gtsam::NoiseModelFactor {
PriorFactor(size_t key, const T& prior, const gtsam::noiseModel::Base* noiseModel);
T prior() const;
@@ -2088,7 +2158,7 @@ virtual class PriorFactor : gtsam::NoiseModelFactor {
#include
-template
+template
virtual class BetweenFactor : gtsam::NoiseModelFactor {
BetweenFactor(size_t key1, size_t key2, const T& relativePose, const gtsam::noiseModel::Base* noiseModel);
T measured() const;
@@ -2099,7 +2169,7 @@ virtual class BetweenFactor : gtsam::NoiseModelFactor {
#include
-template
+template
virtual class NonlinearEquality : gtsam::NoiseModelFactor {
// Constructor - forces exact evaluation
NonlinearEquality(size_t j, const T& feasible);
@@ -2280,7 +2350,7 @@ void writeG2o(const gtsam::NonlinearFactorGraph& graph,
namespace imuBias {
#include
-virtual class ConstantBias : gtsam::Value {
+class ConstantBias {
// Standard Constructor
ConstantBias();
ConstantBias(Vector biasAcc, Vector biasGyro);
@@ -2340,7 +2410,7 @@ virtual class ImuFactor : gtsam::NonlinearFactor {
// Standard Interface
gtsam::ImuFactorPreintegratedMeasurements preintegratedMeasurements() const;
- void Predict(const gtsam::Pose3& pose_i, const gtsam::LieVector& vel_i, gtsam::Pose3& pose_j, gtsam::LieVector& vel_j,
+ void Predict(const gtsam::Pose3& pose_i, const gtsam::Vector3& vel_i, gtsam::Pose3& pose_j, gtsam::Vector3& vel_j,
const gtsam::imuBias::ConstantBias& bias,
const gtsam::ImuFactorPreintegratedMeasurements& preintegratedMeasurements,
Vector gravity, Vector omegaCoriolis) const;
@@ -2383,7 +2453,7 @@ virtual class CombinedImuFactor : gtsam::NonlinearFactor {
// Standard Interface
gtsam::CombinedImuFactorPreintegratedMeasurements preintegratedMeasurements() const;
- void Predict(const gtsam::Pose3& pose_i, const gtsam::LieVector& vel_i, gtsam::Pose3& pose_j, gtsam::LieVector& vel_j,
+ void Predict(const gtsam::Pose3& pose_i, const gtsam::Vector3& vel_i, gtsam::Pose3& pose_j, gtsam::Vector3& vel_j,
const gtsam::imuBias::ConstantBias& bias_i, const gtsam::imuBias::ConstantBias& bias_j,
const gtsam::CombinedImuFactorPreintegratedMeasurements& preintegratedMeasurements,
Vector gravity, Vector omegaCoriolis) const;
diff --git a/gtsam/3rdparty/gtsam_eigen_includes.h.in b/gtsam/3rdparty/gtsam_eigen_includes.h.in
index fa5a51c70..f53e37f07 100644
--- a/gtsam/3rdparty/gtsam_eigen_includes.h.in
+++ b/gtsam/3rdparty/gtsam_eigen_includes.h.in
@@ -27,3 +27,7 @@
#include <@GTSAM_EIGEN_INCLUDE_PREFIX@Eigen/LU>
#include <@GTSAM_EIGEN_INCLUDE_PREFIX@Eigen/SVD>
#include <@GTSAM_EIGEN_INCLUDE_PREFIX@Eigen/Geometry>
+
+
+
+
diff --git a/gtsam/base/GenericValue.h b/gtsam/base/GenericValue.h
index 1c01f7bb7..6c109675d 100644
--- a/gtsam/base/GenericValue.h
+++ b/gtsam/base/GenericValue.h
@@ -70,9 +70,9 @@ struct print {
};
// equals for Matrix types
-template
-struct equals > {
- typedef Eigen::Matrix type;
+template
+struct equals > {
+ typedef Eigen::Matrix type;
typedef bool result_type;
bool operator()(const type& A, const type& B, double tol) {
return equal_with_abs_tol(A, B, tol);
@@ -80,9 +80,9 @@ struct equals > {
};
// print for Matrix types
-template
-struct print > {
- typedef Eigen::Matrix type;
+template
+struct print > {
+ typedef Eigen::Matrix type;
typedef void result_type;
void operator()(const type& A, const std::string& str) {
std::cout << str << ": " << A << std::endl;
diff --git a/gtsam/base/Manifold.h b/gtsam/base/Manifold.h
index ba597ccc6..1d7eda12f 100644
--- a/gtsam/base/Manifold.h
+++ b/gtsam/base/Manifold.h
@@ -253,10 +253,10 @@ struct DefaultChart > {
typedef Eigen::Matrix::value, 1> vector;BOOST_STATIC_ASSERT_MSG((M!=Eigen::Dynamic && N!=Eigen::Dynamic),
"DefaultChart has not been implemented yet for dynamically sized matrices");
static vector local(const T& origin, const T& other) {
- return reshape(other) - reshape(origin);
+ return reshape(other) - reshape(origin);
}
static T retract(const T& origin, const vector& d) {
- return origin + reshape(d);
+ return origin + reshape(d);
}
static int getDimension(const T&origin) {
return origin.rows() * origin.cols();
diff --git a/gtsam/base/Matrix.h b/gtsam/base/Matrix.h
index 00caed44a..132bf79ad 100644
--- a/gtsam/base/Matrix.h
+++ b/gtsam/base/Matrix.h
@@ -294,10 +294,10 @@ void zeroBelowDiagonal(MATRIX& A, size_t cols=0) {
inline Matrix trans(const Matrix& A) { return A.transpose(); }
/// Reshape functor
-template
+template
struct Reshape {
//TODO replace this with Eigen's reshape function as soon as available. (There is a PR already pending : https://bitbucket.org/eigen/eigen/pull-request/41/reshape/diff)
- typedef Eigen::Map > ReshapedType;
+ typedef Eigen::Map > ReshapedType;
static inline ReshapedType reshape(const Eigen::Matrix & in) {
return in.data();
}
@@ -305,7 +305,7 @@ struct Reshape {
/// Reshape specialization that does nothing as shape stays the same (needed to not be ambiguous for square input equals square output)
template
-struct Reshape {
+struct Reshape {
typedef const Eigen::Matrix & ReshapedType;
static inline ReshapedType reshape(const Eigen::Matrix & in) {
return in;
@@ -314,7 +314,7 @@ struct Reshape {
/// Reshape specialization that does nothing as shape stays the same
template
-struct Reshape {
+struct Reshape {
typedef const Eigen::Matrix & ReshapedType;
static inline ReshapedType reshape(const Eigen::Matrix & in) {
return in;
@@ -323,17 +323,17 @@ struct Reshape {
/// Reshape specialization that does transpose
template
-struct Reshape {
+struct Reshape {
typedef typename Eigen::Matrix::ConstTransposeReturnType ReshapedType;
static inline ReshapedType reshape(const Eigen::Matrix & in) {
return in.transpose();
}
};
-template
-inline typename Reshape::ReshapedType reshape(const Eigen::Matrix & m){
+template
+inline typename Reshape::ReshapedType reshape(const Eigen::Matrix & m){
BOOST_STATIC_ASSERT(InM * InN == OutM * OutN);
- return Reshape::reshape(m);
+ return Reshape::reshape(m);
}
/**
diff --git a/gtsam/geometry/Cal3DS2.h b/gtsam/geometry/Cal3DS2.h
index 39f4b7996..95524e115 100644
--- a/gtsam/geometry/Cal3DS2.h
+++ b/gtsam/geometry/Cal3DS2.h
@@ -37,7 +37,7 @@ namespace gtsam {
* k3 (rr + 2 Pn.y^2) + 2*k4 pn.x pn.y ]
* pi = K*pn
*/
-class GTSAM_EXPORT Cal3DS2 : public Cal3DS2_Base, public DerivedValue {
+class GTSAM_EXPORT Cal3DS2 : public Cal3DS2_Base {
typedef Cal3DS2_Base Base;
@@ -87,21 +87,24 @@ public:
/// Return dimensions of calibration manifold object
static size_t Dim() { return 9; } //TODO: make a final dimension variable
- /// @}
private:
+ /// @}
+ /// @name Advanced Interface
+ /// @{
+
/** Serialization function */
friend class boost::serialization::access;
template
void serialize(Archive & ar, const unsigned int version)
{
- ar & boost::serialization::make_nvp("Cal3DS2",
- boost::serialization::base_object(*this));
ar & boost::serialization::make_nvp("Cal3DS2",
boost::serialization::base_object(*this));
}
+ /// @}
+
};
// Define GTSAM traits
diff --git a/gtsam/geometry/Cal3DS2_Base.cpp b/gtsam/geometry/Cal3DS2_Base.cpp
index b8181ab4d..e4397a449 100644
--- a/gtsam/geometry/Cal3DS2_Base.cpp
+++ b/gtsam/geometry/Cal3DS2_Base.cpp
@@ -53,23 +53,23 @@ bool Cal3DS2_Base::equals(const Cal3DS2_Base& K, double tol) const {
}
/* ************************************************************************* */
-static Eigen::Matrix D2dcalibration(double x, double y, double xx,
+static Matrix29 D2dcalibration(double x, double y, double xx,
double yy, double xy, double rr, double r4, double pnx, double pny,
- const Eigen::Matrix& DK) {
- Eigen::Matrix DR1;
+ const Matrix2& DK) {
+ Matrix25 DR1;
DR1 << pnx, 0.0, pny, 1.0, 0.0, 0.0, pny, 0.0, 0.0, 1.0;
- Eigen::Matrix DR2;
+ Matrix24 DR2;
DR2 << x * rr, x * r4, 2 * xy, rr + 2 * xx, //
y * rr, y * r4, rr + 2 * yy, 2 * xy;
- Eigen::Matrix D;
+ Matrix29 D;
D << DR1, DK * DR2;
return D;
}
/* ************************************************************************* */
-static Eigen::Matrix D2dintrinsic(double x, double y, double rr,
+static Matrix2 D2dintrinsic(double x, double y, double rr,
double g, double k1, double k2, double p1, double p2,
- const Eigen::Matrix& DK) {
+ const Matrix2& DK) {
const double drdx = 2. * x;
const double drdy = 2. * y;
const double dgdx = k1 * drdx + k2 * 2. * rr * drdx;
@@ -82,7 +82,7 @@ static Eigen::Matrix D2dintrinsic(double x, double y, double rr,
const double dDydx = 2. * p2 * y + p1 * drdx;
const double dDydy = 2. * p2 * x + p1 * (drdy + 4. * y);
- Eigen::Matrix DR;
+ Matrix2 DR;
DR << g + x * dgdx + dDxdx, x * dgdy + dDxdy, //
y * dgdx + dDydx, g + y * dgdy + dDydy;
@@ -90,8 +90,9 @@ static Eigen::Matrix D2dintrinsic(double x, double y, double rr,
}
/* ************************************************************************* */
-Point2 Cal3DS2_Base::uncalibrate(const Point2& p, boost::optional H1,
- boost::optional H2) const {
+Point2 Cal3DS2_Base::uncalibrate(const Point2& p,
+ boost::optional H1,
+ boost::optional H2) const {
// rr = x^2 + y^2;
// g = (1 + k(1)*rr + k(2)*rr^2);
@@ -110,7 +111,7 @@ Point2 Cal3DS2_Base::uncalibrate(const Point2& p, boost::optional H1,
const double pnx = g * x + dx;
const double pny = g * y + dy;
- Eigen::Matrix DK;
+ Matrix2 DK;
if (H1 || H2) DK << fx_, s_, 0.0, fy_;
// Derivative for calibration
@@ -125,6 +126,26 @@ Point2 Cal3DS2_Base::uncalibrate(const Point2& p, boost::optional H1,
return Point2(fx_ * pnx + s_ * pny + u0_, fy_ * pny + v0_);
}
+/* ************************************************************************* */
+Point2 Cal3DS2_Base::uncalibrate(const Point2& p,
+ boost::optional H1,
+ boost::optional H2) const {
+
+ if (H1 || H2) {
+ Matrix29 D1;
+ Matrix2 D2;
+ Point2 pu = uncalibrate(p, D1, D2);
+ if (H1)
+ *H1 = D1;
+ if (H2)
+ *H2 = D2;
+ return pu;
+
+ } else {
+ return uncalibrate(p);
+ }
+}
+
/* ************************************************************************* */
Point2 Cal3DS2_Base::calibrate(const Point2& pi, const double tol) const {
// Use the following fixed point iteration to invert the radial distortion.
@@ -161,7 +182,7 @@ Matrix Cal3DS2_Base::D2d_intrinsic(const Point2& p) const {
const double rr = xx + yy;
const double r4 = rr * rr;
const double g = (1 + k1_ * rr + k2_ * r4);
- Eigen::Matrix DK;
+ Matrix2 DK;
DK << fx_, s_, 0.0, fy_;
return D2dintrinsic(x, y, rr, g, k1_, k2_, p1_, p2_, DK);
}
@@ -176,7 +197,7 @@ Matrix Cal3DS2_Base::D2d_calibration(const Point2& p) const {
const double dy = 2 * p2_ * xy + p1_ * (rr + 2 * yy);
const double pnx = g * x + dx;
const double pny = g * y + dy;
- Eigen::Matrix DK;
+ Matrix2 DK;
DK << fx_, s_, 0.0, fy_;
return D2dcalibration(x, y, xx, yy, xy, rr, r4, pnx, pny, DK);
}
diff --git a/gtsam/geometry/Cal3DS2_Base.h b/gtsam/geometry/Cal3DS2_Base.h
index 7be5a6874..61c01e349 100644
--- a/gtsam/geometry/Cal3DS2_Base.h
+++ b/gtsam/geometry/Cal3DS2_Base.h
@@ -18,7 +18,6 @@
#pragma once
-#include
#include
namespace gtsam {
@@ -114,9 +113,14 @@ public:
* @param Dp optional 2*2 Jacobian wrpt intrinsic coordinates
* @return point in (distorted) image coordinates
*/
+
Point2 uncalibrate(const Point2& p,
- boost::optional Dcal = boost::none,
- boost::optional Dp = boost::none) const ;
+ boost::optional Dcal = boost::none,
+ boost::optional Dp = boost::none) const ;
+
+ Point2 uncalibrate(const Point2& p,
+ boost::optional Dcal,
+ boost::optional Dp) const ;
/// Convert (distorted) image coordinates uv to intrinsic coordinates xy
Point2 calibrate(const Point2& p, const double tol=1e-5) const;
@@ -127,7 +131,6 @@ public:
/// Derivative of uncalibrate wrpt the calibration parameters
Matrix D2d_calibration(const Point2& p) const ;
-
private:
/// @}
diff --git a/gtsam/geometry/Cal3Unified.cpp b/gtsam/geometry/Cal3Unified.cpp
index e7b408982..6113714a1 100644
--- a/gtsam/geometry/Cal3Unified.cpp
+++ b/gtsam/geometry/Cal3Unified.cpp
@@ -50,6 +50,7 @@ bool Cal3Unified::equals(const Cal3Unified& K, double tol) const {
}
/* ************************************************************************* */
+// todo: make a fixed sized jacobian version of this
Point2 Cal3Unified::uncalibrate(const Point2& p,
boost::optional H1,
boost::optional H2) const {
@@ -70,26 +71,30 @@ Point2 Cal3Unified::uncalibrate(const Point2& p,
// Part2: project NPlane point to pixel plane: use Cal3DS2
Point2 m(x,y);
- Matrix H1base, H2base; // jacobians from Base class
+ Matrix29 H1base;
+ Matrix2 H2base; // jacobians from Base class
Point2 puncalib = Base::uncalibrate(m, H1base, H2base);
// Inlined derivative for calibration
if (H1) {
// part1
- Matrix DU = (Matrix(2,1) << -xs * sqrt_nx * xi_sqrt_nx2,
- -ys * sqrt_nx * xi_sqrt_nx2);
- Matrix DDS2U = H2base * DU;
+ Vector2 DU;
+ DU << -xs * sqrt_nx * xi_sqrt_nx2, //
+ -ys * sqrt_nx * xi_sqrt_nx2;
- *H1 = collect(2, &H1base, &DDS2U);
+ H1->resize(2,10);
+ H1->block<2,9>(0,0) = H1base;
+ H1->block<2,1>(0,9) = H2base * DU;
}
+
// Inlined derivative for points
if (H2) {
// part1
const double denom = 1.0 * xi_sqrt_nx2 / sqrt_nx;
const double mid = -(xi * xs*ys) * denom;
- Matrix DU = (Matrix(2, 2) <<
- (sqrt_nx + xi*(ys*ys + 1)) * denom, mid,
- mid, (sqrt_nx + xi*(xs*xs + 1)) * denom);
+ Matrix2 DU;
+ DU << (sqrt_nx + xi*(ys*ys + 1)) * denom, mid, //
+ mid, (sqrt_nx + xi*(xs*xs + 1)) * denom;
*H2 = H2base * DU;
}
diff --git a/gtsam/geometry/Cal3Unified.h b/gtsam/geometry/Cal3Unified.h
index dc167173e..fb99592f7 100644
--- a/gtsam/geometry/Cal3Unified.h
+++ b/gtsam/geometry/Cal3Unified.h
@@ -40,7 +40,7 @@ namespace gtsam {
* k3 (rr + 2 Pn.y^2) + 2*k4 pn.x pn.y ]
* pi = K*pn
*/
-class GTSAM_EXPORT Cal3Unified : public Cal3DS2_Base, public DerivedValue {
+class GTSAM_EXPORT Cal3Unified : public Cal3DS2_Base {
typedef Cal3Unified This;
typedef Cal3DS2_Base Base;
@@ -129,8 +129,6 @@ private:
template
void serialize(Archive & ar, const unsigned int version)
{
- ar & boost::serialization::make_nvp("Cal3Unified",
- boost::serialization::base_object(*this));
ar & boost::serialization::make_nvp("Cal3Unified",
boost::serialization::base_object(*this));
ar & BOOST_SERIALIZATION_NVP(xi_);
diff --git a/gtsam/geometry/EssentialMatrix.h b/gtsam/geometry/EssentialMatrix.h
index e23b22093..8763504c0 100644
--- a/gtsam/geometry/EssentialMatrix.h
+++ b/gtsam/geometry/EssentialMatrix.h
@@ -58,6 +58,8 @@ public:
return EssentialMatrix(Rot3::Random(rng), Unit3::Random(rng));
}
+ virtual ~EssentialMatrix() {}
+
/// @}
/// @name Testable
diff --git a/gtsam/linear/JacobianFactor.h b/gtsam/linear/JacobianFactor.h
index d1dc7625c..92f43b6be 100644
--- a/gtsam/linear/JacobianFactor.h
+++ b/gtsam/linear/JacobianFactor.h
@@ -187,7 +187,7 @@ namespace gtsam {
/// Return the diagonal of the Hessian for this factor
virtual VectorValues hessianDiagonal() const;
- /* ************************************************************************* */
+ /// Raw memory access version of hessianDiagonal
virtual void hessianDiagonal(double* d) const;
/// Return the block diagonal of the Hessian for this factor
diff --git a/gtsam/linear/Preconditioner.cpp b/gtsam/linear/Preconditioner.cpp
index c180f1160..9af362fba 100644
--- a/gtsam/linear/Preconditioner.cpp
+++ b/gtsam/linear/Preconditioner.cpp
@@ -13,6 +13,7 @@
#include
#include
#include
+#include
#include
#include
@@ -134,30 +135,16 @@ void BlockJacobiPreconditioner::build(
size_t nnz = 0;
for ( size_t i = 0 ; i < n ; ++i ) {
const size_t dim = dims_[i];
- blocks.push_back(Matrix::Zero(dim, dim));
+ // blocks.push_back(Matrix::Zero(dim, dim));
// nnz += (((dim)*(dim+1)) >> 1); // d*(d+1) / 2 ;
nnz += dim*dim;
}
- /* compute the block diagonal by scanning over the factors */
+ /* getting the block diagonals over the factors */
BOOST_FOREACH ( const GaussianFactor::shared_ptr &gf, gfg ) {
- if ( JacobianFactor::shared_ptr jf = boost::dynamic_pointer_cast(gf) ) {
- for ( JacobianFactor::const_iterator it = jf->begin() ; it != jf->end() ; ++it ) {
- const KeyInfoEntry &entry = keyInfo.find(*it)->second;
- const Matrix &Ai = jf->getA(it);
- blocks[entry.index()] += (Ai.transpose() * Ai);
- }
- }
- else if ( HessianFactor::shared_ptr hf = boost::dynamic_pointer_cast(gf) ) {
- for ( HessianFactor::const_iterator it = hf->begin() ; it != hf->end() ; ++it ) {
- const KeyInfoEntry &entry = keyInfo.find(*it)->second;
- const Matrix &Hii = hf->info(it, it).selfadjointView();
- blocks[entry.index()] += Hii;
- }
- }
- else {
- throw invalid_argument("BlockJacobiPreconditioner::build gfg contains a factor that is neither a JacobianFactor nor a HessianFactor.");
- }
+ std::map hessianMap = gf->hessianBlockDiagonal();
+ BOOST_FOREACH ( const Matrix hessian, hessianMap | boost::adaptors::map_values)
+ blocks.push_back(hessian);
}
/* if necessary, allocating the memory for cacheing the factorization results */
diff --git a/gtsam/navigation/ImuFactor.h b/gtsam/navigation/ImuFactor.h
index a8771f995..318c1aa5b 100644
--- a/gtsam/navigation/ImuFactor.h
+++ b/gtsam/navigation/ImuFactor.h
@@ -57,7 +57,7 @@ namespace gtsam {
class PreintegratedMeasurements {
public:
imuBias::ConstantBias biasHat; ///< Acceleration and angular rate bias values used during preintegration
- Matrix measurementCovariance; ///< (Raw measurements uncertainty) Covariance of the vector [integrationError measuredAcc measuredOmega] in R^(9X9)
+ Matrix measurementCovariance; ///< (continuous-time uncertainty) Covariance of the vector [integrationError measuredAcc measuredOmega] in R^(9X9)
Vector3 deltaPij; ///< Preintegrated relative position (does not take into account velocity at time i, see deltap+, in [2]) (in frame i)
Vector3 deltaVij; ///< Preintegrated relative velocity (in global frame)
@@ -216,11 +216,21 @@ namespace gtsam {
H_vel_pos, H_vel_vel, H_vel_angles,
H_angles_pos, H_angles_vel, H_angles_angles;
-
- // first order uncertainty propagation
+ // first order uncertainty propagation:
// the deltaT allows to pass from continuous time noise to discrete time noise
+ // measurementCovariance_discrete = measurementCovariance_contTime * (1/deltaT)
+ // Gt * Qt * G =(approx)= measurementCovariance_discrete * deltaT^2 = measurementCovariance_contTime * deltaT
PreintMeasCov = F * PreintMeasCov * F.transpose() + measurementCovariance * deltaT ;
+ // Extended version, without approximation: Gt * Qt * G =(approx)= measurementCovariance_contTime * deltaT
+ //
+ // Matrix G(9,9);
+ // G << I_3x3 * deltaT, Z_3x3, Z_3x3,
+ // Z_3x3, deltaRij.matrix() * deltaT, Z_3x3,
+ // Z_3x3, Z_3x3, Jrinv_theta_j * Jr_theta_incr * deltaT;
+ //
+ // PreintMeasCov = F * PreintMeasCov * F.transpose() + G * (1/deltaT) * measurementCovariance * G.transpose();
+
// Update preintegrated measurements
/* ----------------------------------------------------------------------------------------------------------------------- */
if(!use2ndOrderIntegration_){
diff --git a/gtsam/nonlinear/Values-inl.h b/gtsam/nonlinear/Values-inl.h
index 0d559cfe6..4595a70ed 100644
--- a/gtsam/nonlinear/Values-inl.h
+++ b/gtsam/nonlinear/Values-inl.h
@@ -215,7 +215,7 @@ namespace gtsam {
Values::Values(const Values::ConstFiltered& view) {
BOOST_FOREACH(const typename ConstFiltered::KeyValuePair& key_value, view) {
Key key = key_value.key;
- insert(key, key_value.value);
+ insert(key, static_cast(key_value.value));
}
}
diff --git a/gtsam/nonlinear/tests/testValues.cpp b/gtsam/nonlinear/tests/testValues.cpp
index 60639e8b7..219cacfd1 100644
--- a/gtsam/nonlinear/tests/testValues.cpp
+++ b/gtsam/nonlinear/tests/testValues.cpp
@@ -383,6 +383,12 @@ TEST(Values, filter) {
expectedSubValues1.insert(3, pose3);
EXPECT(assert_equal(expectedSubValues1, actualSubValues1));
+ // ConstFilter by Key
+ Values::ConstFiltered constfiltered = values.filter(boost::bind(std::greater_equal(), _1, 2));
+ EXPECT_LONGS_EQUAL(2, (long)constfiltered.size());
+ Values fromconstfiltered(constfiltered);
+ EXPECT(assert_equal(expectedSubValues1, fromconstfiltered));
+
// Filter by type
i = 0;
Values::ConstFiltered pose_filtered = values.filter();
diff --git a/gtsam/slam/ImplicitSchurFactor.h b/gtsam/slam/RegularImplicitSchurFactor.h
similarity index 92%
rename from gtsam/slam/ImplicitSchurFactor.h
rename to gtsam/slam/RegularImplicitSchurFactor.h
index 68f23e339..ecf8adfec 100644
--- a/gtsam/slam/ImplicitSchurFactor.h
+++ b/gtsam/slam/RegularImplicitSchurFactor.h
@@ -1,5 +1,5 @@
/**
- * @file ImplicitSchurFactor.h
+ * @file RegularImplicitSchurFactor.h
* @brief A new type of linear factor (GaussianFactor), which is subclass of GaussianFactor
* @author Frank Dellaert
* @author Luca Carlone
@@ -17,13 +17,13 @@
namespace gtsam {
/**
- * ImplicitSchurFactor
+ * RegularImplicitSchurFactor
*/
template //
-class ImplicitSchurFactor: public GaussianFactor {
+class RegularImplicitSchurFactor: public GaussianFactor {
public:
- typedef ImplicitSchurFactor This; ///< Typedef to this class
+ typedef RegularImplicitSchurFactor This; ///< Typedef to this class
typedef boost::shared_ptr shared_ptr; ///< shared_ptr to this class
protected:
@@ -40,11 +40,11 @@ protected:
public:
/// Constructor
- ImplicitSchurFactor() {
+ RegularImplicitSchurFactor() {
}
/// Construct from blcoks of F, E, inv(E'*E), and RHS vector b
- ImplicitSchurFactor(const std::vector& Fblocks, const Matrix& E,
+ RegularImplicitSchurFactor(const std::vector& Fblocks, const Matrix& E,
const Matrix3& P, const Vector& b) :
Fblocks_(Fblocks), PointCovariance_(P), E_(E), b_(b) {
initKeys();
@@ -58,7 +58,7 @@ public:
}
/// Destructor
- virtual ~ImplicitSchurFactor() {
+ virtual ~RegularImplicitSchurFactor() {
}
// Write access, only use for construction!
@@ -87,7 +87,7 @@ public:
/// print
void print(const std::string& s = "",
const KeyFormatter& keyFormatter = DefaultKeyFormatter) const {
- std::cout << " ImplicitSchurFactor " << std::endl;
+ std::cout << " RegularImplicitSchurFactor " << std::endl;
Factor::print(s);
std::cout << " PointCovariance_ \n" << PointCovariance_ << std::endl;
std::cout << " E_ \n" << E_ << std::endl;
@@ -96,7 +96,7 @@ public:
/// equals
bool equals(const GaussianFactor& lf, double tol) const {
- if (!dynamic_cast(&lf))
+ if (!dynamic_cast(&lf))
return false;
else {
return false;
@@ -110,21 +110,21 @@ public:
virtual Matrix augmentedJacobian() const {
throw std::runtime_error(
- "ImplicitSchurFactor::augmentedJacobian non implemented");
+ "RegularImplicitSchurFactor::augmentedJacobian non implemented");
return Matrix();
}
virtual std::pair jacobian() const {
- throw std::runtime_error("ImplicitSchurFactor::jacobian non implemented");
+ throw std::runtime_error("RegularImplicitSchurFactor::jacobian non implemented");
return std::make_pair(Matrix(), Vector());
}
virtual Matrix augmentedInformation() const {
throw std::runtime_error(
- "ImplicitSchurFactor::augmentedInformation non implemented");
+ "RegularImplicitSchurFactor::augmentedInformation non implemented");
return Matrix();
}
virtual Matrix information() const {
throw std::runtime_error(
- "ImplicitSchurFactor::information non implemented");
+ "RegularImplicitSchurFactor::information non implemented");
return Matrix();
}
@@ -210,18 +210,18 @@ public:
}
virtual GaussianFactor::shared_ptr clone() const {
- return boost::make_shared >(Fblocks_,
+ return boost::make_shared >(Fblocks_,
PointCovariance_, E_, b_);
- throw std::runtime_error("ImplicitSchurFactor::clone non implemented");
+ throw std::runtime_error("RegularImplicitSchurFactor::clone non implemented");
}
virtual bool empty() const {
return false;
}
virtual GaussianFactor::shared_ptr negate() const {
- return boost::make_shared >(Fblocks_,
+ return boost::make_shared >(Fblocks_,
PointCovariance_, E_, b_);
- throw std::runtime_error("ImplicitSchurFactor::negate non implemented");
+ throw std::runtime_error("RegularImplicitSchurFactor::negate non implemented");
}
// Raw Vector version of y += F'*alpha*(I - E*P*E')*F*x, for testing
@@ -454,7 +454,7 @@ public:
}
};
-// ImplicitSchurFactor
+// RegularImplicitSchurFactor
}
diff --git a/gtsam/slam/SmartFactorBase.h b/gtsam/slam/SmartFactorBase.h
index 1ecaaf170..8ae26e7cd 100644
--- a/gtsam/slam/SmartFactorBase.h
+++ b/gtsam/slam/SmartFactorBase.h
@@ -21,7 +21,7 @@
#include "JacobianFactorQ.h"
#include "JacobianFactorSVD.h"
-#include "ImplicitSchurFactor.h"
+#include "RegularImplicitSchurFactor.h"
#include "RegularHessianFactor.h"
#include
@@ -73,7 +73,7 @@ public:
/**
* Constructor
- * @param body_P_sensor is the transform from body to sensor frame (default identity)
+ * @param body_P_sensor is the transform from sensor to body frame (default identity)
*/
SmartFactorBase(boost::optional body_P_sensor = boost::none) :
body_P_sensor_(body_P_sensor) {
@@ -271,8 +271,13 @@ public:
Vector bi;
try {
- bi =
- -(cameras[i].project(point, Fi, Ei, Hcali) - this->measured_.at(i)).vector();
+ bi = -(cameras[i].project(point, Fi, Ei, Hcali) - this->measured_.at(i)).vector();
+ if(body_P_sensor_){
+ Pose3 w_Pose_body = (cameras[i].pose()).compose(body_P_sensor_->inverse());
+ Matrix J(6, 6);
+ Pose3 world_P_body = w_Pose_body.compose(*body_P_sensor_, J);
+ Fi = Fi * J;
+ }
} catch (CheiralityException&) {
std::cout << "Cheirality exception " << std::endl;
exit(EXIT_FAILURE);
@@ -624,11 +629,11 @@ public:
}
// ****************************************************************************************************
- boost::shared_ptr > createImplicitSchurFactor(
+ boost::shared_ptr > createRegularImplicitSchurFactor(
const Cameras& cameras, const Point3& point, double lambda = 0.0,
bool diagonalDamping = false) const {
- typename boost::shared_ptr > f(
- new ImplicitSchurFactor());
+ typename boost::shared_ptr > f(
+ new RegularImplicitSchurFactor());
computeJacobians(f->Fblocks(), f->E(), f->PointCovariance(), f->b(),
cameras, point, lambda, diagonalDamping);
f->initKeys();
diff --git a/gtsam/slam/SmartProjectionFactor.h b/gtsam/slam/SmartProjectionFactor.h
index 043528fea..bfd73d9fb 100644
--- a/gtsam/slam/SmartProjectionFactor.h
+++ b/gtsam/slam/SmartProjectionFactor.h
@@ -120,7 +120,7 @@ public:
* @param manageDegeneracy is true, in presence of degenerate triangulation, the factor is converted to a rotation-only constraint,
* otherwise the factor is simply neglected
* @param enableEPI if set to true linear triangulation is refined with embedded LM iterations
- * @param body_P_sensor is the transform from body to sensor frame (default identity)
+ * @param body_P_sensor is the transform from sensor to body frame (default identity)
*/
SmartProjectionFactor(const double rankTol, const double linThreshold,
const bool manageDegeneracy, const bool enableEPI,
@@ -298,7 +298,7 @@ public:
|| (!this->manageDegeneracy_
&& (this->cheiralityException_ || this->degenerate_))) {
if (isDebug) {
- std::cout << "createImplicitSchurFactor: degenerate configuration"
+ std::cout << "createRegularImplicitSchurFactor: degenerate configuration"
<< std::endl;
}
return false;
@@ -409,12 +409,12 @@ public:
}
// create factor
- boost::shared_ptr > createImplicitSchurFactor(
+ boost::shared_ptr > createRegularImplicitSchurFactor(
const Cameras& cameras, double lambda) const {
if (triangulateForLinearize(cameras))
- return Base::createImplicitSchurFactor(cameras, point_, lambda);
+ return Base::createRegularImplicitSchurFactor(cameras, point_, lambda);
else
- return boost::shared_ptr >();
+ return boost::shared_ptr >();
}
/// create factor
@@ -685,7 +685,7 @@ public:
inline bool isPointBehindCamera() const {
return cheiralityException_;
}
- /** return chirality verbosity */
+ /** return cheirality verbosity */
inline bool verboseCheirality() const {
return verboseCheirality_;
}
diff --git a/gtsam/slam/SmartProjectionPoseFactor.h b/gtsam/slam/SmartProjectionPoseFactor.h
index 273102bda..f871ab3aa 100644
--- a/gtsam/slam/SmartProjectionPoseFactor.h
+++ b/gtsam/slam/SmartProjectionPoseFactor.h
@@ -63,7 +63,7 @@ public:
* @param manageDegeneracy is true, in presence of degenerate triangulation, the factor is converted to a rotation-only constraint,
* otherwise the factor is simply neglected
* @param enableEPI if set to true linear triangulation is refined with embedded LM iterations
- * @param body_P_sensor is the transform from body to sensor frame (default identity)
+ * @param body_P_sensor is the transform from sensor to body frame (default identity)
*/
SmartProjectionPoseFactor(const double rankTol = 1,
const double linThreshold = -1, const bool manageDegeneracy = false,
@@ -157,6 +157,9 @@ public:
size_t i=0;
BOOST_FOREACH(const Key& k, this->keys_) {
Pose3 pose = values.at(k);
+ if(Base::body_P_sensor_)
+ pose = pose.compose(*(Base::body_P_sensor_));
+
typename Base::Camera camera(pose, *K_all_[i++]);
cameras.push_back(camera);
}
diff --git a/gtsam/slam/tests/testImplicitSchurFactor.cpp b/gtsam/slam/tests/testRegularImplicitSchurFactor.cpp
similarity index 94%
rename from gtsam/slam/tests/testImplicitSchurFactor.cpp
rename to gtsam/slam/tests/testRegularImplicitSchurFactor.cpp
index 77faaacc1..0161d4fb6 100644
--- a/gtsam/slam/tests/testImplicitSchurFactor.cpp
+++ b/gtsam/slam/tests/testRegularImplicitSchurFactor.cpp
@@ -6,7 +6,7 @@
*/
//#include
-#include
+#include
//#include
#include
//#include "gtsam_unstable/slam/JacobianFactorQR.h"
@@ -38,19 +38,19 @@ const vector > Fblocks = list_of > //
const Vector b = (Vector(6) << 1., 2., 3., 4., 5., 6.);
//*************************************************************************************
-TEST( implicitSchurFactor, creation ) {
+TEST( regularImplicitSchurFactor, creation ) {
// Matrix E = Matrix::Ones(6,3);
Matrix E = zeros(6, 3);
E.block<2,2>(0, 0) = eye(2);
E.block<2,3>(2, 0) = 2 * ones(2, 3);
Matrix3 P = (E.transpose() * E).inverse();
- ImplicitSchurFactor<6> expected(Fblocks, E, P, b);
+ RegularImplicitSchurFactor<6> expected(Fblocks, E, P, b);
Matrix expectedP = expected.getPointCovariance();
EXPECT(assert_equal(expectedP, P));
}
/* ************************************************************************* */
-TEST( implicitSchurFactor, addHessianMultiply ) {
+TEST( regularImplicitSchurFactor, addHessianMultiply ) {
Matrix E = zeros(6, 3);
E.block<2,2>(0, 0) = eye(2);
@@ -81,11 +81,11 @@ TEST( implicitSchurFactor, addHessianMultiply ) {
keys += 0,1,2,3;
Vector x = xvalues.vector(keys);
Vector expected = zero(24);
- ImplicitSchurFactor<6>::multiplyHessianAdd(F, E, P, alpha, x, expected);
+ RegularImplicitSchurFactor<6>::multiplyHessianAdd(F, E, P, alpha, x, expected);
EXPECT(assert_equal(expected, yExpected.vector(keys), 1e-8));
// Create ImplicitSchurFactor
- ImplicitSchurFactor<6> implicitFactor(Fblocks, E, P, b);
+ RegularImplicitSchurFactor<6> implicitFactor(Fblocks, E, P, b);
VectorValues zero = 0 * yExpected;// quick way to get zero w right structure
{ // First Version
@@ -190,7 +190,7 @@ TEST( implicitSchurFactor, addHessianMultiply ) {
}
/* ************************************************************************* */
-TEST(implicitSchurFactor, hessianDiagonal)
+TEST(regularImplicitSchurFactor, hessianDiagonal)
{
/* TESTED AGAINST MATLAB
* F = [ones(2,6) zeros(2,6) zeros(2,6)
@@ -207,7 +207,7 @@ TEST(implicitSchurFactor, hessianDiagonal)
E.block<2,3>(2, 0) << 1,2,3,4,5,6;
E.block<2,3>(4, 0) << 0.5,1,2,3,4,5;
Matrix3 P = (E.transpose() * E).inverse();
- ImplicitSchurFactor<6> factor(Fblocks, E, P, b);
+ RegularImplicitSchurFactor<6> factor(Fblocks, E, P, b);
// hessianDiagonal
VectorValues expected;
diff --git a/gtsam/slam/tests/testSmartProjectionPoseFactor.cpp b/gtsam/slam/tests/testSmartProjectionPoseFactor.cpp
index 6b849ba5a..4e4fde3e4 100644
--- a/gtsam/slam/tests/testSmartProjectionPoseFactor.cpp
+++ b/gtsam/slam/tests/testSmartProjectionPoseFactor.cpp
@@ -292,6 +292,95 @@ TEST( SmartProjectionPoseFactor, 3poses_smart_projection_factor ){
if(isDebugTest) tictoc_print_();
}
+/* *************************************************************************/
+TEST( SmartProjectionPoseFactor, smartFactorWithSensorBodyTransform ){
+
+ // create first camera. Looking along X-axis, 1 meter above ground plane (x-y)
+ Pose3 cameraPose1 = Pose3(Rot3::ypr(-M_PI/2, 0., -M_PI/2), gtsam::Point3(0,0,1)); // body poses
+ Pose3 cameraPose2 = cameraPose1 * Pose3(Rot3(), Point3(1,0,0));
+ Pose3 cameraPose3 = cameraPose1 * Pose3(Rot3(), Point3(0,-1,0));
+
+ SimpleCamera cam1(cameraPose1, *K); // with camera poses
+ SimpleCamera cam2(cameraPose2, *K);
+ SimpleCamera cam3(cameraPose3, *K);
+
+ // create arbitrary body_Pose_sensor (transforms from sensor to body)
+ Pose3 sensor_to_body = Pose3(Rot3::ypr(-M_PI/2, 0., -M_PI/2), gtsam::Point3(1, 1, 1)); // Pose3(); //
+
+ // These are the poses we want to estimate, from camera measurements
+ Pose3 bodyPose1 = cameraPose1.compose(sensor_to_body.inverse());
+ Pose3 bodyPose2 = cameraPose2.compose(sensor_to_body.inverse());
+ Pose3 bodyPose3 = cameraPose3.compose(sensor_to_body.inverse());
+
+ // three landmarks ~5 meters infront of camera
+ Point3 landmark1(5, 0.5, 1.2);
+ Point3 landmark2(5, -0.5, 1.2);
+ Point3 landmark3(5, 0, 3.0);
+
+ vector measurements_cam1, measurements_cam2, measurements_cam3;
+
+ // Project three landmarks into three cameras
+ projectToMultipleCameras(cam1, cam2, cam3, landmark1, measurements_cam1);
+ projectToMultipleCameras(cam1, cam2, cam3, landmark2, measurements_cam2);
+ projectToMultipleCameras(cam1, cam2, cam3, landmark3, measurements_cam3);
+
+ // Create smart factors
+ std::vector views;
+ views.push_back(x1);
+ views.push_back(x2);
+ views.push_back(x3);
+
+ double rankTol = 1;
+ double linThreshold = -1;
+ bool manageDegeneracy = false;
+ bool enableEPI = false;
+
+ SmartFactor::shared_ptr smartFactor1(new SmartFactor(rankTol,linThreshold,manageDegeneracy,enableEPI,sensor_to_body));
+ smartFactor1->add(measurements_cam1, views, model, K);
+
+ SmartFactor::shared_ptr smartFactor2(new SmartFactor(rankTol,linThreshold,manageDegeneracy,enableEPI,sensor_to_body));
+ smartFactor2->add(measurements_cam2, views, model, K);
+
+ SmartFactor::shared_ptr smartFactor3(new SmartFactor(rankTol,linThreshold,manageDegeneracy,enableEPI,sensor_to_body));
+ smartFactor3->add(measurements_cam3, views, model, K);
+
+ const SharedDiagonal noisePrior = noiseModel::Isotropic::Sigma(6, 0.10);
+
+ // Put all factors in factor graph, adding priors
+ NonlinearFactorGraph graph;
+ graph.push_back(smartFactor1);
+ graph.push_back(smartFactor2);
+ graph.push_back(smartFactor3);
+ graph.push_back(PriorFactor(x1, bodyPose1, noisePrior));
+ graph.push_back(PriorFactor(x2, bodyPose2, noisePrior));
+
+ // Check errors at ground truth poses
+ Values gtValues;
+ gtValues.insert(x1, bodyPose1);
+ gtValues.insert(x2, bodyPose2);
+ gtValues.insert(x3, bodyPose3);
+ double actualError = graph.error(gtValues);
+ double expectedError = 0.0;
+ DOUBLES_EQUAL(expectedError, actualError, 1e-7)
+
+ Pose3 noise_pose = Pose3(Rot3::ypr(-M_PI/100, 0., -M_PI/100), gtsam::Point3(0.1,0.1,0.1));
+ Values values;
+ values.insert(x1, bodyPose1);
+ values.insert(x2, bodyPose2);
+ // initialize third pose with some noise, we expect it to move back to original pose3
+ values.insert(x3, bodyPose3*noise_pose);
+
+ LevenbergMarquardtParams params;
+ Values result;
+ LevenbergMarquardtOptimizer optimizer(graph, values, params);
+ result = optimizer.optimize();
+
+ // result.print("results of 3 camera, 3 landmark optimization \n");
+ if(isDebugTest) result.at(x3).print("Smart: Pose3 after optimization: ");
+ EXPECT(assert_equal(bodyPose3,result.at(x3)));
+}
+
+
/* *************************************************************************/
TEST( SmartProjectionPoseFactor, 3poses_iterative_smart_projection_factor ){
// cout << " ************************ SmartProjectionPoseFactor: 3 cams + 3 landmarks **********************" << endl;
diff --git a/gtsam_unstable/gtsam_unstable.h b/gtsam_unstable/gtsam_unstable.h
index 27460bd72..95e635516 100644
--- a/gtsam_unstable/gtsam_unstable.h
+++ b/gtsam_unstable/gtsam_unstable.h
@@ -4,14 +4,15 @@
// specify the classes from gtsam we are using
virtual class gtsam::Value;
-virtual class gtsam::LieScalar;
-virtual class gtsam::LieVector;
-virtual class gtsam::Point2;
-virtual class gtsam::Rot2;
-virtual class gtsam::Pose2;
-virtual class gtsam::Point3;
-virtual class gtsam::Rot3;
-virtual class gtsam::Pose3;
+class gtsam::Vector6;
+class gtsam::LieScalar;
+class gtsam::LieVector;
+class gtsam::Point2;
+class gtsam::Rot2;
+class gtsam::Pose2;
+class gtsam::Point3;
+class gtsam::Rot3;
+class gtsam::Pose3;
virtual class gtsam::noiseModel::Base;
virtual class gtsam::noiseModel::Gaussian;
virtual class gtsam::imuBias::ConstantBias;
@@ -23,8 +24,8 @@ virtual class gtsam::NoiseModelFactor4;
virtual class gtsam::GaussianFactor;
virtual class gtsam::HessianFactor;
virtual class gtsam::JacobianFactor;
-virtual class gtsam::Cal3_S2;
-virtual class gtsam::Cal3DS2;
+class gtsam::Cal3_S2;
+class gtsam::Cal3DS2;
class gtsam::GaussianFactorGraph;
class gtsam::NonlinearFactorGraph;
class gtsam::Ordering;
@@ -48,7 +49,7 @@ class Dummy {
};
#include
-virtual class PoseRTV : gtsam::Value {
+class PoseRTV {
PoseRTV();
PoseRTV(Vector rtv);
PoseRTV(const gtsam::Point3& pt, const gtsam::Rot3& rot, const gtsam::Point3& vel);
@@ -99,7 +100,7 @@ virtual class PoseRTV : gtsam::Value {
};
#include
-virtual class Pose3Upright : gtsam::Value {
+class Pose3Upright {
Pose3Upright();
Pose3Upright(const gtsam::Pose3Upright& x);
Pose3Upright(const gtsam::Rot2& bearing, const gtsam::Point3& t);
@@ -137,7 +138,7 @@ virtual class Pose3Upright : gtsam::Value {
}; // \class Pose3Upright
#include
-virtual class BearingS2 : gtsam::Value {
+class BearingS2 {
BearingS2();
BearingS2(double azimuth, double elevation);
BearingS2(const gtsam::Rot2& azimuth, const gtsam::Rot2& elevation);
@@ -520,14 +521,14 @@ virtual class PendulumFactorPk1 : gtsam::NonlinearFactor {
virtual class Reconstruction : gtsam::NonlinearFactor {
Reconstruction(size_t gKey1, size_t gKey, size_t xiKey, double h);
- Vector evaluateError(const gtsam::Pose3& gK1, const gtsam::Pose3& gK, const gtsam::LieVector& xiK) const;
+ Vector evaluateError(const gtsam::Pose3& gK1, const gtsam::Pose3& gK, const gtsam::Vector6& xiK) const;
};
virtual class DiscreteEulerPoincareHelicopter : gtsam::NonlinearFactor {
DiscreteEulerPoincareHelicopter(size_t xiKey, size_t xiKey_1, size_t gKey,
double h, Matrix Inertia, Vector Fu, double m);
- Vector evaluateError(const gtsam::LieVector& xiK, const gtsam::LieVector& xiK_1, const gtsam::Pose3& gK) const;
+ Vector evaluateError(const gtsam::Vector6& xiK, const gtsam::Vector6& xiK_1, const gtsam::Pose3& gK) const;
};
//*************************************************************************
diff --git a/gtsam_unstable/nonlinear/Expression-inl.h b/gtsam_unstable/nonlinear/Expression-inl.h
index f49b985ba..40fc54751 100644
--- a/gtsam_unstable/nonlinear/Expression-inl.h
+++ b/gtsam_unstable/nonlinear/Expression-inl.h
@@ -64,7 +64,8 @@ public:
}
/// Access via key
VerticalBlockMatrix::Block operator()(Key key) {
- FastVector::const_iterator it = std::find(keys_.begin(),keys_.end(),key);
+ FastVector::const_iterator it = std::find(keys_.begin(), keys_.end(),
+ key);
DenseIndex block = it - keys_.begin();
return Ab_(block);
}
@@ -98,7 +99,7 @@ struct CallRecord {
template
void handleLeafCase(const Eigen::Matrix& dTdA,
JacobianMap& jacobians, Key key) {
- jacobians(key).block < ROWS, COLS > (0, 0) += dTdA; // block makes HUGE difference
+ jacobians(key).block(0, 0) += dTdA; // block makes HUGE difference
}
/// Handle Leaf Case for Dynamic Matrix type (slower)
template<>
@@ -311,9 +312,9 @@ public:
//-----------------------------------------------------------------------------
/// Leaf Expression
-template >
+template >
class LeafExpression: public ExpressionNode {
- typedef ChartValue value_type; // perhaps this can be something else like a std::pair ??
+ typedef ChartValue value_type; // perhaps this can be something else like a std::pair ??
/// The key into values
Key key_;
@@ -347,8 +348,8 @@ public:
}
/// Construct an execution trace for reverse AD
- virtual const value_type& traceExecution(const Values& values, ExecutionTrace& trace,
- char* raw) const {
+ virtual const value_type& traceExecution(const Values& values,
+ ExecutionTrace& trace, char* raw) const {
trace.setLeaf(key_);
return dynamic_cast(values.at(key_));
}
@@ -358,7 +359,7 @@ public:
//-----------------------------------------------------------------------------
/// Leaf Expression, if no chart is given, assume default chart and value_type is just the plain value
template
-class LeafExpression >: public ExpressionNode {
+class LeafExpression > : public ExpressionNode {
typedef T value_type;
/// The key into values
@@ -405,6 +406,7 @@ public:
// C++ Template Metaprogramming: Concepts, Tools, and Techniques from Boost
// and Beyond. Abrahams, David; Gurtovoy, Aleksey. Pearson Education.
// to recursively generate a class, that will be the base for function nodes.
+//
// The class generated, for three arguments A1, A2, and A3 will be
//
// struct Base1 : Argument, FunctionalBase {
@@ -429,6 +431,30 @@ public:
//
// All this magic happens when we generate the Base3 base class of FunctionalNode
// by invoking boost::mpl::fold over the meta-function GenerateFunctionalNode
+//
+// Similarly, the inner Record struct will be
+//
+// struct Record1 : JacobianTrace, CallRecord::value> {
+// ... storage related to A1 ...
+// ... methods that work on A1 ...
+// };
+//
+// struct Record2 : JacobianTrace, Record1 {
+// ... storage related to A2 ...
+// ... methods that work on A2 and (recursively) on A1 ...
+// };
+//
+// struct Record3 : JacobianTrace, Record2 {
+// ... storage related to A3 ...
+// ... methods that work on A3 and (recursively) on A2 and A1 ...
+// };
+//
+// struct Record : Record3 {
+// Provides convenience access to storage in hierarchy by using
+// static_cast &>(*this)
+// }
+//
+
//-----------------------------------------------------------------------------
/// meta-function to generate fixed-size JacobianTA type
@@ -457,6 +483,7 @@ struct FunctionalBase: ExpressionNode {
/// Construct an execution trace for reverse AD
void trace(const Values& values, Record* record, char*& raw) const {
+ // base case: does not do anything
}
};
@@ -562,15 +589,23 @@ struct GenerateFunctionalNode: Argument, Base {
template
struct FunctionalNode {
+ /// The following typedef generates the recursively defined Base class
typedef typename boost::mpl::fold,
GenerateFunctionalNode >::type Base;
+ /**
+ * The type generated by this meta-function derives from Base
+ * and adds access functions as well as the crucial [trace] function
+ */
struct type: public Base {
// Argument types and derived, note these are base 0 !
+ // These are currently not used - useful for Phoenix in future
+#ifdef EXPRESSIONS_PHOENIX
typedef TYPES Arguments;
typedef typename boost::mpl::transform >::type Jacobians;
typedef typename boost::mpl::transform >::type Optionals;
+#endif
/// Reset expression shared pointer
template
@@ -725,7 +760,8 @@ public:
typedef boost::function<
T(const A1&, const A2&, const A3&, typename OptionalJacobian::type,
- typename OptionalJacobian::type, typename OptionalJacobian::type)> Function;
+ typename OptionalJacobian::type,
+ typename OptionalJacobian::type)> Function;
typedef typename FunctionalNode >::type Base;
typedef typename Base::Record Record;
diff --git a/gtsam_unstable/nonlinear/tests/testExpressionMeta.cpp b/gtsam_unstable/nonlinear/tests/testExpressionMeta.cpp
index 84a1ca720..b2cdcdf34 100644
--- a/gtsam_unstable/nonlinear/tests/testExpressionMeta.cpp
+++ b/gtsam_unstable/nonlinear/tests/testExpressionMeta.cpp
@@ -22,6 +22,7 @@
#include
#include
+#include
using namespace std;
using namespace gtsam;
@@ -143,7 +144,7 @@ TEST(ExpressionFactor, Triple) {
// Test out invoke
TEST(ExpressionFactor, Invoke) {
- assert(invoke(add,boost::fusion::make_vector(1,1)) == 2);
+ EXPECT_LONGS_EQUAL(2, invoke(plus(),boost::fusion::make_vector(1,1)));
// Creating a Pose3 (is there another way?)
boost::fusion::vector pair;
diff --git a/gtsam_unstable/partition/FindSeparator-inl.h b/gtsam_unstable/partition/FindSeparator-inl.h
index 849cf7a05..ace13dc41 100644
--- a/gtsam_unstable/partition/FindSeparator-inl.h
+++ b/gtsam_unstable/partition/FindSeparator-inl.h
@@ -17,12 +17,14 @@
#include
#include
+#include "FindSeparator.h"
+
extern "C" {
#include
#include "metislib.h"
}
-#include "FindSeparator.h"
+
namespace gtsam { namespace partition {
diff --git a/gtsam_unstable/slam/BetweenFactorEM.h b/gtsam_unstable/slam/BetweenFactorEM.h
index c147552b3..9082c0101 100644
--- a/gtsam_unstable/slam/BetweenFactorEM.h
+++ b/gtsam_unstable/slam/BetweenFactorEM.h
@@ -25,384 +25,405 @@
namespace gtsam {
+/**
+ * A class for a measurement predicted by "between(config[key1],config[key2])"
+ * @tparam VALUE the Value type
+ * @addtogroup SLAM
+ */
+template
+class BetweenFactorEM: public NonlinearFactor {
+
+public:
+
+ typedef VALUE T;
+
+private:
+
+ typedef BetweenFactorEM This;
+ typedef gtsam::NonlinearFactor Base;
+
+ gtsam::Key key1_;
+ gtsam::Key key2_;
+
+ VALUE measured_; /** The measurement */
+
+ SharedGaussian model_inlier_;
+ SharedGaussian model_outlier_;
+
+ double prior_inlier_;
+ double prior_outlier_;
+
+ bool flag_bump_up_near_zero_probs_;
+
+ /** concept check by type */
+ GTSAM_CONCEPT_LIE_TYPE(T)GTSAM_CONCEPT_TESTABLE_TYPE(T)
+
+public:
+
+ // shorthand for a smart pointer to a factor
+ typedef typename boost::shared_ptr shared_ptr;
+
+ /** default constructor - only use for serialization */
+ BetweenFactorEM() {
+ }
+
+ /** Constructor */
+ BetweenFactorEM(Key key1, Key key2, const VALUE& measured,
+ const SharedGaussian& model_inlier, const SharedGaussian& model_outlier,
+ const double prior_inlier, const double prior_outlier,
+ const bool flag_bump_up_near_zero_probs = false) :
+ Base(cref_list_of<2>(key1)(key2)), key1_(key1), key2_(key2), measured_(
+ measured), model_inlier_(model_inlier), model_outlier_(model_outlier), prior_inlier_(
+ prior_inlier), prior_outlier_(prior_outlier), flag_bump_up_near_zero_probs_(
+ flag_bump_up_near_zero_probs) {
+ }
+
+ virtual ~BetweenFactorEM() {
+ }
+
+ /** implement functions needed for Testable */
+
+ /** print */
+ virtual void print(const std::string& s, const KeyFormatter& keyFormatter =
+ DefaultKeyFormatter) const {
+ std::cout << s << "BetweenFactorEM(" << keyFormatter(key1_) << ","
+ << keyFormatter(key2_) << ")\n";
+ measured_.print(" measured: ");
+ model_inlier_->print(" noise model inlier: ");
+ model_outlier_->print(" noise model outlier: ");
+ std::cout << "(prior_inlier, prior_outlier_) = (" << prior_inlier_ << ","
+ << prior_outlier_ << ")\n";
+ // Base::print(s, keyFormatter);
+ }
+
+ /** equals */
+ virtual bool equals(const NonlinearFactor& f, double tol = 1e-9) const {
+ const This *t = dynamic_cast(&f);
+
+ if (t && Base::equals(f))
+ return key1_ == t->key1_ && key2_ == t->key2_
+ &&
+ // model_inlier_->equals(t->model_inlier_ ) && // TODO: fix here
+ // model_outlier_->equals(t->model_outlier_ ) &&
+ prior_outlier_ == t->prior_outlier_
+ && prior_inlier_ == t->prior_inlier_ && measured_.equals(t->measured_);
+ else
+ return false;
+ }
+
+ /** implement functions needed to derive from Factor */
+
+ /* ************************************************************************* */
+ virtual double error(const gtsam::Values& x) const {
+ return whitenedError(x).squaredNorm();
+ }
+
+ /* ************************************************************************* */
/**
- * A class for a measurement predicted by "between(config[key1],config[key2])"
- * @tparam VALUE the Value type
- * @addtogroup SLAM
+ * Linearize a non-linearFactorN to get a gtsam::GaussianFactor,
+ * \f$ Ax-b \approx h(x+\delta x)-z = h(x) + A \delta x - z \f$
+ * Hence \f$ b = z - h(x) = - \mathtt{error\_vector}(x) \f$
*/
- template
- class BetweenFactorEM: public NonlinearFactor {
+ /* This version of linearize recalculates the noise model each time */
+ virtual boost::shared_ptr linearize(
+ const gtsam::Values& x) const {
+ // Only linearize if the factor is active
+ if (!this->active(x))
+ return boost::shared_ptr();
- public:
+ //std::cout<<"About to linearize"< A(this->size());
+ gtsam::Vector b = -whitenedError(x, A);
+ A1 = A[0];
+ A2 = A[1];
- typedef VALUE T;
+ return gtsam::GaussianFactor::shared_ptr(
+ new gtsam::JacobianFactor(key1_, A1, key2_, A2, b,
+ gtsam::noiseModel::Unit::Create(b.size())));
+ }
- private:
+ /* ************************************************************************* */
+ gtsam::Vector whitenedError(const gtsam::Values& x,
+ boost::optional&> H = boost::none) const {
- typedef BetweenFactorEM This;
- typedef gtsam::NonlinearFactor Base;
+ bool debug = true;
- gtsam::Key key1_;
- gtsam::Key key2_;
+ const T& p1 = x.at(key1_);
+ const T& p2 = x.at(key2_);
- VALUE measured_; /** The measurement */
+ Matrix H1, H2;
- SharedGaussian model_inlier_;
- SharedGaussian model_outlier_;
+ T hx = p1.between(p2, H1, H2); // h(x)
+ // manifold equivalent of h(x)-z -> log(z,h(x))
- double prior_inlier_;
- double prior_outlier_;
+ Vector err = measured_.localCoordinates(hx);
- bool flag_bump_up_near_zero_probs_;
+ // Calculate indicator probabilities (inlier and outlier)
+ Vector p_inlier_outlier = calcIndicatorProb(x);
+ double p_inlier = p_inlier_outlier[0];
+ double p_outlier = p_inlier_outlier[1];
- /** concept check by type */
- GTSAM_CONCEPT_LIE_TYPE(T)
- GTSAM_CONCEPT_TESTABLE_TYPE(T)
+ Vector err_wh_inlier = model_inlier_->whiten(err);
+ Vector err_wh_outlier = model_outlier_->whiten(err);
- public:
+ Matrix invCov_inlier = model_inlier_->R().transpose() * model_inlier_->R();
+ Matrix invCov_outlier = model_outlier_->R().transpose()
+ * model_outlier_->R();
- // shorthand for a smart pointer to a factor
- typedef typename boost::shared_ptr shared_ptr;
+ Vector err_wh_eq;
+ err_wh_eq.resize(err_wh_inlier.rows() * 2);
+ err_wh_eq << sqrt(p_inlier) * err_wh_inlier.array(), sqrt(p_outlier)
+ * err_wh_outlier.array();
- /** default constructor - only use for serialization */
- BetweenFactorEM() {}
+ if (H) {
+ // stack Jacobians for the two indicators for each of the key
- /** Constructor */
- BetweenFactorEM(Key key1, Key key2, const VALUE& measured,
- const SharedGaussian& model_inlier, const SharedGaussian& model_outlier,
- const double prior_inlier, const double prior_outlier, const bool flag_bump_up_near_zero_probs = false) :
- Base(cref_list_of<2>(key1)(key2)), key1_(key1), key2_(key2), measured_(measured),
- model_inlier_(model_inlier), model_outlier_(model_outlier),
- prior_inlier_(prior_inlier), prior_outlier_(prior_outlier), flag_bump_up_near_zero_probs_(flag_bump_up_near_zero_probs){
+ Matrix H1_inlier = sqrt(p_inlier) * model_inlier_->Whiten(H1);
+ Matrix H1_outlier = sqrt(p_outlier) * model_outlier_->Whiten(H1);
+ Matrix H1_aug = gtsam::stack(2, &H1_inlier, &H1_outlier);
+
+ Matrix H2_inlier = sqrt(p_inlier) * model_inlier_->Whiten(H2);
+ Matrix H2_outlier = sqrt(p_outlier) * model_outlier_->Whiten(H2);
+ Matrix H2_aug = gtsam::stack(2, &H2_inlier, &H2_outlier);
+
+ (*H)[0].resize(H1_aug.rows(), H1_aug.cols());
+ (*H)[1].resize(H2_aug.rows(), H2_aug.cols());
+
+ (*H)[0] = H1_aug;
+ (*H)[1] = H2_aug;
}
- virtual ~BetweenFactorEM() {}
+ if (debug) {
+ // std::cout<<"unwhitened error: "<print(" noise model inlier: ");
- model_outlier_->print(" noise model outlier: ");
- std::cout << "(prior_inlier, prior_outlier_) = ("
- << prior_inlier_ << ","
- << prior_outlier_ << ")\n";
- // Base::print(s, keyFormatter);
+ // Matrix Cov_inlier = invCov_inlier.inverse();
+ // Matrix Cov_outlier = invCov_outlier.inverse();
+ // std::cout<<"Cov_inlier: "< (&f);
+ return err_wh_eq;
+ }
- if(t && Base::equals(f))
- return key1_ == t->key1_ && key2_ == t->key2_ &&
- // model_inlier_->equals(t->model_inlier_ ) && // TODO: fix here
- // model_outlier_->equals(t->model_outlier_ ) &&
- prior_outlier_ == t->prior_outlier_ && prior_inlier_ == t->prior_inlier_ && measured_.equals(t->measured_);
- else
- return false;
+ /* ************************************************************************* */
+ gtsam::Vector calcIndicatorProb(const gtsam::Values& x) const {
+
+ bool debug = false;
+
+ Vector err = unwhitenedError(x);
+
+ // Calculate indicator probabilities (inlier and outlier)
+ Vector err_wh_inlier = model_inlier_->whiten(err);
+ Vector err_wh_outlier = model_outlier_->whiten(err);
+
+ Matrix invCov_inlier = model_inlier_->R().transpose() * model_inlier_->R();
+ Matrix invCov_outlier = model_outlier_->R().transpose()
+ * model_outlier_->R();
+
+ double p_inlier = prior_inlier_ * std::sqrt(invCov_inlier.determinant())
+ * exp(-0.5 * err_wh_inlier.dot(err_wh_inlier));
+ double p_outlier = prior_outlier_ * std::sqrt(invCov_outlier.determinant())
+ * exp(-0.5 * err_wh_outlier.dot(err_wh_outlier));
+
+ if (debug) {
+ std::cout << "in calcIndicatorProb. err_unwh: " << err[0] << ", "
+ << err[1] << ", " << err[2] << std::endl;
+ std::cout << "in calcIndicatorProb. err_wh_inlier: " << err_wh_inlier[0]
+ << ", " << err_wh_inlier[1] << ", " << err_wh_inlier[2] << std::endl;
+ std::cout << "in calcIndicatorProb. err_wh_inlier.dot(err_wh_inlier): "
+ << err_wh_inlier.dot(err_wh_inlier) << std::endl;
+ std::cout << "in calcIndicatorProb. err_wh_outlier.dot(err_wh_outlier): "
+ << err_wh_outlier.dot(err_wh_outlier) << std::endl;
+
+ std::cout
+ << "in calcIndicatorProb. p_inlier, p_outlier before normalization: "
+ << p_inlier << ", " << p_outlier << std::endl;
}
- /** implement functions needed to derive from Factor */
+ double sumP = p_inlier + p_outlier;
+ p_inlier /= sumP;
+ p_outlier /= sumP;
- /* ************************************************************************* */
- virtual double error(const gtsam::Values& x) const {
- return whitenedError(x).squaredNorm();
+ if (flag_bump_up_near_zero_probs_) {
+ // Bump up near-zero probabilities (as in linerFlow.h)
+ double minP = 0.05; // == 0.1 / 2 indicator variables
+ if (p_inlier < minP || p_outlier < minP) {
+ if (p_inlier < minP)
+ p_inlier = minP;
+ if (p_outlier < minP)
+ p_outlier = minP;
+ sumP = p_inlier + p_outlier;
+ p_inlier /= sumP;
+ p_outlier /= sumP;
+ }
}
- /* ************************************************************************* */
- /**
- * Linearize a non-linearFactorN to get a gtsam::GaussianFactor,
- * \f$ Ax-b \approx h(x+\delta x)-z = h(x) + A \delta x - z \f$
- * Hence \f$ b = z - h(x) = - \mathtt{error\_vector}(x) \f$
+ return (Vector(2) << p_inlier, p_outlier);
+ }
+
+ /* ************************************************************************* */
+ gtsam::Vector unwhitenedError(const gtsam::Values& x) const {
+
+ const T& p1 = x.at(key1_);
+ const T& p2 = x.at(key2_);
+
+ Matrix H1, H2;
+
+ T hx = p1.between(p2, H1, H2); // h(x)
+
+ return measured_.localCoordinates(hx);
+ }
+
+ /* ************************************************************************* */
+ void set_flag_bump_up_near_zero_probs(bool flag) {
+ flag_bump_up_near_zero_probs_ = flag;
+ }
+
+ /* ************************************************************************* */
+ bool get_flag_bump_up_near_zero_probs() const {
+ return flag_bump_up_near_zero_probs_;
+ }
+
+ /* ************************************************************************* */
+ SharedGaussian get_model_inlier() const {
+ return model_inlier_;
+ }
+
+ /* ************************************************************************* */
+ SharedGaussian get_model_outlier() const {
+ return model_outlier_;
+ }
+
+ /* ************************************************************************* */
+ Matrix get_model_inlier_cov() const {
+ return (model_inlier_->R().transpose() * model_inlier_->R()).inverse();
+ }
+
+ /* ************************************************************************* */
+ Matrix get_model_outlier_cov() const {
+ return (model_outlier_->R().transpose() * model_outlier_->R()).inverse();
+ }
+
+ /* ************************************************************************* */
+ void updateNoiseModels(const gtsam::Values& values,
+ const gtsam::NonlinearFactorGraph& graph) {
+ /* Update model_inlier_ and model_outlier_ to account for uncertainty in robot trajectories
+ * (note these are given in the E step, where indicator probabilities are calculated).
+ *
+ * Principle: R += [H1 H2] * joint_cov12 * [H1 H2]', where H1, H2 are Jacobians of the
+ * unwhitened error w.r.t. states, and R is the measurement covariance (inlier or outlier modes).
+ *
+ * TODO: improve efficiency (info form)
*/
- /* This version of linearize recalculates the noise model each time */
- virtual boost::shared_ptr linearize(const gtsam::Values& x) const {
- // Only linearize if the factor is active
- if (!this->active(x))
- return boost::shared_ptr();
- //std::cout<<"About to linearize"< A(this->size());
- gtsam::Vector b = -whitenedError(x, A);
- A1 = A[0];
- A2 = A[1];
+ // get joint covariance of the involved states
+ std::vector Keys;
+ Keys.push_back(key1_);
+ Keys.push_back(key2_);
+ Marginals marginals(graph, values, Marginals::QR);
+ JointMarginal joint_marginal12 = marginals.jointMarginalCovariance(Keys);
+ Matrix cov1 = joint_marginal12(key1_, key1_);
+ Matrix cov2 = joint_marginal12(key2_, key2_);
+ Matrix cov12 = joint_marginal12(key1_, key2_);
- return gtsam::GaussianFactor::shared_ptr(
- new gtsam::JacobianFactor(key1_, A1, key2_, A2, b, gtsam::noiseModel::Unit::Create(b.size())));
- }
+ updateNoiseModels_givenCovs(values, cov1, cov2, cov12);
+ }
+ /* ************************************************************************* */
+ void updateNoiseModels_givenCovs(const gtsam::Values& values,
+ const Matrix& cov1, const Matrix& cov2, const Matrix& cov12) {
+ /* Update model_inlier_ and model_outlier_ to account for uncertainty in robot trajectories
+ * (note these are given in the E step, where indicator probabilities are calculated).
+ *
+ * Principle: R += [H1 H2] * joint_cov12 * [H1 H2]', where H1, H2 are Jacobians of the
+ * unwhitened error w.r.t. states, and R is the measurement covariance (inlier or outlier modes).
+ *
+ * TODO: improve efficiency (info form)
+ */
- /* ************************************************************************* */
- gtsam::Vector whitenedError(const gtsam::Values& x,
- boost::optional&> H = boost::none) const {
+ const T& p1 = values.at(key1_);
+ const T& p2 = values.at(key2_);
- bool debug = true;
+ Matrix H1, H2;
+ p1.between(p2, H1, H2); // h(x)
- const T& p1 = x.at(key1_);
- const T& p2 = x.at(key2_);
+ Matrix H;
+ H.resize(H1.rows(), H1.rows() + H2.rows());
+ H << H1, H2; // H = [H1 H2]
- Matrix H1, H2;
+ Matrix joint_cov;
+ joint_cov.resize(cov1.rows() + cov2.rows(), cov1.cols() + cov2.cols());
+ joint_cov << cov1, cov12, cov12.transpose(), cov2;
- T hx = p1.between(p2, H1, H2); // h(x)
- // manifold equivalent of h(x)-z -> log(z,h(x))
+ Matrix cov_state = H * joint_cov * H.transpose();
- Vector err = measured_.localCoordinates(hx);
+ // model_inlier_->print("before:");
- // Calculate indicator probabilities (inlier and outlier)
- Vector p_inlier_outlier = calcIndicatorProb(x);
- double p_inlier = p_inlier_outlier[0];
- double p_outlier = p_inlier_outlier[1];
+ // update inlier and outlier noise models
+ Matrix covRinlier =
+ (model_inlier_->R().transpose() * model_inlier_->R()).inverse();
+ model_inlier_ = gtsam::noiseModel::Gaussian::Covariance(
+ covRinlier + cov_state);
- Vector err_wh_inlier = model_inlier_->whiten(err);
- Vector err_wh_outlier = model_outlier_->whiten(err);
+ Matrix covRoutlier =
+ (model_outlier_->R().transpose() * model_outlier_->R()).inverse();
+ model_outlier_ = gtsam::noiseModel::Gaussian::Covariance(
+ covRoutlier + cov_state);
- Matrix invCov_inlier = model_inlier_->R().transpose() * model_inlier_->R();
- Matrix invCov_outlier = model_outlier_->R().transpose() * model_outlier_->R();
+ // model_inlier_->print("after:");
+ // std::cout<<"covRinlier + cov_state: "<Whiten(H1);
- Matrix H1_outlier = sqrt(p_outlier)*model_outlier_->Whiten(H1);
- Matrix H1_aug = gtsam::stack(2, &H1_inlier, &H1_outlier);
+ virtual size_t dim() const {
+ return model_inlier_->R().rows() + model_inlier_->R().cols();
+ }
- Matrix H2_inlier = sqrt(p_inlier)*model_inlier_->Whiten(H2);
- Matrix H2_outlier = sqrt(p_outlier)*model_outlier_->Whiten(H2);
- Matrix H2_aug = gtsam::stack(2, &H2_inlier, &H2_outlier);
+private:
- (*H)[0].resize(H1_aug.rows(),H1_aug.cols());
- (*H)[1].resize(H2_aug.rows(),H2_aug.cols());
+ /** Serialization function */
+ friend class boost::serialization::access;
+ template
+ void serialize(ARCHIVE & ar, const unsigned int version) {
+ ar
+ & boost::serialization::make_nvp("NonlinearFactor",
+ boost::serialization::base_object(*this));
+ ar & BOOST_SERIALIZATION_NVP(measured_);
+ }
+};
+// \class BetweenFactorEM
- (*H)[0] = H1_aug;
- (*H)[1] = H2_aug;
- }
-
- if (debug){
- // std::cout<<"unwhitened error: "<whiten(err);
- Vector err_wh_outlier = model_outlier_->whiten(err);
-
- Matrix invCov_inlier = model_inlier_->R().transpose() * model_inlier_->R();
- Matrix invCov_outlier = model_outlier_->R().transpose() * model_outlier_->R();
-
- double p_inlier = prior_inlier_ * std::sqrt(invCov_inlier.determinant()) * exp( -0.5 * err_wh_inlier.dot(err_wh_inlier) );
- double p_outlier = prior_outlier_ * std::sqrt(invCov_outlier.determinant()) * exp( -0.5 * err_wh_outlier.dot(err_wh_outlier) );
-
- if (debug){
- std::cout<<"in calcIndicatorProb. err_unwh: "<(key1_);
- const T& p2 = x.at(key2_);
-
- Matrix H1, H2;
-
- T hx = p1.between(p2, H1, H2); // h(x)
-
- return measured_.localCoordinates(hx);
- }
-
- /* ************************************************************************* */
- void set_flag_bump_up_near_zero_probs(bool flag) {
- flag_bump_up_near_zero_probs_ = flag;
- }
-
- /* ************************************************************************* */
- bool get_flag_bump_up_near_zero_probs() const {
- return flag_bump_up_near_zero_probs_;
- }
-
- /* ************************************************************************* */
- SharedGaussian get_model_inlier() const {
- return model_inlier_;
- }
-
- /* ************************************************************************* */
- SharedGaussian get_model_outlier() const {
- return model_outlier_;
- }
-
- /* ************************************************************************* */
- Matrix get_model_inlier_cov() const {
- return (model_inlier_->R().transpose()*model_inlier_->R()).inverse();
- }
-
- /* ************************************************************************* */
- Matrix get_model_outlier_cov() const {
- return (model_outlier_->R().transpose()*model_outlier_->R()).inverse();
- }
-
- /* ************************************************************************* */
- void updateNoiseModels(const gtsam::Values& values, const gtsam::NonlinearFactorGraph& graph){
- /* Update model_inlier_ and model_outlier_ to account for uncertainty in robot trajectories
- * (note these are given in the E step, where indicator probabilities are calculated).
- *
- * Principle: R += [H1 H2] * joint_cov12 * [H1 H2]', where H1, H2 are Jacobians of the
- * unwhitened error w.r.t. states, and R is the measurement covariance (inlier or outlier modes).
- *
- * TODO: improve efficiency (info form)
- */
-
- // get joint covariance of the involved states
- std::vector Keys;
- Keys.push_back(key1_);
- Keys.push_back(key2_);
- Marginals marginals( graph, values, Marginals::QR );
- JointMarginal joint_marginal12 = marginals.jointMarginalCovariance(Keys);
- Matrix cov1 = joint_marginal12(key1_, key1_);
- Matrix cov2 = joint_marginal12(key2_, key2_);
- Matrix cov12 = joint_marginal12(key1_, key2_);
-
- updateNoiseModels_givenCovs(values, cov1, cov2, cov12);
- }
-
- /* ************************************************************************* */
- void updateNoiseModels_givenCovs(const gtsam::Values& values, const Matrix& cov1, const Matrix& cov2, const Matrix& cov12){
- /* Update model_inlier_ and model_outlier_ to account for uncertainty in robot trajectories
- * (note these are given in the E step, where indicator probabilities are calculated).
- *
- * Principle: R += [H1 H2] * joint_cov12 * [H1 H2]', where H1, H2 are Jacobians of the
- * unwhitened error w.r.t. states, and R is the measurement covariance (inlier or outlier modes).
- *
- * TODO: improve efficiency (info form)
- */
-
- const T& p1 = values.at(key1_);
- const T& p2 = values.at(key2_);
-
- Matrix H1, H2;
- T hx = p1.between(p2, H1, H2); // h(x)
-
- Matrix H;
- H.resize(H1.rows(), H1.rows()+H2.rows());
- H << H1, H2; // H = [H1 H2]
-
- Matrix joint_cov;
- joint_cov.resize(cov1.rows()+cov2.rows(), cov1.cols()+cov2.cols());
- joint_cov << cov1, cov12,
- cov12.transpose(), cov2;
-
- Matrix cov_state = H*joint_cov*H.transpose();
-
- // model_inlier_->print("before:");
-
- // update inlier and outlier noise models
- Matrix covRinlier = (model_inlier_->R().transpose()*model_inlier_->R()).inverse();
- model_inlier_ = gtsam::noiseModel::Gaussian::Covariance(covRinlier + cov_state);
-
- Matrix covRoutlier = (model_outlier_->R().transpose()*model_outlier_->R()).inverse();
- model_outlier_ = gtsam::noiseModel::Gaussian::Covariance(covRoutlier + cov_state);
-
- // model_inlier_->print("after:");
- // std::cout<<"covRinlier + cov_state: "<