504 lines
16 KiB
C++
504 lines
16 KiB
C++
#pragma once
|
|
|
|
#include "Library.h"
|
|
|
|
#include <glm/gtc/epsilon.hpp>
|
|
|
|
namespace CesiumUtility {
|
|
|
|
/**
|
|
* @brief Mathematical constants and functions
|
|
*/
|
|
class CESIUMUTILITY_API Math final {
|
|
public:
|
|
/** @brief 0.1 */
|
|
static constexpr double Epsilon1 = 1e-1;
|
|
|
|
/** @brief 0.01 */
|
|
static constexpr double Epsilon2 = 1e-2;
|
|
|
|
/** @brief 0.001 */
|
|
static constexpr double Epsilon3 = 1e-3;
|
|
|
|
/** @brief 0.0001 */
|
|
static constexpr double Epsilon4 = 1e-4;
|
|
|
|
/** @brief 0.00001 */
|
|
static constexpr double Epsilon5 = 1e-5;
|
|
|
|
/** @brief 0.000001 */
|
|
static constexpr double Epsilon6 = 1e-6;
|
|
|
|
/** @brief 0.0000001 */
|
|
static constexpr double Epsilon7 = 1e-7;
|
|
|
|
/** @brief 0.00000001 */
|
|
static constexpr double Epsilon8 = 1e-8;
|
|
|
|
/** @brief 0.000000001 */
|
|
static constexpr double Epsilon9 = 1e-9;
|
|
|
|
/** @brief 0.0000000001 */
|
|
static constexpr double Epsilon10 = 1e-10;
|
|
|
|
/** @brief 0.00000000001 */
|
|
static constexpr double Epsilon11 = 1e-11;
|
|
|
|
/** @brief 0.000000000001 */
|
|
static constexpr double Epsilon12 = 1e-12;
|
|
|
|
/** @brief 0.0000000000001 */
|
|
static constexpr double Epsilon13 = 1e-13;
|
|
|
|
/** @brief 0.00000000000001 */
|
|
static constexpr double Epsilon14 = 1e-14;
|
|
|
|
/** @brief 0.000000000000001 */
|
|
static constexpr double Epsilon15 = 1e-15;
|
|
|
|
/** @brief 0.0000000000000001 */
|
|
static constexpr double Epsilon16 = 1e-16;
|
|
|
|
/** @brief 0.00000000000000001 */
|
|
static constexpr double Epsilon17 = 1e-17;
|
|
|
|
/** @brief 0.000000000000000001 */
|
|
static constexpr double Epsilon18 = 1e-18;
|
|
|
|
/** @brief 0.0000000000000000001 */
|
|
static constexpr double Epsilon19 = 1e-19;
|
|
|
|
/** @brief 0.00000000000000000001 */
|
|
static constexpr double Epsilon20 = 1e-20;
|
|
|
|
/** @brief 0.000000000000000000001 */
|
|
static constexpr double Epsilon21 = 1e-21;
|
|
|
|
/**
|
|
* @brief pi
|
|
*/
|
|
static constexpr double OnePi = 3.14159265358979323846;
|
|
|
|
/**
|
|
* @brief two times pi
|
|
*/
|
|
static constexpr double TwoPi = OnePi * 2.0;
|
|
|
|
/**
|
|
* @brief pi divded by two
|
|
*/
|
|
static constexpr double PiOverTwo = OnePi / 2.0;
|
|
|
|
/**
|
|
* @brief Converts a relative to an absolute epsilon, for the epsilon-equality
|
|
* check between two values.
|
|
*
|
|
* @tparam L The length type.
|
|
* @tparam T value value type.
|
|
* @tparam Q The GLM qualifier type.
|
|
*
|
|
* @param a The first value.
|
|
* @param b The second value.
|
|
* @param relativeEpsilon The relative epsilon.
|
|
* @return The absolute epsilon.
|
|
*/
|
|
template <glm::length_t L, typename T, glm::qualifier Q>
|
|
static constexpr glm::vec<L, T, Q> relativeEpsilonToAbsolute(
|
|
const glm::vec<L, T, Q>& a,
|
|
const glm::vec<L, T, Q>& b,
|
|
double relativeEpsilon) noexcept {
|
|
return relativeEpsilon * glm::max(glm::abs(a), glm::abs(b));
|
|
}
|
|
|
|
/**
|
|
* @brief Converts a relative to an absolute epsilon, for the epsilon-equality
|
|
* check between two values.
|
|
*
|
|
* @param a The first value.
|
|
* @param b The second value.
|
|
* @param relativeEpsilon The relative epsilon.
|
|
* @return The absolute epsilon.
|
|
*/
|
|
static constexpr double relativeEpsilonToAbsolute(
|
|
double a,
|
|
double b,
|
|
double relativeEpsilon) noexcept {
|
|
return relativeEpsilon * glm::max(glm::abs(a), glm::abs(b));
|
|
}
|
|
|
|
/**
|
|
* @brief Checks whether two values are equal up to a given relative epsilon.
|
|
*
|
|
* @tparam L The length type.
|
|
* @tparam T value value type.
|
|
* @tparam Q The GLM qualifier type.
|
|
*
|
|
* @param left The first value.
|
|
* @param right The second value.
|
|
* @param relativeEpsilon The relative epsilon.
|
|
* @return Whether the values are epsilon-equal
|
|
*/
|
|
template <glm::length_t L, typename T, glm::qualifier Q>
|
|
static bool constexpr equalsEpsilon(
|
|
const glm::vec<L, T, Q>& left,
|
|
const glm::vec<L, T, Q>& right,
|
|
double relativeEpsilon) noexcept {
|
|
return Math::equalsEpsilon(left, right, relativeEpsilon, relativeEpsilon);
|
|
}
|
|
|
|
/**
|
|
* @brief Checks whether two values are equal up to a given relative epsilon.
|
|
*
|
|
* @param left The first value.
|
|
* @param right The second value.
|
|
* @param relativeEpsilon The relative epsilon.
|
|
* @return Whether the values are epsilon-equal
|
|
*/
|
|
static constexpr bool
|
|
equalsEpsilon(double left, double right, double relativeEpsilon) noexcept {
|
|
return equalsEpsilon(left, right, relativeEpsilon, relativeEpsilon);
|
|
}
|
|
|
|
/**
|
|
* @brief Determines if two values are equal using an absolute or relative
|
|
* tolerance test.
|
|
*
|
|
* This is useful to avoid problems due to roundoff error when comparing
|
|
* floating-point values directly. The values are first compared using an
|
|
* absolute tolerance test. If that fails, a relative tolerance test is
|
|
* performed. Use this test if you are unsure of the magnitudes of left and
|
|
* right.
|
|
*
|
|
* @param left The first value to compare.
|
|
* @param right The other value to compare.
|
|
* @param relativeEpsilon The maximum inclusive delta between `left` and
|
|
* `right` for the relative tolerance test.
|
|
* @param absoluteEpsilon The maximum inclusive delta between `left` and
|
|
* `right` for the absolute tolerance test.
|
|
* @returns `true` if the values are equal within the epsilon; otherwise,
|
|
* `false`.
|
|
*
|
|
* @snippet TestMath.cpp equalsEpsilon
|
|
*/
|
|
static constexpr bool equalsEpsilon(
|
|
double left,
|
|
double right,
|
|
double relativeEpsilon,
|
|
double absoluteEpsilon) noexcept {
|
|
const double diff = glm::abs(left - right);
|
|
return diff <= absoluteEpsilon ||
|
|
diff <= relativeEpsilonToAbsolute(left, right, relativeEpsilon);
|
|
}
|
|
|
|
/**
|
|
* @brief Determines if two values are equal using an absolute or relative
|
|
* tolerance test.
|
|
*
|
|
* This is useful to avoid problems due to roundoff error when comparing
|
|
* floating-point values directly. The values are first compared using an
|
|
* absolute tolerance test. If that fails, a relative tolerance test is
|
|
* performed. Use this test if you are unsure of the magnitudes of left and
|
|
* right.
|
|
*
|
|
* @tparam L The length type.
|
|
* @tparam T value value type.
|
|
* @tparam Q The GLM qualifier type.
|
|
*
|
|
* @param left The first value to compare.
|
|
* @param right The other value to compare.
|
|
* @param relativeEpsilon The maximum inclusive delta between `left` and
|
|
* `right` for the relative tolerance test.
|
|
* @param absoluteEpsilon The maximum inclusive delta between `left` and
|
|
* `right` for the absolute tolerance test.
|
|
* @returns `true` if the values are equal within the epsilon; otherwise,
|
|
* `false`.
|
|
*/
|
|
template <glm::length_t L, typename T, glm::qualifier Q>
|
|
static constexpr bool equalsEpsilon(
|
|
const glm::vec<L, T, Q>& left,
|
|
const glm::vec<L, T, Q>& right,
|
|
double relativeEpsilon,
|
|
double absoluteEpsilon) noexcept {
|
|
const glm::vec<L, T, Q> diff = glm::abs(left - right);
|
|
return glm::lessThanEqual(diff, glm::vec<L, T, Q>(absoluteEpsilon)) ==
|
|
glm::vec<L, bool, Q>(true) ||
|
|
glm::lessThanEqual(
|
|
diff,
|
|
relativeEpsilonToAbsolute(left, right, relativeEpsilon)) ==
|
|
glm::vec<L, bool, Q>(true);
|
|
}
|
|
|
|
/**
|
|
* @brief Returns the sign of the value.
|
|
*
|
|
* This is 1 if the value is positive, -1 if the value is
|
|
* negative, or 0 if the value is 0.
|
|
*
|
|
* @param value The value to return the sign of.
|
|
* @returns The sign of value.
|
|
*/
|
|
static constexpr double sign(double value) noexcept {
|
|
if (value == 0.0 || value != value) {
|
|
// zero or NaN
|
|
return value;
|
|
}
|
|
return value > 0 ? 1 : -1;
|
|
}
|
|
|
|
/**
|
|
* @brief Returns 1.0 if the given value is positive or zero, and -1.0 if it
|
|
* is negative.
|
|
*
|
|
* This is similar to {@link Math::sign} except that returns 1.0 instead of
|
|
* 0.0 when the input value is 0.0.
|
|
*
|
|
* @param value The value to return the sign of.
|
|
* @returns The sign of value.
|
|
*/
|
|
static constexpr double signNotZero(double value) noexcept {
|
|
return value < 0.0 ? -1.0 : 1.0;
|
|
}
|
|
|
|
/**
|
|
* @brief Produces an angle in the range -Pi <= angle <= Pi which is
|
|
* equivalent to the provided angle.
|
|
*
|
|
* @param angle The angle in radians.
|
|
* @returns The angle in the range [`-Math::OnePi`, `Math::OnePi`].
|
|
*/
|
|
static double negativePiToPi(double angle) noexcept {
|
|
if (angle >= -Math::OnePi && angle <= Math::OnePi) {
|
|
// Early exit if the input is already inside the range. This avoids
|
|
// unnecessary math which could introduce floating point error.
|
|
return angle;
|
|
}
|
|
return Math::zeroToTwoPi(angle + Math::OnePi) - Math::OnePi;
|
|
}
|
|
|
|
/**
|
|
* @brief Produces an angle in the range 0 <= angle <= 2Pi which is equivalent
|
|
* to the provided angle.
|
|
*
|
|
* @param angle The angle in radians.
|
|
* @returns The angle in the range [0, `Math::TwoPi`].
|
|
*/
|
|
static double zeroToTwoPi(double angle) noexcept {
|
|
if (angle >= 0 && angle <= Math::TwoPi) {
|
|
// Early exit if the input is already inside the range. This avoids
|
|
// unnecessary math which could introduce floating point error.
|
|
return angle;
|
|
}
|
|
const double mod = Math::mod(angle, Math::TwoPi);
|
|
if (glm::abs(mod) < Math::Epsilon14 && glm::abs(angle) > Math::Epsilon14) {
|
|
return Math::TwoPi;
|
|
}
|
|
return mod;
|
|
}
|
|
|
|
/**
|
|
* @brief The modulo operation that also works for negative dividends.
|
|
*
|
|
* @param m The dividend.
|
|
* @param n The divisor.
|
|
* @returns The remainder.
|
|
*/
|
|
static double mod(double m, double n) noexcept {
|
|
if (Math::sign(m) == Math::sign(n) && glm::abs(m) < glm::abs(n)) {
|
|
// Early exit if the input does not need to be modded. This avoids
|
|
// unnecessary math which could introduce floating point error.
|
|
return m;
|
|
}
|
|
return fmod(fmod(m, n) + n, n);
|
|
}
|
|
|
|
/**
|
|
* @brief Converts degrees to radians.
|
|
*
|
|
* @param angleDegrees The angle to convert in degrees.
|
|
* @returns The corresponding angle in radians.
|
|
*/
|
|
static constexpr double degreesToRadians(double angleDegrees) noexcept {
|
|
return angleDegrees * Math::OnePi / 180.0;
|
|
}
|
|
|
|
/**
|
|
* @brief Converts radians to degrees.
|
|
*
|
|
* @param angleRadians The angle to convert in radians.
|
|
* @returns The corresponding angle in degrees.
|
|
*/
|
|
static constexpr double radiansToDegrees(double angleRadians) noexcept {
|
|
return angleRadians * 180.0 / Math::OnePi;
|
|
}
|
|
|
|
/**
|
|
* @brief Computes the linear interpolation of two values.
|
|
*
|
|
* @param p The start value to interpolate.
|
|
* @param q The end value to interpolate.
|
|
* @param time The time of interpolation generally in the range `[0.0, 1.0]`.
|
|
* @returns The linearly interpolated value.
|
|
*
|
|
* @snippet TestMath.cpp lerp
|
|
*/
|
|
static double lerp(double p, double q, double time) noexcept {
|
|
return glm::mix(p, q, time);
|
|
}
|
|
|
|
/**
|
|
* @brief Constrain a value to lie between two values.
|
|
*
|
|
* @param value The value to constrain.
|
|
* @param min The minimum value.
|
|
* @param max The maximum value.
|
|
* @returns The value clamped so that min <= value <= max.
|
|
*/
|
|
static constexpr double clamp(double value, double min, double max) noexcept {
|
|
return glm::clamp(value, min, max);
|
|
};
|
|
|
|
/**
|
|
* @brief Converts a scalar value in the range [-1.0, 1.0] to a SNORM in the
|
|
* range [0, rangeMaximum]
|
|
*
|
|
* @param value The scalar value in the range [-1.0, 1.0].
|
|
* @param rangeMaximum The maximum value in the mapped range, 255 by default.
|
|
* @returns A SNORM value, where 0 maps to -1.0 and rangeMaximum maps to 1.0.
|
|
*
|
|
* @see Math::fromSNorm
|
|
*/
|
|
static double toSNorm(double value, double rangeMaximum = 255.0) noexcept {
|
|
return glm::round(
|
|
(Math::clamp(value, -1.0, 1.0) * 0.5 + 0.5) * rangeMaximum);
|
|
};
|
|
/**
|
|
* @brief Converts a SNORM value in the range [0, rangeMaximum] to a scalar in
|
|
* the range [-1.0, 1.0].
|
|
*
|
|
* @param value SNORM value in the range [0, rangeMaximum].
|
|
* @param rangeMaximum The maximum value in the SNORM range, 255 by default.
|
|
* @returns Scalar in the range [-1.0, 1.0].
|
|
*
|
|
* @see Math::toSNorm
|
|
*/
|
|
static constexpr double
|
|
fromSNorm(double value, double rangeMaximum = 255.0) noexcept {
|
|
return (Math::clamp(value, 0.0, rangeMaximum) / rangeMaximum) * 2.0 - 1.0;
|
|
}
|
|
|
|
/**
|
|
* Converts a longitude value, in radians, to the range [`-Math::OnePi`,
|
|
* `Math::OnePi`).
|
|
*
|
|
* @param angle The longitude value, in radians, to convert to the range
|
|
* [`-Math::OnePi`, `Math::OnePi`).
|
|
* @returns The equivalent longitude value in the range [`-Math::OnePi`,
|
|
* `Math::OnePi`).
|
|
*
|
|
* @snippet TestMath.cpp convertLongitudeRange
|
|
*/
|
|
static double convertLongitudeRange(double angle) noexcept {
|
|
constexpr double twoPi = Math::TwoPi;
|
|
|
|
const double simplified = angle - glm::floor(angle / twoPi) * twoPi;
|
|
|
|
if (simplified < -Math::OnePi) {
|
|
return simplified + twoPi;
|
|
}
|
|
if (simplified >= Math::OnePi) {
|
|
return simplified - twoPi;
|
|
}
|
|
|
|
return simplified;
|
|
};
|
|
|
|
/**
|
|
* @brief Rounds a value up to the nearest integer, like `ceil`, except
|
|
* that if the value is very close to the lower integer it is rounded down
|
|
* (like `floor`) instead.
|
|
*
|
|
* @param value The value to round.
|
|
* @param tolerance The tolerance. If the value is closer than this to the
|
|
* lower integer, it is rounded down instead.
|
|
* @return The rounded value.
|
|
*/
|
|
static double roundUp(double value, double tolerance) noexcept {
|
|
const double up = glm::ceil(value);
|
|
const double down = glm::floor(value);
|
|
if (value - down < tolerance) {
|
|
return down;
|
|
} else {
|
|
return up;
|
|
}
|
|
}
|
|
|
|
/**
|
|
* @brief Rounds a value down to the nearest integer, like `floor`, except
|
|
* that if the value is very close to the higher integer it is rounded up
|
|
* (like `ceil`) instead.
|
|
*
|
|
* @param value The value to round.
|
|
* @param tolerance The tolerance. If the value is closer than this to the
|
|
* higher integer, it is rounded up instead.
|
|
* @return The rounded value.
|
|
*/
|
|
static double roundDown(double value, double tolerance) noexcept {
|
|
const double up = glm::ceil(value);
|
|
const double down = glm::floor(value);
|
|
if (up - value < tolerance) {
|
|
return up;
|
|
} else {
|
|
return down;
|
|
}
|
|
}
|
|
|
|
/**
|
|
* @brief Construct a vector perpendicular to the argument.
|
|
* @param v The input vector
|
|
* @return A vector perpendicular to the input vector
|
|
*/
|
|
template <typename T, glm::qualifier Q>
|
|
static glm::vec<3, T, Q> perpVec(const glm::vec<3, T, Q>& v) {
|
|
// This constructs a vector whose dot product with v will be 0, hence
|
|
// perpendicular to v. As seen in the "Physically Based Rendering".
|
|
if (std::abs(v.x) > std::abs(v.y)) {
|
|
return glm::vec<3, T, Q>(-v.z, 0, v.x) / std::sqrt(v.x * v.x + v.z * v.z);
|
|
}
|
|
return glm::vec<3, T, Q>(0, v.z, -v.y) / std::sqrt(v.y * v.y + v.z * v.z);
|
|
}
|
|
|
|
/** @brief Compute the rotation between two unit vectors.
|
|
* @param vec1 The first vector.
|
|
* @param vec2 The second vector.
|
|
* @return A quaternion representing the rotation of vec1 to vec2.
|
|
*/
|
|
template <typename T, glm::qualifier Q>
|
|
static glm::qua<T, Q>
|
|
rotation(const glm::vec<3, T, Q>& vec1, const glm::vec<3, T, Q>& vec2) {
|
|
// If we take the dot and cross products of the two vectors and store
|
|
// them in a quaternion, that quaternion represents twice the required
|
|
// rotation. We get the correct quaternion by "averaging" with the zero
|
|
// rotation quaternion, in a way analagous to finding the half vector
|
|
// between two 3D vectors.
|
|
auto cosRot = dot(vec1, vec2);
|
|
auto rotAxis = cross(vec1, vec2);
|
|
auto rotAxisLen2 = dot(rotAxis, rotAxis);
|
|
// Not using epsilon for these tests. If abs(cosRot) < 1.0, we can still
|
|
// create a sensible rotation.
|
|
if (cosRot >= 1 || (rotAxisLen2 == 0 && cosRot > 0)) {
|
|
// zero rotation
|
|
return glm::qua<T, Q>(1, 0, 0, 0);
|
|
}
|
|
if (cosRot <= -1 || (rotAxisLen2 == 0 && cosRot < 0)) {
|
|
auto perpAxis = CesiumUtility::Math::perpVec(vec1);
|
|
// rotation by pi radians
|
|
return glm::qua<T, Q>(0, perpAxis);
|
|
}
|
|
|
|
glm::qua<T, Q> sumQuat(cosRot + 1, rotAxis);
|
|
return normalize(sumQuat);
|
|
}
|
|
};
|
|
|
|
} // namespace CesiumUtility
|