Skip to content

Commit

Permalink
Quternion draft. Untested, must be bugged as hell
Browse files Browse the repository at this point in the history
  • Loading branch information
lhog committed Dec 3, 2024
1 parent 119e917 commit c9f367f
Show file tree
Hide file tree
Showing 7 changed files with 394 additions and 5 deletions.
1 change: 1 addition & 0 deletions rts/System/CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -34,6 +34,7 @@ make_global_var(sources_engine_System_common
"${CMAKE_CURRENT_SOURCE_DIR}/Math/SpringDampers.cpp"
"${CMAKE_CURRENT_SOURCE_DIR}/Math/NURBS.cpp"
"${CMAKE_CURRENT_SOURCE_DIR}/Matrix44f.cpp"
"${CMAKE_CURRENT_SOURCE_DIR}/Quaternion.cpp"
"${CMAKE_CURRENT_SOURCE_DIR}/Misc/RectangleOverlapHandler.cpp"
"${CMAKE_CURRENT_SOURCE_DIR}/Misc/SpringTime.cpp"
"${CMAKE_CURRENT_SOURCE_DIR}/Object.cpp"
Expand Down
28 changes: 27 additions & 1 deletion rts/System/Matrix44f.cpp
Original file line number Diff line number Diff line change
@@ -1,6 +1,7 @@
/* This file is part of the Spring engine (GPL v2 or later), see LICENSE.html */

#include "System/Matrix44f.h"
#include "System/Quaternion.h"
#include "System/SpringMath.h"
#ifndef UNIT_TEST
#include "Rendering/GlobalRendering.h"
Expand Down Expand Up @@ -70,6 +71,22 @@ bool CMatrix44f::IsIdentity() const
return (*this) == IDENTITY;
}

float CMatrix44f::Det3() const
{
// triple product == D
return col[0].dot(col[1].cross(col[2]));
}

float CMatrix44f::Det4() const
{
return
md[0][3] * md[1][2] * md[2][1] * md[3][0] - md[0][2] * md[1][3] * md[2][1] * md[3][0] - md[0][3] * md[1][1] * md[2][2] * md[3][0] + md[0][1] * md[1][3] * md[2][2] * md[3][0] +
md[0][2] * md[1][1] * md[2][3] * md[3][0] - md[0][1] * md[1][2] * md[2][3] * md[3][0] - md[0][3] * md[1][2] * md[2][0] * md[3][1] + md[0][2] * md[1][3] * md[2][0] * md[3][1] +
md[0][3] * md[1][0] * md[2][2] * md[3][1] - md[0][0] * md[1][3] * md[2][2] * md[3][1] - md[0][2] * md[1][0] * md[2][3] * md[3][1] + md[0][0] * md[1][2] * md[2][3] * md[3][1] +
md[0][3] * md[1][1] * md[2][0] * md[3][2] - md[0][1] * md[1][3] * md[2][0] * md[3][2] - md[0][3] * md[1][0] * md[2][1] * md[3][2] + md[0][0] * md[1][3] * md[2][1] * md[3][2] +
md[0][1] * md[1][0] * md[2][3] * md[3][2] - md[0][0] * md[1][1] * md[2][3] * md[3][2] - md[0][2] * md[1][1] * md[2][0] * md[3][3] + md[0][1] * md[1][2] * md[2][0] * md[3][3] +
md[0][2] * md[1][0] * md[2][1] * md[3][3] - md[0][0] * md[1][2] * md[2][1] * md[3][3] - md[0][1] * md[1][0] * md[2][2] * md[3][3] + md[0][0] * md[1][1] * md[2][2] * md[3][3];
}

