commit 8a88c09b55fdafafb0efcc4495408c5889b7a782
parent b8f4e5bf065350d3e6f01383f1941bbd66132da9
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Thu, 13 Sep 2018 10:41:21 +0200
First draft of the binary tree intersector
Diffstat:
4 files changed, 257 insertions(+), 7 deletions(-)
diff --git a/cmake/CMakeLists.txt b/cmake/CMakeLists.txt
@@ -42,6 +42,7 @@ set(VERSION ${VERSION_MAJOR}.${VERSION_MINOR}.${VERSION_PATCH})
set(SVX_FILES_SRC
svx_bintree.c
+ svx_bintree_trace_ray.c
svx_buffer.c
svx_device.c
svx_octree.c
@@ -69,7 +70,7 @@ if(CMAKE_COMPILER_IS_GNUCC)
endif()
add_library(svx SHARED ${SVX_FILES_SRC} ${SVX_FILES_INC} ${SVX_FILES_INC_API})
-target_link_libraries(svx RSys)
+target_link_libraries(svx RSys ${MATH_LIB})
set_target_properties(svx PROPERTIES
DEFINE_SYMBOL SVX_SHARED_BUILD
diff --git a/src/svx.h b/src/svx.h
@@ -229,7 +229,7 @@ svx_bintree_create
(struct svx_device* dev,
const double lower, /* Lower bound of the bintree */
const double upper, /* Upper bound of the bintree */
- const size_t nvoxels, /* # voxels along the range */
+ const size_t nvoxels, /* #voxels along the range */
const enum svx_axis axis, /* Axis along which the binary tree is defined */
const struct svx_voxel_desc* desc, /* Descriptor of a voxel */
struct svx_tree** tree);
@@ -264,9 +264,21 @@ svx_octree_trace_ray
void* context, /* Data sent to the filter functor */
struct svx_hit* hit);
+/* TODO merge with the svx_octree_trace_ray API */
+SVX_API res_T
+svx_bintree_trace_ray
+ (struct svx_tree* tree,
+ const double ray_origin[3],
+ const double ray_direction[3],
+ const double ray_range[2],
+ const svx_hit_challenge_T challenge, /* NULL <=> Traversed up to the leaves */
+ const svx_hit_filter_T filter, /* NULL <=> Stop RT at the 1st hit voxel */
+ void* context, /* Data sent to the filter functor */
+ struct svx_hit* hit);
+
SVX_API res_T
svx_tree_at
- (struct svx_tree* octree,
+ (struct svx_tree* tree,
const double position[3],
svx_at_filter_T filter,
void* context, /* Client data sent as the last argument of the filter func */
diff --git a/src/svx_bintree_trace_ray.c b/src/svx_bintree_trace_ray.c
@@ -0,0 +1,236 @@
+/* Copyright (C) 2018 Université Paul Sabatier, |Meso|Star>
+ *
+ * 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 "svx.h"
+#include "svx_tree.h"
+
+#include <rsys/double3.h>
+
+/*******************************************************************************
+ * Helper function
+ ******************************************************************************/
+static FINLINE void
+setup_hit
+ (struct svx_tree* btree,
+ const struct buffer_index iattr, /* Index toward the voxel attributes */
+ const double distance_min, /* Dst from ray org to the voxel entry point */
+ const double distance_max, /* Dst from ray_org to the voxel exit point */
+ const double lower, /* Lower bound of the current voxel in [0, 1] space */
+ const double scale_exp2, /* Size of the voxel in [0, 1] space */
+ const size_t depth, /* Depth of the voxel in the octree hierarchy */
+ const int is_leaf, /* Define if the voxel is a leaf or not */
+ const int flip, /* Define if the btree was reverted or not */
+ struct svx_hit* hit)
+{
+ double low;
+ double upp;
+ ASSERT(btree && hit);
+ ASSERT(distance_min >= 0 && distance_min <= distance_max);
+
+ hit->distance[0] = distance_min;
+ hit->distance[1] = distance_max;
+ hit->voxel.data = buffer_get_attr(&btree->buffer, iattr);
+ hit->voxel.depth = depth,
+ hit->voxel.id = buffer_absolute_attr_index(&btree->buffer, iattr);
+ hit->voxel.is_leaf = is_leaf;
+
+ /* Transform the voxel aabb in the [0, 1]^3 normalized space and flip if
+ * necessary */
+ low = flip ? 1 - lower - scale_exp2 : lower;
+ upp = low + scale_exp2;
+
+ /* Transform the voxel AABB in world space */
+ d3_splat(hit->voxel.lower, -INF);
+ d3_splat(hit->voxel.upper, INF);
+ hit->voxel.lower[btree->frame[0]] = low;
+ hit->voxel.upper[btree->frame[1]] = upp;
+}
+
+/*******************************************************************************
+ * Exported function
+ ******************************************************************************/
+res_T
+svx_bintree_trace_ray
+ (struct svx_tree* btree,
+ const double ray_org[3],
+ const double ray_dir[3],
+ const double ray_range[2],
+ svx_hit_challenge_T challenge,
+ svx_hit_filter_T filter,
+ void* context,
+ struct svx_hit* hit)
+{
+ #define STACK_LENGTH 23 /* #mantissa bits */
+ ALIGN(64) struct stack_entry {
+ float scale_exp2;
+ uint32_t depth;
+ struct buffer_index inode;
+ } stack[STACK_LENGTH];
+ STATIC_ASSERT(sizeof(struct stack_entry) == 16, Unexpected_sizeof_stack_entry);
+
+ struct buffer_index inode;
+ size_t istack; /* Top of the stack */
+ size_t iaxis; /* Id in [0, 2] of the axis along which the tree is defined */
+ float pos_min, pos_max; /* Min/Max pos along the ray in [0,1] +dir 1D space */
+ float org; /* Ray origin in the [0,1] +dir 1D space */
+ float dir; /* Ray direction in the [0,1] +dir 1D space */
+ float ts; /* 1/(Ray direction) in the [0,1] +dir 1D space */
+ float rcp_btreesz; /* Reciprocal of the binary tree size */
+ float low; /* Lower bound of the traversed node */
+ float scale_exp2; /* Current size of a node in the [0,1] +dir 1D space */
+ uint32_t depth; /* Depth of the traversed nodes */
+ int ichild; /* Traversed child */
+ int flip = 0; /* Is the ray flipped? */
+ int null_dir = 0; /* Is the 1D dir is roughly null? */
+
+ if(!btree || !ray_org || !ray_dir || !ray_range || btree->type != SVX_BINTREE
+ || !d3_is_normalized(ray_dir)) {
+ return RES_BAD_ARG;
+ }
+
+ *hit = SVX_HIT_NULL;
+ if(ray_range[0] >= ray_range[1]) { /* Disabled ray */
+ return RES_OK;
+ }
+
+ iaxis = btree->frame[0];
+ rcp_btreesz = (float)(1.0 / btree->tree_size[0]);
+
+ /* Transform the ray origin in [0, 1] space */
+ org = (float)(ray_org[iaxis] - btree->tree_low[0]) * rcp_btreesz;
+ /* Transform the direction in the normalized bintree space */
+ dir = (float)(ray_dir[iaxis] * rcp_btreesz);
+ /* Define if the 1D ray direction is roughly null */
+ null_dir = eq_epsf(dir, 0, 1.e-6f);
+
+ /* The ray starts outside the binary tree and point outward the bin tree: it
+ * cannot intersect the binary tree */
+ if((org > 1 && (dir > 0 || null_dir))
+ || (org < 0 && (dir < 0 || null_dir)))
+ return RES_OK;
+
+ /* Mirror rays with negative direction */
+ if(dir < 0) {
+ flip = 1;
+ org = 1 - org;
+ dir = -dir;
+ }
+
+ /* Let a ray r defined as O + tD with O the ray origin and D the direction;
+ * and X an axis aligned plane. r intersects X at a distance t computed as
+ * below :
+ * t = (X-O) / D
+ * Note that one can transform r in r' with a transform M without any impact
+ * on the t parameter:
+ * M.X = M.O * t(M.D)
+ * t = (M.X-M.O)/(M.D) = (M.X-M.O)*ts; with ts = 1/(M.D) */
+ ts = null_dir ? (float)INF : 1.f / dir;
+
+ /* Compute the range in [0, 1] of the ray/binary tree intersection */
+ pos_min = (float)MMAX(dir * ray_range[0] + org, 0);
+ pos_max = (float)MMIN(dir * ray_range[1] + org, 1);
+
+ /* Defone the first traversed child and set its lower bound */
+ if(org < 0.5) {
+ ichild = 0;
+ low = 0.0;
+ } else {
+ ichild = 1;
+ low = 0.5;
+ }
+
+ /* Init the traversal */
+ depth = 1;
+ istack = 0;
+ scale_exp2 = 0.5;
+ inode = btree->root;
+
+ /* Here we go */
+ while(pos_min < pos_max) {
+ const struct buffer_xnode* node = buffer_get_node(&btree->buffer, inode);
+ const int ichild_adjusted = ichild ^ flip;
+ const int ichild_flag = BIT(ichild_adjusted);
+
+ if(node->is_valid & ichild_flag) {
+ const int is_leaf = (node->is_leaf & ichild_flag) != 0;
+ int go_deeper = 1;
+
+ if(is_leaf || challenge) {
+ struct svx_hit hit_tmp;
+ const float upp = low + scale_exp2;
+ const double t_min = null_dir ? ray_range[0] : (pos_min - org) * ts;
+ const double t_max = null_dir ? ray_range[1] : (upp - org) * ts;
+ const struct buffer_index iattr = buffer_get_child_attr_index
+ (&btree->buffer, inode, ichild_adjusted);
+ ASSERT(t_min < t_max);
+
+ setup_hit(btree, iattr, t_min, t_max, low, scale_exp2, depth, is_leaf,
+ flip, &hit_tmp);
+
+ if(is_leaf || challenge(&hit_tmp, context)) {
+ go_deeper = 0;
+ /* Stop the traversal if no filter is defined or if the filter
+ * function returns 0 */
+ if(!filter /* By default, i.e. with no filter, stop the traversal */
+ || !filter(&hit_tmp, ray_org, ray_dir, ray_range, context)) {
+ *hit = hit_tmp;
+ break;
+ }
+ }
+ }
+
+ if(go_deeper) {
+ const float scale_exp2_child = scale_exp2 * 0.5f;
+ const float child_mid = low + scale_exp2_child;
+
+ /* Push parent node if necessary */
+ if(ichild == 0 && !null_dir) {
+ stack[istack].inode = inode;
+ stack[istack].scale_exp2 = scale_exp2;
+ stack[istack].depth = depth;
+ ++istack;
+ }
+
+ if(pos_min > child_mid) low = child_mid;
+ scale_exp2 = scale_exp2_child;
+ ++depth;
+ }
+ }
+
+ /* Purse traversal */
+ ASSERT(pos_min >= low && pos_min < low + scale_exp2);
+ low += scale_exp2; /* Lower bound of the next node to traverse */
+ pos_min = low; /* Snap the pos_min to the next node lower bound */
+
+ if(ichild == 1) { /* No more child to traverse in this node => Pop node */
+ if(istack == 0) break; /* No more node to traverse */
+
+ /* Pop node */
+ --istack;
+ inode = stack[istack].inode;
+ scale_exp2 = stack[istack].scale_exp2;
+ depth = stack[istack].depth;
+
+ /* A node is pushed on the stack if one of its child was not traversed.
+ * Furthermore, this intersector ensures that the ray dir is alwas
+ * positive along the axis of the binary tree. As a consequence, the
+ * child 0 is always traversed first and the remaining child to traverse
+ * of a pushed node is thus always the child 1 */
+ ichild = 1;
+ }
+ }
+ return RES_OK;
+}
+
diff --git a/src/svx_octree_trace_ray.c b/src/svx_octree_trace_ray.c
@@ -61,7 +61,7 @@ setup_hit
const double distance_min, /* Dst from ray org to the voxel entry point */
const double distance_max, /* Dst from ray_org to the voxel exit point */
const float corner[3], /* Corner of the current voxel in [1, 2]^3 space */
- const double scale_exp2, /* Size of the voxel in [1, 2]^3] space */
+ const double scale_exp2, /* Size of the voxel in [1, 2]^3 space */
const size_t depth, /* Depth of the voxel in the octree hierarchy */
const int is_leaf, /* Define if the voxel is a leaf or not */
const uint32_t octant_mask, /* bitmask of reverted dimensions */
@@ -140,7 +140,7 @@ setup_ray
dir_adjusted[1] = dir[1] * rcp_ocsz[1];
dir_adjusted[2] = dir[2] * rcp_ocsz[2];
- /* Let a ray defines as org + t*dir and X the coordinate of an axis aligned
+ /* Let a ray defined as org + t*dir and X the coordinate of an axis aligned
* plane. The ray intersects X at t = (X - org)/dir = (X - org) * ts; with ts
* = 1/dir. Note that one assume that dir is always negative. */
ray->ts[0] = 1.0 / -fabs(dir_adjusted[0]);
@@ -373,7 +373,7 @@ trace_ray
}
/*******************************************************************************
- * Exported functions
+ * Exported function
******************************************************************************/
res_T
svx_octree_trace_ray
@@ -389,7 +389,7 @@ svx_octree_trace_ray
struct ray ray;
res_T res = RES_OK;
- if(!oct || !org || !range || !hit || oct->type != SVX_OCTREE) {
+ if(!oct || !org || !dir || !range || !hit || oct->type != SVX_OCTREE) {
res = RES_BAD_ARG;
goto error;
}
@@ -411,3 +411,4 @@ error:
goto exit;
}
+
