star-vx

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

commit e3b3bd37057cdac1e232ad86ce94d2b7e8c6a3ad
parent b17ab1b376a237491fa1f342f396a9e767d57689
Author: Vincent Forest <vincent.forest@meso-star.com>
Date:   Thu,  8 Mar 2018 16:19:56 +0100

Write a first draft of the svx_octree_trace_ray function

Diffstat:
Asrc/svx_octree_trace_ray.c | 330+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Dsrc/svx_scene_trace_ray.c | 228-------------------------------------------------------------------------------
2 files changed, 330 insertions(+), 228 deletions(-)

diff --git a/src/svx_octree_trace_ray.c b/src/svx_octree_trace_ray.c @@ -0,0 +1,330 @@ +/* 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_scene.h" + +struct ray { + double org[3]; + double range[2]; + double ts[3]; /* 1 / -abs(dir) */ + uint32_t octant_mask; +}; + +/******************************************************************************* + * Helper functions + ******************************************************************************/ +static FINLINE uint32_t +ftoui(const float f) +{ + union { uint32_t ui; float f; } ucast; + ucast.f = f; + return ucast.ui; +} + +static FINLINE float +uitof(const uint32_t ui) +{ + union { uint32_t ui; float f; } ucast; + ucast.ui = ui; + return ucast.f; +} + +/* Return RES_BAD_OP if the ray cannot intersect the scene */ +static res_T +setup_ray + (const struct scene* scn, + const double org[3], + const double dir[3], + const double range[2] + struct ray* ray) +{ + double rcp_ocsz[3]; /* Reciprocal size of the World space octree AABB */ + ASSERT(scn); + + if(range[0] >= range[1]) return RES_BAD_OP; /* Disabled ray */ + + /* Ray paralelle to an axis and that does not intersect the scene AABB */ + if((!dir[0] && (org[0] < scn->lower[0] || org[0] > scn->upper[0])) + || (!dir[1] && (org[1] < scn->lower[1] || org[1] > scn->upper[1])) + || (!dir[2] && (org[2] < scn->lower[2] || org[2] > scn->upper[2]))) { + return RES_BAD_OP; + } + + /* Compute reciprocal size of the world space octree AABB */ + rcp_ocsz[0] = 1.0 / (scn->ocupp[0] - scn->oclow[0]); + rcp_ocsz[1] = 1.0 / (scn->ocupp[1] - scn->oclow[1]); + rcp_ocsz[2] = 1.0 / (scn->ocupp[2] - scn->oclow[2]); + + /* Transform the ray origin in the [1, 2] space */ + ray->org[0] = (org[0] - scn->oclow[0]) * rcp_ocsz[0] + 1.0; + ray->org[1] = (org[1] - scn->oclow[1]) * rcp_ocsz[1] + 1.0; + ray->org[2] = (org[2] - scn->oclow[2]) * rcp_ocsz[2] + 1.0; + + /* Transform the range in the normalized octree space */ + ray->range[0] = range[0] * rcp_ocsz[0]; + ray->range[1] = range[1] * rcp_ocsz[1]; + + /* Transform the direction in the normalized octree space */ + ray->dir[0] = dir[0] * rcp_ocsz[0]; + ray->dir[1] = dir[1] * rcp_ocsz[1]; + ray->dir[2] = dir[2] * rcp_ocsz[2]; + + /* Let a ray defines 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 / -abs(dir[0]); + ray->ts[1] = 1.0 / -abs(dir[1]); + ray->ts[2] = 1.0 / -abs(dir[2]); + + /* Mirror rays with position directions */ + ray->octant_mask = 0; + if(ray_dir[0] > 0) { ray->octant_mask ^= 4; ray->org[0] = 3.0 - ray->org[0]; } + if(ray_dir[1] > 0) { ray->octant_mask ^= 2; ray->org[0] = 3.0 - ray->org[1]; } + if(ray_dir[2] > 0) { ray->octant_mask ^= 1; ray->org[0] = 3.0 - ray->org[2]; } + + return RES_OK; +} + +static res_T +trace_ray(const struct svx_octree* oct, const struct* ray, struct svx_hit* hit) +{ + #define SCALE_MAX 23 /* #mantisse bits */ + struct stack_entry { + struct octree_index inode; + double t_max; + } stack[SCALE_MAX + 1/*Dummy entry use to avoid invalid read*/]; + struct octree_index inode; + double t_min, tmax; + float low[3], upp[3]; /* AABB of the current node in normalized space */ + float mid[3]; /* Middle point of the current node in normalized space */ + float corner[3]; + float scale_exp2; + uint32_t scale; /* stack index */ + uint32_t ichild; + ASSERT(oct && ray && hit); + + hit = SVX_HIT_NONE; /* Initialise the hit to "no intersection" */ + + f3_splat(low, 1); + f3_splat(upp, 2); + + /* Compute the min/max ray intersection with the octree AABB in normalized + * space. Note that in this space, the octree AABB is in [1, 2]^3 */ + t_min = (2 - ray->org[0]) * ray->ts[0]; + t_min = MMAX((2 - ray->org[1]) * ray->ts[1], t_min); + t_min = MMAX((2 - ray->org[2]) * ray->ts[2], t_min); + t_max = (1 - ray->org[0]) * ray->ts[0]; + t_max = MMIN((1 - ray->org[1]) * ray->ts[1], t_max); + t_max = MMIN((1 - ray->org[2]) * ray->ts[2], t_max); + t_min = MMAX(ray->range[0], t_min); + t_max = MMIN(ray->range[1], t_max); + if(t_min >= t_max) return RES_OK; /* No intersection */ + + /* Traverrsal initialisation */ + inode = oct->root; + scale_exp2 = 0.5f; + scale = SCALE_MAX - 1; + + /* Define the first child id and the position of its lower left corner in the + * normalized octree space, i.e. in [1, 2]^3 */ + ichild = 0; + corner[0] = 1.f; + corner[1] = 1.f; + corner[2] = 1.f; + if((1.5 - ray->org[0])*ray->ts[0] > t_min) { ichild ^= 4; corner[0] = 1.5f; } + if((1.5 - ray->org[1])*ray->ts[1] > t_min) { ichild ^= 2; corner[1] = 1.5f; } + if((1.5 - ray->org[2])*ray->ts[2] > t_min) { ichild ^= 1; corner[2] = 1.5f; } + + /* Octree traversal */ + scale_max = scale + 1; + while(scale < scale_max) { + const struct octree_xnode* node = octree_buffer_get_node(&oct->buffer, inode); + double t_corner[3]; + double t_max_corner; + double t_max_child; + uint32_t ichild_adjusted = ichild ^ ray->octant_mask; + uint32_t ichild_flag = (uint32_t)BIT(ichild_adjusted); + uint32_t istep; + + /* Compute the exit point of the ray in the current child node */ + t_corner[0] = (corner[0] - ray->org[0])*ray->ts[0]; + t_corner[1] = (corner[1] - ray->org[1])*ray->ts[1]; + t_corner[2] = (corner[2] - ray->org[2])*ray->ts[2]; + t_max_corner = MMIN(MMIN(t_corner[0], t_corner[1]), t_corner[2]); + + /* Traverse the current child */ + if((node->is_valid & ichild_flag) + && t_min <= (t_max_child = MMIN(t_max, t_max_corner))) { + double t_max_parent; + double t_center[3]; + float scale_exp2_child; + + t_max_parent = t_max; + t_max = t_max_child; + + if(node->is_leaf & ichild_flag) break; /* Leaf node */ + /* TODO invoke filter function if necessary */ + + scale_exp2_child = scale_exp2 * 0.5f; + + /* center = corner - scale_exp2_child => + * t_center = ts*(corner + scale_exp2_child - org) + * t_center = t_corner + ts*scale_exp2_child + * Anyway we perforrm the whole computation to avoid numerical issues */ + t_center[0] = (corner[0] - ray->org[0] + scale_exp2_child) * ray->ts[0]; + t_center[1] = (corner[1] - ray->org[1] + scale_exp2_child) * ray->ts[1]; + t_center[2] = (corner[2] - ray->org[2] + scale_exp2_child) * ray->ts[2]; + + /* Push the parent node */ + stack[scale].t_max = t_max_parent; + stack[scale].inode = inode; + + /* Define the id and the lower left corner of the first child */ + ichild = 0; + if(t_center[0] > t_min) { ichild ^= 4; corner[0] += scale_exp2_child; } + if(t_center[1] > t_min) { ichild ^= 2; corner[1] += scale_exp2_child; } + if(t_center[2] > t_min) { ichild ^= 1; corner[2] += scale_exp2_child; } + + --scale; + scale_exp2 = scale_exp2_child; + continue; + } + + /* Define the id and the lower left corner of the next child */ + istep = 0; + if(t_corner[0] <= t_max_corner) { istep ^= 4; corner[0] -= scale_exp2; } + if(t_corner[1] <= t_max_corner) { istep ^= 2; corner[1] -= scale_exp2; } + if(t_corner[2] <= t_max_corner) { istep ^= 1; corner[2] -= scale_exp2; } + ichild ^= istep; + + t_min = t_max_corner; /* Adjust the ray range */ + + if((ichild & istep) != 0) { /* The ray exits the child. Pop the stack. */ + uint32_t diff_bits = 0; + uint32_t shift[3]; + + /* The IEEE-754 encoding of a single precision floating point number `f' + * whose binary representation is `F' is defined as: + * f = (-1)^S * 2^(E-127) + * * (1 + For(i in [0..22]) { M & BIT(22-i) ? 2^i : 0 }) + * with S = F / 2^31; E = F / 2^23 and M = F & (2^23 - 1). + * + * We transformed the SVO in a normalized translated space of [1, 2]^3. + * As a consequence, the coordinates of the lower left `corner' of a node + * have always a null exponent (i.e. E = 127). In addition, note that for + * each dimension, the i^th bit of the mantissa M is set if the corner is + * greater or equal to the median split of the node at the (23-i)^th + * level. + * + * For instance, considering the mantissa M=0x480000 of a X coordinate of + * a node lower left corner. According to its binary encoding, i.e. + * `100 1000 0000 0000 0000 0000', one can assume that: + * - X >= 1 + 2^-1 + * - X < 1 + 2^-1 + 2^-2 + * - X < 1 + 2^-1 + 2^-3 + * - X >= 1 + 2^-1 + 2^-4 + * + * Note that we ensure that the traversal along each dimension is from 2 + * to 1, i.e. the ray direction are always negative. To define if the + * median split of a node is traversed by a ray, it is thus sufficient to + * track the update of its corresponding bit from 1 to 0. + * + * In the following code we use this property to find the highest level + * into the node hierarchy whose median split was traversed by the ray. + * This is particularly usefull since, thanks to this information, one + * can pop the traversal stack directly up to this level. This remove the + * recursive popping and thus drastically reduced the computation cost of + * the 'stack pop' procedure. */ + if(istep & 4) diff_bits |= ftoui(corner[0]) ^ ftoui(corner[0]+scale_exp2); + if(istep & 2) diff_bits |= ftoui(corner[1]) ^ ftoui(corner[1]+scale_exp2); + if(istep & 1) diff_bits |= ftoui(corner[2]) ^ ftoui(corner[2]+scale_exp2); + scale = (ftoui((float)diff_bits) >> 23) - 127; + scale_exp2 = uitof(((scale - SCALE_MAX) + 127) << 23); + + inode = stack[scale].inode; + t_max = stack[scale].t_max; + + /* Compute the lower corner of the popped node */ + shift[0] = ftoui(corner[0]) >> scale; + shift[1] = ftoui(corner[1]) >> scale; + shift[2] = ftoui(corner[2]) >> scale; + corner[0] = uitof(shift[0] << scale); + corner[1] = uitof(shift[1] << scale); + corner[2] = uitof(shift[2] << scale); + + /* Define the index of the popped node */ + ichild = ((shift[0] & 1) << 2) | ((shift[1] & 1) << 1) | (shift[2] & 1); + } + } + + if(scale < scale_max) { + const uint8_t ichild_adjusted = (uint8_t)(ichild ^ ray->octant_mask); + const struct octree_xnode* node = octree_buffer_get_node(&oct->buffer, inode); + const struct octree_index iattr; + + /* Retrieve the node attribs */ + iattr = octree_buffer_get_child_attr_index(&oct->buffer, &inode, ichild_adjusted); + + /* Setu the hit */ + hit->distance = t_min < ray->range[0] ? t_max : t_min; + hit->voxel.data = octree_buffer_get_attr(&oct->buffer, &iattr); + hit->voxel.depth = SCALE_MAX - scale; + hit->voxel.id = absolute_attr_index(&oct->buffer, iattr); + hit->voxel.is_leaf = (node->is_leaf & ichild_adjusted) != 0; + hit->voxel.lower[0] = corner[0]; + hit->voxel.lower[1] = corner[1]; + hit->voxel.lower[2] = corner[2]; + hit->voxel.upper[0] = corner[0] + scale_exp2; + hit->voxel.upper[1] = corner[1] + scale_exp2; + hit->voxel.upper[2] = corner[2] + scale_exp2; + } + return RES_OK; +} + +/******************************************************************************* + * Exported functions + ******************************************************************************/ +res_T +svx_octree_trace_ray + (struct svx_octree* oct, + const double org[3], + const double dir[3], + const double range[2], + svx_hit_filter_function_T filter, + void* context, + struct svx_hit* hit) +{ + struct ray ray; + res_T res = RES_OK; + + if(!scn || !org || !range || !hit) { + res = RES_BAD_ARG; + goto error; + } + + *hit = SVX_HIT_NULL; + + res = setup_ray(scn, org, dir, range, &ray); + if(res == RES_BAD_OP) { /* The ray cannot intersect the scene. */ + res = RES_OK; + goto exit; + } + +exit: + return res; +error: + goto exit; +} diff --git a/src/svx_scene_trace_ray.c b/src/svx_scene_trace_ray.c @@ -1,228 +0,0 @@ -/* 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_scene.h" - -struct ray { - float org[3]; - float range[2]; - float ts[3]; /* 1 / -abs(dir) */ - uint32_t octant_mask; -}; - -/******************************************************************************* - * Helper functions - ******************************************************************************/ -/* Return RES_BAD_OP if the ray cannot intersect the scene */ -static res_T -setup_ray - (const struct scene* scn, - const double org[3], - const double dir[3], - const double range[2] - struct ray* ray) -{ - double rcp_ocsz[3]; /* Reciprocal size of the World space octree AABB */ - ASSERT(scn); - - if(range[0] >= range[1]) return RES_BAD_OP; /* Disabled ray */ - - /* Ray paralelle to an axis and that does not intersect the scene AABB */ - if((!dir[0] && (org[0] < scn->lower[0] || org[0] > scn->upper[0])) - || (!dir[1] && (org[1] < scn->lower[1] || org[1] > scn->upper[1])) - || (!dir[2] && (org[2] < scn->lower[2] || org[2] > scn->upper[2]))) { - return RES_BAD_OP; - } - - /* Compute reciprocal size of the world space octree AABB */ - rcp_ocsz[0] = 1.0 / (scn->ocupp[0] - scn->oclow[0]); - rcp_ocsz[1] = 1.0 / (scn->ocupp[1] - scn->oclow[1]); - rcp_ocsz[2] = 1.0 / (scn->ocupp[2] - scn->oclow[2]); - - /* Transform the ray origin in the [1, 2] space */ - ray->org[0] = (float)((org[0] - scn->oclow[0]) * rcp_ocsz[0] + 1.f); - ray->org[1] = (float)((org[1] - scn->oclow[1]) * rcp_ocsz[1] + 1.f); - ray->org[2] = (float)((org[2] - scn->oclow[2]) * rcp_ocsz[2] + 1.f); - - /* Transform the range in the normalized octree space */ - ray->range[0] = (float)(range[0] * rcp_ocsz[0]); - ray->range[1] = (float)(range[1] * rcp_ocsz[1]); - - /* Transform the direction in the normalized octree space */ - ray->dir[0] = (float)(dir[0] * rcp_ocsz[0]); - ray->dir[1] = (float)(dir[1] * rcp_ocsz[1]); - ray->dir[2] = (float)(dir[2] * rcp_ocsz[2]); - - /* Let a ray defines 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] = (float)(1.0 / -fabs(dir[0])); - ray->ts[1] = (float)(1.0 / -fabs(dir[1])); - ray->ts[2] = (float)(1.0 / -fabs(dir[2])); - - /* Mirror rays with position directions */ - ray->octant_mask = 0; - if(ray_dir[0] > 0) { ray->octant_mask ^= 4; ray->org[0] = 3.f - ray->org[0]; } - if(ray_dir[1] > 0) { ray->octant_mask ^= 2; ray->org[0] = 3.f - ray->org[1]; } - if(ray_dir[2] > 0) { ray->octant_mask ^= 1; ray->org[0] = 3.f - ray->org[2]; } - - return RES_OK; -} - -static res_T -trace_ray(const struct scene* scn, const struct* ray, struct svx_hit* hit) -{ - #define SCALE_MAX 23 /* #mantisse bits */ - struct octree_index inode; - float scale_exp2; - float t_min, tmax; - float low[3], upp[3]; /* AABB of the current node in normalized space */ - float mid[3]; /* Middle point of the current node in normalized space */ - float corner[3]; - uint32_t scale; /* stack index */ - uint32_t ichild; - ASSERT(scn && ray && hit); - - hit = SVX_HIT_NONE; /* Initialise the hit to "no intersection" */ - - f3_splat(lower, 1); - f3_splat(upper, 2); - - /* Compute the min/max ray intersection with the octree AABB in normalized - * space. Note that in this space, the octree AABB is in [1, 2]^3 */ - t_min = (2 - ray->org[0]) * ray->ts[0]; - t_min = MMAX((2 - ray->org[1]) * ray->ts[1], t_min); - t_min = MMAX((2 - ray->org[2]) * ray->ts[2], t_min); - t_max = (1 - ray->org[0]) * ray->ts[0]; - t_max = MMIN((1 - ray->org[1]) * ray->ts[1], t_max); - t_max = MMIN((1 - ray->org[2]) * ray->ts[2], t_max); - t_min = MMAX(ray->range[0], t_min); - t_max = MMIN(ray->range[1], t_max); - if(t_min >= t_max) return RES_OK; /* No intersection */ - - /* Traverrsal initialisation */ - inode = scn->root; - scale_exp2 = 0.5f; - scale = SCALE_MAX - 1; - - /* Define the first child id and the position of its lower left corner in the - * normalized octree space, i.e. in [1, 2]^3 */ - ichild = 0; - corner[0] = 1.f; - corner[1] = 1.f; - corner[2] = 1.f; - if((1.5f - ray->org[0])*ray->ts[0] > t_min) { ichild ^= 4; corner[0] = 1.5f; } - if((1.5f - ray->org[1])*ray->ts[1] > t_min) { ichild ^= 2; corner[1] = 1.5f; } - if((1.5f - ray->org[2])*ray->ts[2] > t_min) { ichild ^= 1; corner[2] = 1.5f; } - - /* Octree traversal */ - scale_max = scale + 1; - while(scale < scale_max) { - const struct octree_xnode* node; - float t_corner[3]; - float t_max_corner; - float t_max_child; - uint32_t ichild_adjusted = ichild ^ ray->octant_mask; - uint32_t istep; - - inode = octree_buffer_get_child_index - (&scn->buffer, inode, (int)ichild_adjusted); - node = octree_buffer_get_node(&scn->buffer, inode); - - /* Compute the exit point of the ray in the current child node */ - t_corner[0] = (corner[0] - ray->org[0])*ray->ts[0]; - t_corner[1] = (corner[1] - ray->org[1])*ray->ts[1]; - t_corner[2] = (corner[2] - ray->org[2])*ray->ts[2]; - t_max_corner = MMIN(MMIN(t_corner[0], t_corner[1]), t_corner[2]); - - /* Traverse the current child */ - if(!OCTREE_XNODE_IS_EMPTY(node) - && t_min <= (t_max_child = MMIN(t_max, t_max_corner))) { - float scale_exp2_child; - float t_max_parent; - float t_center[3]; - - t_max_parent = t_max; - t_max = t_max_child; - - if(OCTREE_XNODE_IS_LEAF(node)) break; /* Leaf node */ - - scale_exp2_child = scale_exp2 * 0.5f; - - /* center = corner - scale_exp2_child => - * t_center = ts*(corner + scale_exp2_child - org) - * t_center = t_corner + ts*scale_exp2_child - * Anyway we perforrm the whole computation to avoid numerical issues */ - t_center[0] = (corner[0] - ray->org[0] + scale_exp2_child) * ray->ts[0]; - t_center[1] = (corner[1] - ray->org[1] + scale_exp2_child) * ray->ts[1]; - t_center[2] = (corner[2] - ray->org[2] + scale_exp2_child) * ray->ts[2]; - - /* Push the parent node */ - stack[scale].t_max = t_max_parent; - stack[scale].inode = inode; - - /* Define the id and the lower left corner of the first child */ - ichild = 0; - if(t_center[0] > t_min) { ichild ^= 4; corner[0] += scale_exp2_child; } - if(t_center[1] > t_min) { ichild ^= 2; corner[1] += scale_exp2_child; } - if(t_center[2] > t_min) { ichild ^= 1; corner[2] += scale_exp2_child; } - - --scale; - scale_exp2 = scale_exp2_child; - continue; - } - - /* Define the id and the lower left corner of the next child */ - istep = 0; - if(t_corner[0] <= t_max_corner) { istep ^= 4; corner[0] -= scale_exp2; } - if(t_corner[1] <= t_max_corner) { istep ^= 2; corner[1] -= scale_exp2; } - if(t_corner[2] <= t_max_corner) { istep ^= 1; corner[2] -= scale_exp2; } - ichild ^= istep; - } -} - -/******************************************************************************* - * Exported functions - ******************************************************************************/ -res_T -svx_scene_trace_ray - (struct svx_scene* scn, - const double org[3], - const double dir[3], - const double range[2], - struct svx_hit* hit) -{ - struct ray ray; - res_T res = RES_OK; - - if(!scn || !org || !range || !hit) { - res = RES_BAD_ARG; - goto error; - } - - *hit = SVX_HIT_NULL; - - res = setup_ray(scn, org, dir, range, &ray); - if(res == RES_BAD_OP) { /* The ray cannot intersect the scene. */ - res = RES_OK; - goto exit; - } - -exit: - return res; -error: - goto exit; -}