commit 82bc303fbdd6a1ee854e9ca7e8205669f5d233c6
parent 95a8cd9eaa5e383b3f788fe1e1153d0e0d071997
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Wed, 20 Jan 2021 13:56:49 +0100
Add cache support to the setup_octrees procedure
Diffstat:
9 files changed, 767 insertions(+), 732 deletions(-)
diff --git a/cmake/CMakeLists.txt b/cmake/CMakeLists.txt
@@ -54,19 +54,20 @@ set(VERSION ${VERSION_MAJOR}.${VERSION_MINOR}.${VERSION_PATCH})
set(ATRSTM_FILES_SRC
atrstm.c
- atrstm_build_octrees.c
atrstm_cache.c
atrstm_log.c
atrstm_optprops.c
atrstm_partition.c
+ atrstm_setup_octrees.c
atrstm_setup_uvm.c)
set(ATRSTM_FILES_INC
- atrstm_build_octrees.h
atrstm_c.h
atrstm_cache.h
atrstm_log.h
atrstm_optprops.h
- atrstm_partition.h)
+ atrstm_partition.h
+ atrstm_setup_octrees.h)
+
set(ATRSTM_FILES_INC_API
atrstm.h)
diff --git a/src/atrstm.c b/src/atrstm.c
@@ -16,10 +16,10 @@
#define _POSIX_C_SOURCE 200112L /* snprintf support */
#include "atrstm.h"
-#include "atrstm_build_octrees.h"
#include "atrstm_c.h"
#include "atrstm_cache.h"
#include "atrstm_log.h"
+#include "atrstm_setup_octrees.h"
#include <astoria/atrri.h>
#include <astoria/atrtp.h>
@@ -84,7 +84,6 @@ atrstm_create
const struct atrstm_args* args,
struct atrstm** out_atrstm)
{
- struct build_octrees_args build_octree_args = BUILD_OCTREES_ARGS_NULL;
struct atrstm* atrstm = NULL;
struct mem_allocator* allocator = NULL;
int nthreads_max;
@@ -120,6 +119,11 @@ atrstm_create
atrstm->spectral_type = args->spectral_type;
atrstm->wlen_range[0] = args->wlen_range[0];
atrstm->wlen_range[1] = args->wlen_range[1];
+ atrstm->grid_max_definition[0] = args->grid_max_definition[0];
+ atrstm->grid_max_definition[1] = args->grid_max_definition[1];
+ atrstm->grid_max_definition[2] = args->grid_max_definition[2];
+ atrstm->optical_thickness = args->optical_thickness;
+
str_init(allocator, &atrstm->name);
if(logger) {
@@ -191,11 +195,7 @@ atrstm_create
}
/* Build the octree */
- build_octree_args.grid_max_definition[0] = args->grid_max_definition[0];
- build_octree_args.grid_max_definition[1] = args->grid_max_definition[1];
- build_octree_args.grid_max_definition[2] = args->grid_max_definition[2];
- build_octree_args.optical_thickness_threshold = args->optical_thickness;
- res = build_octrees(atrstm, &build_octree_args);
+ res = setup_octrees(atrstm);
if(res != RES_OK) goto error;
exit:
diff --git a/src/atrstm_build_octrees.c b/src/atrstm_build_octrees.c
@@ -1,669 +0,0 @@
-/* Copyright (C) 2020 CNRS
- *
- * This program is free software: you can redistribute it and/or modify
- * it under the terms of the GNU General Public License as published by
- * the Free Software Foundation, either version 3 of the License, or
- * (at your option) any later version.
- *
- * This program is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- * GNU General Public License for more details.
- *
- * You should have received a copy of the GNU General Public License
- * along with this program. If not, see <http://www.gnu.org/licenses/>. */
-
-#define _POSIX_C_SOURCE 200112L /* nextafterf */
-
-#include "atrstm_c.h"
-#include "atrstm_build_octrees.h"
-#include "atrstm_log.h"
-#include "atrstm_optprops.h"
-#include "atrstm_partition.h"
-#include "atrstm_svx.h"
-
-#include <astoria/atrri.h>
-
-#include <star/suvm.h>
-#include <star/svx.h>
-
-#include <rsys/clock_time.h>
-#include <rsys/cstr.h>
-#include <rsys/dynamic_array_size_t.h>
-#include <rsys/morton.h>
-
-#include <math.h>
-#include <omp.h>
-
-struct build_octree_context {
- struct atrstm* atrstm;
- struct pool* pools; /* Per thread pool of voxel partitions */
-
- /* Optical thickness threshold criteria for the merge process */
- double tau_threshold;
-};
-#define BUILD_OCTREE_CONTEXT_NULL__ { NULL, NULL, 0 }
-static const struct build_octree_context BUILD_OCTREE_CONTEXT_NULL =
- BUILD_OCTREE_CONTEXT_NULL__;
-
-/* Define the darray of list of primitive ids */
-#define DARRAY_NAME prims_list
-#define DARRAY_DATA struct darray_size_t
-#define DARRAY_FUNCTOR_INIT darray_size_t_init
-#define DARRAY_FUNCTOR_RELEASE darray_size_t_release
-#define DARRAY_FUNCTOR_COPY darray_size_t_copy
-#define DARRAY_FUNCTOR_COPY_AND_RELEASE darray_size_t_copy_and_release
-#include <rsys/dynamic_array.h>
-
-/*******************************************************************************
- * Helper functions
- ******************************************************************************/
-static INLINE int
-check_build_octrees_args(const struct build_octrees_args* args)
-{
- return args
- && args->grid_max_definition[0]
- && args->grid_max_definition[1]
- && args->grid_max_definition[2]
- && args->optical_thickness_threshold >= 0
- && (!args->cache.stream || args->cache.name);
-}
-
-static INLINE void
-register_primitive
- (const struct suvm_primitive* prim,
- const double low[3],
- const double upp[3],
- void* context)
-{
- struct darray_size_t* prims = context;
- ASSERT(prim && low && upp && context);
- ASSERT(low[0] < upp[0]);
- ASSERT(low[1] < upp[1]);
- ASSERT(low[2] < upp[2]);
- (void)low, (void)upp;
- CHK(darray_size_t_push_back(prims, &prim->iprim) == RES_OK);
-}
-
-static void
-voxel_commit_optprops_range
- (float* vox,
- const enum atrstm_component cpnt,
- const struct optprops* optprops_min,
- const struct optprops* optprops_max)
-{
- double ka[2], ks[2], kext[2];
- float vox_ka[2], vox_ks[2], vox_kext[2];
- ASSERT(vox);
- ASSERT(optprops_min);
- ASSERT(optprops_max);
-
- ka[0] = optprops_min->ka;
- ka[1] = optprops_max->ka;
- ks[0] = optprops_min->ks;
- ks[1] = optprops_max->ks;
- kext[0] = optprops_min->kext;
- kext[1] = optprops_max->kext;
-
- /* Ensure that the single precision bounds include their double precision
- * version. */
- if(ka[0] != (float)ka[0]) ka[0] = nextafterf((float)ka[0],-FLT_MAX);
- if(ka[1] != (float)ka[1]) ka[1] = nextafterf((float)ka[1], FLT_MAX);
- if(ks[0] != (float)ks[0]) ks[0] = nextafterf((float)ks[0],-FLT_MAX);
- if(ks[1] != (float)ks[1]) ks[1] = nextafterf((float)ks[1], FLT_MAX);
- if(kext[0] != (float)kext[0]) kext[0] = nextafterf((float)kext[0],-FLT_MAX);
- if(kext[1] != (float)kext[1]) kext[1] = nextafterf((float)kext[1], FLT_MAX);
-
- /* Fetch the range of the voxel optical properties */
- vox_ka[0] = vox_get(vox, cpnt, ATRSTM_Ka, ATRSTM_SVX_MIN);
- vox_ka[1] = vox_get(vox, cpnt, ATRSTM_Ka, ATRSTM_SVX_MAX);
- vox_ks[0] = vox_get(vox, cpnt, ATRSTM_Ks, ATRSTM_SVX_MIN);
- vox_ks[1] = vox_get(vox, cpnt, ATRSTM_Ks, ATRSTM_SVX_MAX);
- vox_kext[0] = vox_get(vox, cpnt, ATRSTM_Kext, ATRSTM_SVX_MIN);
- vox_kext[1] = vox_get(vox, cpnt, ATRSTM_Kext, ATRSTM_SVX_MAX);
-
- /* Update the range of the voxel optical properties */
- vox_ka[0] = MMIN(vox_ka[0], (float)ka[0]);
- vox_ka[1] = MMAX(vox_ka[1], (float)ka[1]);
- vox_ks[0] = MMIN(vox_ks[0], (float)ks[0]);
- vox_ks[1] = MMAX(vox_ks[1], (float)ks[1]);
- vox_kext[0] = MMIN(vox_kext[0], (float)kext[0]);
- vox_kext[1] = MMAX(vox_kext[1], (float)kext[1]);
-
- /* Commit the upd to the voxel */
- vox_set(vox, cpnt, ATRSTM_Ka, ATRSTM_SVX_MIN, vox_ka[0]);
- vox_set(vox, cpnt, ATRSTM_Ka, ATRSTM_SVX_MAX, vox_ka[1]);
- vox_set(vox, cpnt, ATRSTM_Ks, ATRSTM_SVX_MIN, vox_ks[0]);
- vox_set(vox, cpnt, ATRSTM_Ks, ATRSTM_SVX_MAX, vox_ks[1]);
- vox_set(vox, cpnt, ATRSTM_Kext, ATRSTM_SVX_MIN, vox_kext[0]);
- vox_set(vox, cpnt, ATRSTM_Kext, ATRSTM_SVX_MAX, vox_kext[1]);
-}
-
-static res_T
-voxelize_partition
- (struct atrstm* atrstm,
- const struct darray_size_t* prims, /* Primitives intersecting the partition */
- const double part_low[3], /* Lower bound of the partition */
- const double part_upp[3], /* Upper bound of the partition */
- const double vxsz[3], /* Size of a partition voxel */
- const struct atrri_refractive_index* refract_id, /* Refractive index */
- struct part* part)
-{
- size_t iprim;
- res_T res = RES_OK;
- ASSERT(atrstm && prims && part_low && part_upp && vxsz && part && refract_id);
- ASSERT(part_low[0] < part_upp[0]);
- ASSERT(part_low[1] < part_upp[1]);
- ASSERT(part_low[2] < part_upp[2]);
- ASSERT(vxsz[0] > 0 && vxsz[1] > 0 && vxsz[2] > 0);
-
- /* Clean the partition voxels */
- part_clear_voxels(part);
-
- FOR_EACH(iprim, 0, darray_size_t_size_get(prims)) {
- struct suvm_primitive prim = SUVM_PRIMITIVE_NULL;
- struct suvm_polyhedron poly;
- double poly_low[3];
- double poly_upp[3];
- uint32_t ivxl_low[3];
- uint32_t ivxl_upp[3];
- uint32_t ivxl[3];
- size_t prim_id;
-
- prim_id = darray_size_t_cdata_get(prims)[iprim];
-
- /* Retrieve the primitive data and setup its polyhedron */
- SUVM(volume_get_primitive(atrstm->volume, prim_id, &prim));
- SUVM(primitive_setup_polyhedron(&prim, &poly));
- ASSERT(poly.lower[0] <= part_upp[0] && poly.upper[0] >= part_low[0]);
- ASSERT(poly.lower[1] <= part_upp[1] && poly.upper[1] >= part_low[1]);
- ASSERT(poly.lower[2] <= part_upp[2] && poly.upper[2] >= part_low[2]);
-
- /* Clamp the poly AABB to the partition boundaries */
- poly_low[0] = MMAX(poly.lower[0], part_low[0]);
- poly_low[1] = MMAX(poly.lower[1], part_low[1]);
- poly_low[2] = MMAX(poly.lower[2], part_low[2]);
- poly_upp[0] = MMIN(poly.upper[0], part_upp[0]);
- poly_upp[1] = MMIN(poly.upper[1], part_upp[1]);
- poly_upp[2] = MMIN(poly.upper[2], part_upp[2]);
-
- /* Transform the polyhedron AABB in partition voxel space */
- ivxl_low[0] = (uint32_t)((poly_low[0] - part_low[0]) / vxsz[0]);
- ivxl_low[1] = (uint32_t)((poly_low[1] - part_low[1]) / vxsz[1]);
- ivxl_low[2] = (uint32_t)((poly_low[2] - part_low[2]) / vxsz[2]);
- ivxl_upp[0] = (uint32_t)ceil((poly_upp[0] - part_low[0]) / vxsz[0]);
- ivxl_upp[1] = (uint32_t)ceil((poly_upp[1] - part_low[1]) / vxsz[1]);
- ivxl_upp[2] = (uint32_t)ceil((poly_upp[2] - part_low[2]) / vxsz[2]);
- ASSERT(ivxl_upp[0] <= PARTITION_DEFINITION);
- ASSERT(ivxl_upp[0] <= PARTITION_DEFINITION);
- ASSERT(ivxl_upp[0] <= PARTITION_DEFINITION);
-
- /* Iterate over the partition voxels intersected by the primitive AABB */
- FOR_EACH(ivxl[2], ivxl_low[2], ivxl_upp[2]) {
- float vxl_low[3]; /* Voxel lower bound */
- float vxl_upp[3]; /* Voxel upper bound */
- uint64_t mcode[3]; /* Morton code cache */
-
- vxl_low[2] = (float)((double)ivxl[2] * vxsz[2] + part_low[2]);
- vxl_upp[2] = vxl_low[2] + (float)vxsz[2];
- mcode[2] = morton3D_encode_u21(ivxl[2]) << 0;
-
- FOR_EACH(ivxl[1], ivxl_low[1], ivxl_upp[1]) {
- vxl_low[1] = (float)((double)ivxl[1] * vxsz[1] + part_low[1]);
- vxl_upp[1] = vxl_low[1] + (float)vxsz[1];
- mcode[1] = (morton3D_encode_u21(ivxl[1]) << 1) | mcode[2];
-
- FOR_EACH(ivxl[0], ivxl_low[0], ivxl_upp[0]) {
- enum suvm_intersection_type intersect;
- vxl_low[0] = (float)((double)ivxl[0] * vxsz[0] + part_low[0]);
- vxl_upp[0] = vxl_low[0] + (float)vxsz[0];
-
- /* NOTE mcode[0] store the morton index of the voxel */
- mcode[0] = (morton3D_encode_u21(ivxl[0]) << 2) | mcode[1];
-
- intersect = suvm_polyhedron_intersect_aabb(&poly, vxl_low, vxl_upp);
- if(intersect != SUVM_INTERSECT_NONE) {
- struct optprops optprops_min;
- struct optprops optprops_max;
- float* vox = part_get_voxel(part, mcode[2]);
-
- /* Currently, gas is not handled */
- optprops_min = OPTPROPS_NULL;
- optprops_max = OPTPROPS_NULL;
- voxel_commit_optprops_range
- (vox, ATRSTM_CPNT_GAS, &optprops_min, &optprops_max);
-
- res = primitive_compute_optprops_range
- (atrstm, refract_id, &prim, &optprops_min, &optprops_max);
- if(UNLIKELY(res != RES_OK)) goto error;
- voxel_commit_optprops_range
- (vox, ATRSTM_CPNT_SOOT, &optprops_min, &optprops_max);
- }
- }
- }
- }
- }
-
-exit:
- return res;
-error:
- goto exit;
-}
-
-static res_T
-voxelize_volumetric_mesh
- (struct atrstm* atrstm,
- const struct build_octrees_args* args,
- const struct atrri_refractive_index* refract_id,
- struct pool* pools)
-{
- struct darray_prims_list prims_list; /* Per thread list of primitives */
- const size_t part_def = PARTITION_DEFINITION;
- size_t nparts[3]; /* #partitions along the 3 axis */
- size_t nparts_adjusted;
- double vol_low[3];
- double vol_upp[3];
- double vxsz[3];
- uint64_t ipart;
- ATOMIC res = RES_OK;
- ASSERT(atrstm && args && pools);
-
- darray_prims_list_init(atrstm->allocator, &prims_list);
-
- /* Allocate the per thread primitive lists */
- res = darray_prims_list_resize(&prims_list, atrstm->nthreads);
- if(res != RES_OK) goto error;
-
- /* Fetch the volume AABB and compute the size of one voxel */
- SUVM(volume_get_aabb(atrstm->volume, vol_low, vol_upp));
- vxsz[0] = (vol_upp[0] - vol_low[0]) / (double)args->grid_max_definition[0];
- vxsz[1] = (vol_upp[1] - vol_low[1]) / (double)args->grid_max_definition[1];
- vxsz[2] = (vol_upp[2] - vol_low[2]) / (double)args->grid_max_definition[2];
-
- /* Compute the number of partitions */
- nparts[0] = (args->grid_max_definition[0] + (part_def-1)/*ceil*/) / part_def;
- nparts[1] = (args->grid_max_definition[1] + (part_def-1)/*ceil*/) / part_def;
- nparts[2] = (args->grid_max_definition[2] + (part_def-1)/*ceil*/) / part_def;
-
- /* Compute the overall #partitions allowing Morton indexing */
- nparts_adjusted = MMAX(nparts[0], MMAX(nparts[1], nparts[2]));
- nparts_adjusted = round_up_pow2(nparts_adjusted);
- nparts_adjusted = nparts_adjusted*nparts_adjusted*nparts_adjusted;
-
- /* Iterate over the partition according to their Morton order and voxelize
- * the primitives that intersect it */
- omp_set_num_threads((int)atrstm->nthreads);
- #pragma omp parallel for schedule(static, 1/*chunk size*/)
- for(ipart = 0; ipart < nparts_adjusted; ++ipart) {
- struct darray_size_t* prims = NULL;
- struct pool* pool = NULL;
- double part_low[3];
- double part_upp[3];
- uint32_t part_ids[3];
- const int ithread = omp_get_thread_num();
- struct part* part = NULL;
- res_T res_local = RES_OK;
-
- /* Handle error */
- if(ATOMIC_GET(&res) != RES_OK) continue;
-
- /* Fetch the per thread list of primitives and partition pool */
- prims = darray_prims_list_data_get(&prims_list)+ithread;
- darray_size_t_clear(prims);
- pool = pools + ithread;
-
- morton_xyz_decode_u21(ipart, part_ids);
-
- /* Check that the partition is not out of bound due to Morton indexing */
- if(part_ids[0] >= nparts[0]
- || part_ids[1] >= nparts[1]
- || part_ids[2] >= nparts[2])
- continue;
-
- /* Compute the partition AABB */
- part_low[0] = (double)(part_ids[0] * part_def) * vxsz[0] + vol_low[0];
- part_low[1] = (double)(part_ids[1] * part_def) * vxsz[1] + vol_low[1];
- part_low[2] = (double)(part_ids[2] * part_def) * vxsz[2] + vol_low[2];
- part_upp[0] = part_low[0] + (double)part_def * vxsz[0];
- part_upp[1] = part_low[1] + (double)part_def * vxsz[1];
- part_upp[2] = part_low[2] + (double)part_def * vxsz[2];
-
- /* Find the list of primitives that intersects the partition */
- res_local = suvm_volume_intersect_aabb(atrstm->volume, part_low, part_upp,
- register_primitive, prims);
- if(res_local != RES_OK) {
- ATOMIC_SET(&res, res_local);
- continue;
- }
-
- part = pool_new_partition(pool, ipart);
- if(!part) { /* An error occurs */
- ATOMIC_SET(&res, res_local);
- continue;
- }
- res = voxelize_partition
- (atrstm, prims, part_low, part_upp, vxsz, refract_id, part);
- if(res != RES_OK) {
- pool_free_partition(pool, part);
- ATOMIC_SET(&res, res_local);
- continue;
- }
- pool_commit_partition(pool, part);
- }
-
-exit:
- darray_prims_list_release(&prims_list);
- return (res_T)res;
-error:
- goto exit;
-}
-
-static void
-voxel_get(const size_t xyz[3], const uint64_t mcode, void* dst, void* context)
-{
- const struct build_octree_context* ctx = context;
- struct pool* pool = NULL;
- struct part* part = NULL;
- float* vox = NULL;
- uint64_t ivox, ipart, ipool;
- ASSERT(xyz && dst && ctx);
- (void)xyz;
-
- /* Retrieve the partition morton id and the per partition voxel morton id
- * from the overall voxel morton code. As voxels, partitions are sorted in
- * morton order and thus the partition morton ID is encoded in the MSB of the
- * morton code while the voxels morton ID is stored in its LSB. */
- ipart = (mcode >> (LOG2_PARTITION_DEFINITION*3));
- ivox = (mcode & (BIT_U64(LOG2_PARTITION_DEFINITION*3)-1));
-
- /* Compute the pool index containing the partition. Partitions are
- * alternatively stored into the per thread pool. Consequentlu the i^th
- * partition is stored in the (i % #thread)^th pool. */
- ipool = ipart % ctx->atrstm->nthreads;
-
- /* Fetch the pool storing the partition */
- pool = ctx->pools + ipool;
-
- /* Fetch the partition storing the voxel */
- part = pool_fetch_partition(pool, ipart);
-
- if(!part) { /* An error occurs */
- memset(dst, 0, NFLOATS_PER_VOXEL * sizeof(float));
- } else {
- vox = part_get_voxel(part, ivox); /* Fetch the voxel data */
-
- /* Copy voxel data */
- memcpy(dst, vox, NFLOATS_PER_VOXEL * sizeof(float));
-
- /* Free the partition if the fetch voxel was its last one */
- if(ivox == PARTITION_NVOXELS - 1) {
- pool_free_partition(pool, part);
- }
- }
-}
-
-static void
-voxels_merge_component
- (void* dst,
- const enum atrstm_component cpnt,
- const void* voxels[],
- const size_t nvoxs)
-{
- float ka[2] = {FLT_MAX,-FLT_MAX};
- float ks[2] = {FLT_MAX,-FLT_MAX};
- float kext[2] = {FLT_MAX,-FLT_MAX};
- size_t ivox;
- ASSERT(dst && voxels && nvoxs);
-
- FOR_EACH(ivox, 0, nvoxs) {
- const float* vox = voxels[ivox];
- ka[0] = MMIN(ka[0], vox_get(vox, cpnt, ATRSTM_Ka, ATRSTM_SVX_MIN));
- ka[1] = MMAX(ka[1], vox_get(vox, cpnt, ATRSTM_Ka, ATRSTM_SVX_MAX));
- ks[0] = MMIN(ks[0], vox_get(vox, cpnt, ATRSTM_Ks, ATRSTM_SVX_MIN));
- ks[1] = MMAX(ks[1], vox_get(vox, cpnt, ATRSTM_Ks, ATRSTM_SVX_MAX));
- kext[0] = MMIN(kext[0], vox_get(vox, cpnt, ATRSTM_Kext, ATRSTM_SVX_MIN));
- kext[1] = MMAX(kext[1], vox_get(vox, cpnt, ATRSTM_Kext, ATRSTM_SVX_MAX));
- }
- vox_set(dst, cpnt, ATRSTM_Ka, ATRSTM_SVX_MIN, ka[0]);
- vox_set(dst, cpnt, ATRSTM_Ka, ATRSTM_SVX_MAX, ka[1]);
- vox_set(dst, cpnt, ATRSTM_Ks, ATRSTM_SVX_MIN, ks[0]);
- vox_set(dst, cpnt, ATRSTM_Ks, ATRSTM_SVX_MAX, ks[1]);
- vox_set(dst, cpnt, ATRSTM_Kext, ATRSTM_SVX_MIN, kext[0]);
- vox_set(dst, cpnt, ATRSTM_Kext, ATRSTM_SVX_MAX, kext[1]);
-}
-
-static void
-voxels_merge
- (void* dst,
- const void* voxels[],
- const size_t nvoxs,
- void* ctx)
-{
- (void)ctx;
- voxels_merge_component(dst, ATRSTM_CPNT_GAS, voxels, nvoxs);
- voxels_merge_component(dst, ATRSTM_CPNT_SOOT, voxels, nvoxs);
-}
-
-static INLINE void
-voxels_component_compute_kext_range
- (const enum atrstm_component cpnt,
- const struct svx_voxel voxels[],
- const size_t nvoxs,
- float kext[2])
-{
- size_t ivox;
- ASSERT(voxels && nvoxs && kext);
-
- kext[0] = FLT_MAX;
- kext[1] =-FLT_MAX;
-
- /* Compute the Kext range of the submitted voxels */
- FOR_EACH(ivox, 0, nvoxs) {
- const float* vox = voxels[ivox].data;
- kext[0] = MMIN(kext[0], vox_get(vox, cpnt, ATRSTM_Kext, ATRSTM_SVX_MIN));
- kext[1] = MMAX(kext[1], vox_get(vox, cpnt, ATRSTM_Kext, ATRSTM_SVX_MAX));
- }
-}
-
-static int
-voxels_challenge_merge
- (const struct svx_voxel voxels[],
- const size_t nvoxs,
- void* context)
-{
- const struct build_octree_context* ctx = context;
- double lower[3] = { DBL_MAX, DBL_MAX, DBL_MAX};
- double upper[3] = {-DBL_MAX,-DBL_MAX,-DBL_MAX};
- double sz_max;
- size_t ivox;
- float kext_gas[2] = {0, 0};
- float kext_soot[2] = {0, 0};
- int merge_gas = 0;
- int merge_soot = 0;
- ASSERT(voxels && nvoxs && context);
-
- /* Compute the AABB of the group of voxels and its maximum size */
- FOR_EACH(ivox, 0, nvoxs) {
- lower[0] = MMIN(voxels[ivox].lower[0], lower[0]);
- lower[1] = MMIN(voxels[ivox].lower[1], lower[1]);
- lower[2] = MMIN(voxels[ivox].lower[2], lower[2]);
- upper[0] = MMAX(voxels[ivox].upper[0], upper[0]);
- upper[1] = MMAX(voxels[ivox].upper[1], upper[1]);
- upper[2] = MMAX(voxels[ivox].upper[2], upper[2]);
- }
- sz_max = upper[0] - lower[0];
- sz_max = MMAX(upper[1] - lower[1], sz_max);
- sz_max = MMAX(upper[2] - lower[2], sz_max);
-
- /* Compute the Kext range for each component */
- voxels_component_compute_kext_range(ATRSTM_CPNT_GAS, voxels, nvoxs, kext_gas);
- voxels_component_compute_kext_range(ATRSTM_CPNT_SOOT, voxels, nvoxs, kext_soot);
-
- if(kext_gas[0] > kext_gas[1]) { /* Empty voxels */
- /* Always merge empty voxels */
- merge_gas = 1;
- merge_soot = 1;
- } else {
- /* Check if the voxels are candidates to the merge process by evaluating its
- * optical thickness regarding the maximum size of their AABB and a user
- * defined threshold */
- ASSERT(kext_gas[0] <= kext_gas[1]);
- ASSERT(kext_soot[0] <= kext_soot[1]);
- merge_gas = (kext_gas[1] - kext_gas[0])*sz_max <= ctx->tau_threshold;
- merge_soot = (kext_soot[1] - kext_soot[0])*sz_max <= ctx->tau_threshold;
- }
- return merge_gas && merge_soot;
-}
-
-static res_T
-build_octree
- (struct atrstm* atrstm,
- const struct build_octrees_args* args,
- struct pool* pools)
-{
- struct svx_voxel_desc vox_desc = SVX_VOXEL_DESC_NULL;
- struct build_octree_context ctx = BUILD_OCTREE_CONTEXT_NULL;
- double lower[3];
- double upper[3];
- size_t definition[3];
- res_T res = RES_OK;
- ASSERT(atrstm && args && pools);
-
- /* Setup the build octree context */
- ctx.atrstm = atrstm;
- ctx.pools = pools;
- ctx.tau_threshold = args->optical_thickness_threshold;
-
- /* Setup the voxel descriptor. TODO in shortwave, the NFLOATS_PER_VOXEL
- * could be divided by 2 since there is no gas to handle */
- vox_desc.get = voxel_get;
- vox_desc.merge = voxels_merge;
- vox_desc.challenge_merge = voxels_challenge_merge;
- vox_desc.context = &ctx;
- vox_desc.size = NFLOATS_PER_VOXEL * sizeof(float);
-
- SUVM(volume_get_aabb(atrstm->volume, lower, upper));
-
- /* Create the octree */
- definition[0] = (size_t)args->grid_max_definition[0];
- definition[1] = (size_t)args->grid_max_definition[1];
- definition[2] = (size_t)args->grid_max_definition[2];
- res = svx_octree_create
- (atrstm->svx, lower, upper, definition, &vox_desc, &atrstm->octree);
- if(res != RES_OK) goto error;
-
-exit:
- return res;
-error:
- goto exit;
-}
-
-/*******************************************************************************
- * Local functions
- ******************************************************************************/
-res_T
-build_octrees(struct atrstm* atrstm, const struct build_octrees_args* args)
-{
- char buf[128];
- struct atrri_refractive_index refract_id = ATRRI_REFRACTIVE_INDEX_NULL;
- struct time t0, t1;
- struct pool* pools = NULL;
- size_t i;
- double wlen;
- ATOMIC res = RES_OK;
- ASSERT(atrstm && check_build_octrees_args(args));
-
- /* Currently only shortwave is handled and in this case, radiative
- * transfert is monochromatic */
- ASSERT(atrstm->spectral_type == ATRSTM_SPECTRAL_SW);
- ASSERT(atrstm->wlen_range[0] == atrstm->wlen_range[1]);
- wlen = atrstm->wlen_range[0];
-
- /* Fetch the submitted wavelength */
- res = atrri_fetch_refractive_index(atrstm->atrri, wlen, &refract_id);
- if(res != RES_OK) goto error;
-
- /* Allocate the per thread pool of partition */
- pools = MEM_CALLOC(atrstm->allocator, atrstm->nthreads, sizeof(*pools));
- if(!pools) { res = RES_MEM_ERR; goto error; }
- FOR_EACH(i, 0, atrstm->nthreads) {
- res = pool_init(atrstm->allocator, pools+i);
- if(res != RES_OK) {
- log_err(atrstm,
- "Error initializing the pool of voxel partitions -- %s.\n",
- res_to_cstr((res_T)res));
- goto error;
- }
- }
-
- log_info(atrstm,
- "Evaluating and partitionning the field of optical properties.\n");
- time_current(&t0);
-
-#if 0
- {
- FILE* fp = fopen("optprops.txt", "w");
- CHK(fp);
- CHK(dump_optprops(atrstm, &refract_id, fp) == RES_OK);
- fclose(fp);
- }
-#endif
-
- omp_set_nested(1); /* Enable nested threads for voxelize_volumetric_mesh */
- #pragma omp parallel sections num_threads(2)
- {
- #pragma omp section
- {
- const res_T res_local = voxelize_volumetric_mesh
- (atrstm, args, &refract_id, pools);
- if(res_local != RES_OK) {
- size_t ipool;
- log_err(atrstm, "Error voxelizing the volumetric mesh -- %s\n",
- res_to_cstr(res_local));
- FOR_EACH(ipool, 0, atrstm->nthreads) {
- pool_invalidate(pools+ipool);
- }
- ATOMIC_SET(&res, res_local);
- }
- }
-
- #pragma omp section
- {
- const res_T res_local = build_octree(atrstm, args, pools);
- if(res_local != RES_OK) {
- size_t ipool;
- log_err(atrstm, "Error building the octree -- %s\n",
- res_to_cstr(res_local));
- FOR_EACH(ipool, 0, atrstm->nthreads) {
- pool_invalidate(pools+ipool);
- }
- ATOMIC_SET(&res, res_local);
- }
- }
- }
- if(res != RES_OK) goto error;
-
- time_sub(&t0, time_current(&t1), &t0);
- time_dump(&t0, TIME_ALL, NULL, buf, sizeof(buf));
- log_info(atrstm, "Setup the partitionning data structures in %s\n", buf);
-
- i = MEM_ALLOCATED_SIZE(&atrstm->svx_allocator);
- dump_memory_size(i, NULL, buf, sizeof(buf));
- log_info(atrstm, "Star-VoXel memory usage: %s\n", buf);
-
-exit:
- if(pools) {
- FOR_EACH(i, 0, atrstm->nthreads) pool_release(pools+i);
- MEM_RM(atrstm->allocator, pools);
- }
- return (res_T)res;
-error:
- goto exit;
-}
-
-void
-octrees_clean(struct atrstm* atrstm)
-{
- ASSERT(atrstm);
- if(atrstm->octree) SVX(tree_ref_put(atrstm->octree));
-}
diff --git a/src/atrstm_build_octrees.h b/src/atrstm_build_octrees.h
@@ -1,51 +0,0 @@
-/* Copyright (C) 2020 CNRS
- *
- * This program is free software: you can redistribute it and/or modify
- * it under the terms of the GNU General Public License as published by
- * the Free Software Foundation, either version 3 of the License, or
- * (at your option) any later version.
- *
- * This program is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- * GNU General Public License for more details.
- *
- * You should have received a copy of the GNU General Public License
- * along with this program. If not, see <http://www.gnu.org/licenses/>. */
-
-#ifndef ATRSTM_BUILD_OCTREES_H
-#define ATRSTM_BUILD_OCTREES_H
-
-#include <rsys/rsys.h>
-
-/* Forward declaration of external data type */
-struct sth;
-
-struct cache {
- FILE* stream; /* Stream where data are stored */
- const char* name; /* Name of the stream where data are stored */
- int force_update; /* Force overwrite of the cached data */
-};
-#define CACHE_NULL__ {NULL, NULL, 0}
-static const struct cache CACHE_NULL = CACHE_NULL__;
-
-struct build_octrees_args {
- unsigned grid_max_definition[3]; /* Max definition of the Octrees */
- double optical_thickness_threshold; /* Threshold used to build the octrees */
- struct cache cache; /* Cache of the acceleration data structures */
-};
-#define BUILD_OCTREES_ARGS_NULL__ {{0,0,0}, 0, CACHE_NULL__}
-static const struct build_octrees_args BUILD_OCTREES_ARGS_NULL =
- BUILD_OCTREES_ARGS_NULL__;
-
-extern LOCAL_SYM res_T
-build_octrees
- (struct atrstm* atstmm,
- const struct build_octrees_args* args);
-
-extern LOCAL_SYM void
-octrees_clean
- (struct atrstm* atrstm);
-
-#endif /* ATRSTM_BUILD_OCTREES_H */
-
diff --git a/src/atrstm_c.h b/src/atrstm_c.h
@@ -49,7 +49,9 @@ struct atrstm {
enum atrstm_spectral_type spectral_type;
double wlen_range[2]; /* Spectral range in nm */
-
+ unsigned grid_max_definition[3];
+ double optical_thickness;
+
struct cache* cache; /* Cache of the SVX data structures */
struct mem_allocator* allocator;
diff --git a/src/atrstm_cache.c b/src/atrstm_cache.c
@@ -36,6 +36,7 @@
struct cache {
FILE* stream;
struct atrstm* atrstm;
+ int empty;
ref_T ref;
};
@@ -162,6 +163,8 @@ write_cache_header(struct cache* cache, const char* filename)
WRITE(&cache->atrstm->fractal_dimension, 1);
WRITE(&cache->atrstm->spectral_type, 1);
WRITE(cache->atrstm->wlen_range, 2);
+ WRITE(cache->atrstm->grid_max_definition, 3);
+ WRITE(&cache->atrstm->optical_thickness, 1);
#undef WRITE
res = hash_compute(cache->atrstm, &hash);
@@ -183,7 +186,9 @@ read_cache_header(struct cache* cache, const char* filename)
double gyration_radius_prefactor = 0;
double fractal_dimension = 0;
enum atrstm_spectral_type spectral_type = ATRSTM_SPECTRAL_TYPES_COUNT__;
- double wlen_range[2] = {0, 0};
+ double wlen_range[2] = {0,0};
+ unsigned grid_max_definition[3] = {0,0,0};
+ double optical_thickness = 0;
int cache_version = 0;
res_T res = RES_OK;
ASSERT(cache);
@@ -241,6 +246,18 @@ read_cache_header(struct cache* cache, const char* filename)
CHK_VAR(wlen_range[1], cache->atrstm->wlen_range[1],
"spectral upper bound", "%g");
+ READ(grid_max_definition, 3);
+ CHK_VAR(grid_max_definition[0], cache->atrstm->grid_max_definition[0],
+ "grid X definition", "%u");
+ CHK_VAR(grid_max_definition[1], cache->atrstm->grid_max_definition[1],
+ "grid Y definition", "%u");
+ CHK_VAR(grid_max_definition[2], cache->atrstm->grid_max_definition[2],
+ "grid Z definition", "%u");
+
+ READ(&optical_thickness, 1);
+ CHK_VAR(optical_thickness, cache->atrstm->optical_thickness,
+ "optical thickness", "%g");
+
#undef CHK_VAR
res = hash_read(&cached_hash, cache->stream);
@@ -335,9 +352,11 @@ cache_create
}
if(cache_stat.st_size == 0) { /* Empty cache */
+ cache->empty = 1;
res = write_cache_header(cache, filename);
if(res != RES_OK) goto error;
} else { /* Cache already exists */
+ cache->empty = 0;
res = read_cache_header(cache, filename);
if(res != RES_OK) goto error;
}
@@ -371,3 +390,10 @@ cache_get_stream(struct cache* cache)
ASSERT(cache);
return cache->stream;
}
+
+int
+cache_is_empty(const struct cache* cache)
+{
+ ASSERT(cache);
+ return cache->empty;
+}
diff --git a/src/atrstm_cache.h b/src/atrstm_cache.h
@@ -41,4 +41,8 @@ extern LOCAL_SYM FILE*
cache_get_stream
(struct cache* cache);
+extern LOCAL_SYM int
+cache_is_empty
+ (const struct cache* cache);
+
#endif /* ATRSTM_CACHE_H */
diff --git a/src/atrstm_setup_octrees.c b/src/atrstm_setup_octrees.c
@@ -0,0 +1,689 @@
+/* Copyright (C) 2020 CNRS
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. */
+
+#define _POSIX_C_SOURCE 200112L /* nextafterf */
+
+#include "atrstm_c.h"
+#include "atrstm_cache.h"
+#include "atrstm_log.h"
+#include "atrstm_optprops.h"
+#include "atrstm_partition.h"
+#include "atrstm_setup_octrees.h"
+#include "atrstm_svx.h"
+
+#include <astoria/atrri.h>
+
+#include <star/suvm.h>
+#include <star/svx.h>
+
+#include <rsys/clock_time.h>
+#include <rsys/cstr.h>
+#include <rsys/dynamic_array_size_t.h>
+#include <rsys/morton.h>
+
+#include <math.h>
+#include <omp.h>
+
+struct build_octree_context {
+ struct atrstm* atrstm;
+ struct pool* pools; /* Per thread pool of voxel partitions */
+
+ /* Optical thickness threshold criteria for the merge process */
+ double tau_threshold;
+};
+#define BUILD_OCTREE_CONTEXT_NULL__ { NULL, NULL, 0 }
+static const struct build_octree_context BUILD_OCTREE_CONTEXT_NULL =
+ BUILD_OCTREE_CONTEXT_NULL__;
+
+/* Define the darray of list of primitive ids */
+#define DARRAY_NAME prims_list
+#define DARRAY_DATA struct darray_size_t
+#define DARRAY_FUNCTOR_INIT darray_size_t_init
+#define DARRAY_FUNCTOR_RELEASE darray_size_t_release
+#define DARRAY_FUNCTOR_COPY darray_size_t_copy
+#define DARRAY_FUNCTOR_COPY_AND_RELEASE darray_size_t_copy_and_release
+#include <rsys/dynamic_array.h>
+
+/*******************************************************************************
+ * Helper functions
+ ******************************************************************************/
+static INLINE int
+check_octrees_parameters(const struct atrstm* atrstm)
+{
+ return atrstm
+ && atrstm->grid_max_definition[0]
+ && atrstm->grid_max_definition[1]
+ && atrstm->grid_max_definition[2]
+ && atrstm->optical_thickness >= 0;
+}
+
+static INLINE void
+register_primitive
+ (const struct suvm_primitive* prim,
+ const double low[3],
+ const double upp[3],
+ void* context)
+{
+ struct darray_size_t* prims = context;
+ ASSERT(prim && low && upp && context);
+ ASSERT(low[0] < upp[0]);
+ ASSERT(low[1] < upp[1]);
+ ASSERT(low[2] < upp[2]);
+ (void)low, (void)upp;
+ CHK(darray_size_t_push_back(prims, &prim->iprim) == RES_OK);
+}
+
+static void
+voxel_commit_optprops_range
+ (float* vox,
+ const enum atrstm_component cpnt,
+ const struct optprops* optprops_min,
+ const struct optprops* optprops_max)
+{
+ double ka[2], ks[2], kext[2];
+ float vox_ka[2], vox_ks[2], vox_kext[2];
+ ASSERT(vox);
+ ASSERT(optprops_min);
+ ASSERT(optprops_max);
+
+ ka[0] = optprops_min->ka;
+ ka[1] = optprops_max->ka;
+ ks[0] = optprops_min->ks;
+ ks[1] = optprops_max->ks;
+ kext[0] = optprops_min->kext;
+ kext[1] = optprops_max->kext;
+
+ /* Ensure that the single precision bounds include their double precision
+ * version. */
+ if(ka[0] != (float)ka[0]) ka[0] = nextafterf((float)ka[0],-FLT_MAX);
+ if(ka[1] != (float)ka[1]) ka[1] = nextafterf((float)ka[1], FLT_MAX);
+ if(ks[0] != (float)ks[0]) ks[0] = nextafterf((float)ks[0],-FLT_MAX);
+ if(ks[1] != (float)ks[1]) ks[1] = nextafterf((float)ks[1], FLT_MAX);
+ if(kext[0] != (float)kext[0]) kext[0] = nextafterf((float)kext[0],-FLT_MAX);
+ if(kext[1] != (float)kext[1]) kext[1] = nextafterf((float)kext[1], FLT_MAX);
+
+ /* Fetch the range of the voxel optical properties */
+ vox_ka[0] = vox_get(vox, cpnt, ATRSTM_Ka, ATRSTM_SVX_MIN);
+ vox_ka[1] = vox_get(vox, cpnt, ATRSTM_Ka, ATRSTM_SVX_MAX);
+ vox_ks[0] = vox_get(vox, cpnt, ATRSTM_Ks, ATRSTM_SVX_MIN);
+ vox_ks[1] = vox_get(vox, cpnt, ATRSTM_Ks, ATRSTM_SVX_MAX);
+ vox_kext[0] = vox_get(vox, cpnt, ATRSTM_Kext, ATRSTM_SVX_MIN);
+ vox_kext[1] = vox_get(vox, cpnt, ATRSTM_Kext, ATRSTM_SVX_MAX);
+
+ /* Update the range of the voxel optical properties */
+ vox_ka[0] = MMIN(vox_ka[0], (float)ka[0]);
+ vox_ka[1] = MMAX(vox_ka[1], (float)ka[1]);
+ vox_ks[0] = MMIN(vox_ks[0], (float)ks[0]);
+ vox_ks[1] = MMAX(vox_ks[1], (float)ks[1]);
+ vox_kext[0] = MMIN(vox_kext[0], (float)kext[0]);
+ vox_kext[1] = MMAX(vox_kext[1], (float)kext[1]);
+
+ /* Commit the upd to the voxel */
+ vox_set(vox, cpnt, ATRSTM_Ka, ATRSTM_SVX_MIN, vox_ka[0]);
+ vox_set(vox, cpnt, ATRSTM_Ka, ATRSTM_SVX_MAX, vox_ka[1]);
+ vox_set(vox, cpnt, ATRSTM_Ks, ATRSTM_SVX_MIN, vox_ks[0]);
+ vox_set(vox, cpnt, ATRSTM_Ks, ATRSTM_SVX_MAX, vox_ks[1]);
+ vox_set(vox, cpnt, ATRSTM_Kext, ATRSTM_SVX_MIN, vox_kext[0]);
+ vox_set(vox, cpnt, ATRSTM_Kext, ATRSTM_SVX_MAX, vox_kext[1]);
+}
+
+static res_T
+voxelize_partition
+ (struct atrstm* atrstm,
+ const struct darray_size_t* prims, /* Primitives intersecting the partition */
+ const double part_low[3], /* Lower bound of the partition */
+ const double part_upp[3], /* Upper bound of the partition */
+ const double vxsz[3], /* Size of a partition voxel */
+ const struct atrri_refractive_index* refract_id, /* Refractive index */
+ struct part* part)
+{
+ size_t iprim;
+ res_T res = RES_OK;
+ ASSERT(atrstm && prims && part_low && part_upp && vxsz && part && refract_id);
+ ASSERT(part_low[0] < part_upp[0]);
+ ASSERT(part_low[1] < part_upp[1]);
+ ASSERT(part_low[2] < part_upp[2]);
+ ASSERT(vxsz[0] > 0 && vxsz[1] > 0 && vxsz[2] > 0);
+
+ /* Clean the partition voxels */
+ part_clear_voxels(part);
+
+ FOR_EACH(iprim, 0, darray_size_t_size_get(prims)) {
+ struct suvm_primitive prim = SUVM_PRIMITIVE_NULL;
+ struct suvm_polyhedron poly;
+ double poly_low[3];
+ double poly_upp[3];
+ uint32_t ivxl_low[3];
+ uint32_t ivxl_upp[3];
+ uint32_t ivxl[3];
+ size_t prim_id;
+
+ prim_id = darray_size_t_cdata_get(prims)[iprim];
+
+ /* Retrieve the primitive data and setup its polyhedron */
+ SUVM(volume_get_primitive(atrstm->volume, prim_id, &prim));
+ SUVM(primitive_setup_polyhedron(&prim, &poly));
+ ASSERT(poly.lower[0] <= part_upp[0] && poly.upper[0] >= part_low[0]);
+ ASSERT(poly.lower[1] <= part_upp[1] && poly.upper[1] >= part_low[1]);
+ ASSERT(poly.lower[2] <= part_upp[2] && poly.upper[2] >= part_low[2]);
+
+ /* Clamp the poly AABB to the partition boundaries */
+ poly_low[0] = MMAX(poly.lower[0], part_low[0]);
+ poly_low[1] = MMAX(poly.lower[1], part_low[1]);
+ poly_low[2] = MMAX(poly.lower[2], part_low[2]);
+ poly_upp[0] = MMIN(poly.upper[0], part_upp[0]);
+ poly_upp[1] = MMIN(poly.upper[1], part_upp[1]);
+ poly_upp[2] = MMIN(poly.upper[2], part_upp[2]);
+
+ /* Transform the polyhedron AABB in partition voxel space */
+ ivxl_low[0] = (uint32_t)((poly_low[0] - part_low[0]) / vxsz[0]);
+ ivxl_low[1] = (uint32_t)((poly_low[1] - part_low[1]) / vxsz[1]);
+ ivxl_low[2] = (uint32_t)((poly_low[2] - part_low[2]) / vxsz[2]);
+ ivxl_upp[0] = (uint32_t)ceil((poly_upp[0] - part_low[0]) / vxsz[0]);
+ ivxl_upp[1] = (uint32_t)ceil((poly_upp[1] - part_low[1]) / vxsz[1]);
+ ivxl_upp[2] = (uint32_t)ceil((poly_upp[2] - part_low[2]) / vxsz[2]);
+ ASSERT(ivxl_upp[0] <= PARTITION_DEFINITION);
+ ASSERT(ivxl_upp[0] <= PARTITION_DEFINITION);
+ ASSERT(ivxl_upp[0] <= PARTITION_DEFINITION);
+
+ /* Iterate over the partition voxels intersected by the primitive AABB */
+ FOR_EACH(ivxl[2], ivxl_low[2], ivxl_upp[2]) {
+ float vxl_low[3]; /* Voxel lower bound */
+ float vxl_upp[3]; /* Voxel upper bound */
+ uint64_t mcode[3]; /* Morton code cache */
+
+ vxl_low[2] = (float)((double)ivxl[2] * vxsz[2] + part_low[2]);
+ vxl_upp[2] = vxl_low[2] + (float)vxsz[2];
+ mcode[2] = morton3D_encode_u21(ivxl[2]) << 0;
+
+ FOR_EACH(ivxl[1], ivxl_low[1], ivxl_upp[1]) {
+ vxl_low[1] = (float)((double)ivxl[1] * vxsz[1] + part_low[1]);
+ vxl_upp[1] = vxl_low[1] + (float)vxsz[1];
+ mcode[1] = (morton3D_encode_u21(ivxl[1]) << 1) | mcode[2];
+
+ FOR_EACH(ivxl[0], ivxl_low[0], ivxl_upp[0]) {
+ enum suvm_intersection_type intersect;
+ vxl_low[0] = (float)((double)ivxl[0] * vxsz[0] + part_low[0]);
+ vxl_upp[0] = vxl_low[0] + (float)vxsz[0];
+
+ /* NOTE mcode[0] store the morton index of the voxel */
+ mcode[0] = (morton3D_encode_u21(ivxl[0]) << 2) | mcode[1];
+
+ intersect = suvm_polyhedron_intersect_aabb(&poly, vxl_low, vxl_upp);
+ if(intersect != SUVM_INTERSECT_NONE) {
+ struct optprops optprops_min;
+ struct optprops optprops_max;
+ float* vox = part_get_voxel(part, mcode[2]);
+
+ /* Currently, gas is not handled */
+ optprops_min = OPTPROPS_NULL;
+ optprops_max = OPTPROPS_NULL;
+ voxel_commit_optprops_range
+ (vox, ATRSTM_CPNT_GAS, &optprops_min, &optprops_max);
+
+ res = primitive_compute_optprops_range
+ (atrstm, refract_id, &prim, &optprops_min, &optprops_max);
+ if(UNLIKELY(res != RES_OK)) goto error;
+ voxel_commit_optprops_range
+ (vox, ATRSTM_CPNT_SOOT, &optprops_min, &optprops_max);
+ }
+ }
+ }
+ }
+ }
+
+exit:
+ return res;
+error:
+ goto exit;
+}
+
+static res_T
+voxelize_volumetric_mesh
+ (struct atrstm* atrstm,
+ const struct atrri_refractive_index* refract_id,
+ struct pool* pools)
+{
+ struct darray_prims_list prims_list; /* Per thread list of primitives */
+ const size_t part_def = PARTITION_DEFINITION;
+ size_t nparts[3]; /* #partitions along the 3 axis */
+ size_t nparts_adjusted;
+ double vol_low[3];
+ double vol_upp[3];
+ double vxsz[3];
+ uint64_t ipart;
+ ATOMIC res = RES_OK;
+ ASSERT(atrstm && pools);
+
+ darray_prims_list_init(atrstm->allocator, &prims_list);
+
+ /* Allocate the per thread primitive lists */
+ res = darray_prims_list_resize(&prims_list, atrstm->nthreads);
+ if(res != RES_OK) goto error;
+
+ /* Fetch the volume AABB and compute the size of one voxel */
+ SUVM(volume_get_aabb(atrstm->volume, vol_low, vol_upp));
+ vxsz[0] = (vol_upp[0] - vol_low[0]) / (double)atrstm->grid_max_definition[0];
+ vxsz[1] = (vol_upp[1] - vol_low[1]) / (double)atrstm->grid_max_definition[1];
+ vxsz[2] = (vol_upp[2] - vol_low[2]) / (double)atrstm->grid_max_definition[2];
+
+ /* Compute the number of partitions */
+ nparts[0] = (atrstm->grid_max_definition[0] + (part_def-1)/*ceil*/) / part_def;
+ nparts[1] = (atrstm->grid_max_definition[1] + (part_def-1)/*ceil*/) / part_def;
+ nparts[2] = (atrstm->grid_max_definition[2] + (part_def-1)/*ceil*/) / part_def;
+
+ /* Compute the overall #partitions allowing Morton indexing */
+ nparts_adjusted = MMAX(nparts[0], MMAX(nparts[1], nparts[2]));
+ nparts_adjusted = round_up_pow2(nparts_adjusted);
+ nparts_adjusted = nparts_adjusted*nparts_adjusted*nparts_adjusted;
+
+ /* Iterate over the partition according to their Morton order and voxelize
+ * the primitives that intersect it */
+ omp_set_num_threads((int)atrstm->nthreads);
+ #pragma omp parallel for schedule(static, 1/*chunk size*/)
+ for(ipart = 0; ipart < nparts_adjusted; ++ipart) {
+ struct darray_size_t* prims = NULL;
+ struct pool* pool = NULL;
+ double part_low[3];
+ double part_upp[3];
+ uint32_t part_ids[3];
+ const int ithread = omp_get_thread_num();
+ struct part* part = NULL;
+ res_T res_local = RES_OK;
+
+ /* Handle error */
+ if(ATOMIC_GET(&res) != RES_OK) continue;
+
+ /* Fetch the per thread list of primitives and partition pool */
+ prims = darray_prims_list_data_get(&prims_list)+ithread;
+ darray_size_t_clear(prims);
+ pool = pools + ithread;
+
+ morton_xyz_decode_u21(ipart, part_ids);
+
+ /* Check that the partition is not out of bound due to Morton indexing */
+ if(part_ids[0] >= nparts[0]
+ || part_ids[1] >= nparts[1]
+ || part_ids[2] >= nparts[2])
+ continue;
+
+ /* Compute the partition AABB */
+ part_low[0] = (double)(part_ids[0] * part_def) * vxsz[0] + vol_low[0];
+ part_low[1] = (double)(part_ids[1] * part_def) * vxsz[1] + vol_low[1];
+ part_low[2] = (double)(part_ids[2] * part_def) * vxsz[2] + vol_low[2];
+ part_upp[0] = part_low[0] + (double)part_def * vxsz[0];
+ part_upp[1] = part_low[1] + (double)part_def * vxsz[1];
+ part_upp[2] = part_low[2] + (double)part_def * vxsz[2];
+
+ /* Find the list of primitives that intersects the partition */
+ res_local = suvm_volume_intersect_aabb(atrstm->volume, part_low, part_upp,
+ register_primitive, prims);
+ if(res_local != RES_OK) {
+ ATOMIC_SET(&res, res_local);
+ continue;
+ }
+
+ part = pool_new_partition(pool, ipart);
+ if(!part) { /* An error occurs */
+ ATOMIC_SET(&res, res_local);
+ continue;
+ }
+ res = voxelize_partition
+ (atrstm, prims, part_low, part_upp, vxsz, refract_id, part);
+ if(res != RES_OK) {
+ pool_free_partition(pool, part);
+ ATOMIC_SET(&res, res_local);
+ continue;
+ }
+ pool_commit_partition(pool, part);
+ }
+
+exit:
+ darray_prims_list_release(&prims_list);
+ return (res_T)res;
+error:
+ goto exit;
+}
+
+static void
+voxel_get(const size_t xyz[3], const uint64_t mcode, void* dst, void* context)
+{
+ const struct build_octree_context* ctx = context;
+ struct pool* pool = NULL;
+ struct part* part = NULL;
+ float* vox = NULL;
+ uint64_t ivox, ipart, ipool;
+ ASSERT(xyz && dst && ctx);
+ (void)xyz;
+
+ /* Retrieve the partition morton id and the per partition voxel morton id
+ * from the overall voxel morton code. As voxels, partitions are sorted in
+ * morton order and thus the partition morton ID is encoded in the MSB of the
+ * morton code while the voxels morton ID is stored in its LSB. */
+ ipart = (mcode >> (LOG2_PARTITION_DEFINITION*3));
+ ivox = (mcode & (BIT_U64(LOG2_PARTITION_DEFINITION*3)-1));
+
+ /* Compute the pool index containing the partition. Partitions are
+ * alternatively stored into the per thread pool. Consequentlu the i^th
+ * partition is stored in the (i % #thread)^th pool. */
+ ipool = ipart % ctx->atrstm->nthreads;
+
+ /* Fetch the pool storing the partition */
+ pool = ctx->pools + ipool;
+
+ /* Fetch the partition storing the voxel */
+ part = pool_fetch_partition(pool, ipart);
+
+ if(!part) { /* An error occurs */
+ memset(dst, 0, NFLOATS_PER_VOXEL * sizeof(float));
+ } else {
+ vox = part_get_voxel(part, ivox); /* Fetch the voxel data */
+
+ /* Copy voxel data */
+ memcpy(dst, vox, NFLOATS_PER_VOXEL * sizeof(float));
+
+ /* Free the partition if the fetch voxel was its last one */
+ if(ivox == PARTITION_NVOXELS - 1) {
+ pool_free_partition(pool, part);
+ }
+ }
+}
+
+static void
+voxels_merge_component
+ (void* dst,
+ const enum atrstm_component cpnt,
+ const void* voxels[],
+ const size_t nvoxs)
+{
+ float ka[2] = {FLT_MAX,-FLT_MAX};
+ float ks[2] = {FLT_MAX,-FLT_MAX};
+ float kext[2] = {FLT_MAX,-FLT_MAX};
+ size_t ivox;
+ ASSERT(dst && voxels && nvoxs);
+
+ FOR_EACH(ivox, 0, nvoxs) {
+ const float* vox = voxels[ivox];
+ ka[0] = MMIN(ka[0], vox_get(vox, cpnt, ATRSTM_Ka, ATRSTM_SVX_MIN));
+ ka[1] = MMAX(ka[1], vox_get(vox, cpnt, ATRSTM_Ka, ATRSTM_SVX_MAX));
+ ks[0] = MMIN(ks[0], vox_get(vox, cpnt, ATRSTM_Ks, ATRSTM_SVX_MIN));
+ ks[1] = MMAX(ks[1], vox_get(vox, cpnt, ATRSTM_Ks, ATRSTM_SVX_MAX));
+ kext[0] = MMIN(kext[0], vox_get(vox, cpnt, ATRSTM_Kext, ATRSTM_SVX_MIN));
+ kext[1] = MMAX(kext[1], vox_get(vox, cpnt, ATRSTM_Kext, ATRSTM_SVX_MAX));
+ }
+ vox_set(dst, cpnt, ATRSTM_Ka, ATRSTM_SVX_MIN, ka[0]);
+ vox_set(dst, cpnt, ATRSTM_Ka, ATRSTM_SVX_MAX, ka[1]);
+ vox_set(dst, cpnt, ATRSTM_Ks, ATRSTM_SVX_MIN, ks[0]);
+ vox_set(dst, cpnt, ATRSTM_Ks, ATRSTM_SVX_MAX, ks[1]);
+ vox_set(dst, cpnt, ATRSTM_Kext, ATRSTM_SVX_MIN, kext[0]);
+ vox_set(dst, cpnt, ATRSTM_Kext, ATRSTM_SVX_MAX, kext[1]);
+}
+
+static void
+voxels_merge
+ (void* dst,
+ const void* voxels[],
+ const size_t nvoxs,
+ void* ctx)
+{
+ (void)ctx;
+ voxels_merge_component(dst, ATRSTM_CPNT_GAS, voxels, nvoxs);
+ voxels_merge_component(dst, ATRSTM_CPNT_SOOT, voxels, nvoxs);
+}
+
+static INLINE void
+voxels_component_compute_kext_range
+ (const enum atrstm_component cpnt,
+ const struct svx_voxel voxels[],
+ const size_t nvoxs,
+ float kext[2])
+{
+ size_t ivox;
+ ASSERT(voxels && nvoxs && kext);
+
+ kext[0] = FLT_MAX;
+ kext[1] =-FLT_MAX;
+
+ /* Compute the Kext range of the submitted voxels */
+ FOR_EACH(ivox, 0, nvoxs) {
+ const float* vox = voxels[ivox].data;
+ kext[0] = MMIN(kext[0], vox_get(vox, cpnt, ATRSTM_Kext, ATRSTM_SVX_MIN));
+ kext[1] = MMAX(kext[1], vox_get(vox, cpnt, ATRSTM_Kext, ATRSTM_SVX_MAX));
+ }
+}
+
+static int
+voxels_challenge_merge
+ (const struct svx_voxel voxels[],
+ const size_t nvoxs,
+ void* context)
+{
+ const struct build_octree_context* ctx = context;
+ double lower[3] = { DBL_MAX, DBL_MAX, DBL_MAX};
+ double upper[3] = {-DBL_MAX,-DBL_MAX,-DBL_MAX};
+ double sz_max;
+ size_t ivox;
+ float kext_gas[2] = {0, 0};
+ float kext_soot[2] = {0, 0};
+ int merge_gas = 0;
+ int merge_soot = 0;
+ ASSERT(voxels && nvoxs && context);
+
+ /* Compute the AABB of the group of voxels and its maximum size */
+ FOR_EACH(ivox, 0, nvoxs) {
+ lower[0] = MMIN(voxels[ivox].lower[0], lower[0]);
+ lower[1] = MMIN(voxels[ivox].lower[1], lower[1]);
+ lower[2] = MMIN(voxels[ivox].lower[2], lower[2]);
+ upper[0] = MMAX(voxels[ivox].upper[0], upper[0]);
+ upper[1] = MMAX(voxels[ivox].upper[1], upper[1]);
+ upper[2] = MMAX(voxels[ivox].upper[2], upper[2]);
+ }
+ sz_max = upper[0] - lower[0];
+ sz_max = MMAX(upper[1] - lower[1], sz_max);
+ sz_max = MMAX(upper[2] - lower[2], sz_max);
+
+ /* Compute the Kext range for each component */
+ voxels_component_compute_kext_range(ATRSTM_CPNT_GAS, voxels, nvoxs, kext_gas);
+ voxels_component_compute_kext_range(ATRSTM_CPNT_SOOT, voxels, nvoxs, kext_soot);
+
+ if(kext_gas[0] > kext_gas[1]) { /* Empty voxels */
+ /* Always merge empty voxels */
+ merge_gas = 1;
+ merge_soot = 1;
+ } else {
+ /* Check if the voxels are candidates to the merge process by evaluating its
+ * optical thickness regarding the maximum size of their AABB and a user
+ * defined threshold */
+ ASSERT(kext_gas[0] <= kext_gas[1]);
+ ASSERT(kext_soot[0] <= kext_soot[1]);
+ merge_gas = (kext_gas[1] - kext_gas[0])*sz_max <= ctx->tau_threshold;
+ merge_soot = (kext_soot[1] - kext_soot[0])*sz_max <= ctx->tau_threshold;
+ }
+ return merge_gas && merge_soot;
+}
+
+static res_T
+build_octree
+ (struct atrstm* atrstm,
+ struct pool* pools)
+{
+ struct svx_voxel_desc vox_desc = SVX_VOXEL_DESC_NULL;
+ struct build_octree_context ctx = BUILD_OCTREE_CONTEXT_NULL;
+ double lower[3];
+ double upper[3];
+ size_t definition[3];
+ res_T res = RES_OK;
+ ASSERT(atrstm && pools);
+
+ /* Setup the build octree context */
+ ctx.atrstm = atrstm;
+ ctx.pools = pools;
+ ctx.tau_threshold = atrstm->optical_thickness;
+
+ /* Setup the voxel descriptor. TODO in shortwave, the NFLOATS_PER_VOXEL
+ * could be divided by 2 since there is no gas to handle */
+ vox_desc.get = voxel_get;
+ vox_desc.merge = voxels_merge;
+ vox_desc.challenge_merge = voxels_challenge_merge;
+ vox_desc.context = &ctx;
+ vox_desc.size = NFLOATS_PER_VOXEL * sizeof(float);
+
+ SUVM(volume_get_aabb(atrstm->volume, lower, upper));
+
+ /* Create the octree */
+ definition[0] = (size_t)atrstm->grid_max_definition[0];
+ definition[1] = (size_t)atrstm->grid_max_definition[1];
+ definition[2] = (size_t)atrstm->grid_max_definition[2];
+ res = svx_octree_create
+ (atrstm->svx, lower, upper, definition, &vox_desc, &atrstm->octree);
+ if(res != RES_OK) goto error;
+
+exit:
+ return res;
+error:
+ goto exit;
+}
+
+static res_T
+build_octrees(struct atrstm* atrstm)
+{
+ char buf[128];
+ struct atrri_refractive_index refract_id = ATRRI_REFRACTIVE_INDEX_NULL;
+ struct time t0, t1;
+ struct pool* pools = NULL;
+ size_t i;
+ double wlen;
+ ATOMIC res = RES_OK;
+ ASSERT(atrstm && check_octrees_parameters(atrstm));
+
+ /* Currently only shortwave is handled and in this case, radiative
+ * transfert is monochromatic */
+ ASSERT(atrstm->spectral_type == ATRSTM_SPECTRAL_SW);
+ ASSERT(atrstm->wlen_range[0] == atrstm->wlen_range[1]);
+ wlen = atrstm->wlen_range[0];
+
+ /* Fetch the submitted wavelength */
+ res = atrri_fetch_refractive_index(atrstm->atrri, wlen, &refract_id);
+ if(res != RES_OK) goto error;
+
+ /* Allocate the per thread pool of partition */
+ pools = MEM_CALLOC(atrstm->allocator, atrstm->nthreads, sizeof(*pools));
+ if(!pools) { res = RES_MEM_ERR; goto error; }
+ FOR_EACH(i, 0, atrstm->nthreads) {
+ res = pool_init(atrstm->allocator, pools+i);
+ if(res != RES_OK) {
+ log_err(atrstm,
+ "Error initializing the pool of voxel partitions -- %s.\n",
+ res_to_cstr((res_T)res));
+ goto error;
+ }
+ }
+
+ log_info(atrstm,
+ "Evaluating and partitionning the field of optical properties.\n");
+ time_current(&t0);
+
+ omp_set_nested(1); /* Enable nested threads for voxelize_volumetric_mesh */
+ #pragma omp parallel sections num_threads(2)
+ {
+ #pragma omp section
+ {
+ const res_T res_local = voxelize_volumetric_mesh
+ (atrstm, &refract_id, pools);
+ if(res_local != RES_OK) {
+ size_t ipool;
+ log_err(atrstm, "Error voxelizing the volumetric mesh -- %s\n",
+ res_to_cstr(res_local));
+ FOR_EACH(ipool, 0, atrstm->nthreads) {
+ pool_invalidate(pools+ipool);
+ }
+ ATOMIC_SET(&res, res_local);
+ }
+ }
+
+ #pragma omp section
+ {
+ const res_T res_local = build_octree(atrstm, pools);
+ if(res_local != RES_OK) {
+ size_t ipool;
+ log_err(atrstm, "Error building the octree -- %s\n",
+ res_to_cstr(res_local));
+ FOR_EACH(ipool, 0, atrstm->nthreads) {
+ pool_invalidate(pools+ipool);
+ }
+ ATOMIC_SET(&res, res_local);
+ }
+ }
+ }
+ if(res != RES_OK) goto error;
+
+ time_sub(&t0, time_current(&t1), &t0);
+ time_dump(&t0, TIME_ALL, NULL, buf, sizeof(buf));
+ log_info(atrstm, "Setup the partitionning data structures in %s\n", buf);
+
+ if(atrstm->cache) {
+ ASSERT(cache_is_empty(atrstm->cache));
+ res = svx_tree_write(atrstm->octree, cache_get_stream(atrstm->cache));
+ if(res != RES_OK) {
+ log_err(atrstm, "Could not write the octrees to the cache -- %s.\n",
+ res_to_cstr((res_T)res));
+ goto error;
+ }
+ }
+
+ i = MEM_ALLOCATED_SIZE(&atrstm->svx_allocator);
+ dump_memory_size(i, NULL, buf, sizeof(buf));
+ log_info(atrstm, "Star-VoXel memory usage: %s\n", buf);
+
+exit:
+ if(pools) {
+ FOR_EACH(i, 0, atrstm->nthreads) pool_release(pools+i);
+ MEM_RM(atrstm->allocator, pools);
+ }
+ return (res_T)res;
+error:
+ goto exit;
+}
+
+static res_T
+load_octrees(struct atrstm* atrstm)
+{
+ FILE* stream;
+ ASSERT(atrstm && atrstm->cache && !cache_is_empty(atrstm->cache));
+
+ stream = cache_get_stream(atrstm->cache);
+ return svx_tree_create_from_stream(atrstm->svx, stream, &atrstm->octree);
+}
+
+/*******************************************************************************
+ * Local functions
+ ******************************************************************************/
+res_T
+setup_octrees(struct atrstm* atrstm)
+{
+ ASSERT(atrstm);
+ if(atrstm->cache && !cache_is_empty(atrstm->cache)) {
+ return load_octrees(atrstm);;
+ } else {
+ return build_octrees(atrstm);
+ }
+}
+
+void
+octrees_clean(struct atrstm* atrstm)
+{
+ ASSERT(atrstm);
+ if(atrstm->octree) SVX(tree_ref_put(atrstm->octree));
+}
diff --git a/src/atrstm_setup_octrees.h b/src/atrstm_setup_octrees.h
@@ -0,0 +1,33 @@
+/* Copyright (C) 2020 CNRS
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. */
+
+#ifndef ATRSTM_SETUP_OCTREES_H
+#define ATRSTM_SETUP_OCTREES_H
+
+#include <rsys/rsys.h>
+
+/* Forward declaration of external data type */
+struct sth;
+
+extern LOCAL_SYM res_T
+setup_octrees
+ (struct atrstm* atrstm);
+
+extern LOCAL_SYM void
+octrees_clean
+ (struct atrstm* atrstm);
+
+#endif /* ATRSTM_SETUP_OCTREES_H */
+
