148 lines
		
	
	
		
			5.2 KiB
		
	
	
	
		
			C++
		
	
	
			
		
		
	
	
			148 lines
		
	
	
		
			5.2 KiB
		
	
	
	
		
			C++
		
	
	
/* ----------------------------------------------------------------------------
 | 
						|
 | 
						|
 * GTSAM Copyright 2010, Georgia Tech Research Corporation,
 | 
						|
 * Atlanta, Georgia 30332-0415
 | 
						|
 * All Rights Reserved
 | 
						|
 * Authors: Frank Dellaert, et al. (see THANKS for the full author list)
 | 
						|
 | 
						|
 * See LICENSE for the license information
 | 
						|
 | 
						|
 * -------------------------------------------------------------------------- */
 | 
						|
 | 
						|
/**
 | 
						|
 * @file ImuFactorExample2
 | 
						|
 * @brief Test example for using GTSAM ImuFactor and ImuCombinedFactor with ISAM2.
 | 
						|
 * @author Robert Truax
 | 
						|
 */
 | 
						|
 | 
						|
#include <gtsam/geometry/PinholeCamera.h>
 | 
						|
#include <gtsam/geometry/Cal3_S2.h>
 | 
						|
#include <gtsam/inference/Symbol.h>
 | 
						|
#include <gtsam/navigation/ImuBias.h>
 | 
						|
#include <gtsam/navigation/ImuFactor.h>
 | 
						|
#include <gtsam/navigation/Scenario.h>
 | 
						|
#include <gtsam/nonlinear/ISAM2.h>
 | 
						|
#include <gtsam/slam/BetweenFactor.h>
 | 
						|
 | 
						|
#include <vector>
 | 
						|
 | 
						|
using namespace std;
 | 
						|
using namespace gtsam;
 | 
						|
 | 
						|
// Shorthand for velocity and pose variables
 | 
						|
using symbol_shorthand::V;
 | 
						|
using symbol_shorthand::X;
 | 
						|
 | 
						|
const double kGravity = 9.81;
 | 
						|
 | 
						|
/* ************************************************************************* */
 | 
						|
