Representation and interpolation of rotations in complex numbers.

main.cpp

// C++14 or later.
#include <iostream>
#include <cstdlib>
#include <complex>

namespace
{

template <class T>
constexpr T pi = static_cast<T>(3.14159265358979323846);    

template <class T>
T deg_to_rad(T angle)
{
    return angle / static_cast<T>(180) * pi<T>;
}

template <class T>
T rad_to_deg(T angle)
{
    return angle / pi<T> * static_cast<T>(180);
}

template<typename T>
T dot(const std::complex<T>& z0, const std::complex<T>& z1)
{
    return z0.real() * z1.real() + z0.imag() * z1.imag();
}
    
}   // namespace

// standard functions.

template<typename T>
std::complex<T>
slerp(const std::complex<T>& p, const std::complex<T>& q, T h)
{
    constexpr auto tolerance = static_cast<T>(0.01);
    const auto cosine = dot(p, q);

    T coeff0, coeff1;
    
    if((static_cast<T>(1) - std::abs(cosine)) > tolerance)
    {
        const auto angle = std::acos(cosine);
        const auto cosec = static_cast<T>(1) / std::sin(angle);
        coeff0 = std::sin((static_cast<T>(1) - h) * angle) * cosec;
        coeff1 = std::sin(h * angle) * cosec;
    }
    else
    {   // // Switch to linear interpolation.
        coeff0 = static_cast<T>(1) - h;
        coeff1 = h;
    }
    return coeff0 * p + coeff1 * q;
}

template<typename T>
std::complex<T>
slerp_longest(const std::complex<T>& p, const std::complex<T>& q, T h)
{
    constexpr auto tolerance = static_cast<T>(0.01);
    constexpr auto two_pi = static_cast<T>(2) * pi<T>; 
    const auto cosine = dot(p, q);

    T coeff0, coeff1;
    
    if((static_cast<T>(1) - std::abs(cosine)) > tolerance)
    {
        auto phi = std::arg(p);
        if(phi < static_cast<T>(0))
            phi += two_pi;

        auto psi = std::arg(p);
        if(psi < static_cast<T>(0))
            psi += two_pi;

        const auto det = p.real() * q.imag() - p.imag() * q.real();

        T angle;
        if(psi - phi < static_cast<T>(0))
        {
            angle = (det < static_cast<T>(0))? std::acos(cosine) : two_pi - std::acos(cosine);
        }
        else
        {
            angle = (det < static_cast<T>(0))? two_pi - std::acos(cosine) : std::acos(cosine);
        }

        const auto cosec = static_cast<T>(1) / std::sin(angle);
        coeff0 = std::sin((static_cast<T>(1) - h) * angle) * cosec;
        coeff1 = std::sin(h * angle) * cosec;
    }
    else
    {   // // Switch to linear interpolation.
        coeff0 = static_cast<T>(1) - h;
        coeff1 = h;
    }
    return coeff0 * p + coeff1 * q;
}

// spinor functions.

template<typename T>
std::complex<T>
make_spinor(T angle)
{
    const auto half_angle = angle / static_cast<T>(2);
    return std::complex<T>(std::cos(half_angle), std::sin(half_angle));
}

template<typename T>
std::complex<T>
rotate_by_spinor(const std::complex<T>& p, const std::complex<T>& v)
{
      return p * v * p;
}

template<typename T>
std::complex<T>
spinor_slerp_shortest(const std::complex<T>& p, const std::complex<T>& q, T h)
{
    const auto cosine = dot(p, q);
    return slerp(p, cosine < static_cast<T>(0)? -q : q, h);
}

int main()
{
    const auto v = std::complex<double>(1.0, 0.0);

    std::cout << "=== Test for rotation." << std::endl;
    {
        std::cout << "[Angle] : " << "[Standard], [Spinor]" << std::endl;
        for(int i = 0; i < 9; i++)
        {
            auto angle = static_cast<double>(i * 45);
            auto ret0 = std::polar(1.0, deg_to_rad(angle)) * v;
            auto ret1 = rotate_by_spinor(make_spinor(deg_to_rad(angle)), v);
            std::cout << "Deg " << angle << ": " << ret0 << ", " << ret1 << std::endl;
        }
    }
    std::cout << std::endl;
    std::cout << "=== Test for SLERP." << std::endl;
    {
        const auto angle0 = deg_to_rad(0.0);
        const auto angle1 = deg_to_rad(270.0);

        std::cout << "[Standard]" << std::endl;
        {
            const auto p = std::polar(1.0, angle0);
            const auto q = std::polar(1.0, angle1);
            
            std::cout << "p: " << p << std::endl;
            std::cout << "q: " << q << std::endl;

            std::cout << "[h]: [default], [longest]" << std::endl;
            std::cout << "0.0: " << slerp(p, q, 0.0) << ", " << slerp_longest(p, q, 0.0) << std::endl;
            std::cout << "0.5: " << slerp(p, q, 0.5) << ", " << slerp_longest(p, q, 0.5) << std::endl;
            std::cout << "1.0: " << slerp(p, q, 1.0) << ", " << slerp_longest(p, q, 1.0) << std::endl;
        }
        std::cout << std::endl;
        std::cout << "[Spinor]" << std::endl;
        {
            const auto p = make_spinor(angle0);
            const auto q = make_spinor(angle1);

            std::cout << "p: " << p << std::endl;
            std::cout << "q: " << q << std::endl;

            std::cout << "[h]: [default], [shortest]" << std::endl;
            std::cout << "0.0: " << slerp(p, q, 0.0) << ", " << spinor_slerp_shortest(p, q, 0.0) << std::endl;
            std::cout << "0.5: " << slerp(p, q, 0.5) << ", " << spinor_slerp_shortest(p, q, 0.5) << std::endl;
            std::cout << "1.0: " << slerp(p, q, 1.0) << ", " << spinor_slerp_shortest(p, q, 1.0) << std::endl;
        }
    }
    return 1;
}