star-vx

Structuring voxels for ray-tracing
git clone git://git.meso-star.fr/star-vx.git
Log | Files | Refs | README | LICENSE

commit 6ab54cf71641413bb95e1f987021d0ee40d7aedd
parent bdcf2435076972f31f252793fc870ce6f023565e
Author: Vincent Forest <vincent.forest@meso-star.com>
Date:   Thu,  1 Mar 2018 16:09:04 +0100

Rename the scene in octree

Diffstat:
Mcmake/CMakeLists.txt | 10+++++-----
Msrc/htvox.h | 58+++++++++++++++++++++++++++++++++-------------------------
Asrc/htvox_octree.c | 783+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Asrc/htvox_octree.h | 40++++++++++++++++++++++++++++++++++++++++
Msrc/htvox_octree_buffer.h | 4++--
Dsrc/htvox_scene.c | 765-------------------------------------------------------------------------------
Dsrc/htvox_scene.h | 37-------------------------------------
Asrc/test_htvox_octree.c | 381+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Dsrc/test_htvox_scene.c | 381-------------------------------------------------------------------------------
9 files changed, 1244 insertions(+), 1215 deletions(-)

diff --git a/cmake/CMakeLists.txt b/cmake/CMakeLists.txt @@ -42,12 +42,12 @@ set(VERSION ${VERSION_MAJOR}.${VERSION_MINOR}.${VERSION_PATCH}) set(HTVOX_FILES_SRC htvox_device.c - htvox_octree_buffer.c - htvox_scene.c) + htvox_octree.c + htvox_octree_buffer.c) set(SSOL_FILES_INC htvox_device.h - htvox_octree_buffer.h - htvox_scene.h) + htvox_octree.h + htvox_octree_buffer.h) set(SSOL_FILES_INC_API htvox.h) @@ -86,7 +86,7 @@ if(NOT NO_TEST) endfunction() new_test(test_htvox_device) - new_test(test_htvox_scene) + new_test(test_htvox_octree) endif() ################################################################################ diff --git a/src/htvox.h b/src/htvox.h @@ -82,6 +82,20 @@ struct htvox_voxel_desc { static const struct htvox_voxel_desc HTVOX_VOXEL_DESC_NULL = HTVOX_VOXEL_DESC_NULL__; +struct htvox_octree_desc { + /* Submitted Axis Aligned Bounding Box */ + double lower[3], upper[3]; + + size_t nleaves; /* #leaves */ + size_t nvoxels; /* #voxels, i.e. #leaves + #parents */ + size_t depth; /* Maximum depth of the octree */ +}; + +#define HTVOX_OCTREE_DESC_NULL__ \ + {{DBL_MAX, DBL_MAX, DBL_MAX}, {-DBL_MAX,-DBL_MAX,-DBL_MAX}, 0, 0, 0} +static const struct htvox_octree_desc HTVOX_OCTREE_DESC_NULL = + HTVOX_OCTREE_DESC_NULL__; + struct htvox_hit { double distance; /* Distance from the ray origin to the impacted voxel */ struct htvox_voxel voxel; /* Intersected voxel */ @@ -98,7 +112,7 @@ struct mem_allocator; /* Forward declaration of opaque data types */ struct htvox_device; -struct htvox_scene; +struct htvox_octree; BEGIN_DECLS @@ -121,39 +135,33 @@ htvox_device_ref_put (struct htvox_device* htvox); /******************************************************************************* - * Scene + * Octree ******************************************************************************/ HTVOX_API res_T -htvox_scene_create +htvox_octree_create (struct htvox_device* dev, - const double lower[3], /* Lower bound of the scene */ - const double upper[3], /* Upper bound of the scene */ + const double lower[3], /* Lower bound of the octree */ + const double upper[3], /* Upper bound of the octree */ const size_t nvoxels[3], /* # voxels along the 3 axis */ const struct htvox_voxel_desc* desc, /* Descriptor of a voxel */ - struct htvox_scene** scn); - -HTVOX_API res_T -htvox_scene_ref_get - (struct htvox_scene* scn); + struct htvox_octree** octree); HTVOX_API res_T -htvox_scene_ref_put - (struct htvox_scene* scn); +htvox_octree_ref_get + (struct htvox_octree* octree); HTVOX_API res_T -htvox_scene_get_leaves_count - (const struct htvox_scene* scn, - size_t* nleaves); +htvox_octree_ref_put + (struct htvox_octree* octree); HTVOX_API res_T -htvox_scene_get_aabb - (const struct htvox_scene* scn, - double lower[3], - double upper[3]); +htvox_octree_get_desc + (const struct htvox_octree* octree, + struct htvox_octree_desc* desc); HTVOX_API res_T -htvox_scene_for_each_leaf - (struct htvox_scene* scn, +htvox_octree_for_each_leaf + (struct htvox_octree* octree, void (*functor) (const void* val, /* Value of the voxel */ const size_t ileaf, /* Identifier of the leaf in [0, #leafs[ */ @@ -163,16 +171,16 @@ htvox_scene_for_each_leaf void* context); /* Client data sent as the last argument of the callback */ HTVOX_API res_T -htvox_scene_trace_ray - (struct htvox_scene* scn, +htvox_octree_trace_ray + (struct htvox_octree* octree, const double ray_origin[3], const double ray_direction[3], const double ray_range[2], struct htvox_hit* hit); HTVOX_API res_T -htvox_scene_at - (struct htvox_scene* scn, +htvox_octree_at + (struct htvox_octree* octree, const double position[3], struct htvox_voxel* voxel); diff --git a/src/htvox_octree.c b/src/htvox_octree.c @@ -0,0 +1,783 @@ +/* Copyright (C) CNRS 2018 + * + * 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/>. */ + +#include "htvox.h" +#include "htvox_c.h" +#include "htvox_device.h" +#include "htvox_octree.h" + +#include <rsys/double3.h> +#include <rsys/mem_allocator.h> + +#define OCTREE_DEPTH_MAX 16 /* Maximum depth of an octree */ + +struct voxel { + uint64_t mcode; /* Morton code of the voxel */ + void* data; /* Data of the voxel */ +}; +static const struct voxel VOXEL_NULL = {0, NULL}; + +struct octree_node { + struct octree_index ichild_node; /* Index of the 1st child node */ + struct octree_index ichild_attr; /* Index of the 1st child attr */ + uint8_t is_valid; /* Mask defining whether the children are valid or not */ + uint8_t is_leaf; /* Mask defining whether the children are leaves or not */ + ALIGN(16) char data[8][HTVOX_MAX_SIZEOF_VOXEL]; /* Data of the leaves */ +}; + +/* Stacked children of an octree node */ +struct stack { + struct octree_node nodes[8]; /* List of registered children nodes */ + uint8_t mask; /* Mask of valid children nodes (0 = empty) */ +}; + +struct octree_builder { + struct stack stacks[OCTREE_DEPTH_MAX]; + struct octree_buffer* buffer; + + const struct htvox_voxel_desc* desc; + size_t nleaves; /* Number of emitted leaves */ + + int octree_depth; /* Maximum octree depth */ + int non_empty_lvl; /* Index of the 1st non empty level */ + + uint64_t mcode; /* Morton code of the last registered voxels */ +}; + +/******************************************************************************* + * Helper function + ******************************************************************************/ +#ifndef NDEBUG +static void +check_octree(struct octree_buffer* buf, const struct octree_index root) +{ + const struct octree_xnode* node; + int ichild; + ASSERT(buf); + + node = octree_buffer_get_node(buf, root); + FOR_EACH(ichild, 0, 8) { + const int ichild_flag = BIT(ichild); + if((node->is_valid & ichild_flag) == 0) continue; + + if(node->is_leaf & ichild_flag) { + struct octree_index iattr; + iattr = octree_buffer_get_child_attr_index(buf, root, ichild); + ASSERT(octree_buffer_get_attr(buf, iattr) != NULL); + } else { + struct octree_index child; + child = octree_buffer_get_child_node_index(buf, root, ichild); + check_octree(buf, child); + } + } +} +#endif + +static INLINE int +check_htvox_voxel_desc(const struct htvox_voxel_desc* desc) +{ + return desc + && desc->get + && desc->merge + && desc->challenge_merge + && desc->size > 0 + && desc->size <= HTVOX_MAX_SIZEOF_VOXEL; +} + +static INLINE void +stack_clear(struct stack* stack) +{ + int inode; + FOR_EACH(inode, 0, 8) { + int ichild; + stack->nodes[inode].is_leaf = 0; + stack->nodes[inode].is_valid = 0; + stack->nodes[inode].ichild_node = OCTREE_INDEX_NULL; + stack->nodes[inode].ichild_attr = OCTREE_INDEX_NULL; + FOR_EACH(ichild, 0, 8) { + memset(stack->nodes[inode].data[ichild], 0, HTVOX_MAX_SIZEOF_VOXEL); + } + } + stack->mask = 0; +} + +/* Build a parent octree node from the registered stack nodes */ +static INLINE void +stack_setup_node + (struct stack* stack, + struct octree_node* node, + const struct htvox_voxel_desc* desc) +{ + int ichild; + ASSERT(stack && node && check_htvox_voxel_desc(desc)); + + node->ichild_node = OCTREE_INDEX_NULL; + node->ichild_attr = OCTREE_INDEX_NULL; + node->is_valid = stack->mask; + node->is_leaf = 0; + + if(stack->mask == 0) return; /* Empty stack */ + + /* Try to merge the child's leaves */ + FOR_EACH(ichild, 0, 8) { + const void* data[8]; + const uint8_t ichild_flag = (uint8_t)BIT(ichild); + struct octree_node* child = stack->nodes + ichild; + int igrandchild; + + /* Fetch the grandchildren data */ + FOR_EACH(igrandchild, 0, 8) { + data[igrandchild] = child->data[igrandchild]; + } + + desc->merge(node->data[ichild], data, 8, desc->context); + + if(child->is_leaf==0xFF && desc->challenge_merge(data, 8, desc->context)) { + /* The node becomes a leaf : the children does not exist anymore */ + node->is_leaf |= ichild_flag; + stack->mask ^= ichild_flag; + } + } +} + +static res_T +stack_write + (struct stack* stack, /* Node to write */ + struct octree_buffer* buf, /* Buffer where nodes are written */ + struct octree_index* out_index, /* Index of the first written node */ + size_t* out_nleaves) /* #writen leaves */ +{ + struct octree_index nodes_id = OCTREE_INDEX_NULL; + struct octree_node* node = NULL; + size_t nleaves = 0; + int inode; + res_T res = RES_OK; + ASSERT(stack && buf && out_index && out_nleaves); + + /* No registered nodes, this means that the nodes were merged in an higher + * level */ + if(!stack->mask) goto exit; + + /* Write the attrib of the children */ + FOR_EACH(inode, 0, 8) { + char* data = NULL; + size_t nattrs = 0; + size_t nvoxs = 0; + int ichild = 0; + + if((stack->mask & BIT(inode)) == 0) continue; /* Empty node */ + + node = stack->nodes + inode; + + nattrs = (size_t)popcount(node->is_valid); + ASSERT(nattrs > 0); + + res = octree_buffer_alloc_attrs(buf, nattrs, &node->ichild_attr); + if(res != RES_OK) goto error; + + data = octree_buffer_get_attr(buf, node->ichild_attr); + nvoxs = 0; + FOR_EACH(ichild, 0, 8) { + if(!(node->is_valid & BIT(ichild))) continue; + memcpy(data + nvoxs*buf->voxsize, node->data[ichild], buf->voxsize); + ++nvoxs; + } + ASSERT(nvoxs == nattrs); + nleaves += (size_t)popcount(node->is_leaf); + } + + do { + struct octree_index index = OCTREE_INDEX_NULL; + struct octree_xnode* xnodes = NULL; + const size_t nnodes = (size_t)popcount(stack->mask); + size_t ixnode = 0; + + /* Alloc the octree nodes */ + res = octree_buffer_alloc_nodes(buf, (size_t)nnodes, &nodes_id); + if(res != RES_OK) goto error; + xnodes = octree_buffer_get_node(buf, nodes_id); + + FOR_EACH(inode, 0, 8) { + uint16_t attr_offset = UINT16_MAX; + uint16_t node_offset = UINT16_MAX; + + if((stack->mask & BIT(inode)) == 0) continue; /* Empty node */ + node = stack->nodes + inode; + + /* Setup the offset toward the children and children attribs */ + if(node->ichild_node.ipage == nodes_id.ipage) { + node_offset = node->ichild_node.inode; + attr_offset = node->ichild_attr.inode; + } + + /* The page id of the children is not the same as that of node */ + if(node_offset > OCTREE_XNODE_MAX_CHILDREN_OFFSET) { + res = octree_buffer_alloc_far_index(buf, &index); + if(res != RES_OK) break; + *octree_buffer_get_far_index(buf, index) = node->ichild_node; + node_offset = OCTREE_XNODE_FLAG_FAR_INDEX | index.inode; + } + + /* The page id of the attribs is not tthe same as that of node */ + if(attr_offset > OCTREE_XNODE_FLAG_FAR_INDEX) { + res = octree_buffer_alloc_far_index(buf, &index); + if(res != RES_OK) break; + *octree_buffer_get_far_index(buf, index) = node->ichild_attr; + attr_offset = OCTREE_XNODE_FLAG_FAR_INDEX | index.inode; + } + + xnodes[ixnode].node_offset = node_offset; + xnodes[ixnode].attr_offset = attr_offset; + xnodes[ixnode].is_valid = node->is_valid; + xnodes[ixnode].is_leaf = node->is_leaf; + ++ixnode; + } + /* inode < 8 <=> not enough memory in the current page. A far index could not + * be stored in the same page of its associated node. The write process was + * stoped. Rewrite the whole stacked nodes in a new page. */ + } while(inode < 8); + + +exit: + /* Return the index toward the first writen nodes */ + *out_index = nodes_id; + *out_nleaves = nleaves; + return res; +error: + goto exit; +} + +static res_T +octree_builder_init + (struct octree_builder* bldr, + const size_t definition, + const struct htvox_voxel_desc* desc, + struct octree_buffer* buffer) +{ + int ilvl; + res_T res = RES_OK; + ASSERT(bldr && IS_POW2(definition) && check_htvox_voxel_desc(desc)); + memset(bldr, 0, sizeof(struct octree_builder)); + + /* Compute the maximum depth of the octree */ + bldr->octree_depth = log2i((int)definition); + if(bldr->octree_depth > OCTREE_DEPTH_MAX) { + res = RES_MEM_ERR; + goto error; + } + + /* Init the per octree level stack */ + FOR_EACH(ilvl, 0, bldr->octree_depth) { + stack_clear(&bldr->stacks[ilvl]); + } + + octree_buffer_clear(buffer); + bldr->nleaves = 0; + bldr->desc = desc; + bldr->buffer = buffer; + bldr->non_empty_lvl = bldr->octree_depth - 1; + +exit: + return res; +error: + goto exit; +} + +static res_T +octree_builder_add_voxel + (struct octree_builder* bldr, + const struct voxel* vox) +{ + uint64_t mcode_xor; + int inode; + int ichild; + uint8_t ichild_flag; + res_T res = RES_OK; + ASSERT(bldr && vox && vox->mcode >= bldr->mcode); + + /* Define if the bits in [4 .. 63] of the previous and the next Morton + * codes are the same */ + mcode_xor = bldr->mcode ^ vox->mcode; + + /* The next voxel is not in the current node */ + if(mcode_xor >= 8) { + size_t ilvl; + + inode = (bldr->mcode >> 3) & 7; + bldr->stacks[0].mask |= (uint8_t)BIT(inode); + + /* Flush the stack of the octree level that does not contain the next + * voxel. The last octree level is actually the level that contains the + * root node while the penultimate describes the root node itself. These 2 + * levels contain all the voxels and can be skipped */ + FOR_EACH(ilvl, 0, bldr->octree_depth-2/*The 2 last leaves contain all voxels*/) { + struct octree_node* stack_node; + uint64_t mcode_max_lvl; + size_t nleaves; + + /* Compute the maximum morton code value for the current octree level */ + mcode_max_lvl = 8lu/*#children*/ * (1lu<<(3*(ilvl+1)))/*#voxels per children*/; + + if(mcode_xor < mcode_max_lvl) break; + + /* Retrieve the node index of the next level */ + inode = (bldr->mcode >> (3*(ilvl+2))) & 7; + + /* The next voxel is not in the ilvl^th stack. Setup the parent node of the + * nodes registered into the stack */ + stack_node = &bldr->stacks[ilvl+1].nodes[inode]; + stack_setup_node(&bldr->stacks[ilvl], stack_node, bldr->desc); + bldr->stacks[ilvl+1].mask |= (uint8_t)BIT(inode); + + /* Write the nodes of the stack of the current octree level into the buf */ + res = stack_write + (&bldr->stacks[ilvl], bldr->buffer, &stack_node->ichild_node, &nleaves); + if(res != RES_OK) goto error; + + bldr->nleaves += nleaves; + if(nleaves) bldr->non_empty_lvl = MMIN(bldr->non_empty_lvl, (int)ilvl); + + /* Reset the current stack */ + stack_clear(&bldr->stacks[ilvl]); + } + } + + /* Retrieve the index of the current voxel and of its parent */ + ichild = vox->mcode & 7; + inode = (vox->mcode >> 3) & 7; + ichild_flag = (uint8_t)BIT(ichild); + + /* Register the voxel */ + memcpy(bldr->stacks[0].nodes[inode].data[ichild], vox->data, bldr->desc->size); + bldr->stacks[0].nodes[inode].is_valid |= ichild_flag; + bldr->stacks[0].nodes[inode].is_leaf |= ichild_flag; + + /* Update morton code of the last registered voxel */ + bldr->mcode = vox->mcode; + +exit: + return res; +error: + goto exit; +} + +static INLINE res_T +octree_builder_finalize + (struct octree_builder* bldr, + struct octree_index* root_id, + void* root_data) +{ + const void* data[8]; + size_t inode; + size_t nleaves; + int ilvl; + int ichild; + res_T res = RES_OK; + ASSERT(bldr); + + inode = (bldr->mcode >> 3) & 7; + bldr->stacks[0].mask |= (uint8_t)BIT(inode); + + /* Flush the stacked nodes */ + FOR_EACH(ilvl, 0, bldr->octree_depth-1) { + struct octree_node* parent_node; + + if(bldr->stacks[ilvl].mask == 0) continue; + + /* Retrieve the node index of the next level */ + inode = (bldr->mcode >> (3*(ilvl+2))) & 7; + + /* Setup the parent node of the nodes registered into the current stack */ + parent_node = &bldr->stacks[ilvl+1].nodes[inode]; /* Fetch the parent node */ + stack_setup_node(&bldr->stacks[ilvl], parent_node, bldr->desc); + bldr->stacks[ilvl+1].mask |= (uint8_t)BIT(inode); + + /* Write the stacked nodes of the current level */ + res = stack_write + (&bldr->stacks[ilvl], bldr->buffer, &parent_node->ichild_node, &nleaves); + if(res != RES_OK) goto error; + bldr->nleaves += nleaves; + } + + ilvl = bldr->octree_depth-1; /* Root level */ + + /* Write the root node */ + res = stack_write(&bldr->stacks[ilvl], bldr->buffer, root_id, &nleaves); + if(res != RES_OK) goto error; + bldr->nleaves += nleaves; + + /* Setup the root attribs */ + ASSERT(bldr->stacks[ilvl].mask == 1); /* Only the root node is active */ + FOR_EACH(ichild, 0, 8) { + data[ichild] = bldr->stacks[ilvl].nodes[0].data[ichild]; + } + bldr->desc->merge(root_data, data, 8, bldr->desc->context); + +exit: + return res; +error: + goto exit; +} + +static void +octree_release(ref_T* ref) +{ + struct htvox_octree* oct; + struct htvox_device* dev; + ASSERT(ref); + oct = CONTAINER_OF(ref, struct htvox_octree, ref); + octree_buffer_release(&oct->buffer); + dev = oct->dev; + MEM_RM(dev->allocator, oct); + HTVOX(device_ref_put(dev)); +} + +/******************************************************************************* + * Exported functions + ******************************************************************************/ +res_T +htvox_octree_create + (struct htvox_device* dev, + const double lower[3], /* Lower bound of the octree */ + const double upper[3], /* Upper bound of the octree */ + const size_t nvoxels[3], /* # voxels along the 3 axis */ + const struct htvox_voxel_desc* desc, /* Descriptor of a voxel */ + struct htvox_octree** out_oct) +{ + struct htvox_octree* oct = NULL; + double vox_sz[3]; /* World space size of a voxel */ + struct octree_builder bldr; + struct voxel vox = VOXEL_NULL; + uint64_t mcode_max; + uint64_t mcode; + res_T res = RES_OK; + + if(!dev || !lower || !upper || !nvoxels + || !check_htvox_voxel_desc(desc) || !out_oct) { + res = RES_BAD_ARG; + goto error; + } + if(lower[0] >= upper[0] + || lower[1] >= upper[1] + || lower[2] >= upper[2]) { + log_err(dev, + "%s: the submitted AABB is degenerated.\n" + "\tlower={%g, %g, %g}, upper={%g, %g, %g}.\n", + FUNC_NAME, SPLIT3(lower), SPLIT3(upper)); + res = RES_BAD_ARG; + goto error; + } + if(!nvoxels[0] || !nvoxels[1] || !nvoxels[2]) { + log_err(dev, + "%s: the number of voxels along each axis cannot be null.\n" + "\t#voxels XYZ = {%lu, %lu, %lu}.\n", + FUNC_NAME, + (unsigned long)nvoxels[0], + (unsigned long)nvoxels[1], + (unsigned long)nvoxels[2]); + res = RES_BAD_ARG; + goto error; + } + oct = MEM_CALLOC(dev->allocator, 1, sizeof(struct htvox_octree)); + if(!oct) { + res = RES_MEM_ERR; + goto error; + } + ref_init(&oct->ref); + HTVOX(device_ref_get(dev)); + octree_buffer_init(dev->allocator, desc->size, &oct->buffer); + oct->dev = dev; + + /* Compute the octree definition */ + oct->definition = MMAX(nvoxels[0], 2); + oct->definition = MMAX(nvoxels[1], oct->definition); + oct->definition = MMAX(nvoxels[2], oct->definition); + oct->definition = round_up_pow2(oct->definition); + + /* Intialize the octree builder */ + res = octree_builder_init(&bldr, oct->definition, desc, &oct->buffer); + if(res != RES_OK) goto error; + + vox.data = MEM_CALLOC(dev->allocator, 1, desc->size); + if(!vox.data) { + res = RES_MEM_ERR; + goto error; + } + + mcode_max = oct->definition * oct->definition * oct->definition; + FOR_EACH(mcode, 0, mcode_max) { + size_t xyz[3]; + uint32_t ui3[3]; + + morton_xyz_decode_u21(mcode, ui3); + + /* Out of bound voxels */ + if(ui3[0] >= nvoxels[0] + || ui3[1] >= nvoxels[1] + || ui3[2] >= nvoxels[2]) + continue; + + /* Retrieve the voxel data from the caller */ + xyz[0] = (size_t)ui3[0]; + xyz[1] = (size_t)ui3[1]; + xyz[2] = (size_t)ui3[2]; + desc->get(xyz, vox.data, desc->context); + vox.mcode = mcode; + + /* Register the voxel against the octree */ + res = octree_builder_add_voxel(&bldr, &vox); + if(res != RES_OK) goto error; + } + + res = octree_builder_finalize(&bldr, &oct->root, oct->root_attr); + if(res != RES_OK) goto error; + + oct->nleaves = bldr.nleaves; + oct->depth = (size_t)(bldr.octree_depth - bldr.non_empty_lvl) + + 1 /* leaf level */; + ASSERT(bldr.octree_depth > bldr.non_empty_lvl); + +#ifndef NDEBUG + check_octree(&oct->buffer, oct->root); +#endif + /* Setup the octree AABB in world space */ + d3_set(oct->lower, lower); + d3_set(oct->upper, upper); + + /* Compute the voxel size in world space */ + vox_sz[0] = (upper[0] - lower[0]) / (double)nvoxels[0]; + vox_sz[1] = (upper[1] - lower[1]) / (double)nvoxels[1]; + vox_sz[2] = (upper[2] - lower[2]) / (double)nvoxels[2]; + + /* Compute the octree AABB in world space */ + d3_set(oct->oclow, lower); + oct->ocupp[0] = (double)oct->definition * vox_sz[0]; + oct->ocupp[1] = (double)oct->definition * vox_sz[1]; + oct->ocupp[2] = (double)oct->definition * vox_sz[2]; + + /* Adjust the octree definition with respect to the octree depth since finest + * levels might be removed by the build process due to the merging process */ + oct->definition = oct->definition / (1u << bldr.non_empty_lvl); + +exit: + if(vox.data) MEM_RM(dev->allocator, vox.data); + if(out_oct) *out_oct = oct; + return res; +error: + if(oct) { + HTVOX(octree_ref_put(oct)); + oct = NULL; + } + goto exit; +} + +res_T +htvox_octree_ref_get(struct htvox_octree* oct) +{ + if(!oct) return RES_BAD_ARG; + ref_get(&oct->ref); + return RES_OK; +} + +res_T +htvox_octree_ref_put(struct htvox_octree* oct) +{ + if(!oct) return RES_BAD_ARG; + ref_put(&oct->ref, octree_release); + return RES_OK; +} + +res_T +htvox_octree_get_desc + (const struct htvox_octree* oct, struct htvox_octree_desc* desc) +{ + size_t nvoxs_per_page; + if(!oct || !desc) return RES_BAD_ARG; + nvoxs_per_page = oct->buffer.pagesize / oct->buffer.voxsize; + d3_set(desc->lower, oct->lower); + d3_set(desc->upper, oct->upper); + desc->nleaves = oct->nleaves; + desc->nvoxels = oct->buffer.attr_head.ipage * nvoxs_per_page + + oct->buffer.attr_head.inode + + 1; /* Root node */ + desc->depth = oct->depth; + return RES_OK; +} + +res_T +htvox_octree_for_each_leaf + (struct htvox_octree* oct, + void (*func) + (const void* val, + const size_t ileaf, + const double low[3], + const double upp[3], + void* ctx), + void* ctx) +{ + struct stack_entry { + struct octree_index inode; + double low[3]; + double upp[3]; + } stack[OCTREE_DEPTH_MAX*8]; + int istack; + size_t ileaf = 0; + + if(!oct || !func) return RES_BAD_ARG; + + stack[0].inode = oct->root; + stack[0].low[0] = oct->oclow[0]; + stack[0].low[1] = oct->oclow[1]; + stack[0].low[2] = oct->oclow[2]; + stack[0].upp[0] = oct->ocupp[0]; + stack[0].upp[1] = oct->ocupp[1]; + stack[0].upp[2] = oct->ocupp[2]; + istack = 1; + + do { + const struct stack_entry entry = stack[--istack]; + struct octree_xnode* node; + int ichild; + + node = octree_buffer_get_node(&oct->buffer, entry.inode); + + FOR_EACH(ichild, 0, 8) { + const uint8_t ichild_flag = (uint8_t)BIT(ichild); + double half_sz[3]; /* Half size of the current node */ + double mid[3]; /* Middle point of the current node */ + double upp[3]; + double low[3]; + + if((node->is_valid & ichild_flag) == 0) continue; /* Empty node */ + + half_sz[0] = (entry.upp[0] - entry.low[0])*0.5; + half_sz[1] = (entry.upp[1] - entry.low[1])*0.5; + half_sz[2] = (entry.upp[2] - entry.low[2])*0.5; + mid[0] = entry.low[0] + half_sz[0]; + mid[1] = entry.low[1] + half_sz[1]; + mid[2] = entry.low[2] + half_sz[2]; + low[0] = ichild&4 ? mid[0] : entry.low[0]; + low[1] = ichild&2 ? mid[1] : entry.low[1]; + low[2] = ichild&1 ? mid[2] : entry.low[2]; + upp[0] = low[0] + half_sz[0]; + upp[1] = low[1] + half_sz[1]; + upp[2] = low[2] + half_sz[2]; + + if(node->is_leaf & ichild_flag) { + struct octree_index iattr; + const void* val; + + iattr = octree_buffer_get_child_attr_index + (&oct->buffer, entry.inode, ichild); + val = octree_buffer_get_attr(&oct->buffer, iattr); + + ASSERT(upp[0] <= oct->upper[0]); + ASSERT(upp[1] <= oct->upper[1]); + ASSERT(upp[2] <= oct->upper[2]); + ASSERT(low[0] >= oct->lower[0]); + ASSERT(low[1] >= oct->lower[1]); + ASSERT(low[2] >= oct->lower[2]); + + func(val, ileaf, low, upp, ctx); + ileaf++; + } else { + struct stack_entry* top = stack + istack; + + top->inode = octree_buffer_get_child_node_index + (&oct->buffer, entry.inode, ichild); + top->low[0] = low[0]; + top->low[1] = low[1]; + top->low[2] = low[2]; + top->upp[0] = upp[0]; + top->upp[1] = upp[1]; + top->upp[2] = upp[2]; + ++istack; + } + } + } while(istack); + + return RES_OK; +} + +res_T +htvox_octree_at + (struct htvox_octree* oct, + const double position[3], + struct htvox_voxel* voxel) +{ + struct octree_index inode; + double scale_exp2; + double low[3]; + double pos[3]; + + if(!oct || !position || !voxel) return RES_BAD_ARG; + + *voxel = HTVOX_VOXEL_NULL; + + if(position[0] > oct->upper[0] || position[0] < oct->lower[0] + || position[1] > oct->upper[1] || position[1] < oct->lower[1] + || position[2] > oct->upper[2] || position[2] < oct->lower[2]) { + return RES_OK; + } + + /* Transform the position in the normalized octree space, + * i.e. octree lies in [0, 1]^2 */ + pos[0] = (position[0] - oct->oclow[0]) / (oct->ocupp[0] - oct->oclow[0]); + pos[1] = (position[1] - oct->oclow[1]) / (oct->ocupp[1] - oct->oclow[1]); + pos[2] = (position[2] - oct->oclow[2]) / (oct->ocupp[2] - oct->oclow[2]); + + /* Initialized the lower left corner of the current node */ + low[0] = 0; + low[1] = 0; + low[2] = 0; + + scale_exp2 = 0.5; + inode = oct->root; + *voxel = HTVOX_VOXEL_NULL; + for(;;) { + struct octree_xnode* node = octree_buffer_get_node(&oct->buffer, inode); + int ichild; + uint8_t ichild_flag; + double mid[3]; + + /* Compute the middle point of the node */ + mid[0] = low[0] + scale_exp2; + mid[1] = low[1] + scale_exp2; + mid[2] = low[2] + scale_exp2; + + /* Define the child of the current node into which pos `lies' and compute + * its lower left corner */ + ichild = 0; + if(pos[0] > mid[0]) { ichild |= 4; low[0] = mid[0]; } + if(pos[1] > mid[1]) { ichild |= 2; low[1] = mid[1]; } + if(pos[2] > mid[2]) { ichild |= 1; low[2] = mid[2]; } + + ichild_flag = (uint8_t)BIT(ichild); + if((node->is_valid & ichild_flag) == 0) { /* Empty node */ + break; + } else if(node->is_leaf & ichild_flag) { /* Leaf node */ + inode = octree_buffer_get_child_attr_index(&oct->buffer, inode, ichild); + voxel->data = octree_buffer_get_attr(&oct->buffer, inode); + voxel->id = inode.ipage * oct->buffer.pagesize / oct->buffer.voxsize + + inode.inode; + break; + } else { /* Child node */ + inode = octree_buffer_get_child_node_index(&oct->buffer, inode, ichild); + scale_exp2 *= 0.5; + } + } + return RES_OK; +} + diff --git a/src/htvox_octree.h b/src/htvox_octree.h @@ -0,0 +1,40 @@ +/* Copyright (C) CNRS 2018 + * + * 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 HTVOX_OCTREE_H +#define HTVOX_OCTREE_H + +#include "htvox_octree_buffer.h" +#include <rsys/ref_count.h> + +struct htvox_octree { + size_t definition; /* #voxels along the X, Y and Z axis */ + + double lower[3], upper[3]; /* Submitted World space AABB of the octree */ + double oclow[3], ocupp[3]; /* Adjusted World space AABB of the octree */ + + struct octree_buffer buffer; + struct octree_index root; /* Index toward the children of the root */ + ALIGN(16) char root_attr[HTVOX_MAX_SIZEOF_VOXEL]; /* Attribute of the root */ + + size_t nleaves; /* #leaves */ + size_t depth; /* Maximum depth of the octree */ + + struct htvox_device* dev; + ref_T ref; +}; + +#endif /* HTVOX_OCTREE_H */ + diff --git a/src/htvox_octree_buffer.h b/src/htvox_octree_buffer.h @@ -30,8 +30,8 @@ * stored into the same page, that defines the absolute position of its first * valid child into the whole list of node pages * - * The data of the nodes are stored in their own of memory pages. The attribs - * of the children of a node are stored consecutively into a page. If the page + * The data of the nodes are stored in their own memory pages. The attribs of + * the children of a node are stored consecutively into a page. If the page * identifier of the attribs is the same of the page into which their parent * node lies, then the node saves the index toward the first valid attrib into * the page of attribs. In the other case, the node references a `struct diff --git a/src/htvox_scene.c b/src/htvox_scene.c @@ -1,765 +0,0 @@ -/* Copyright (C) CNRS 2018 - * - * 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/>. */ - -#include "htvox.h" -#include "htvox_c.h" -#include "htvox_device.h" -#include "htvox_scene.h" - -#include <rsys/double3.h> -#include <rsys/mem_allocator.h> - -#define OCTREE_DEPTH_MAX 16 /* Maximum depth of an octree */ - -struct voxel { - uint64_t mcode; /* Morton code of the voxel */ - void* data; /* Data of the voxel */ -}; -static const struct voxel VOXEL_NULL = {0, NULL}; - -struct octree_node { - struct octree_index ichild_node; /* Index of the 1st child node */ - struct octree_index ichild_attr; /* Index of the 1st child attr */ - uint8_t is_valid; /* Mask defining whether the children are valid or not */ - uint8_t is_leaf; /* Mask defining whether the children are leaves or not */ - ALIGN(16) char data[8][HTVOX_MAX_SIZEOF_VOXEL]; /* Data of the leaves */ -}; - -/* Stacked children of an octree node */ -struct stack { - struct octree_node nodes[8]; /* List of registered children nodes */ - uint8_t mask; /* Mask of valid children nodes (0 = empty) */ -}; - -struct octree_builder { - struct stack stacks[OCTREE_DEPTH_MAX]; - struct octree_buffer* buffer; - - const struct htvox_voxel_desc* desc; - size_t nleaves; - - int octree_depth; - - uint64_t mcode; /* Morton code of the last registered voxels */ -}; - -/******************************************************************************* - * Helper function - ******************************************************************************/ -#ifndef NDEBUG -static void -check_octree(struct octree_buffer* buf, const struct octree_index root) -{ - const struct octree_xnode* node; - int ichild; - ASSERT(buf); - - node = octree_buffer_get_node(buf, root); - FOR_EACH(ichild, 0, 8) { - const int ichild_flag = BIT(ichild); - if((node->is_valid & ichild_flag) == 0) continue; - - if(node->is_leaf & ichild_flag) { - struct octree_index iattr; - iattr = octree_buffer_get_child_attr_index(buf, root, ichild); - ASSERT(octree_buffer_get_attr(buf, iattr) != NULL); - } else { - struct octree_index child; - child = octree_buffer_get_child_node_index(buf, root, ichild); - check_octree(buf, child); - } - } -} -#endif - -static INLINE int -check_htvox_voxel_desc(const struct htvox_voxel_desc* desc) -{ - return desc - && desc->get - && desc->merge - && desc->challenge_merge - && desc->size > 0 - && desc->size <= HTVOX_MAX_SIZEOF_VOXEL; -} - -static INLINE void -stack_clear(struct stack* stack) -{ - int inode; - FOR_EACH(inode, 0, 8) { - int ichild; - stack->nodes[inode].is_leaf = 0; - stack->nodes[inode].is_valid = 0; - stack->nodes[inode].ichild_node = OCTREE_INDEX_NULL; - stack->nodes[inode].ichild_attr = OCTREE_INDEX_NULL; - FOR_EACH(ichild, 0, 8) { - memset(stack->nodes[inode].data[ichild], 0, HTVOX_MAX_SIZEOF_VOXEL); - } - } - stack->mask = 0; -} - -/* Build a parent octree node from the registered stack nodes */ -static INLINE void -stack_setup_node - (struct stack* stack, - struct octree_node* node, - const struct htvox_voxel_desc* desc) -{ - int ichild; - ASSERT(stack && node && check_htvox_voxel_desc(desc)); - - node->ichild_node = OCTREE_INDEX_NULL; - node->ichild_attr = OCTREE_INDEX_NULL; - node->is_valid = stack->mask; - node->is_leaf = 0; - - if(stack->mask == 0) return; /* Empty stack */ - - /* Try to merge the child's leaves */ - FOR_EACH(ichild, 0, 8) { - const void* data[8]; - const uint8_t ichild_flag = (uint8_t)BIT(ichild); - struct octree_node* child = stack->nodes + ichild; - int igrandchild; - - /* Fetch the grandchildren data */ - FOR_EACH(igrandchild, 0, 8) { - data[igrandchild] = child->data[igrandchild]; - } - - desc->merge(node->data[ichild], data, 8, desc->context); - - if(child->is_leaf==0xFF && desc->challenge_merge(data, 8, desc->context)) { - /* The node becomes a leaf : the children does not exist anymore */ - node->is_leaf |= ichild_flag; - stack->mask ^= ichild_flag; - } - } -} - -static res_T -stack_write - (struct stack* stack, /* Node to write */ - struct octree_buffer* buf, /* Buffer where nodes are written */ - struct octree_index* out_index, /* Index of the first written node */ - size_t* out_nleaves) /* #writen leaves */ -{ - struct octree_index nodes_id = OCTREE_INDEX_NULL; - struct octree_node* node = NULL; - size_t nleaves = 0; - int inode; - res_T res = RES_OK; - ASSERT(stack && buf && out_index && out_nleaves); - - /* No registered nodes, this means that the nodes were merged in an higher - * level */ - if(!stack->mask) goto exit; - - /* Write the attrib of the children */ - FOR_EACH(inode, 0, 8) { - char* data = NULL; - size_t nattrs = 0; - size_t nvoxs = 0; - int ichild = 0; - - if((stack->mask & BIT(inode)) == 0) continue; /* Empty node */ - - node = stack->nodes + inode; - - nattrs = (size_t)popcount(node->is_valid); - ASSERT(nattrs > 0); - - res = octree_buffer_alloc_attrs(buf, nattrs, &node->ichild_attr); - if(res != RES_OK) goto error; - - data = octree_buffer_get_attr(buf, node->ichild_attr); - nvoxs = 0; - FOR_EACH(ichild, 0, 8) { - if(!(node->is_valid & BIT(ichild))) continue; - memcpy(data + nvoxs*buf->voxsize, node->data[ichild], buf->voxsize); - ++nvoxs; - } - ASSERT(nvoxs == nattrs); - nleaves += (size_t)popcount(node->is_leaf); - } - - do { - struct octree_index index = OCTREE_INDEX_NULL; - struct octree_xnode* xnodes = NULL; - const size_t nnodes = (size_t)popcount(stack->mask); - size_t ixnode = 0; - - /* Alloc the octree nodes */ - res = octree_buffer_alloc_nodes(buf, (size_t)nnodes, &nodes_id); - if(res != RES_OK) goto error; - xnodes = octree_buffer_get_node(buf, nodes_id); - - FOR_EACH(inode, 0, 8) { - uint16_t attr_offset = UINT16_MAX; - uint16_t node_offset = UINT16_MAX; - - if((stack->mask & BIT(inode)) == 0) continue; /* Empty node */ - node = stack->nodes + inode; - - /* Setup the offset toward the children and children attribs */ - if(node->ichild_node.ipage == nodes_id.ipage) { - node_offset = node->ichild_node.inode; - attr_offset = node->ichild_attr.inode; - } - - /* The page id of the children is not the same as that of node */ - if(node_offset > OCTREE_XNODE_MAX_CHILDREN_OFFSET) { - res = octree_buffer_alloc_far_index(buf, &index); - if(res != RES_OK) break; - *octree_buffer_get_far_index(buf, index) = node->ichild_node; - node_offset = OCTREE_XNODE_FLAG_FAR_INDEX | index.inode; - } - - /* The page id of the attribs is not tthe same as that of node */ - if(attr_offset > OCTREE_XNODE_FLAG_FAR_INDEX) { - res = octree_buffer_alloc_far_index(buf, &index); - if(res != RES_OK) break; - *octree_buffer_get_far_index(buf, index) = node->ichild_attr; - attr_offset = OCTREE_XNODE_FLAG_FAR_INDEX | index.inode; - } - - xnodes[ixnode].node_offset = node_offset; - xnodes[ixnode].attr_offset = attr_offset; - xnodes[ixnode].is_valid = node->is_valid; - xnodes[ixnode].is_leaf = node->is_leaf; - ++ixnode; - } - /* inode < 8 <=> not enough memory in the current page. A far index could not - * be stored in the same page of its associated node. The write process was - * stoped. Rewrite the whole stacked nodes in a new page. */ - } while(inode < 8); - - -exit: - /* Return the index toward the first writen nodes */ - *out_index = nodes_id; - *out_nleaves = nleaves; - return res; -error: - goto exit; -} - -static res_T -octree_builder_init - (struct octree_builder* bldr, - const size_t definition, - const struct htvox_voxel_desc* desc, - struct octree_buffer* buffer) -{ - int ilvl; - res_T res = RES_OK; - ASSERT(bldr && IS_POW2(definition) && check_htvox_voxel_desc(desc)); - memset(bldr, 0, sizeof(struct octree_builder)); - - /* Compute the maximum depth of the octree */ - bldr->octree_depth = log2i((int)definition); - if(bldr->octree_depth > OCTREE_DEPTH_MAX) { - res = RES_MEM_ERR; - goto error; - } - - /* Init the per octree level stack */ - FOR_EACH(ilvl, 0, bldr->octree_depth) { - stack_clear(&bldr->stacks[ilvl]); - } - - octree_buffer_clear(buffer); - bldr->nleaves = 0; - bldr->desc = desc; - bldr->buffer = buffer; - -exit: - return res; -error: - goto exit; -} - -static res_T -octree_builder_add_voxel - (struct octree_builder* bldr, - const struct voxel* vox) -{ - uint64_t mcode_xor; - int inode; - int ichild; - uint8_t ichild_flag; - res_T res = RES_OK; - ASSERT(bldr && vox && vox->mcode >= bldr->mcode); - - /* Define if the bits in [4 .. 63] of the previous and the next Morton - * codes are the same */ - mcode_xor = bldr->mcode ^ vox->mcode; - - /* The next voxel is not in the current node */ - if(mcode_xor >= 8) { - size_t ilvl; - - inode = (bldr->mcode >> 3) & 7; - bldr->stacks[0].mask |= (uint8_t)BIT(inode); - - /* Flush the stack of the octree level that does not contain the next - * voxel. The last octree level is actually the level that contains the - * root node while the penultimate describes the root node itself. These 2 - * levels contain all the voxels and can be skipped */ - FOR_EACH(ilvl, 0, bldr->octree_depth-2/*The 2 last leaves contain all voxels*/) { - struct octree_node* stack_node; - uint64_t mcode_max_lvl; - size_t nleaves; - - /* Compute the maximum morton code value for the current octree level */ - mcode_max_lvl = 8lu/*#children*/ * (1lu<<(3*(ilvl+1)))/*#voxels per children*/; - - if(mcode_xor < mcode_max_lvl) break; - - /* Retrieve the node index of the next level */ - inode = (bldr->mcode >> (3*(ilvl+2))) & 7; - - /* The next voxel is not in the ilvl^th stack. Setup the parent node of the - * nodes registered into the stack */ - stack_node = &bldr->stacks[ilvl+1].nodes[inode]; - stack_setup_node(&bldr->stacks[ilvl], stack_node, bldr->desc); - bldr->stacks[ilvl+1].mask |= (uint8_t)BIT(inode); - - /* Write the nodes of the stack of the current octree level into the buf */ - res = stack_write - (&bldr->stacks[ilvl], bldr->buffer, &stack_node->ichild_node, &nleaves); - if(res != RES_OK) goto error; - - bldr->nleaves += nleaves; - - /* Reset the current stack */ - stack_clear(&bldr->stacks[ilvl]); - } - } - - /* Retrieve the index of the current voxel and of its parent */ - ichild = vox->mcode & 7; - inode = (vox->mcode >> 3) & 7; - ichild_flag = (uint8_t)BIT(ichild); - - /* Register the voxel */ - memcpy(bldr->stacks[0].nodes[inode].data[ichild], vox->data, bldr->desc->size); - bldr->stacks[0].nodes[inode].is_valid |= ichild_flag; - bldr->stacks[0].nodes[inode].is_leaf |= ichild_flag; - - /* Update morton code of the last registered voxel */ - bldr->mcode = vox->mcode; - -exit: - return res; -error: - goto exit; -} - -static INLINE res_T -octree_builder_finalize - (struct octree_builder* bldr, - struct octree_index* root_id) -{ - size_t inode; - size_t nleaves; - int ilvl; - res_T res = RES_OK; - ASSERT(bldr); - - inode = (bldr->mcode >> 3) & 7; - bldr->stacks[0].mask |= (uint8_t)BIT(inode); - - /* Flush the stacked nodes */ - FOR_EACH(ilvl, 0, bldr->octree_depth-1) { - struct octree_node* parent_node; - - if(bldr->stacks[ilvl].mask == 0) continue; - - /* Retrieve the node index of the next level */ - inode = (bldr->mcode >> (3*(ilvl+2))) & 7; - - /* Setup the parent node of the nodes registered into the current stack */ - parent_node = &bldr->stacks[ilvl+1].nodes[inode]; /* Fetch the parent node */ - stack_setup_node(&bldr->stacks[ilvl], parent_node, bldr->desc); - bldr->stacks[ilvl+1].mask |= (uint8_t)BIT(inode); - - /* Write the stacked nodes of the current level */ - res = stack_write - (&bldr->stacks[ilvl], bldr->buffer, &parent_node->ichild_node, &nleaves); - if(res != RES_OK) goto error; - bldr->nleaves += nleaves; - } - - /* Write the root node */ - ilvl = bldr->octree_depth-1; /* Root level */ - res = stack_write(&bldr->stacks[ilvl], bldr->buffer, root_id, &nleaves); - if(res != RES_OK) goto error; - bldr->nleaves += nleaves; - -exit: - return res; -error: - goto exit; -} - -static void -scene_release(ref_T* ref) -{ - struct htvox_scene* scn; - struct htvox_device* dev; - ASSERT(ref); - scn = CONTAINER_OF(ref, struct htvox_scene, ref); - octree_buffer_release(&scn->buffer); - dev = scn->dev; - MEM_RM(dev->allocator, scn); - HTVOX(device_ref_put(dev)); -} - -/******************************************************************************* - * Exported functions - ******************************************************************************/ -res_T -htvox_scene_create - (struct htvox_device* dev, - const double lower[3], /* Lower bound of the scene */ - const double upper[3], /* Upper bound of the scene */ - const size_t nvoxels[3], /* # voxels along the 3 axis */ - const struct htvox_voxel_desc* desc, /* Descriptor of a voxel */ - struct htvox_scene** out_scn) -{ - struct htvox_scene* scn = NULL; - double vox_sz[3]; /* World space size of a voxel */ - struct octree_builder bldr; - struct voxel vox = VOXEL_NULL; - uint64_t mcode_max; - uint64_t mcode; - res_T res = RES_OK; - - if(!dev || !lower || !upper || !nvoxels - || !check_htvox_voxel_desc(desc) || !out_scn) { - res = RES_BAD_ARG; - goto error; - } - if(lower[0] >= upper[0] - || lower[1] >= upper[1] - || lower[2] >= upper[2]) { - log_err(dev, - "%s: the scene AABB is degenerated.\n" - "\tlower={%g, %g, %g}, upper={%g, %g, %g}.\n", - FUNC_NAME, SPLIT3(lower), SPLIT3(upper)); - res = RES_BAD_ARG; - goto error; - } - if(!nvoxels[0] || !nvoxels[1] || !nvoxels[2]) { - log_err(dev, - "%s: the number of voxels along each axis cannot be null.\n" - "\t#voxels XYZ = {%lu, %lu, %lu}.\n", - FUNC_NAME, - (unsigned long)nvoxels[0], - (unsigned long)nvoxels[1], - (unsigned long)nvoxels[2]); - res = RES_BAD_ARG; - goto error; - } - scn = MEM_CALLOC(dev->allocator, 1, sizeof(struct htvox_scene)); - if(!scn) { - res = RES_MEM_ERR; - goto error; - } - ref_init(&scn->ref); - HTVOX(device_ref_get(dev)); - octree_buffer_init(dev->allocator, desc->size, &scn->buffer); - scn->dev = dev; - - /* Compute the octree definition */ - scn->definition = MMAX(nvoxels[0], 2); - scn->definition = MMAX(nvoxels[1], scn->definition); - scn->definition = MMAX(nvoxels[2], scn->definition); - scn->definition = round_up_pow2(scn->definition); - - /* Intialize the octree builder */ - res = octree_builder_init(&bldr, scn->definition, desc, &scn->buffer); - if(res != RES_OK) goto error; - - vox.data = MEM_CALLOC(dev->allocator, 1, desc->size); - if(!vox.data) { - res = RES_MEM_ERR; - goto error; - } - - mcode_max = scn->definition * scn->definition * scn->definition; - FOR_EACH(mcode, 0, mcode_max) { - size_t xyz[3]; - uint32_t ui3[3]; - - morton_xyz_decode_u21(mcode, ui3); - - /* Out of bound voxels */ - if(ui3[0] >= nvoxels[0] - || ui3[1] >= nvoxels[1] - || ui3[2] >= nvoxels[2]) - continue; - - /* Retrieve the voxel data from the caller */ - xyz[0] = (size_t)ui3[0]; - xyz[1] = (size_t)ui3[1]; - xyz[2] = (size_t)ui3[2]; - desc->get(xyz, vox.data, desc->context); - vox.mcode = mcode; - - /* Register the voxel against the octree */ - res = octree_builder_add_voxel(&bldr, &vox); - if(res != RES_OK) goto error; - } - - res = octree_builder_finalize(&bldr, &scn->root); - if(res != RES_OK) goto error; - - scn->nleaves = bldr.nleaves; - -#ifndef NDEBUG - check_octree(&scn->buffer, scn->root); -#endif - /* Setup the scene AABB in world space */ - d3_set(scn->lower, lower); - d3_set(scn->upper, upper); - - /* Compute the voxel size in world space */ - vox_sz[0] = (upper[0] - lower[0]) / (double)nvoxels[0]; - vox_sz[1] = (upper[1] - lower[1]) / (double)nvoxels[1]; - vox_sz[2] = (upper[2] - lower[2]) / (double)nvoxels[2]; - - /* Compute the octree AABB in world space */ - d3_set(scn->oclow, lower); - scn->ocupp[0] = (double)scn->definition * vox_sz[0]; - scn->ocupp[1] = (double)scn->definition * vox_sz[1]; - scn->ocupp[2] = (double)scn->definition * vox_sz[2]; - -exit: - if(vox.data) MEM_RM(dev->allocator, vox.data); - if(out_scn) *out_scn = scn; - return res; -error: - if(scn) { - HTVOX(scene_ref_put(scn)); - scn = NULL; - } - goto exit; -} - -res_T -htvox_scene_ref_get(struct htvox_scene* scn) -{ - if(!scn) return RES_BAD_ARG; - ref_get(&scn->ref); - return RES_OK; -} - -res_T -htvox_scene_ref_put(struct htvox_scene* scn) -{ - if(!scn) return RES_BAD_ARG; - ref_put(&scn->ref, scene_release); - return RES_OK; -} - -res_T -htvox_scene_get_aabb - (const struct htvox_scene* scn, - double lower[3], - double upper[3]) -{ - if(!scn || !lower || !upper) return RES_BAD_ARG; - d3_set(lower, scn->lower); - d3_set(upper, scn->upper); - return RES_OK; -} - -res_T -htvox_scene_get_leaves_count(const struct htvox_scene* scn, size_t* nleaves) -{ - if(!scn || !nleaves) return RES_BAD_ARG; - *nleaves = scn->nleaves; - return RES_OK; -} - -res_T -htvox_scene_for_each_leaf - (struct htvox_scene* scn, - void (*func) - (const void* val, - const size_t ileaf, - const double low[3], - const double upp[3], - void* ctx), - void* ctx) -{ - struct stack_entry { - struct octree_index inode; - double low[3]; - double upp[3]; - } stack[OCTREE_DEPTH_MAX*8]; - int istack; - size_t ileaf = 0; - - if(!scn || !func) return RES_BAD_ARG; - - stack[0].inode = scn->root; - stack[0].low[0] = scn->oclow[0]; - stack[0].low[1] = scn->oclow[1]; - stack[0].low[2] = scn->oclow[2]; - stack[0].upp[0] = scn->ocupp[0]; - stack[0].upp[1] = scn->ocupp[1]; - stack[0].upp[2] = scn->ocupp[2]; - istack = 1; - - do { - const struct stack_entry entry = stack[--istack]; - struct octree_xnode* node; - int ichild; - - node = octree_buffer_get_node(&scn->buffer, entry.inode); - - FOR_EACH(ichild, 0, 8) { - const uint8_t ichild_flag = (uint8_t)BIT(ichild); - double half_sz[3]; /* Half size of the current node */ - double mid[3]; /* Middle point of the current node */ - double upp[3]; - double low[3]; - - if((node->is_valid & ichild_flag) == 0) continue; /* Empty node */ - - half_sz[0] = (entry.upp[0] - entry.low[0])*0.5; - half_sz[1] = (entry.upp[1] - entry.low[1])*0.5; - half_sz[2] = (entry.upp[2] - entry.low[2])*0.5; - mid[0] = entry.low[0] + half_sz[0]; - mid[1] = entry.low[1] + half_sz[1]; - mid[2] = entry.low[2] + half_sz[2]; - low[0] = ichild&4 ? mid[0] : entry.low[0]; - low[1] = ichild&2 ? mid[1] : entry.low[1]; - low[2] = ichild&1 ? mid[2] : entry.low[2]; - upp[0] = low[0] + half_sz[0]; - upp[1] = low[1] + half_sz[1]; - upp[2] = low[2] + half_sz[2]; - - if(node->is_leaf & ichild_flag) { - struct octree_index iattr; - const void* val; - - iattr = octree_buffer_get_child_attr_index - (&scn->buffer, entry.inode, ichild); - val = octree_buffer_get_attr(&scn->buffer, iattr); - - ASSERT(upp[0] <= scn->upper[0]); - ASSERT(upp[1] <= scn->upper[1]); - ASSERT(upp[2] <= scn->upper[2]); - ASSERT(low[0] >= scn->lower[0]); - ASSERT(low[1] >= scn->lower[1]); - ASSERT(low[2] >= scn->lower[2]); - - func(val, ileaf, low, upp, ctx); - ileaf++; - } else { - struct stack_entry* top = stack + istack; - - top->inode = octree_buffer_get_child_node_index - (&scn->buffer, entry.inode, ichild); - top->low[0] = low[0]; - top->low[1] = low[1]; - top->low[2] = low[2]; - top->upp[0] = upp[0]; - top->upp[1] = upp[1]; - top->upp[2] = upp[2]; - ++istack; - } - } - } while(istack); - - return RES_OK; -} - -res_T -htvox_scene_at - (struct htvox_scene* scn, - const double position[3], - struct htvox_voxel* voxel) -{ - struct octree_index inode; - double scale_exp2; - double low[3]; - double pos[3]; - - if(!scn || !position || !voxel) return RES_BAD_ARG; - - *voxel = HTVOX_VOXEL_NULL; - - if(position[0] > scn->upper[0] || position[0] < scn->lower[0] - || position[1] > scn->upper[1] || position[1] < scn->lower[1] - || position[2] > scn->upper[2] || position[2] < scn->lower[2]) { - return RES_OK; - } - - /* Transform the position in the normalized octree space, - * i.e. octree lies in [0, 1]^2 */ - pos[0] = (position[0] - scn->oclow[0]) / (scn->ocupp[0] - scn->oclow[0]); - pos[1] = (position[1] - scn->oclow[1]) / (scn->ocupp[1] - scn->oclow[1]); - pos[2] = (position[2] - scn->oclow[2]) / (scn->ocupp[2] - scn->oclow[2]); - - /* Initialized the lower left corner of the current node */ - low[0] = 0; - low[1] = 0; - low[2] = 0; - - scale_exp2 = 0.5; - inode = scn->root; - *voxel = HTVOX_VOXEL_NULL; - for(;;) { - struct octree_xnode* node = octree_buffer_get_node(&scn->buffer, inode); - int ichild; - uint8_t ichild_flag; - double mid[3]; - - /* Compute the middle point of the node */ - mid[0] = low[0] + scale_exp2; - mid[1] = low[1] + scale_exp2; - mid[2] = low[2] + scale_exp2; - - /* Define the child of the current node into which pos `lies' and compute - * its lower left corner */ - ichild = 0; - if(pos[0] > mid[0]) { ichild |= 4; low[0] = mid[0]; } - if(pos[1] > mid[1]) { ichild |= 2; low[1] = mid[1]; } - if(pos[2] > mid[2]) { ichild |= 1; low[2] = mid[2]; } - - ichild_flag = (uint8_t)BIT(ichild); - if((node->is_valid & ichild_flag) == 0) { /* Empty node */ - break; - } else if(node->is_leaf & ichild_flag) { /* Leaf node */ - inode = octree_buffer_get_child_attr_index(&scn->buffer, inode, ichild); - voxel->data = octree_buffer_get_attr(&scn->buffer, inode); - voxel->id = inode.ipage * scn->buffer.pagesize / scn->buffer.voxsize - + inode.inode; - break; - } else { /* Child node */ - inode = octree_buffer_get_child_node_index(&scn->buffer, inode, ichild); - scale_exp2 *= 0.5; - } - } - return RES_OK; -} - diff --git a/src/htvox_scene.h b/src/htvox_scene.h @@ -1,37 +0,0 @@ -/* Copyright (C) CNRS 2018 - * - * 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 HTVOX_SCENE_H -#define HTVOX_SCENE_H - -#include "htvox_octree_buffer.h" -#include <rsys/ref_count.h> - -struct htvox_scene { - size_t definition; /* Definition of the octree */ - - double lower[3], upper[3]; /* World space AABB of the scene */ - double oclow[3], ocupp[3]; /* World space AABB of the octree */ - - struct octree_buffer buffer; - struct octree_index root; /* Index toward the root node of the octree */ - size_t nleaves; /* #leaves */ - - struct htvox_device* dev; - ref_T ref; -}; - -#endif /* HTVOX_SCENE_H */ - diff --git a/src/test_htvox_octree.c b/src/test_htvox_octree.c @@ -0,0 +1,381 @@ +/* Copyright (C) CNRS 2018 + * + * 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/>. */ + +#include "htvox.h" +#include "htvox_c.h" /* For morton_xyz_encode_u21 */ +#include "test_htvox_utils.h" + +#include <rsys/double3.h> + +struct check_context { + double* lower; + double* upper; + size_t* nvoxels; +}; + +static int +no_merge(const void* voxels[], const size_t nvoxels, void* ctx) +{ + CHK(voxels != NULL); + CHK(nvoxels != 0); + CHK((intptr_t)ctx == 0xDECAFBAD); + return 0; /* Merge nothing */ +} + +static int +merge_level0(const void** voxels, const size_t nvoxels, void* ctx) +{ + double min_val = DBL_MAX; + double max_val =-DBL_MAX; + size_t i; + CHK(voxels != NULL); + CHK(nvoxels != 0); + CHK((intptr_t)ctx == 0xDECAFBAD); + + FOR_EACH(i, 0, nvoxels) { + const double* val = voxels[i]; + min_val = MMIN(min_val, *val); + max_val = MMAX(max_val, *val); + } + + return (max_val - min_val) < 8; +} + +static void +keep_max(void* dst, const void* voxels[], const size_t nvoxels, void* ctx) +{ + double* vox_dst = dst; + double max_val = -DBL_MAX; + size_t i; + + CHK(dst != NULL); + CHK(voxels != NULL); + CHK(nvoxels != 0); + CHK(ctx != NULL); + + FOR_EACH(i, 0, nvoxels) { + const double* val = voxels[i]; + max_val = MMAX(max_val, *val); + } + *vox_dst = max_val; +} + +static void +get(const size_t xyz[3], void* dst, void* ctx) +{ + uint32_t ui3[3]; + uint64_t mcode; + double* val = dst; + CHK(xyz != NULL); + CHK(val != NULL); + CHK((intptr_t)ctx == 0xDECAFBAD); + + ui3[0] = (uint32_t)xyz[0]; + ui3[1] = (uint32_t)xyz[1]; + ui3[2] = (uint32_t)xyz[2]; + + mcode = morton_xyz_encode_u21(ui3); + *val = (double)mcode; +} + +static void +check_voxel + (const void* val, + const size_t ivoxel, + const double low[3], + const double upp[3], + void* context) +{ + const double* dbl = val; + struct check_context* ctx = context; + uint64_t mcode; + uint32_t xyz[3]; + double lower[3]; + double delta[3]; + + CHK(val != NULL); + CHK(low != NULL); + CHK(upp != NULL); + CHK(ctx != NULL); + CHK(low[0] < upp[0]); + CHK(low[1] < upp[1]); + CHK(low[2] < upp[2]); + CHK(ivoxel < ctx->nvoxels[0]*ctx->nvoxels[1]*ctx->nvoxels[2]); + CHK(*dbl >= 0); + + mcode = (uint64_t)(*dbl); + CHK(*dbl == (double)mcode); + + delta[0] = (ctx->upper[0] - ctx->lower[0]) / (double)ctx->nvoxels[0]; + delta[1] = (ctx->upper[1] - ctx->lower[1]) / (double)ctx->nvoxels[1]; + delta[2] = (ctx->upper[2] - ctx->lower[2]) / (double)ctx->nvoxels[2]; + + morton_xyz_decode_u21(mcode, xyz); + lower[0] = xyz[0] * delta[0]; + lower[1] = xyz[1] * delta[1]; + lower[2] = xyz[2] * delta[2]; + + CHK(eq_eps(lower[0], low[0], 1.e-6)); + CHK(eq_eps(lower[1], low[1], 1.e-6)); + CHK(eq_eps(lower[2], low[2], 1.e-6)); + CHK(eq_eps(lower[0] + delta[0], upp[0], 1.e-6)); + CHK(eq_eps(lower[1] + delta[1], upp[1], 1.e-6)); + CHK(eq_eps(lower[2] + delta[2], upp[2], 1.e-6)); +} + +static void +write_points + (const void* val, + const size_t ivoxel, + const double low[3], + const double upp[3], + void* context) +{ + FILE* stream = context; + (void)val, (void)ivoxel; + CHK(stream != NULL); + CHK(low != NULL); + CHK(upp != NULL); + fprintf(stream, + "%g %g %g\n%g %g %g\n%g %g %g\n%g %g %g\n" + "%g %g %g\n%g %g %g\n%g %g %g\n%g %g %g\n", + low[0], low[1], low[2], + upp[0], low[1], low[2], + low[0], upp[1], low[2], + upp[0], upp[1], low[2], + low[0], low[1], upp[2], + upp[0], low[1], upp[2], + low[0], upp[1], upp[2], + upp[0], upp[1], upp[2]); +} + +static void +write_cells + (const void* val, + const size_t ivoxel, + const double low[3], + const double upp[3], + void* context) +{ + FILE* stream = context; + (void)ivoxel, (void)val; + CHK(stream != NULL); + CHK(low != NULL); + CHK(upp != NULL); + fprintf(stream, "8 %lu %lu %lu %lu %lu %lu %lu %lu\n", + (unsigned long)(ivoxel*8 + 0), + (unsigned long)(ivoxel*8 + 1), + (unsigned long)(ivoxel*8 + 2), + (unsigned long)(ivoxel*8 + 3), + (unsigned long)(ivoxel*8 + 4), + (unsigned long)(ivoxel*8 + 5), + (unsigned long)(ivoxel*8 + 6), + (unsigned long)(ivoxel*8 + 7)); +} + +static void +write_scalars + (const void* val, + const size_t ivoxel, + const double low[3], + const double upp[3], + void* context) +{ + FILE* stream = context; + (void)ivoxel; + CHK(stream != NULL); + CHK(low != NULL); + CHK(upp != NULL); + fprintf(stream, "%g\n", *(double*)val); +} + +static void +dump_data(FILE* stream, struct htvox_octree* oct) +{ + struct htvox_octree_desc desc; + size_t ileaf; + + CHK(stream != NULL); + CHK(oct != NULL); + + fprintf(stream, "# vtk DataFile Version 2.0\n"); + fprintf(stream, "Volume\n"); + fprintf(stream, "ASCII\n"); + fprintf(stream, "DATASET UNSTRUCTURED_GRID\n"); + + CHK(htvox_octree_get_desc(oct, &desc) == RES_OK); + fprintf(stream, "POINTS %lu float\n", (unsigned long)(desc.nleaves* 8)); + CHK(htvox_octree_for_each_leaf(oct, write_points, stream) == RES_OK); + + fprintf(stream, "CELLS %lu %lu\n", + (unsigned long)desc.nleaves, + (unsigned long)(desc.nleaves*(8/*#verts per leaf*/ + 1/*1st field of a cell*/))); + CHK(htvox_octree_for_each_leaf(oct, write_cells, stream) == RES_OK); + + fprintf(stream, "CELL_TYPES %lu\n", (unsigned long)desc.nleaves); + FOR_EACH(ileaf, 0, desc.nleaves) fprintf(stream, "11\n"); + + fprintf(stream, "CELL_DATA %lu\n", (unsigned long)desc.nleaves); + fprintf(stream, "SCALARS K float 1\n"); + fprintf(stream, "LOOKUP_TABLE default\n"); + CHK(htvox_octree_for_each_leaf(oct, write_scalars, stream) == RES_OK); +} + +int +main(int argc, char** argv) +{ + struct htvox_device* dev = NULL; + struct htvox_octree* oct = NULL; + struct mem_allocator allocator; + struct htvox_voxel_desc voxdesc = HTVOX_VOXEL_DESC_NULL; + struct htvox_octree_desc octdesc = HTVOX_OCTREE_DESC_NULL; + double low[3]; + double upp[3]; + double pos[3]; + double delta[3]; + size_t nvxls[3]; + struct check_context ctx; + void* ptr = (void*)0xDECAFBAD; + size_t x, y, z; + (void)argc, (void)argv; + + CHK(mem_init_proxy_allocator(&allocator, &mem_default_allocator) == RES_OK); + + CHK(htvox_device_create(NULL, &allocator, 1, &dev) == RES_OK); + + d3_splat(low, 0.0); + d3_splat(upp, 1.0); + nvxls[0] = nvxls[1] = nvxls[2] = 5; + + ctx.lower = low; + ctx.upper = upp; + ctx.nvoxels = nvxls; + + #define NEW_SCN htvox_octree_create + + voxdesc.get = get; + voxdesc.merge = keep_max; + voxdesc.challenge_merge = no_merge; + voxdesc.context = ptr; + voxdesc.size = sizeof(double); + CHK(NEW_SCN(dev, low, upp, nvxls, &voxdesc, &oct) == RES_OK); + + CHK(htvox_octree_ref_get(NULL) == RES_BAD_ARG); + CHK(htvox_octree_ref_get(oct) == RES_OK); + CHK(htvox_octree_ref_put(NULL) == RES_BAD_ARG); + CHK(htvox_octree_ref_put(oct) == RES_OK); + CHK(htvox_octree_ref_put(oct) == RES_OK); + + upp[0] = low[0]; + CHK(NEW_SCN(dev, low, upp, nvxls, &voxdesc, &oct) == RES_BAD_ARG); + upp[0] = 1.0; + + nvxls[2] = 0; + CHK(NEW_SCN(dev, low, upp, nvxls, &voxdesc, &oct) == RES_BAD_ARG); + nvxls[2] = nvxls[0]; + + CHK(NEW_SCN(NULL, low, upp, nvxls, &voxdesc, &oct) == RES_BAD_ARG); + CHK(NEW_SCN(dev, NULL, upp, nvxls, &voxdesc, &oct) == RES_BAD_ARG); + CHK(NEW_SCN(dev, low, NULL, nvxls, &voxdesc, &oct) == RES_BAD_ARG); + CHK(NEW_SCN(dev, low, upp, NULL, &voxdesc, &oct) == RES_BAD_ARG); + CHK(NEW_SCN(dev, low, upp, nvxls, &voxdesc, NULL) == RES_BAD_ARG); + + voxdesc.get = NULL; + CHK(NEW_SCN(dev, low, upp, nvxls, &voxdesc, &oct) == RES_BAD_ARG); + voxdesc.get = get; + + voxdesc.merge = NULL; + CHK(NEW_SCN(dev, low, upp, nvxls, &voxdesc, &oct) == RES_BAD_ARG); + voxdesc.merge = keep_max; + + voxdesc.challenge_merge = NULL; + CHK(NEW_SCN(dev, low, upp, nvxls, &voxdesc, &oct) == RES_BAD_ARG); + voxdesc.challenge_merge = no_merge; + + voxdesc.size = 0; + CHK(NEW_SCN(dev, low, upp, nvxls, &voxdesc, &oct) == RES_BAD_ARG); + voxdesc.size = sizeof(double); + + voxdesc.size = HTVOX_MAX_SIZEOF_VOXEL + 1; + CHK(NEW_SCN(dev, low, upp, nvxls, &voxdesc, &oct) == RES_BAD_ARG); + voxdesc.size = sizeof(double); + + CHK(NEW_SCN(dev, low, upp, nvxls, &voxdesc, &oct) == RES_OK); + + CHK(htvox_octree_for_each_leaf(oct, check_voxel, &ctx) == RES_OK); + + CHK(htvox_octree_get_desc(NULL, &octdesc) == RES_BAD_ARG); + CHK(htvox_octree_get_desc(oct, NULL) == RES_BAD_ARG); + CHK(htvox_octree_get_desc(oct, &octdesc) == RES_OK); + CHK(octdesc.nleaves == 5*5*5); + CHK(octdesc.nvoxels == octdesc.nleaves + 36/*#parents*/); + CHK(octdesc.depth == 4); + + d3_splat(low, 0); + d3_splat(upp, 1); + CHK(d3_eq(octdesc.lower, low)); + CHK(d3_eq(octdesc.upper, upp)); + + delta[0] = (upp[0] - low[0])/(double)nvxls[0]; + delta[1] = (upp[1] - low[1])/(double)nvxls[1]; + delta[2] = (upp[2] - low[2])/(double)nvxls[2]; + + FOR_EACH(x, 0, nvxls[0]) { + pos[0] = low[0] + ((double)x + 0.5) * delta[0]; + FOR_EACH(y, 0, nvxls[1]) { + pos[1] = low[1] + ((double)y + 0.5) * delta[1]; + FOR_EACH(z, 0, nvxls[2]) { + struct htvox_voxel vox; + uint32_t ui3[3]; + uint64_t mcode; + + pos[2] = low[2] + ((double)z + 0.5) * delta[2]; + + ui3[0] = (uint32_t)x; + ui3[1] = (uint32_t)y; + ui3[2] = (uint32_t)z; + + CHK(htvox_octree_at(oct, pos, &vox) == RES_OK); + CHK(!HTVOX_VOXEL_NONE(&vox)); + + mcode = morton_xyz_encode_u21(ui3); + CHK(*((double*)vox.data) == mcode); + } + } + } + + CHK(htvox_octree_ref_put(oct) == RES_OK); + + nvxls[0] = nvxls[1] = nvxls[2] = 32; + voxdesc.challenge_merge = merge_level0; + CHK(NEW_SCN(dev, low, upp, nvxls, &voxdesc, &oct) == RES_OK); + CHK(htvox_octree_get_desc(oct, &octdesc) == RES_OK); + CHK(octdesc.nleaves == nvxls[0]*nvxls[1]*nvxls[2] / 8); + CHK(octdesc.nvoxels == (octdesc.nleaves*8 - 1) / 7); + CHK(octdesc.depth == (size_t)log2i((int)(nvxls[0]/2))+1); + + dump_data(stdout, oct); + + #undef NEW_SCN + + CHK(htvox_device_ref_put(dev) == RES_OK); + CHK(htvox_octree_ref_put(oct) == RES_OK); + + check_memory_allocator(&allocator); + mem_shutdown_proxy_allocator(&allocator); + CHK(mem_allocated_size() == 0); + return 0; +} + diff --git a/src/test_htvox_scene.c b/src/test_htvox_scene.c @@ -1,381 +0,0 @@ -/* Copyright (C) CNRS 2018 - * - * 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/>. */ - -#include "htvox.h" -#include "htvox_c.h" /* For morton_xyz_encode_u21 */ -#include "test_htvox_utils.h" - -#include <rsys/double3.h> - -struct check_context { - double* lower; - double* upper; - size_t* nvoxels; -}; - -static int -no_merge(const void* voxels[], const size_t nvoxels, void* ctx) -{ - CHK(voxels != NULL); - CHK(nvoxels != 0); - CHK((intptr_t)ctx == 0xDECAFBAD); - return 0; /* Merge nothing */ -} - -static int -merge_level0(const void** voxels, const size_t nvoxels, void* ctx) -{ - double min_val = DBL_MAX; - double max_val =-DBL_MAX; - size_t i; - CHK(voxels != NULL); - CHK(nvoxels != 0); - CHK((intptr_t)ctx == 0xDECAFBAD); - - FOR_EACH(i, 0, nvoxels) { - const double* val = voxels[i]; - min_val = MMIN(min_val, *val); - max_val = MMAX(max_val, *val); - } - - return (max_val - min_val) < 8; -} - -static void -keep_max(void* dst, const void* voxels[], const size_t nvoxels, void* ctx) -{ - double* vox_dst = dst; - double max_val = -DBL_MAX; - size_t i; - - CHK(dst != NULL); - CHK(voxels != NULL); - CHK(nvoxels != 0); - CHK(ctx != NULL); - - FOR_EACH(i, 0, nvoxels) { - const double* val = voxels[i]; - max_val = MMAX(max_val, *val); - } - *vox_dst = max_val; -} - -static void -get(const size_t xyz[3], void* dst, void* ctx) -{ - uint32_t ui3[3]; - uint64_t mcode; - double* val = dst; - CHK(xyz != NULL); - CHK(val != NULL); - CHK((intptr_t)ctx == 0xDECAFBAD); - - ui3[0] = (uint32_t)xyz[0]; - ui3[1] = (uint32_t)xyz[1]; - ui3[2] = (uint32_t)xyz[2]; - - mcode = morton_xyz_encode_u21(ui3); - *val = (double)mcode; -} - -static void -check_voxel - (const void* val, - const size_t ivoxel, - const double low[3], - const double upp[3], - void* context) -{ - const double* dbl = val; - struct check_context* ctx = context; - uint64_t mcode; - uint32_t xyz[3]; - double lower[3]; - double delta[3]; - - CHK(val != NULL); - CHK(low != NULL); - CHK(upp != NULL); - CHK(ctx != NULL); - CHK(low[0] < upp[0]); - CHK(low[1] < upp[1]); - CHK(low[2] < upp[2]); - CHK(ivoxel < ctx->nvoxels[0]*ctx->nvoxels[1]*ctx->nvoxels[2]); - CHK(*dbl >= 0); - - mcode = (uint64_t)(*dbl); - CHK(*dbl == (double)mcode); - - delta[0] = (ctx->upper[0] - ctx->lower[0]) / (double)ctx->nvoxels[0]; - delta[1] = (ctx->upper[1] - ctx->lower[1]) / (double)ctx->nvoxels[1]; - delta[2] = (ctx->upper[2] - ctx->lower[2]) / (double)ctx->nvoxels[2]; - - morton_xyz_decode_u21(mcode, xyz); - lower[0] = xyz[0] * delta[0]; - lower[1] = xyz[1] * delta[1]; - lower[2] = xyz[2] * delta[2]; - - CHK(eq_eps(lower[0], low[0], 1.e-6)); - CHK(eq_eps(lower[1], low[1], 1.e-6)); - CHK(eq_eps(lower[2], low[2], 1.e-6)); - CHK(eq_eps(lower[0] + delta[0], upp[0], 1.e-6)); - CHK(eq_eps(lower[1] + delta[1], upp[1], 1.e-6)); - CHK(eq_eps(lower[2] + delta[2], upp[2], 1.e-6)); -} - -static void -write_points - (const void* val, - const size_t ivoxel, - const double low[3], - const double upp[3], - void* context) -{ - FILE* stream = context; - (void)val, (void)ivoxel; - CHK(stream != NULL); - CHK(low != NULL); - CHK(upp != NULL); - fprintf(stream, - "%g %g %g\n%g %g %g\n%g %g %g\n%g %g %g\n" - "%g %g %g\n%g %g %g\n%g %g %g\n%g %g %g\n", - low[0], low[1], low[2], - upp[0], low[1], low[2], - low[0], upp[1], low[2], - upp[0], upp[1], low[2], - low[0], low[1], upp[2], - upp[0], low[1], upp[2], - low[0], upp[1], upp[2], - upp[0], upp[1], upp[2]); -} - -static void -write_cells - (const void* val, - const size_t ivoxel, - const double low[3], - const double upp[3], - void* context) -{ - FILE* stream = context; - (void)ivoxel, (void)val; - CHK(stream != NULL); - CHK(low != NULL); - CHK(upp != NULL); - fprintf(stream, "8 %lu %lu %lu %lu %lu %lu %lu %lu\n", - (unsigned long)(ivoxel*8 + 0), - (unsigned long)(ivoxel*8 + 1), - (unsigned long)(ivoxel*8 + 2), - (unsigned long)(ivoxel*8 + 3), - (unsigned long)(ivoxel*8 + 4), - (unsigned long)(ivoxel*8 + 5), - (unsigned long)(ivoxel*8 + 6), - (unsigned long)(ivoxel*8 + 7)); -} - -static void -write_scalars - (const void* val, - const size_t ivoxel, - const double low[3], - const double upp[3], - void* context) -{ - FILE* stream = context; - (void)ivoxel; - CHK(stream != NULL); - CHK(low != NULL); - CHK(upp != NULL); - fprintf(stream, "%g\n", *(double*)val); -} - -static void -dump_data(FILE* stream, struct htvox_scene* scn) -{ - size_t ivxl; - size_t nvxls; - - CHK(stream != NULL); - CHK(scn != NULL); - - fprintf(stream, "# vtk DataFile Version 2.0\n"); - fprintf(stream, "Volume\n"); - fprintf(stream, "ASCII\n"); - fprintf(stream, "DATASET UNSTRUCTURED_GRID\n"); - - CHK(htvox_scene_get_leaves_count(scn, &nvxls) == RES_OK); - fprintf(stream, "POINTS %lu float\n", (unsigned long)(nvxls * 8)); - CHK(htvox_scene_for_each_leaf(scn, write_points, stream) == RES_OK); - - fprintf(stream, "CELLS %lu %lu\n", - (unsigned long)nvxls, - (unsigned long)(nvxls*(8/*#verts per voxel*/ + 1/*1st field of a cell*/))); - CHK(htvox_scene_for_each_leaf(scn, write_cells, stream) == RES_OK); - - fprintf(stream, "CELL_TYPES %lu\n", (unsigned long)nvxls); - FOR_EACH(ivxl, 0, nvxls) fprintf(stream, "11\n"); - - fprintf(stream, "CELL_DATA %lu\n", (unsigned long)nvxls); - fprintf(stream, "SCALARS K float 1\n"); - fprintf(stream, "LOOKUP_TABLE default\n"); - CHK(htvox_scene_for_each_leaf(scn, write_scalars, stream) == RES_OK); -} - -int -main(int argc, char** argv) -{ - struct htvox_device* dev = NULL; - struct htvox_scene* scn = NULL; - struct mem_allocator allocator; - struct htvox_voxel_desc desc = HTVOX_VOXEL_DESC_NULL; - double low[3]; - double upp[3]; - double pos[3]; - double delta[3]; - size_t nvxls[3]; - size_t nvoxels; - struct check_context ctx; - void* ptr = (void*)0xDECAFBAD; - size_t x, y, z; - (void)argc, (void)argv; - - CHK(mem_init_proxy_allocator(&allocator, &mem_default_allocator) == RES_OK); - - CHK(htvox_device_create(NULL, &allocator, 1, &dev) == RES_OK); - - d3_splat(low, 0.0); - d3_splat(upp, 1.0); - nvxls[0] = nvxls[1] = nvxls[2] = 5; - - ctx.lower = low; - ctx.upper = upp; - ctx.nvoxels = nvxls; - - #define NEW_SCN htvox_scene_create - - desc.get = get; - desc.merge = keep_max; - desc.challenge_merge = no_merge; - desc.context = ptr; - desc.size = sizeof(double); - CHK(NEW_SCN(dev, low, upp, nvxls, &desc, &scn) == RES_OK); - - CHK(htvox_scene_ref_get(NULL) == RES_BAD_ARG); - CHK(htvox_scene_ref_get(scn) == RES_OK); - CHK(htvox_scene_ref_put(NULL) == RES_BAD_ARG); - CHK(htvox_scene_ref_put(scn) == RES_OK); - CHK(htvox_scene_ref_put(scn) == RES_OK); - - upp[0] = low[0]; - CHK(NEW_SCN(dev, low, upp, nvxls, &desc, &scn) == RES_BAD_ARG); - upp[0] = 1.0; - - nvxls[2] = 0; - CHK(NEW_SCN(dev, low, upp, nvxls, &desc, &scn) == RES_BAD_ARG); - nvxls[2] = nvxls[0]; - - CHK(NEW_SCN(NULL, low, upp, nvxls, &desc, &scn) == RES_BAD_ARG); - CHK(NEW_SCN(dev, NULL, upp, nvxls, &desc, &scn) == RES_BAD_ARG); - CHK(NEW_SCN(dev, low, NULL, nvxls, &desc, &scn) == RES_BAD_ARG); - CHK(NEW_SCN(dev, low, upp, NULL, &desc, &scn) == RES_BAD_ARG); - CHK(NEW_SCN(dev, low, upp, nvxls, &desc, NULL) == RES_BAD_ARG); - - desc.get = NULL; - CHK(NEW_SCN(dev, low, upp, nvxls, &desc, &scn) == RES_BAD_ARG); - desc.get = get; - - desc.merge = NULL; - CHK(NEW_SCN(dev, low, upp, nvxls, &desc, &scn) == RES_BAD_ARG); - desc.merge = keep_max; - - desc.challenge_merge = NULL; - CHK(NEW_SCN(dev, low, upp, nvxls, &desc, &scn) == RES_BAD_ARG); - desc.challenge_merge = no_merge; - - desc.size = 0; - CHK(NEW_SCN(dev, low, upp, nvxls, &desc, &scn) == RES_BAD_ARG); - desc.size = sizeof(double); - - desc.size = HTVOX_MAX_SIZEOF_VOXEL + 1; - CHK(NEW_SCN(dev, low, upp, nvxls, &desc, &scn) == RES_BAD_ARG); - desc.size = sizeof(double); - - CHK(NEW_SCN(dev, low, upp, nvxls, &desc, &scn) == RES_OK); - - CHK(htvox_scene_for_each_leaf(scn, check_voxel, &ctx) == RES_OK); - - CHK(htvox_scene_get_leaves_count(NULL, &nvoxels) == RES_BAD_ARG); - CHK(htvox_scene_get_leaves_count(scn, NULL) == RES_BAD_ARG); - CHK(htvox_scene_get_leaves_count(scn, &nvoxels) == RES_OK); - CHK(nvoxels == nvxls[0]*nvxls[1]*nvxls[2]); - - delta[0] = (upp[0] - low[0])/(double)nvxls[0]; - delta[1] = (upp[1] - low[1])/(double)nvxls[1]; - delta[2] = (upp[2] - low[2])/(double)nvxls[2]; - - FOR_EACH(x, 0, nvxls[0]) { - pos[0] = low[0] + ((double)x + 0.5) * delta[0]; - FOR_EACH(y, 0, nvxls[1]) { - pos[1] = low[1] + ((double)y + 0.5) * delta[1]; - FOR_EACH(z, 0, nvxls[2]) { - struct htvox_voxel vox; - uint32_t ui3[3]; - uint64_t mcode; - - pos[2] = low[2] + ((double)z + 0.5) * delta[2]; - - ui3[0] = (uint32_t)x; - ui3[1] = (uint32_t)y; - ui3[2] = (uint32_t)z; - - CHK(htvox_scene_at(scn, pos, &vox) == RES_OK); - CHK(!HTVOX_VOXEL_NONE(&vox)); - - mcode = morton_xyz_encode_u21(ui3); - CHK(*((double*)vox.data) == mcode); - } - } - } - - d3_splat(low, DBL_MAX); - d3_splat(upp,-DBL_MAX); - CHK(htvox_scene_get_aabb(NULL, low, upp) == RES_BAD_ARG); - CHK(htvox_scene_get_aabb(scn, NULL, upp) == RES_BAD_ARG); - CHK(htvox_scene_get_aabb(scn, low, NULL) == RES_BAD_ARG); - CHK(htvox_scene_get_aabb(scn, low, upp) == RES_OK); - CHK(low[0] == 0 && low[1] == 0 && low[2] == 0); - CHK(upp[0] == 1 && upp[1] == 1 && upp[2] == 1); - - CHK(htvox_scene_ref_put(scn) == RES_OK); - - nvxls[0] = nvxls[1] = nvxls[2] = 32; - desc.challenge_merge = merge_level0; - CHK(NEW_SCN(dev, low, upp, nvxls, &desc, &scn) == RES_OK); - CHK(htvox_scene_get_leaves_count(scn, &nvoxels) == RES_OK); - CHK(nvoxels == nvxls[0]*nvxls[1]*nvxls[2] / 8); - - dump_data(stdout, scn); - - #undef NEW_SCN - - CHK(htvox_device_ref_put(dev) == RES_OK); - CHK(htvox_scene_ref_put(scn) == RES_OK); - - check_memory_allocator(&allocator); - mem_shutdown_proxy_allocator(&allocator); - CHK(mem_allocated_size() == 0); - return 0; -} -