commit 8d06f72773e62c5e259e5087b4ed0fb59bc94a35
parent 8718e005b3c2b149203a3b845607014ff437d3be
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Thu, 8 Feb 2018 14:39:18 +0100
First draft of the rt routine
Diffstat:
4 files changed, 242 insertions(+), 16 deletions(-)
diff --git a/src/htvox.h b/src/htvox.h
@@ -52,6 +52,11 @@ struct htvox_hit {
struct htvox_voxel voxel; /* Intersected voxel */
};
+#define HTVOX_HIT_NULL__ {DBL_MAX, HTVOX_VOXEL_NULL__}
+static const struct htvox_hit HTVOX_HIT_NULL = HTVOX_HIT_NULL__;
+
+#define HTVOX_HIT_NONE(Hit) ((Hit)->distance >= DBL_MAX)
+
/* Forward declaration of external data types */
struct logger;
struct mem_allocator;
@@ -129,6 +134,7 @@ htvox_scene_trace_ray
(struct htvox_scene* scn,
const double ray_origin[3],
const double ray_direction[3],
+ const double ray_range[2],
struct htvox_hit* hit);
HTVOX_API res_T
diff --git a/src/htvox_scene.c b/src/htvox_scene.c
@@ -442,20 +442,6 @@ htvox_scene_create
octree_buffer_init(dev->allocator, &scn->buffer);
scn->dev = dev;
- /* 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];
- vox_sz_max = MMAX(MMAX(vox_sz[0], vox_sz[1]), vox_sz[2]);
-
- /* Compute the scale factor of the voxels to transform them in cube */
- scn->vox_scale[0] = 1.0;
- scn->vox_scale[1] = 1.0;
- scn->vox_scale[2] = 1.0;
- if(vox_sz[0] != vox_sz_max) scn->vox_scale[0] = vox_sz_max / vox_sz[0];
- if(vox_sz[1] != vox_sz_max) scn->vox_scale[1] = vox_sz_max / vox_sz[1];
- if(vox_sz[2] != vox_sz_max) scn->vox_scale[2] = vox_sz_max / vox_sz[2];
-
/* Compute the octree definition */
scn->definition = MMAX(nvoxels[0], 1);
scn->definition = MMAX(nvoxels[1], scn->definition);
@@ -502,9 +488,16 @@ htvox_scene_create
#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];
diff --git a/src/htvox_scene.h b/src/htvox_scene.h
@@ -20,7 +20,6 @@
#include <rsys/ref_count.h>
struct htvox_scene {
- double vox_scale[3]; /* Scale factor of the octree */
size_t definition; /* Definition of the octree */
size_t nvoxels; /* Overall #voxels effectively stored in the octree */
diff --git a/src/htvox_scene_trace_ray.c b/src/htvox_scene_trace_ray.c
@@ -0,0 +1,228 @@
+/* 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_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 htvox_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 = HTVOX_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
+htvox_scene_trace_ray
+ (struct htvox_scene* scn,
+ const double org[3],
+ const double dir[3],
+ const double range[2],
+ struct htvox_hit* hit)
+{
+ struct ray ray;
+ res_T res = RES_OK;
+
+ if(!scn || !org || !range || !hit) {
+ res = RES_BAD_ARG;
+ goto error;
+ }
+
+ *hit = HTVOX_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;
+}
