commit f79261fed9077859666f05c97d27bb64b9ce3894
parent 21b67b362f5bf7750699bdbcf4fc627f729c3b1f
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Wed, 9 Dec 2020 17:54:09 +0100
First draft of the aabb/tetrahedron intersection test
Diffstat:
3 files changed, 335 insertions(+), 9 deletions(-)
diff --git a/src/suvm.h b/src/suvm.h
@@ -151,11 +151,9 @@ suvm_volume_at
struct suvm_primitive* prim, /* Geometric primitive where `pos' lies */
double barycentric_coords[4]); /* `pos' into the primitive */
-/* Iterate over the volumetric primitives that intersect the submitted Axis
- * Aligned Bounding Box and invoke the `clbk' function onto them. A primitive
- * is included in the AABB if at least one of its vertices is _strictly_
- * included into the submitted boundaries; a vertex lying exactly on `low' or
- * `upp' are considered outside the AABB */
+/* Iterate over the volumetric primitives whose axis aligned bounding box
+ * intersects the submitted Axis Aligned Bounding Box and invoke the `clbk'
+ * function onto them. */
SUVM_API res_T
suvm_volume_intersect_aabb
(struct suvm_volume* volume,
diff --git a/src/suvm_volume.h b/src/suvm_volume.h
@@ -23,6 +23,22 @@
#include <embree3/rtcore.h>
+/*
+ * Tetrahedron geometric layout
+ *
+ * 3 (top vertex) facet 0 = {0,1,2}
+ * + facet 1 = {0,3,1}
+ * /|`. facet 2 = {1,3,2}
+ * / | `. facet 3 = {2,3,0}
+ * / | `.
+ * / | `.
+ * 0 +....|.........+ 2
+ * `. | _,-'
+ * `,|_,-'
+ * +
+ * 1
+ */
+
/* Forward declarations */
struct suvm_device;
struct mem_allocator;
@@ -68,7 +84,7 @@ struct ALIGN(16) node_leaf {
struct suvm_volume {
struct darray_u32 indices; /* List of uint32_t[4] */
- struct darray_float positions; /* List of double[3] */
+ struct darray_float positions; /* List of float[3] */
struct darray_float normals; /* List of per tetrahedron facet normals */
struct buffer prim_data; /* Per primitive data */
struct buffer vert_data; /* Per vertex data */
@@ -105,6 +121,22 @@ volume_get_vertices_count(const struct suvm_volume* vol)
return sz / 3;
}
+static FINLINE float*
+volume_get_primitive_vertex_position
+ (const struct suvm_volume* vol,
+ const size_t iprim,
+ const size_t ivert, /* In [0, 3] */
+ float pos[3])
+{
+ const uint32_t* ids;
+ ASSERT(vol && iprim < volume_get_primitives_count(vol) && pos);
+ ids = darray_u32_cdata_get(&vol->indices) + iprim*4/*#vertex per tetra*/;
+ pos[0] = darray_float_cdata_get(&vol->positions)[ids[ivert]*3 + 0];
+ pos[1] = darray_float_cdata_get(&vol->positions)[ids[ivert]*3 + 1];
+ pos[2] = darray_float_cdata_get(&vol->positions)[ids[ivert]*3 + 2];
+ return pos;
+}
+
static INLINE float*
volume_compute_tetrahedron_normal
(const struct suvm_volume* vol,
@@ -145,16 +177,16 @@ static FINLINE float*
volume_get_tetrahedron_normal
(const struct suvm_volume* vol,
const size_t itetra,
- const int ifacet /* In [0, 3] */,
+ const int ifacet, /* In [0, 3] */
float normal[3])
{
- const size_t id =
- (itetra*4/*#facets per tetra*/ + (size_t)ifacet)*3/*#coords per normal*/;
ASSERT(vol && itetra < volume_get_primitives_count(vol) && normal);
if(!darray_float_size_get(&vol->normals)) {
volume_compute_tetrahedron_normal(vol, itetra, ifacet, normal);
} else {
+ const size_t id =
+ (itetra*4/*#facets per tetra*/ + (size_t)ifacet)*3/*#coords per normal*/;
ASSERT(id + 3 <= darray_float_size_get(&vol->normals));
normal[0] = darray_float_cdata_get(&vol->normals)[id+0];
normal[1] = darray_float_cdata_get(&vol->normals)[id+1];
diff --git a/src/suvm_volume_intersect_aabb.c b/src/suvm_volume_intersect_aabb.c
@@ -17,9 +17,305 @@
#include "suvm_device.h"
#include "suvm_volume.h"
+#include <rsys/float2.h>
+#include <rsys/float3.h>
+
+enum { X, Y, Z, AXES_COUNT }; /* Syntactic sugar */
+
+struct tetrahedron {
+ float v[4][AXES_COUNT]; /* Tetrahedron vertices */
+ float N[4][AXES_COUNT]; /* Tetrahedron normals */
+ float D[4]; /* Slope parameter of the plane equation Ax + By + Cz + D = 0 */
+
+ /* 2D equation of the silhouette edges projected along the X, Y and Z axis of
+ * the form Ax + By + C = 0, with (A, B) the normals of the edge */
+ float Ep[AXES_COUNT][4/*#edges max*/][3];
+
+ /* Number of silhouette edges in the ZY, XZ and XY planes (3 or 4) */
+ int nEp[AXES_COUNT];
+
+ /* Tetrahedron axis aligned bounding box */
+ float lower[AXES_COUNT];
+ float upper[AXES_COUNT];
+};
+
/*******************************************************************************
* Helper functions
******************************************************************************/
+static INLINE void
+tetrahedron_setup
+ (const struct suvm_volume* vol,
+ const float low[3],
+ const float upp[3],
+ const size_t itetra,
+ struct tetrahedron* tetra)
+{
+#ifndef NDEBUG
+ float center[AXES_COUNT]; /* Center of the tetrahedron */
+#endif
+ float e[6][AXES_COUNT];
+ float (*v)[AXES_COUNT];
+ float (*N)[AXES_COUNT];
+ float (*Ep)[4][3];
+ int i;
+ ASSERT(vol && low && upp && tetra);
+
+ /* Fetch tetrahedron vertices */
+ volume_get_primitive_vertex_position(vol, itetra, 0, tetra->v[0]);
+ volume_get_primitive_vertex_position(vol, itetra, 1, tetra->v[1]);
+ volume_get_primitive_vertex_position(vol, itetra, 2, tetra->v[2]);
+ volume_get_primitive_vertex_position(vol, itetra, 4, tetra->v[3]);
+ v = tetra->v;
+
+#ifndef NDEBUG
+ /* Compute the center of the tetrahedron */
+ center[X] = (v[0][X] + v[1][X] + v[2][X] + v[3][X]) * 0.25f;
+ center[Y] = (v[0][Y] + v[1][Y] + v[2][Y] + v[3][Y]) * 0.25f;
+ center[Z] = (v[0][Z] + v[1][Z] + v[2][Z] + v[3][Z]) * 0.25f;
+#endif
+
+ /* Check argument consistency */
+ ASSERT(MMIN(MMIN(v[0][X], v[1][X]), MMIN(v[2][X], v[3][X])) == low[X]);
+ ASSERT(MMIN(MMIN(v[0][Y], v[1][Y]), MMIN(v[2][Y], v[3][Y])) == low[Y]);
+ ASSERT(MMIN(MMIN(v[0][Z], v[1][Z]), MMIN(v[2][Z], v[3][Z])) == low[Z]);
+ ASSERT(MMAX(MMAX(v[0][X], v[1][X]), MMAX(v[2][X], v[3][X])) == upp[X]);
+ ASSERT(MMAX(MMAX(v[0][Y], v[1][Y]), MMAX(v[2][Y], v[3][Y])) == upp[Y]);
+ ASSERT(MMAX(MMAX(v[0][Z], v[1][Z]), MMAX(v[2][Z], v[3][Z])) == upp[Z]);
+
+ /* Setup tetrahedron AABB */
+ f3_set(tetra->lower, low);
+ f3_set(tetra->upper, upp);
+
+ /* Fetch tetrahedron normals */
+ volume_get_tetrahedron_normal(vol, itetra, 0, tetra->N[0]);
+ volume_get_tetrahedron_normal(vol, itetra, 1, tetra->N[1]);
+ volume_get_tetrahedron_normal(vol, itetra, 2, tetra->N[2]);
+ volume_get_tetrahedron_normal(vol, itetra, 3, tetra->N[3]);
+ N = tetra->N;
+
+ /* Compute the slope of the planes */
+ tetra->D[0] = -f3_dot(tetra->N[0], v[0]);
+ tetra->D[1] = -f3_dot(tetra->N[1], v[1]);
+ tetra->D[2] = -f3_dot(tetra->N[2], v[2]);
+ tetra->D[3] = -f3_dot(tetra->N[3], v[3]);
+
+ /* Compute tetrahedron edges */
+ f3_sub(e[0], v[1], v[0]);
+ f3_sub(e[1], v[2], v[1]);
+ f3_sub(e[2], v[0], v[2]);
+ f3_sub(e[3], v[0], v[3]);
+ f3_sub(e[4], v[1], v[3]);
+ f3_sub(e[5], v[2], v[3]);
+
+ /* Detect the silhouette edges of the tetrahedron once projected in the YZ,
+ * ZX and XY planes, and compute their 2D equation in the projected plane. The
+ * variable 'i' defines the projection axis which is successively X (YZ
+ * plane), Y (ZX plane) and Z (XY plane). The axes of the corresponding 2D
+ * repairs are defined by the 'j' and 'k' variables. */
+ Ep = tetra->Ep;
+ FOR_EACH(i, 0, AXES_COUNT) {
+ /* On 1st iteration 'j' and 'k' define the YZ plane
+ * On 2nd iteration 'j' and 'k' define the ZX plane
+ * On 3rd ietration 'j' and 'k' define the XY plane */
+ const int j = (i + 1) % AXES_COUNT;
+ const int k = (j + 1) % AXES_COUNT;
+
+ /* Register the number of detected silhouette edges */
+ int n = 0;
+
+ /* To detect the silhouette edges, check the sign of the normals of two
+ * adjacent facets for the coordinate of the projection axis. If the signs
+ * are the same, the facets look at the same direction regarding the
+ * projection axis. The other case means that one facet points toward the
+ * projection axis whole the other one points to the opposite direction:
+ * this is the definition of a silhouette edge.
+ *
+ * Once detected, we compute the equation of the silhouette edge:
+ *
+ * Ax + By + C = 0
+ *
+ * with A and B the coordinates of the edge normals and C the slope of the
+ * edge. Computing the normal is actually as simple as swizzling the
+ * coordinates of the edges projected in 2D and sign tuning. Let the
+ * projection axis X (i==X) and thus the YZ projection plane (j==Y &&
+ * k==Z). Assuming that the first edge 'e[0]' is a silhouette edge, ie the
+ * edge between the facet0 and the facet1 (refers to the suvm_volume.h for
+ * the geometric layout of a tetrahedron) once projected in the YZ plane
+ * the edge is {0, e[0][Y], e[0][Z]}. Its normal N can is defined as the
+ * cross product between this projected edge ands the projection axis:
+ *
+ * | 0 | | 1 | | 0 |
+ * N = | e[0][Y] | X | 0 | = | e[0][Z] |
+ * | e[0][Z] | | 0 | |-e[0][Y] |
+ *
+ * In the projection plane, the _unormalized_ normal of the {e[0][Y],
+ * e[0][Z]} edge is thus {e[0][Z], -e[0][Y]}. More generally, the normal of
+ * an edge 'I' {e[I][j], e[I][k]} projected along the 'i' axis in the 'jk'
+ * plane is {e[I][k],-e[I]j]}. Anyway, one has to revert the normal
+ * orientation to ensure a specific convention wrt the tetrahedron
+ * footprint. In our case we ensure that the normal points toward the
+ * footprint. Finalle the edge slope 'C' can be computed as:
+ * | A |
+ * C =- | B | . | Pj, Pk |
+ *
+ * with {Pj, Pk} a point belonging to the edge as for instance one of its
+ * vertices. */
+ if(signf(N[0][i]) != signf(N[1][i])) { /* The edge 0 is silhouette */
+ f2_normalize(Ep[i][n], f2(Ep[i][n], e[0][k], -e[0][j]));
+ if(N[0][i] > 0) f2_minus(Ep[i][n], Ep[i][n]);
+ Ep[i][n][2] = -f2_dot(Ep[i][n], v[0]);
+ ++n;
+ }
+ if(signf(N[0][i]) != signf(N[2][i])) { /* The edge 1 is silhouette */
+ f2_normalize(Ep[i][n], f2(Ep[i][n], e[1][k], -e[1][j]));
+ if(N[0][i] > 0) f2_minus(Ep[i][n], Ep[i][n]);
+ Ep[i][n][2] = -f2_dot(Ep[i][n], v[1]);
+ ++n;
+ }
+ if(signf(N[0][i]) != signf(N[3][i])) { /* The edge 2 is silhouette */
+ f2_normalize(Ep[i][n], f2(Ep[i][n], e[2][k], -e[2][j]));
+ if(N[0][i] > 0) f2_minus(Ep[i][n], Ep[i][n]);
+ Ep[i][n][2] = -f2_dot(Ep[i][n], v[2]);
+ ++n;
+ }
+ if(signf(N[1][i]) != signf(N[3][i])) { /* The edge 3 is silhouette */
+ f2_normalize(Ep[i][n], f2(Ep[i][n], e[3][k], -e[3][j]));
+ if(N[1][i] > 0) f2_minus(Ep[i][n], Ep[i][n]);
+ Ep[i][n][2] = -f2_dot(Ep[i][n], v[0]);
+ ++n;
+ }
+ if(signf(N[1][i]) != signf(N[2][i])) { /* The edge 4 is silhouette */
+ f2_normalize(Ep[X][n], f2(Ep[i][n], e[4][k], -e[4][j]));
+ if(N[2][i] > 0) f2_minus(Ep[i][n], Ep[i][n]);
+ Ep[i][n][2] = -f2_dot(Ep[i][n], v[1]);
+ ++n;
+ }
+ if(signf(N[2][i]) != signf(N[3][i])) { /* The edge 5 is silhouette */
+ f2_normalize(Ep[i][n], f2(Ep[i][n], e[5][k], -e[5][j]));
+ if(N[3][i] > 0) f2_minus(Ep[i][n], Ep[i][n]);
+ Ep[i][n][2] = -f2_dot(Ep[i][n], v[2]);
+ ++n;
+ }
+ /* A tetrahedron can have 3 or 4 silhouette edges once projected in 2D */
+ ASSERT(n == 3 || n == 4);
+ tetra->nEp[i] = n; /* Register the #silouhette edges for this project axis */
+
+#ifndef NDEBUG
+ /* Assert that the normals points toward the expected half space */
+ {
+ int iedge;
+ FOR_EACH(iedge, 0, n) {
+ /* Evaluate the edge equation regarding the projected polyhedron center
+ * and check that its sign is positive */
+ const float dst =
+ Ep[i][iedge][0] * center[j]
+ + Ep[i][iedge][1] * center[k]
+ + Ep[i][iedge][2];
+ ASSERT(dst > 0);
+ }
+#endif
+ }
+}
+
+/* Define if an axis aligned bounding box and a tetrahedron are intersected.
+ * Its implementation follows the algorithms proposed by N. Green in "Detecting
+ * intersection of a rectangular solid and a convex polyhedron" (Graphics Gems
+ * IV, 1994, p74--82) and by H. Ratschek and J. Rokne in "Test for intersection
+ * between box and tetrahedron" (Internation Journal of Computer Mathematics,
+ * 1997, p191--204) */
+static INLINE int
+aabb_intersect_tetrahedron
+ (const float low[3],
+ const float upp[3],
+ const struct tetrahedron* tetra)
+{
+ float intersect_aabb_low[3];
+ float intersect_aabb_upp[3];
+ int i;
+ int nplanes_including_aabb;
+ ASSERT(low && upp && tetra);
+ ASSERT(low[X] < upp[X]);
+ ASSERT(low[Y] < upp[Y]);
+ ASSERT(low[Z] < upp[Z]);
+
+ /* Check the intersection between the AABB and the tetra AABB */
+ FOR_EACH(i, 0, AXES_COUNT) {
+ intersect_aabb_low[i] = MMAX(low[i], tetra->lower[i]);
+ intersect_aabb_upp[i] = MMIN(upp[i], tetra->upper[i]);
+ if(intersect_aabb_low[i] > intersect_aabb_upp[i]) /* Do not intersect */
+ return 0;
+ }
+
+ /* Check if the tetrahedron is included into the aabb */
+ if(intersect_aabb_low[X] == low[X]
+ && intersect_aabb_low[Y] == low[Y]
+ && intersect_aabb_low[Z] == low[Z]) {
+ return 1;
+ }
+
+ /* #planes that totally includes the aabb on its positive side */
+ nplanes_including_aabb = 0;
+
+ /* Check the tetrahedron planes against the aabb */
+ FOR_EACH(i, 0, 4) {
+ float p[3], n[3];
+
+ /* Define the aabb vertices that is the farthest in the positive/negative
+ * direction of the face normal. */
+ f3_set(p, upp);
+ f3_set(n, low);
+ if(tetra->N[i][X] < 0) SWAP(float, p[X], n[X]);
+ if(tetra->N[i][Y] < 0) SWAP(float, p[Y], n[Y]);
+ if(tetra->N[i][Z] < 0) SWAP(float, p[Z], n[Z]);
+
+ /* Check that the box is totally outside the plane, To do this, check that
+ * 'p', aka the farthest aabb vertex in the positive direction of the plane
+ * normal, is not in the positive half space. In this case, the aabb is
+ * entirely outside the tetrahedron */
+ if((f3_dot(tetra->N[i], p) + tetra->D[i]) < 0) return 0;
+
+ /* Check if the box is totally inside the given plane. To do this, check
+ * that 'n', aka the farthest aabb vertex in the negative direction of the
+ * plane normal, is inside the positive half space */
+ if((f3_dot(tetra->N[i], n) + tetra->D[i]) > 0) {
+ /* Register this plane as a plane including the aabb */
+ nplanes_including_aabb += 1;
+ }
+ }
+
+ /* Check if the aabb is entirely included into the tetrahedron */
+ if(nplanes_including_aabb == 4) return 1;
+
+ /* For each silhouette edge in each projection plane, check if it totally
+ * exclude the projected aabb */
+ FOR_EACH(i, 0, AXES_COUNT) {
+ /* Compute the remaining axes.
+ * On 1st iteration 'j' and 'k' define the YZ plane
+ * On 2nd iteration 'j' and 'k' define the ZX plane
+ * On 3rd ietration 'j' and 'k' define the XY plane */
+ const int j = (i + 1) % AXES_COUNT;
+ const int k = (j + 1) % AXES_COUNT;
+
+ int iedge;
+ FOR_EACH(iedge, 0, tetra->nEp[i]) {
+ float p[2];
+
+ /* Define the aabb vertices that is the farthest in the positive/negative
+ * direction of the normal of the silhouette edge. */
+ p[0] = tetra->Ep[i][iedge][0] > 0 ? upp[j] : low[j];
+ p[1] = tetra->Ep[i][iedge][1] > 0 ? upp[k] : low[k];
+
+ /* Check that the projected aabb along the 'i'th axis is totally outside
+ * the current silhouette edge. That means that the aabb is entirely
+ * outside the tetrahedron and thus that they do not intersect */
+ if((f2_dot(tetra->Ep[i][iedge], p) + tetra->Ep[i][iedge][2]) < 0) return 0;
+ }
+ }
+
+ /* The tetrahedron and the aabb are partially intersecting */
+ return 1;
+}
+
static INLINE int
aabb_intersect_aabb
(const float low0[3],
