Rotor Swing-Twist Decomposition - Multivector Element Access Clarity

[CaseySanchez]

Hello again Dr. Fernandes,

I have yet another example for you. In the snippet below I perform a swing-twist decomposition of a rotor about an axis. I then constrain these two resulting swing and twist rotors by a 45-degree angle. I reconstitute the two rotors into a single rotor, yielding a fully-constrained rotor.

#include <gatl/ga3e.hpp>

using namespace ga3e;

template <typename T, typename = std::enable_if_t<std::is_floating_point<T>::value>>
using Rotor = decltype(full_kvector_t<T, 3, 0>() + full_kvector_t<T, 3, 2>());

template <typename T, typename = std::enable_if_t<std::is_floating_point<T>::value>>
using Translator = decltype(full_kvector_t<T, 3, 1>());

template <typename T, typename = std::enable_if_t<std::is_floating_point<T>::value>>
using Motor = decltype(Rotor<T>() * Translator<T>());

template <typename T, typename = std::enable_if_t<std::is_floating_point<T>::value>>
void decompose(Rotor<T> const &rotor, Vector<T> const &vector, Rotor<T> &swing, Rotor<T> &twist)
{
    Vector<T> const rotor_vector = *std::next(std::cbegin(rotor.values()), 1) * e1 + *std::next(std::cbegin(rotor.values()), 2) * e2 + *std::next(std::cbegin(rotor.values()), 3) * e3;
    Vector<T> const rotor_project = project(rotor_vector, vector);

    twist = *std::cbegin(rotor.values()) + *std::cbegin(rotor_project.values()) * (e1 ^ e2) + *std::next(std::cbegin(rotor_project.values()), 1) * (e1 ^ e3) + *std::next(std::cbegin(rotor_project.values()), 2) * (e2 ^ e3);

    T const twist_norm = std::sqrt(*std::cbegin((twist * ~twist).values()));

    twist = (static_cast<T>(1.0) / twist_norm) * twist;
    swing = rotor * inv(twist) * inv(rotor);
}

template <typename T, typename = std::enable_if_t<std::is_floating_point<T>::value>>
Rotor<T> constrain(T const &theta, Rotor<T> const &rotor)
{
    T const magnitude = std::sin(theta * static_cast<T>(0.5));
    T const sqr_magnitude = magnitude * magnitude;

    Translator<T> const vector = *std::next(std::cbegin(rotor.values()), 1) * e1 + *std::next(std::cbegin(rotor.values()), 2) * e2 + *std::next(std::cbegin(rotor.values()), 3) * e3;

    T const vector_sqr_magnitude = *std::cbegin((vector * ~vector).values());

    if (vector_sqr_magnitude > sqr_magnitude) {
        Translator<T> const vector_normal = (static_cast<T>(1.0) / std::sqrt(vector_sqr_magnitude)) * vector;

        T const sign = std::signbit(*std::cbegin(rotor.values())) ? static_cast<T>(-1.0) : static_cast<T>(+1.0);

        return std::sqrt(static_cast<T>(1.0) - sqr_magnitude) * sign + *std::cbegin(vector_normal.values()) * (e1 ^ e2) + *std::next(std::cbegin(vector_normal.values()), 1) * (e1 ^ e3) + *std::next(std::cbegin(vector_normal.values()), 2) * (e2 ^ e3);
    }

    return rotor;
}

int main(int argc, char **argv)
{
    rotor_t<float> const rotor = std::cos(-3.141592f * 0.5f * 0.5f) + std::sin(-3.141592f * 0.5f * 0.5f) * (e1 ^ e2) + 0.0f * (e1 ^ e3) + 0.0f * (e2 ^ e3);
    translator_t<float> const translator = std::sqrt(3.0f) / 3.0f * e1 + std::sqrt(3.0f) / 3.0f * e2 + std::sqrt(3.0f) / 3.0f * e3;

    rotor_t<float> swing;
    rotor_t<float> twist;

    decompose<float>(rotor, translator, swing, twist);

    constrain<float>(3.141592f * 0.25f, swing);
    constrain<float>(3.141592f * 0.25f, twist);

    rotor_t<float> const rotor_constrained = swing * twist;

    std::cout << "Rotor: ";

    for (float const &value : rotor.values()) {
        std::cout << value << " ";
    }

    std::cout << std::endl;

    std::cout << "Rotor constrained: ";

    for (float const &value : rotor_constrained.values()) {
        std::cout << value << " ";
    }

    std::cout << std::endl;

    return 0;
}

As you can see, the solution in its current state is less-than-elegant. Is there any better way of accessing the individual elements of the multivectors aside from manipulating the iterator with std::cbegin and std::next?

Regards, Casey

[CaseySanchez]

