util/math.h
| Line | Branch | Exec | Source |
|---|---|---|---|
| 1 | // SPDX-License-Identifier: BSD-3-Clause | ||
| 2 | /// @file math.h | ||
| 3 | /// @brief Small scalar functions shared across Ernie's numerics. | ||
| 4 | #pragma once | ||
| 5 | |||
| 6 | #include <glm/glm.hpp> | ||
| 7 | |||
| 8 | #include <algorithm> | ||
| 9 | #include <array> | ||
| 10 | #include <cmath> | ||
| 11 | #include <cstddef> | ||
| 12 | #include <numbers> | ||
| 13 | #include <span> | ||
| 14 | #include <utility> | ||
| 15 | |||
| 16 | namespace ernie { | ||
| 17 | |||
| 18 | /// @brief The fractional part of a number. | ||
| 19 | /// @param x The number. | ||
| 20 | /// @return x minus its floor, in [0, 1). | ||
| 21 | 82283488 | inline double Fract(double x) { | |
| 22 | 82283488 | return x - std::floor(x); | |
| 23 | } | ||
| 24 | |||
| 25 | /// @brief Floored modulo, which wraps negative numbers into [0, m). | ||
| 26 | /// @param x The number. | ||
| 27 | /// @param m The modulus. | ||
| 28 | /// @return x modulo m, with the sign of m. | ||
| 29 | 10283448 | inline double Mod(double x, double m) { | |
| 30 | 10283448 | return x - m * std::floor(x / m); | |
| 31 | } | ||
| 32 | |||
| 33 | /// @brief Hermite interpolation from 0 to 1 between two edges. | ||
| 34 | /// @param edge0 Where the result starts rising from 0. | ||
| 35 | /// @param edge1 Where the result reaches 1. | ||
| 36 | /// @param x The value. | ||
| 37 | /// @return 0 below edge0, 1 above edge1, and a smooth cubic between. | ||
| 38 | 49156 | inline double SmoothStep(double edge0, double edge1, double x) { | |
| 39 | 49156 | const double t = std::clamp((x - edge0) / (edge1 - edge0), 0.0, 1.0); | |
| 40 | 49156 | return t * t * (3.0 - 2.0 * t); | |
| 41 | } | ||
| 42 | |||
| 43 | /// @brief Trilinear interpolation between the eight corners of a cell. | ||
| 44 | /// @param corners The corner values, indexed x + 2y + 4z. | ||
| 45 | /// @param t The position within the cell, each component in [0, 1]. | ||
| 46 | /// @return The interpolated value. | ||
| 47 | 10294452 | inline double Trilinear(const std::array<double, 8>& corners, const glm::dvec3& t) { | |
| 48 | 10294452 | const double z0 = std::lerp(std::lerp(corners[0], corners[1], t.x), | |
| 49 | 10294452 | std::lerp(corners[2], corners[3], t.x), t.y); | |
| 50 | 10294452 | const double z1 = std::lerp(std::lerp(corners[4], corners[5], t.x), | |
| 51 | 10294452 | std::lerp(corners[6], corners[7], t.x), t.y); | |
| 52 | 10294452 | return std::lerp(z0, z1, t.z); | |
| 53 | } | ||
| 54 | |||
| 55 | /// @brief Cosine spacing: sample positions in [0, 1] that crowd towards both ends. | ||
| 56 | /// @param i The sample. | ||
| 57 | /// @param count How many samples, at least 2. | ||
| 58 | /// @return (1 - cos(pi i / (count - 1))) / 2. | ||
| 59 | 97918 | inline double CosineSpacing(std::size_t i, std::size_t count) { | |
| 60 | 97922 | const double u = static_cast<double>(i) / static_cast<double>(count - 1); | |
| 61 | 97918 | return 0.5 * (1.0 - std::cos(std::numbers::pi * u)); | |
| 62 | } | ||
| 63 | |||
| 64 | /// @brief The interval of a table holding a position, and the fraction across it. | ||
| 65 | /// @param xs The sample positions, increasing, at least two. | ||
| 66 | /// @param x The position. | ||
| 67 | /// @return The interval's first index and the fraction across it, clamped to [0, 1]; 0 across an | ||
| 68 | /// interval of no width. | ||
| 69 | 67990380 | inline std::pair<std::size_t, double> Bracket(std::span<const double> xs, double x) { | |
| 70 | 67990380 | const std::ptrdiff_t upper = std::upper_bound(xs.begin(), xs.end(), x) - xs.begin(); | |
| 71 | 67990380 | const auto i = static_cast<std::size_t>( | |
| 72 | 67990380 | std::clamp<std::ptrdiff_t>(upper - 1, 0, static_cast<std::ptrdiff_t>(xs.size()) - 2)); | |
| 73 | 67990380 | const double width = xs[i + 1] - xs[i]; | |
| 74 | 135980759 | return {i, width > 0.0 ? std::clamp((x - xs[i]) / width, 0.0, 1.0) : 0.0}; | |
| 75 | } | ||
| 76 | |||
| 77 | /// @brief Piecewise-linear interpolation of a table, constant beyond its ends. | ||
| 78 | /// @param xs The sample positions, increasing. | ||
| 79 | /// @param values The value at each position. | ||
| 80 | /// @param x Where to interpolate. | ||
| 81 | /// @return The interpolated value, or a default value for an empty table. | ||
| 82 | template <typename T> | ||
| 83 | 67861246 | T Interpolate(std::span<const double> xs, std::span<const T> values, double x) { | |
| 84 | 67861246 | if (xs.empty()) { | |
| 85 | ✗ | return T{}; | |
| 86 | } | ||
| 87 | 67861244 | if (x <= xs.front()) { | |
| 88 | 15 | return values.front(); | |
| 89 | } | ||
| 90 | 67861229 | if (x >= xs.back()) { | |
| 91 | 539 | return values.back(); | |
| 92 | } | ||
| 93 | 67860690 | const auto [i, t] = Bracket(xs, x); | |
| 94 | 67860690 | return values[i] + (values[i + 1] - values[i]) * t; | |
| 95 | } | ||
| 96 | |||
| 97 | /// @return How many groups of `size` cover `count`, as workgroups cover threads. | ||
| 98 | 123 | template <typename T> constexpr T DivideRoundUp(T count, T size) { | |
| 99 | 123 | return (count + size - 1) / size; | |
| 100 | } | ||
| 101 | |||
| 102 | } // namespace ernie | ||
| 103 |