material/uplift.cpp - GCC Code Coverage Report

material/uplift.cpp

Line	Branch	Exec	Source
1			// SPDX-License-Identifier: BSD-3-Clause
2			#include "chrono_batch_render/material/uplift.h"
3
4			#include <algorithm>
5			#include <chrono>
6			#include <cmath>
7			#include <cstring>
8			#include <optional>
9			#include <string>
10			#include <thread>
11
12			#include <fmt/format.h>
13
14			#include "chrono_batch_render/util/log.h"
15
16			namespace chrono::batch_render {
17
18			namespace {
19
20			constexpr uint32_t kSamples = 95;
21			constexpr float kStep = (kSigmoidLastNm - kSigmoidFirstNm) / float(kSamples - 1);
22			constexpr float kWhiteKelvin = 6504;
23
24		665	float Lobe(float nanometers, float mean, float below, float above) {
25	2/2 ✓ Branch 0 taken 458 times. ✓ Branch 1 taken 207 times.	665	const float t = (nanometers - mean) / (nanometers < mean ? below : above);
26		665	return std::exp(-0.5f * t * t);
27			}
28
29		95	std::array<float, 3> Observer(float nanometers) {
30		95	const float x = 1.056f * Lobe(nanometers, 599.8f, 37.9f, 31.0f) +
31		95	0.362f * Lobe(nanometers, 442.0f, 16.0f, 26.7f) -
32		95	0.065f * Lobe(nanometers, 501.1f, 20.4f, 26.2f);
33		95	const float y = 0.821f * Lobe(nanometers, 568.8f, 46.9f, 40.5f) +
34		95	0.286f * Lobe(nanometers, 530.9f, 16.3f, 31.1f);
35		95	const float z = 1.217f * Lobe(nanometers, 437.0f, 11.8f, 36.0f) +
36		95	0.681f * Lobe(nanometers, 459.0f, 26.0f, 13.8f);
37		95	return {x, y, z};
38			}
39
40		95	float Planck(float nanometers, float kelvin) {
41		95	const double meters = double(nanometers) * 1e-9;
42		95	return float(1 / (std::pow(meters, 5) * (std::exp(1.438776877e-2 / (meters * kelvin)) - 1)));
43			}
44
45			struct Weights {
46			std::array<std::array<float, 3>, kSamples> rgb{};
47			std::array<float, kSamples> nanometers{};
48			};
49
50		2548169	const Weights& RgbWeights() {
51		1	static const Weights weights = [] {
52		1	Weights made;
53		1	std::array<std::array<float, 3>, kSamples> xyz{};
54		1	std::array<double, 3> white{};
55	2/2 ✓ Branch 0 taken 95 times. ✓ Branch 1 taken 1 time.	96	for (uint32_t i = 0; i < kSamples; ++i) {
56		95	const float nanometers = kSigmoidFirstNm + float(i) * kStep;
57		95	made.nanometers[i] = nanometers;
58		95	const std::array<float, 3> observer = Observer(nanometers);
59		95	const float light = Planck(nanometers, kWhiteKelvin);
60	2/2 ✓ Branch 1 taken 285 times. ✓ Branch 2 taken 95 times.	380	for (int c = 0; c < 3; ++c) {
61		285	xyz[i][c] = observer[c] * light;
62			}
63			}
64		1	constexpr float kToRgb[3][3] = {{3.2404542f, -1.5371385f, -0.4985314f},
65			{-0.9692660f, 1.8760108f, 0.0415560f},
66			{0.0556434f, -0.2040259f, 1.0572252f}};
67	2/2 ✓ Branch 0 taken 95 times. ✓ Branch 1 taken 1 time.	96	for (uint32_t i = 0; i < kSamples; ++i) {
68	2/2 ✓ Branch 0 taken 285 times. ✓ Branch 1 taken 95 times.	380	for (int c = 0; c < 3; ++c) {
69		285	made.rgb[i][c] =
70		285	kToRgb[c][0] * xyz[i][0] + kToRgb[c][1] * xyz[i][1] + kToRgb[c][2] * xyz[i][2];
71		285	white[c] += made.rgb[i][c];
72			}
73			}
74	2/2 ✓ Branch 0 taken 95 times. ✓ Branch 1 taken 1 time.	96	for (uint32_t i = 0; i < kSamples; ++i) {
75	2/2 ✓ Branch 0 taken 285 times. ✓ Branch 1 taken 95 times.	380	for (int c = 0; c < 3; ++c) {
76		285	made.rgb[i][c] = float(made.rgb[i][c] / white[c]);
77			}
78			}
79		1	return made;
80	3/4 ✓ Branch 0 taken 1 time. ✓ Branch 1 taken 2548168 times. ✓ Branch 3 taken 1 time. ✗ Branch 4 not taken.	2548169	}();
81		2548169	return weights;
82			}
83
84		242076059	float Sigmoid(float x) {
85	1/4 ✗ Branch 0 not taken. ✓ Branch 1 taken 242076059 times. ✗ Branch 2 not taken. ✗ Branch 3 not taken.	242076059	return std::isinf(x) ? (x > 0 ? 1.0f : 0.0f) : 0.5f + x / (2 * std::sqrt(1 + x * x));
86			}
87
88		242076059	float Normalized(float nanometers) {
89		242076059	return (nanometers - kSigmoidFirstNm) / (kSigmoidLastNm - kSigmoidFirstNm);
90			}
91
92		2548169	std::array<float, 3> Residual(const std::array<float, 3>& c, const std::array<float, 3>& rgb) {
93		2548169	const Weights& weights = RgbWeights();
94		2548169	std::array<float, 3> sum{};
95	2/2 ✓ Branch 0 taken 242076055 times. ✓ Branch 1 taken 2548169 times.	244624224	for (uint32_t i = 0; i < kSamples; ++i) {
96		242076055	const float value = SigmoidValue(c, weights.nanometers[i]);
97	2/2 ✓ Branch 1 taken 726228165 times. ✓ Branch 2 taken 242076055 times.	968304220	for (int k = 0; k < 3; ++k) {
98		726228165	sum[k] += value * weights.rgb[i][k];
99			}
100			}
101		2548169	return {sum[0] - rgb[0], sum[1] - rgb[1], sum[2] - rgb[2]};
102			}
103
104			template <size_t N>
105		530875	std::optional<std::array<double, N>> Solve(std::array<std::array<double, N + 1>, N> m) {
106	2/4 std::optional<std::array<double, 3ul> > chrono::batch_render::(anonymous namespace)::Solve<3ul>(std::array<std::array<double, (3ul)+(1)>, 3ul>): ✓ Branch 0 taken 1592481 times. ✓ Branch 1 taken 530782 times. std::optional<std::array<double, 4ul> > chrono::batch_render::(anonymous namespace)::Solve<4ul>(std::array<std::array<double, (4ul)+(1)>, 4ul>): ✗ Branch 0 not taken. ✗ Branch 1 not taken.	2123263	for (size_t pivot = 0; pivot < N; ++pivot) {
107		1592481	size_t best = pivot;
108	2/4 std::optional<std::array<double, 3ul> > chrono::batch_render::(anonymous namespace)::Solve<3ul>(std::array<std::array<double, (3ul)+(1)>, 3ul>): ✓ Branch 0 taken 1592565 times. ✓ Branch 1 taken 1592481 times. std::optional<std::array<double, 4ul> > chrono::batch_render::(anonymous namespace)::Solve<4ul>(std::array<std::array<double, (4ul)+(1)>, 4ul>): ✗ Branch 0 not taken. ✗ Branch 1 not taken.	3185046	for (size_t r = pivot + 1; r < N; ++r) {
109	2/4 std::optional<std::array<double, 3ul> > chrono::batch_render::(anonymous namespace)::Solve<3ul>(std::array<std::array<double, (3ul)+(1)>, 3ul>): ✓ Branch 0 taken 1240914 times. ✓ Branch 1 taken 351651 times. std::optional<std::array<double, 4ul> > chrono::batch_render::(anonymous namespace)::Solve<4ul>(std::array<std::array<double, (4ul)+(1)>, 4ul>): ✗ Branch 0 not taken. ✗ Branch 1 not taken.	2833479	best = std::abs(m[r][pivot]) > std::abs(m[best][pivot]) ? r : best;
110			}
111	2/4 std::optional<std::array<double, 3ul> > chrono::batch_render::(anonymous namespace)::Solve<3ul>(std::array<std::array<double, (3ul)+(1)>, 3ul>): ✓ Branch 0 taken 1592388 times. ✓ Branch 1 taken 93 times. std::optional<std::array<double, 4ul> > chrono::batch_render::(anonymous namespace)::Solve<4ul>(std::array<std::array<double, (4ul)+(1)>, 4ul>): ✗ Branch 0 not taken. ✗ Branch 1 not taken.	3184962	std::swap(m[pivot], m[best]);
112	2/4 std::optional<std::array<double, 3ul> > chrono::batch_render::(anonymous namespace)::Solve<3ul>(std::array<std::array<double, (3ul)+(1)>, 3ul>): ✓ Branch 0 taken 1592388 times. ✓ Branch 1 taken 93 times. std::optional<std::array<double, 4ul> > chrono::batch_render::(anonymous namespace)::Solve<4ul>(std::array<std::array<double, (4ul)+(1)>, 4ul>): ✗ Branch 0 not taken. ✗ Branch 1 not taken.	1592481	if (std::abs(m[pivot][pivot]) < 1e-30) {
113		93	return std::nullopt;
114			}
115	2/4 std::optional<std::array<double, 3ul> > chrono::batch_render::(anonymous namespace)::Solve<3ul>(std::array<std::array<double, (3ul)+(1)>, 3ul>): ✓ Branch 0 taken 4777164 times. ✓ Branch 1 taken 1592388 times. std::optional<std::array<double, 4ul> > chrono::batch_render::(anonymous namespace)::Solve<4ul>(std::array<std::array<double, (4ul)+(1)>, 4ul>): ✗ Branch 0 not taken. ✗ Branch 1 not taken.	6369552	for (size_t r = 0; r < N; ++r) {
116	2/4 std::optional<std::array<double, 3ul> > chrono::batch_render::(anonymous namespace)::Solve<3ul>(std::array<std::array<double, (3ul)+(1)>, 3ul>): ✓ Branch 0 taken 1592388 times. ✓ Branch 1 taken 3184776 times. std::optional<std::array<double, 4ul> > chrono::batch_render::(anonymous namespace)::Solve<4ul>(std::array<std::array<double, (4ul)+(1)>, 4ul>): ✗ Branch 0 not taken. ✗ Branch 1 not taken.	4777164	if (r == pivot) {
117		1592388	continue;
118			}
119		3184776	const double factor = m[r][pivot] / m[pivot][pivot];
120	2/4 std::optional<std::array<double, 3ul> > chrono::batch_render::(anonymous namespace)::Solve<3ul>(std::array<std::array<double, (3ul)+(1)>, 3ul>): ✓ Branch 0 taken 9554394 times. ✓ Branch 1 taken 3184776 times. std::optional<std::array<double, 4ul> > chrono::batch_render::(anonymous namespace)::Solve<4ul>(std::array<std::array<double, (4ul)+(1)>, 4ul>): ✗ Branch 0 not taken. ✗ Branch 1 not taken.	12739170	for (size_t c = pivot; c <= N; ++c) {
121		9554394	m[r][c] -= factor * m[pivot][c];
122			}
123			}
124			}
125		530782	std::array<double, N> x{};
126	2/4 std::optional<std::array<double, 3ul> > chrono::batch_render::(anonymous namespace)::Solve<3ul>(std::array<std::array<double, (3ul)+(1)>, 3ul>): ✓ Branch 0 taken 1592346 times. ✓ Branch 1 taken 530782 times. std::optional<std::array<double, 4ul> > chrono::batch_render::(anonymous namespace)::Solve<4ul>(std::array<std::array<double, (4ul)+(1)>, 4ul>): ✗ Branch 0 not taken. ✗ Branch 1 not taken.	2123128	for (size_t r = 0; r < N; ++r) {
127		1592346	x[r] = m[r][N] / m[r][r];
128			}
129		530782	return x;
130			}
131
132		98304	size_t TableIndex(uint32_t l, uint32_t zi, uint32_t y, uint32_t x) {
133		98304	constexpr size_t n = kSigmoidResolution;
134		98304	return (((l * n + zi) * n + y) * n + x) * 3;
135			}
136
137		98304	std::array<float, 3> TableColor(uint32_t l, float z, uint32_t y, uint32_t x) {
138		98304	constexpr float last = float(kSigmoidResolution - 1);
139		98304	std::array<float, 3> rgb{};
140		98304	rgb[l] = z;
141		98304	rgb[(l + 1) % 3] = float(x) / last * z;
142		98304	rgb[(l + 2) % 3] = float(y) / last * z;
143		98304	return rgb;
144			}
145
146			}
147
148		98308	std::array<float, 3> FitSigmoid(std::array<float, 3> rgb, std::array<float, 3> c) {
149	6/6 ✓ Branch 0 taken 5056 times. ✓ Branch 1 taken 93252 times. ✓ Branch 2 taken 3103 times. ✓ Branch 3 taken 1953 times. ✓ Branch 4 taken 3072 times. ✓ Branch 5 taken 31 times.	98308	if (rgb[0] <= 0 && rgb[1] <= 0 && rgb[2] <= 0) {
150		3072	return {0, 0, -1e4f};
151			}
152		1563225	const auto norm = [](const std::array<float, 3>& r) {
153		1467989	return r[0] * r[0] + r[1] * r[1] + r[2] * r[2];
154			};
155		95236	std::array<float, 3> r = Residual(c, rgb);
156	4/4 ✓ Branch 1 taken 607684 times. ✓ Branch 2 taken 8 times. ✓ Branch 3 taken 530875 times. ✓ Branch 4 taken 76809 times.	607692	for (int iteration = 0; iteration < 50 && norm(r) > 1e-14f; ++iteration) {
157		530875	std::array<std::array<double, 4>, 3> system{};
158	2/2 ✓ Branch 0 taken 1592625 times. ✓ Branch 1 taken 530875 times.	2123500	for (int k = 0; k < 3; ++k) {
159		1592625	std::array<float, 3> moved = c;
160	2/2 ✓ Branch 0 taken 1456442 times. ✓ Branch 1 taken 136183 times.	1592625	const float delta = 1e-3f * std::max(1.0f, std::abs(c[k]));
161		1592625	moved[k] += delta;
162		1592625	const std::array<float, 3> shifted = Residual(moved, rgb);
163	2/2 ✓ Branch 1 taken 4777875 times. ✓ Branch 2 taken 1592625 times.	6370500	for (int row = 0; row < 3; ++row) {
164		4777875	system[row][k] = (shifted[row] - r[row]) / delta;
165		4777875	system[row][3] = r[row];
166			}
167			}
168		530875	const std::optional<std::array<double, 3>> solved = Solve<3>(system);
169	2/2 ✓ Branch 0 taken 530782 times. ✓ Branch 1 taken 93 times.	530875	if (!solved) {
170			break;
171			}
172		530782	const std::array<float, 3> step = {float((solved)[0]), float((solved)[1]),
173		530782	float((*solved)[2])};
174		530782	bool improved = false;
175	2/2 ✓ Branch 0 taken 860305 times. ✓ Branch 1 taken 18326 times.	878631	for (float scale = 1; scale > 1e-4f; scale *= 0.5f) {
176		860305	const std::array<float, 3> tried = {c[0] - scale * step[0], c[1] - scale * step[1],
177		860305	c[2] - scale * step[2]};
178		860305	const std::array<float, 3> residual = Residual(tried, rgb);
179	2/2 ✓ Branch 0 taken 512456 times. ✓ Branch 1 taken 347849 times.	860305	if (norm(residual) < norm(r)) {
180		512456	c = tried;
181		512456	r = residual;
182		512456	improved = true;
183		512456	break;
184			}
185			}
186		512456	if (!improved) {
187			break;
188			}
189			}
190		95236	return c;
191			}
192
193		242076059	float SigmoidValue(const std::array<float, 3>& c, float nanometers) {
194		242076059	const float t = Normalized(nanometers);
195		242076059	return Sigmoid((c[0] * t + c[1]) * t + c[2]);
196			}
197
198		3	std::array<float, 3> SigmoidRgb(const std::array<float, 3>& c) {
199		3	return Residual(c, {0, 0, 0});
200			}
201
202		✗	std::vector<std::array<float, 3>> LinearUplift() {
203		✗	const SigmoidTable& table = Sigmoids();
204		✗	constexpr uint32_t n = kSigmoidResolution;
205		✗	std::vector<std::array<float, 3>> colours;
206		✗	std::vector<std::array<float, 3>> coefficients;
207	0/2 ✗ Branch 0 not taken. ✗ Branch 1 not taken.	✗	for (uint32_t l = 0; l < 3; ++l) {
208	0/2 ✗ Branch 0 not taken. ✗ Branch 1 not taken.	✗	for (uint32_t zi = 0; zi < n; ++zi) {
209	0/2 ✗ Branch 0 not taken. ✗ Branch 1 not taken.	✗	for (uint32_t y = 0; y < n; ++y) {
210	0/2 ✗ Branch 0 not taken. ✗ Branch 1 not taken.	✗	for (uint32_t x = 0; x < n; ++x) {
211		✗	const size_t at = TableIndex(l, zi, y, x);
212	0/1 ✗ Branch 1 not taken.	✗	colours.push_back(TableColor(l, table.scales[zi], y, x));
213	0/1 ✗ Branch 1 not taken.	✗	coefficients.push_back({table.coefficients[at], table.coefficients[at + 1],
214	0/1 ✗ Branch 1 not taken.	✗	table.coefficients[at + 2]});
215			}
216			}
217			}
218			}
219		✗	std::vector<std::array<float, 3>> weights;
220	0/2 ✗ Branch 0 not taken. ✗ Branch 1 not taken.	✗	for (uint32_t k = 0; k < kLinearUpliftSamples; ++k) {
221		✗	const float nanometers = kLinearUpliftFirstNm + float(k) * kLinearUpliftStepNm;
222		✗	std::array<std::array<double, 5>, 4> m{};
223	0/2 ✗ Branch 0 not taken. ✗ Branch 1 not taken.	✗	for (size_t i = 0; i < colours.size(); ++i) {
224		✗	const double value = SigmoidValue(coefficients[i], nanometers);
225	0/2 ✗ Branch 0 not taken. ✗ Branch 1 not taken.	✗	for (int r = 0; r < 3; ++r) {
226	0/2 ✗ Branch 0 not taken. ✗ Branch 1 not taken.	✗	for (int c = 0; c < 3; ++c) {
227		✗	m[r][c] += 2.0 * colours[i][r] * colours[i][c];
228			}
229		✗	m[r][4] += 2.0 * colours[i][r] * value;
230			}
231			}
232	0/2 ✗ Branch 0 not taken. ✗ Branch 1 not taken.	✗	for (int r = 0; r < 3; ++r) {
233		✗	m[r][3] = 1;
234		✗	m[3][r] = 1;
235			}
236		✗	m[3][4] = 1;
237	0/2 ✗ Branch 1 not taken. ✗ Branch 2 not taken.	✗	const std::array<double, 4> solved = Solve<4>(m).value_or(std::array<double, 4>{});
238	0/1 ✗ Branch 1 not taken.	✗	weights.push_back({float(solved[0]), float(solved[1]), float(solved[2])});
239			}
240		✗	return weights;
241		✗	}
242
243			namespace {
244
245		1	SigmoidTable FitSigmoidTable() {
246		1	const auto start = std::chrono::steady_clock::now();
247		1	constexpr uint32_t n = kSigmoidResolution;
248		1	SigmoidTable made;
249	1/1 ✓ Branch 1 taken 1 time.	1	made.scales.resize(n);
250		32	const auto smoothstep = [](float x) { return x * x * (3 - 2 * x); };
251	2/2 ✓ Branch 0 taken 32 times. ✓ Branch 1 taken 1 time.	33	for (uint32_t k = 0; k < n; ++k) {
252		32	made.scales[k] = smoothstep(smoothstep(float(k) / float(n - 1)));
253			}
254	1/1 ✓ Branch 1 taken 1 time.	1	made.coefficients.assign(size_t(3) * n * n * n * 3, 0);
255		1	std::vector<std::thread> workers;
256	2/2 ✓ Branch 0 taken 3 times. ✓ Branch 1 taken 1 time.	4	for (uint32_t l = 0; l < 3; ++l) {
257	1/1 ✓ Branch 1 taken 3 times.	3	workers.emplace_back([&made, l] {
258	2/2 ✓ Branch 0 taken 96 times. ✓ Branch 1 taken 3 times.	99	for (uint32_t y = 0; y < n; ++y) {
259	2/2 ✓ Branch 0 taken 3072 times. ✓ Branch 1 taken 96 times.	3168	for (uint32_t x = 0; x < n; ++x) {
260		101376	const auto fit = [&](uint32_t zi, std::array<float, 3> c) {
261		98304	c = FitSigmoid(TableColor(l, made.scales[zi], y, x), c);
262		196608	std::copy(c.begin(), c.end(),
263		98304	made.coefficients.begin() + ptrdiff_t(TableIndex(l, zi, y, x)));
264		98304	return c;
265		3072	};
266		3072	const uint32_t middle = n / 5;
267		3072	std::array<float, 3> c{};
268	2/2 ✓ Branch 0 taken 79872 times. ✓ Branch 1 taken 3072 times.	82944	for (uint32_t zi = middle; zi < n; ++zi) {
269		79872	c = fit(zi, c);
270			}
271		3072	c = {};
272	2/2 ✓ Branch 0 taken 18432 times. ✓ Branch 1 taken 3072 times.	21504	for (uint32_t zi = middle; zi-- > 0;) {
273		18432	c = fit(zi, c);
274			}
275			}
276			}
277		3	});
278			}
279	2/2 ✓ Branch 0 taken 3 times. ✓ Branch 1 taken 1 time.	4	for (std::thread& worker : workers) {
280	1/1 ✓ Branch 1 taken 3 times.	3	worker.join();
281			}
282	1/1 ✓ Branch 1 taken 1 time.	1	Log(LogLevel::Debug, "fitted the {}-cubed RGB sigmoid uplift table in {:.1f} s", n,
283	1/1 ✓ Branch 3 taken 1 time.	1	std::chrono::duration<double>(std::chrono::steady_clock::now() - start).count());
284		1	return made;
285		1	}
286
287			}
288
289		8	SigmoidTable CachedSigmoids(const DiskCache& cache) {
290		8	constexpr size_t n = kSigmoidResolution;
291		8	constexpr size_t scaleBytes = n * sizeof(float);
292		8	constexpr size_t coefficientBytes = 3 * n * n * n * 3 * sizeof(float);
293		8	const std::string key =
294		8	fmt::format("sigmoid uplift table {}", CHRONO_BATCH_RENDER_UPLIFT_DIGEST);
295	1/1 ✓ Branch 1 taken 8 times.	8	if (const std::optional<std::string> bytes = cache.Load(key);
296	3/4 ✓ Branch 0 taken 7 times. ✓ Branch 1 taken 1 time. ✓ Branch 2 taken 7 times. ✗ Branch 3 not taken.	8	bytes && bytes->size() == scaleBytes + coefficientBytes) {
297		7	SigmoidTable loaded;
298	1/1 ✓ Branch 1 taken 7 times.	7	loaded.scales.resize(n);
299	1/1 ✓ Branch 1 taken 7 times.	7	loaded.coefficients.resize(coefficientBytes / sizeof(float));
300		7	std::memcpy(loaded.scales.data(), bytes->data(), scaleBytes);
301		7	std::memcpy(loaded.coefficients.data(), bytes->data() + scaleBytes, coefficientBytes);
302		7	return loaded;
303		7	}
304	1/1 ✓ Branch 1 taken 1 time.	1	SigmoidTable made = FitSigmoidTable();
305	2/2 ✓ Branch 1 taken 1 time. ✓ Branch 4 taken 1 time.	1	std::string bytes(scaleBytes + coefficientBytes, '\0');
306	1/1 ✓ Branch 1 taken 1 time.	1	std::memcpy(bytes.data(), made.scales.data(), scaleBytes);
307		1	std::memcpy(bytes.data() + scaleBytes, made.coefficients.data(), coefficientBytes);
308	1/1 ✓ Branch 1 taken 1 time.	1	cache.Store(key, bytes);
309		1	return made;
310		10	}
311
312		275	const SigmoidTable& Sigmoids() {
313	7/8 ✓ Branch 0 taken 6 times. ✓ Branch 1 taken 269 times. ✓ Branch 3 taken 6 times. ✗ Branch 4 not taken. ✓ Branch 6 taken 6 times. ✓ Branch 9 taken 6 times. ✓ Branch 12 taken 6 times. ✓ Branch 15 taken 6 times.	305	static const SigmoidTable table = CachedSigmoids(DiskCache(UserCacheDirectory() / "uplift"));
314		275	return table;
315			}
316			}
317