Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
85 changes: 34 additions & 51 deletions examples/AbcEquivariantFilterExample.cpp
Original file line number Diff line number Diff line change
@@ -1,32 +1,31 @@
/**
* @file AbcEquivariantFilterExample.cpp
* @brief Demonstration of the full Attitude-Bias-Calibration Equivariant Filter
* @brief Demonstration of the Attitude-Bias-Calibration Equivariant Filter
*
* This demo shows the Equivariant Filter (EqF) for attitude estimation
* with both gyroscope bias and sensor extrinsic calibration, based on the
* paper: "Overcoming Bias: Equivariant Filter Design for Biased Attitude
* Estimation with Online Calibration" by Fornasier et al.
*
* This example uses the simplified AbcEquivariantFilter class which
* provides a clean interface for predict/update operations without requiring
* manual computation of Jacobian matrices and innovation functions.
*
* @author Darshan Rajasekaran
* @author Jennifer Oum
* @author Rohan Bansal
* @author Frank Dellaert
* @date 2025
*/
#include <gtsam/navigation/EquivariantFilter.h>
#include <gtsam/slam/dataset.h>
#include <gtsam_unstable/geometry/ABC.h>
#include <gtsam_unstable/geometry/ABCEquivariantFilter.h>

// Use namespace for convenience
using namespace gtsam;
constexpr size_t n = 1; // Number of calibration states
using M = abc::State<n>;
using G = abc::Group<n>;
using Symmetry = abc::Symmetry<n>;
using EqFilter = gtsam::EquivariantFilter<M, Symmetry>;
using Lift = abc::Lift<n>;
using InputOrbit = typename abc::InputAction<n>::Orbit;
using Innovation = abc::Innovation<n>;
using AbcFilter = abc::AbcEquivariantFilter<n>;