I suppose I should re-write my question thusly: how would you yourself approach this problem @laffernandes? The way I wrote it clearly relies upon explicitly having a scalar part and the three bivector parts of a rotor, which is not representative of all possible rotor combinations as you have mentioned (e.g. the identity rotor being only a scalar with a value of 1, or a 90-degree rotation of e1 ^ e2 being sqrt(2) * 0.5 + sqrt(2) * 0.5 * (e1 ^ e2)). This goes against your philosophy of GATL; I would much prefer to abide by your philosophy.

EDIT: Would it perhaps be most wise to accept any arbitrary multivector type and thereafter perform a trivial_copy to the Rotor<T> in order to extract the relevant components?

EDIT2: As described in EDIT, does this seem better to you? Is there room for more improvement in your opinion?

#pragma once

#include <gatl/ga3e.hpp>

using namespace ga3e;

template <typename T, typename = std::enable_if_t<std::is_floating_point<T>::value>>
using Rotor = decltype(full_kvector_t<T, 3, 0>() + full_kvector_t<T, 3, 2>());

template <typename T, typename = std::enable_if_t<std::is_floating_point<T>::value>>
using Translator = decltype(full_kvector_t<T, 3, 1>());

template <typename T, typename = std::enable_if_t<std::is_floating_point<T>::value>>
using Vector = decltype(full_kvector_t<T, 3, 1>());

template <typename T, typename = std::enable_if_t<std::is_floating_point<T>::value>>
using Motor = decltype(Rotor<T>() * Translator<T>());

template <typename T, typename RotorExpression, typename VectorExpression, typename = std::enable_if_t<std::is_floating_point<T>::value && is_clifford_expression<RotorExpression>::value && is_clifford_expression<VectorExpression>::value>>
void decompose(RotorExpression const &rotor_expression, VectorExpression const &vector_expression, Rotor<T> &swing, Rotor<T> &twist)
{
    Rotor<T> rotor;
    Vector<T> vector;

    checked_trivial_copy(rotor_expression, rotor);
    checked_trivial_copy(vector_expression, vector);

    Vector<T> const rotor_vector = *std::next(std::cbegin(rotor.values()), 1) * e1 + *std::next(std::cbegin(rotor.values()), 2) * e2 + *std::next(std::cbegin(rotor.values()), 3) * e3;
    Vector<T> const rotor_project = project(rotor_vector, vector);

    twist = *std::cbegin(rotor.values()) + *std::cbegin(rotor_project.values()) * (e1 ^ e2) + *std::next(std::cbegin(rotor_project.values()), 1) * (e1 ^ e3) + *std::next(std::cbegin(rotor_project.values()), 2) * (e2 ^ e3);

    T const twist_norm = std::sqrt(*std::cbegin((twist * ~twist).values()));

    twist = (static_cast<T>(1.0) / twist_norm) * twist;
    swing = rotor * inv(twist) * inv(rotor);
}

template <typename T, typename RotorExpression, typename = std::enable_if_t<std::is_floating_point<T>::value && is_clifford_expression<RotorExpression>::value>>
Rotor<T> constrain(T const &theta, RotorExpression const &rotor_expression)
{
    Rotor<T> rotor;

    checked_trivial_copy(rotor_expression, rotor);

    T const norm = std::sin(theta * static_cast<T>(0.5));
    T const sqr_norm = norm * norm;

    Translator<T> const vector = *std::next(std::cbegin(rotor.values()), 1) * e1 + *std::next(std::cbegin(rotor.values()), 2) * e2 + *std::next(std::cbegin(rotor.values()), 3) * e3;

    T const vector_sqr_norm = *std::cbegin((vector * ~vector).values());

    if (vector_sqr_norm > sqr_norm) {
        Translator<T> const vector_normal = (static_cast<T>(1.0) / std::sqrt(vector_sqr_norm)) * vector;

        T const sign = std::signbit(*std::cbegin(rotor.values())) ? static_cast<T>(-1.0) : static_cast<T>(+1.0);

        return std::sqrt(static_cast<T>(1.0) - sqr_norm) * sign + *std::cbegin(vector_normal.values()) * (e1 ^ e2) + *std::next(std::cbegin(vector_normal.values()), 1) * (e1 ^ e3) + *std::next(std::cbegin(vector_normal.values()), 2) * (e2 ^ e3);
    }
    else {
        return rotor;
    }
}

int main(int argc, char **argv)
{
    auto const rotor = std::cos(-3.141592f * 0.5f * 0.5f) + std::sin(-3.141592f * 0.5f * 0.5f) * (e1 ^ e2);
    auto const translator = std::sqrt(3.0f) / 3.0f * e1 + std::sqrt(3.0f) / 3.0f * e2 + std::sqrt(3.0f) / 3.0f * e3;

    Rotor<float> swing;
    Rotor<float> twist;

    decompose<float>(rotor, translator, swing, twist);

    constrain<float>(3.141592f * 0.25f, swing);
    constrain<float>(3.141592f * 0.25f, twist);

    Rotor<float> const rotor_constrained = swing * twist;

    std::cout << "Rotor: ";

    for (float const &value : rotor.values()) {
        std::cout << value << " ";
    }

    std::cout << std::endl;

    std::cout << "Rotor constrained: ";

    for (float const &value : rotor_constrained.values()) {
        std::cout << value << " ";
    }

    std::cout << std::endl;

    return 0;
}

