3#ifndef TATOOINE_RENDERING_DIRECT_ISOSURFACE_H
4#define TATOOINE_RENDERING_DIRECT_ISOSURFACE_H
22template <camera Camera, arithmetic IsoReal,
typename Dim0,
typename Dim1,
23 typename Dim2,
typename Field,
typename Shader>
26 Field&&
field, std::vector<IsoReal>
const& isovalues,
37 constexpr auto use_indices =
38 std::is_invocable_v<Field, std::size_t, std::size_t, std::size_t>;
44 using cam_real_type =
typename Camera::real_type;
46 using color_type =
invoke_result<Shader, pos_type, IsoReal, pos_type,
49 using alpha_type =
typename color_type::value_type;
53 "Shader must return a vector with 3 or 4 components.");
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();
68 rendered_image.template vertex_property<rgb_type>(
"rendered_isosurface");
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};
74 [&](
auto const... is) {
75 rendering(is...) = bg_color;
76 auto r = cam.ray(
static_cast<cam_real_type
>(is)...);
79 auto lock = std::lock_guard {mutex};
86 auto cell_data = std::make_unique<double[]>(8);
88 auto accumulated_color = rgb_type::zeros();
89 if constexpr (color_type::num_components() == 3) {
90 accumulated_color = bg_color;
92 auto accumulated_alpha = alpha_type{};
93 auto const& [r, t, pixel_coord] = rays[i];
95 auto entry_point = r(t);
96 for (std::size_t i = 0; i < 3; ++i) {
97 if (entry_point(i) <
aabb.
min(i)) {
100 if (entry_point(i) >
aabb.
max(i)) {
105 auto cell_pos = pos_type{};
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) {
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) {
122 cell_pos[dim] =
static_cast<grid_real_type
>(ci) + t;
124 assert(cell_pos[dim] >= 0 &&
126 static_cast<grid_real_type
>(g.size(dim) - 1));
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;
139 plane_indices_to_check[dim] =
140 static_cast<std::size_t
>(cell_pos[dim]) - 1;
143 if (r.direction(dim) > 0) {
144 plane_indices_to_check[dim] =
static_cast<std::size_t
>(gcem::ceil(cell_pos[dim]));
146 plane_indices_to_check[dim] =
static_cast<std::size_t
>(gcem::floor(cell_pos[dim]));
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);
159 ? std::numeric_limits<cam_real_type>::max()
160 : (dim0[plane_indices_to_check[0]] - r.origin(0)) /
164 ? std::numeric_limits<cam_real_type>::max()
165 : (dim1[plane_indices_to_check[1]] - r.origin(1)) /
169 ? std::numeric_limits<cam_real_type>::max()
170 : (dim2[plane_indices_to_check[2]] - r.origin(2)) /
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]);
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]));
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]);
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]));
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)};
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;
230 std::vector<std::tuple<grid_real_type, IsoReal, pos_type>>
232 for (
auto const isovalue : isovalues) {
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))) {
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)] -
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)];
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)];
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);
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;
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);
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)]) *
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)]) *
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)]) *
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)]) *
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)]) *
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)]) *
385 cell_data[indexing(0, 0, 1)] -
386 cell_data[indexing(0, 0, 0)]};
387 if constexpr (color_type::num_components() == 3) {
389 shader(x_iso, isovalue, gradient, r.direction(), pixel_coord);
391 }
else if constexpr (color_type::num_components() == 4) {
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) {
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);
418 return rendered_image;
421template <
typename IsoReal,
typename Dim0,
typename Dim1,
typename Dim2,
422 typename Field,
typename FieldReal,
typename Shader>
426 IsoReal
const isovalue, Shader&& shader) {
428 std::forward<Shader>(shader));
441template <camera Camera,
typename GridVertexProperty,
typename Shader,
448 Iso
const isovalue, Shader&& shader)
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));
463template <
typename DistOnRay,
typename AABBReal,
typename DataEvaluator,
464 typename Isovalue,
typename DomainCheck,
typename Shader>
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;
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>;
479 using alpha_type =
typename color_type::value_type;
481 static_assert(is_floating_point<value_type>,
482 "DataEvaluator must return scalar type.");
484 "ColorScale must return scalar type or tatooine::vec.");
486 color_type::num_components() == 3 || color_type::num_components() == 4,
487 "ColorScale must return scalar type or tatooine::vec.");
488 auto rendered_image =
493 cam.plane_height()}};
495 rendered_image.template vertex_property<rgb_type>(
"rendered_isosurface");
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};
502 [&](
auto const... is) {
503 rendering(is...) = bg_color;
504 auto r = cam.ray(is...);
507 std::lock_guard lock{mutex};
508 rays.push_back(std::tuple{r, i->t, is...});
514 auto const [r, t, x, y] = rays[i];
515 auto accumulated_alpha = alpha_type{};
516 auto accumulated_color = rgb_type{};
520 for (std::size_t i = 0; i < 3; ++i) {
528 auto sample0 = data_evaluator(x0);
530 auto t1 = t0 + distance_on_ray;
532 auto sample1 = sample0;
536 if (domain_check(x0) && domain_check(x1)) {
537 sample1 = data_evaluator(x1);
538 if ((sample0 <= isovalue && sample1 > isovalue) ||
539 (sample0 >= isovalue && sample1 < isovalue)) {
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)) {
550 cur_sample1 = sample_center;
553 cur_sample0 = sample_center;
557 (isovalue - cur_sample0) / (cur_sample1 - cur_sample0);
558 auto const iso_pos = r(t0 + t_iso * distance_on_ray);
560 if constexpr (color_type::num_components() == 3) {
561 accumulated_color = shader(iso_pos, r.direction());
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) {
577 sample0 = std::move(sample1);
578 t1 += distance_on_ray;
581 rendering(x, y) = accumulated_color * accumulated_alpha +
582 bg_color * (1 - accumulated_alpha);
585 return rendered_image;
Definition: rectilinear_grid.h:38
common_type< typename Dimensions::value_type... > real_type
Definition: rectilinear_grid.h:46
Definition: concepts.h:33
Definition: concepts.h:121
Definition: tensor_concepts.h:23
static constexpr parallel_t parallel
Definition: tags.h:60
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
auto begin(Range &&range)
Definition: iterator_facade.h:318
typename value_type_impl< T >::type value_type
Definition: type_traits.h:280
auto end(Range &&range)
Definition: iterator_facade.h:322
std::invoke_result_t< F, Args... > invoke_result
Definition: concepts.h:127
auto solve(polynomial< Real, 1 > const &p) -> std::vector< Real >
solve a + b*x
Definition: polynomial.h:187
constexpr auto min(A &&a, B &&b)
Definition: math.h:15
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
Definition: axis_aligned_bounding_box.h:103
auto constexpr is_inside(pos_type const &p) const
Definition: axis_aligned_bounding_box.h:191
auto constexpr max() const -> auto const &
Definition: axis_aligned_bounding_box.h:156
auto constexpr min() const -> auto const &
Definition: axis_aligned_bounding_box.h:151
auto check_intersection(ray_type const &r, Real const =0) const -> optional_intersection_type override
Definition: axis_aligned_bounding_box.h:33
Definition: vertex_property_sampler.h:271
Definition: interpolation.h:16
Definition: linspace.h:26
Definition: polynomial.h:14
Definition: static_multidim_size.h:81