Skip to content
Closed
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
126 changes: 126 additions & 0 deletions gtsam/geometry/SOT3.cpp
Original file line number Diff line number Diff line change
@@ -0,0 +1,126 @@
/* ----------------------------------------------------------------------------

* 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 SOT3.cpp
* @brief The scaled orthogonal transforms SOT(3) = SO(3) x R>0
* @author Rohan Bansal
*/

#include <gtsam/geometry/SOT3.h>

#include <cmath>
#include <iostream>

namespace gtsam {

void SOT3::print(const std::string& s) const {
if (!s.empty()) std::cout << s << "\n";
std::cout << "SOT3:\n";
R_.print(" R: ");
std::cout << " c: " << c_ << "\n";
}

bool SOT3::equals(const SOT3& other, double tol) const {
return R_.equals(other.R_, tol) && std::abs(c_ - other.c_) < tol;
}

// MatrixLieGroup interface

SOT3::MatrixNN SOT3::matrix() const {
MatrixNN M = MatrixNN::Zero();
M.topLeftCorner<3, 3>() = R_.matrix();
M(3, 3) = c_;
return M;
}

SOT3::MatrixNN SOT3::Hat(const TangentVector& xi) {
// (Omega, s) -> [[Omega^x, 0], [0, s]]
MatrixNN X = MatrixNN::Zero();
X.topLeftCorner<3, 3>() = SO3::Hat(Vector3(xi.head<3>()));
X(3, 3) = xi(3);
return X;
}

SOT3::TangentVector SOT3::Vee(const MatrixNN& X) {
// [[Omega^x, 0], [0, s]] -> (Omega, s)
TangentVector xi;
xi.head<3>() = SO3::Vee(Matrix3(X.topLeftCorner<3, 3>()));
xi(3) = X(3, 3);
return xi;
}

//******************************************************************************
// Lie Group

SOT3 SOT3::Expmap(const TangentVector& xi, ChartJacobian H) {
// SOT(3) = SO(3) x R>0 is direct product, so:
// exp(Omega, s) = (SO3::Expmap(Omega), exp(s))
const Vector3 Omega = xi.head<3>();
const double s = xi(3);
const double c = std::exp(s);

if (H) {
// block-diagonal Jacobian where top-left 3x3: right Jacobian of
// SO3::Expmap at Omega, and bottom-right: identity, since
// Local(exp(s), exp(s + ds)) = ds
Matrix3 H_R;
const SO3 R = SO3::Expmap(Omega, H ? &H_R : nullptr);
H->setZero();
H->topLeftCorner<3, 3>() = H_R;
(*H)(3, 3) = 1.0;
return SOT3(R, c);
}

return SOT3(SO3::Expmap(Omega), c);
}

SOT3::TangentVector SOT3::Logmap(const SOT3& Q, ChartJacobian H) {
// log(R, c) = (SO3::Logmap(R), log(c))
const double s = std::log(Q.c_);

if (H) {
// block-diagonal Jacobian where top-left 3x3: derivative of
// SO3::Logmap at R, and bottom-right: identity, since
// Local(exp(s), exp(s + ds)) = ds
Matrix3 H_R;
const Vector3 Omega = SO3::Logmap(Q.R_, H ? &H_R : nullptr);
H->setZero();
H->topLeftCorner<3, 3>() = H_R;
(*H)(3, 3) = 1.0;
TangentVector xi;
xi.head<3>() = Omega;
xi(3) = s;
return xi;
}

TangentVector xi;
xi.head<3>() = SO3::Logmap(Q.R_);
xi(3) = s;
return xi;
}

SOT3::Jacobian SOT3::AdjointMap() const {
Jacobian Ad = Jacobian::Identity();
Ad.topLeftCorner<3, 3>() = R_.matrix();
return Ad;
}

SOT3 SOT3::ChartAtOrigin::Retract(const TangentVector& xi, ChartJacobian H) {
return SOT3::Expmap(xi, H);
}

SOT3::TangentVector SOT3::ChartAtOrigin::Local(const SOT3& Q,
ChartJacobian H) {
return SOT3::Logmap(Q, H);
}

} // namespace gtsam
153 changes: 153 additions & 0 deletions gtsam/geometry/SOT3.h
Original file line number Diff line number Diff line change
@@ -0,0 +1,153 @@
/* ----------------------------------------------------------------------------

* 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 SOT3.h
* @brief The scaled orthogonal transforms SOT(3) = SO(3) x R>0
*
* This file implements the SOT(3) group as defined by eqVIO (van Goor et al., https://arxiv.org/pdf/2205.01980)
*
* @author Rohan Bansal
* @date 2026
*/

