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 9091cf25808a0c164c1ea2555ddf5eef895f985f
parent be3a7a49dedfc98ef4fb3d42837a715e2ab72257
Author: Christophe Coustet <christophe.coustet@meso-star.com>
Date:   Tue, 12 Jan 2021 15:57:52 +0100

Improve closest_point accuracy

Diffstat:
Msrc/s3d_scene_view_closest_point.c | 160+++++++++++++++++++++++++++++++++++++++++--------------------------------------
1 file changed, 83 insertions(+), 77 deletions(-)

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]; - float E0[3], E1[3], Ng[3]; - float vec[3]; - float uv[2]; + 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]; + double E0[3], E1[3], Ng[3]; + 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) { @@ -235,23 +239,23 @@ 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); + d3_sub(E0, v2, v0); + d3_sub(E1, v1, v0); + d3_cross(Ng, E0, E1); /* Flip the geometric normal wrt the flip surface flag */ flip_surface ^= geom->flip_surface; - if(flip_surface) f3_minus(Ng, Ng); + if(flip_surface) d3_minus(Ng, Ng); /* Setup the hit */ hit.prim.shape__ = geom; hit.prim.inst__ = inst; hit.distance = dst; - hit.uv[0] = uv[0]; - hit.uv[1] = uv[1]; - hit.normal[0] = Ng[0]; - hit.normal[1] = Ng[1]; - hit.normal[2] = Ng[2]; + hit.uv[0] = (float)uv[0]; + hit.uv[1] = (float)uv[1]; + hit.normal[0] = (float)Ng[0]; + hit.normal[1] = (float)Ng[1]; + hit.normal[2] = (float)Ng[2]; hit.prim.prim_id = args->primID; hit.prim.geom_id = geom->name; hit.prim.inst_id = inst ? inst->name : S3D_INVALID_ID; @@ -260,12 +264,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! */ }