commit 92ab8fe9483159c087a6a222e1097284caec480f
parent 21b67b362f5bf7750699bdbcf4fc627f729c3b1f
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Thu, 7 Jan 2021 12:08:20 +0100
Merge branch 'feature_tetrahedron_aabb_intersection'
Diffstat:
12 files changed, 1499 insertions(+), 153 deletions(-)
diff --git a/cmake/CMakeLists.txt b/cmake/CMakeLists.txt
@@ -26,6 +26,7 @@ option(NO_TEST "Do not build tests" OFF)
find_package(Embree 3.6 REQUIRED)
find_package(RCMake 0.4 REQUIRED)
find_package(RSys 0.10 REQUIRED)
+find_package(StarTetraHedra)
set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} ${RCMAKE_SOURCE_DIR})
include(rcmake)
@@ -47,6 +48,7 @@ set(VERSION ${VERSION_MAJOR}.${VERSION_MINOR}.${VERSION_PATCH})
set(SUVM_FILES_SRC
suvm_device.c
+ suvm_primitive.c
suvm_volume.c
suvm_volume_at.c
suvm_volume_intersect_aabb.c)
@@ -95,6 +97,12 @@ if(NOT NO_TEST)
new_test(test_suvm_device)
new_test(test_suvm_volume)
+ new_test(test_suvm_primitive_intersection)
+endif()
+
+if(StarTetraHedra_FOUND)
+ add_executable(suvm_voxelize ${SUVM_SOURCE_DIR}/suvm_voxelize_sth.c)
+ target_link_libraries(suvm_voxelize RSys StarTetraHedra suvm)
endif()
################################################################################
diff --git a/src/suvm.h b/src/suvm.h
@@ -17,6 +17,7 @@
#define SUVM_H
#include <rsys/rsys.h>
+#include <float.h>
/* Library symbol management */
#if defined(SUVM_SHARED_BUILD) /* Build shared library */
@@ -44,19 +45,48 @@ struct suvm_data {
size_t size; /* Size of the data in bytes */
size_t alignment; /* Alignment of the data */
};
-static const struct suvm_data SUVM_DATA_NULL;
+#define SUVM_DATA_NULL__ {NULL, 0, 0}
+static const struct suvm_data SUVM_DATA_NULL = SUVM_DATA_NULL__;
struct suvm_primitive {
const void* data; /* Data of the primitive */
const void* vertex_data[SUVM_PRIMITIVE_MAX_VERTICES_COUNT]; /* Vertex data */
size_t iprim; /* Identifier of the primitive */
size_t nvertices; /* #vertices of the primitive */
+
+ /* Internal data */
+ const struct suvm_volume* volume__;
};
-#define SUVM_PRIMITIVE_NULL__ {NULL, {NULL}, SIZE_MAX, SIZE_MAX}
+#define SUVM_PRIMITIVE_NULL__ { \
+ NULL, /* Primitive data */ \
+ {NULL}, /* Vertex data */ \
+ SIZE_MAX, /* Primitive id */ \
+ SIZE_MAX, /* #vertices */ \
+ NULL /* Pointer toward its associated volume */ \
+}
static const struct suvm_primitive SUVM_PRIMITIVE_NULL = SUVM_PRIMITIVE_NULL__;
#define SUVM_PRIMITIVE_NONE(Prim) ((Prim)->iprim == SIZE_MAX)
+/* Precomputed data from a suvm_primitive used to speed up AABB/primitive
+ * intersection tests */
+struct suvm_polyhedron {
+ float v[4/*#vertices*/][3/*#coords*/]; /* Vertices */
+ float N[4/*#vertices*/][3/*#coords*/]; /* Normals */
+ float D[4/*#facets*/]; /* Slope parameter of the plane (D = -(Ax+By+Cz)) */
+
+ /* 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[3/*#axes*/][4/*#edges max*/][3/*#equation parameters*/];
+
+ /* Number of silhouette edges in the ZY, XZ and XY planes (3 or 4) */
+ int nEp[3/*#coords*/];
+
+ /* Tetrahedron axis aligned bounding box */
+ float lower[3/*#coords*/];
+ float upper[3/*#coords*/];
+};
+
struct suvm_tetrahedral_mesh_args {
size_t ntetrahedra; /* #tetrahedra */
size_t nvertices; /* #vertices */
@@ -78,7 +108,18 @@ struct suvm_tetrahedral_mesh_args {
void* context; /* Client data set as the last param of the callbacks */
};
-static const struct suvm_tetrahedral_mesh_args SUVM_TETRAHEDRAL_MESH_ARGS_NULL;
+#define SUVM_TETRAHEDRAL_MESH_ARGS_NULL__ { \
+ 0, 0, NULL, NULL, SUVM_DATA_NULL__, SUVM_DATA_NULL__, 0, NULL \
+}
+static const struct suvm_tetrahedral_mesh_args SUVM_TETRAHEDRAL_MESH_ARGS_NULL =
+ SUVM_TETRAHEDRAL_MESH_ARGS_NULL__;
+
+enum suvm_intersection_type {
+ SUVM_INTERSECT_NONE,
+ SUVM_INTERSECT_INCLUDE,
+ SUVM_INTERSECT_IS_INCLUDED,
+ SUVM_INTERSECT_PARTIAL
+};
/* Callback invoked on suvm_volume_intersect_aabb invocation on each primitives
* intersected by the submitted Axis Aligned Bounding Box */
@@ -151,11 +192,8 @@ 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 intersecting 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,
@@ -164,6 +202,31 @@ suvm_volume_intersect_aabb
suvm_primitive_intersect_aabb_T clbk,
void* context); /* User data sent as the last argument of clbk */
+SUVM_API res_T
+suvm_volume_get_primitives_count
+ (const struct suvm_volume* volume,
+ size_t* nprims);
+
+SUVM_API res_T
+suvm_volume_get_primitive
+ (const struct suvm_volume* volume,
+ const size_t iprim, /* In [0, suvm_volume_get_primitives_count[ */
+ struct suvm_primitive* prim);
+
+/*******************************************************************************
+ * Primitive API
+ ******************************************************************************/
+SUVM_API res_T
+suvm_primitive_setup_polyhedron
+ (const struct suvm_primitive* prim,
+ struct suvm_polyhedron* poly);
+
+SUVM_API enum suvm_intersection_type
+suvm_polyhedron_intersect_aabb
+ (const struct suvm_polyhedron* poly,
+ const float low[3],
+ const float upp[3]);
+
END_DECLS
#endif /* SUVM_H */
diff --git a/src/suvm_device.c b/src/suvm_device.c
@@ -83,7 +83,7 @@ suvm_device_create
dev->logger = logger;
dev->verbose = verbose;
- sz = snprintf(embree_opts, sizeof(embree_opts), "verbose=%d", verbose);
+ sz = snprintf(embree_opts, sizeof(embree_opts), "verbose=%d", verbose > 1 ? 1 : 0);
if((size_t)sz >= sizeof(embree_opts)) {
log_err(dev, "Could not setup the Embree option string.\n");
res = RES_MEM_ERR;
diff --git a/src/suvm_primitive.c b/src/suvm_primitive.c
@@ -0,0 +1,335 @@
+/* Copyright (C) 2020 |Meso|Star> (contact@meso-star.com)
+ *
+ * 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 "suvm.h"
+#include "suvm_volume.h"
+
+#include <rsys/float2.h>
+#include <rsys/float3.h>
+
+/*******************************************************************************
+ * Exported functions
+ ******************************************************************************/
+res_T
+suvm_primitive_setup_polyhedron
+ (const struct suvm_primitive* prim,
+ struct suvm_polyhedron* poly)
+{
+ if(!prim || !poly) return RES_BAD_ARG;
+ tetrahedron_setup(prim->volume__, NULL, NULL, prim->iprim, poly);
+ return RES_OK;
+}
+
+/* Define if an axis aligned bounding box and a tetrahedron are intersecting.
+ * Its implementation follows the algorithms proposed by N. Greene in
+ * "Detecting intersection of a rectangular solid and a convex polyhedron"
+ * (Graphics Gems IV, 1994, p74--82) */
+enum suvm_intersection_type
+suvm_polyhedron_intersect_aabb
+ (const struct suvm_polyhedron* tetra,
+ const float low[3],
+ const float upp[3])
+{
+ enum { X, Y, Z, AXES_COUNT }; /* Syntactic sugar */
+ float intersect_aabb_low[AXES_COUNT];
+ float intersect_aabb_upp[AXES_COUNT];
+ int nplanes_including_aabb;
+ int i;
+
+ ASSERT(tetra && low && upp);
+ ASSERT(low[0] < upp[0]);
+ ASSERT(low[1] < upp[1]);
+ ASSERT(low[2] < upp[2]);
+
+ /* 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 SUVM_INTERSECT_NONE;
+ }
+
+ /* Check if the tetrahedron is included into the aabb */
+ if(intersect_aabb_low[X] == tetra->lower[X]
+ && intersect_aabb_low[Y] == tetra->lower[Y]
+ && intersect_aabb_low[Z] == tetra->lower[Z]
+ && intersect_aabb_upp[X] == tetra->upper[X]
+ && intersect_aabb_upp[Y] == tetra->upper[Y]
+ && intersect_aabb_upp[Z] == tetra->upper[Z]) {
+ return SUVM_INTERSECT_IS_INCLUDED;
+ }
+
+ /* #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 SUVM_INTERSECT_NONE;
+ }
+
+ /* 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 SUVM_INTERSECT_INCLUDE;
+ }
+
+ /* 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 SUVM_INTERSECT_NONE;
+ }
+ }
+ }
+
+ /* The tetrahedron and the aabb are partially intersecting */
+ return SUVM_INTERSECT_PARTIAL;
+}
+
+
+/*******************************************************************************
+ * Local functions
+ ******************************************************************************/
+void
+tetrahedron_setup
+ (const struct suvm_volume* vol,
+ const float lower[3],
+ const float upper[3],
+ const size_t itetra,
+ struct suvm_polyhedron* tetra)
+{
+ enum { X, Y, Z, AXES_COUNT }; /* Syntactic sugar */
+ const float* low = lower;
+ const float* upp = upper;
+ float low__[3];
+ float upp__[3];
+ float center[AXES_COUNT]; /* Center of the tetrahedron */
+ float e[6][AXES_COUNT];
+ float (*v)[AXES_COUNT];
+ float (*N)[AXES_COUNT];
+ float (*Ep)[4][3];
+ int i;
+ ASSERT(vol && tetra);
+
+ /* Fetch tetrahedron vertices */
+ volume_primitive_get_vertex_position(vol, itetra, 0, tetra->v[0]);
+ volume_primitive_get_vertex_position(vol, itetra, 1, tetra->v[1]);
+ volume_primitive_get_vertex_position(vol, itetra, 2, tetra->v[2]);
+ volume_primitive_get_vertex_position(vol, itetra, 3, tetra->v[3]);
+ v = tetra->v;
+
+ /* 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;
+
+ /* Define the primitive AABB if necessary */
+ if(!low) {
+ low__[0] = MMIN(MMIN(v[0][X], v[1][X]), MMIN(v[2][X], v[3][X]));
+ low__[1] = MMIN(MMIN(v[0][Y], v[1][Y]), MMIN(v[2][Y], v[3][Y]));
+ low__[2] = MMIN(MMIN(v[0][Z], v[1][Z]), MMIN(v[2][Z], v[3][Z]));
+ low = low__;
+ }
+ if(!upp) {
+ upp__[0] = MMAX(MMAX(v[0][X], v[1][X]), MMAX(v[2][X], v[3][X]));
+ upp__[1] = MMAX(MMAX(v[0][Y], v[1][Y]), MMAX(v[2][Y], v[3][Y]));
+ upp__[2] = MMAX(MMAX(v[0][Z], v[1][Z]), MMAX(v[2][Z], v[3][Z]));
+ upp = upp__;
+ }
+
+#ifndef NDEBUG
+ /* Check argument consistency */
+ if(lower) {
+ ASSERT(MMIN(MMIN(v[0][X], v[1][X]), MMIN(v[2][X], v[3][X])) == lower[X]);
+ ASSERT(MMIN(MMIN(v[0][Y], v[1][Y]), MMIN(v[2][Y], v[3][Y])) == lower[Y]);
+ ASSERT(MMIN(MMIN(v[0][Z], v[1][Z]), MMIN(v[2][Z], v[3][Z])) == lower[Z]);
+ }
+ if(upper) {
+ ASSERT(MMAX(MMAX(v[0][X], v[1][X]), MMAX(v[2][X], v[3][X])) == upper[X]);
+ ASSERT(MMAX(MMAX(v[0][Y], v[1][Y]), MMAX(v[2][Y], v[3][Y])) == upper[Y]);
+ ASSERT(MMAX(MMAX(v[0][Z], v[1][Z]), MMAX(v[2][Z], v[3][Z])) == upper[Z]);
+ }
+#endif
+
+ /* Setup tetrahedron AABB */
+ f3_set(tetra->lower, low);
+ f3_set(tetra->upper, upp);
+
+ /* Fetch tetrahedron normals */
+ volume_primitive_get_facet_normal(vol, itetra, 0, tetra->N[0]);
+ volume_primitive_get_facet_normal(vol, itetra, 1, tetra->N[1]);
+ volume_primitive_get_facet_normal(vol, itetra, 2, tetra->N[2]);
+ volume_primitive_get_facet_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) {
+ float c[2]; /* Projected tetrahedron center */
+
+ /* 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;
+
+ int iedge;
+
+ /* Project the tetrahedron center */
+ c[0] = center[j];
+ c[1] = center[k];
+
+ /* 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]));
+ Ep[i][n][2] = -(Ep[i][n][0]*v[0][j] + Ep[i][n][1]*v[0][k]);
+ ++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]));
+ Ep[i][n][2] = -(Ep[i][n][0]*v[1][j] + Ep[i][n][1]*v[1][k]);
+ ++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]));
+ Ep[i][n][2] = -(Ep[i][n][0]*v[2][j] + Ep[i][n][1]*v[2][k]);
+ ++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]));
+ Ep[i][n][2] = -(Ep[i][n][0]*v[0][j] + Ep[i][n][1]*v[0][k]);
+ ++n;
+ }
+ if(signf(N[1][i]) != signf(N[2][i])) { /* The edge 4 is silhouette */
+ f2_normalize(Ep[i][n], f2(Ep[i][n], e[4][k], -e[4][j]));
+ Ep[i][n][2] = -(Ep[i][n][0]*v[1][j] + Ep[i][n][1]*v[1][k]);
+ ++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]));
+ Ep[i][n][2] = -(Ep[i][n][0]*v[2][j] + Ep[i][n][1]*v[2][k]);
+ ++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 */
+
+ /* Ensure that the edge normals point toward the tetrahedron center */
+ FOR_EACH(iedge, 0, n) {
+ if(f2_dot(Ep[i][iedge], c) + Ep[i][iedge][2] < 0)
+ f3_minus(Ep[i][iedge], Ep[i][iedge]);
+ }
+ }
+}
+
diff --git a/src/suvm_volume.c b/src/suvm_volume.c
@@ -179,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 + 2)*3;
float* n3 = darray_float_data_get(&vol->normals) + (itetra*4 + 3)*3;
- 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);
+ volume_primitive_compute_facet_normal(vol, itetra, 0, n0);
+ volume_primitive_compute_facet_normal(vol, itetra, 1, n1);
+ volume_primitive_compute_facet_normal(vol, itetra, 2, n2);
+ volume_primitive_compute_facet_normal(vol, itetra, 3, n3);
}
exit:
@@ -426,65 +426,92 @@ error:
}
static res_T
-check_tetrahedra_normals(const struct suvm_volume* vol)
+fixup_tetrahedra_normals(struct suvm_volume* vol)
{
size_t itetra, ntetra;
- res_T res = RES_OK;
+ int fixup = 0;
ASSERT(vol);
ntetra = volume_get_primitives_count(vol);
FOR_EACH(itetra, 0, ntetra) {
- const uint32_t* ids = NULL;
- const float* verts[4];
- float normals[4][3];
- float pt[3];
- float v0[3], v1[3];
-
+ uint32_t* ids = NULL;
+ const float* vert0 = NULL;
+ const float* vert3 = NULL;
+ float normal0[3];
+ float v0[3];
/* Fetch tetrahedron indices */
- ids = darray_u32_cdata_get(&vol->indices) + itetra*4;
+ ids = darray_u32_data_get(&vol->indices) + itetra*4;
ASSERT(ids[0] < volume_get_vertices_count(vol));
ASSERT(ids[1] < volume_get_vertices_count(vol));
ASSERT(ids[2] < volume_get_vertices_count(vol));
ASSERT(ids[3] < volume_get_vertices_count(vol));
/* Fetch the tetrahedron vertices */
- verts[0] = darray_float_cdata_get(&vol->positions) + ids[0]*3/*#coords*/;
- verts[1] = darray_float_cdata_get(&vol->positions) + ids[1]*3/*#coords*/;
- verts[2] = darray_float_cdata_get(&vol->positions) + ids[2]*3/*#coords*/;
- verts[3] = darray_float_cdata_get(&vol->positions) + ids[3]*3/*#coords*/;
+ vert0 = darray_float_cdata_get(&vol->positions) + ids[0]*3/*#coords*/;
+ vert3 = darray_float_cdata_get(&vol->positions) + ids[3]*3/*#coords*/;
+
+ /* Fetch the normal of the 1st facet */
+ volume_primitive_get_facet_normal(vol, itetra, 0, normal0);
+
+ /* Check that the normal of the 1st facet looks the fourth vertex */
+ if(f3_dot(normal0, f3_sub(v0, vert3, vert0)) < 0) {
+ fixup = 1;
+
+ /* Revert tetrahedron orientation, i.e. the vertices of the 1st facet*/
+ SWAP(uint32_t, ids[1], ids[2]);
+
+ /* Revert precomputed normals if any */
+ if(darray_float_size_get(&vol->normals)) {
+ size_t ifacet;
+ FOR_EACH(ifacet, 0, 4) {
+ const size_t id = (itetra*4/*#facets*/ + ifacet)*3/*#coords*/;
+ float* N = darray_float_data_get(&vol->normals) + id;
+ N[0] = -N[0];
+ N[1] = -N[1];
+ N[2] = -N[2];
+ }
+ }
+ }
- /* Fetch the tetrahedron normals */
- volume_get_tetrahedron_normal(vol, itetra, 0, normals[0]);
- volume_get_tetrahedron_normal(vol, itetra, 1, normals[1]);
- volume_get_tetrahedron_normal(vol, itetra, 2, normals[2]);
- volume_get_tetrahedron_normal(vol, itetra, 3, normals[3]);
-
- /* Compute the position at the center of the tetrahedron */
- pt[0] = (verts[0][0] + verts[1][0] + verts[2][0] + verts[3][0]) * 0.25f;
- pt[1] = (verts[0][1] + verts[1][1] + verts[2][1] + verts[3][1]) * 0.25f;
- pt[2] = (verts[0][2] + verts[1][2] + verts[2][2] + verts[3][2]) * 0.25f;
-
- /* Check that normals look toward the tetrahedron center */
- f3_sub(v0, pt, verts[0]);
- f3_sub(v1, pt, verts[3]);
- if(f3_dot(normals[0], v0) < 0
- || f3_dot(normals[1], v1) < 0
- || f3_dot(normals[2], v1) < 0
- || f3_dot(normals[3], v1) < 0) {
- log_err(vol->dev,
- "Invalid vertex orientation for the tetrahedron %lu.\n",
- (unsigned long)itetra);
- res = RES_BAD_ARG;
- goto error;
+#ifndef NDEBUG
+ {
+ const float* verts[4];
+ float normals[4][3];
+ float pt[3], v1[3];
+
+ verts[0] = vert0;
+ verts[1] = darray_float_cdata_get(&vol->positions) + ids[1]*3/*#coords*/;
+ verts[2] = darray_float_cdata_get(&vol->positions) + ids[2]*3/*#coords*/;
+ verts[3] = vert3;
+
+ /* Compute the position at the center of the tetrahedron */
+ pt[0] = (verts[0][0] + verts[1][0] + verts[2][0] + verts[3][0]) * 0.25f;
+ pt[1] = (verts[0][1] + verts[1][1] + verts[2][1] + verts[3][1]) * 0.25f;
+ pt[2] = (verts[0][2] + verts[1][2] + verts[2][2] + verts[3][2]) * 0.25f;
+
+ /* Fetch the tetrahedron normals */
+ volume_primitive_get_facet_normal(vol, itetra, 0, normals[0]);
+ volume_primitive_get_facet_normal(vol, itetra, 1, normals[1]);
+ volume_primitive_get_facet_normal(vol, itetra, 2, normals[2]);
+ volume_primitive_get_facet_normal(vol, itetra, 3, normals[3]);
+
+ /* Check normals orientation */
+ f3_sub(v0, pt, verts[0]);
+ f3_sub(v1, pt, verts[3]);
+ ASSERT(f3_dot(normals[0], v0) >= 0);
+ ASSERT(f3_dot(normals[1], v1) >= 0);
+ ASSERT(f3_dot(normals[2], v1) >= 0);
+ ASSERT(f3_dot(normals[3], v1) >= 0);
}
+#endif
}
-exit:
- return res;
-error:
- goto exit;
+ if(fixup) {
+ log_warn(vol->dev, "Tetrahedrals were not correctly orientated.\n");
+ }
+ return RES_OK;
}
static void
@@ -540,7 +567,7 @@ suvm_tetrahedral_mesh_create
/* Ensure that the vertices of the tetrahedra are well ordered regarding the
* normal convention */
- res = check_tetrahedra_normals(vol);
+ res = fixup_tetrahedra_normals(vol);
if(res != RES_OK) goto error;
/* Build the BVH of the volumetric mesh */
@@ -607,3 +634,23 @@ suvm_volume_get_aabb
upper[2] = volume->upp[2];
return RES_OK;
}
+
+res_T
+suvm_volume_get_primitives_count(const struct suvm_volume* vol, size_t* nprims)
+{
+ if(!vol || !nprims) return RES_BAD_ARG;
+ *nprims = volume_get_primitives_count(vol);
+ return RES_OK;
+}
+
+res_T
+suvm_volume_get_primitive
+ (const struct suvm_volume* vol,
+ const size_t iprim,
+ struct suvm_primitive* prim)
+{
+ if(!vol|| !prim || iprim >= volume_get_primitives_count(vol))
+ return RES_BAD_ARG;
+ volume_primitive_setup(vol, iprim, prim);
+ return RES_OK;
+}
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,10 +121,26 @@ volume_get_vertices_count(const struct suvm_volume* vol)
return sz / 3;
}
+static FINLINE float*
+volume_primitive_get_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 && ivert < 4);
+ 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
+volume_primitive_compute_facet_normal
(const struct suvm_volume* vol,
- const size_t itetra,
+ const size_t iprim,
const int ifacet/*In [0,3]*/,
float normal[3])
{
@@ -122,10 +154,11 @@ volume_compute_tetrahedron_normal
/* 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);
+ ASSERT(vol && iprim < volume_get_primitives_count(vol) && normal);
+ ASSERT(ifacet < 4);
/* Fetch tetrahedron indices */
- ids = darray_u32_cdata_get(&vol->indices) + itetra*4;
+ ids = darray_u32_cdata_get(&vol->indices) + iprim*4;
/* Fetch facet vertices */
v0 = darray_float_cdata_get(&vol->positions) + ids[facet[ifacet][0]]*3;
@@ -142,19 +175,20 @@ volume_compute_tetrahedron_normal
}
static FINLINE float*
-volume_get_tetrahedron_normal
+volume_primitive_get_facet_normal
(const struct suvm_volume* vol,
- const size_t itetra,
- const int ifacet /* In [0, 3] */,
+ const size_t iprim,
+ 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);
+ ASSERT(vol && iprim < volume_get_primitives_count(vol) && normal);
+ ASSERT(ifacet < 4);
if(!darray_float_size_get(&vol->normals)) {
- volume_compute_tetrahedron_normal(vol, itetra, ifacet, normal);
+ volume_primitive_compute_facet_normal(vol, iprim, ifacet, normal);
} else {
+ const size_t id =
+ (iprim*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];
@@ -174,7 +208,7 @@ buffer_at(const struct buffer* buf, const size_t i)
}
static INLINE const void*
-volume_get_primitive_data(const struct suvm_volume* vol, const size_t iprim)
+volume_primitive_get_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);
@@ -193,5 +227,34 @@ volume_primitive_get_vertex_data
return buffer_at(&vol->vert_data, ids[ivert]);
}
+static INLINE void
+volume_primitive_setup
+ (const struct suvm_volume* vol,
+ const size_t iprim,
+ struct suvm_primitive* prim)
+{
+ ASSERT(vol && prim && iprim < volume_get_primitives_count(vol));
+ prim->nvertices = 4;
+ prim->iprim = iprim;
+ prim->volume__ = vol;
+ if(vol->has_prim_data) {
+ prim->data = volume_primitive_get_data(vol, prim->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);
+ }
+}
+
+extern LOCAL_SYM void
+tetrahedron_setup
+ (const struct suvm_volume* vol,
+ const float low[3], /* NULL <=> compute on the fly */
+ const float upp[3], /* NULL <=> compute on the fly */
+ const size_t iprim,
+ struct suvm_polyhedron* polyhedron);
+
#endif /* SUVM_VOLUME_H */
diff --git a/src/suvm_volume_at.c b/src/suvm_volume_at.c
@@ -56,7 +56,7 @@ intersect_leaf
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);
+ volume_primitive_get_facet_normal(vol, itetra, 0, n0);
d0 = f3_dot(n0, f3_sub(p, pos, v0));
if(d0 < 0) return 0;
@@ -65,13 +65,13 @@ intersect_leaf
f3_sub(p, pos, v3);
/* Compute the distance from pos to the 'side' facets */
- volume_get_tetrahedron_normal(vol, itetra, 1, n1);
+ volume_primitive_get_facet_normal(vol, itetra, 1, n1);
d1 = f3_dot(n1, p);
if(d1 < 0) return 0;
- volume_get_tetrahedron_normal(vol, itetra, 2, n2);
+ volume_primitive_get_facet_normal(vol, itetra, 2, n2);
d2 = f3_dot(n2, p);
if(d2 < 0) return 0;
- volume_get_tetrahedron_normal(vol, itetra, 3, n3);
+ volume_primitive_get_facet_normal(vol, itetra, 3, n3);
d3 = f3_dot(n3, p);
if(d3 < 0) return 0;
@@ -184,17 +184,7 @@ suvm_volume_at
/* 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);
- }
+ volume_primitive_setup(vol, iprim, prim);
}
exit:
diff --git a/src/suvm_volume_intersect_aabb.c b/src/suvm_volume_intersect_aabb.c
@@ -17,6 +17,9 @@
#include "suvm_device.h"
#include "suvm_volume.h"
+#include <rsys/float2.h>
+#include <rsys/float3.h>
+
/*******************************************************************************
* Helper functions
******************************************************************************/
@@ -93,21 +96,19 @@ suvm_volume_intersect_aabb
size_t istack;
if(node->type == NODE_LEAF) {
+ struct suvm_polyhedron tetra;
+ enum suvm_intersection_type intersect;
+
leaf = CONTAINER_OF(node, struct node_leaf, node);
- if(aabb_intersect_aabb(lowf, uppf, leaf->low, leaf->upp)) {
+
+ /* Check the intersection between the tetrahedron and the AABB */
+ tetrahedron_setup(vol, leaf->low, leaf->upp, leaf->prim_id, &tetra);
+ intersect = suvm_polyhedron_intersect_aabb(&tetra, lowf, uppf);
+
+ if(intersect != SUVM_INTERSECT_NONE) {
struct suvm_primitive prim = SUVM_PRIMITIVE_NULL;
/* Setup the intersected primitive */
- prim.nvertices = 4;
- prim.iprim = leaf->prim_id;
- if(vol->has_prim_data) {
- prim.data = volume_get_primitive_data(vol, prim.iprim);
- }
- if(vol->has_vert_data) {
- prim.vertex_data[0] = volume_primitive_get_vertex_data(vol, prim.iprim, 0);
- prim.vertex_data[1] = volume_primitive_get_vertex_data(vol, prim.iprim, 1);
- prim.vertex_data[2] = volume_primitive_get_vertex_data(vol, prim.iprim, 2);
- prim.vertex_data[3] = volume_primitive_get_vertex_data(vol, prim.iprim, 3);
- }
+ volume_primitive_setup(vol, leaf->prim_id, &prim);
/* Invoke the user callback on the intersected primitive */
clbk(&prim, low, upp, context);
}
diff --git a/src/suvm_voxelize_sth.c b/src/suvm_voxelize_sth.c
@@ -0,0 +1,423 @@
+/* Copyright (C) 2020 |Meso|Star> (contact@meso-star.com)
+ *
+ * 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/>. */
+
+#define _POSIX_C_SOURCE 200112L /* getopt/close support */
+
+#include "suvm.h"
+#include <star/sth.h>
+#include <rsys/clock_time.h>
+#include <rsys/cstr.h>
+#include <rsys/mem_allocator.h>
+
+#include <string.h>
+#include <unistd.h> /* getopt & close functions */
+
+struct args {
+ const char* input_filename;
+ const char* output_filename;
+ unsigned def[3];
+ int quit;
+ int verbose;
+};
+
+static const struct args ARGS_DEFAULT = {
+ NULL, /* Input file */
+ NULL, /* Output file */
+ {64, 64, 64}, /* Definition */
+ 0, /* Quit */
+ 0 /* Verbose */
+};
+
+/*******************************************************************************
+ * Helper functions
+ ******************************************************************************/
+static void
+print_help(const char* cmd)
+{
+ ASSERT(cmd);
+ printf(
+"Usage: %s [options] [file]\n"
+"Voxelize the submitted Star-TetraHedra geometry and save the result\n"
+"into a VTK file.\n",
+ cmd);
+ printf("\n");
+ printf(
+" -d X:Y:Z voxelisation definition along the 3 axes.\n"
+" The default definition is %u:%u:%u.\n",
+ ARGS_DEFAULT.def[0],
+ ARGS_DEFAULT.def[1],
+ ARGS_DEFAULT.def[2]);
+ printf(
+" -h display this help and exit.\n");
+ printf(
+" -o <output> filename of the output VTK file. If not defined, write\n"
+" results to standard ouput\n");
+ printf(
+" -v make the program verobse.\n");
+ printf("\n");
+ printf(
+"Copyright (C) 2020 |Meso|Star> <contact@meso-star.com>.\n"
+"This is free software released under the GNU GPL license, version 3 or\n"
+"later. You are free to change or redistribute it under certain\n"
+"conditions <http://gnu.org.licenses/gpl.html>\n");
+}
+
+static void
+args_release(struct args* args)
+{
+ ASSERT(args);
+ *args = ARGS_DEFAULT;
+}
+
+static res_T
+args_init(struct args* args, const int argc, char** argv)
+{
+ size_t n;
+ int opt;
+ res_T res = RES_OK;
+ ASSERT(args);
+
+ *args = ARGS_DEFAULT;
+
+ while((opt = getopt(argc, argv, "d:ho:v")) != -1) {
+ switch(opt) {
+ case 'd':
+ res = cstr_to_list_uint(optarg, ':', args->def, &n, 3);
+ if(res == RES_OK && n != 3) res = RES_BAD_ARG;
+ if(!args->def[0] || !args->def[1] || !args->def[2]) res = RES_BAD_ARG;
+ break;
+ case 'h':
+ print_help(argv[0]);
+ args_release(args);
+ args->quit = 1;
+ break;
+ case 'o':
+ args->output_filename = optarg;
+ break;
+ case 'v': args->verbose = 1; break;
+ default: res = RES_BAD_ARG; break;
+ }
+ if(res != RES_OK) {
+ if(optarg) {
+ fprintf(stderr, "%s: invalid option argument '%s' -- '%c'\n",
+ argv[0], optarg, opt);
+ }
+ goto error;
+ }
+ }
+
+ if(optind < argc) {
+ args->input_filename = argv[optind];
+ }
+
+exit:
+ optind = 1;
+ return res;
+error:
+ args_release(args);
+ goto exit;
+}
+
+static void
+get_indices(const size_t itetra, size_t ids[4], void* ctx)
+{
+ const uint64_t* ui64;
+ struct sth_desc* desc = ctx;
+ ASSERT(desc);
+ ui64 = sth_desc_get_tetrahedron_indices(desc, itetra);
+ ids[0] = (size_t)ui64[0];
+ ids[1] = (size_t)ui64[1];
+ ids[2] = (size_t)ui64[2];
+ ids[3] = (size_t)ui64[3];
+}
+
+static void
+get_position(const size_t ivert, double pos[3], void* ctx)
+{
+ const double* dbl;
+ struct sth_desc* desc = ctx;
+ ASSERT(desc);
+ dbl = sth_desc_get_vertex_position(desc, ivert);
+ pos[0] = dbl[0];
+ pos[1] = dbl[1];
+ pos[2] = dbl[2];
+}
+
+static void
+write_voxels
+ (FILE* stream,
+ const int* vxls,
+ const unsigned def[3],
+ const double vxl_sz[3],
+ const double org[3])
+{
+ size_t ix, iy, iz;
+ size_t ivxl;
+ size_t i;
+ ASSERT(stream && vxls && def && def[0] && def[1] && def[2]);
+
+ fprintf(stream, "# vtk DataFile Version 2.0\n");
+ fprintf(stream, "Volumetric grid\n");
+ fprintf(stream, "ASCII\n");
+
+ fprintf(stream, "DATASET RECTILINEAR_GRID\n");
+ fprintf(stream, "DIMENSIONS %u %u %u\n", def[0]+1, def[1]+1, def[2]+1);
+ fprintf(stream, "X_COORDINATES %u float\n", def[0]+1);
+ FOR_EACH(i, 0, def[0]+1) fprintf(stream, "%g ", (double)i*vxl_sz[0] + org[0]);
+ fprintf(stream, "\n");
+ fprintf(stream, "Y_COORDINATES %u float\n", def[1]+1);
+ FOR_EACH(i, 0, def[1]+1) fprintf(stream, "%g ", (double)i*vxl_sz[1] + org[1]);
+ fprintf(stream, "\n");
+ fprintf(stream, "Z_COORDINATES %u float\n", def[2]+1);
+ FOR_EACH(i, 0, def[2]+1) fprintf(stream, "%g ", (double)i*vxl_sz[2] + org[2]);
+ fprintf(stream, "\n");
+
+ fprintf(stream, "CELL_DATA %u\n", def[0]*def[1]*def[2]);
+ fprintf(stream, "SCALARS intersect_type int 1\n");
+ fprintf(stream, "LOOKUP_TABLE default\n");
+
+ ivxl = 0;
+ FOR_EACH(iz, 0, def[2]) {
+ FOR_EACH(iy, 0, def[1]) {
+ FOR_EACH(ix, 0, def[0]) {
+ fprintf(stream, "%i\n", vxls[ivxl]);
+ ++ivxl;
+ }
+ }
+ }
+}
+
+static res_T
+voxelize_suvm_volume
+ (const struct suvm_volume* vol,
+ const unsigned def[3],
+ int* vxls)
+{
+ double low[3], upp[3];
+ double vxl_sz[3];
+ size_t iprim;
+ size_t nprims;
+ int progress = 0;
+ res_T res = RES_OK;
+ ASSERT(vol && def && def[0] && def[1] && def[2] && vxls);
+
+ /* Compute the voxel size */
+ res = suvm_volume_get_aabb(vol, low, upp);
+ if(res != RES_OK) goto error;
+ vxl_sz[0] = (upp[0] - low[0]) / (double)def[0];
+ vxl_sz[1] = (upp[1] - low[1]) / (double)def[1];
+ vxl_sz[2] = (upp[2] - low[2]) / (double)def[2];
+
+ res = suvm_volume_get_primitives_count(vol, &nprims);
+ if(res != RES_OK) goto error;
+
+ FOR_EACH(iprim, 0, nprims) {
+ struct suvm_primitive prim = SUVM_PRIMITIVE_NULL;
+ struct suvm_polyhedron poly;
+ size_t ivxl_low[3];
+ size_t ivxl_upp[3];
+ size_t ivxl[3];
+ size_t i = 0;
+ int pcent;
+
+ CHK(suvm_volume_get_primitive(vol, iprim, &prim) == RES_OK);
+ CHK(suvm_primitive_setup_polyhedron(&prim, &poly) == RES_OK);
+
+ /* Transform the polyhedron AABB in voxel space */
+ ivxl_low[0] = (size_t)((poly.lower[0] - low[0]) / vxl_sz[0]);
+ ivxl_low[1] = (size_t)((poly.lower[1] - low[1]) / vxl_sz[1]);
+ ivxl_low[2] = (size_t)((poly.lower[2] - low[2]) / vxl_sz[2]);
+ ivxl_upp[0] = (size_t)ceil((poly.upper[0] - low[0]) / vxl_sz[0]);
+ ivxl_upp[1] = (size_t)ceil((poly.upper[1] - low[1]) / vxl_sz[1]);
+ ivxl_upp[2] = (size_t)ceil((poly.upper[2] - low[2]) / vxl_sz[2]);
+ CHK(ivxl_low[0] < def[0] && ivxl_upp[0] <= def[0]);
+ CHK(ivxl_low[1] < def[1] && ivxl_upp[1] <= def[1]);
+ CHK(ivxl_low[2] < def[2] && ivxl_upp[2] <= def[2]);
+
+ FOR_EACH(ivxl[2], ivxl_low[2], ivxl_upp[2]) {
+ float vxl_low[3];
+ float vxl_upp[3];
+ vxl_low[0] = (float)low[0];
+ vxl_low[1] = (float)low[1];
+ vxl_low[2] = (float)((double)ivxl[2] * vxl_sz[2] + low[2]);
+ vxl_upp[0] = (float)upp[0];
+ vxl_upp[1] = (float)upp[1];
+ vxl_upp[2] = vxl_low[2] + (float)vxl_sz[2];
+
+ FOR_EACH(ivxl[1], ivxl_low[1], ivxl_upp[1]) {
+ vxl_low[1] = (float)((double)ivxl[1] * vxl_sz[1] + low[1]);
+ vxl_upp[1] = vxl_low[1] + (float)vxl_sz[1];
+ FOR_EACH(ivxl[0], ivxl_low[0], ivxl_upp[0]) {
+ vxl_low[0] = (float)((double)ivxl[0] * vxl_sz[0] + low[0]);
+ vxl_upp[0] = vxl_low[0] + (float)vxl_sz[0];
+
+ i = ivxl[0] + ivxl[1]*def[0] + ivxl[2]*def[0]*def[1];
+ vxls[i] += (int)suvm_polyhedron_intersect_aabb(&poly, vxl_low, vxl_upp);
+ }
+ }
+ }
+ pcent = (int)((double)iprim * 100 / (double)nprims + 0.5/*round*/);
+ if(pcent > progress) {
+ progress = pcent;
+ fprintf(stdout, "Voxelizing: %3d%%\r", progress);
+ fflush(stdout);
+ }
+ }
+
+ printf("Voxelizing: %3d%%\n", progress);
+
+exit:
+ return res;
+error:
+ goto exit;
+}
+
+/*******************************************************************************
+ * Main functions
+ ******************************************************************************/
+int
+main(int argc, char** argv)
+{
+ char dump[128];
+ struct args args = ARGS_DEFAULT;
+ FILE* stream_out = stdout;
+ FILE* stream_in = stdin;
+ const char* stream_out_name = "stdout";
+ const char* stream_in_name = "stdin";
+ struct sth* sth = NULL;
+ struct sth_desc desc = STH_DESC_NULL;
+ struct suvm_device* suvm = NULL;
+ struct suvm_volume* vol = NULL;
+ struct suvm_tetrahedral_mesh_args msh_args = SUVM_TETRAHEDRAL_MESH_ARGS_NULL;
+ struct time t0, t1;
+ int* vxls = NULL;
+ double low[3];
+ double upp[3];
+ double vxl_sz[3];
+ res_T res = RES_OK;
+ int err = 0;
+
+ res = args_init(&args, argc, argv);
+ if(res != RES_OK) goto error;
+ if(args.quit) goto exit;
+
+ /* Setup input stream */
+ if(args.input_filename) {
+ stream_in = fopen(args.input_filename, "r");
+ if(!stream_in) {
+ fprintf(stderr, "Could not open the input file '%s' -- %s.\n",
+ args.input_filename, strerror(errno));
+ goto error;
+ }
+ stream_in_name = args.input_filename;
+ }
+
+ /* Setup output stream */
+ if(args.output_filename) {
+ stream_out = fopen(args.output_filename, "w");
+ if(!stream_out) {
+ fprintf(stderr, "Could not open the output file '%s' -- %s.\n",
+ args.output_filename, strerror(errno));
+ goto error;
+ }
+ stream_out_name = args.output_filename;
+ (void)stream_out_name;
+ }
+
+
+ /* Load the submitted file */
+ res = sth_create(NULL, NULL, args.verbose, &sth);
+ if(res != RES_OK) goto error;
+ res = sth_load_stream(sth, stream_in, stream_in_name);
+ if(res != RES_OK) goto error;
+ res = sth_get_desc(sth, &desc);
+ if(res != RES_OK) goto error;
+ if(args.verbose) {
+ printf("#nodes: %lu; #tetrahedra: %lu\n", desc.nvertices, desc.ntetrahedra);
+ }
+
+ res = suvm_device_create(NULL, NULL, args.verbose, &suvm);
+ if(res != RES_OK) goto error;
+
+ msh_args.nvertices = desc.nvertices;
+ msh_args.ntetrahedra = desc.ntetrahedra;
+ msh_args.get_indices = get_indices;
+ msh_args.get_position = get_position;
+ msh_args.context = &desc;
+
+ /* Setup the volume from the loaded data */
+ if(args.verbose) {
+ printf("Partitionning the tetrahedral mesh '%s'\n", stream_in_name);
+ }
+ time_current(&t0);
+ res = suvm_tetrahedral_mesh_create(suvm, &msh_args, &vol);
+ if(res != RES_OK) goto error;
+ time_sub(&t0, time_current(&t1), &t0);
+ time_dump(&t0, TIME_ALL, NULL, dump, sizeof(dump));
+ if(args.verbose) {
+ printf("Build acceleration structure in %s\n", dump);
+ }
+
+ /* Allocate the list of voxels */
+ vxls = mem_calloc(args.def[0]*args.def[1]*args.def[2], sizeof(*vxls));
+ if(!vxls) {
+ fprintf(stderr, "Could not allocate the list of voxels.\n");
+ res = RES_MEM_ERR;
+ goto error;
+ }
+
+ /* Voxelize the volume */
+ time_current(&t0);
+ res = voxelize_suvm_volume(vol, args.def, vxls);
+ if(res != RES_OK) goto error;
+ time_sub(&t0, time_current(&t1), &t0);
+ time_dump(&t0, TIME_ALL, NULL, dump, sizeof(dump));
+ if(args.verbose) {
+ printf("Volume voxelized in %s\n", dump);
+ }
+
+ /* Write the voxels */
+ if(args.verbose) {
+ printf("Writing voxels in '%s'\n", stream_out_name);
+ }
+ SUVM(volume_get_aabb(vol, low, upp));
+ vxl_sz[0] = (upp[0] - low[0]) / (double)args.def[0];
+ vxl_sz[1] = (upp[1] - low[1]) / (double)args.def[1];
+ vxl_sz[2] = (upp[2] - low[2]) / (double)args.def[2];
+ time_current(&t0);
+ write_voxels(stream_out, vxls, args.def, vxl_sz, low);
+ time_sub(&t0, time_current(&t1), &t0);
+ time_dump(&t0, TIME_ALL, NULL, dump, sizeof(dump));
+ if(args.verbose) {
+ printf("Voxels written in %s\n", dump);
+ }
+
+exit:
+ if(stream_out && stream_out != stdout) fclose(stream_out);
+ if(stream_in && stream_in != stdin) fclose(stream_in);
+ if(suvm) SUVM(device_ref_put(suvm));
+ if(vol) SUVM(volume_ref_put(vol));
+ if(sth) STH(ref_put(sth));
+ if(vxls) mem_rm(vxls);
+ args_release(&args);
+ if(mem_allocated_size() != 0) {
+ fprintf(stderr, "Memory leaks: %lu Bytes.\n",
+ (unsigned long)mem_allocated_size());
+ }
+ return err;
+error:
+ err = -1;
+ goto exit;
+}
diff --git a/src/test_suvm_box.h b/src/test_suvm_box.h
@@ -0,0 +1,52 @@
+/* Copyright (C) 2020 |Meso|Star> (contact@meso-star.com)
+ *
+ * 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/>. */
+
+#ifndef TEST_SUVM_BOX_H
+#define TEST_SUVM_BOX_H
+
+static const double box_vertices[9/*#vertices*/*3/*#coords per vertex*/] = {
+ 0.0, 0.0, 0.0,
+ 1.0, 0.0, 0.0,
+ 0.0, 1.0, 0.0,
+ 1.0, 1.0, 0.0,
+ 0.0, 0.0, 1.0,
+ 1.0, 0.0, 1.0,
+ 0.0, 1.0, 1.0,
+ 1.0, 1.0, 1.0,
+ 0.5, 0.5, 0.5
+};
+static const size_t box_nverts = sizeof(box_vertices) / sizeof(double[3]);
+
+/* The following array lists the indices toward the 3D vertices of each
+ * boundary triangle. The index 8 references the vertex at the center of the
+ * box.
+ * ,2---,3 ,2----3
+ * ,' | ,'/| ,'/| \ | Y
+ * 6----7' / | 6' / | \ | |
+ * |', | / ,1 | / ,0---,1 o--X
+ * | ',|/,' |/,' | ,' /
+ * 4----5' 4----5' Z */
+static const size_t box_indices[12/*#tetras*/*4/*#indices per tetra*/] = {
+ 0, 1, 2, 8, 1, 3, 2, 8, /* -Z */
+ 0, 2, 4, 8, 2, 6, 4, 8, /* -X */
+ 4, 6, 5, 8, 6, 7, 5, 8, /* +Z */
+ 3, 5, 7, 8, 5, 3, 1, 8, /* +X */
+ 2, 7, 6, 8, 7, 2, 3, 8, /* +Y */
+ 0, 5, 1, 8, 5, 0, 4, 8 /* -Y */
+};
+static const size_t box_ntetras = sizeof(box_indices) / sizeof(size_t[4]);
+
+
+#endif /* TEST_SUVM_BOX_H */
diff --git a/src/test_suvm_primitive_intersection.c b/src/test_suvm_primitive_intersection.c
@@ -0,0 +1,310 @@
+/* Copyright (C) 2020 |Meso|Star> (contact@meso-star.com)
+ *
+ * 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 "suvm.h"
+#include "test_suvm_ball.h"
+#include "test_suvm_box.h"
+#include "test_suvm_utils.h"
+
+#include <rsys/float3.h>
+#include <string.h>
+
+struct mesh {
+ const double* vertices; /* List of double[3] */
+ size_t nvertices;
+ const size_t* tetrahedra; /* List of size_t[4] */
+ size_t ntetrahedra;
+};
+static const struct mesh MESH_NULL;
+
+/*******************************************************************************
+ * Helper functions
+ *************v*****************************************************************/
+static void
+get_indices(const size_t itetra, size_t ids[4], void* ctx)
+{
+ struct mesh* msh = ctx;
+ CHK(itetra < msh->ntetrahedra);
+ CHK(ids != NULL);
+ ids[0] = msh->tetrahedra[itetra*4+0];
+ ids[1] = msh->tetrahedra[itetra*4+1];
+ ids[2] = msh->tetrahedra[itetra*4+2];
+ ids[3] = msh->tetrahedra[itetra*4+3];
+}
+
+static void
+get_position(const size_t ivert, double pos[3], void* ctx)
+{
+ struct mesh* msh = ctx;
+ CHK(ctx != NULL);
+ CHK(pos != NULL);
+ CHK(ivert < msh->nvertices);
+ pos[0] = msh->vertices[ivert*3+0];
+ pos[1] = msh->vertices[ivert*3+1];
+ pos[2] = msh->vertices[ivert*3+2];
+}
+
+static void
+write_voxels
+ (FILE* stream,
+ const int* vxls,
+ const size_t def[3],
+ const double vxl_sz[3])
+{
+ size_t ix, iy, iz;
+ size_t ivxl;
+ size_t i;
+ CHK(stream && vxls && def && def[0] && def[1] && def[2]);
+
+ fprintf(stream, "# vtk DataFile Version 2.0\n");
+ fprintf(stream, "nothing\n");
+ fprintf(stream, "ASCII\n");
+
+ fprintf(stream, "DATASET RECTILINEAR_GRID\n");
+ fprintf(stream, "DIMENSIONS %lu %lu %lu\n", def[0]+1, def[1]+1, def[2]+1);
+ fprintf(stream, "X_COORDINATES %lu float\n", def[0]+1);
+ FOR_EACH(i, 0, def[0]+1) fprintf(stream, "%g ", (double)i*vxl_sz[0]);
+ fprintf(stream, "\n");
+ fprintf(stream, "Y_COORDINATES %lu float\n", def[1]+1);
+ FOR_EACH(i, 0, def[1]+1) fprintf(stream, "%g ", (double)i*vxl_sz[1]);
+ fprintf(stream, "\n");
+ fprintf(stream, "Z_COORDINATES %lu float\n", def[2]+1);
+ FOR_EACH(i, 0, def[2]+1) fprintf(stream, "%g ", (double)i*vxl_sz[2]);
+ fprintf(stream, "\n");
+
+ fprintf(stream, "CELL_DATA %lu\n", def[0]*def[1]*def[2]);
+ fprintf(stream, "SCALARS intersect_type int 1\n");
+ fprintf(stream, "LOOKUP_TABLE default\n");
+
+ ivxl = 0;
+ FOR_EACH(iz, 0, def[2]) {
+ FOR_EACH(iy, 0, def[1]) {
+ FOR_EACH(ix, 0, def[0]) {
+ fprintf(stream, "%i\n", vxls[ivxl]);
+ ++ivxl;
+ }
+ }
+ }
+}
+
+static void
+voxelise_volume
+ (const struct suvm_volume* vol,
+ int* vxls,
+ const size_t def[3])
+{
+ double low[3];
+ double upp[3];
+ double vxl_sz[3];
+ size_t iprim;
+ size_t nprims;
+
+ CHK(suvm_volume_get_aabb(vol, low, upp) == RES_OK);
+ CHK(suvm_volume_get_primitives_count(vol, &nprims) == RES_OK);
+
+ memset(vxls, 0, sizeof(*vxls)*def[0]*def[1]*def[2]);
+
+ vxl_sz[0] = (upp[0] - low[0]) / (double)def[0];
+ vxl_sz[1] = (upp[1] - low[1]) / (double)def[1];
+ vxl_sz[2] = (upp[2] - low[2]) / (double)def[2];
+
+ FOR_EACH(iprim, 0, nprims) {
+ struct suvm_primitive prim = SUVM_PRIMITIVE_NULL;
+ struct suvm_polyhedron poly;
+ size_t ivxl_low[3];
+ size_t ivxl_upp[3];
+ size_t ivxl[3];
+ size_t i = 0;
+
+ CHK(suvm_volume_get_primitive(vol, iprim, &prim) == RES_OK);
+ CHK(suvm_primitive_setup_polyhedron(&prim, &poly) == RES_OK);
+
+ /* Transform the polyhedron AABB in voxel space */
+ ivxl_low[0] = (size_t)((poly.lower[0] - low[0]) / vxl_sz[0]);
+ ivxl_low[1] = (size_t)((poly.lower[1] - low[1]) / vxl_sz[1]);
+ ivxl_low[2] = (size_t)((poly.lower[2] - low[2]) / vxl_sz[2]);
+ ivxl_upp[0] = (size_t)ceil((poly.upper[0] - low[0]) / vxl_sz[0]);
+ ivxl_upp[1] = (size_t)ceil((poly.upper[1] - low[1]) / vxl_sz[1]);
+ ivxl_upp[2] = (size_t)ceil((poly.upper[2] - low[2]) / vxl_sz[2]);
+ CHK(ivxl_low[0] < def[0] && ivxl_upp[0] <= def[0]);
+ CHK(ivxl_low[1] < def[1] && ivxl_upp[1] <= def[1]);
+ CHK(ivxl_low[2] < def[2] && ivxl_upp[2] <= def[2]);
+
+ FOR_EACH(ivxl[2], ivxl_low[2], ivxl_upp[2]) {
+ float vxl_low[3];
+ float vxl_upp[3];
+ vxl_low[2] = (float)((double)ivxl[2] * vxl_sz[2] + low[2]);
+ vxl_upp[2] = vxl_low[2] + (float)vxl_sz[2];
+ FOR_EACH(ivxl[1], ivxl_low[1], ivxl_upp[1]) {
+ vxl_low[1] = (float)((double)ivxl[1] * vxl_sz[1] + low[1]);
+ vxl_upp[1] = vxl_low[1] + (float)vxl_sz[1];
+ FOR_EACH(ivxl[0], ivxl_low[0], ivxl_upp[0]) {
+ vxl_low[0] = (float)((double)ivxl[0] * vxl_sz[0] + low[0]);
+ vxl_upp[0] = vxl_low[0] + (float)vxl_sz[0];
+
+ i = ivxl[0] + ivxl[1]*def[0] + ivxl[2]*def[0]*def[1];
+ vxls[i] += (int)suvm_polyhedron_intersect_aabb(&poly, vxl_low, vxl_upp);
+ }
+ }
+ }
+ }
+}
+
+static void
+test_mesh_voxelization
+ (struct suvm_device* dev,
+ struct mem_allocator* allocator,
+ struct mesh* msh,
+ const size_t def[3],
+ const char* filename) /* NULL <=> do not write the result onto disk */
+{
+ struct suvm_volume* vol = NULL;
+ struct suvm_tetrahedral_mesh_args args = SUVM_TETRAHEDRAL_MESH_ARGS_NULL;
+ double low[3], upp[3];
+ double vxl_sz[3];
+ int* vxls = NULL;
+
+ args.ntetrahedra = msh->ntetrahedra;
+ args.nvertices = msh->nvertices;
+ args.get_indices = get_indices;
+ args.get_position = get_position;
+ args.context = msh;
+
+ CHK(suvm_tetrahedral_mesh_create(dev, &args, &vol) == RES_OK);
+ CHK(suvm_volume_get_aabb(vol, low, upp) == RES_OK);
+ vxl_sz[0] = (upp[0] - low[0]) / (double)def[0];
+ vxl_sz[1] = (upp[1] - low[1]) / (double)def[1];
+ vxl_sz[2] = (upp[2] - low[2]) / (double)def[2];
+
+ CHK(allocator != NULL);
+ vxls = MEM_CALLOC(allocator, def[0]*def[1]*def[2], sizeof(*vxls));
+ CHK(vxls != NULL);
+
+ voxelise_volume(vol, vxls, def);
+
+ if(filename) {
+ FILE* fp = fopen(filename, "w");
+ CHK(fp != NULL);
+ write_voxels(fp, vxls, def, vxl_sz);
+ CHK(fclose(fp) == 0);
+ }
+
+ MEM_RM(allocator, vxls);
+ CHK(suvm_volume_ref_put(vol) == RES_OK);
+}
+
+
+static void
+test_tetra_aabb_intersection
+ (struct suvm_device* dev,
+ struct mem_allocator* allocator)
+{
+ const double vertices[] = {
+ 0.0, 0.0, 0.0,
+ 0.0, 0.0, 1.0,
+ 1.0, 0.0, 1.0,
+ 0.0, 1.0, 1.0
+ };
+ const size_t tetra[] = { 0, 1, 2, 3 };
+ struct mesh msh = MESH_NULL;
+ struct suvm_tetrahedral_mesh_args args = SUVM_TETRAHEDRAL_MESH_ARGS_NULL;
+ struct suvm_polyhedron poly;
+ struct suvm_primitive prim = SUVM_PRIMITIVE_NULL;
+ struct suvm_volume* vol= NULL;
+ float low[3];
+ float upp[3];
+ double vxl_sz[3];
+ const size_t def[3] = {16, 16, 16};
+ size_t nprims;
+ int* vxls = NULL;
+ (void)vxl_sz;
+
+ msh.vertices = vertices;
+ msh.nvertices = 4;
+ msh.tetrahedra = tetra;
+ msh.ntetrahedra = 1;
+
+ args.ntetrahedra = 1;
+ args.nvertices = 4;
+ args.get_indices = get_indices;
+ args.get_position = get_position;
+ args.context = &msh;
+
+ CHK(suvm_tetrahedral_mesh_create(dev, &args, &vol) == RES_OK);
+
+ CHK(suvm_volume_get_primitives_count(vol, &nprims) == RES_OK);
+ CHK(nprims == 1);
+
+ CHK(suvm_volume_get_primitive(vol, 0, &prim) == RES_OK);
+ CHK(suvm_primitive_setup_polyhedron(&prim, &poly) == RES_OK);
+
+ f3(low, 0.f, 0.f, 0.f);
+ f3(upp, 1.f, 1.f, 1.f);
+ CHK(suvm_polyhedron_intersect_aabb(&poly, low, upp) == SUVM_INTERSECT_IS_INCLUDED);
+
+ f3(low, 0.f, 0.f, 0.5f);
+ f3(upp, 0.5f, 0.5f, 1.f);
+ CHK(suvm_polyhedron_intersect_aabb(&poly, low, upp) == SUVM_INTERSECT_PARTIAL);
+
+ f3(low, 0.f, 0.f, 0.66667f /* ~1-1/3 */);
+ f3(upp, 0.33332f/*~1/3*/, 0.33332f/*~1-1/3*/, 1.f);
+ CHK(suvm_polyhedron_intersect_aabb(&poly, low, upp) == SUVM_INTERSECT_INCLUDE);
+
+ f3(low, 0.33334f/*~1.3*/, 0.33334f/* ~1/3 */, 0);
+ f3(upp, 1.f, 1.f, 0.66665f /*~ 1-1/3 */);
+ CHK(suvm_polyhedron_intersect_aabb(&poly, low, upp) == SUVM_INTERSECT_NONE);
+
+ CHK(allocator != NULL);
+ vxls = MEM_CALLOC(allocator, def[0]*def[1]*def[2], sizeof(*vxls));
+ CHK(vxls != NULL);
+
+ voxelise_volume(vol, vxls, def);
+
+ MEM_RM(allocator, vxls);
+ CHK(suvm_volume_ref_put(vol) == RES_OK);
+}
+
+/*******************************************************************************
+ * Main function
+ ******************************************************************************/
+int
+main(int argc, char** argv)
+{
+ struct mesh msh = MESH_NULL;
+ struct suvm_device* dev = NULL;
+ size_t def[3] = {64, 64, 64};
+ (void)argc, (void)argv;
+
+ CHK(suvm_device_create(NULL, &mem_default_allocator, 1, &dev) == RES_OK);
+
+ test_tetra_aabb_intersection(dev, &mem_default_allocator);
+
+ msh.vertices = box_vertices;
+ msh.nvertices = box_nverts;
+ msh.tetrahedra = box_indices;
+ msh.ntetrahedra = box_ntetras;
+ test_mesh_voxelization(dev, &mem_default_allocator, &msh, def, "box.vtk");
+
+ msh.vertices = ball_vertices;
+ msh.nvertices = ball_nvertices;
+ msh.tetrahedra = ball_tetrahedra;
+ msh.ntetrahedra = ball_ntetrahedra;
+ test_mesh_voxelization(dev, &mem_default_allocator, &msh, def, "ball.vtk");
+
+ CHK(suvm_device_ref_put(dev) == RES_OK);
+ check_memory_allocator(&mem_default_allocator);
+ CHK(mem_allocated_size() == 0);
+ return 0;
+}
diff --git a/src/test_suvm_volume.c b/src/test_suvm_volume.c
@@ -16,46 +16,16 @@
#include "suvm.h"
#include "test_suvm_utils.h"
#include "test_suvm_ball.h"
+#include "test_suvm_box.h"
#include <rsys/dynamic_array.h>
#include <rsys/double3.h>
+#include <rsys/float3.h>
#define DARRAY_NAME prim
#define DARRAY_DATA struct suvm_primitive
#include <rsys/dynamic_array.h>
-static const double box_vertices[9/*#vertices*/*3/*#coords per vertex*/] = {
- 0.0, 0.0, 0.0,
- 1.0, 0.0, 0.0,
- 0.0, 1.0, 0.0,
- 1.0, 1.0, 0.0,
- 0.0, 0.0, 1.0,
- 1.0, 0.0, 1.0,
- 0.0, 1.0, 1.0,
- 1.0, 1.0, 1.0,
- 0.5, 0.5, 0.5
-};
-static const size_t box_nverts = sizeof(box_vertices) / sizeof(double[3]);
-
-/* The following array lists the indices toward the 3D vertices of each
- * boundary triangle. The index 8 references the vertex at the center of the
- * box.
- * ,2---,3 ,2----3
- * ,' | ,'/| ,'/| \ | Y
- * 6----7' / | 6' / | \ | |
- * |', | / ,1 | / ,0---,1 o--X
- * | ',|/,' |/,' | ,' /
- * 4----5' 4----5' Z */
-static const size_t box_indices[12/*#tetras*/*4/*#indices per tetra*/] = {
- 0, 1, 2, 8, 1, 3, 2, 8, /* -Z */
- 0, 2, 4, 8, 2, 6, 4, 8, /* -X */
- 4, 6, 5, 8, 6, 7, 5, 8, /* +Z */
- 3, 5, 7, 8, 5, 3, 1, 8, /* +X */
- 2, 7, 6, 8, 7, 2, 3, 8, /* +Y */
- 0, 5, 1, 8, 5, 0, 4, 8 /* -Y */
-};
-static const size_t box_ntetras = sizeof(box_indices) / sizeof(size_t[4]);
-
struct mesh {
const double* vertices; /* List of double[3] */
size_t nvertices;
@@ -138,11 +108,75 @@ cmp_prim(const void* a, const void* b)
}
static void
+check_volume_polyhedra(const struct suvm_volume* vol, struct mesh* msh)
+{
+ size_t iprim;
+ size_t nprims;
+
+ CHK(suvm_volume_get_primitives_count(vol, &nprims) == RES_OK);
+ CHK(nprims == msh->ntetrahedra);
+
+ FOR_EACH(iprim, 0, nprims) {
+ struct suvm_primitive prim = SUVM_PRIMITIVE_NULL;
+ struct suvm_polyhedron poly;
+ double v[4][3];
+ double E0[3];
+ double E1[3];
+ double N[3];
+ size_t ids[4];
+ float f[3];
+ CHK(suvm_volume_get_primitive(vol, iprim, &prim) == RES_OK);
+ CHK(!SUVM_PRIMITIVE_NONE(&prim));
+ CHK(prim.iprim == iprim);
+ CHK(prim.nvertices == 4);
+ CHK(suvm_primitive_setup_polyhedron(&prim, &poly) == RES_OK);
+
+ get_indices(iprim, ids, msh);
+ get_position(ids[0], v[0], msh);
+ get_position(ids[1], v[1], msh);
+ get_position(ids[2], v[2], msh);
+ get_position(ids[3], v[3], msh);
+
+ CHK(f3_eq(poly.v[0], f3_set_d3(f, v[0])));
+ CHK(f3_eq(poly.v[1], f3_set_d3(f, v[1])));
+ CHK(f3_eq(poly.v[2], f3_set_d3(f, v[2])));
+ CHK(f3_eq(poly.v[3], f3_set_d3(f, v[3])));
+
+ d3_sub(E0, v[1], v[0]);
+ d3_sub(E1, v[2], v[0]);
+ d3_cross(N, E0, E1);
+ f3_normalize(f, f3_set_d3(f, N));
+ CHK(f3_eq_eps(poly.N[0], f, 1.e-4f));
+
+ d3_sub(E0, v[3], v[0]);
+ d3_sub(E1, v[1], v[0]);
+ d3_cross(N, E0, E1);
+ f3_normalize(f, f3_set_d3(f, N));
+ CHK(f3_eq_eps(poly.N[1], f, 1.e-4f));
+
+ d3_sub(E0, v[3], v[1]);
+ d3_sub(E1, v[2], v[1]);
+ d3_cross(N, E0, E1);
+ f3_normalize(f, f3_set_d3(f, N));
+ CHK(f3_eq_eps(poly.N[2], f, 1.e-4f));
+
+ d3_sub(E0, v[3], v[2]);
+ d3_sub(E1, v[0], v[2]);
+ d3_cross(N, E0, E1);
+ f3_normalize(f, f3_set_d3(f, N));
+ CHK(f3_eq_eps(poly.N[3], f, 1.e-4f));
+ }
+}
+
+static void
test_tetrahedral_mesh_creation(struct suvm_device* dev)
{
struct mesh msh = MESH_NULL;
struct suvm_tetrahedral_mesh_args args = SUVM_TETRAHEDRAL_MESH_ARGS_NULL;
struct suvm_volume* vol = NULL;
+ struct suvm_primitive prim = SUVM_PRIMITIVE_NULL;
+ struct suvm_polyhedron poly;
+ size_t nprims;
args.ntetrahedra = box_ntetras;
args.nvertices = box_nverts;
@@ -221,11 +255,31 @@ test_tetrahedral_mesh_creation(struct suvm_device* dev)
args.vertex_data.alignment = 32;
msh.vertex_data_alignment = 32;
CHK(suvm_tetrahedral_mesh_create(dev, &args, &vol) == RES_OK);
+
+ CHK(suvm_volume_get_primitives_count(NULL, &nprims) == RES_BAD_ARG);
+ CHK(suvm_volume_get_primitives_count(vol, NULL) == RES_BAD_ARG);
+ CHK(suvm_volume_get_primitives_count(vol, &nprims) == RES_OK);
+ CHK(nprims == msh.ntetrahedra);
+
+ CHK(suvm_volume_get_primitive(NULL, 0, &prim) == RES_BAD_ARG);
+ CHK(suvm_volume_get_primitive(vol, nprims, &prim) == RES_BAD_ARG);
+ CHK(suvm_volume_get_primitive(vol, 0, NULL) == RES_BAD_ARG);
+ CHK(suvm_volume_get_primitive(vol, 0, &prim) == RES_OK);
+
+ CHK(suvm_primitive_setup_polyhedron(NULL, &poly) == RES_BAD_ARG);
+ CHK(suvm_primitive_setup_polyhedron(&prim, NULL) == RES_BAD_ARG);
+ CHK(suvm_primitive_setup_polyhedron(&prim, &poly) == RES_OK);
+
+ check_volume_polyhedra(vol, &msh);
+
CHK(suvm_volume_ref_put(vol) == RES_OK);
}
static void
-check_prim(const struct suvm_primitive* prim, struct mesh* msh)
+check_prim
+ (const struct suvm_primitive* prim,
+ const struct suvm_volume* vol,
+ struct mesh* msh)
{
size_t ids[4];
double verts[4][3];
@@ -234,6 +288,8 @@ check_prim(const struct suvm_primitive* prim, struct mesh* msh)
/* Check primitive data */
CHK(prim != NULL);
+ CHK(vol != NULL);
+ CHK(prim->volume__ == vol);
CHK(!SUVM_PRIMITIVE_NONE(prim));
CHK(prim->data != NULL);
CHK(prim->data != NULL);
@@ -271,6 +327,7 @@ check_prim(const struct suvm_primitive* prim, struct mesh* msh)
static void
check_prim_at
(const struct suvm_primitive* prim,
+ const struct suvm_volume* vol,
const double pos[3],
const double bcoords[4],
struct mesh* msh)
@@ -283,7 +340,7 @@ check_prim_at
double p[3];
const double* vert_data[4];
- check_prim(prim, msh);
+ check_prim(prim, vol, msh);
/* Fetch tetrahedron vertices */
get_indices(prim->iprim, ids, msh);
@@ -345,6 +402,7 @@ check_prim_at
static void
check_prims_intersect_aabb
(struct darray_prim* primitives,
+ const struct suvm_volume* vol,
const double low[3],
const double upp[3],
struct mesh* msh)
@@ -357,13 +415,15 @@ check_prims_intersect_aabb
double vtx[4][3];
double prim_low[3];
double prim_upp[3];
+ float lowf[3];
+ float uppf[3];
CHK(primitives);
prims = darray_prim_data_get(primitives);
nprims = darray_prim_size_get(primitives);
FOR_EACH(iprim, 0, nprims) {
- check_prim(prims + iprim, msh);
+ check_prim(prims + iprim, vol, msh);
/* Fetch tetrahedron vertices */
get_indices(prims[iprim].iprim, ids, msh);
@@ -393,26 +453,19 @@ check_prims_intersect_aabb
* that the input primitives are effectively the same of the ones detected by
* this brute force procedure */
iprim_challenge = 0;
+ f3_set_d3(lowf, low);
+ f3_set_d3(uppf, upp);
FOR_EACH(iprim, 0, msh->ntetrahedra) {
- /* Fetch tetrahedron vertices */
- get_indices(iprim, ids, msh);
- get_position(ids[0], vtx[0], msh);
- get_position(ids[1], vtx[1], msh);
- get_position(ids[2], vtx[2], msh);
- get_position(ids[3], vtx[3], msh);
+ struct suvm_primitive prim;
+ struct suvm_polyhedron tetra;
+ enum suvm_intersection_type intersect;
- /* Compute the primitive AABB */
- prim_low[0] = MMIN(MMIN(vtx[0][0], vtx[1][0]), MMIN(vtx[2][0], vtx[3][0]));
- prim_low[1] = MMIN(MMIN(vtx[0][1], vtx[1][1]), MMIN(vtx[2][1], vtx[3][1]));
- prim_low[2] = MMIN(MMIN(vtx[0][2], vtx[1][2]), MMIN(vtx[2][2], vtx[3][2]));
- prim_upp[0] = MMAX(MMAX(vtx[0][0], vtx[1][0]), MMAX(vtx[2][0], vtx[3][0]));
- prim_upp[1] = MMAX(MMAX(vtx[0][1], vtx[1][1]), MMAX(vtx[2][1], vtx[3][1]));
- prim_upp[2] = MMAX(MMAX(vtx[0][2], vtx[1][2]), MMAX(vtx[2][2], vtx[3][2]));
+ CHK(suvm_volume_get_primitive(vol, iprim, &prim) == RES_OK);
+ CHK(suvm_primitive_setup_polyhedron(&prim, &tetra) == RES_OK);
- /* Check that the primitives intersect the AABB */
- if(prim_low[0] < upp[0] && prim_upp[0] > low[0]
- && prim_low[1] < upp[1] && prim_upp[1] > low[1]
- && prim_low[2] < upp[2] && prim_upp[2] > low[2]) {
+ /* Check that the primitive intersects the AABB */
+ intersect = suvm_polyhedron_intersect_aabb(&tetra, lowf, uppf);
+ if(intersect != SUVM_INTERSECT_NONE) {
CHK(iprim == prims[iprim_challenge].iprim);
++iprim_challenge;
}
@@ -485,7 +538,7 @@ test_volume_at_cube(struct suvm_device* dev)
CHK(suvm_volume_at(vol, pos, &prim, bcoords) == RES_OK);
CHK(!SUVM_PRIMITIVE_NONE(&prim));
CHK(prim.iprim == 2);
- check_prim_at(&prim, pos, bcoords, &msh);
+ check_prim_at(&prim, vol, pos, bcoords, &msh);
CHK(suvm_volume_get_aabb(NULL, low, upp) == RES_BAD_ARG);
CHK(suvm_volume_get_aabb(vol, NULL, upp) == RES_BAD_ARG);
@@ -499,7 +552,7 @@ test_volume_at_cube(struct suvm_device* dev)
samp[1] = low[1] + rand_canonic()*(upp[1] - low[1]);
samp[2] = low[2] + rand_canonic()*(upp[2] - low[2]);
CHK(suvm_volume_at(vol, samp, &prim, bcoords) == RES_OK);
- check_prim_at(&prim, samp, bcoords, &msh);
+ check_prim_at(&prim, vol, samp, bcoords, &msh);
}
pos[0] = upp[0] + 0.1;
@@ -523,7 +576,7 @@ test_volume_at_cube(struct suvm_device* dev)
(vol, low, upp, prim_intersect_aabb, &prims) == RES_OK);
CHK(darray_prim_size_get(&prims) == box_ntetras);
- check_prims_intersect_aabb(&prims, low, upp, &msh);
+ check_prims_intersect_aabb(&prims, vol, low, upp, &msh);
CHK(suvm_volume_intersect_aabb
(vol, upp, low, prim_intersect_aabb, &prims) == RES_BAD_ARG);
@@ -548,8 +601,7 @@ test_volume_at_cube(struct suvm_device* dev)
CHK(suvm_volume_intersect_aabb
(vol, low, upp, prim_intersect_aabb, &prims) == RES_OK);
CHK(darray_prim_size_get(&prims) == 10);
- check_prims_intersect_aabb(&prims, low, upp, &msh);
-
+ check_prims_intersect_aabb(&prims, vol, low, upp, &msh);
CHK(suvm_volume_ref_put(vol) == RES_OK);
darray_prim_release(&prims);
@@ -610,6 +662,8 @@ test_volume_at_ball(struct suvm_device* dev)
CHK(suvm_tetrahedral_mesh_create(dev, &args, &vol) == RES_OK);
+ check_volume_polyhedra(vol, &msh);
+
/* Check volume AABB */
CHK(suvm_volume_get_aabb(vol, low, upp) == RES_OK);
CHK(d3_eq_eps(aabb_low, low, aabb_sz*1.e-6));
@@ -626,7 +680,7 @@ test_volume_at_ball(struct suvm_device* dev)
if(!SUVM_PRIMITIVE_NONE(&prim)) {
check_at += 1;
- check_prim_at(&prim, samp, bcoords, &msh);
+ check_prim_at(&prim, vol, samp, bcoords, &msh);
}
}
CHK(check_at != 0);
@@ -657,7 +711,7 @@ test_volume_at_ball(struct suvm_device* dev)
(vol, aabb_low, aabb_upp, prim_intersect_aabb, &prims) == RES_OK);
check_intersect_aabb = darray_prim_size_get(&prims) != 0;
- check_prims_intersect_aabb(&prims, aabb_low, aabb_upp, &msh);
+ check_prims_intersect_aabb(&prims, vol, aabb_low, aabb_upp, &msh);
darray_prim_clear(&prims);
}
CHK(check_intersect_aabb != 0);
