atrstm_setup_octrees.c (25821B)
1 /* Copyright (C) 2022, 2023 |Méso|Star> (contact@meso-star.com) 2 * Copyright (C) 2020, 2021 Centre National de la Recherche Scientifique 3 * 4 * This program is free software: you can redistribute it and/or modify 5 * it under the terms of the GNU General Public License as published by 6 * the Free Software Foundation, either version 3 of the License, or 7 * (at your option) any later version. 8 * 9 * This program is distributed in the hope that it will be useful, 10 * but WITHOUT ANY WARRANTY; without even the implied warranty of 11 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 12 * GNU General Public License for more details. 13 * 14 * You should have received a copy of the GNU General Public License 15 * along with this program. If not, see <http://www.gnu.org/licenses/>. */ 16 17 #define _POSIX_C_SOURCE 200112L /* nextafterf */ 18 19 #include "atrstm_c.h" 20 #include "atrstm_cache.h" 21 #include "atrstm_log.h" 22 #include "atrstm_radcoefs.h" 23 #include "atrstm_partition.h" 24 #include "atrstm_setup_octrees.h" 25 #include "atrstm_svx.h" 26 27 #ifdef ATRSTM_USE_SIMD 28 #include "atrstm_radcoefs_simd4.h" 29 #endif 30 31 #include <astoria/atrri.h> 32 33 #include <star/suvm.h> 34 #include <star/svx.h> 35 36 #include <rsys/clock_time.h> 37 #include <rsys/cstr.h> 38 #include <rsys/dynamic_array_size_t.h> 39 #include <rsys/morton.h> 40 41 #include <math.h> 42 #include <omp.h> 43 44 struct build_octree_context { 45 struct atrstm* atrstm; 46 struct pool* pool; /* Pool of voxel partitions */ 47 struct part* part; /* Current partition */ 48 49 /* Optical thickness threshold criteria for the merge process */ 50 double tau_threshold; 51 }; 52 #define BUILD_OCTREE_CONTEXT_NULL__ { NULL, NULL, NULL, 0 } 53 static const struct build_octree_context BUILD_OCTREE_CONTEXT_NULL = 54 BUILD_OCTREE_CONTEXT_NULL__; 55 56 /* Define the darray of list of primitive ids */ 57 #define DARRAY_NAME prims_list 58 #define DARRAY_DATA struct darray_size_t 59 #define DARRAY_FUNCTOR_INIT darray_size_t_init 60 #define DARRAY_FUNCTOR_RELEASE darray_size_t_release 61 #define DARRAY_FUNCTOR_COPY darray_size_t_copy 62 #define DARRAY_FUNCTOR_COPY_AND_RELEASE darray_size_t_copy_and_release 63 #include <rsys/dynamic_array.h> 64 65 /******************************************************************************* 66 * Helper functions 67 ******************************************************************************/ 68 static INLINE int 69 check_octrees_parameters(const struct atrstm* atrstm) 70 { 71 return atrstm 72 && atrstm->grid_max_definition[0] 73 && atrstm->grid_max_definition[1] 74 && atrstm->grid_max_definition[2] 75 && atrstm->optical_thickness >= 0; 76 } 77 78 static INLINE void 79 register_primitive 80 (const struct suvm_primitive* prim, 81 const double low[3], 82 const double upp[3], 83 void* context) 84 { 85 struct darray_size_t* prims = context; 86 ASSERT(prim && low && upp && context); 87 ASSERT(low[0] < upp[0]); 88 ASSERT(low[1] < upp[1]); 89 ASSERT(low[2] < upp[2]); 90 (void)low, (void)upp; 91 CHK(darray_size_t_push_back(prims, &prim->iprim) == RES_OK); 92 } 93 94 static void 95 voxel_commit_radcoefs_range 96 (float* vox, 97 const enum atrstm_component cpnt, 98 const struct radcoefs* radcoefs_min, 99 const struct radcoefs* radcoefs_max) 100 { 101 double ka[2], ks[2], kext[2]; 102 float vox_ka[2], vox_ks[2], vox_kext[2]; 103 ASSERT(vox); 104 ASSERT(radcoefs_min); 105 ASSERT(radcoefs_max); 106 107 ka[0] = radcoefs_min->ka; 108 ka[1] = radcoefs_max->ka; 109 ks[0] = radcoefs_min->ks; 110 ks[1] = radcoefs_max->ks; 111 kext[0] = radcoefs_min->kext; 112 kext[1] = radcoefs_max->kext; 113 114 /* Ensure that the single precision bounds include their double precision 115 * version. */ 116 if(ka[0] != (float)ka[0]) ka[0] = nextafterf((float)ka[0],-FLT_MAX); 117 if(ka[1] != (float)ka[1]) ka[1] = nextafterf((float)ka[1], FLT_MAX); 118 if(ks[0] != (float)ks[0]) ks[0] = nextafterf((float)ks[0],-FLT_MAX); 119 if(ks[1] != (float)ks[1]) ks[1] = nextafterf((float)ks[1], FLT_MAX); 120 if(kext[0] != (float)kext[0]) kext[0] = nextafterf((float)kext[0],-FLT_MAX); 121 if(kext[1] != (float)kext[1]) kext[1] = nextafterf((float)kext[1], FLT_MAX); 122 123 /* Fetch the range of the voxel optical properties */ 124 vox_ka[0] = vox_get(vox, cpnt, ATRSTM_RADCOEF_Ka, ATRSTM_SVX_OP_MIN); 125 vox_ka[1] = vox_get(vox, cpnt, ATRSTM_RADCOEF_Ka, ATRSTM_SVX_OP_MAX); 126 vox_ks[0] = vox_get(vox, cpnt, ATRSTM_RADCOEF_Ks, ATRSTM_SVX_OP_MIN); 127 vox_ks[1] = vox_get(vox, cpnt, ATRSTM_RADCOEF_Ks, ATRSTM_SVX_OP_MAX); 128 vox_kext[0] = vox_get(vox, cpnt, ATRSTM_RADCOEF_Kext, ATRSTM_SVX_OP_MIN); 129 vox_kext[1] = vox_get(vox, cpnt, ATRSTM_RADCOEF_Kext, ATRSTM_SVX_OP_MAX); 130 131 /* Update the range of the voxel optical properties */ 132 vox_ka[0] = MMIN(vox_ka[0], (float)ka[0]); 133 vox_ka[1] = MMAX(vox_ka[1], (float)ka[1]); 134 vox_ks[0] = MMIN(vox_ks[0], (float)ks[0]); 135 vox_ks[1] = MMAX(vox_ks[1], (float)ks[1]); 136 vox_kext[0] = MMIN(vox_kext[0], (float)kext[0]); 137 vox_kext[1] = MMAX(vox_kext[1], (float)kext[1]); 138 139 /* Commit the upd to the voxel */ 140 vox_set(vox, cpnt, ATRSTM_RADCOEF_Ka, ATRSTM_SVX_OP_MIN, vox_ka[0]); 141 vox_set(vox, cpnt, ATRSTM_RADCOEF_Ka, ATRSTM_SVX_OP_MAX, vox_ka[1]); 142 vox_set(vox, cpnt, ATRSTM_RADCOEF_Ks, ATRSTM_SVX_OP_MIN, vox_ks[0]); 143 vox_set(vox, cpnt, ATRSTM_RADCOEF_Ks, ATRSTM_SVX_OP_MAX, vox_ks[1]); 144 vox_set(vox, cpnt, ATRSTM_RADCOEF_Kext, ATRSTM_SVX_OP_MIN, vox_kext[0]); 145 vox_set(vox, cpnt, ATRSTM_RADCOEF_Kext, ATRSTM_SVX_OP_MAX, vox_kext[1]); 146 } 147 148 static res_T 149 voxelize_partition 150 (struct atrstm* atrstm, 151 const struct darray_size_t* prims, /* Primitives intersecting the partition */ 152 const double part_low[3], /* Lower bound of the partition */ 153 const double part_upp[3], /* Upper bound of the partition */ 154 const double vxsz[3], /* Size of a partition voxel */ 155 const struct atrri_refractive_index* refract_id, /* Refractive index */ 156 struct part* part) 157 { 158 size_t iprim; 159 res_T res = RES_OK; 160 ASSERT(atrstm && prims && part_low && part_upp && vxsz && part && refract_id); 161 ASSERT(part_low[0] < part_upp[0]); 162 ASSERT(part_low[1] < part_upp[1]); 163 ASSERT(part_low[2] < part_upp[2]); 164 ASSERT(vxsz[0] > 0 && vxsz[1] > 0 && vxsz[2] > 0); 165 166 /* Clean the partition voxels */ 167 part_clear_voxels(part); 168 169 FOR_EACH(iprim, 0, darray_size_t_size_get(prims)) { 170 struct suvm_primitive prim = SUVM_PRIMITIVE_NULL; 171 struct suvm_polyhedron poly; 172 double poly_low[3]; 173 double poly_upp[3]; 174 uint32_t ivxl_low[3]; 175 uint32_t ivxl_upp[3]; 176 uint32_t ivxl[3]; 177 size_t prim_id; 178 179 prim_id = darray_size_t_cdata_get(prims)[iprim]; 180 181 /* Retrieve the primitive data and setup its polyhedron */ 182 SUVM(volume_get_primitive(atrstm->volume, prim_id, &prim)); 183 SUVM(primitive_setup_polyhedron(&prim, &poly)); 184 ASSERT(poly.lower[0] <= part_upp[0] && poly.upper[0] >= part_low[0]); 185 ASSERT(poly.lower[1] <= part_upp[1] && poly.upper[1] >= part_low[1]); 186 ASSERT(poly.lower[2] <= part_upp[2] && poly.upper[2] >= part_low[2]); 187 188 /* Clamp the poly AABB to the partition boundaries */ 189 poly_low[0] = MMAX(poly.lower[0], part_low[0]); 190 poly_low[1] = MMAX(poly.lower[1], part_low[1]); 191 poly_low[2] = MMAX(poly.lower[2], part_low[2]); 192 poly_upp[0] = MMIN(poly.upper[0], part_upp[0]); 193 poly_upp[1] = MMIN(poly.upper[1], part_upp[1]); 194 poly_upp[2] = MMIN(poly.upper[2], part_upp[2]); 195 196 /* Transform the polyhedron AABB in partition voxel space */ 197 ivxl_low[0] = (uint32_t)((poly_low[0] - part_low[0]) / vxsz[0]); 198 ivxl_low[1] = (uint32_t)((poly_low[1] - part_low[1]) / vxsz[1]); 199 ivxl_low[2] = (uint32_t)((poly_low[2] - part_low[2]) / vxsz[2]); 200 ivxl_upp[0] = (uint32_t)ceil((poly_upp[0] - part_low[0]) / vxsz[0]); 201 ivxl_upp[1] = (uint32_t)ceil((poly_upp[1] - part_low[1]) / vxsz[1]); 202 ivxl_upp[2] = (uint32_t)ceil((poly_upp[2] - part_low[2]) / vxsz[2]); 203 ASSERT(ivxl_upp[0] <= PARTITION_DEFINITION); 204 ASSERT(ivxl_upp[1] <= PARTITION_DEFINITION); 205 ASSERT(ivxl_upp[2] <= PARTITION_DEFINITION); 206 207 /* Iterate over the partition voxels intersected by the primitive AABB */ 208 FOR_EACH(ivxl[2], ivxl_low[2], ivxl_upp[2]) { 209 float vxl_low[3]; /* Voxel lower bound */ 210 float vxl_upp[3]; /* Voxel upper bound */ 211 uint64_t mcode[3]; /* Morton code cache */ 212 213 vxl_low[2] = (float)((double)ivxl[2] * vxsz[2] + part_low[2]); 214 vxl_upp[2] = vxl_low[2] + (float)vxsz[2]; 215 mcode[2] = morton3D_encode_u21(ivxl[2]) << 0; 216 217 FOR_EACH(ivxl[1], ivxl_low[1], ivxl_upp[1]) { 218 vxl_low[1] = (float)((double)ivxl[1] * vxsz[1] + part_low[1]); 219 vxl_upp[1] = vxl_low[1] + (float)vxsz[1]; 220 mcode[1] = (morton3D_encode_u21(ivxl[1]) << 1) | mcode[2]; 221 222 FOR_EACH(ivxl[0], ivxl_low[0], ivxl_upp[0]) { 223 enum suvm_intersection_type intersect; 224 vxl_low[0] = (float)((double)ivxl[0] * vxsz[0] + part_low[0]); 225 vxl_upp[0] = vxl_low[0] + (float)vxsz[0]; 226 227 /* NOTE mcode[0] store the morton index of the voxel */ 228 mcode[0] = (morton3D_encode_u21(ivxl[0]) << 2) | mcode[1]; 229 230 intersect = suvm_polyhedron_intersect_aabb(&poly, vxl_low, vxl_upp); 231 if(intersect != SUVM_INTERSECT_NONE) { 232 struct radcoefs radcoefs_min; 233 struct radcoefs radcoefs_max; 234 float* vox = part_get_voxel(part, mcode[0]); 235 236 /* Currently, gas is not handled */ 237 radcoefs_min = RADCOEFS_NULL; 238 radcoefs_max = RADCOEFS_NULL; 239 voxel_commit_radcoefs_range 240 (vox, ATRSTM_CPNT_GAS, &radcoefs_min, &radcoefs_max); 241 242 /* Compute the range of the radiative coefficient for the current 243 * primitive */ 244 #ifndef ATRSTM_USE_SIMD 245 { 246 res = primitive_compute_radcoefs_range 247 (atrstm, refract_id, &prim, &radcoefs_min, &radcoefs_max); 248 ASSERT(atrstm->use_simd == 0); 249 } 250 #else 251 { 252 if(atrstm->use_simd) { 253 res = primitive_compute_radcoefs_range_simd4 254 (atrstm, refract_id, &prim, &radcoefs_min, &radcoefs_max); 255 } else { 256 res = primitive_compute_radcoefs_range 257 (atrstm, refract_id, &prim, &radcoefs_min, &radcoefs_max); 258 } 259 } 260 #endif 261 if(UNLIKELY(res != RES_OK)) goto error; 262 voxel_commit_radcoefs_range 263 (vox, ATRSTM_CPNT_SOOT, &radcoefs_min, &radcoefs_max); 264 } 265 } 266 } 267 } 268 } 269 270 exit: 271 return res; 272 error: 273 goto exit; 274 } 275 276 static res_T 277 voxelize_volumetric_mesh 278 (struct atrstm* atrstm, 279 const struct atrri_refractive_index* refract_id, 280 struct pool* pool) 281 { 282 struct darray_prims_list prims_list; /* Per thread list of primitives */ 283 const size_t part_def = PARTITION_DEFINITION; 284 size_t nparts[3]; /* #partitions along the 3 axis */ 285 size_t nparts_adjusted; 286 double vol_low[3]; 287 double vol_upp[3]; 288 double vxsz[3]; 289 int64_t i; 290 int progress = 0; 291 ATOMIC nparts_voxelized = 0; 292 ATOMIC res = RES_OK; 293 ASSERT(atrstm && pool); 294 295 darray_prims_list_init(atrstm->allocator, &prims_list); 296 297 /* Allocate the per thread primitive lists */ 298 res = darray_prims_list_resize(&prims_list, atrstm->nthreads); 299 if(res != RES_OK) goto error; 300 301 /* Fetch the volume AABB and compute the size of one voxel */ 302 SUVM(volume_get_aabb(atrstm->volume, vol_low, vol_upp)); 303 vxsz[0] = (vol_upp[0] - vol_low[0]) / (double)atrstm->grid_max_definition[0]; 304 vxsz[1] = (vol_upp[1] - vol_low[1]) / (double)atrstm->grid_max_definition[1]; 305 vxsz[2] = (vol_upp[2] - vol_low[2]) / (double)atrstm->grid_max_definition[2]; 306 307 /* Compute the number of partitions */ 308 nparts[0] = (atrstm->grid_max_definition[0] + (part_def-1)/*ceil*/) / part_def; 309 nparts[1] = (atrstm->grid_max_definition[1] + (part_def-1)/*ceil*/) / part_def; 310 nparts[2] = (atrstm->grid_max_definition[2] + (part_def-1)/*ceil*/) / part_def; 311 312 /* Compute the overall #partitions allowing Morton indexing */ 313 nparts_adjusted = MMAX(nparts[0], MMAX(nparts[1], nparts[2])); 314 nparts_adjusted = round_up_pow2(nparts_adjusted); 315 nparts_adjusted = nparts_adjusted*nparts_adjusted*nparts_adjusted; 316 317 #define LOG_MSG "Voxelizing the volumetric mesh: %3d%%\r" 318 log_info(atrstm, LOG_MSG, 0); 319 320 /* Iterate over the partition according to their Morton order and voxelize 321 * the primitives that intersect it */ 322 omp_set_num_threads((int)atrstm->nthreads); 323 #pragma omp parallel for schedule(static, 1/*chunk size*/) 324 for(i = 0; i < (int64_t)nparts_adjusted; ++i) { 325 struct darray_size_t* prims = NULL; 326 double part_low[3]; 327 double part_upp[3]; 328 uint32_t part_ids[3]; 329 const int ithread = omp_get_thread_num(); 330 struct part* part = NULL; 331 int pcent; 332 size_t n; 333 res_T res_local = RES_OK; 334 335 /* Handle error */ 336 if(ATOMIC_GET(&res) != RES_OK) continue; 337 338 /* Fetch the per thread list of primitives and partition pool */ 339 prims = darray_prims_list_data_get(&prims_list)+ithread; 340 darray_size_t_clear(prims); 341 342 part = pool_next_partition(pool); 343 if(!part) { /* An error occurs */ 344 ATOMIC_SET(&res, RES_UNKNOWN_ERR); 345 continue; 346 } 347 348 morton_xyz_decode_u21((uint64_t)part->id, part_ids); 349 350 /* Check that the partition is not out of bound due to Morton indexing */ 351 if(part_ids[0] >= nparts[0] 352 || part_ids[1] >= nparts[1] 353 || part_ids[2] >= nparts[2]) { 354 pool_free_partition(pool, part); 355 continue; 356 } 357 358 /* Compute the partition AABB */ 359 part_low[0] = (double)(part_ids[0] * part_def) * vxsz[0] + vol_low[0]; 360 part_low[1] = (double)(part_ids[1] * part_def) * vxsz[1] + vol_low[1]; 361 part_low[2] = (double)(part_ids[2] * part_def) * vxsz[2] + vol_low[2]; 362 part_upp[0] = part_low[0] + (double)part_def * vxsz[0]; 363 part_upp[1] = part_low[1] + (double)part_def * vxsz[1]; 364 part_upp[2] = part_low[2] + (double)part_def * vxsz[2]; 365 366 /* Find the list of primitives that intersects the partition */ 367 res_local = suvm_volume_intersect_aabb(atrstm->volume, part_low, part_upp, 368 register_primitive, prims); 369 if(res_local != RES_OK) { 370 ATOMIC_SET(&res, res_local); 371 continue; 372 } 373 374 res = voxelize_partition 375 (atrstm, prims, part_low, part_upp, vxsz, refract_id, part); 376 if(res != RES_OK) { 377 pool_free_partition(pool, part); 378 ATOMIC_SET(&res, res_local); 379 continue; 380 } 381 382 pool_commit_partition(pool, part); 383 384 /* Update progress message */ 385 n = (size_t)ATOMIC_INCR(&nparts_voxelized); 386 pcent = (int)((n * 100) / (nparts[0]*nparts[1]*nparts[2])); 387 388 #pragma omp critical 389 if(pcent > progress) { 390 progress = pcent; 391 log_info(atrstm, LOG_MSG, pcent); 392 } 393 } 394 395 if(res != RES_OK) goto error; 396 log_info(atrstm, LOG_MSG"\n", 100); 397 398 exit: 399 darray_prims_list_release(&prims_list); 400 return (res_T)res; 401 error: 402 goto exit; 403 } 404 405 static void 406 voxel_get(const size_t xyz[3], const uint64_t mcode, void* dst, void* context) 407 { 408 struct build_octree_context* ctx = context; 409 float* vox = NULL; 410 uint64_t ivox, ipart; 411 ASSERT(xyz && dst && ctx); 412 (void)xyz; 413 414 /* Retrieve the partition morton id and the per partition voxel morton id 415 * from the overall voxel morton code. As voxels, partitions are sorted in 416 * morton order and thus the partition morton ID is encoded in the MSB of the 417 * morton code while the voxels morton ID is stored in its LSB. */ 418 ipart = (mcode >> (LOG2_PARTITION_DEFINITION*3)); 419 ivox = (mcode & (BIT_U64(LOG2_PARTITION_DEFINITION*3)-1)); 420 421 if(!ctx->part || ctx->part->id != ipart) { 422 if(ctx->part) { /* Free the previous partition */ 423 pool_free_partition(ctx->pool, ctx->part); 424 } 425 /* Fetch the partition storing the voxel */ 426 ctx->part = pool_fetch_partition(ctx->pool, ipart); 427 } 428 429 if(!ctx->part) { /* An error occurs */ 430 memset(dst, 0, NFLOATS_PER_VOXEL * sizeof(float)); 431 } else { 432 int cpnt; 433 vox = part_get_voxel(ctx->part, ivox); /* Fetch the voxel data */ 434 435 /* Setup destination data */ 436 FOR_EACH(cpnt, 0, ATRSTM_CPNTS_COUNT__) { 437 int radcoef; 438 FOR_EACH(radcoef, 0, ATRSTM_RADCOEFS_COUNT__) { 439 const float k_min = vox_get(vox, cpnt, radcoef, ATRSTM_SVX_OP_MIN); 440 const float k_max = vox_get(vox, cpnt, radcoef, ATRSTM_SVX_OP_MAX); 441 ASSERT(ATRSTM_SVX_OPS_COUNT__ == 2); 442 443 if(k_min <= k_max) { 444 vox_set(dst, cpnt, radcoef, ATRSTM_SVX_OP_MIN, k_min); 445 vox_set(dst, cpnt, radcoef, ATRSTM_SVX_OP_MAX, k_max); 446 } else { 447 /* The voxel was not overlaped by any tetrahedron. Set its radiative 448 * coefficients to 0 */ 449 vox_set(dst, cpnt, radcoef, ATRSTM_SVX_OP_MIN, 0); 450 vox_set(dst, cpnt, radcoef, ATRSTM_SVX_OP_MAX, 0); 451 } 452 } 453 } 454 } 455 } 456 457 static void 458 voxels_merge_component 459 (void* dst, 460 const enum atrstm_component cpnt, 461 const void* voxels[], 462 const size_t nvoxs) 463 { 464 float ka[2] = {FLT_MAX,-FLT_MAX}; 465 float ks[2] = {FLT_MAX,-FLT_MAX}; 466 float kext[2] = {FLT_MAX,-FLT_MAX}; 467 size_t ivox; 468 ASSERT(dst && voxels && nvoxs); 469 470 FOR_EACH(ivox, 0, nvoxs) { 471 const float* vox = voxels[ivox]; 472 ka[0] = MMIN(ka[0], vox_get(vox, cpnt, ATRSTM_RADCOEF_Ka, ATRSTM_SVX_OP_MIN)); 473 ka[1] = MMAX(ka[1], vox_get(vox, cpnt, ATRSTM_RADCOEF_Ka, ATRSTM_SVX_OP_MAX)); 474 ks[0] = MMIN(ks[0], vox_get(vox, cpnt, ATRSTM_RADCOEF_Ks, ATRSTM_SVX_OP_MIN)); 475 ks[1] = MMAX(ks[1], vox_get(vox, cpnt, ATRSTM_RADCOEF_Ks, ATRSTM_SVX_OP_MAX)); 476 kext[0] = MMIN(kext[0], vox_get(vox, cpnt, ATRSTM_RADCOEF_Kext, ATRSTM_SVX_OP_MIN)); 477 kext[1] = MMAX(kext[1], vox_get(vox, cpnt, ATRSTM_RADCOEF_Kext, ATRSTM_SVX_OP_MAX)); 478 } 479 vox_set(dst, cpnt, ATRSTM_RADCOEF_Ka, ATRSTM_SVX_OP_MIN, ka[0]); 480 vox_set(dst, cpnt, ATRSTM_RADCOEF_Ka, ATRSTM_SVX_OP_MAX, ka[1]); 481 vox_set(dst, cpnt, ATRSTM_RADCOEF_Ks, ATRSTM_SVX_OP_MIN, ks[0]); 482 vox_set(dst, cpnt, ATRSTM_RADCOEF_Ks, ATRSTM_SVX_OP_MAX, ks[1]); 483 vox_set(dst, cpnt, ATRSTM_RADCOEF_Kext, ATRSTM_SVX_OP_MIN, kext[0]); 484 vox_set(dst, cpnt, ATRSTM_RADCOEF_Kext, ATRSTM_SVX_OP_MAX, kext[1]); 485 } 486 487 static void 488 voxels_merge 489 (void* dst, 490 const void* voxels[], 491 const size_t nvoxs, 492 void* ctx) 493 { 494 (void)ctx; 495 voxels_merge_component(dst, ATRSTM_CPNT_GAS, voxels, nvoxs); 496 voxels_merge_component(dst, ATRSTM_CPNT_SOOT, voxels, nvoxs); 497 } 498 499 static INLINE void 500 voxels_component_compute_kext_range 501 (const enum atrstm_component cpnt, 502 const struct svx_voxel voxels[], 503 const size_t nvoxs, 504 float kext[2]) 505 { 506 size_t ivox; 507 ASSERT(voxels && nvoxs && kext); 508 509 kext[0] = FLT_MAX; 510 kext[1] =-FLT_MAX; 511 512 /* Compute the Kext range of the submitted voxels */ 513 FOR_EACH(ivox, 0, nvoxs) { 514 const float* vox = voxels[ivox].data; 515 kext[0] = MMIN(kext[0], vox_get(vox, cpnt, ATRSTM_RADCOEF_Kext, ATRSTM_SVX_OP_MIN)); 516 kext[1] = MMAX(kext[1], vox_get(vox, cpnt, ATRSTM_RADCOEF_Kext, ATRSTM_SVX_OP_MAX)); 517 } 518 } 519 520 static int 521 voxels_challenge_merge 522 (const struct svx_voxel voxels[], 523 const size_t nvoxs, 524 void* context) 525 { 526 const struct build_octree_context* ctx = context; 527 double lower[3] = { DBL_MAX, DBL_MAX, DBL_MAX}; 528 double upper[3] = {-DBL_MAX,-DBL_MAX,-DBL_MAX}; 529 double sz_max; 530 size_t ivox; 531 float kext_gas[2] = {0, 0}; 532 float kext_soot[2] = {0, 0}; 533 int merge_gas = 0; 534 int merge_soot = 0; 535 ASSERT(voxels && nvoxs && context); 536 537 /* Compute the AABB of the group of voxels and its maximum size */ 538 FOR_EACH(ivox, 0, nvoxs) { 539 lower[0] = MMIN(voxels[ivox].lower[0], lower[0]); 540 lower[1] = MMIN(voxels[ivox].lower[1], lower[1]); 541 lower[2] = MMIN(voxels[ivox].lower[2], lower[2]); 542 upper[0] = MMAX(voxels[ivox].upper[0], upper[0]); 543 upper[1] = MMAX(voxels[ivox].upper[1], upper[1]); 544 upper[2] = MMAX(voxels[ivox].upper[2], upper[2]); 545 } 546 sz_max = upper[0] - lower[0]; 547 sz_max = MMAX(upper[1] - lower[1], sz_max); 548 sz_max = MMAX(upper[2] - lower[2], sz_max); 549 550 /* Compute the Kext range for each component */ 551 voxels_component_compute_kext_range(ATRSTM_CPNT_GAS, voxels, nvoxs, kext_gas); 552 voxels_component_compute_kext_range(ATRSTM_CPNT_SOOT, voxels, nvoxs, kext_soot); 553 554 if(kext_gas[0] > kext_gas[1]) { /* Empty voxels */ 555 /* Always merge empty voxels */ 556 merge_gas = 1; 557 merge_soot = 1; 558 } else { 559 /* Check if the voxels are candidates to the merge process by evaluating its 560 * optical thickness regarding the maximum size of their AABB and a user 561 * defined threshold */ 562 ASSERT(kext_gas[0] <= kext_gas[1]); 563 ASSERT(kext_soot[0] <= kext_soot[1]); 564 merge_gas = (kext_gas[1] - kext_gas[0])*sz_max <= ctx->tau_threshold; 565 merge_soot = (kext_soot[1] - kext_soot[0])*sz_max <= ctx->tau_threshold; 566 } 567 return merge_gas && merge_soot; 568 } 569 570 static res_T 571 build_octree 572 (struct atrstm* atrstm, 573 struct pool* pool) 574 { 575 struct svx_voxel_desc vox_desc = SVX_VOXEL_DESC_NULL; 576 struct build_octree_context ctx = BUILD_OCTREE_CONTEXT_NULL; 577 double lower[3]; 578 double upper[3]; 579 size_t definition[3]; 580 res_T res = RES_OK; 581 ASSERT(atrstm && pool); 582 583 /* Setup the build octree context */ 584 ctx.atrstm = atrstm; 585 ctx.pool = pool; 586 ctx.tau_threshold = atrstm->optical_thickness; 587 588 /* Setup the voxel descriptor. TODO in shortwave, the NFLOATS_PER_VOXEL 589 * could be divided by 2 since there is no gas to handle */ 590 vox_desc.get = voxel_get; 591 vox_desc.merge = voxels_merge; 592 vox_desc.challenge_merge = voxels_challenge_merge; 593 vox_desc.context = &ctx; 594 vox_desc.size = NFLOATS_PER_VOXEL * sizeof(float); 595 596 SUVM(volume_get_aabb(atrstm->volume, lower, upper)); 597 598 /* Create the octree */ 599 definition[0] = (size_t)atrstm->grid_max_definition[0]; 600 definition[1] = (size_t)atrstm->grid_max_definition[1]; 601 definition[2] = (size_t)atrstm->grid_max_definition[2]; 602 res = svx_octree_create 603 (atrstm->svx, lower, upper, definition, &vox_desc, &atrstm->octree); 604 if(res != RES_OK) goto error; 605 606 exit: 607 return res; 608 error: 609 goto exit; 610 } 611 612 static res_T 613 build_octrees(struct atrstm* atrstm) 614 { 615 struct atrri_refractive_index refract_id = ATRRI_REFRACTIVE_INDEX_NULL; 616 struct pool pool; 617 double wlen; 618 int pool_is_init = 0; 619 ATOMIC res = RES_OK; 620 ASSERT(atrstm && check_octrees_parameters(atrstm)); 621 622 /* Currently only shortwave is handled and in this case, radiative 623 * transfert is monochromatic */ 624 ASSERT(atrstm->spectral_type == ATRSTM_SPECTRAL_SW); 625 ASSERT(atrstm->wlen_range[0] == atrstm->wlen_range[1]); 626 wlen = atrstm->wlen_range[0]; 627 628 /* Fetch the submitted wavelength */ 629 res = atrri_fetch_refractive_index(atrstm->atrri, wlen, &refract_id); 630 if(res != RES_OK) goto error; 631 632 /* Initialise the pool of partitions */ 633 res = pool_init(atrstm->allocator, 16 * atrstm->nthreads, &pool); 634 if(res != RES_OK) { 635 log_err(atrstm, 636 "Error initializing the pool of voxel partitions -- %s.\n", 637 res_to_cstr((res_T)res)); 638 goto error; 639 } 640 pool_is_init = 1; 641 642 log_info(atrstm, 643 "Evaluate and partition the field of optical properties.\n"); 644 log_info(atrstm, 645 "Grid definition: %ux%ux%u.\n", SPLIT3(atrstm->grid_max_definition)); 646 647 omp_set_nested(1); /* Enable nested threads for voxelize_volumetric_mesh */ 648 #pragma omp parallel sections num_threads(2) 649 { 650 #pragma omp section 651 { 652 const res_T res_local = voxelize_volumetric_mesh 653 (atrstm, &refract_id, &pool); 654 if(res_local != RES_OK) { 655 log_err(atrstm, "Error voxelizing the volumetric mesh -- %s\n", 656 res_to_cstr(res_local)); 657 pool_invalidate(&pool); 658 ATOMIC_SET(&res, res_local); 659 } 660 } 661 662 #pragma omp section 663 { 664 const res_T res_local = build_octree(atrstm, &pool); 665 if(res_local != RES_OK) { 666 log_err(atrstm, "Error building the octree -- %s\n", 667 res_to_cstr(res_local)); 668 pool_invalidate(&pool); 669 ATOMIC_SET(&res, res_local); 670 } 671 } 672 } 673 if(res != RES_OK) goto error; 674 675 if(atrstm->cache) { 676 ASSERT(cache_is_empty(atrstm->cache)); 677 678 log_info(atrstm, "Write acceleration data structure to '%s'.\n", 679 cache_get_name(atrstm->cache)); 680 681 res = svx_tree_write(atrstm->octree, cache_get_stream(atrstm->cache)); 682 if(res != RES_OK) { 683 log_err(atrstm, "Could not write the octrees to the cache -- %s.\n", 684 res_to_cstr((res_T)res)); 685 goto error; 686 } 687 688 fflush(cache_get_stream(atrstm->cache)); 689 } 690 691 exit: 692 if(pool_is_init) pool_release(&pool); 693 return (res_T)res; 694 error: 695 goto exit; 696 } 697 698 static res_T 699 load_octrees(struct atrstm* atrstm) 700 { 701 FILE* stream = NULL; 702 res_T res = RES_OK; 703 ASSERT(atrstm && atrstm->cache && !cache_is_empty(atrstm->cache)); 704 705 log_info(atrstm, "Load the acceleration data structures of the optical " 706 "properties from '%s'.\n", cache_get_name(atrstm->cache)); 707 708 stream = cache_get_stream(atrstm->cache); 709 res = svx_tree_create_from_stream(atrstm->svx, stream, &atrstm->octree); 710 if(res != RES_OK) { 711 log_err(atrstm, "Could not load the acceleration data structures -- %s.\n", 712 res_to_cstr(res)); 713 goto error; 714 } 715 716 exit: 717 return res; 718 error: 719 goto exit; 720 } 721 722 /******************************************************************************* 723 * Local functions 724 ******************************************************************************/ 725 res_T 726 setup_octrees(struct atrstm* atrstm) 727 { 728 char buf[128]; 729 struct time t0, t1; 730 size_t sz; 731 res_T res = RES_OK; 732 ASSERT(atrstm); 733 734 time_current(&t0); 735 736 if(atrstm->cache && !cache_is_empty(atrstm->cache)) { 737 res = load_octrees(atrstm);; 738 if(res != RES_OK) goto error; 739 } else { 740 res = build_octrees(atrstm); 741 if(res != RES_OK) goto error; 742 } 743 744 time_sub(&t0, time_current(&t1), &t0); 745 time_dump(&t0, TIME_ALL, NULL, buf, sizeof(buf)); 746 log_info(atrstm, "Setup the partitionning data structures in %s.\n", buf); 747 748 sz = MEM_ALLOCATED_SIZE(&atrstm->svx_allocator); 749 dump_memory_size(sz, NULL, buf, sizeof(buf)); 750 log_info(atrstm, "Star-VoXel memory usage: %s.\n", buf); 751 752 exit: 753 return res; 754 error: 755 goto exit; 756 } 757 758 void 759 octrees_clean(struct atrstm* atrstm) 760 { 761 ASSERT(atrstm); 762 if(atrstm->octree) SVX(tree_ref_put(atrstm->octree)); 763 }