Tweak vector & quaternion slerp

This commit is contained in:
Nikita Lisitsa 2026-04-14 21:04:56 +03:00
parent 2d9355f829
commit b3cb60371e
2 changed files with 26 additions and 34 deletions

View file

@ -258,28 +258,25 @@ namespace psemek::math
template <typename T>
quaternion<T> slerp(quaternion<T> const & q0, quaternion<T> const & q1, T const & t)
{
// threshold is chosen so that for abs(x) < threshold the second term in
// sin(x) Taylor series is less than the minimum value representable by T
static auto const threshold = std::pow(6 * std::numeric_limits<T>::min(), T{1}/T{3});
using std::sin;
using std::acos;
using std::abs;
auto const d = clamp(dot(normalized(q0.coords), normalized(q1.coords)), {T(-1), T(1)});
auto const omega = std::acos(std::abs(d));
auto const d = dot(q0.coords, q1.coords);
// prevent division by zero
if (d >= T{1})
return quaternion<T>{lerp(q0.coords, q1.coords, t)};
auto const omega = acos(abs(d));
// NB: the case of omega ~ pi is ambiguous and isn't handled in any special way
if (std::abs(omega) < threshold)
{
// prevent division by zero
return quaternion<T>{normalized(lerp(q0.coords, q1.coords, t))};
}
else
{
auto const s = std::sin(omega);
auto const w0 = std::sin((1 - t) * omega) / s;
auto const w1 = std::sin(t * omega) / s * ((d > 0) ? 1 : -1);
auto const s = sin(omega);
auto const w0 = sin((1 - t) * omega) / s;
auto const w1 = sin(t * omega) / s * ((d > T{0}) ? T{1} : -T{1});
return quaternion<T>{q0.coords * w0 + q1.coords * w1};
}
}
template <typename T>
quaternion<T> exp(quaternion<T> const & q)

View file

@ -473,31 +473,26 @@ namespace psemek::math
return v0 * (T(1) - t) + v1 * t;
}
// NB: v0 and v1 are assumed to be normalized
template <typename T, std::size_t N>
vector<T, N> slerp(vector<T, N> const & v0, vector<T, N> const & v1, T const & t)
{
using std::pow;
using std::abs;
using std::sin;
using std::acos;
// threshold is chosen so that for abs(x) < threshold the second term in
// sin(x) Taylor series is less than the minimum value representable by T
static auto const threshold = pow(T{6} * std::numeric_limits<T>::min(), T{1}/T{3});
auto const d = dot(v0, v1);
auto const omega = angle(normalized(v0), normalized(v1));
// prevent division by zero
if (d >= T{1})
return lerp(v0, v1, t);
auto const omega = acos(d);
// NB: the case of omega ~ pi is ambiguous and isn't handled in any special way
if (abs(omega) < threshold)
{
// prevent division by zero
return normalized(lerp(v0, v1, t));
}
else
{
auto const s = sin(omega);
return v0 * sin((1 - t) * omega) / s + v1 * sin(t * omega) / s;
}
auto const w0 = std::sin((1 - t) * omega) / s;
auto const w1 = (sin(t * omega) / s);
return w0 * v0 + w1 * v1;
}
// Return vector orthogonal to n