#pragma once

#include <gtsam/base/MatrixLieGroup.h>
#include <gtsam/dllexport.h>
#include <gtsam/geometry/SO3.h>

#include <string>

namespace gtsam {

/**
* @class SOT3
* @brief Scaled orthogonal transforms: SOT(3) = SO(3) x R>0.
*
* Represents group elements Q = (R_Q, c_Q) with R_Q in SO(3) and c_Q > 0.
* The 4x4 matrix representation is block-diagonal: [[R, 0], [0, c]].
* The Lie algebra has dimension 4: tangent vector xi = (Omega, s) in R^4.
*/
class GTSAM_EXPORT SOT3 : public MatrixLieGroup<SOT3, 4, 4> {
public:
/// @name Type definitions
/// @{

static constexpr int dimension = 4;
using MatrixNN = Eigen::Matrix4d;
using LieAlgebra = MatrixNN;
using TangentVector = Eigen::Matrix<double, 4, 1>;
using ChartJacobian = OptionalJacobian<4, 4>;
using Jacobian = Eigen::Matrix<double, 4, 4>;

/// @}
/// @name Constructors
/// @{

/// Default constructor: identity element (R = I_3, c = 1).
SOT3() : R_(), c_(1.0) {}

/// Construct from rotation and positive scale.
SOT3(const SO3& R, double c) : R_(R), c_(c) {}

Comment on lines +59 to +61
Copy link

Copilot AI Mar 4, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

SOT3 is defined with c>0, but the public constructor currently accepts any c (including 0/negative). This allows creating invalid group elements that will later cause undefined results (e.g., inverse() divides by 0 and Logmap() takes log(c_)). Consider enforcing the invariant at construction time (e.g., throw std::invalid_argument if c <= 0) and/or guarding in inverse()/Logmap() to fail fast with a clear error.

Copilot uses AI. Check for mistakes.
/// Construct 4x4 [[R, 0], [0, c]].
explicit SOT3(const MatrixNN& M)
: R_(SO3(M.topLeftCorner<3, 3>())), c_(M(3, 3)) {}

Comment on lines +63 to +65
Copy link

Copilot AI Mar 4, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The 4x4-matrix constructor assumes M is exactly [[R,0],[0,c]] but does not validate the zero blocks or that c is positive. This can silently accept arbitrary matrices and produce an invalid SOT3 instance. Suggest adding structural checks (similar to Gal3::Gal3(const Matrix5&) throwing std::invalid_argument for malformed matrices) and rejecting M(3,3) <= 0.

Copilot uses AI. Check for mistakes.
/// @}
/// @name Testable
/// @{

void print(const std::string& s = "") const;
bool equals(const SOT3& other, double tol = 1e-9) const;

/// @}
/// @name Group
/// @{

/// Identity element.
static SOT3 Identity() { return SOT3(); }

/// Inverse element
SOT3 inverse() const { return SOT3(R_.inverse(), 1.0 / c_); }

/// Group multiplication: (R1,c1)*(R2,c2) = (R1*R2, c1*c2).
SOT3 operator*(const SOT3& other) const {
return SOT3(R_ * other.R_, c_ * other.c_);
Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

If this is true, is SOT3 then not just a product group? :-)
It could be that we need things in the EqVio filter that they are not provided by product group, but if not, it's just a typedef :-/

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Could still be the case that we need to extend it, in which case we need to convert the product group to a CRPT, as I did for ExtendedPose3

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I also noticed that it can be represented as a product group, but thought it might be beneficial to stay faithful to van Goor's SOT(3) implementation from his library just in case it is needed down the road during my port.

That being said, maybe I can club this implementation together with the VIOGroup implementation, at which point I will know for sure whether it needs to be extended?

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Well, how is Van Goor's implementation different from product group? Will read answer in the morning :-)

Copy link
Copy Markdown
Contributor Author

@rohan-bansal rohan-bansal Mar 4, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

You are correct that the group itself is mathematically a product group.

What I meant by staying faithful to LiePP was that it is not implemented as a composition of groups, but as a dedicated type that exposes wedge, vee, asMatrix, adjoint as all explicit 4x4 matrix operations.

In order to expose those same functions in a product group setting, correct me if I am wrong, but I assume we would have to make SOT(3) a wrapper around ProductLieGroup<SO3, double> and manually define them.

However, after tracing through the code, I am also noticing now that wedge and vee are only defined for the purposes of van Goor's library, and aren't used for the rest of the eqVIO implementation. Maybe in this case it would make more sense to just define it as a product group in GTSAM.

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

If we need those, 2 ways to go: CRTP+extend, or actually define ProductMatrixLieGroup - let’s chat offline. Try with ProductLieGroup for now.

}

/// @}
/// @name MatrixLieGroup interface
/// @{

/// Return 4x4 [[R, 0], [0, c]].
MatrixNN matrix() const;

/// Hat: (Omega, s) -> [[Omega^x, 0], [0, s]].
static MatrixNN Hat(const TangentVector& xi);

/// Vee: inverse of Hat.
static TangentVector Vee(const MatrixNN& X);

/// Exponential map at identity.
static SOT3 Expmap(const TangentVector& xi, ChartJacobian H = {});

/// Logarithmic map at identity.
static TangentVector Logmap(const SOT3& Q, ChartJacobian H = {});

/// Adjoint map Ad_Q: block-diagonal [[R, 0], [0, 1]].
Jacobian AdjointMap() const;

// Chart at origin uses Expmap/Logmap directly.
struct ChartAtOrigin {
static SOT3 Retract(const TangentVector& xi, ChartJacobian H = {});
static TangentVector Local(const SOT3& Q, ChartJacobian H = {});
};

using LieGroup<SOT3, 4>::inverse;

/// @}
/// @name Group action on R^3
/// @{

/// Apply Q to point p in R^3: Qp = c * R * p.
Vector3 operator*(const Vector3& p) const { return c_ * (R_.matrix() * p); }

/// Apply inverse of Q to point p: Q^{-1}p = (1/c) * R^T * p.
Vector3 applyInverse(const Vector3& p) const {
return (1.0 / c_) * (R_.matrix().transpose() * p);
}

/// @}
/// @name Accessors
/// @{

/// Return R_Q, the SO(3) rotation component.
const SO3& rotation() const { return R_; }

/// Return c_Q, the positive scale component.
double scalar() const { return c_; }

/// @}

private:
SO3 R_; ///< Rotation component R_Q in SO(3).
double c_; ///< Scale component c_Q > 0.
};
Comment on lines +142 to +145
Copy link

Copilot AI Mar 4, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Most geometry Lie-group types in this repo provide Boost-serialization support behind #if GTSAM_ENABLE_BOOST_SERIALIZATION (e.g., Pose2 has a serialize() method in gtsam/geometry/Pose2.h:335-343, Similarity3 in gtsam/geometry/Similarity3.h:222-231). SOT3 currently has no serialization path, which can limit interoperability with existing serialization utilities/tests. Consider adding a serialize() method (for R_ and c_) and, if appropriate, extending testSerializationGeometry.cpp to cover SOT3.

Copilot uses AI. Check for mistakes.

template <>
struct traits<SOT3> : public internal::MatrixLieGroup<SOT3, 4> {};

template <>
struct traits<const SOT3> : public internal::MatrixLieGroup<SOT3, 4> {};

} // namespace gtsam
Loading
Loading