CMatrix44f& CMatrix44f::RotateX(float angle)
{
Expand Down Expand Up @@ -272,7 +289,7 @@ CMatrix44f& CMatrix44f::RotateEulerZYX(const float3 angles)
// [ 0 0 sz 0]
// [ 0 0 0 1]
//
CMatrix44f& CMatrix44f::Scale(const float3 scales)
CMatrix44f& CMatrix44f::Scale(const float3& scales)
{
m[ 0] *= scales.x;
m[ 1] *= scales.x;
Expand All @@ -291,6 +308,15 @@ CMatrix44f& CMatrix44f::Scale(const float3 scales)
return *this;
}

CMatrix44f& CMatrix44f::FromTQS(const float3& pos, const CQuaternion& quat, const float3& scale)
{
*this = quat.ToRotMatrix();
this->Scale(scale);
this->Translate(pos);

return *this;
}

CMatrix44f& CMatrix44f::Translate(const float x, const float y, const float z)
{
m[12] += (x*m[0] + y*m[4] + z*m[ 8]);
Expand Down
11 changes: 8 additions & 3 deletions rts/System/Matrix44f.h
Original file line number Diff line number Diff line change
Expand Up @@ -9,6 +9,7 @@
#include "System/float3.h"
#include "System/float4.h"

class CQuaternion;
class CMatrix44f
{
public:
Expand All @@ -30,6 +31,10 @@ class CMatrix44f

bool IsOrthoNormal() const;
bool IsIdentity() const;
bool IsRotMatrix() const { return IsOrthoNormal() && math::fabs(1.0f - Det4()) <= float3::cmp_eps(); }
bool IsRotOrRotTranMatrix() const { return IsOrthoNormal() && math::fabs(1.0f - Det3()) <= float3::cmp_eps(); }
float Det3() const;
float Det4() const;

CMatrix44f& LoadIdentity() { return (*this = CMatrix44f()); }

Expand All @@ -43,9 +48,9 @@ class CMatrix44f
CMatrix44f& RotateEulerZXY(const float3 angles); // executes Rotate{Z,X,Y}
CMatrix44f& RotateEulerZYX(const float3 angles); // executes Rotate{Z,Y,X}
CMatrix44f& Translate(const float x, const float y, const float z);
CMatrix44f& Translate(const float3 pos) { return Translate(pos.x, pos.y, pos.z); }
CMatrix44f& Scale(const float3 scales);
CMatrix44f& Scale(float scaleX, float scaleY, float scaleZ) { return Scale(float3{ scaleX, scaleY, scaleZ }); }
CMatrix44f& Translate(const float3& pos) { return Translate(pos.x, pos.y, pos.z); }
CMatrix44f& Scale(const float3& scales);
CMatrix44f& FromTQS(const float3& pos, const CQuaternion& quat, const float3& scale);

void SetPos(const float3 pos) { m[12] = pos.x; m[13] = pos.y; m[14] = pos.z; }
void SetX (const float3 dir) { m[ 0] = dir.x; m[ 1] = dir.y; m[ 2] = dir.z; }
Expand Down
275 changes: 275 additions & 0 deletions rts/System/Quaternion.cpp
Original file line number Diff line number Diff line change
@@ -0,0 +1,275 @@
#include "Quaternion.h"
/* This file is part of the Spring engine (GPL v2 or later), see LICENSE.html */

#include "System/Quaternion.h"
#include "System/SpringMath.h"


//contains some code from
// https://github.com/ilmola/gml/blob/master/include/gml/quaternion.hpp
// https://github.com/ilmola/gml/blob/master/include/gml/mat.hpp
// https://github.com/g-truc/glm/blob/master/glm/ext/quaternion_common.inl

CR_BIND(CQuaternion, )
CR_REG_METADATA(CQuaternion, CR_MEMBER(q))

/// <summary>
/// Quaternion from Euler angles
/// </summary>
CQuaternion CQuaternion::MakeFrom(const float3& euler)
{
const float a1 = 0.5f * euler.x;
const float a2 = 0.5f * euler.y;
const float a3 = 0.5f * euler.z;

const float sx = math::sin(a1);
const float cx = math::cos(a1);
const float sy = math::sin(a2);
const float cy = math::cos(a2);
const float sz = math::sin(a3);
const float cz = math::cos(a3);

return CQuaternion(
cx * cy * cz + sx * sy * sz,
float3(
cy * cz * sx - cx * sy * sz,
cx * cz * sy + cy * sx * sz,
-cz * sx * sy + cx * cy * sz
)
);
}

/// <summary>
/// Quaternion from rotation angle and axis
/// </summary>
CQuaternion CQuaternion::MakeFrom(float angle, const float3& axis)
{
assert(axis.Normalized());

const float a = 0.5f * angle;
return CQuaternion(math::cos(a), axis * math::sin(a));
}

/// <summary>
/// Quaternion to rotate from v1 to v2
/// </summary>
CQuaternion CQuaternion::MakeFrom(const float3& v1, const float3& v2)
{
float dp = v1.dot(v2);
if unlikely(v1.same(v2)) {
return CQuaternion(0.0f, v1);
}
else if unlikely(v1.same(-v2)) {
float3 v;
if (v1.x > -float3::cmp_eps() && v1.x < float3::cmp_eps()) // if x ~= 0
v = { 1.0f, 0.0f, 0.0f };
else if (v1.y > -float3::cmp_eps() && v1.y < float3::cmp_eps()) // if y ~= 0
v = { 0.0f, 1.0f, 0.0f };
else // if z ~= 0
v = { 0.0f, 0.0f, 1.0f };

return CQuaternion(math::HALFPI, v);
}
else {
float3 u1 = v1; u1.Normalize();
float3 u2 = v1; u2.Normalize();

float3 v = u1.cross(u2); // compute rotation axis
float angle = math::acosf(u1.dot(u2)); // rotation angle
return CQuaternion(angle * 0.5f, v); // half angle
}
}

/// <summary>
/// Quaternion from a rotation matrix
/// Should only be called on R or T * R matrices
/// </summary>
CQuaternion CQuaternion::MakeFrom(const CMatrix44f& mat)
{
assert(mat.IsRotOrRotTranMatrix());
const float trace = mat.md[0][0] + mat.md[1][1] + mat.md[2][2];

if (trace > 0.0f) {
const float s = 0.5f * math::isqrt(trace + 1.0f);

return CQuaternion(
0.25f / s,
s * float3(
(mat.md[1][2] - mat.md[2][1]),
(mat.md[2][0] - mat.md[0][2]),
(mat.md[0][1] - mat.md[1][0])
)
);
}
else if (mat.md[0][0] > mat.md[1][1] && mat.md[0][0] > mat.md[2][2]) {
const float s = 2.0f * math::sqrt(1.0f + mat.md[0][0] - mat.md[1][1] - mat.md[2][2]);

return CQuaternion(
(mat.md[1][2] - mat.md[2][1]) / s,
float3(
0.25f * s,
(mat.md[1][0] + mat.md[0][1]) / s,
(mat.md[2][0] + mat.md[0][2]) / s
)
);
}
else if (mat.md[1][1] > mat.md[2][2]) {
const float s = 2.0f * math::sqrt(1.0f + mat.md[1][1] - mat.md[0][0] - mat.md[2][2]);

return CQuaternion(
(mat.md[2][0] - mat.md[0][2]) / s,
float3(
(mat.md[1][0] + mat.md[0][1]) / s,
0.25f * s,
(mat.md[2][1] + mat.md[1][2]) / s
)
);
}
else {
const float s = 2.0f * math::sqrt(1.0f + mat.md[2][2] - mat.md[0][0] - mat.md[1][1]);

return CQuaternion(
(mat.md[0][1] - mat.md[1][0]) / s,
float3(
(mat.md[2][0] + mat.md[0][2]) / s,
(mat.md[2][1] + mat.md[1][2]) / s,
0.25f * s
)
);
}
}

/// <summary>
/// Decompose a transformation matrix into translate, rotation (Quaternion), scale components
/// </summary>
std::tuple<float3, CQuaternion, float3> CQuaternion::DecomposeIntoTRS(const CMatrix44f& mat)
{
CMatrix44f tmpMat = mat;
float4& t0 = tmpMat.col[0];
float4& t1 = tmpMat.col[1];
float4& t2 = tmpMat.col[2];

const float4& c0 = mat.col[0];
const float4& c1 = mat.col[1];
const float4& c2 = mat.col[2];

//triple product == determinant
const float d = c0.dot(c1.cross(c2));

const float s = Sign(d);

float3 scaling {s * c0.Length(), c1.Length(), c2.Length()};


t0[0] /= scaling[0];
t1[1] /= scaling[1];
t2[2] /= scaling[2];

return std::make_tuple(
float3(mat.col[3]), //translate
CQuaternion::MakeFrom(tmpMat), //rotate (quat)
scaling //scale
);
}


CQuaternion& CQuaternion::Normalize()
{
const float sqn = SqNorm();
if unlikely(sqn < float3::nrm_eps())
return *this;

q *= math::isqrt(sqn);

return *this;
}

/// <summary>
/// Find axis and angle equivalent rotation from a quaternion
/// </summary>
float4 CQuaternion::ToAxisAndAngle() const
{
return float4(
static_cast<float3>(q) * math::isqrt(std::max(0.0f, 1.0f - q.w * q.w)),
2.0f * math::acos(std::clamp(q.w, -1.0f, 1.0f))
);
}

/// <summary>
/// Converts a quaternion to rotational matrix
/// </summary>
CMatrix44f CQuaternion::ToRotMatrix() const
{
const float qxx = q.x * q.x;
const float qyy = q.y * q.y;
const float qzz = q.z * q.z;
const float qxz = q.x * q.z;
const float qxy = q.x * q.y;
const float qyz = q.y * q.z;
const float qwx = q.w * q.x;
const float qwy = q.w * q.y;
const float qwz = q.w * q.z;

return CMatrix44f(
1.0f - 2.0f * (qyy + qzz), 2.0f * (qxy - qwz) , 2.0f * (qxz + qwy) , 0.0f,
2.0f * (qxy + qwz) , 1.0f - 2.0f * (qxx + qzz), 2.0f * (qyz - qwx) , 0.0f,
2.0f * (qxz - qwy) , 2.0f * (qyz + qwx) , 1.0f - 2.0f * (qxx + qyy), 0.0f,
0.0f , 0.0f , 0.0f , 1.0f
);
}

CQuaternion& CQuaternion::Inverse()
{
const float sqn = SqNorm();
if unlikely(sqn < float3::nrm_eps())
return *this;

Conjugate() / SqNorm();
return *this;
};

CQuaternion CQuaternion::operator*(const CQuaternion& rhs) const
{
const float3& limag = q;
const float3& rimag = rhs.q;

return CQuaternion(
q.w * rhs.q.w - limag.dot(rimag),
q.w * rimag + rhs.q.w * limag + limag.cross(rimag)
);
}

float CQuaternion::SqNorm() const {
return q.SqLength();
}

CQuaternion CQuaternion::Lerp(const CQuaternion& q1, const CQuaternion& q2, const float a) {
return q1 * (1.0f - a) + (q2 * a);
}

CQuaternion CQuaternion::SLerp(const CQuaternion& q1, const CQuaternion& q2_, const float a) {
if (a == 0.0f)
return q1;
else if (a == 1.0f)
return q2_;

float cosTheta = q1.dot(q2_);

const float s = Sign(cosTheta);

CQuaternion q2 = q2_ * s;
cosTheta *= s;

if unlikely(cosTheta > 1.0f - float3::cmp_eps()) {
// Linear interpolation
return CQuaternion(
mix(q1.q, q2.q, a)
);
} else {
// Essential Mathematics, page 467
const float angle = math::acos(cosTheta);
return CQuaternion(
(math::sin((1.0f - a) * angle) * q1.q + math::sin(a * angle) * q2.q) / math::sin(angle)
);
}
}
Loading

0 comments on commit c9f367f

Please sign in to comment.