commit da5c391fd93d5b0ee2997ecac9d8b3abac1a8cd2
parent 54d6871c4061ad7ec9beab96b8221b9fb97d3c8c
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Fri, 25 Sep 2020 17:00:07 +0200
Implement the suvm_volume_at function
Diffstat:
5 files changed, 324 insertions(+), 48 deletions(-)
diff --git a/cmake/CMakeLists.txt b/cmake/CMakeLists.txt
@@ -47,7 +47,8 @@ set(VERSION ${VERSION_MAJOR}.${VERSION_MINOR}.${VERSION_PATCH})
set(SUVM_FILES_SRC
suvm_device.c
- suvm_volume.c)
+ suvm_volume.c
+ suvm_volume_at.c)
set(SUVM_FILES_INC
suvm_device.h
suvm_volume.h)
diff --git a/src/suvm.h b/src/suvm.h
@@ -125,7 +125,7 @@ suvm_volume_ref_put
SUVM_API res_T
suvm_volume_get_aabb
- (struct suvm_volume* volume,
+ (const struct suvm_volume* volume,
double lower[3],
double upper[3]);
diff --git a/src/suvm_volume.c b/src/suvm_volume.c
@@ -21,11 +21,8 @@
#include <rsys/cstr.h>
#include <rsys/dynamic_array_double.h>
#include <rsys/dynamic_array_size_t.h>
-#include <rsys/float3.h>
#include <rsys/ref_count.h>
-#include <embree3/rtcore_builder.h>
-
/* Generate the dynamic array of RTCBuildPrimitive */
#define DARRAY_NAME rtc_prim
#define DARRAY_DATA struct RTCBuildPrimitive
@@ -163,42 +160,6 @@ error:
goto exit;
}
-static INLINE float*
-compute_tetrahedron_normal
- (const struct suvm_volume* vol,
- const size_t itetra,
- const int ifacet/*In [0,3]*/,
- float normal[3])
-{
- const uint32_t* ids;
- const float* v0;
- const float* v1;
- const float* v2;
-
- float E0[3], E1[3];
-
- /* Facet indices, Ensure CCW ordering */
- const int facet[4][3] = {{0,1,2}, {0,3,1}, {1,3,2}, {2,3,0}};
-
- ASSERT(vol && itetra < volume_get_primitives_count(vol));
-
- /* Fetch tetrahedron indices */
- ids = darray_u32_cdata_get(&vol->indices) + itetra*4;
-
- /* Fetch facet vertices */
- v0 = darray_float_cdata_get(&vol->positions) + ids[facet[ifacet][0]]*3;
- v1 = darray_float_cdata_get(&vol->positions) + ids[facet[ifacet][1]]*3;
- v2 = darray_float_cdata_get(&vol->positions) + ids[facet[ifacet][2]]*3;
-
- /* Compute the normal */
- f3_sub(E0, v1, v0);
- f3_sub(E1, v2, v0);
- f3_cross(normal, E0, E1);
- f3_normalize(normal, normal);
-
- return normal;
-}
-
static res_T
compute_tetrahedral_mesh_normals(struct suvm_volume* vol)
{
@@ -218,10 +179,10 @@ compute_tetrahedral_mesh_normals(struct suvm_volume* vol)
float* n1 = darray_float_data_get(&vol->normals) + (itetra*4 + 1)*3;
float* n2 = darray_float_data_get(&vol->normals) + (itetra*4 + 3)*3;
float* n3 = darray_float_data_get(&vol->normals) + (itetra*4 + 2)*3;
- compute_tetrahedron_normal(vol, itetra, 0, n0);
- compute_tetrahedron_normal(vol, itetra, 1, n1);
- compute_tetrahedron_normal(vol, itetra, 2, n2);
- compute_tetrahedron_normal(vol, itetra, 3, n3);
+ volume_compute_tetrahedron_normal(vol, itetra, 0, n0);
+ volume_compute_tetrahedron_normal(vol, itetra, 1, n1);
+ volume_compute_tetrahedron_normal(vol, itetra, 2, n2);
+ volume_compute_tetrahedron_normal(vol, itetra, 3, n3);
}
exit:
@@ -254,6 +215,7 @@ setup_tetrahedral_mesh
res = buffer_init_from_data(vol->dev->allocator, &vol->prim_data,
&args->tetrahedron_data, args->ntetrahedra, args->context);
if(res != RES_OK) goto error;
+ vol->has_prim_data = 1;
}
/* Store the per vertex data */
@@ -261,6 +223,7 @@ setup_tetrahedral_mesh
res = buffer_init_from_data(vol->dev->allocator, &vol->vert_data,
&args->vertex_data, args->nvertices, args->context);
if(res != RES_OK) goto error;
+ vol->has_vert_data = 1;
}
exit:
diff --git a/src/suvm_volume.h b/src/suvm_volume.h
@@ -18,9 +18,10 @@
#include <rsys/dynamic_array_u32.h>
#include <rsys/dynamic_array_float.h>
+#include <rsys/float3.h>
#include <rsys/ref_count.h>
-#include <embree3/rtcore_device.h>
+#include <embree3/rtcore.h>
/* Forward declarations */
struct suvm_device;
@@ -54,9 +55,9 @@ struct ALIGN(16) node_inner {
/* Leaf node */
struct ALIGN(16) node_leaf {
float low[3];
- unsigned int geom_id;
+ uint32_t geom_id;
float upp[3];
- unsigned int prim_id;
+ uint32_t prim_id;
struct node node;
};
@@ -67,6 +68,9 @@ struct suvm_volume {
struct buffer prim_data; /* Per primitive data */
struct buffer vert_data; /* Per vertex data */
+ int has_prim_data;
+ int has_vert_data;
+
float low[3], upp[3]; /* Volume AABB */
RTCBVH bvh;
@@ -86,5 +90,103 @@ volume_get_primitives_count(const struct suvm_volume* vol)
return sz / 4;
}
+static FINLINE size_t
+volume_get_vertices_count(const struct suvm_volume* vol)
+{
+ size_t sz = 0;
+ ASSERT(vol);
+ sz = darray_float_size_get(&vol->positions);
+ ASSERT(sz % 3/*#coords per vertex*/ == 0);
+ return sz / 3;
+}
+
+static INLINE float*
+volume_compute_tetrahedron_normal
+ (const struct suvm_volume* vol,
+ const size_t itetra,
+ const int ifacet/*In [0,3]*/,
+ float normal[3])
+{
+ const uint32_t* ids;
+ const float* v0;
+ const float* v1;
+ const float* v2;
+
+ float E0[3], E1[3];
+
+ /* Facet indices, Ensure CCW ordering */
+ const int facet[4][3] = {{0,1,2}, {0,3,1}, {1,3,2}, {2,3,0}};
+
+ ASSERT(vol && itetra < volume_get_primitives_count(vol) && normal);
+
+ /* Fetch tetrahedron indices */
+ ids = darray_u32_cdata_get(&vol->indices) + itetra*4;
+
+ /* Fetch facet vertices */
+ v0 = darray_float_cdata_get(&vol->positions) + ids[facet[ifacet][0]]*3;
+ v1 = darray_float_cdata_get(&vol->positions) + ids[facet[ifacet][1]]*3;
+ v2 = darray_float_cdata_get(&vol->positions) + ids[facet[ifacet][2]]*3;
+
+ /* Compute the normal */
+ f3_sub(E0, v1, v0);
+ f3_sub(E1, v2, v0);
+ f3_cross(normal, E0, E1);
+ f3_normalize(normal, normal);
+
+ return normal;
+}
+
+static FINLINE float*
+volume_get_tetrahedron_normal
+ (const struct suvm_volume* vol,
+ const size_t itetra,
+ 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 {
+ 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];
+ normal[2] = darray_float_cdata_get(&vol->normals)[id+2];
+ }
+ return normal;
+}
+
+static INLINE const void*
+buffer_at(const struct buffer* buf, const size_t i)
+{
+ void* data;
+ ASSERT(buf && i < buf->size);
+ data = (char*)buf->mem + i * buf->elmt_stride;
+ ASSERT(IS_ALIGNED(data, buf->elmt_alignment));
+ return data;
+}
+
+static INLINE const void*
+volume_get_primitive_data(const struct suvm_volume* vol, const size_t iprim)
+{
+ ASSERT(vol && vol->has_prim_data && iprim < volume_get_primitives_count(vol));
+ return buffer_at(&vol->prim_data, iprim);
+}
+
+static INLINE const void*
+volume_primitive_get_vertex_data
+ (const struct suvm_volume* vol,
+ const size_t iprim,
+ const size_t ivert /* In [0, 3] */)
+{
+ const uint32_t* ids;
+ ASSERT(vol && vol->has_vert_data && iprim < volume_get_primitives_count(vol));
+ ASSERT(ivert < 4);
+ ids = darray_u32_cdata_get(&vol->indices) + iprim*4/*#vertex per tetra*/;
+ return buffer_at(&vol->vert_data, ids[ivert]);
+}
+
#endif /* SUVM_VOLUME_H */
diff --git a/src/suvm_volume_at.c b/src/suvm_volume_at.c
@@ -0,0 +1,210 @@
+/* Copyright (C) 2013-2020 Vincent Forest (vaplv@free.fr)
+ *
+ * The RSys library is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * The RSys library 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 Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with the RSys library. If not, see <http://www.gnu.org/licenses/>. */
+
+#include "suvm.h"
+#include "suvm_device.h"
+#include "suvm_volume.h"
+
+/* Generate the dynamic array of BVH node */
+#define DARRAY_NAME node
+#define DARRAY_DATA struct node*
+#include <rsys/dynamic_array.h>
+
+/*******************************************************************************
+ * Helper functions
+ ******************************************************************************/
+static FINLINE int
+pos_in_aabb
+ (const float pos[3],
+ const float low[3],
+ const float upp[3])
+{
+ return low[0] <= pos[0] && upp[0] >= pos[0]
+ && low[1] <= pos[1] && upp[1] >= pos[1]
+ && low[2] <= pos[2] && upp[2] >= pos[2];
+}
+
+/* Return 1 if pos lies is included into the leaf or 0 otherwise */
+static int
+intersect_leaf
+ (const struct suvm_volume* vol,
+ const struct node_leaf* leaf,
+ const float pos[3],
+ float bcoords[4]) /* Barycentric coordinates of pos into the leaf */
+{
+ const uint32_t* tetra;
+ const float *v0, *v1, *v2, *v3;
+ float n0[3], n1[3], n2[3], n3[3];
+ float p[3];
+ float d0, d1, d2, d3;
+ float h0, h1, h2, h3;
+ size_t itetra;
+ ASSERT(vol && leaf && pos && bcoords);
+ ASSERT(leaf->prim_id < volume_get_primitives_count(vol));
+
+ itetra = leaf->prim_id;
+ tetra = darray_u32_cdata_get(&vol->indices) + itetra*4/*#vertices per tetra*/;
+
+ /* Fetch the first tetrahedron vertex */
+ v0 = darray_float_cdata_get(&vol->positions) + tetra[0]*3;
+
+ /* Compute the distance from pos to the bottom facet */
+ volume_get_tetrahedron_normal(vol, itetra, 0, n0);
+ d0 = f3_dot(n0, f3_sub(p, pos, v0));
+ if(d0 < 0) return 0;
+
+ /* Fetch the top tetrahedron vertex */
+ v3 = darray_float_cdata_get(&vol->positions) + tetra[3]*3;
+ f3_sub(p, pos, v3);
+
+ /* Compute the distance from pos to the 'side' facets */
+ volume_get_tetrahedron_normal(vol, itetra, 1, n1);
+ d1 = f3_dot(n1, p);
+ if(d1 < 0) return 0;
+ volume_get_tetrahedron_normal(vol, itetra, 2, n2);
+ d2 = f3_dot(n2, p);
+ if(d2 < 0) return 0;
+ volume_get_tetrahedron_normal(vol, itetra, 3, n3);
+ d3 = f3_dot(n3, p);
+ if(d3 < 0) return 0;
+
+ /* Fetch the two remaining vertices (i.e. the second and third ones) */
+ v1 = darray_float_cdata_get(&vol->positions) + tetra[1]*3;
+ v2 = darray_float_cdata_get(&vol->positions) + tetra[1]*3;
+
+ /* Distance of tetrahedron corners to their opposite face */
+ h0 = f3_dot(n0, f3_sub(p, v3, v0));
+ h1 = f3_dot(n1, f3_sub(p, v2, v1));
+ h2 = f3_dot(n2, f3_sub(p, v0, v2));
+ h3 = f3_dot(n3, f3_sub(p, v1, v3));
+
+ /* Compute the barycentric coordinates, i.e. ratio of distances */
+ bcoords[0] = d2 / h2;
+ bcoords[1] = d3 / h3;
+ bcoords[2] = d1 / h1;
+ bcoords[3] = d0 / h0;
+
+ return 1;
+}
+
+/*******************************************************************************
+ * Exported function
+ ******************************************************************************/
+res_T
+suvm_volume_at
+ (struct suvm_volume* vol,
+ const double pos_in[3],
+ struct suvm_primitive* prim,
+ double uvw[4])
+{
+ struct darray_node stack;
+ struct node* node = NULL;
+ struct node_leaf* leaf = NULL;
+ struct node_inner* inner = NULL;
+ float pos[3];
+ size_t iprim = SIZE_MAX; /* Index of the hit primitive */
+ res_T res = RES_OK;
+
+ darray_node_init(vol->dev->allocator, &stack);
+
+ if(!vol || !pos_in || !prim || !uvw) {
+ res = RES_BAD_ARG;
+ goto error;
+ }
+
+ *prim = SUVM_PRIMITIVE_NULL;
+ pos[0] = (float)pos_in[0];
+ pos[1] = (float)pos_in[1];
+ pos[2] = (float)pos_in[2];
+
+ node = vol->bvh_root;
+ ASSERT(node);
+
+ /* Is the position included in the AABB of the volume */
+ if(!pos_in_aabb(pos, vol->low, vol->upp)) {
+ goto exit;
+ }
+
+ res = darray_node_reserve(&stack, 32);
+ if(res != RES_OK) goto error;
+
+ for(;;) {
+ size_t istack;
+
+ if(node->type == NODE_LEAF) {
+ float bcoords[4];
+ leaf = CONTAINER_OF(node, struct node_leaf, node);
+ if(intersect_leaf(vol, leaf, pos, bcoords)) {
+ iprim = leaf->prim_id;
+ uvw[0] = (double)bcoords[0];
+ uvw[1] = (double)bcoords[1];
+ uvw[2] = (double)bcoords[2];
+ break;
+ }
+
+ } else {
+ int traverse_child[2] = {0,0};
+ inner = CONTAINER_OF(node, struct node_inner, node);
+
+ /* Check pos against the children's AABBs */
+ traverse_child[0] = pos_in_aabb(pos, inner->low[0], inner->upp[0]);
+ traverse_child[1] = pos_in_aabb(pos, inner->low[1], inner->upp[1]);
+
+ if(traverse_child[0] && traverse_child[1]) { /* Traverse both children */
+ res = darray_node_push_back(&stack, &inner->children[1]);
+ if(UNLIKELY(res != RES_OK)) goto error;
+ node = inner->children[0];
+ continue;
+
+ } else if(traverse_child[0]) { /* Traverse only child 0 */
+ node = inner->children[1];
+ continue;
+
+ } else if(traverse_child[1]) { /* Traverse only child 1 */
+ node = inner->children[1];
+ continue;
+ }
+ }
+
+ istack = darray_node_size_get(&stack);
+ if(!istack) break; /* No more stack entry: stop traversal */
+
+ /* Pop the stack */
+ node = darray_node_data_get(&stack)[istack-1];
+ darray_node_pop_back(&stack);
+ }
+
+ /* pos lies into a tetrahedron */
+ if(iprim != SIZE_MAX) {
+ prim->iprim = iprim;
+ prim->nvertices = 4;
+ if(vol->has_prim_data) {
+ prim->data = volume_get_primitive_data(vol, iprim);
+ }
+ if(vol->has_vert_data) {
+ prim->vertex_data[0] = volume_primitive_get_vertex_data(vol, iprim, 0);
+ prim->vertex_data[1] = volume_primitive_get_vertex_data(vol, iprim, 1);
+ prim->vertex_data[2] = volume_primitive_get_vertex_data(vol, iprim, 2);
+ prim->vertex_data[3] = volume_primitive_get_vertex_data(vol, iprim, 3);
+ }
+ }
+
+exit:
+ darray_node_release(&stack);
+ return res;
+error:
+ goto exit;
+}
+