/// Measurement struct
struct Measurement {
Expand Down Expand Up @@ -72,7 +71,7 @@ std::vector<Data> loadDataFromCSV(const std::string& filename, int startRow = 0,
int maxRows = -1, int downsample = 1);

/// Process data with EqF and print summary results
void processDataWithEqF(EqFilter& filter, const std::vector<Data>& data_list,
void processDataWithEqF(AbcFilter& filter, const std::vector<Data>& data_list,
int printInterval = 10);

//========================================================================
Expand Down Expand Up @@ -235,7 +234,7 @@ std::vector<Data> loadDataFromCSV(const std::string& filename, int startRow,
}

/// Takes in the data and runs an EqF on it and reports the results
void processDataWithEqF(EqFilter& filter, const std::vector<Data>& data_list,
void processDataWithEqF(AbcFilter& filter, const std::vector<Data>& data_list,
int printInterval) {
if (data_list.empty()) {
std::cerr << "No data to process" << std::endl;
Expand Down Expand Up @@ -267,18 +266,7 @@ void processDataWithEqF(EqFilter& filter, const std::vector<Data>& data_list,

for (size_t i = 0; i < data_list.size(); i++) {
const Data& data = data_list[i];
Matrix Q = abc::inputProcessNoise<n>(data.inputCovariance);
// Propagate filter with current input and time step
Vector6 u = abc::toInputVector(data.omega);
Lift lift_u(u);
InputOrbit psi_u(u);

// Use Explicit Matrices API
G X_hat = filter.groupEstimate();
Matrix A = abc::stateMatrixA<n>(psi_u, X_hat);
Matrix B = abc::inputMatrixB<n>(X_hat);
Matrix Qc = B * Q * B.transpose(); // continuous-time manifold covariance
filter.predictWithJacobian<2>(lift_u, A, Qc, data.dt);
filter.predict(data.omega, data.inputCovariance, data.dt);

// Process all measurements
for (const auto& measurement : data.measurements) {
Expand All @@ -294,30 +282,26 @@ void processDataWithEqF(EqFilter& filter, const std::vector<Data>& data_list,
}

try {
Innovation innovation(measurement.y, measurement.d,
measurement.cal_idx);
// Use Explicit Matrices API
X_hat = filter.groupEstimate();
const Matrix3 D = abc::outputMatrixD<n>(X_hat, measurement.cal_idx);
const Matrix3 R = D * measurement.R * D.transpose();
filter.update<Vector3>(innovation, Z_3x1, R);

filter.update(measurement.y, measurement.d, measurement.R,
measurement.cal_idx);
validMeasurements++;
} catch (const std::exception& e) {
std::cerr << "Error updating at t=" << data.t << ": " << e.what()
<< std::endl;
}
}

// Get current state estimate
M estimate = filter.state();
// Get current estimates using accessor methods
Rot3 att_est = filter.attitude();
Vector3 bias_est = filter.bias();
Rot3 cal_est = filter.calibration(0);

// Calculate errors
Vector3 att_error = Rot3::Logmap(data.xi.R.between(estimate.R));
Vector3 bias_error = estimate.b - data.xi.b;
Vector3 att_error = Rot3::Logmap(data.xi.R.between(att_est));
Vector3 bias_error = bias_est - data.xi.b;
Vector3 cal_error = Z_3x1;
if (!data.xi.S.empty() && !estimate.S.empty()) {
cal_error = Rot3::Logmap(data.xi.S[0].between(estimate.S[0]));
if (!data.xi.S.empty()) {
cal_error = Rot3::Logmap(data.xi.S[0].between(cal_est));
}

// Store errors
Expand Down Expand Up @@ -353,14 +337,16 @@ void processDataWithEqF(EqFilter& filter, const std::vector<Data>& data_list,

// Calculate final errors from last data point
const Data& final_data = data_list.back();
M final_estimate = filter.state();
Rot3 final_att_est = filter.attitude();
Vector3 final_bias_est = filter.bias();
Rot3 final_cal_est = filter.calibration(0);

Vector3 final_att_error =
Rot3::Logmap(final_data.xi.R.between(final_estimate.R));
Vector3 final_bias_error = final_estimate.b - final_data.xi.b;
Rot3::Logmap(final_data.xi.R.between(final_att_est));
Vector3 final_bias_error = final_bias_est - final_data.xi.b;
Vector3 final_cal_error = Z_3x1;
if (!final_data.xi.S.empty() && !final_estimate.S.empty()) {
final_cal_error =
Rot3::Logmap(final_data.xi.S[0].between(final_estimate.S[0]));
if (!final_data.xi.S.empty()) {
final_cal_error = Rot3::Logmap(final_data.xi.S[0].between(final_cal_est));
}

// Print summary statistics
Expand Down Expand Up @@ -388,16 +374,15 @@ void processDataWithEqF(EqFilter& filter, const std::vector<Data>& data_list,
// Print a brief comparison of final estimate vs ground truth
std::cout << "\n-- Final State vs Ground Truth --" << std::endl;
std::cout << "Attitude (RPY) - Estimate: "
<< (final_estimate.R.rpy() * RAD_TO_DEG).transpose()
<< (final_att_est.rpy() * RAD_TO_DEG).transpose()
<< "° | Truth: " << (final_data.xi.R.rpy() * RAD_TO_DEG).transpose()
<< "°" << std::endl;
std::cout << "Bias - Estimate: " << final_estimate.b.transpose()
std::cout << "Bias - Estimate: " << final_bias_est.transpose()
<< " | Truth: " << final_data.xi.b.transpose() << std::endl;

if (!final_estimate.S.empty() && !final_data.xi.S.empty()) {
if (!final_data.xi.S.empty()) {
std::cout << "Calibration (RPY) - Estimate: "
<< (final_estimate.S[0].rpy() * RAD_TO_DEG).transpose()
<< "° | Truth: "
<< (final_cal_est.rpy() * RAD_TO_DEG).transpose() << "° | Truth: "
<< (final_data.xi.S[0].rpy() * RAD_TO_DEG).transpose() << "°"
<< std::endl;
}
Expand Down Expand Up @@ -456,10 +441,8 @@ int main(int argc, char* argv[]) {
initialSigma.diagonal().tail<3>() =
Vector3::Constant(0.1); // Calibration uncertainty

M initialState = M::identity();

// Create filter
EqFilter filter(initialState, initialSigma);
// Create filter with initial covariance (starts at identity state)
AbcFilter filter(initialSigma);

// Process data
processDataWithEqF(filter, data);
Expand Down
2 changes: 1 addition & 1 deletion gtsam/navigation/EquivariantFilter.h
Original file line number Diff line number Diff line change
Expand Up @@ -93,7 +93,7 @@ class EquivariantFilter : public ManifoldEKF<M> {
using Base::state;

/// errorCovariance that returns P_, on the equivariant filter error
Matrix errorCovariance() const { return this->P_; }
const typename Base::Covariance& errorCovariance() const { return this->P_; }

/// Covariance in the tangent space at the current state.
CovarianceM covariance() const {
Expand Down
21 changes: 21 additions & 0 deletions gtsam/navigation/navigation.i
Original file line number Diff line number Diff line change
Expand Up @@ -651,6 +651,27 @@ virtual class InvariantEKF : gtsam::LeftLinearEKF<G> {
void predict(const gtsam::Vector& u, double dt, gtsam::Matrix Q);
};

// ---------------------------------------------------------------------------
// ABC Equivariant Filter
#include <gtsam_unstable/geometry/ABCEquivariantFilter.h>
namespace abc {
template <N = {1, 2, 3}>
class AbcEquivariantFilter {
// Constructors
AbcEquivariantFilter();
AbcEquivariantFilter(gtsam::Matrix Sigma0);

// Predict and update methods
void predict(const gtsam::Vector3& omega, const gtsam::Matrix6& inputCovariance, double dt);
void update(const gtsam::Unit3& y, const gtsam::Unit3& d, const gtsam::Matrix3& R, int cal_idx);

// Accessors
gtsam::Rot3 attitude() const;
gtsam::Vector3 bias() const;
gtsam::Rot3 calibration(size_t i) const;
};
} // namespace abc

// Specialized NavState IMU EKF
#include <gtsam/navigation/NavStateImuEKF.h>
class NavStateImuEKF : gtsam::LeftLinearEKF<gtsam::NavState> {
Expand Down
136 changes: 136 additions & 0 deletions gtsam_unstable/geometry/ABCEquivariantFilter.h
Original file line number Diff line number Diff line change
@@ -0,0 +1,136 @@
/**
* @file ABCEquivariantFilter.h
* @brief Attitude-Bias-Calibration Equivariant Filter for state estimation
* @author Rohan Bansal
* @date 2026
*
* This file extends the EquivariantFilter class to provide a more user-friendly
* interface for the ABC Equivariant Filter. This class templates on the number
* of calibrated sensors N, abstracting away the details of the ABC system.
*/

#pragma once

#include <gtsam/base/Matrix.h>
#include <gtsam/base/Vector.h>
#include <gtsam/geometry/Rot3.h>
#include <gtsam/geometry/Unit3.h>
#include <gtsam/navigation/EquivariantFilter.h>
#include <gtsam_unstable/geometry/ABC.h>

namespace gtsam {
namespace abc {

/**
* @class AbcEquivariantFilter
* @brief Equivariant Filter for Attitude-Bias-Calibration (ABC) estimation
*
* This class implements an equivariant filter for estimating:
* - Attitude (rotation): The orientation of the body frame relative to a
* reference frame
* - Bias: Gyroscope bias correction vector (3D)
* - Calibration: N sensor calibration rotation matrices
*
* The filter uses the ABC Lie group structure and equivariant dynamics to
* provide consistent state estimation. It inherits from EquivariantFilter
* and provides a simplified interface for ABC-specific operations.
*
* @tparam N Number of calibrated sensors (typically 1 or more)
*/
template <size_t N>
class AbcEquivariantFilter
: public gtsam::EquivariantFilter<State<N>, Symmetry<N>> {
public:
using gtsam::EquivariantFilter<State<N>, Symmetry<N>>::update;

/**
* @brief Default constructor with identity initial covariance
*
* Initializes the filter with an identity state and identity covariance
* matrix of dimension (6 + 3*N) x (6 + 3*N), where 6 accounts for the
* attitude and bias parameters, and 3*N for N calibration parameters.
*/
AbcEquivariantFilter()
: AbcEquivariantFilter(Matrix::Identity(6 + 3 * N, 6 + 3 * N)) {}

/**
* @brief Construct filter with custom initial covariance
*
* Initializes the filter at the identity state with a specified initial
* covariance matrix on the manifold tangent space.
*
* @param Sigma0 Initial covariance matrix (6+3*N) x (6+3*N)
*/
explicit AbcEquivariantFilter(const Matrix& Sigma0)
: gtsam::EquivariantFilter<State<N>, Symmetry<N>>(State<N>::identity(),
Sigma0) {}

/**
* @brief Prediction step using gyroscope measurements
*
* Propagates the filter state forward in time using angular velocity
* measurements and process noise. This method uses the explicit Jacobian
* matrices (A and B) computed from the ABC dynamics.
*
* @param omega Angular velocity measurement (body frame, rad/s)
* @param inputCovariance Process noise covariance (6x6) for the input
* @param dt Time step (seconds)
*/
void predict(const Vector3& omega, const Matrix6& inputCovariance,
double dt) {
const Matrix Q = inputProcessNoise<N>(inputCovariance);
const Vector6 u = toInputVector(omega);
const Lift<N> lift_u(u);
const typename InputAction<N>::Orbit psi_u(u);

const Group<N> X_hat = this->groupEstimate();
const Matrix A = stateMatrixA<N>(psi_u, X_hat);
const Matrix B = inputMatrixB<N>(X_hat);
const Matrix Qc = B * Q * B.transpose();

this->template predictWithJacobian<2>(lift_u, A, Qc, dt);
}

/**
* @brief Measurement update using a direction observation
*
* Corrects the filter estimate using a direction measurement from a
* calibrated sensor. The measurement model assumes the sensor observes
* a known reference direction in the body frame.
*
* @param y Measured direction (unit vector in body frame)
* @param d Reference direction (unit vector in reference frame)
* @param R Measurement noise covariance (3x3)
* @param cal_idx Index of the calibrated sensor (0 to N-1), or -1 for
* uncalibrated measurements
*/
void update(const Unit3& y, const Unit3& d, const Matrix3& R, int cal_idx) {
const Innovation<N> innovation(y, d, cal_idx);
const Group<N> X_hat = this->groupEstimate();
const Matrix3 D = outputMatrixD<N>(X_hat, cal_idx);
const Matrix3 R_adjusted = D * R * D.transpose();
this->template update<Vector3>(innovation, Z_3x1, R_adjusted);
}

/**
* @brief Get current attitude estimate
* @return Current rotation estimate (body frame to reference frame)
*/
Rot3 attitude() const { return this->state().R; }

/**
* @brief Get current gyroscope bias estimate
* @return Current bias vector (rad/s)
*/
Vector3 bias() const { return this->state().b; }

/**
* @brief Get calibration estimate for a specific sensor
* @param i Sensor index (0 to N-1)
* @return Calibration rotation for sensor i
*/
Rot3 calibration(size_t i) const { return this->state().S[i]; }
};

} // namespace abc
} // namespace gtsam
Loading
Loading