int main(int argc, char* argv[]) {
 | 
						|
  auto params = PreintegrationParams::MakeSharedU(kGravity);
 | 
						|
  params->setAccelerometerCovariance(I_3x3 * 0.1);
 | 
						|
  params->setGyroscopeCovariance(I_3x3 * 0.1);
 | 
						|
  params->setIntegrationCovariance(I_3x3 * 0.1);
 | 
						|
  params->setUse2ndOrderCoriolis(false);
 | 
						|
  params->setOmegaCoriolis(Vector3(0, 0, 0));
 | 
						|
 | 
						|
  Pose3 delta(Rot3::Rodrigues(-0.1, 0.2, 0.25), Point3(0.05, -0.10, 0.20));
 | 
						|
 | 
						|
  // Start with a camera on x-axis looking at origin
 | 
						|
  double radius = 30;
 | 
						|
  const Point3 up(0, 0, 1), target(0, 0, 0);
 | 
						|
  const Point3 position(radius, 0, 0);
 | 
						|
  const auto camera = PinholeCamera<Cal3_S2>::Lookat(position, target, up);
 | 
						|
  const auto pose_0 = camera.pose();
 | 
						|
 | 
						|
  // Now, create a constant-twist scenario that makes the camera orbit the
 | 
						|
  // origin
 | 
						|
  double angular_velocity = M_PI,  // rad/sec
 | 
						|
      delta_t = 1.0 / 18;          // makes for 10 degrees per step
 | 
						|
  Vector3 angular_velocity_vector(0, -angular_velocity, 0);
 | 
						|
  Vector3 linear_velocity_vector(radius * angular_velocity, 0, 0);
 | 
						|
  auto scenario = ConstantTwistScenario(angular_velocity_vector,
 | 
						|
                                        linear_velocity_vector, pose_0);
 | 
						|
 | 
						|
  // Create a factor graph
 | 
						|
  NonlinearFactorGraph newgraph;
 | 
						|
 | 
						|
  // Create (incremental) ISAM2 solver
 | 
						|
  ISAM2 isam;
 | 
						|
 | 
						|
  // Create the initial estimate to the solution
 | 
						|
  // Intentionally initialize the variables off from the ground truth
 | 
						|
  Values initialEstimate, totalEstimate, result;
 | 
						|
 | 
						|
  // Add a prior on pose x0. This indirectly specifies where the origin is.
 | 
						|
  // 0.1 rad std on roll, pitch, yaw, 30cm std on x,y,z.
 | 
						|
  auto noise = noiseModel::Diagonal::Sigmas(
 | 
						|
      (Vector(6) << Vector3::Constant(0.1), Vector3::Constant(0.3)).finished());
 | 
						|
  newgraph.addPrior(X(0), pose_0, noise);
 | 
						|
 | 
						|
  // Add imu priors
 | 
						|
  Key biasKey = Symbol('b', 0);
 | 
						|
  auto biasnoise = noiseModel::Diagonal::Sigmas(Vector6::Constant(0.1));
 | 
						|
  newgraph.addPrior(biasKey, imuBias::ConstantBias(), biasnoise);
 | 
						|
  initialEstimate.insert(biasKey, imuBias::ConstantBias());
 | 
						|
  auto velnoise = noiseModel::Diagonal::Sigmas(Vector3(0.1, 0.1, 0.1));
 | 
						|
 | 
						|
  Vector n_velocity(3);
 | 
						|
  n_velocity << 0, angular_velocity * radius, 0;
 | 
						|
  newgraph.addPrior(V(0), n_velocity, velnoise);
 | 
						|
 | 
						|
  initialEstimate.insert(V(0), n_velocity);
 | 
						|
 | 
						|
  // IMU preintegrator
 | 
						|
  PreintegratedImuMeasurements accum(params);
 | 
						|
 | 
						|
  // Simulate poses and imu measurements, adding them to the factor graph
 | 
						|
  for (size_t i = 0; i < 36; ++i) {
 | 
						|
    double t = i * delta_t;
 | 
						|
    if (i == 0) {  // First time add two poses
 | 
						|
      auto pose_1 = scenario.pose(delta_t);
 | 
						|
      initialEstimate.insert(X(0), pose_0.compose(delta));
 | 
						|
      initialEstimate.insert(X(1), pose_1.compose(delta));
 | 
						|
    } else if (i >= 2) {  // Add more poses as necessary
 | 
						|
      auto pose_i = scenario.pose(t);
 | 
						|
      initialEstimate.insert(X(i), pose_i.compose(delta));
 | 
						|
    }
 | 
						|
 | 
						|
    if (i > 0) {
 | 
						|
      // Add Bias variables periodically
 | 
						|
      if (i % 5 == 0) {
 | 
						|
        biasKey++;
 | 
						|
        Symbol b1 = biasKey - 1;
 | 
						|
        Symbol b2 = biasKey;
 | 
						|
        Vector6 covvec;
 | 
						|
        covvec << 0.1, 0.1, 0.1, 0.1, 0.1, 0.1;
 | 
						|
        auto cov = noiseModel::Diagonal::Variances(covvec);
 | 
						|
        auto f = std::make_shared<BetweenFactor<imuBias::ConstantBias> >(
 | 
						|
            b1, b2, imuBias::ConstantBias(), cov);
 | 
						|
        newgraph.add(f);
 | 
						|
        initialEstimate.insert(biasKey, imuBias::ConstantBias());
 | 
						|
      }
 | 
						|
      // Predict acceleration and gyro measurements in (actual) body frame
 | 
						|
      Vector3 measuredAcc = scenario.acceleration_b(t) -
 | 
						|
                            scenario.rotation(t).transpose() * params->n_gravity;
 | 
						|
      Vector3 measuredOmega = scenario.omega_b(t);
 | 
						|
      accum.integrateMeasurement(measuredAcc, measuredOmega, delta_t);
 | 
						|
 | 
						|
      // Add Imu Factor
 | 
						|
      ImuFactor imufac(X(i - 1), V(i - 1), X(i), V(i), biasKey, accum);
 | 
						|
      newgraph.add(imufac);
 | 
						|
 | 
						|
      // insert new velocity, which is wrong
 | 
						|
      initialEstimate.insert(V(i), n_velocity);
 | 
						|
      accum.resetIntegration();
 | 
						|
    }
 | 
						|
 | 
						|
    // Incremental solution
 | 
						|
    isam.update(newgraph, initialEstimate);
 | 
						|
    result = isam.calculateEstimate();
 | 
						|
    newgraph = NonlinearFactorGraph();
 | 
						|
    initialEstimate.clear();
 | 
						|
  }
 | 
						|
  GTSAM_PRINT(result);
 | 
						|
  return 0;
 | 
						|
}
 | 
						|
/* ************************************************************************* */
 |