Math.h

#pragma once
 
#include "Library.h"
 
#include <glm/gtc/epsilon.hpp>
 
namespace CesiumUtility {
 
class CESIUMUTILITY_API Math final {
public:
  static constexpr double Epsilon1 = 1e-1;
 
  static constexpr double Epsilon2 = 1e-2;
 
  static constexpr double Epsilon3 = 1e-3;
 
  static constexpr double Epsilon4 = 1e-4;
 
  static constexpr double Epsilon5 = 1e-5;
 
  static constexpr double Epsilon6 = 1e-6;
 
  static constexpr double Epsilon7 = 1e-7;
 
  static constexpr double Epsilon8 = 1e-8;
 
  static constexpr double Epsilon9 = 1e-9;
 
  static constexpr double Epsilon10 = 1e-10;
 
  static constexpr double Epsilon11 = 1e-11;
 
  static constexpr double Epsilon12 = 1e-12;
 
  static constexpr double Epsilon13 = 1e-13;
 
  static constexpr double Epsilon14 = 1e-14;
 
  static constexpr double Epsilon15 = 1e-15;
 
  static constexpr double Epsilon16 = 1e-16;
 
  static constexpr double Epsilon17 = 1e-17;
 
  static constexpr double Epsilon18 = 1e-18;
 
  static constexpr double Epsilon19 = 1e-19;
 
  static constexpr double Epsilon20 = 1e-20;
 
  static constexpr double Epsilon21 = 1e-21;
 
  static constexpr double OnePi = 3.14159265358979323846;
 
  static constexpr double TwoPi = OnePi * 2.0;
 
  static constexpr double PiOverTwo = OnePi / 2.0;
 
  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));
  }
 
  static constexpr double relativeEpsilonToAbsolute(
      double a,
      double b,
      double relativeEpsilon) noexcept {
    return relativeEpsilon * glm::max(glm::abs(a), glm::abs(b));
  }
 
  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);
  }
 
  static constexpr bool
  equalsEpsilon(double left, double right, double relativeEpsilon) noexcept {
    return equalsEpsilon(left, right, relativeEpsilon, relativeEpsilon);
  }
 
  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);
  }
 
  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);
  }
 
  static constexpr double sign(double value) noexcept {
    if (value == 0.0 || value != value) {
      // zero or NaN
      return value;
    }
    return value > 0 ? 1 : -1;
  }
 
  static constexpr double signNotZero(double value) noexcept {
    return value < 0.0 ? -1.0 : 1.0;
  }
 
  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;
  }
 
  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;
  }
 
  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);
  }
 
  static constexpr double degreesToRadians(double angleDegrees) noexcept {
    return angleDegrees * Math::OnePi / 180.0;
  }
 
  static constexpr double radiansToDegrees(double angleRadians) noexcept {
    return angleRadians * 180.0 / Math::OnePi;
  }
 
  static double lerp(double p, double q, double time) noexcept {
    return glm::mix(p, q, time);
  }
 
  static constexpr double clamp(double value, double min, double max) noexcept {
    return glm::clamp(value, min, max);
  };
 
  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);
  };
  static constexpr double
  fromSNorm(double value, double rangeMaximum = 255.0) noexcept {
    return (Math::clamp(value, 0.0, rangeMaximum) / rangeMaximum) * 2.0 - 1.0;
  }
 
  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;
  };
 
  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;
    }
  }
 
  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;
    }
  }
 
  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);
  }
 
  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