physics/spectrum.h - GCC Code Coverage Report

physics/spectrum.h

Line	Branch	Exec	Source
1			// SPDX-License-Identifier: BSD-3-Clause
2			/// @file spectrum.h
3			/// @brief Spectra, spectral bands, black bodies and colour matching.
4			/// Wavelengths are in nm, irradiance in W m^-2 nm^-1 and radiance in
5			/// W m^-2 sr^-1 nm^-1.
6			#pragma once
7
8			#include "ernie/physics/constants.h"
9
10			#include <glm/glm.hpp>
11
12			#include <cstddef>
13			#include <functional>
14			#include <span>
15			#include <utility>
16			#include <vector>
17
18			namespace ernie {
19
20			/// @brief A camera's spectral sampling: a range split into equal
21			/// bands.
22			/// Band i covers [Lower(i), Lower(i) + Width()).
23			struct SpectralBands {
24			/// Lower end of the range, in nm.
25			double minNm = 380.0;
26			/// Upper end of the range, in nm.
27			double maxNm = 780.0;
28			/// Number of bands.
29			int count = 16;
30			/// RGB mode: three bands carry linear sRGB blue, green and red,
31			/// so colours match a 3-channel renderer's instead of shifting.
32			bool rgb = false;
33
34			/// @return RGB mode's bands: 400-700 nm in three.
35		1	static SpectralBands Rgb() {
36		1	return {400.0, 700.0, 3, true};
37			}
38
39			/// @return A band's width, in nm.
40		7977	double Width() const {
41	3/6 <unknown function 35>: ✓ Branch 1 taken 1 time. ✗ Branch 2 not taken. <unknown function 40>: ✓ Branch 0 taken 3 times. ✗ Branch 1 not taken. ✓ Branch 1 taken 144 times. ✗ Branch 2 not taken.	160	return (maxNm - minNm) / count;
42			}
43
44			/// @return Band `band`'s lower edge, in nm.
45		7829	double Lower(int band) const {
46	2/4 ✓ Branch 1 taken 150 times. ✗ Branch 2 not taken. ✓ Branch 4 taken 2 times. ✗ Branch 5 not taken.	7333	return minNm + band * Width();
47			}
48
49			/// @return Band `band`'s centre, in nm.
50		508	double Centre(int band) const {
51	1/2 ✓ Branch 1 taken 114 times. ✗ Branch 2 not taken.	508	return Lower(band) + 0.5 * Width();
52			}
53
54			/// @brief Validates the bands.
55			/// @throws std::invalid_argument Unless 0 < minNm < maxNm,
56			/// 1 <= count <= kMaxBands, and count is 3 in RGB mode.
57			void Check() const;
58
59			/// The most bands a camera may have.
60			static constexpr int kMaxBands = 64;
61
62			bool operator==(const SpectralBands&) const = default;
63			};
64
65			/// @brief A spectrum's value in one band.
66			/// @param bands The bands.
67			/// @param band The band index.
68			/// @param spectrum The spectrum, by wavelength in nm.
69			/// @return Its band average; in RGB mode, the band's linear sRGB
70			/// channel relative to equal-energy white, at least 0.
71			double BandValue(const SpectralBands& bands, int band,
72			const std::function<double(double)>& spectrum);
73
74			/// @brief A tabulated spectrum, linear between samples.
75			/// Constant beyond the ends; an empty spectrum is zero.
76	4/8 ✓ Branch 3 taken 5 times. ✗ Branch 4 not taken. ✗ Branch 23 not taken. ✓ Branch 24 taken 1 time. ✗ Branch 27 not taken. ✓ Branch 28 taken 1 time. ✗ Branch 31 not taken. ✓ Branch 32 taken 1 time.	149	class Spectrum {
77			public:
78	1/2 ✓ Branch 1 taken 1 time. ✗ Branch 2 not taken.	1	Spectrum() = default;
79
80			/// @brief Tabulates a spectrum.
81			/// @param wavelengthsNm Strictly increasing wavelengths, in nm.
82			/// @param values The value at each.
83			/// @throws std::invalid_argument If the wavelengths are not
84			/// strictly increasing or the sizes differ.
85			Spectrum(std::vector<double> wavelengthsNm, std::vector<double> values);
86
87			/// @param value The value.
88			/// @return A spectrum equal to `value` everywhere.
89			static Spectrum Constant(double value);
90
91			/// @brief Builds from interleaved (wavelength, value) pairs, the
92			/// form of `ernie:` spectrum attributes.
93			/// @param pairs The pairs.
94			/// @return The spectrum.
95			/// @throws std::invalid_argument For an odd count, or as the
96			/// constructor.
97			static Spectrum FromPairs(std::span<const float> pairs);
98
99			/// @param wavelengthNm The wavelength, in nm.
100			/// @return The interpolated value.
101			double operator()(double wavelengthNm) const;
102
103			/// @brief Exact mean of the interpolant over a range.
104			/// @param lowNm The range's start, in nm.
105			/// @param highNm The range's end, in nm.
106			/// @return The mean.
107			double Average(double lowNm, double highNm) const;
108
109			/// @param bands The bands.
110			/// @return The average over each band.
111			std::vector<float> Bin(const SpectralBands& bands) const;
112
113			/// @return Whether it has no samples.
114			bool Empty() const {
115			return wavelengths_.empty();
116			}
117
118			/// @return The sample wavelengths, in nm.
119		21	const std::vector<double>& Wavelengths() const {
120	2/4 ✓ Branch 1 taken 19 times. ✗ Branch 2 not taken. ✓ Branch 4 taken 19 times. ✗ Branch 5 not taken.	41	return wavelengths_;
121			}
122
123			/// @return The sample values.
124			const std::vector<double>& Values() const {
125		19	return values_;
126			}
127
128			private:
129			double Integral(double lowNm, double highNm) const;
130
131			std::vector<double> wavelengths_;
132			std::vector<double> values_;
133			};
134
135			/// @brief Planck's law.
136			/// @param wavelengthNm The wavelength, in nm.
137			/// @param kelvin The temperature, in kelvin.
138			/// @return Black-body spectral radiance, in W m^-2 sr^-1 nm^-1.
139			double PlanckRadiance(double wavelengthNm, double kelvin);
140
141			/// @brief CIE standard illuminant D (daylight).
142			/// Interpolated from the CIE's 10 nm basis over 300-830 nm; 0
143			/// outside.
144			/// @param wavelengthNm The wavelength, in nm.
145			/// @param kelvin Correlated colour temperature, 4000 to 25000 K
146			/// (D65 is 6504 K).
147			/// @return Relative spectral power, 100 at 560 nm.
148			/// @throws std::invalid_argument For a temperature out of range.
149			///
150			/// References:
151			/// - CIE 015:2018, "Colorimetry", 4th edition, section 4.1.2 and
152			/// table 6 (the S0, S1 and S2 components of daylight).
153			double CieDaylight(double wavelengthNm, double kelvin);
154
155			/// @brief A black body's irradiance averaged over each band.
156			/// @param bands The bands.
157			/// @param kelvin The temperature, in kelvin.
158			/// @param totalIrradiance Its integral over all wavelengths, in
159			/// W/m^2.
160			/// @return Spectral irradiance per band, in W m^-2 nm^-1.
161			std::vector<float> BlackbodyIrradiance(const SpectralBands& bands, double kelvin,
162			double totalIrradiance);
163
164			/// @brief Thermal emission of a surface, by Kirchhoff's law.
165			/// @param bands The bands.
166			/// @param kelvin The temperature; at or below 0 K gives zeros.
167			/// @param emissivity Emissivity by wavelength in nm.
168			/// @return Emissivity times Planck radiance, averaged per band, in
169			/// W m^-2 sr^-1 nm^-1.
170			///
171			/// References:
172			/// - M. Planck, "The Theory of Heat Radiation", 1914 (Kirchhoff's
173			/// law and the black-body spectrum).
174			std::vector<float> ThermalEmission(const SpectralBands& bands, double kelvin,
175			const std::function<double(double)>& emissivity);
176
177			/// @brief CIE 1931 2-degree colour matching functions, by Wyman et
178			/// al.'s multi-lobe Gaussian fit.
179			/// @param wavelengthNm The wavelength, in nm.
180			/// @return x-bar, y-bar and z-bar; y-bar peaks at 1 near 555 nm.
181			///
182			/// References:
183			/// - C. Wyman, P.-P. Sloan and P. Shirley, "Simple Analytic
184			/// Approximations to the CIE XYZ Color Matching Functions", Journal
185			/// of Computer Graphics Techniques 2(2):1-11, 2013.
186			glm::dvec3 CieXyz(double wavelengthNm);
187
188			/// @brief Colour matching functions integrated over each band.
189			/// XYZ = sum_i value_i * weights_i integrates a spectrum.
190			/// @param bands The bands.
191			/// @return One XYZ weight per band.
192			std::vector<glm::dvec3> CieWeights(const SpectralBands& bands);
193
194			/// @brief The weights a renderer turns band values into XYZ with.
195			/// `CieWeights` for 16 bands or more. For fewer, a least-squares fit
196			/// over sigmoid reflectances sampled at band centres, greys exact;
197			/// plain averages would shift colours (green reads yellow with 3).
198			/// @param bands The bands.
199			/// @return One XYZ weight per band.
200			///
201			/// References:
202			/// - W. Jakob and J. Hanika, "A Low-Dimensional Function Space for
203			/// Efficient Spectral Upsampling", Computer Graphics Forum
204			/// 38(2):147-155, 2019 (the reflectances fitted over).
205			std::vector<glm::dvec3> ColourWeights(const SpectralBands& bands);
206
207			/// @brief CIE XYZ of band values, without the luminous efficacy.
208			/// @param bands The bands.
209			/// @param values One value per band.
210			/// @return The sum of value times `CieWeights`.
211			glm::dvec3 ToXyz(const SpectralBands& bands, std::span<const float> values);
212
213			/// @brief Converts CIE XYZ to linear sRGB (D65, IEC 61966-2-1).
214			/// @param xyz The colour.
215			/// @return Linear sRGB.
216			glm::dvec3 XyzToLinearSrgb(glm::dvec3 xyz);
217
218			/// @brief Converts XYZ relative to illuminant E to linear sRGB,
219			/// Bradford-adapted to D65 so flat spectra map to greys.
220			/// @param xyz The colour, white at a flat spectrum's XYZ.
221			/// @return Linear sRGB.
222			///
223			/// References:
224			/// - K. M. Lam, "Metamerism and Colour Constancy", PhD thesis,
225			/// University of Bradford, 1985; matrices as given by B. Lindbloom,
226			/// http://www.brucelindbloom.com/Eqn_ChromAdapt.html
227			glm::dvec3 EqualEnergyXyzToLinearSrgb(glm::dvec3 xyz);
228
229			/// @brief A smooth reflectance in Jakob and Hanika's
230			/// sigmoid-polynomial form.
231			/// s = 1/2 + x / (2 sqrt(1 + x^2)), x = c0 t^2 + c1 t + c2, with t
232			/// the wavelength mapped from 360-830 nm to 0-1. Always in [0, 1].
233			///
234			/// References:
235			/// - W. Jakob and J. Hanika, "A Low-Dimensional Function Space for
236			/// Efficient Spectral Upsampling", Computer Graphics Forum (Proc.
237			/// Eurographics) 38(2):147-155, 2019.
238			/// https://doi.org/10.1111/cgf.13626
239			struct SigmoidSpectrum {
240			/// c0, c1 and c2; flat when c0 = c1 = 0.
241			glm::dvec3 coefficients{0.0};
242
243			/// @param wavelengthNm The wavelength, in nm.
244			/// @return The reflectance.
245			double operator()(double wavelengthNm) const;
246
247			/// @return The spectrum every 5 nm over 360-830 nm.
248			Spectrum Tabulate() const;
249			};
250
251			/// @brief Upsamples a linear sRGB colour to a sigmoid reflectance.
252			/// Gauss-Newton as Jakob and Hanika, but minimising in linear sRGB
253			/// rather than CIELAB. Greys give flat spectra exactly; unreachable
254			/// colours end at the nearest spectrum found.
255			/// @param rgb The colour, each component clamped to [0, 1].
256			/// @return The spectrum whose colour under E, adapted to D65, is
257			/// `rgb`.
258			SigmoidSpectrum RgbToSigmoid(glm::dvec3 rgb);
259
260			/// @brief A reflectance's colour under illuminant E, adapted to D65.
261			/// @param reflectance The spectrum.
262			/// @return Linear sRGB.
263			glm::dvec3 ReflectanceToLinearSrgb(const Spectrum& reflectance);
264
265			/// @brief `RgbToSigmoid` tabulated as a `Spectrum`.
266			/// @param rgb The colour.
267			/// @return The spectrum.
268			Spectrum RgbToSpectrum(glm::dvec3 rgb);
269
270			} // namespace ernie
271