Tatooine
direct_isosurface.h
Go to the documentation of this file.
1
3#ifndef TATOOINE_RENDERING_DIRECT_ISOSURFACE_H
4#define TATOOINE_RENDERING_DIRECT_ISOSURFACE_H
5//==============================================================================
6#include <tatooine/demangling.h>
7#include <tatooine/field.h>
8#include <tatooine/rectilinear_grid.h>
9#include <tatooine/rendering/camera.h>
10//==============================================================================
11namespace tatooine::rendering {
12//==============================================================================
22template <camera Camera, arithmetic IsoReal, typename Dim0, typename Dim1,
23 typename Dim2, typename Field, typename Shader>
24auto direct_isosurface(Camera const& cam,
25 rectilinear_grid<Dim0, Dim1, Dim2> const& g,
26 Field&& field, std::vector<IsoReal> const& isovalues,
27 Shader&& shader)
28 //
29 requires invocable<
30 Shader, vec<typename rectilinear_grid<Dim0, Dim1, Dim2>::real_type, 3>,
31 IsoReal, vec<typename rectilinear_grid<Dim0, Dim1, Dim2>::real_type, 3>,
32 vec<typename Camera::real_type, 3>, vec<std::size_t, 2>>
33//
34{
35 using grid_real_type = typename rectilinear_grid<Dim0, Dim1, Dim2>::real_type;
36 using pos_type = vec<grid_real_type, 3>;
37 constexpr auto use_indices =
38 std::is_invocable_v<Field, std::size_t, std::size_t, std::size_t>;
39 // using value_type =
40 // std::conditional_t<use_indices,
41 // std::invoke_result_t<Field, std::size_t, std::size_t,
42 // std::size_t>, std::invoke_result_t<Field, pos_type>>;
43 // static_assert(is_floating_point<value_type>);
44 using cam_real_type = typename Camera::real_type;
45 using viewdir_type = vec<cam_real_type, 3>;
46 using color_type = invoke_result<Shader, pos_type, IsoReal, pos_type,
47 viewdir_type, vec<std::size_t, 2>>;
48 using rgb_type = vec<typename color_type::value_type, 3>;
49 using alpha_type = typename color_type::value_type;
50 using ray_type = ray<cam_real_type, 3>;
51 using pixel_pos = Vec2<std::size_t>;
52 static_assert(static_vec<color_type>,
53 "Shader must return a vector with 3 or 4 components.");
54 static_assert(
55 color_type::num_components() == 3 || color_type::num_components() == 4,
56 "Shader must return a vector with 3 or 4 components.");
57 auto const& dim0 = g.template dimension<0>();
58 auto const& dim1 = g.template dimension<1>();
59 auto const& dim2 = g.template dimension<2>();
60 auto const aabb = g.bounding_box();
61 auto rendered_image =
62 rectilinear_grid<linspace<cam_real_type>, linspace<cam_real_type>>{
63 linspace<cam_real_type>{0.0, cam.plane_width() - 1,
64 cam.plane_width()},
65 linspace<cam_real_type>{0.0, cam.plane_height() - 1,
66 cam.plane_height()}};
67 auto& rendering =
68 rendered_image.template vertex_property<rgb_type>("rendered_isosurface");
69 //
70 auto rays = std::vector<std::tuple<ray_type, double, pixel_pos>>{};
71 auto mutex = std::mutex{};
72 auto const bg_color = rgb_type{1, 1, 1};
73 for_loop(
74 [&](auto const... is) {
75 rendering(is...) = bg_color;
76 auto r = cam.ray(static_cast<cam_real_type>(is)...);
77 r.normalize();
78 if (auto const i = aabb.check_intersection(r); i) {
79 auto lock = std::lock_guard {mutex};
80 rays.push_back(std::tuple{r, i->t, Vec2<std::size_t>{is...}});
81 }
82 },
83 execution_policy::parallel, cam.plane_width(), cam.plane_height());
84 for_loop(
85 [&](auto const i) {
86 auto cell_data = std::make_unique<double[]>(8);
87 auto indexing = static_multidim_size<x_fastest, 2, 2, 2>{};
88 auto accumulated_color = rgb_type::zeros();
89 if constexpr (color_type::num_components() == 3) {
90 accumulated_color = bg_color;
91 }
92 auto accumulated_alpha = alpha_type{};
93 auto const& [r, t, pixel_coord] = rays[i];
94
95 auto entry_point = r(t);
96 for (std::size_t i = 0; i < 3; ++i) {
97 if (entry_point(i) < aabb.min(i)) {
98 entry_point(i) = aabb.min(i);
99 }
100 if (entry_point(i) > aabb.max(i)) {
101 entry_point(i) = aabb.max(i);
102 }
103 }
104
105 auto cell_pos = pos_type{};
106 auto done = false;
107 auto update_cell_pos = [&](auto const& r) {
108 auto cells = g.cell_index(entry_point);
109 for (std::size_t dim = 0; dim < 3; ++dim) {
110 auto const [ci, t] = cells[dim];
111 if (gcem::abs(t) < 1e-7) {
112 cell_pos[dim] = static_cast<grid_real_type>(ci);
113 if (r.direction(dim) < 0 && ci == 0) {
114 done = true;
115 }
116 } else if (gcem::abs(t - 1) < 1e-7) {
117 cell_pos[dim] = static_cast<grid_real_type>(ci + 1);
118 if (r.direction(dim) > 0 && ci + 1 == g.size(dim) - 1) {
119 done = true;
120 }
121 } else {
122 cell_pos[dim] = static_cast<grid_real_type>(ci) + t;
123 }
124 assert(cell_pos[dim] >= 0 &&
125 cell_pos[dim] <=
126 static_cast<grid_real_type>(g.size(dim) - 1));
127 }
128 };
129 update_cell_pos(r);
130
131 while (!done) {
132 auto plane_indices_to_check = make_array<std::size_t, 3>();
133 for (std::size_t dim = 0; dim < 3; ++dim) {
134 if (cell_pos[dim] - gcem::floor(cell_pos[dim]) == 0) {
135 if (r.direction(dim) > 0) {
136 plane_indices_to_check[dim] =
137 static_cast<std::size_t>(cell_pos[dim]) + 1;
138 } else {
139 plane_indices_to_check[dim] =
140 static_cast<std::size_t>(cell_pos[dim]) - 1;
141 }
142 } else {
143 if (r.direction(dim) > 0) {
144 plane_indices_to_check[dim] = static_cast<std::size_t>(gcem::ceil(cell_pos[dim]));
145 } else {
146 plane_indices_to_check[dim] = static_cast<std::size_t>(gcem::floor(cell_pos[dim]));
147 }
148 }
149 }
150
151 assert(plane_indices_to_check[0] < dim0.size() &&
152 plane_indices_to_check[0] != std::size_t(0) - 1);
153 assert(plane_indices_to_check[1] < dim1.size() &&
154 plane_indices_to_check[1] != std::size_t(0) - 1);
155 assert(plane_indices_to_check[2] < dim2.size() &&
156 plane_indices_to_check[2] != std::size_t(0) - 1);
157 auto const t0 =
158 r.direction(0) == 0
159 ? std::numeric_limits<cam_real_type>::max()
160 : (dim0[plane_indices_to_check[0]] - r.origin(0)) /
161 r.direction(0);
162 auto const t1 =
163 r.direction(1) == 0
164 ? std::numeric_limits<cam_real_type>::max()
165 : (dim1[plane_indices_to_check[1]] - r.origin(1)) /
166 r.direction(1);
167 auto const t2 =
168 r.direction(2) == 0
169 ? std::numeric_limits<cam_real_type>::max()
170 : (dim2[plane_indices_to_check[2]] - r.origin(2)) /
171 r.direction(2);
172
173 std::array<std::size_t, 3> i0{0, 0, 0};
174 std::array<std::size_t, 3> i1{0, 0, 0};
175 for (std::size_t dim = 0; dim < 3; ++dim) {
176 if (cell_pos[dim] - gcem::floor(cell_pos[dim]) == 0) {
177 if (r.direction(dim) > 0) {
178 i0[dim] = static_cast<std::size_t>(cell_pos[dim]);
179 i1[dim] = static_cast<std::size_t>(cell_pos[dim]) + 1;
180 } else if (r.direction(dim) < 0) {
181 i0[dim] = static_cast<std::size_t>(cell_pos[dim]) - 1;
182 i1[dim] = static_cast<std::size_t>(cell_pos[dim]);
183 }
184 } else {
185 i0[dim] = static_cast<std::size_t>(gcem::floor(cell_pos[dim]));
186 i1[dim] = static_cast<std::size_t>(gcem::ceil(cell_pos[dim]));
187 }
188 }
189 if constexpr (use_indices) {
190 cell_data[indexing(0, 0, 0)] = field(i0[0], i0[1], i0[2]);
191 cell_data[indexing(1, 0, 0)] = field(i1[0], i0[1], i0[2]);
192 cell_data[indexing(0, 1, 0)] = field(i0[0], i1[1], i0[2]);
193 cell_data[indexing(1, 1, 0)] = field(i1[0], i1[1], i0[2]);
194 cell_data[indexing(0, 0, 1)] = field(i0[0], i0[1], i1[2]);
195 cell_data[indexing(1, 0, 1)] = field(i1[0], i0[1], i1[2]);
196 cell_data[indexing(0, 1, 1)] = field(i0[0], i1[1], i1[2]);
197 cell_data[indexing(1, 1, 1)] = field(i1[0], i1[1], i1[2]);
198 } else {
199 cell_data[indexing(0, 0, 0)] =
200 field(g.vertex_at(i0[0], i0[1], i0[2]));
201 cell_data[indexing(1, 0, 0)] =
202 field(g.vertex_at(i1[0], i0[1], i0[2]));
203 cell_data[indexing(0, 1, 0)] =
204 field(g.vertex_at(i0[0], i1[1], i0[2]));
205 cell_data[indexing(1, 1, 0)] =
206 field(g.vertex_at(i1[0], i1[1], i0[2]));
207 cell_data[indexing(0, 0, 1)] =
208 field(g.vertex_at(i0[0], i0[1], i1[2]));
209 cell_data[indexing(1, 0, 1)] =
210 field(g.vertex_at(i1[0], i0[1], i1[2]));
211 cell_data[indexing(0, 1, 1)] =
212 field(g.vertex_at(i0[0], i1[1], i1[2]));
213 cell_data[indexing(1, 1, 1)] =
214 field(g.vertex_at(i1[0], i1[1], i1[2]));
215 }
216
217 auto const x0 = g.vertex_at(i0[0], i0[1], i0[2]);
218 auto const x1 = g.vertex_at(i1[0], i1[1], i1[2]);
219 auto const& xa = r.origin();
220 auto const& xb = r.direction();
221 auto const cell_extent = x1 - x0;
222 auto const inv_cell_extent = pos_type{
223 1 / cell_extent(0), 1 / cell_extent(1), 1 / cell_extent(2)};
224 // create rays in different spaces
225 auto const a0 = (x1 - xa) * inv_cell_extent;
226 auto const b0 = xb * inv_cell_extent;
227 auto const a1 = (xa - x0) * inv_cell_extent;
228 auto const b1 = -xb * inv_cell_extent;
229
230 std::vector<std::tuple<grid_real_type, IsoReal, pos_type>>
231 found_surfaces;
232 for (auto const isovalue : isovalues) {
233 // check if isosurface is present in current cell
234 if (!((cell_data[indexing(0, 0, 0)] > isovalue &&
235 cell_data[indexing(0, 0, 1)] > isovalue &&
236 cell_data[indexing(0, 1, 0)] > isovalue &&
237 cell_data[indexing(0, 1, 1)] > isovalue &&
238 cell_data[indexing(1, 0, 0)] > isovalue &&
239 cell_data[indexing(1, 0, 1)] > isovalue &&
240 cell_data[indexing(1, 1, 0)] > isovalue &&
241 cell_data[indexing(1, 1, 1)] > isovalue) ||
242 (cell_data[indexing(0, 0, 0)] < isovalue &&
243 cell_data[indexing(0, 0, 1)] < isovalue &&
244 cell_data[indexing(0, 1, 0)] < isovalue &&
245 cell_data[indexing(0, 1, 1)] < isovalue &&
246 cell_data[indexing(1, 0, 0)] < isovalue &&
247 cell_data[indexing(1, 0, 1)] < isovalue &&
248 cell_data[indexing(1, 1, 0)] < isovalue &&
249 cell_data[indexing(1, 1, 1)] < isovalue))) {
250 // construct coefficients of cubic polynomial A + B*t + C*t*t +
251 // D*t*t*t
252 auto const A =
253 a0(0) * a0(1) * a0(2) * cell_data[indexing(0, 0, 0)] +
254 a0(0) * a0(1) * a1(2) * cell_data[indexing(0, 0, 1)] +
255 a0(0) * a1(1) * a0(2) * cell_data[indexing(0, 1, 0)] +
256 a0(0) * a1(1) * a1(2) * cell_data[indexing(0, 1, 1)] +
257 a1(0) * a0(1) * a0(2) * cell_data[indexing(1, 0, 0)] +
258 a1(0) * a0(1) * a1(2) * cell_data[indexing(1, 0, 1)] +
259 a1(0) * a1(1) * a0(2) * cell_data[indexing(1, 1, 0)] +
260 a1(0) * a1(1) * a1(2) * cell_data[indexing(1, 1, 1)] -
261 isovalue;
262 auto const B = (b0(0) * a0(1) * a0(2) + a0(0) * b0(1) * a0(2) +
263 a0(0) * a0(1) * b0(2)) *
264 cell_data[indexing(0, 0, 0)] +
265 (b0(0) * a0(1) * a1(2) + a0(0) * b0(1) * a1(2) +
266 a0(0) * a0(1) * b1(2)) *
267 cell_data[indexing(0, 0, 1)] +
268 (b0(0) * a1(1) * a0(2) + a0(0) * b1(1) * a0(2) +
269 a0(0) * a1(1) * b0(2)) *
270 cell_data[indexing(0, 1, 0)] +
271 (b0(0) * a1(1) * a1(2) + a0(0) * b1(1) * a1(2) +
272 a0(0) * a1(1) * b1(2)) *
273 cell_data[indexing(0, 1, 1)] +
274 (b1(0) * a0(1) * a0(2) + a1(0) * b0(1) * a0(2) +
275 a1(0) * a0(1) * b0(2)) *
276 cell_data[indexing(1, 0, 0)] +
277 (b1(0) * a0(1) * a1(2) + a1(0) * b0(1) * a1(2) +
278 a1(0) * a0(1) * b1(2)) *
279 cell_data[indexing(1, 0, 1)] +
280 (b1(0) * a1(1) * a0(2) + a1(0) * b1(1) * a0(2) +
281 a1(0) * a1(1) * b0(2)) *
282 cell_data[indexing(1, 1, 0)] +
283 (b1(0) * a1(1) * a1(2) + a1(0) * b1(1) * a1(2) +
284 a1(0) * a1(1) * b1(2)) *
285 cell_data[indexing(1, 1, 1)];
286 auto const C = (a0(0) * b0(1) * b0(2) + b0(0) * a0(1) * b0(2) +
287 b0(0) * b0(1) * a0(2)) *
288 cell_data[indexing(0, 0, 0)] +
289 (a0(0) * b0(1) * b1(2) + b0(0) * a0(1) * b1(2) +
290 b0(0) * b0(1) * a1(2)) *
291 cell_data[indexing(0, 0, 1)] +
292 (a0(0) * b1(1) * b0(2) + b0(0) * a1(1) * b0(2) +
293 b0(0) * b1(1) * a0(2)) *
294 cell_data[indexing(0, 1, 0)] +
295 (a0(0) * b1(1) * b1(2) + b0(0) * a1(1) * b1(2) +
296 b0(0) * b1(1) * a1(2)) *
297 cell_data[indexing(0, 1, 1)] +
298 (a1(0) * b0(1) * b0(2) + b1(0) * a0(1) * b0(2) +
299 b1(0) * b0(1) * a0(2)) *
300 cell_data[indexing(1, 0, 0)] +
301 (a1(0) * b0(1) * b1(2) + b1(0) * a0(1) * b1(2) +
302 b1(0) * b0(1) * a1(2)) *
303 cell_data[indexing(1, 0, 1)] +
304 (a1(0) * b1(1) * b0(2) + b1(0) * a1(1) * b0(2) +
305 b1(0) * b1(1) * a0(2)) *
306 cell_data[indexing(1, 1, 0)] +
307 (a1(0) * b1(1) * b1(2) + b1(0) * a1(1) * b1(2) +
308 b1(0) * b1(1) * a1(2)) *
309 cell_data[indexing(1, 1, 1)];
310 auto const D =
311 b0(0) * b0(1) * b0(2) * cell_data[indexing(0, 0, 0)] +
312 b0(0) * b0(1) * b1(2) * cell_data[indexing(0, 0, 1)] +
313 b0(0) * b1(1) * b0(2) * cell_data[indexing(0, 1, 0)] +
314 b0(0) * b1(1) * b1(2) * cell_data[indexing(0, 1, 1)] +
315 b1(0) * b0(1) * b0(2) * cell_data[indexing(1, 0, 0)] +
316 b1(0) * b0(1) * b1(2) * cell_data[indexing(1, 0, 1)] +
317 b1(0) * b1(1) * b0(2) * cell_data[indexing(1, 1, 0)] +
318 b1(0) * b1(1) * b1(2) * cell_data[indexing(1, 1, 1)];
319
320 auto const s = solve(polynomial{A, B, C, D});
321
322 if (!s.empty()) {
323 for (auto const t : s) {
324 constexpr auto eps = 1e-10;
325 if (auto x = a0 + t * b0; //
326 - eps <= x(0) && x(0) <= 1 + eps && //
327 -eps <= x(1) && x(1) <= 1 + eps && //
328 -eps <= x(2) && x(2) <= 1 + eps) {
329 found_surfaces.emplace_back(t, isovalue, x);
330 }
331 }
332 }
333 }
334 }
335 std::sort(begin(found_surfaces), end(found_surfaces),
336 [](auto const& s0, auto const& s1) {
337 auto const& [t0, iso0, uvw0] = s0;
338 auto const& [t1, iso1, uvw1] = s1;
339 return t0 < t1;
340 });
341 for (auto const& [t, isovalue, uvw1] : found_surfaces) {
342 auto const uvw0 = pos_type{1 - uvw1(0), 1 - uvw1(1), 1 - uvw1(2)};
343 auto const x_iso = uvw0 * cell_extent + x0;
344 assert(uvw0(0) >= 0 && uvw0(0) <= 1);
345 assert(uvw0(1) >= 0 && uvw0(1) <= 1);
346 assert(uvw0(2) >= 0 && uvw0(2) <= 1);
347 assert(uvw1(0) >= 0 && uvw1(0) <= 1);
348 assert(uvw1(1) >= 0 && uvw1(1) <= 1);
349 assert(uvw1(2) >= 0 && uvw1(2) <= 1);
350 auto const k =
351 cell_data[indexing(1, 1, 1)] - cell_data[indexing(0, 1, 1)] -
352 cell_data[indexing(1, 0, 1)] + cell_data[indexing(0, 0, 1)] -
353 cell_data[indexing(1, 1, 0)] + cell_data[indexing(0, 1, 0)] +
354 cell_data[indexing(1, 0, 0)] - cell_data[indexing(0, 0, 0)];
355 auto const gradient = vec{
356 (k * uvw0(1) + cell_data[indexing(1, 0, 1)] -
357 cell_data[indexing(0, 0, 1)] - cell_data[indexing(1, 0, 0)] +
358 cell_data[indexing(0, 0, 0)]) *
359 uvw0(2) +
360 (cell_data[indexing(1, 1, 0)] -
361 cell_data[indexing(0, 1, 0)] -
362 cell_data[indexing(1, 0, 0)] +
363 cell_data[indexing(0, 0, 0)]) *
364 uvw0(1) +
365 cell_data[indexing(1, 0, 0)] - cell_data[indexing(0, 0, 0)],
366 (k * uvw0(0) + cell_data[indexing(0, 1, 1)] -
367 cell_data[indexing(0, 0, 1)] - cell_data[indexing(0, 1, 0)] +
368 cell_data[indexing(0, 0, 0)]) *
369 uvw0(2) +
370 (cell_data[indexing(1, 1, 0)] -
371 cell_data[indexing(0, 1, 0)] -
372 cell_data[indexing(1, 0, 0)] +
373 cell_data[indexing(0, 0, 0)]) *
374 uvw0(0) +
375 cell_data[indexing(0, 1, 0)] - cell_data[indexing(0, 0, 0)],
376 (k * uvw0(0) + cell_data[indexing(0, 1, 1)] -
377 cell_data[indexing(0, 0, 1)] - cell_data[indexing(0, 1, 0)] +
378 cell_data[indexing(0, 0, 0)]) *
379 uvw0(1) +
380 (cell_data[indexing(1, 0, 1)] -
381 cell_data[indexing(0, 0, 1)] -
382 cell_data[indexing(1, 0, 0)] +
383 cell_data[indexing(0, 0, 0)]) *
384 uvw0(0) +
385 cell_data[indexing(0, 0, 1)] -
386 cell_data[indexing(0, 0, 0)]};
387 if constexpr (color_type::num_components() == 3) {
388 accumulated_color =
389 shader(x_iso, isovalue, gradient, r.direction(), pixel_coord);
390 done = true;
391 } else if constexpr (color_type::num_components() == 4) {
392 auto const rgba =
393 shader(x_iso, isovalue, gradient, r.direction(), pixel_coord);
394 auto const rgb = vec{rgba(0), rgba(1), rgba(2)};
395 auto const alpha = rgba(3);
396 accumulated_color += (1 - accumulated_alpha) * alpha * rgb;
397 accumulated_alpha += (1 - accumulated_alpha) * alpha;
398 if (accumulated_alpha >= 0.95) {
399 done = true;
400 }
401 }
402 }
403
404 if (!done) {
405 entry_point = r(tatooine::min(t0, t1, t2));
406 update_cell_pos(r);
407 }
408 }
409 if constexpr (color_type::num_components() == 3) {
410 rendering(pixel_coord.x(), pixel_coord.y()) = accumulated_color;
411 } else if constexpr (color_type::num_components() == 4) {
412 rendering(pixel_coord.x(), pixel_coord.y()) =
413 accumulated_color * accumulated_alpha +
414 bg_color * (1 - accumulated_alpha);
415 }
416 },
417 execution_policy::parallel, rays.size());
418 return rendered_image;
419}
420//==============================================================================
421template <typename IsoReal, typename Dim0, typename Dim1, typename Dim2,
422 typename Field, typename FieldReal, typename Shader>
423auto direct_isosurface(camera auto const& cam,
424 rectilinear_grid<Dim0, Dim1, Dim2> const& g,
425 scalarfield<Field, FieldReal, 3> const& field,
426 IsoReal const isovalue, Shader&& shader) {
427 return direct_isosurface(cam, g, field, std::vector{isovalue},
428 std::forward<Shader>(shader));
429}
430//==============================================================================
441template <camera Camera, typename GridVertexProperty, typename Shader,
442 arithmetic Iso>
443auto direct_isosurface(
444 Camera const& cam,
445 tatooine::detail::rectilinear_grid::vertex_property_sampler<
446 GridVertexProperty, interpolation::linear, interpolation::linear,
447 interpolation::linear> const& linear_field,
448 Iso const isovalue, Shader&& shader)
449 //
450 requires invocable<
451 Shader, vec<typename GridVertexProperty::grid_type::real_type, 3>,
452 Iso, vec<typename GridVertexProperty::grid_type::real_type, 3>,
453 vec<typename Camera::real_type, 3>, vec<std::size_t, 2>>
454//
455{
456 return direct_isosurface(
457 cam, linear_field.grid(),
458 [&](std::size_t const ix, std::size_t const iy, std::size_t const iz)
459 -> auto const& { return linear_field.data_at(ix, iy, iz); },
460 std::vector{isovalue}, std::forward<Shader>(shader));
461}
462//------------------------------------------------------------------------------
463template <typename DistOnRay, typename AABBReal, typename DataEvaluator,
464 typename Isovalue, typename DomainCheck, typename Shader>
465auto direct_isosurface(camera auto const& cam,
466 axis_aligned_bounding_box<AABBReal, 3> const& aabb,
467 DataEvaluator&& data_evaluator,
468 DomainCheck&& domain_check, Isovalue isovalue,
469 DistOnRay const distance_on_ray, Shader&& shader) {
470 using cam_real_type = typename std::decay_t<decltype(cam)>::real_type;
471 using pos_type = vec<cam_real_type, 3>;
472 using viewdir_type = vec<cam_real_type, 3>;
473 static_assert(
474 std::is_invocable_v<Shader, pos_type, viewdir_type, vec<std::size_t, 2>>,
475 "Shader must be invocable with position and view direction.");
476 using value_type = std::invoke_result_t<DataEvaluator, pos_type>;
477 using color_type = std::invoke_result_t<Shader, pos_type, viewdir_type>;
478 using rgb_type = vec<typename color_type::value_type, 3>;
479 using alpha_type = typename color_type::value_type;
480 using ray_type = ray<cam_real_type, 3>;
481 static_assert(is_floating_point<value_type>,
482 "DataEvaluator must return scalar type.");
483 static_assert(static_vec<color_type>,
484 "ColorScale must return scalar type or tatooine::vec.");
485 static_assert(
486 color_type::num_components() == 3 || color_type::num_components() == 4,
487 "ColorScale must return scalar type or tatooine::vec.");
488 auto rendered_image =
489 rectilinear_grid<linspace<cam_real_type>, linspace<cam_real_type>>{
490 linspace<cam_real_type>{0.0, cam.plane_width() - 1,
491 cam.plane_width()},
492 linspace<cam_real_type>{0.0, cam.plane_height() - 1,
493 cam.plane_height()}};
494 auto& rendering =
495 rendered_image.template vertex_property<rgb_type>("rendered_isosurface");
496
497 auto rays = std::vector<
498 std::tuple<ray_type, AABBReal, std::size_t, std::size_t>>{};
499 auto mutex = std::mutex{};
500 auto const bg_color = rgb_type{1, 1, 1};
501 for_loop(
502 [&](auto const... is) {
503 rendering(is...) = bg_color;
504 auto r = cam.ray(is...);
505 r.normalize();
506 if (auto const i = aabb.check_intersection(r); i) {
507 std::lock_guard lock{mutex};
508 rays.push_back(std::tuple{r, i->t, is...});
509 }
510 },
511 execution_policy::parallel, cam.plane_width(), cam.plane_height());
512 for_loop(
513 [&](auto const i) {
514 auto const [r, t, x, y] = rays[i];
515 auto accumulated_alpha = alpha_type{};
516 auto accumulated_color = rgb_type{};
517
518 auto t0 = t;
519 auto x0 = r(t0);
520 for (std::size_t i = 0; i < 3; ++i) {
521 if (x0(i) < aabb.min(i)) {
522 x0(i) = aabb.min(i);
523 }
524 if (x0(i) > aabb.max(i)) {
525 x0(i) = aabb.max(i);
526 }
527 }
528 auto sample0 = data_evaluator(x0);
529
530 auto t1 = t0 + distance_on_ray;
531 auto x1 = r(t1);
532 auto sample1 = sample0;
533
534 auto done = false;
535 while (!done && aabb.is_inside(x1)) {
536 if (domain_check(x0) && domain_check(x1)) {
537 sample1 = data_evaluator(x1);
538 if ((sample0 <= isovalue && sample1 > isovalue) ||
539 (sample0 >= isovalue && sample1 < isovalue)) {
540 auto cur_x0 = x0;
541 auto cur_x1 = x1;
542 auto cur_sample0 = sample0;
543 auto cur_sample1 = sample1;
544 for (std::size_t i = 0; i < 100; ++i) {
545 auto x_center = (cur_x0 + cur_x1) / 2;
546 auto sample_center = data_evaluator(x_center);
547 if ((cur_sample0 <= isovalue && sample_center > isovalue) ||
548 (cur_sample0 >= isovalue && sample_center < isovalue)) {
549 cur_x1 = x_center;
550 cur_sample1 = sample_center;
551 } else {
552 cur_x0 = x_center;
553 cur_sample0 = sample_center;
554 }
555 }
556 auto const t_iso =
557 (isovalue - cur_sample0) / (cur_sample1 - cur_sample0);
558 auto const iso_pos = r(t0 + t_iso * distance_on_ray);
559
560 if constexpr (color_type::num_components() == 3) {
561 accumulated_color = shader(iso_pos, r.direction());
562 done = true;
563 } else if constexpr (color_type::num_components() == 4) {
564 auto const rgba = shader(iso_pos, r.direction());
565 auto const rgb = vec{rgba(0), rgba(1), rgba(2)};
566 auto const alpha = rgba(3);
567 accumulated_color += (1 - accumulated_alpha) * alpha * rgb;
568 accumulated_alpha += (1 - accumulated_alpha) * alpha;
569 if (accumulated_alpha >= 0.95) {
570 done = true;
571 }
572 }
573 }
574 }
575 t0 = t1;
576 x0 = std::move(x1);
577 sample0 = std::move(sample1);
578 t1 += distance_on_ray;
579 x1 = r(t1);
580 }
581 rendering(x, y) = accumulated_color * accumulated_alpha +
582 bg_color * (1 - accumulated_alpha);
583 },
584 execution_policy::parallel, rays.size());
585 return rendered_image;
586}
587//==============================================================================
588} // namespace tatooine::rendering
589//==============================================================================
590#endif
real_type
tatooine::rectilinear_grid
Definition: rectilinear_grid.h:38
tatooine::rectilinear_grid::real_type
common_type< typename Dimensions::value_type... > real_type
Definition: rectilinear_grid.h:46
tatooine::arithmetic
Definition: concepts.h:33
tatooine::invocable
Definition: concepts.h:121
tatooine::rendering::camera
Definition: camera.h:312
tatooine::static_vec
Definition: tensor_concepts.h:23
rectilinear_grid.h
tatooine::execution_policy::parallel
static constexpr parallel_t parallel
Definition: tags.h:60
tatooine::rendering
Definition: camera.h:12
tatooine::rendering::direct_isosurface
auto direct_isosurface(Camera const &cam, rectilinear_grid< Dim0, Dim1, Dim2 > const &g, Field &&field, std::vector< IsoReal > const &isovalues, Shader &&shader)
Definition: direct_isosurface.h:24
tatooine::begin
auto begin(Range &&range)
Definition: iterator_facade.h:318
tatooine::value_type
typename value_type_impl< T >::type value_type
Definition: type_traits.h:280
tatooine::end
auto end(Range &&range)
Definition: iterator_facade.h:322
tatooine::invoke_result
std::invoke_result_t< F, Args... > invoke_result
Definition: concepts.h:127
tatooine::solve
auto solve(polynomial< Real, 1 > const &p) -> std::vector< Real >
solve a + b*x
Definition: polynomial.h:187
tatooine::min
constexpr auto min(A &&a, B &&b)
Definition: math.h:15
tatooine::for_loop
constexpr auto for_loop(Iteration &&iteration, execution_policy::sequential_t, Ranges(&&... ranges)[2]) -> void
Use this function for creating a sequential nested loop.
Definition: for_loop.h:336
tatooine::axis_aligned_bounding_box
Definition: axis_aligned_bounding_box.h:103
tatooine::axis_aligned_bounding_box::is_inside
auto constexpr is_inside(pos_type const &p) const
Definition: axis_aligned_bounding_box.h:191
tatooine::axis_aligned_bounding_box::max
auto constexpr max() const -> auto const &
Definition: axis_aligned_bounding_box.h:156
tatooine::axis_aligned_bounding_box::min
auto constexpr min() const -> auto const &
Definition: axis_aligned_bounding_box.h:151
tatooine::detail::aabb_ray_intersectable_parent::check_intersection
auto check_intersection(ray_type const &r, Real const =0) const -> optional_intersection_type override
Definition: axis_aligned_bounding_box.h:33
tatooine::detail::rectilinear_grid::vertex_property_sampler
Definition: vertex_property_sampler.h:271
tatooine::field
Definition: field.h:134
tatooine::interpolation::linear
Definition: interpolation.h:16
tatooine::linspace
Definition: linspace.h:26
tatooine::polynomial
Definition: polynomial.h:14
tatooine::ray
Definition: ray.h:10
tatooine::static_multidim_size
Definition: static_multidim_size.h:81
tatooine::vec
Definition: vec.h:12
Generated by doxygen 1.9.5