star-3d

Surface structuring for efficient 3D geometric queries
git clone git://git.meso-star.fr/star-3d.git
Log | Files | Refs | README | LICENSE

commit ae294c92192df43eb5441bc2df6bfdc3a8674f47
parent dff27e2ea9869abde08d775cd52772c9c319e2d6
Author: Vincent Forest <vincent.forest@meso-star.com>
Date:   Tue, 15 Jun 2021 10:12:42 +0200

Merge branch 'release_0.7.4'

Diffstat:
MREADME.md | 6++++++
Mcmake/CMakeLists.txt | 2+-
Msrc/s3d.h | 4++--
Msrc/s3d_device.c | 6++----
Msrc/s3d_scene_view_closest_point.c | 155++++++++++++++++++++++++++++++++++++++++++-------------------------------------
Msrc/s3d_scene_view_trace_ray.c | 8++++----
Msrc/test_s3d_accel_struct_conf.c | 6++++--
Msrc/test_s3d_closest_point.c | 131++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++-
Msrc/test_s3d_seams.c | 2+-
Msrc/test_s3d_trace_ray.c | 120+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++--
10 files changed, 350 insertions(+), 90 deletions(-)

diff --git a/README.md b/README.md @@ -120,6 +120,12 @@ with `<STAR3D_INSTALL_DIR>` the install directory of Star-3D and ## Release notes +### Version 0.7.4 + +- Fix the barycentric coordinates of the intersection of the ray onto the + triangle: the coordinates could lie outside the triangle. +- Fix compilation warnings. + ### Version 0.7.3 - Fix the `s3d_scene_view_closest_point` function on the scenes containing diff --git a/cmake/CMakeLists.txt b/cmake/CMakeLists.txt @@ -58,7 +58,7 @@ endif() ################################################################################ set(VERSION_MAJOR 0) set(VERSION_MINOR 7) -set(VERSION_PATCH 3) +set(VERSION_PATCH 4) set(VERSION ${VERSION_MAJOR}.${VERSION_MINOR}.${VERSION_PATCH}) set(S3D_FILES_SRC diff --git a/src/s3d.h b/src/s3d.h @@ -144,7 +144,7 @@ static const struct s3d_vertex_data S3D_VERTEX_DATA_NULL = S3D_VERTEX_DATA_NULL_ /* Intersection point */ struct s3d_hit { struct s3d_primitive prim; /* Intersected primitive */ - float normal[3]; /* Un-normalized geometry normal */ + float normal[3]; /* Un-normalized geometry normal (left hand convention) */ float uv[2]; /* Barycentric coordinates of the hit onto `prim' */ float distance; /* Hit distance from the query origin */ }; @@ -382,7 +382,7 @@ S3D_API res_T s3d_scene_view_closest_point (struct s3d_scene_view* scnview, const float pos[3], /* Position to query */ - const float radius, /* Search distance in ]0, radius[ */ + const float radius, /* Search distance in [0, radius[ */ void* query_data, /* User data sent to the hit filter func. May be NULL */ struct s3d_hit* hit); diff --git a/src/s3d_device.c b/src/s3d_device.c @@ -45,10 +45,8 @@ static INLINE void rtc_error_func(void* context, enum RTCError err, const char* str) { - char msg[128]; - (void)str, (void)err, (void)context; - snprintf(msg, sizeof(msg), "Embree:error: %s\n", rtc_error_string(err)); - FATAL(msg); + (void)str, (void)context; + VFATAL("Embree:error: %s\n", ARG1(rtc_error_string(err))); } static INLINE void diff --git a/src/s3d_scene_view_closest_point.c b/src/s3d_scene_view_closest_point.c @@ -38,7 +38,10 @@ #include "s3d_scene_view_c.h" #include "s3d_sphere.h" +#include <rsys/float2.h> #include <rsys/float3.h> +#include <rsys/double2.h> +#include <rsys/double3.h> #include <rsys/float33.h> struct point_query_context { @@ -50,103 +53,103 @@ struct point_query_context { /******************************************************************************* * Helper functions ******************************************************************************/ -static INLINE float* +static INLINE double* closest_point_triangle - (const float p[3], /* Point */ - const float a[3], /* 1st triangle vertex */ - const float b[3], /* 2nd triangle vertex */ - const float c[3], /* 3rd triangle vertex */ - float closest_pt[3], /* Closest position */ - float uv[2]) /* UV of the closest position */ + (const double p[3], /* Point */ + const double a[3], /* 1st triangle vertex */ + const double b[3], /* 2nd triangle vertex */ + const double c[3], /* 3rd triangle vertex */ + double closest_pt[3], /* Closest position */ + double uv[2]) /* UV of the closest position */ { - float ab[3], ac[3], ap[3], bp[3], cp[3]; - float d1, d2, d3, d4, d5, d6; - float va, vb, vc; - float rcp_triangle_area; - float v, w; + double ab[3], ac[3], ap[3], bp[3], cp[3]; + double d1, d2, d3, d4, d5, d6; + double va, vb, vc; + double rcp_triangle_area; + double v, w; ASSERT(p && a && b && c && closest_pt && uv); - f3_sub(ab, b, a); - f3_sub(ac, c, a); + d3_sub(ab, b, a); + d3_sub(ac, c, a); /* Check if the closest point is the triangle vertex 'a' */ - f3_sub(ap, p, a); - d1 = f3_dot(ab, ap); - d2 = f3_dot(ac, ap); - if(d1 <= 0.f && d2 <= 0.f) { - uv[0] = 1.f; - uv[1] = 0.f; - return f3_set(closest_pt, a); + d3_sub(ap, p, a); + d1 = d3_dot(ab, ap); + d2 = d3_dot(ac, ap); + if(d1 <= 0 && d2 <= 0) { + uv[0] = 1; + uv[1] = 0; + return d3_set(closest_pt, a); } /* Check if the closest point is the triangle vertex 'b' */ - f3_sub(bp, p, b); - d3 = f3_dot(ab, bp); - d4 = f3_dot(ac, bp); - if(d3 >= 0.f && d4 <= d3) { - uv[0] = 0.f; - uv[1] = 1.f; - return f3_set(closest_pt, b); + d3_sub(bp, p, b); + d3 = d3_dot(ab, bp); + d4 = d3_dot(ac, bp); + if(d3 >= 0 && d4 <= d3) { + uv[0] = 0; + uv[1] = 1; + return d3_set(closest_pt, b); } /* Check if the closest point is the triangle vertex 'c' */ - f3_sub(cp, p, c); - d5 = f3_dot(ab, cp); - d6 = f3_dot(ac, cp); - if(d6 >= 0.f && d5 <= d6) { - uv[0] = 0.f; - uv[1] = 0.f; - return f3_set(closest_pt, c); + d3_sub(cp, p, c); + d5 = d3_dot(ab, cp); + d6 = d3_dot(ac, cp); + if(d6 >= 0 && d5 <= d6) { + uv[0] = 0; + uv[1] = 0; + return d3_set(closest_pt, c); } /* Check if the closest point is on the triangle edge 'ab' */ vc = d1*d4 - d3*d2; - if(vc <= 0.f && d1 >= 0.f && d3 <= 0.f) { - const float s = d1 / (d1 - d3); + if(vc <= 0 && d1 >= 0 && d3 <= 0) { + const double s = d1 / (d1 - d3); closest_pt[0] = a[0] + s*ab[0]; closest_pt[1] = a[1] + s*ab[1]; closest_pt[2] = a[2] + s*ab[2]; - uv[0] = 1.f - s; + uv[0] = 1 - s; uv[1] = s; return closest_pt; } /* Check if the closest point is on the triangle edge 'ac' */ vb = d5*d2 - d1*d6; - if(vb <= 0.f && d2 >= 0.f && d6 <= 0.f) { - const float s = d2 / (d2 - d6); + if(vb <= 0 && d2 >= 0 && d6 <= 0) { + const double s = d2 / (d2 - d6); closest_pt[0] = a[0] + s*ac[0]; closest_pt[1] = a[1] + s*ac[1]; closest_pt[2] = a[2] + s*ac[2]; - uv[0] = 1.f - s; - uv[1] = 0.f; + uv[0] = 1 - s; + uv[1] = 0; return closest_pt; } /* Check if the closest point is on the triangle edge 'bc' */ va = d3*d6 - d5*d4; - if(va <= 0.f && (d4 - d3) >= 0.f && (d5 - d6) >= 0.f) { - const float s = (d4 - d3) / ((d4 - d3) + (d5 - d6)); + if(va <= 0 && (d4 - d3) >= 0 && (d5 - d6) >= 0) { + const double s = (d4 - d3) / ((d4 - d3) + (d5 - d6)); closest_pt[0] = b[0] + s*(c[0] - b[0]); closest_pt[1] = b[1] + s*(c[1] - b[1]); closest_pt[2] = b[2] + s*(c[2] - b[2]); - uv[0] = 0.f; - uv[1] = 1.f - s; + uv[0] = 0; + uv[1] = 1 - s; return closest_pt; } /* The closest point lies in the triangle: compute its barycentric * coordinates */ - rcp_triangle_area = 1.f / (va + vb + vc); + rcp_triangle_area = 1 / (va + vb + vc); v = vb * rcp_triangle_area; w = vc * rcp_triangle_area; /* Save the uv barycentric coordinates */ - uv[0] = 1.f - v - w; + uv[0] = 1 - v - w; uv[1] = v; - ASSERT(eq_epsf(uv[0] + uv[1] + w, 1.f, 1.e-4f)); + ASSERT(eq_eps(uv[0] + uv[1] + w, 1, 1.e-4)); - if(uv[0] < 0.f) { /* Handle precision issues */ + if(uv[0] < 0) { /* Handle precision issues */ if(uv[1] > w) uv[1] += uv[0]; uv[0] = 0; } @@ -170,14 +173,15 @@ closest_point_mesh struct s3d_hit* out_hit = NULL; struct hit_filter* filter = NULL; const uint32_t* ids = NULL; - float closest_point[3]; - float query_pos_ws[3]; /* World space query position */ - float query_pos_ls[3]; /* Local space query position */ - float v0[3], v1[3], v2[3]; + double closest_point[3]; + double query_pos_ws[3]; /* World space query position */ + double query_pos_ls[3]; /* Local space query position */ + double v0[3], v1[3], v2[3]; float E0[3], E1[3], Ng[3]; - float vec[3]; - float uv[2]; + double vec[3]; + double uv[2]; float dst; + float pos[3], dir[3]; int flip_surface = 0; ASSERT(args && geom && geom->type == GEOM_MESH); ASSERT(args->primID < mesh_get_ntris(geom->data.mesh)); @@ -187,9 +191,9 @@ closest_point_mesh /* Fetch triangle vertices */ ASSERT(geom->data.mesh->attribs_type[S3D_POSITION] == S3D_FLOAT3); - f3_set(v0, mesh_get_pos(geom->data.mesh) + ids[0]*3/*#coords per vertex*/); - f3_set(v1, mesh_get_pos(geom->data.mesh) + ids[1]*3/*#coords per vertex*/); - f3_set(v2, mesh_get_pos(geom->data.mesh) + ids[2]*3/*#coords per vertex*/); + d3_set_f3(v0, mesh_get_pos(geom->data.mesh) + ids[0]*3/*#coords per vertex*/); + d3_set_f3(v1, mesh_get_pos(geom->data.mesh) + ids[1]*3/*#coords per vertex*/); + d3_set_f3(v2, mesh_get_pos(geom->data.mesh) + ids[2]*3/*#coords per vertex*/); /* Local copy of the query position */ query_pos_ws[0] = args->query->x; @@ -202,16 +206,16 @@ closest_point_mesh query_pos_ls[2] = query_pos_ws[2]; } else { const float* world2inst; - float a[3], b[3], c[3]; + double a[3], b[3], c[3], tmp[3]; ASSERT(args->context->instStackSize == 1); ASSERT(inst && inst->type == GEOM_INSTANCE); world2inst = args->context->world2inst[0]; /* Transform the query position in instance space */ - f3_mulf(a, world2inst+0, query_pos_ws[0]); - f3_mulf(b, world2inst+4, query_pos_ws[1]); - f3_mulf(c, world2inst+8, query_pos_ws[2]); + d3_muld(a, d3_set_f3(tmp, world2inst+0), query_pos_ws[0]); + d3_muld(b, d3_set_f3(tmp, world2inst+4), query_pos_ws[1]); + d3_muld(c, d3_set_f3(tmp, world2inst+8), query_pos_ws[2]); query_pos_ls[0] = a[0] + b[0] + c[0] + world2inst[12]; query_pos_ls[1] = a[1] + b[1] + c[1] + world2inst[13]; query_pos_ls[2] = a[2] + b[2] + c[2] + world2inst[14]; @@ -223,8 +227,8 @@ closest_point_mesh closest_point_triangle(query_pos_ls, v0, v1, v2, closest_point, uv); /* Compute the distance */ - f3_sub(vec, closest_point, query_pos_ls); - dst = f3_len(vec); + d3_sub(vec, closest_point, query_pos_ls); + dst = (float)d3_len(vec); /* Transform the distance in world space */ if(args->context->instStackSize != 0) { @@ -234,10 +238,13 @@ closest_point_mesh if(dst >= args->query->radius) return 0; - /* Compute the triangle normal in world space */ - f3_sub(E0, v2, v0); - f3_sub(E1, v1, v0); - f3_cross(Ng, E0, E1); + /* Compute the triangle normal in world space (left hand convention). Keep + * it in float to avoid double-cast accuracy loss wrt user computed result */ + f3_sub(E0, mesh_get_pos(geom->data.mesh) + ids[1] * 3, + mesh_get_pos(geom->data.mesh) + ids[0] * 3); + f3_sub(E1, mesh_get_pos(geom->data.mesh) + ids[2] * 3, + mesh_get_pos(geom->data.mesh) + ids[0] * 3); + f3_cross(Ng, E1, E0); /* Flip the geometric normal wrt the flip surface flag */ flip_surface ^= geom->flip_surface; @@ -247,8 +254,8 @@ closest_point_mesh hit.prim.shape__ = geom; hit.prim.inst__ = inst; hit.distance = dst; - hit.uv[0] = uv[0]; - hit.uv[1] = uv[1]; + hit.uv[0] = (float)uv[0]; + hit.uv[1] = (float)uv[1]; hit.normal[0] = Ng[0]; hit.normal[1] = Ng[1]; hit.normal[2] = Ng[2]; @@ -260,12 +267,14 @@ closest_point_mesh + geom->scene_prim_id_offset + (inst ? inst->scene_prim_id_offset : 0); - /* `vec' is the direction along which the closest point was found. We thus + /* `dir' is the direction along which the closest point was found. We thus * submit it to the filter function as the direction corresponding to the * computed hit */ + f3_set_d3(dir, vec); + f3_set_d3(pos, query_pos_ws); filter = &geom->data.mesh->filter; if(filter->func - && filter->func(&hit, query_pos_ws, vec, query_data, filter->data)) { + && filter->func(&hit, pos, dir, query_data, filter->data)) { return 0; /* This point is filtered. Discard it! */ } diff --git a/src/s3d_scene_view_trace_ray.c b/src/s3d_scene_view_trace_ray.c @@ -116,13 +116,13 @@ hit_setup ASSERT(hit->prim.shape__); ASSERT(((struct geometry*)hit->prim.shape__)->type == GEOM_MESH ||((struct geometry*)hit->prim.shape__)->type == GEOM_SPHERE); - - hit->uv[0] = ray_hit->hit.u; - hit->uv[1] = ray_hit->hit.v; + + /* Handle Embree returning uv out of range */ + hit->uv[0] = CLAMP(ray_hit->hit.u, 0, 1); + hit->uv[1] = CLAMP(ray_hit->hit.v, 0, 1); if(((struct geometry*)hit->prim.shape__)->type == GEOM_MESH) { w = 1.f - hit->uv[0] - hit->uv[1]; - ASSERT(w <= 1.f); /* This may not occurs */ if(w < 0.f) { /* Handle precision error */ if(hit->uv[0] > hit->uv[1]) hit->uv[0] += w; else hit->uv[1] += w; diff --git a/src/test_s3d_accel_struct_conf.c b/src/test_s3d_accel_struct_conf.c @@ -67,8 +67,10 @@ mesh_init_sphere sphere->ntris = 2*nthetas * (nphis-2)/*#contour tris*/ + 2*nthetas/*#polar tris*/; sphere->allocator = allocator; - CHK(sphere->pos = MEM_CALLOC(allocator, sphere->nverts, sizeof(double[3]))); - CHK(sphere->ids = MEM_CALLOC(allocator, sphere->ntris, sizeof(size_t[3]))); + sphere->pos = MEM_CALLOC(allocator, sphere->nverts, sizeof(double[3])); + CHK(sphere->pos); + sphere->ids = MEM_CALLOC(allocator, sphere->ntris, sizeof(size_t[3])); + CHK(sphere->ids); /* Build the contour vertices */ i = 0; diff --git a/src/test_s3d_closest_point.c b/src/test_s3d_closest_point.c @@ -34,6 +34,7 @@ #include "test_s3d_cbox.h" #include "test_s3d_utils.h" +#include <rsys/float2.h> #include <rsys/float3.h> #include <rsys/float33.h> #include <limits.h> @@ -853,7 +854,7 @@ test_single_triangle(struct s3d_device* dev) float closest_pos[3] = {0,0,0}; float low[3], upp[3], mid[3]; union { float f; uint32_t ui32; } ucast; - size_t i; + size_t a, i; f3(vertices+0, -0.5f, -0.3f, 0.1f); f3(vertices+3, -0.4f, 0.2f, 0.3f); @@ -950,6 +951,134 @@ test_single_triangle(struct s3d_device* dev) CHK(s3d_shape_ref_put(msh) == RES_OK); CHK(s3d_scene_view_ref_put(view) == RES_OK); CHK(s3d_scene_ref_put(scn) == RES_OK); + + /* Check accuracy on a configuration whose analytic distance is known */ + FOR_EACH(a, 0, 16) { + const float amplitude = exp2f((float)a); + const float eps = 5e-6f * amplitude; + FOR_EACH(i, 0, 1000) { + float A[3], B[3], C[3], AB[3], AC[3], BC[3], N[3], hit_N[3]; + int j, n; + + /* Randomly generate a triangle ABC */ + FOR_EACH(n, 0, 3) + A[n] = (rand_canonic() - 0.5f) * amplitude; + do { + FOR_EACH(n, 0, 3) B[n] = (rand_canonic() - 0.5f) * amplitude; + } while (f3_eq_eps(A, B, eps)); + do { + FOR_EACH(n, 0, 3) C[n] = (rand_canonic() - 0.5f) * amplitude; + } while (f3_eq_eps(A, C, eps) || f3_eq_eps(B, C, eps)); + + f3_sub(AB, B, A); + f3_sub(AC, C, A); + f3_sub(BC, C, B); + f3_cross(N, AC, AB); /* Left hand convention */ + f3_normalize(N, N); + + f3_set(vertices + 0, A); + f3_set(vertices + 3, B); + f3_set(vertices + 6, C); + + CHK(s3d_scene_create(dev, &scn) == RES_OK); + CHK(s3d_shape_create_mesh(dev, &msh) == RES_OK); + CHK(s3d_scene_attach_shape(scn, msh) == RES_OK); + + vdata.usage = S3D_POSITION; + vdata.type = S3D_FLOAT3; + vdata.get = triangle_get_pos; + CHK(s3d_mesh_setup_indexed_vertices + (msh, 1, triangle_get_ids, 3, &vdata, 1, vertices) == RES_OK); + + CHK(s3d_scene_view_create(scn, S3D_TRACE, &view) == RES_OK); + + FOR_EACH(j, 0, 1000) { + float proj[3]; /* Projection of pos on the line */ + float AP[3], BP[3], CP[3], closest[3], tmp[3]; + float u, v, w, h, x, dist, d; + + /* Randomly generate a pos not on the triangle + * with know position wrt the problem: pos = A + u.AB + v.AC + k.N */ + u = 3 * rand_canonic() - 1; + v = 3 * rand_canonic() - 1; + w = 1 - u - v; + h = (2 * rand_canonic() - 1) * amplitude; + f3_add(proj, A, f3_add(proj, f3_mulf(proj, AB, u), f3_mulf(tmp, AC, v))); + f3_add(pos, proj, f3_mulf(pos, N, h)); + f3_sub(AP, proj, A); + f3_sub(BP, proj, B); + f3_sub(CP, proj, C); + + /* Compute closest point */ + CHK(s3d_scene_view_closest_point(view, pos, (float)INF, NULL, &hit) + == RES_OK); + CHK(!S3D_HIT_NONE(&hit)); + CHK(s3d_primitive_get_attrib(&hit.prim, S3D_POSITION, hit.uv, &attr) + == RES_OK); + + /* Check result + * Due to known uv lack of accuracy we mainly check distance */ + if(u >= 0 && v >= 0 && w >= 0) { + /* proj is inside the triangle and is the closest point */ + f3_set(closest, proj); + dist = fabsf(h); + } else { + /* proj is outside the triangle */ + float lab2 = f3_dot(AB, AB); + float lac2 = f3_dot(AC, AC); + float lbc2 = f3_dot(BC, BC); + if(w >= 0 && u < 0) { + /* proj is closest to either AB or AC */ + x = f3_dot(AP, AB); + if(v < 0 && x > 0) { + /* proj is closest to AB */ + f3_add(closest, A, f3_mulf(tmp, AB, MMIN(1, x / lab2))); + } else { + /* proj is closest to AC */ + f3_add(closest, A, + f3_mulf(tmp, AC, MMIN(1, MMAX(0, f3_dot(AP, AC) / lac2)))); + } + } + else if(u >= 0 && v < 0) { + /* proj is closest to either BC or BA */ + x = f3_dot(BP, BC); + if(w < 0 && x > 0) { + /* proj is closest to BC */ + f3_add(closest, B, f3_mulf(tmp, BC, MMIN(1, x / lbc2))); + } else { + /* proj is closest to BA */ + f3_add(closest, B, + f3_mulf(tmp, AB, -MMIN(1, MMAX(0, -f3_dot(BP, AB) / lab2)))); + } + } + else if(v >= 0 && w < 0) { + /* proj is closest to either CA or CB */ + x = -f3_dot(CP, AC); + if(u < 0 && x > 0) { + /* proj is closest to CA */ + f3_add(closest, C, f3_mulf(tmp, AC, -MMIN(1, x / lac2))); + } else { + /* proj is closest to CB */ + f3_add(closest, C, + f3_mulf(tmp, BC, -MMIN(1, MMAX(0, -f3_dot(CP, BC) / lbc2)))); + } + } + else { ASSERT(0); } + dist = f3_len(f3_sub(tmp, pos, closest)); + } + CHK(eq_epsf(hit.distance, dist, eps)); + /* Intersection-point's position is less accurate than hit distance */ + d = f3_len(f3_sub(tmp, closest, attr.value)); + CHK(d <= 10 * eps); + f3_normalize(hit_N, hit.normal); + CHK(f3_eq_eps(N, hit_N, FLT_EPSILON)); + } + + CHK(s3d_shape_ref_put(msh) == RES_OK); + CHK(s3d_scene_view_ref_put(view) == RES_OK); + CHK(s3d_scene_ref_put(scn) == RES_OK); + } + } } static void diff --git a/src/test_s3d_seams.c b/src/test_s3d_seams.c @@ -161,7 +161,7 @@ check_ray(int use_double) CHK(s3d_scene_view_create(scn2, S3D_TRACE, &scnview) == RES_OK); CHK(s3d_scene_view_trace_ray(scnview, org, dir, range, NULL, &hit) == RES_OK); printf("\nRaytrace using %s: ", use_double ? "double" : "float"); - if (!S3D_HIT_NONE(&hit)) { + if(!S3D_HIT_NONE(&hit)) { f3_add(pos, org, f3_mulf(pos, dir, hit.distance)); printf("Hit at [%g %g %g]\n",SPLIT3(pos)); } else { diff --git a/src/test_s3d_trace_ray.c b/src/test_s3d_trace_ray.c @@ -60,6 +60,32 @@ filter_func return hit->prim.prim_id % 2 == 0; } +static void +triangle_get_ids(const unsigned itri, unsigned ids[3], void* ctx) +{ + (void)ctx; + CHK(itri == 0); + CHK(ids); + ids[0] = 0; + ids[1] = 1; + ids[2] = 2; +} + +static void +triangle_get_pos(const unsigned ivert, float pos[3], void* ctx) +{ + float* vertices = ctx; + CHK(ctx); + CHK(ivert < 3); + CHK(pos); + switch (ivert) { /* Setup a random triangle */ + case 0: f3_set(pos, vertices + 0); break; + case 1: f3_set(pos, vertices + 3); break; + case 2: f3_set(pos, vertices + 6); break; + default: FATAL("Unreachable code\n"); break; + } +} + int main(int argc, char** argv) { @@ -93,7 +119,7 @@ main(int argc, char** argv) unsigned walls_id; unsigned tall_block_id; unsigned short_block_id; - size_t i; + size_t a, i; char filter = 0; mem_init_proxy_allocator(&allocator, &mem_default_allocator); @@ -354,7 +380,6 @@ main(int argc, char** argv) CHK(image_write_ppm_stream(&img, 0, stdout) == RES_OK); image_release(&img); - CHK(s3d_device_ref_put(dev) == RES_OK); CHK(s3d_shape_ref_put(inst) == RES_OK); CHK(s3d_shape_ref_put(short_block) == RES_OK); CHK(s3d_shape_ref_put(tall_block) == RES_OK); @@ -362,6 +387,97 @@ main(int argc, char** argv) CHK(s3d_scene_ref_put(scn) == RES_OK); CHK(s3d_scene_ref_put(scn2) == RES_OK); + /* Check accuracy on a configuration whose analytic distance is known */ + FOR_EACH(a, 0, 16) { + const float amplitude = exp2f((float)a); + const float eps = 5e-6f * amplitude; + float vertices[9]; + struct s3d_vertex_data vdata = S3D_VERTEX_DATA_NULL; + struct s3d_scene_view* view = NULL; + struct s3d_shape* msh = NULL; + FOR_EACH(i, 0, 1000) { + float A[3], B[3], C[3], AB[3], AC[3], N[3]; + int j, n; + + /* Randomly generate a triangle ABC */ + FOR_EACH(n, 0, 3) + A[n] = (rand_canonic() - 0.5f) * amplitude; + do { + FOR_EACH(n, 0, 3) B[n] = (rand_canonic() - 0.5f) * amplitude; + } while (f3_eq_eps(A, B, eps)); + do { + FOR_EACH(n, 0, 3) C[n] = (rand_canonic() - 0.5f) * amplitude; + } while (f3_eq_eps(A, C, eps) || f3_eq_eps(B, C, eps)); + + f3_sub(AB, B, A); + f3_sub(AC, C, A); + f3_cross(N, AB, AC); + f3_normalize(N, N); + + f3_set(vertices + 0, A); + f3_set(vertices + 3, B); + f3_set(vertices + 6, C); + + CHK(s3d_scene_create(dev, &scn) == RES_OK); + CHK(s3d_shape_create_mesh(dev, &msh) == RES_OK); + CHK(s3d_scene_attach_shape(scn, msh) == RES_OK); + + vdata.usage = S3D_POSITION; + vdata.type = S3D_FLOAT3; + vdata.get = triangle_get_pos; + CHK(s3d_mesh_setup_indexed_vertices + (msh, 1, triangle_get_ids, 3, &vdata, 1, vertices) == RES_OK); + + CHK(s3d_scene_view_create(scn, S3D_TRACE, &view) == RES_OK); + + FOR_EACH(j, 0, 1000) { + float proj[3]; /* Projection of pos on the line */ + float tmp[3]; + float u, v, w, h; + + /* Randomly generate a pos not on the triangle + * with know position wrt the problem: pos = A + u.AB + v.AC + k.N */ + u = 3 * rand_canonic() - 1; + v = 3 * rand_canonic() - 1; + w = 1 - u - v; + h = (2 * rand_canonic() - 1) * amplitude; + f3_add(proj, A, f3_add(proj, f3_mulf(proj, AB, u), f3_mulf(tmp, AC, v))); + f3_add(pos, proj, f3_mulf(pos, N, h)); + + /* Raytrace from pos towards proj */ + f3_mulf(dir, N, (h > 0 ? -1.f : 1.f)); + f3_normalize(dir, dir); + CHK(s3d_scene_view_trace_ray(view, pos, dir, range, NULL, &hit) + == RES_OK); + + /* Check result */ + if(u < 0 || v < 0 || w < 0) { + if(!S3D_HIT_NONE(&hit)) + CHK(u >= -FLT_EPSILON && v >= -FLT_EPSILON && w >= -FLT_EPSILON); + } else { + if(S3D_HIT_NONE(&hit)) + CHK(u <= FLT_EPSILON || v <= FLT_EPSILON || w <= FLT_EPSILON); + } + if(!S3D_HIT_NONE(&hit)) { + struct s3d_attrib attr; + float d; + CHK(eq_epsf(hit.distance, fabsf(h), eps)); + CHK(s3d_primitive_get_attrib(&hit.prim, S3D_POSITION, hit.uv, &attr) + == RES_OK); + /* Intersection-point's position is less accurate than hit distance */ + d = f3_len(f3_sub(tmp, attr.value, proj)); + CHK(d <= 10 * eps); + } + } + + CHK(s3d_shape_ref_put(msh) == RES_OK); + CHK(s3d_scene_view_ref_put(view) == RES_OK); + CHK(s3d_scene_ref_put(scn) == RES_OK); + } + } + + CHK(s3d_device_ref_put(dev) == RES_OK); + check_memory_allocator(&allocator); mem_shutdown_proxy_allocator(&allocator); CHK(mem_allocated_size() == 0);