physics/spectrum.cpp
| Line | Branch | Exec | Source |
|---|---|---|---|
| 1 | // SPDX-License-Identifier: BSD-3-Clause | ||
| 2 | /// @file spectrum.cpp | ||
| 3 | /// @brief Spectra, black bodies, colour matching and RGB upsampling. | ||
| 4 | #include "ernie/physics/spectrum.h" | ||
| 5 | |||
| 6 | #include "ernie/util/linear_algebra.h" | ||
| 7 | #include "ernie/util/math.h" | ||
| 8 | |||
| 9 | #include <algorithm> | ||
| 10 | #include <array> | ||
| 11 | #include <cmath> | ||
| 12 | #include <map> | ||
| 13 | #include <mutex> | ||
| 14 | #include <numbers> | ||
| 15 | #include <optional> | ||
| 16 | #include <stdexcept> | ||
| 17 | #include <string> | ||
| 18 | #include <tuple> | ||
| 19 | #include <type_traits> | ||
| 20 | #include <utility> | ||
| 21 | |||
| 22 | namespace ernie { | ||
| 23 | |||
| 24 | namespace { | ||
| 25 | |||
| 26 | /// The visible range RGB mode's channels are integrated over, in nm. | ||
| 27 | constexpr double kVisibleMinNm = 380.0; | ||
| 28 | constexpr double kVisibleMaxNm = 780.0; | ||
| 29 | |||
| 30 | /// The range sigmoid spectra are defined and fitted over, in nm. | ||
| 31 | constexpr double kSigmoidMinNm = 360.0; | ||
| 32 | constexpr double kSigmoidMaxNm = 830.0; | ||
| 33 | |||
| 34 | /// Sub-samples per band when averaging analytic functions, enough for 1e-4 relative accuracy with | ||
| 35 | /// bands up to 100 nm. | ||
| 36 | constexpr int kSubSamples = 32; | ||
| 37 | |||
| 38 | /// Below this many bands, colour weights are fitted rather than integrated. | ||
| 39 | constexpr int kFitBelow = 16; | ||
| 40 | |||
| 41 | // ---------------------------------------------------------------------------------------------- | ||
| 42 | // Colour matching | ||
| 43 | |||
| 44 | /// Piecewise Gaussian of Wyman et al.: width `below` left of `mean` and `above` right of it. | ||
| 45 | 1261568 | double Lobe(double wavelength, double mean, double below, double above) { | |
| 46 | 1261568 | const double t = (wavelength - mean) / (wavelength < mean ? below : above); | |
| 47 | 1261568 | return std::exp(-0.5 * t * t); | |
| 48 | } | ||
| 49 | |||
| 50 | /// A spectrum's CIE XYZ, integrated in 1 nm steps at their midpoints over [lowNm, highNm). | ||
| 51 | template <typename Function> | ||
| 52 | 128 | glm::dvec3 IntegrateXyz(Function spectrum, double lowNm, double highNm) { | |
| 53 | 128 | glm::dvec3 xyz(0.0); | |
| 54 | 51328 | for (double w = lowNm + 0.5; w < highNm; w += 1.0) { | |
| 55 | 51200 | xyz += spectrum(w) * CieXyz(w); | |
| 56 | } | ||
| 57 | 128 | return xyz; | |
| 58 | } | ||
| 59 | |||
| 60 | /// The colour matching functions every 1 nm over the sigmoid range, computed once. | ||
| 61 | 6852 | const std::vector<glm::dvec3>& SampledCmfs() { | |
| 62 | 3 | static const std::vector<glm::dvec3> cmfs = [] { | |
| 63 | 3 | std::vector<glm::dvec3> values; | |
| 64 | 1416 | for (double w = kSigmoidMinNm; w <= kSigmoidMaxNm; w += 1.0) { | |
| 65 | 1413 | values.push_back(CieXyz(w)); | |
| 66 | } | ||
| 67 | 3 | return values; | |
| 68 | 6852 | }(); | |
| 69 | 6852 | return cmfs; | |
| 70 | } | ||
| 71 | |||
| 72 | /// A reflectance's XYZ under the equal-energy illuminant, normalised by the integral of y. | ||
| 73 | 6852 | template <typename Function> glm::dvec3 EqualEnergyXyz(Function reflectance) { | |
| 74 | 6852 | const std::vector<glm::dvec3>& cmfs = SampledCmfs(); | |
| 75 | 6852 | glm::dvec3 xyz(0.0); | |
| 76 | 6852 | double y = 0.0; | |
| 77 | 3234144 | for (std::size_t i = 0; i < cmfs.size(); ++i) { | |
| 78 | 3227292 | xyz += reflectance(kSigmoidMinNm + static_cast<double>(i)) * cmfs[i]; | |
| 79 | 3227292 | y += cmfs[i].y; | |
| 80 | } | ||
| 81 | 6852 | return xyz / y; | |
| 82 | } | ||
| 83 | |||
| 84 | /// The Bradford adaptation from the equal-energy white to D65, as one matrix. | ||
| 85 | 3 | glm::dmat3 EqualEnergyToD65() { | |
| 86 | // The cone response matrix, rows written out and transposed into GLM's columns. | ||
| 87 | 3 | const glm::dmat3 bradford = glm::transpose(glm::dmat3(0.8951, 0.2664, -0.1614, // | |
| 88 | -0.7502, 1.7135, 0.0367, // | ||
| 89 | 0.0389, -0.0685, 1.0296)); | ||
| 90 | 3 | const glm::dvec3 source = bradford * EqualEnergyXyz([](double) { return 1.0; }); | |
| 91 | 3 | const glm::dvec3 target = bradford * glm::dvec3(0.95047, 1.0, 1.08883); | |
| 92 | 3 | const glm::dmat3 scale(glm::dvec3(target.x / source.x, 0.0, 0.0), | |
| 93 | 3 | glm::dvec3(0.0, target.y / source.y, 0.0), | |
| 94 | 3 | glm::dvec3(0.0, 0.0, target.z / source.z)); | |
| 95 | 3 | return glm::inverse(bradford) * scale * bradford; | |
| 96 | } | ||
| 97 | |||
| 98 | // ---------------------------------------------------------------------------------------------- | ||
| 99 | // Bands and RGB mode | ||
| 100 | |||
| 101 | /// The equal-energy white's linear sRGB, which RGB mode's channels are relative to. | ||
| 102 | 1 | glm::dvec3 RgbWhite() { | |
| 103 | 1 | static const glm::dvec3 white = EqualEnergyXyzToLinearSrgb( | |
| 104 | 2 | IntegrateXyz([](double) { return 1.0; }, kVisibleMinNm, kVisibleMaxNm)); | |
| 105 | 1 | return white; | |
| 106 | } | ||
| 107 | |||
| 108 | /// The linear sRGB channel RGB-mode band `band` carries of spectrum `f` (blue, green then red, by | ||
| 109 | /// rising wavelength), relative to the equal-energy white's, so a flat spectrum has its value in | ||
| 110 | /// every channel. At least 0, for lights more saturated than sRGB holds. | ||
| 111 | ✗ | template <typename Function> double RgbChannel(int band, Function f) { | |
| 112 | ✗ | const glm::dvec3 xyz = IntegrateXyz(f, kVisibleMinNm, kVisibleMaxNm); | |
| 113 | ✗ | const glm::dvec3 rgb = EqualEnergyXyzToLinearSrgb(xyz) / RgbWhite(); | |
| 114 | ✗ | return std::max(rgb[2 - band], 0.0); | |
| 115 | } | ||
| 116 | |||
| 117 | /// The mean of `f` over a band, from midpoint sub-samples. In RGB mode a scalar spectrum's band | ||
| 118 | /// value is its RGB channel instead. | ||
| 119 | 224 | template <typename Function> auto BandAverage(const SpectralBands& bands, int band, Function f) { | |
| 120 | using Result = decltype(f(0.0)); | ||
| 121 | if constexpr (std::is_same_v<Result, double>) { | ||
| 122 | 80 | if (bands.rgb) { | |
| 123 | ✗ | return RgbChannel(band, f); | |
| 124 | } | ||
| 125 | } | ||
| 126 | 144 | Result sum{}; | |
| 127 | 7392 | for (int k = 0; k < kSubSamples; ++k) { | |
| 128 | 7168 | sum += f(bands.Lower(band) + (k + 0.5) * bands.Width() / kSubSamples); | |
| 129 | } | ||
| 130 | 224 | return sum / static_cast<double>(kSubSamples); | |
| 131 | } | ||
| 132 | |||
| 133 | // ---------------------------------------------------------------------------------------------- | ||
| 134 | // Colour weights | ||
| 135 | |||
| 136 | /// RGB mode's weights: each channel's XYZ, inverting `RgbChannel`'s projection, so band values of | ||
| 137 | /// 1 give the equal-energy white's XYZ. | ||
| 138 | 1 | std::vector<glm::dvec3> RgbModeColourWeights() { | |
| 139 | 2 | const glm::dmat3 toRgb(EqualEnergyXyzToLinearSrgb({1.0, 0.0, 0.0}), | |
| 140 | 1 | EqualEnergyXyzToLinearSrgb({0.0, 1.0, 0.0}), | |
| 141 | 1 | EqualEnergyXyzToLinearSrgb({0.0, 0.0, 1.0})); | |
| 142 | 1 | const glm::dmat3 toXyz = glm::inverse(toRgb); | |
| 143 | 1 | const glm::dvec3 white = RgbWhite(); | |
| 144 | 1 | std::vector<glm::dvec3> weights(3); | |
| 145 | 4 | for (int band = 0; band < 3; ++band) { | |
| 146 | 3 | glm::dvec3 channel(0.0); | |
| 147 | 3 | channel[2 - band] = white[2 - band]; | |
| 148 | 3 | weights[static_cast<std::size_t>(band)] = toXyz * channel; | |
| 149 | } | ||
| 150 | 1 | return weights; | |
| 151 | } | ||
| 152 | |||
| 153 | /// Colours of every hue and saturation at full brightness; the fit is linear, so darker ones | ||
| 154 | /// follow. Greys are left out: white is added separately with a large weight. | ||
| 155 | 1 | std::vector<glm::dvec3> SaturatedColours() { | |
| 156 | 1 | constexpr int kSteps = 6; | |
| 157 | 1 | std::vector<glm::dvec3> colours; | |
| 158 | 8 | for (int r = 0; r <= kSteps; ++r) { | |
| 159 | 56 | for (int g = 0; g <= kSteps; ++g) { | |
| 160 | 392 | for (int b = 0; b <= kSteps; ++b) { | |
| 161 | 343 | const glm::dvec3 rgb = glm::dvec3(r, g, b) / static_cast<double>(kSteps); | |
| 162 | 343 | const double top = std::max({rgb.r, rgb.g, rgb.b}); | |
| 163 | 343 | if (top == 1.0 && top != std::min({rgb.r, rgb.g, rgb.b})) { | |
| 164 | 126 | colours.push_back(rgb); | |
| 165 | } | ||
| 166 | } | ||
| 167 | } | ||
| 168 | } | ||
| 169 | 1 | return colours; | |
| 170 | ✗ | } | |
| 171 | |||
| 172 | /// Weights w minimising the error of sum_i w_i a_i over many reflectances, a being a reflectance | ||
| 173 | /// sampled at the band centres as shaders sample it, against its exact XYZ. | ||
| 174 | 1 | std::vector<glm::dvec3> FitColourWeights(const SpectralBands& bands) { | |
| 175 | 1 | const auto n = static_cast<std::size_t>(bands.count); | |
| 176 | 1 | std::vector<double> normal(n * n, 0.0); | |
| 177 | 1 | std::vector<glm::dvec3> right(n, glm::dvec3(0.0)); | |
| 178 | 128 | const auto add = [&](const SigmoidSpectrum& spectrum, double weight) { | |
| 179 | 127 | std::vector<double> sampled(n); | |
| 180 | 508 | for (std::size_t i = 0; i < n; ++i) { | |
| 181 | 381 | sampled[i] = spectrum(bands.Centre(static_cast<int>(i))); | |
| 182 | } | ||
| 183 | 127 | const glm::dvec3 xyz = IntegrateXyz(spectrum, bands.minNm, bands.maxNm); | |
| 184 | 508 | for (std::size_t i = 0; i < n; ++i) { | |
| 185 | 1524 | for (std::size_t j = 0; j < n; ++j) { | |
| 186 | 1143 | normal[i * n + j] += weight * sampled[i] * sampled[j]; | |
| 187 | } | ||
| 188 | 381 | right[i] += weight * sampled[i] * xyz; | |
| 189 | } | ||
| 190 | 128 | }; | |
| 191 | 127 | for (const glm::dvec3& rgb : SaturatedColours()) { | |
| 192 | 126 | add(RgbToSigmoid(rgb), 1.0); | |
| 193 | 1 | } | |
| 194 | 1 | add(SigmoidSpectrum{}, 1e3); // White, weighted to keep greys exact. | |
| 195 | 2 | return SolveLinearSystem(std::move(normal), std::move(right)); | |
| 196 | 1 | } | |
| 197 | |||
| 198 | /// `FitColourWeights`, computed once per band layout since renderers ask every frame. | ||
| 199 | 1 | const std::vector<glm::dvec3>& CachedColourFit(const SpectralBands& bands) { | |
| 200 | 1 | static std::mutex mutex; | |
| 201 | 1 | static std::map<std::tuple<double, double, int>, std::vector<glm::dvec3>> fitted; | |
| 202 | 1 | const std::lock_guard lock(mutex); | |
| 203 | 1 | const auto key = std::tuple(bands.minNm, bands.maxNm, bands.count); | |
| 204 | 1 | if (const auto it = fitted.find(key); it != fitted.end()) { | |
| 205 | ✗ | return it->second; | |
| 206 | } | ||
| 207 | 1 | return fitted.emplace(key, FitColourWeights(bands)).first->second; | |
| 208 | 1 | } | |
| 209 | |||
| 210 | // ---------------------------------------------------------------------------------------------- | ||
| 211 | // CIE daylight | ||
| 212 | |||
| 213 | /// CIE daylight components S0, S1 and S2 from 300 to 830 nm every 10 nm (CIE 015:2018, table 6). | ||
| 214 | constexpr double kDaylightFirstNm = 300.0; | ||
| 215 | constexpr double kDaylightStepNm = 10.0; | ||
| 216 | constexpr std::array<double, 54> kDaylightS0{ | ||
| 217 | 0.04, 6.0, 29.6, 55.3, 57.3, 61.8, 61.5, 68.8, 63.4, 65.8, 94.8, 104.8, 105.9, 96.8, | ||
| 218 | 113.9, 125.6, 125.5, 121.3, 121.3, 113.5, 113.1, 110.8, 106.5, 108.8, 105.3, 104.4, 100.0, 96.0, | ||
| 219 | 95.1, 89.1, 90.5, 90.3, 88.4, 84.0, 85.1, 81.9, 82.6, 84.9, 81.3, 71.9, 74.3, 76.4, | ||
| 220 | 63.3, 71.7, 77.0, 65.2, 47.7, 68.6, 65.0, 66.0, 61.0, 53.3, 58.9, 61.9}; | ||
| 221 | constexpr std::array<double, 54> kDaylightS1{ | ||
| 222 | 0.02, 4.5, 22.4, 42.0, 40.6, 41.6, 38.0, 42.4, 38.5, 35.0, 43.4, 46.3, 43.9, 37.1, | ||
| 223 | 36.7, 35.9, 32.6, 27.9, 24.3, 20.1, 16.2, 13.2, 8.6, 6.1, 4.2, 1.9, 0.0, -1.6, | ||
| 224 | -3.5, -3.5, -5.8, -7.2, -8.6, -9.5, -10.9, -10.7, -12.0, -14.0, -13.6, -12.0, -13.3, -12.9, | ||
| 225 | -10.6, -11.6, -12.2, -10.2, -7.8, -11.2, -10.4, -10.6, -9.7, -8.3, -9.3, -9.8}; | ||
| 226 | constexpr std::array<double, 54> kDaylightS2{ | ||
| 227 | 0.0, 2.0, 4.0, 8.5, 7.8, 6.7, 5.3, 6.1, 3.0, 1.2, -1.1, -0.5, -0.7, -1.2, | ||
| 228 | -2.6, -2.9, -2.8, -2.6, -2.6, -1.8, -1.5, -1.3, -1.2, -1.0, -0.5, -0.3, 0.0, 0.2, | ||
| 229 | 0.5, 2.1, 3.2, 4.1, 4.7, 5.1, 6.7, 7.3, 8.6, 9.8, 10.2, 8.3, 9.6, 8.5, | ||
| 230 | 7.0, 7.6, 8.0, 6.7, 5.2, 7.4, 6.8, 7.0, 6.4, 5.5, 6.1, 6.5}; | ||
| 231 | |||
| 232 | /// The daylight locus: the chromaticity (x, y) of daylight at a correlated colour temperature. | ||
| 233 | 804 | glm::dvec2 DaylightChromaticity(double kelvin) { | |
| 234 | 804 | const double t = kelvin; | |
| 235 | 1608 | const double x = t <= 7000.0 | |
| 236 | 804 | ? -4.6070e9 / (t * t * t) + 2.9678e6 / (t * t) + 0.09911e3 / t + 0.244063 | |
| 237 | ✗ | : -2.0064e9 / (t * t * t) + 1.9018e6 / (t * t) + 0.24748e3 / t + 0.237040; | |
| 238 | 804 | return {x, -3.000 * x * x + 2.870 * x - 0.275}; | |
| 239 | } | ||
| 240 | |||
| 241 | /// The weights M1 and M2 of components S1 and S2 for daylight of chromaticity `xy`. | ||
| 242 | 804 | glm::dvec2 DaylightComponentWeights(glm::dvec2 xy) { | |
| 243 | 804 | const double m = 0.0241 + 0.2562 * xy.x - 0.7341 * xy.y; | |
| 244 | 804 | return {(-1.3515 - 1.7703 * xy.x + 5.9114 * xy.y) / m, | |
| 245 | 804 | (0.0300 - 31.4424 * xy.x + 30.0717 * xy.y) / m}; | |
| 246 | } | ||
| 247 | |||
| 248 | // ---------------------------------------------------------------------------------------------- | ||
| 249 | // Sigmoid spectra | ||
| 250 | |||
| 251 | 3275493 | double Sigmoid(double x) { | |
| 252 | 3275493 | return 0.5 + x / (2.0 * std::sqrt(1.0 + x * x)); | |
| 253 | } | ||
| 254 | |||
| 255 | /// Inverse of `Sigmoid` for values strictly inside (0, 1). | ||
| 256 | 145 | double InverseSigmoid(double y) { | |
| 257 | 145 | return (2.0 * y - 1.0) / (2.0 * std::sqrt(y * (1.0 - y))); | |
| 258 | } | ||
| 259 | |||
| 260 | /// A sigmoid spectrum's colour under illuminant E, adapted to D65. | ||
| 261 | 6843 | glm::dvec3 SigmoidColour(const SigmoidSpectrum& spectrum) { | |
| 262 | 6843 | return EqualEnergyXyzToLinearSrgb(EqualEnergyXyz(spectrum)); | |
| 263 | } | ||
| 264 | |||
| 265 | /// The colour's derivatives with respect to the three coefficients, by forward differences. | ||
| 266 | 837 | glm::dmat3 SigmoidColourJacobian(const SigmoidSpectrum& spectrum) { | |
| 267 | 837 | glm::dmat3 jacobian; | |
| 268 | 3348 | for (int k = 0; k < 3; ++k) { | |
| 269 | 2511 | SigmoidSpectrum nudged = spectrum; | |
| 270 | 2511 | const double h = 1e-5 * std::max(1.0, std::abs(spectrum.coefficients[k])); | |
| 271 | 2511 | nudged.coefficients[k] += h; | |
| 272 | 2511 | jacobian[k] = (SigmoidColour(nudged) - SigmoidColour(spectrum)) / h; | |
| 273 | } | ||
| 274 | 837 | return jacobian; | |
| 275 | } | ||
| 276 | |||
| 277 | /// The Gauss-Newton step `step` from `spectrum`, halved up to ten times until it brings the colour | ||
| 278 | /// closer to `rgb` than `error`; nothing if no length of it does. | ||
| 279 | 837 | std::optional<SigmoidSpectrum> DampedStep(const SigmoidSpectrum& spectrum, const glm::dvec3& step, | |
| 280 | const glm::dvec3& rgb, double error) { | ||
| 281 | 837 | double damping = 1.0; | |
| 282 | 842 | for (int halving = 0; halving < 10; ++halving) { | |
| 283 | 842 | SigmoidSpectrum next = spectrum; | |
| 284 | 842 | next.coefficients = spectrum.coefficients - damping * step; | |
| 285 | 842 | if (glm::length(SigmoidColour(next) - rgb) < error) { | |
| 286 | 837 | return next; | |
| 287 | } | ||
| 288 | 5 | damping *= 0.5; | |
| 289 | } | ||
| 290 | ✗ | return std::nullopt; | |
| 291 | } | ||
| 292 | |||
| 293 | } // namespace | ||
| 294 | |||
| 295 | // ---------------------------------------------------------------------------------------------- | ||
| 296 | // Spectral bands | ||
| 297 | |||
| 298 | 14 | void SpectralBands::Check() const { | |
| 299 | 14 | if (!(minNm > 0.0) || !(maxNm > minNm) || count < 1 || count > kMaxBands) { | |
| 300 | 6 | throw std::invalid_argument("A spectral range needs 0 < min < max and 1 to " + | |
| 301 | 12 | std::to_string(kMaxBands) + " bands"); | |
| 302 | } | ||
| 303 | 11 | if (rgb && count != 3) { | |
| 304 | ✗ | throw std::invalid_argument("RGB mode has 3 bands"); | |
| 305 | } | ||
| 306 | 11 | } | |
| 307 | |||
| 308 | ✗ | double BandValue(const SpectralBands& bands, int band, | |
| 309 | const std::function<double(double)>& spectrum) { | ||
| 310 | ✗ | return BandAverage(bands, band, spectrum); | |
| 311 | } | ||
| 312 | |||
| 313 | // ---------------------------------------------------------------------------------------------- | ||
| 314 | // Tabulated spectra | ||
| 315 | |||
| 316 | 54 | Spectrum::Spectrum(std::vector<double> wavelengthsNm, std::vector<double> values) | |
| 317 | 54 | : wavelengths_(std::move(wavelengthsNm)), values_(std::move(values)) { | |
| 318 | 54 | if (wavelengths_.size() != values_.size()) { | |
| 319 | 1 | throw std::invalid_argument("A spectrum needs one value per wavelength"); | |
| 320 | } | ||
| 321 | 20908 | for (std::size_t i = 1; i < wavelengths_.size(); ++i) { | |
| 322 | 20856 | if (!(wavelengths_[i] > wavelengths_[i - 1])) { | |
| 323 | 1 | throw std::invalid_argument("Spectrum wavelengths must be strictly increasing"); | |
| 324 | } | ||
| 325 | } | ||
| 326 | 54 | } | |
| 327 | |||
| 328 | 2 | Spectrum Spectrum::Constant(double value) { | |
| 329 | 4 | return Spectrum({1.0}, {value}); | |
| 330 | } | ||
| 331 | |||
| 332 | 2 | Spectrum Spectrum::FromPairs(std::span<const float> pairs) { | |
| 333 | 2 | if (pairs.size() % 2 != 0) { | |
| 334 | 1 | throw std::invalid_argument("A spectrum needs (wavelength, value) pairs"); | |
| 335 | } | ||
| 336 | 1 | std::vector<double> wavelengths; | |
| 337 | 1 | std::vector<double> values; | |
| 338 | 3 | for (std::size_t i = 0; i < pairs.size(); i += 2) { | |
| 339 | 2 | wavelengths.push_back(pairs[i]); | |
| 340 | 2 | values.push_back(pairs[i + 1]); | |
| 341 | } | ||
| 342 | 2 | return Spectrum(std::move(wavelengths), std::move(values)); | |
| 343 | 1 | } | |
| 344 | |||
| 345 | 10482 | double Spectrum::operator()(double wavelengthNm) const { | |
| 346 | 10482 | return Interpolate<double>(wavelengths_, values_, wavelengthNm); | |
| 347 | } | ||
| 348 | |||
| 349 | 9 | double Spectrum::Integral(double lowNm, double highNm) const { | |
| 350 | // The interpolant is linear between samples, so the trapezoid rule over the range's ends and | ||
| 351 | // every sample inside it is exact. | ||
| 352 | 9 | std::vector<double> knots{lowNm, highNm}; | |
| 353 | 15256 | for (const double w : wavelengths_) { | |
| 354 | 15247 | if (w > lowNm && w < highNm) { | |
| 355 | 3796 | knots.push_back(w); | |
| 356 | } | ||
| 357 | } | ||
| 358 | 9 | std::sort(knots.begin(), knots.end()); | |
| 359 | 9 | double sum = 0.0; | |
| 360 | 3814 | for (std::size_t i = 1; i < knots.size(); ++i) { | |
| 361 | 3805 | sum += 0.5 * ((*this)(knots[i - 1]) + (*this)(knots[i])) * (knots[i] - knots[i - 1]); | |
| 362 | } | ||
| 363 | 9 | return sum; | |
| 364 | 9 | } | |
| 365 | |||
| 366 | 9 | double Spectrum::Average(double lowNm, double highNm) const { | |
| 367 | 9 | if (!(highNm > lowNm)) { | |
| 368 | ✗ | return (*this)(lowNm); | |
| 369 | } | ||
| 370 | 9 | return Integral(lowNm, highNm) / (highNm - lowNm); | |
| 371 | } | ||
| 372 | |||
| 373 | 1 | std::vector<float> Spectrum::Bin(const SpectralBands& bands) const { | |
| 374 | 1 | bands.Check(); | |
| 375 | 1 | std::vector<float> binned(static_cast<std::size_t>(bands.count)); | |
| 376 | 3 | for (int i = 0; i < bands.count; ++i) { | |
| 377 | 4 | binned[static_cast<std::size_t>(i)] = | |
| 378 | 4 | static_cast<float>(bands.rgb ? RgbChannel(i, [this](double w) { return (*this)(w); }) | |
| 379 | 2 | : Average(bands.Lower(i), bands.Lower(i) + bands.Width())); | |
| 380 | } | ||
| 381 | 1 | return binned; | |
| 382 | ✗ | } | |
| 383 | |||
| 384 | // ---------------------------------------------------------------------------------------------- | ||
| 385 | // Black bodies and daylight | ||
| 386 | |||
| 387 | 127800 | double PlanckRadiance(double wavelengthNm, double kelvin) { | |
| 388 | 127800 | const double lambda = wavelengthNm * 1e-9; | |
| 389 | 127800 | const double exponent = kPlanck * kSpeedOfLight / (lambda * kBoltzmann * kelvin); | |
| 390 | 127800 | const double perMetre = 2.0 * kPlanck * kSpeedOfLight * kSpeedOfLight / | |
| 391 | 127800 | (std::pow(lambda, 5.0) * std::expm1(exponent)); | |
| 392 | 127800 | return perMetre * 1e-9; | |
| 393 | } | ||
| 394 | |||
| 395 | 805 | double CieDaylight(double wavelengthNm, double kelvin) { | |
| 396 | 805 | if (!(kelvin >= 4000.0 && kelvin <= 25000.0)) { | |
| 397 | 1 | throw std::invalid_argument("CIE illuminant D needs a temperature from 4000 to 25000 K"); | |
| 398 | } | ||
| 399 | 804 | const glm::dvec2 m = DaylightComponentWeights(DaylightChromaticity(kelvin)); | |
| 400 | 804 | const double position = (wavelengthNm - kDaylightFirstNm) / kDaylightStepNm; | |
| 401 | 804 | if (position < 0.0 || position > static_cast<double>(kDaylightS0.size() - 1)) { | |
| 402 | return 0.0; | ||
| 403 | } | ||
| 404 | 803 | const auto i = std::min(static_cast<std::size_t>(position), kDaylightS0.size() - 2); | |
| 405 | 803 | const double f = position - static_cast<double>(i); | |
| 406 | 2409 | const auto daylight = [&](std::size_t k) { | |
| 407 | 1606 | return kDaylightS0[k] + m.x * kDaylightS1[k] + m.y * kDaylightS2[k]; | |
| 408 | 803 | }; | |
| 409 | 803 | return (1.0 - f) * daylight(i) + f * daylight(i + 1); | |
| 410 | } | ||
| 411 | |||
| 412 | 3 | std::vector<float> BlackbodyIrradiance(const SpectralBands& bands, double kelvin, | |
| 413 | double totalIrradiance) { | ||
| 414 | 3 | bands.Check(); | |
| 415 | 3 | if (!(kelvin > 0.0)) { | |
| 416 | 1 | throw std::invalid_argument("A black body needs a positive temperature"); | |
| 417 | } | ||
| 418 | 2 | const double totalRadiance = kStefanBoltzmann * std::pow(kelvin, 4.0) / std::numbers::pi; | |
| 419 | 2 | std::vector<float> irradiance(static_cast<std::size_t>(bands.count)); | |
| 420 | 82 | for (int i = 0; i < bands.count; ++i) { | |
| 421 | 80 | const double radiance = | |
| 422 | 2640 | BandAverage(bands, i, [kelvin](double w) { return PlanckRadiance(w, kelvin); }); | |
| 423 | 80 | irradiance[static_cast<std::size_t>(i)] = | |
| 424 | 80 | static_cast<float>(totalIrradiance * radiance / totalRadiance); | |
| 425 | } | ||
| 426 | 2 | return irradiance; | |
| 427 | ✗ | } | |
| 428 | |||
| 429 | ✗ | std::vector<float> ThermalEmission(const SpectralBands& bands, double kelvin, | |
| 430 | const std::function<double(double)>& emissivity) { | ||
| 431 | ✗ | bands.Check(); | |
| 432 | ✗ | std::vector<float> radiance(static_cast<std::size_t>(bands.count), 0.0f); | |
| 433 | ✗ | if (!(kelvin > 0.0)) { | |
| 434 | return radiance; | ||
| 435 | } | ||
| 436 | ✗ | for (int i = 0; i < bands.count; ++i) { | |
| 437 | ✗ | radiance[static_cast<std::size_t>(i)] = static_cast<float>(BandAverage( | |
| 438 | ✗ | bands, i, [&](double w) { return emissivity(w) * PlanckRadiance(w, kelvin); })); | |
| 439 | } | ||
| 440 | return radiance; | ||
| 441 | ✗ | } | |
| 442 | |||
| 443 | // ---------------------------------------------------------------------------------------------- | ||
| 444 | // Colour | ||
| 445 | |||
| 446 | 180224 | glm::dvec3 CieXyz(double w) { | |
| 447 | 180224 | const double x = 1.056 * Lobe(w, 599.8, 37.9, 31.0) + 0.362 * Lobe(w, 442.0, 16.0, 26.7) - | |
| 448 | 180224 | 0.065 * Lobe(w, 501.1, 20.4, 26.2); | |
| 449 | 180224 | const double y = 0.821 * Lobe(w, 568.8, 46.9, 40.5) + 0.286 * Lobe(w, 530.9, 16.3, 31.1); | |
| 450 | 180224 | const double z = 1.217 * Lobe(w, 437.0, 11.8, 36.0) + 0.681 * Lobe(w, 459.0, 26.0, 13.8); | |
| 451 | 180224 | return {x, y, z}; | |
| 452 | } | ||
| 453 | |||
| 454 | 3 | std::vector<glm::dvec3> CieWeights(const SpectralBands& bands) { | |
| 455 | 3 | bands.Check(); | |
| 456 | 3 | std::vector<glm::dvec3> weights(static_cast<std::size_t>(bands.count)); | |
| 457 | 147 | for (int i = 0; i < bands.count; ++i) { | |
| 458 | 144 | weights[static_cast<std::size_t>(i)] = BandAverage(bands, i, CieXyz) * bands.Width(); | |
| 459 | } | ||
| 460 | 3 | return weights; | |
| 461 | ✗ | } | |
| 462 | |||
| 463 | 3 | std::vector<glm::dvec3> ColourWeights(const SpectralBands& bands) { | |
| 464 | 3 | if (bands.rgb) { | |
| 465 | 1 | return RgbModeColourWeights(); | |
| 466 | } | ||
| 467 | 2 | if (bands.count >= kFitBelow) { | |
| 468 | 1 | return CieWeights(bands); | |
| 469 | } | ||
| 470 | 1 | return CachedColourFit(bands); | |
| 471 | } | ||
| 472 | |||
| 473 | 2 | glm::dvec3 ToXyz(const SpectralBands& bands, std::span<const float> values) { | |
| 474 | 2 | const std::vector<glm::dvec3> weights = CieWeights(bands); | |
| 475 | 2 | if (values.size() != weights.size()) { | |
| 476 | 1 | throw std::invalid_argument("ToXyz needs one value per band"); | |
| 477 | } | ||
| 478 | 1 | glm::dvec3 xyz(0.0); | |
| 479 | 65 | for (std::size_t i = 0; i < values.size(); ++i) { | |
| 480 | 64 | xyz += static_cast<double>(values[i]) * weights[i]; | |
| 481 | } | ||
| 482 | 1 | return xyz; | |
| 483 | 2 | } | |
| 484 | |||
| 485 | 6870 | glm::dvec3 XyzToLinearSrgb(glm::dvec3 xyz) { | |
| 486 | 6870 | return {3.2404542 * xyz.x - 1.5371385 * xyz.y - 0.4985314 * xyz.z, | |
| 487 | 6870 | -0.9692660 * xyz.x + 1.8760108 * xyz.y + 0.0415560 * xyz.z, | |
| 488 | 6870 | 0.0556434 * xyz.x - 0.2040259 * xyz.y + 1.0572252 * xyz.z}; | |
| 489 | } | ||
| 490 | |||
| 491 | 6869 | glm::dvec3 EqualEnergyXyzToLinearSrgb(glm::dvec3 xyz) { | |
| 492 | 6869 | static const glm::dmat3 adaptation = EqualEnergyToD65(); | |
| 493 | 6869 | return XyzToLinearSrgb(adaptation * xyz); | |
| 494 | } | ||
| 495 | |||
| 496 | // ---------------------------------------------------------------------------------------------- | ||
| 497 | // RGB upsampling | ||
| 498 | |||
| 499 | 3275493 | double SigmoidSpectrum::operator()(double wavelengthNm) const { | |
| 500 | 3275493 | const double t = (wavelengthNm - kSigmoidMinNm) / (kSigmoidMaxNm - kSigmoidMinNm); | |
| 501 | 3275493 | return Sigmoid((coefficients.x * t + coefficients.y) * t + coefficients.z); | |
| 502 | } | ||
| 503 | |||
| 504 | 7 | Spectrum SigmoidSpectrum::Tabulate() const { | |
| 505 | 7 | std::vector<double> wavelengths; | |
| 506 | 7 | std::vector<double> values; | |
| 507 | 672 | for (double w = kSigmoidMinNm; w <= kSigmoidMaxNm; w += 5.0) { | |
| 508 | 665 | wavelengths.push_back(w); | |
| 509 | 665 | values.push_back((*this)(w)); | |
| 510 | } | ||
| 511 | 14 | return Spectrum(std::move(wavelengths), std::move(values)); | |
| 512 | 7 | } | |
| 513 | |||
| 514 | 6 | glm::dvec3 ReflectanceToLinearSrgb(const Spectrum& reflectance) { | |
| 515 | 2832 | return EqualEnergyXyzToLinearSrgb(EqualEnergyXyz([&](double w) { return reflectance(w); })); | |
| 516 | } | ||
| 517 | |||
| 518 | 145 | SigmoidSpectrum RgbToSigmoid(glm::dvec3 rgb) { | |
| 519 | 145 | rgb = glm::clamp(rgb, glm::dvec3(0.0), glm::dvec3(1.0)); | |
| 520 | // Start from the flat spectrum of the colour's mean, which is exact for greys. | ||
| 521 | 145 | const double mean = (rgb.r + rgb.g + rgb.b) / 3.0; | |
| 522 | 290 | SigmoidSpectrum spectrum{{0.0, 0.0, InverseSigmoid(std::clamp(mean, 1e-6, 1.0 - 1e-6))}}; | |
| 523 | 145 | if (rgb.r == rgb.g && rgb.g == rgb.b) { | |
| 524 | return spectrum; | ||
| 525 | } | ||
| 526 | 979 | for (int iteration = 0; iteration < 40; ++iteration) { | |
| 527 | 979 | const glm::dvec3 residual = SigmoidColour(spectrum) - rgb; | |
| 528 | 979 | const double error = glm::length(residual); | |
| 529 | 979 | if (error < 1e-7) { | |
| 530 | break; | ||
| 531 | } | ||
| 532 | 837 | const glm::dvec3 step = Solve(SigmoidColourJacobian(spectrum), residual); | |
| 533 | 837 | const std::optional<SigmoidSpectrum> next = DampedStep(spectrum, step, rgb, error); | |
| 534 | 837 | if (!next) { | |
| 535 | break; | ||
| 536 | } | ||
| 537 | 837 | spectrum = *next; | |
| 538 | } | ||
| 539 | return spectrum; | ||
| 540 | } | ||
| 541 | |||
| 542 | 2 | Spectrum RgbToSpectrum(glm::dvec3 rgb) { | |
| 543 | 2 | return RgbToSigmoid(rgb).Tabulate(); | |
| 544 | } | ||
| 545 | |||
| 546 | } // namespace ernie | ||
| 547 |