From e554ea30be39d332b331e55aafff7b2e655e669e Mon Sep 17 00:00:00 2001 From: Alexei Kotov Date: Wed, 23 Jul 2025 01:34:57 +0300 Subject: [PATCH] Implement quadratic/TCB interpolation for quaternions (#2379) --- components/misc/mathutil.hpp | 23 ++++++++++++ components/nif/nifkey.hpp | 60 +++++++++++++++++++++++++++++--- components/nifosg/controller.hpp | 11 +++++- 3 files changed, 89 insertions(+), 5 deletions(-) diff --git a/components/misc/mathutil.hpp b/components/misc/mathutil.hpp index e0c52eaf0c..0529240fd5 100644 --- a/components/misc/mathutil.hpp +++ b/components/misc/mathutil.hpp @@ -94,6 +94,29 @@ namespace Misc } return 1 << depth; } + + inline osg::Quat quatexp(const osg::Quat& q) + { + const float imagNorm2 = q.asVec3().length2(); + const float e = std::exp(q.w()); + if (imagNorm2 < 1e-6f) + return { 0.f, 0.f, 0.f, e }; + + const float imagNorm = std::sqrt(imagNorm2); + const float f = e * std::sin(imagNorm) / imagNorm; + return { q.x() * f, q.y() * f, q.z() * f, e * std::cos(imagNorm) }; + } + + inline osg::Quat quatlog(const osg::Quat& q) + { + const float imagNorm2 = q.asVec3().length2(); + const float norm = q.length(); + if (imagNorm2 < 1e-6f) + return { 0.f, 0.f, 0.f, std::log(norm) }; + + const float f = std::acos(q.w() / norm) / std::sqrt(imagNorm2); + return { q.x() * f, q.y() * f, q.z() * f, std::log(norm) }; + } } #endif diff --git a/components/nif/nifkey.hpp b/components/nif/nifkey.hpp index e32ef76d95..775b732a4f 100644 --- a/components/nif/nifkey.hpp +++ b/components/nif/nifkey.hpp @@ -6,6 +6,8 @@ #include #include +#include + #include "exception.hpp" #include "niffile.hpp" #include "nifstream.hpp" @@ -28,7 +30,7 @@ namespace Nif { T mValue; T mInTan; // Only for Quadratic interpolation, and never for QuaternionKeyList - T mOutTan; // Only for Quadratic interpolation, and never for QuaternionKeyList + T mOutTan; // For quaternions these are generated }; template @@ -107,6 +109,7 @@ namespace Nif readQuadratic(*nif, key); mKeys.emplace_back(time, key); } + generateQuadraticTangents(mKeys); } else if (mInterpolationType == InterpolationType_TCB) { @@ -155,6 +158,30 @@ namespace Nif static void readQuadratic(NIFStream& nif, KeyT& key) { readValue(nif, key); } + template + static void generateQuadraticTangents(std::vector>>& keys) + { + // These are predefined for all types except quaternions + } + + static void generateQuadraticTangents(std::vector>>& keys) + { + if (keys.size() <= 1) + return; + + for (std::size_t i = 0; i < keys.size(); ++i) + { + KeyT& curr = keys[i].second; + const std::size_t prevIndex = i == 0 ? 0 : i - 1; + const std::size_t nextIndex = i == keys.size() - 1 ? i : i + 1; + const osg::Quat inv = curr.mValue.inverse(); + const osg::Quat& prev = keys[prevIndex].second.mValue; + const osg::Quat& next = keys[nextIndex].second.mValue; + const osg::Quat len = Misc::quatlog(inv * prev) + Misc::quatlog(inv * next); + curr.mInTan = curr.mOutTan = curr.mValue * Misc::quatexp(len * -0.25f); + } + } + template static void generateTCBTangents(std::vector>& keys) { @@ -176,8 +203,8 @@ namespace Nif const float w = (1.f - curr.mTension) * (1.f - curr.mContinuity) * (1.f - curr.mBias); const U prevDelta = prev != nullptr ? curr.mValue - prev->mValue : next->mValue - curr.mValue; const U nextDelta = next != nullptr ? next->mValue - curr.mValue : curr.mValue - prev->mValue; - curr.mInTan = (prevDelta * x + nextDelta * y) * prevLen / (prevLen + nextLen); - curr.mOutTan = (prevDelta * z + nextDelta * w) * nextLen / (prevLen + nextLen); + curr.mInTan = (prevDelta * x + nextDelta * y) * (prevLen / (prevLen + nextLen)); + curr.mOutTan = (prevDelta * z + nextDelta * w) * (nextLen / (prevLen + nextLen)); } } @@ -188,7 +215,32 @@ namespace Nif static void generateTCBTangents(std::vector>& keys) { - // TODO: implement TCB interpolation for quaternions + if (keys.size() <= 1) + return; + + for (std::size_t i = 0; i < keys.size(); ++i) + { + TCBKey& curr = keys[i]; + const TCBKey* prev = (i == 0) ? nullptr : &keys[i - 1]; + const TCBKey* next = (i == keys.size() - 1) ? nullptr : &keys[i + 1]; + const float prevLen = prev != nullptr && next != nullptr ? curr.mTime - prev->mTime : 1.f; + const float nextLen = prev != nullptr && next != nullptr ? next->mTime - curr.mTime : 1.f; + if (prevLen + nextLen == 0.f) + continue; + const float x = (1.f - curr.mTension) * (1.f - curr.mContinuity) * (1.f + curr.mBias); + const float y = (1.f - curr.mTension) * (1.f + curr.mContinuity) * (1.f - curr.mBias); + const float z = (1.f - curr.mTension) * (1.f + curr.mContinuity) * (1.f + curr.mBias); + const float w = (1.f - curr.mTension) * (1.f - curr.mContinuity) * (1.f - curr.mBias); + const osg::Quat inv = curr.mValue.inverse(); + const osg::Quat prevDelta = prev != nullptr ? Misc::quatlog(prev->mValue.inverse() * curr.mValue) + : Misc::quatlog(inv * next->mValue); + const osg::Quat nextDelta = next != nullptr ? Misc::quatlog(inv * next->mValue) + : Misc::quatlog(prev->mValue.inverse() * curr.mValue); + const osg::Quat t1 = (prevDelta * x + nextDelta * y) * (nextLen / (prevLen + nextLen)); + const osg::Quat t2 = (prevDelta * z + nextDelta * w) * (prevLen / (prevLen + nextLen)); + curr.mInTan = curr.mValue * Misc::quatexp((prevDelta - t1) * 0.5f); + curr.mOutTan = curr.mValue * Misc::quatexp((t2 - nextDelta) * 0.5f); + } } }; diff --git a/components/nifosg/controller.hpp b/components/nifosg/controller.hpp index cceec279b2..872bcc04c4 100644 --- a/components/nifosg/controller.hpp +++ b/components/nifosg/controller.hpp @@ -165,7 +165,16 @@ namespace NifOsg { case Nif::InterpolationType_Constant: return fraction > 0.5f ? b.mValue : a.mValue; - // TODO: Implement Quadratic and TBC interpolation + case Nif::InterpolationType_Quadratic: + case Nif::InterpolationType_TCB: + { + // Spherical quadrangle interpolation + osg::Quat from, to, result; + from.slerp(fraction, a.mValue, b.mValue); + to.slerp(fraction, a.mOutTan, b.mInTan); + result.slerp(2.f * fraction * (1.f - fraction), from, to); + return result; + } default: { osg::Quat result;