EDIT3: Would it be possible to implement user-defined literals such as multivector_e12 or even something more complicated as multivector_(e1 ^ e2) (I'm not even sure if this particular example is possible) to extract these coefficients?

[laffernandes]

The first code you presented can be rewritten this way:

#include <gatl/ga3e.hpp>

using namespace ga3e;

template<typename T1, typename RotorExpression, typename T2, typename TranslatorExpression>
auto decompose(clifford_expression<T1, RotorExpression> const &rotor, clifford_expression<T2, TranslatorExpression> const &translator)
{
    T1 x = sp(rotor, -e1^e2); // Write mathematical expressions to access particular coefficients of a multivector without using iterators. It is for free due to how GATL handles those expressions in compile time. Here, `sp` is the scalar product of blades, and the minus sign in e1^e2 prevents the change of orientation of the result (this is typical in Geometric Algebra... you will get used to it).
    T1 y = sp(rotor, -e1^e3);
    T1 z = sp(rotor, -e2^e3);
    auto rotor_vector = vector(x, y, z);

    auto rotor_project = project(rotor_vector, translator);

    auto twist = unit(sp(rotor, c<1>) + sp(rotor_project, e1) * (e1 ^ e2) + sp(rotor_project, e2) * (e1 ^ e3) + sp(rotor_project, e3) * (e2 ^ e3)); // In GATL, c<integral_value> means a compile-time defined constant integral value. In this case, the value 1. The `unit` function return the unit versor under reverse norm. 
    auto swing = apply_rotor(rotor, inv(twist));

    return std::make_tuple(swing, twist);
}

template<typename T1, typename T2, typename RotorExpression>
auto constrain(T1 const theta, clifford_expression<T2, RotorExpression> const &rotor)
{
    T1 magnitude = std::sin(theta * 0.5f);
    T1 sqr_magnitude = magnitude * magnitude;

    T2 x = sp(rotor, -e1^e2); // Write mathematical expressions to access particular coefficients of a multivector without using iterators. It is for free due to how GATL handles those expressions in compile time. Here, `sp` is the scalar product of blades, and the minus sign in e1^e2 prevents the change of orientation of the result (this is typical in Geometric Algebra... you will get used to it).
    T2 y = sp(rotor, -e1^e3);
    T2 z = sp(rotor, -e2^e3);
    auto rotor_vector = vector(x, y, z);

    T2 rotor_vector_sqr_magnitude = rnorm_sqr(rotor_vector);

    if (rotor_vector_sqr_magnitude > sqr_magnitude) {
        auto vector_normal = rotor_vector / sqrt(rotor_vector_sqr_magnitude);
        int sign = sp(rotor, c<1>) < 0.0f ? -1 : 1;
        return std::sqrt(1.0 - sqr_magnitude) * sign + sp(vector_normal, e1) * (e1 ^ e2) + sp(vector_normal, e2) * (e1 ^ e3) + sp(vector_normal, e3) * (e2 ^ e3);
    }

    return rotor;
}

int main(int argc, char **argv)
{
    auto rotor = exp((-3.141592f * 0.5f * 0.5f) * (e1 ^ e2));
    auto translator = unit(vector(1.0f, 1.0f, 1.0f));

    auto [swing, twist] = decompose(rotor, translator); // The structured binding in C++17 is convenient! It looks like Python.

    constrain(3.141592f * 0.25f, swing);
    constrain(3.141592f * 0.25f, twist);

    auto rotor_constrained = swing * twist;

    std::cout << "Rotor: " << rotor << std::endl;
    std::cout << "Rotor constrained: " << rotor_constrained << std::endl;

    return EXIT_SUCCESS;
}

What I find strange about the mathematical expressions in your implementation (and in my re-implementation above following your ideas) is that you have to look at the multivector components. More general solutions in Geometrical Algebra do not usually do this. Although it is common for some papers to cherry-pick coefficients, elegant solutions do not. For example, instead of selecting each component of the bivector portion of the rotor, wouldn't it be better to use take_grade to take the two-dimensional portion of the rotor (the rotation plane) and then take its dual (using the dual function) to have the vector normal to that plane? This would work on any number of dimensions, whereas the implemented solution only works in 3-D.

As selecting specific coefficients is not a good practice in Geometric Algebra, I will not include user-defined literal to avoid encouraging this practice.

[CaseySanchez]

This is excellent! Thank you for sharing.

Regarding the concepts of swing-twist decomposition of quaternions/rotors, I must say that I am a fraction of the mathematician that you are, so perhaps you would better understand it than I do. All that I know is that it works!

I had originally seen it described here: https://www.euclideanspace.com/maths/geometry/rotations/for/decomposition/

I'm going to close this issue, hopefully it will serve as useful reference for users in the future.

Thank you Dr. Fernandes.