commit 46d767820a86fea21ed0c1a77da29ecee91c7cb4
parent e9d897180d10179bcfb585e4bd57294fffa560e3
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Mon, 16 Apr 2018 14:53:15 +0200
Add and test the sdis_scene_boundary_project_pos function
This function computes the parametric coordinates of a world space
position projected onto a given primitive.
Diffstat:
4 files changed, 337 insertions(+), 30 deletions(-)
diff --git a/src/sdis.h b/src/sdis.h
@@ -404,9 +404,43 @@ SDIS_API res_T
sdis_scene_get_boundary_position
(const struct sdis_scene* scn,
const size_t iprim, /* Primitive index */
- const double uv[2], /* Parametric coordinate onto the pimitive */
+ const double uv[2], /* Parametric coordinate onto the primitive */
double pos[3]); /* World space position */
+/* Project a world space position onto a primitive wrt its normal and compute
+ * the parametric coordinates of the projected point onto the primitive. This
+ * function may help to define the probe position onto a boundary as expected
+ * by the sdis_solve_probe_boundary function.
+ *
+ * Note that the projected point can lie outside the submitted primitive. In
+ * this case, the parametric coordinates are clamped against the primitive
+ * boundaries in order to ensure that the returned parametric coordinates are
+ * valid according to the primitive. To ensure this, in 2D, the parametric
+ * coordinate is simply clamped to [0, 1]. In 3D, the `uv' coordinates are
+ * clamped against the triangle edges. For instance, let the
+ * following triangle whose vertices are `a', `b' and `c':
+ * , ,
+ * , B ,
+ * , ,
+ * b E1
+ * E0 / \ ,P
+ * / \,*^
+ * / \
+ * ....a-------c......
+ * ' '
+ * A ' E2 ' C
+ * ' '
+ * The projected point `P' is orthogonally wrapped to the edge `ab', `bc' or
+ * `ca' if it lies in the `E0', `E1' or `E2' region, respectively. If `P' is in
+ * the `A', `B' or `C' region, then it is taken back to the `a', `b' or `c'
+ * vertex, respectively. */
+SDIS_API res_T
+sdis_scene_boundary_project_position
+ (const struct sdis_scene* scn,
+ const size_t iprim,
+ const double pos[3],
+ double uv[]);
+
/*******************************************************************************
* An estimator stores the state of a simulation
******************************************************************************/
diff --git a/src/sdis_scene.c b/src/sdis_scene.c
@@ -39,6 +39,65 @@ struct geometry_context {
/*******************************************************************************
* Helper function
******************************************************************************/
+static void
+project_position
+ (const double V0[3],
+ const double E0[3],
+ const double N[3],
+ const double NxE1[3],
+ const double rcp_det,
+ const double pos[3],
+ double uvw[3])
+{
+ double T[3], Q[3], k;
+ ASSERT(V0 && E0 && N && NxE1 && pos && uvw);
+
+ /* Use Moller/Trumbore intersection test the compute the parametric
+ * coordinates of the intersection between the triangle and the ray
+ * `r = pos + N*d' */
+ d3_sub(T, pos, V0);
+ uvw[0] = d3_dot(T, NxE1) * rcp_det;
+ d3_cross(Q, T, E0);
+ uvw[1] = d3_dot(Q, N) * rcp_det;
+ uvw[2] = 1.0 - uvw[0] - uvw[1];
+
+ if(uvw[0] >= 0 && uvw[1] >= 0 && uvw[2] >= 0) {/* The ray hits the triangle */
+ ASSERT(eq_eps(uvw[0] + uvw[1] + uvw[2], 1.0, 1.e-6));
+ return;
+ }
+
+ /* Clamp barycentric coordinates to triangle edges */
+ if(uvw[0] >= 0) {
+ if(uvw[1] >= 0) {
+ k = 1.0 / (uvw[0] + uvw[1]);
+ uvw[0] *= k;
+ uvw[1] *= k;
+ uvw[2] = 0;
+ } else if( uvw[2] >= 0) {
+ k = 1.0 / (uvw[0] + uvw[2]);
+ uvw[0] *= k;
+ uvw[1] = 0;
+ uvw[2] *= k;
+ } else {
+ ASSERT(uvw[0] >= 1.f);
+ d3(uvw, 1, 0, 0);
+ }
+ } else if(uvw[1] >= 0) {
+ if(uvw[2] >= 0) {
+ k = 1.0 / (uvw[1] + uvw[2]);
+ uvw[0] = 0;
+ uvw[1] *= k;
+ uvw[2] *= k;
+ } else {
+ ASSERT(uvw[1] >= 1);
+ d3(uvw, 0, 1, 0);
+ }
+ } else {
+ ASSERT(uvw[2] >= 1);
+ d3(uvw, 0, 0, 1);
+ }
+}
+
/* Check that `hit' roughly lies on a vertex. For segments, a simple but
* approximative way is to test that its position have at least one barycentric
* coordinate roughly equal to 0 or 1. */
@@ -639,6 +698,69 @@ sdis_scene_get_boundary_position
return RES_OK;
}
+res_T
+sdis_scene_boundary_project_position
+ (const struct sdis_scene* scn,
+ const size_t iprim,
+ const double pos[],
+ double uv[])
+{
+ if(!scn || !pos || !uv) return RES_BAD_ARG;
+ if(iprim >= scene_get_primitives_count(scn)) return RES_BAD_ARG;
+
+ if(scene_is_2d(scn)) {
+ struct s2d_primitive prim;
+ struct s2d_attrib a;
+ double V[2][2]; /* Vertices */
+ double E[2][3]; /* V0->V1 and V0->pos */
+ double proj;
+
+ /* Retrieve the segment vertices */
+ S2D(scene_view_get_primitive(scn->s2d_view, (unsigned int)iprim, &prim));
+ S2D(primitive_get_attrib(&prim, S2D_POSITION, 0, &a)); d2_set_f2(V[0], a.value);
+ S2D(primitive_get_attrib(&prim, S2D_POSITION, 1, &a)); d2_set_f2(V[1], a.value);
+
+ /* Compute the parametric coordinate of the project of `pos' onto the
+ * segment.*/
+ d2_sub(E[0], V[1], V[0]);
+ d2_normalize(E[0], E[0]);
+ d2_sub(E[1], pos, V[0]);
+ proj = d2_dot(E[0], E[1]);
+
+ uv[0] = CLAMP(proj, 0, 1); /* Clamp the parametric coordinate in [0, 1] */
+
+ } else {
+ struct s3d_primitive prim;
+ struct s3d_attrib a;
+ double V[3][3]; /* Vertices */
+ double E[2][3]; /* V0->V1 and V0->V2 edges */
+ double N[3]; /* Normal */
+ double NxE1[3], rcp_det; /* Muller/Trumboer triangle parameters */
+ double uvw[3];
+
+ S3D(scene_view_get_primitive(scn->s3d_view, (unsigned int)iprim, &prim));
+ S3D(triangle_get_vertex_attrib(&prim, 0, S3D_POSITION, &a)); d3_set_f3(V[0], a.value);
+ S3D(triangle_get_vertex_attrib(&prim, 1, S3D_POSITION, &a)); d3_set_f3(V[1], a.value);
+ S3D(triangle_get_vertex_attrib(&prim, 2, S3D_POSITION, &a)); d3_set_f3(V[2], a.value);
+ d3_sub(E[0], V[1], V[0]);
+ d3_sub(E[1], V[2], V[0]);
+ d3_cross(N, E[0], E[1]);
+
+ /* Muller/Trumbore triangle parameters */
+ d3_cross(NxE1, N, E[1]);
+ rcp_det = 1.0 / d3_dot(NxE1, E[0]);
+
+ /* Use the Muller/Trumbore intersection test to project `pos' onto the
+ * triangle and to retrieve the parametric coordinates of the projection
+ * point */
+ project_position(V[0], E[0], N, NxE1, rcp_det, pos, uvw);
+
+ uv[0] = uvw[2];
+ uv[1] = uvw[0];
+ }
+ return RES_OK;
+}
+
/*******************************************************************************
* Local miscellaneous function
******************************************************************************/
diff --git a/src/test_sdis_scene.c b/src/test_sdis_scene.c
@@ -15,6 +15,9 @@
#include "sdis.h"
#include "test_sdis_utils.h"
+
+#include <rsys/double2.h>
+#include <rsys/double3.h>
#include <rsys/math.h>
struct context {
@@ -23,8 +26,14 @@ struct context {
struct sdis_interface* interf;
};
+static INLINE double
+rand_canonic()
+{
+ return (double)rand()/(double)((double)RAND_MAX+1);
+}
+
static INLINE void
-get_indices(const size_t itri, size_t ids[3], void* context)
+get_indices_3d(const size_t itri, size_t ids[3], void* context)
{
struct context* ctx = context;
CHK(ctx != NULL);
@@ -35,7 +44,17 @@ get_indices(const size_t itri, size_t ids[3], void* context)
}
static INLINE void
-get_position(const size_t ivert, double pos[3], void* context)
+get_indices_2d(const size_t iseg, size_t ids[2], void* context)
+{
+ struct context* ctx = context;
+ CHK(ctx != NULL);
+ CHK(iseg < square_nsegments);
+ ids[0] = ctx->indices[iseg*2+0];
+ ids[1] = ctx->indices[iseg*2+1];
+}
+
+static INLINE void
+get_position_3d(const size_t ivert, double pos[3], void* context)
{
struct context* ctx = context;
CHK(ctx != NULL);
@@ -46,6 +65,16 @@ get_position(const size_t ivert, double pos[3], void* context)
}
static INLINE void
+get_position_2d(const size_t ivert, double pos[2], void* context)
+{
+ struct context* ctx = context;
+ CHK(ctx != NULL);
+ CHK(ivert < square_nvertices);
+ pos[0] = ctx->positions[ivert*2+0];
+ pos[1] = ctx->positions[ivert*2+1];
+}
+
+static INLINE void
get_interface(const size_t itri, struct sdis_interface** bound, void* context)
{
struct context* ctx = context;
@@ -55,30 +84,15 @@ get_interface(const size_t itri, struct sdis_interface** bound, void* context)
*bound = ctx->interf;
}
-int
-main(int argc, char** argv)
+static void
+test_scene_3d(struct sdis_device* dev, struct sdis_interface* interf)
{
- struct mem_allocator allocator;
- struct sdis_device* dev = NULL;
- struct sdis_medium* solid = NULL;
- struct sdis_medium* fluid = NULL;
- struct sdis_interface* interf = NULL;
struct sdis_scene* scn = NULL;
- struct sdis_fluid_shader fluid_shader = DUMMY_FLUID_SHADER;
- struct sdis_solid_shader solid_shader = DUMMY_SOLID_SHADER;
- struct sdis_interface_shader interface_shader = DUMMY_INTERFACE_SHADER;
double lower[3], upper[3];
+ double uv0[2], uv1[2], pos[3], pos1[3];
struct context ctx;
size_t ntris, npos;
- (void)argc, (void)argv;
-
- CHK(mem_init_proxy_allocator(&allocator, &mem_default_allocator) == RES_OK);
- CHK(sdis_device_create(NULL, &allocator, 1, 0, &dev) == RES_OK);
-
- CHK(sdis_fluid_create(dev, &fluid_shader, NULL, &fluid) == RES_OK);
- CHK(sdis_solid_create(dev, &solid_shader, NULL, &solid) == RES_OK);
- CHK(sdis_interface_create
- (dev, solid, fluid, &interface_shader, NULL, &interf) == RES_OK);
+ size_t i;
ctx.positions = box_vertices;
ctx.indices = box_indices;
@@ -87,8 +101,8 @@ main(int argc, char** argv)
npos = box_nvertices;
#define CREATE sdis_scene_create
- #define IDS get_indices
- #define POS get_position
+ #define IDS get_indices_3d
+ #define POS get_position_3d
#define IFA get_interface
CHK(CREATE(NULL, 0, NULL, NULL, 0, NULL, &ctx, NULL) == RES_BAD_ARG);
@@ -115,17 +129,158 @@ main(int argc, char** argv)
CHK(eq_eps(upper[1], 1, 1.e-6));
CHK(eq_eps(upper[2], 1, 1.e-6));
+ uv0[0] = 0.3;
+ uv0[1] = 0.3;
+
+ CHK(sdis_scene_get_boundary_position(NULL, 6, uv0, pos) == RES_BAD_ARG);
+ CHK(sdis_scene_get_boundary_position(scn, 12, uv0, pos) == RES_BAD_ARG);
+ CHK(sdis_scene_get_boundary_position(scn, 6, NULL, pos) == RES_BAD_ARG);
+ CHK(sdis_scene_get_boundary_position(scn, 6, uv0, NULL) == RES_BAD_ARG);
+ CHK(sdis_scene_get_boundary_position(scn, 6, uv0, pos) == RES_OK);
+
+ CHK(sdis_scene_boundary_project_position(NULL, 6, pos, uv1) == RES_BAD_ARG);
+ CHK(sdis_scene_boundary_project_position(scn, 12, pos, uv1) == RES_BAD_ARG);
+ CHK(sdis_scene_boundary_project_position(scn, 6, NULL, uv1) == RES_BAD_ARG);
+ CHK(sdis_scene_boundary_project_position(scn, 6, pos, NULL) == RES_BAD_ARG);
+ CHK(sdis_scene_boundary_project_position(scn, 6, pos, uv1) == RES_OK);
+
+ CHK(d2_eq_eps(uv0, uv1, 1.e-6));
+
+ FOR_EACH(i, 0, 64) {
+ uv0[0] = rand_canonic();
+ uv0[1] = rand_canonic() * (1 - uv0[0]);
+
+ CHK(sdis_scene_get_boundary_position(scn, 4, uv0, pos) == RES_OK);
+ CHK(sdis_scene_boundary_project_position(scn, 4, pos, uv1) == RES_OK);
+ CHK(d2_eq_eps(uv0, uv1, 1.e-6));
+ }
+
+ pos[0] = 10;
+ pos[1] = 0.1;
+ pos[2] = 0.5;
+ CHK(sdis_scene_boundary_project_position(scn, 6, pos, uv1) == RES_OK);
+ CHK(sdis_scene_get_boundary_position(scn, 6, uv1, pos1) == RES_OK);
+ CHK(!d3_eq_eps(pos1, pos, 1.e-6));
+
CHK(sdis_scene_ref_get(NULL) == RES_BAD_ARG);
CHK(sdis_scene_ref_get(scn) == RES_OK);
CHK(sdis_scene_ref_put(NULL) == RES_BAD_ARG);
CHK(sdis_scene_ref_put(scn) == RES_OK);
CHK(sdis_scene_ref_put(scn) == RES_OK);
+}
+
+static void
+test_scene_2d(struct sdis_device* dev, struct sdis_interface* interf)
+{
+ struct sdis_scene* scn = NULL;
+ double lower[2], upper[2];
+ double u0, u1, pos[2];
+ struct context ctx;
+ size_t nsegs, npos;
+ size_t i;
+
+ ctx.positions = square_vertices;
+ ctx.indices = square_indices;
+ ctx.interf = interf;
+ nsegs = square_nsegments;
+ npos = square_nvertices;
+
+ #define CREATE sdis_scene_2d_create
+ #define IDS get_indices_2d
+ #define POS get_position_2d
+ #define IFA get_interface
+
+ CHK(CREATE(NULL, 0, NULL, NULL, 0, NULL, &ctx, NULL) == RES_BAD_ARG);
+ CHK(CREATE(dev, 0, IDS, IFA, npos, POS, &ctx, &scn) == RES_BAD_ARG);
+ CHK(CREATE(dev, nsegs, NULL, IFA, npos, POS, &ctx, &scn) == RES_BAD_ARG);
+ CHK(CREATE(dev, nsegs, IDS, NULL, npos, POS, &ctx, &scn) == RES_BAD_ARG);
+ CHK(CREATE(dev, nsegs, IDS, IFA, 0, POS, &ctx, &scn) == RES_BAD_ARG);
+ CHK(CREATE(dev, nsegs, IDS, IFA, npos, NULL, &ctx, &scn) == RES_BAD_ARG);
+ CHK(CREATE(dev, nsegs, IDS, IFA, npos, POS, &ctx, &scn) == RES_OK);
+
+ #undef CREATE
+ #undef IDS
+ #undef POS
+ #undef IFA
+
+ CHK(sdis_scene_get_aabb(NULL, lower, upper) == RES_BAD_ARG);
+ CHK(sdis_scene_get_aabb(scn, NULL, upper) == RES_BAD_ARG);
+ CHK(sdis_scene_get_aabb(scn, lower, NULL) == RES_BAD_ARG);
+ CHK(sdis_scene_get_aabb(scn, lower, upper) == RES_OK);
+ CHK(eq_eps(lower[0], 0, 1.e-6));
+ CHK(eq_eps(lower[1], 0, 1.e-6));
+ CHK(eq_eps(upper[0], 1, 1.e-6));
+ CHK(eq_eps(upper[1], 1, 1.e-6));
+
+ u0 = 0.5;
+
+ CHK(sdis_scene_get_boundary_position(NULL, 1, &u0, pos) == RES_BAD_ARG);
+ CHK(sdis_scene_get_boundary_position(scn, 4, &u0, pos) == RES_BAD_ARG);
+ CHK(sdis_scene_get_boundary_position(scn, 1, NULL, pos) == RES_BAD_ARG);
+ CHK(sdis_scene_get_boundary_position(scn, 1, &u0, NULL) == RES_BAD_ARG);
+ CHK(sdis_scene_get_boundary_position(scn, 1, &u0, pos) == RES_OK);
+
+ CHK(sdis_scene_boundary_project_position(NULL, 1, pos, &u1) == RES_BAD_ARG);
+ CHK(sdis_scene_boundary_project_position(scn, 4, pos, &u1) == RES_BAD_ARG);
+ CHK(sdis_scene_boundary_project_position(scn, 1, NULL, &u1) == RES_BAD_ARG);
+ CHK(sdis_scene_boundary_project_position(scn, 1, pos, NULL) == RES_BAD_ARG);
+ CHK(sdis_scene_boundary_project_position(scn, 1, pos, &u1) == RES_OK);
+
+ CHK(eq_eps(u0, u1, 1.e-6));
+
+ FOR_EACH(i, 0, 64) {
+ u0 = rand_canonic();
+
+ CHK(sdis_scene_get_boundary_position(scn, 2, &u0, pos) == RES_OK);
+ CHK(sdis_scene_boundary_project_position(scn, 2, pos, &u1) == RES_OK);
+ CHK(eq_eps(u0, u1, 1.e-6));
+ }
+
+ d2(pos, 5, 0.5);
+ CHK(sdis_scene_boundary_project_position(scn, 3, pos, &u0) == RES_OK);
+ CHK(eq_eps(u0, 0.5, 1.e-6));
+
+ d2(pos, 1, 2);
+ CHK(sdis_scene_boundary_project_position(scn, 3, pos, &u0) == RES_OK);
+ CHK(eq_eps(u0, 0, 1.e-6));
+
+ d2(pos, 1, -1);
+ CHK(sdis_scene_boundary_project_position(scn, 3, pos, &u0) == RES_OK);
+ CHK(eq_eps(u0, 1, 1.e-6));
+
+ CHK(sdis_scene_ref_put(scn) == RES_OK);
+}
+
+int
+main(int argc, char** argv)
+{
+ struct mem_allocator allocator;
+ struct sdis_device* dev = NULL;
+ struct sdis_medium* solid = NULL;
+ struct sdis_medium* fluid = NULL;
+ struct sdis_interface* interf = NULL;
+ struct sdis_fluid_shader fluid_shader = DUMMY_FLUID_SHADER;
+ struct sdis_solid_shader solid_shader = DUMMY_SOLID_SHADER;
+ struct sdis_interface_shader interface_shader = DUMMY_INTERFACE_SHADER;
+ (void)argc, (void)argv;
+
+ CHK(mem_init_proxy_allocator(&allocator, &mem_default_allocator) == RES_OK);
+ CHK(sdis_device_create(NULL, &allocator, 1, 0, &dev) == RES_OK);
+
+ CHK(sdis_fluid_create(dev, &fluid_shader, NULL, &fluid) == RES_OK);
+ CHK(sdis_solid_create(dev, &solid_shader, NULL, &solid) == RES_OK);
+ CHK(sdis_interface_create
+ (dev, solid, fluid, &interface_shader, NULL, &interf) == RES_OK);
- CHK(sdis_device_ref_put(dev) == RES_OK);
- CHK(sdis_interface_ref_put(interf) == RES_OK);
CHK(sdis_medium_ref_put(solid) == RES_OK);
CHK(sdis_medium_ref_put(fluid) == RES_OK);
+ test_scene_3d(dev, interf);
+ test_scene_2d(dev, interf);
+
+ CHK(sdis_device_ref_put(dev) == RES_OK);
+ CHK(sdis_interface_ref_put(interf) == RES_OK);
+
check_memory_allocator(&allocator);
mem_shutdown_proxy_allocator(&allocator);
CHK(mem_allocated_size() == 0);
diff --git a/src/test_sdis_solve_probe_boundary.c b/src/test_sdis_solve_probe_boundary.c
@@ -355,10 +355,6 @@ main(int argc, char** argv)
CHK(sdis_estimator_get_temperature(estimator, &T) == RES_OK);
CHK(sdis_estimator_ref_put(estimator) == RES_OK);
- CHK(sdis_scene_get_boundary_position(NULL, iprim, uv, pos) == RES_BAD_ARG);
- CHK(sdis_scene_get_boundary_position(box_scn, 12, uv, pos) == RES_BAD_ARG);
- CHK(sdis_scene_get_boundary_position(box_scn, iprim, NULL, pos) == RES_BAD_ARG);
- CHK(sdis_scene_get_boundary_position(box_scn, iprim, uv, NULL) == RES_BAD_ARG);
CHK(sdis_scene_get_boundary_position(box_scn, iprim, uv, pos) == RES_OK);
ref = (H*Tf + LAMBDA * Tb) / (H + LAMBDA);
