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 a8d1e50b58b70de5b71e3aaefb6288d1b2102d6c
parent 8e782af7a7bef563c739a67806fec29bb643d92a
Author: Vincent Forest <vincent.forest@meso-star.com>
Date:   Mon, 30 Sep 2019 11:56:28 +0200

Add the s3d_scene_view_point_query function

Diffstat:
Mcmake/CMakeLists.txt | 1+
Msrc/s3d.h | 29++++++++++++++++++++++-------
Msrc/s3d_geometry.c | 20++++----------------
Msrc/s3d_geometry.h | 2+-
Asrc/s3d_scene_view_point_query.c | 439+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Msrc/s3d_sphere.h | 22++++++++++++++++++++++
6 files changed, 489 insertions(+), 24 deletions(-)

diff --git a/cmake/CMakeLists.txt b/cmake/CMakeLists.txt @@ -69,6 +69,7 @@ set(S3D_FILES_SRC s3d_primitive.c s3d_scene.c s3d_scene_view.c + s3d_scene_view_point_query.c s3d_shape.c s3d_sphere.c) set(S3D_FILES_INC_API s3d.h) diff --git a/src/s3d.h b/src/s3d.h @@ -146,7 +146,7 @@ struct s3d_hit { struct s3d_primitive prim; /* Intersected primitive */ float normal[3]; /* Unnormalized geometry normal */ float uv[2]; /* Barycentric coordinates of the hit onto `prim' */ - float distance; /* Hit distance from the ray origin */ + float distance; /* Hit distance from the query origin */ }; /* Constant defining a NULL intersection. Should be used to initialize a hit */ @@ -295,6 +295,11 @@ s3d_scene_view_get_mask (struct s3d_scene_view* scnview, int* mask); +/* Trace a ray into the scene and return the closest intersection along it. The + * ray is defined by `origin' + t*`direction' = 0 with t in [`range[0]', + * `range[1]'). Note that if a range is degenerated (i.e. `range[0]' >= + * `range[1]') then the ray is not traced and `hit' is set to S3D_HIT_NULL. Can + * be called only if the scnview was created with the S3D_TRACE flag. */ S3D_API res_T s3d_scene_view_trace_ray (struct s3d_scene_view* scnview, @@ -304,12 +309,8 @@ s3d_scene_view_trace_ray void* ray_data, /* User ray data sent to the hit filter func. May be NULL */ struct s3d_hit* hit); -/* Trace a bundle of rays into the scene and return the closest intersection - * along them. The rays are defined by `origin' + t*`direction' = 0 with t in - * [`range[0]', `range[1]'). Note that if a range is degenerated (i.e. - * `range[0]' >= `range[1]') then its associated ray is not traced and `hit' is - * set to S3D_HIT_NULL. Can be called only if the scnview was created with the - * S3D_TRACE flag. */ +/* Trace a bundle of rays into the scene. Can be called only if the scnview was + * created with the S3D_TRACE flag. */ S3D_API res_T s3d_scene_view_trace_rays (struct s3d_scene_view* scnview, @@ -322,6 +323,20 @@ s3d_scene_view_trace_rays const size_t sizeof_ray_data, /* Size in Bytes of *one* ray data */ struct s3d_hit* hits); +/* Return the point onto the scene surfaces that is the closest of the + * submitted `pos'. The `radius' parameter defines the maximum search distance + * around `pos'. Each candidate position are internally filtered by the + * hit_filter_function attached to the corresponding shape; the user can thus + * reject a candidate position according to its own criteria. Cab ne called + * only if the scenview was created with the S3D_TRACE flag */ +S3D_API res_T +s3d_scene_view_point_query + (struct s3d_scene_view* scnview, + const float pos[3], /* Position to query */ + const float radius, /* Maxium search distance around pos */ + void* query_data, /* User data sent to the hit filter func. May be NULL */ + struct s3d_hit* hit); + /* Uniformly sample the scene and return the sampled primitive and its sample * uv position. Can be called only if the scnview was created with the * S3D_SAMPLE flag */ diff --git a/src/s3d_geometry.c b/src/s3d_geometry.c @@ -51,9 +51,8 @@ sphere_ray_hit_setup struct RTCHitN* hitN; struct RTCHit hit; struct RTCRay ray; - float cos_theta; float Ng[3]; - float u, v; + float uv[2]; size_t i; ASSERT(args && args->primID == 0 && args->N == 1 && args->valid[0] != 0); @@ -71,24 +70,13 @@ sphere_ray_hit_setup Ng[2] = ray.dir_z*tfar + ray.org_z - geom->data.sphere->pos[2]; f3_normalize(Ng, Ng); - cos_theta = Ng[2]; - - v = (1.f - cos_theta) * 0.5f; - if(absf(cos_theta) == 1) { - u = 0; - } else if(eq_epsf(Ng[0], 0.f, 1.e-6f)) { - u = Ng[1] > 0 ? 0.25f : 0.75f; - } else { - double phi = atan2f(Ng[1], Ng[0]); /* phi in [-PI, PI] */ - if(phi < 0) phi = 2*PI + phi; /* phi in [0, 2PI] */ - u = (float)(phi / (2*PI)); - } + sphere_normal_to_uv(Ng, uv); hit.Ng_x = Ng[0]; hit.Ng_y = Ng[1]; hit.Ng_z = Ng[2]; - hit.u = u; - hit.v = v; + hit.u = uv[0]; + hit.v = uv[1]; hit.primID = 0; hit.geomID = geom->rtc_id; FOR_EACH(i, 0, RTC_MAX_INSTANCE_LEVEL_COUNT) { diff --git a/src/s3d_geometry.h b/src/s3d_geometry.h @@ -35,6 +35,7 @@ #include "s3d.h" #include "s3d_backend.h" +#include <rsys/float3.h> #include <rsys/ref_count.h> enum geometry_type { @@ -99,4 +100,3 @@ geometry_rtc_sphere_intersect (const struct RTCIntersectFunctionNArguments* args); #endif /* S3D_GEOMETRY_H */ - diff --git a/src/s3d_scene_view_point_query.c b/src/s3d_scene_view_point_query.c @@ -0,0 +1,439 @@ +/* Copyright (C) 2015-2019 |Meso|Star> (contact@meso-star.com) + * + * This software is a computer program whose purpose is to describe a + * virtual 3D environment that can be ray-traced and sampled both robustly + * and efficiently. + * + * This software is governed by the CeCILL license under French law and + * abiding by the rules of distribution of free software. You can use, + * modify and/or redistribute the software under the terms of the CeCILL + * license as circulated by CEA, CNRS and INRIA at the following URL + * "http://www.cecill.info". + * + * As a counterpart to the access to the source code and rights to copy, + * modify and redistribute granted by the license, users are provided only + * with a limited warranty and the software's author, the holder of the + * economic rights, and the successive licensors have only limited + * liability. + * + * In this respect, the user's attention is drawn to the risks associated + * with loading, using, modifying and/or developing or reproducing the + * software by the user in light of its specific status of free software, + * that may mean that it is complicated to manipulate, and that also + * therefore means that it is reserved for developers and experienced + * professionals having in-depth computer knowledge. Users are therefore + * encouraged to load and test the software's suitability as regards their + * requirements in conditions enabling the security of their systems and/or + * data to be ensured and, more generally, to use and operate it in the + * same conditions as regards security. + * + * The fact that you are presently reading this means that you have had + * knowledge of the CeCILL license and that you accept its terms. */ + +#include "s3d.h" +#include "s3d_device_c.h" +#include "s3d_instance.h" +#include "s3d_geometry.h" +#include "s3d_mesh.h" +#include "s3d_scene_view_c.h" +#include "s3d_sphere.h" + +#include <rsys/float3.h> +#include <rsys/float33.h> + +struct point_query_context { + struct RTCPointQueryContext rtc; + struct s3d_scene_view* scnview; + void* data; /* Per point query defined data */ +}; + +/******************************************************************************* + * Helper functions + ******************************************************************************/ +static INLINE float* +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 */ +{ + float ab[3], ac[3], ap[3], bp[3], cp[3]; + float d1, d2, d3, d4, d5, d6; + float va, vb, vc; + float triangle_area; + float v, w; + ASSERT(p && a && b && c && closest_pt && uv); + + f3_sub(ab, b, a); + f3_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); + } + + /* 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); + } + + /* 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); + } + + /* 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); + 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[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); + 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; + 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)); + 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; + return closest_pt; + } + + /* The closest point lies in the triangle: compute its barycentric + * coordinates */ + triangle_area = va + vb + vc; + v = vb * triangle_area; + w = vc * triangle_area; + + /* Save the uv barycentric coordinates */ + uv[0] = 1.f - v - w; + uv[1] = v; + + /* Use the barycentric coordinates to compute the world space position of the + * closest point onto the triangle */ + closest_pt[0] = a[0] + v*ab[0] + w*ac[0]; + closest_pt[1] = a[1] + v*ab[1] + w*ac[1]; + closest_pt[2] = a[2] + v*ab[2] + w*ac[2]; + return closest_pt; +} + +static bool +closest_point_mesh + (struct RTCPointQueryFunctionArguments* args, + struct geometry* geom, + struct geometry* inst, /* Can be NULL */ + void* query_data) +{ + struct s3d_hit hit = S3D_HIT_NULL; + struct s3d_hit* out_hit = NULL; + struct hit_filter* filter = NULL; + const uint32_t* ids = NULL; + float closest_point[3]; + float query_pos[3]; + float v0[3], v1[3], v2[3]; + float E0[3], E1[3], Ng[3]; + float vec[3]; + float uv[2]; + float dst; + int flip_surface = 0; + ASSERT(args && geom && geom->type == GEOM_MESH); + ASSERT(args->primID < mesh_get_ntris(geom->data.mesh)); + + /* Fetch triangle indices */ + ids = mesh_get_ids(geom->data.mesh) + args->primID*3/*#indices per triangle*/; + + /* 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*/); + + /* Local copy of the query position */ + query_pos[0] = args->query->x; + query_pos[1] = args->query->y; + query_pos[2] = args->query->z; + + if(args->context->instStackSize) { /* The mesh is instantiated */ + const float* transform; + transform = inst->data.instance->transform; + ASSERT(args->context->instStackSize == 1); + ASSERT(args->similarityScale == 1); + ASSERT(inst && inst->type == GEOM_INSTANCE); + ASSERT(f3_eq_eps(transform+0, args->context->inst2world[0]+0, 1.e-6f)); + ASSERT(f3_eq_eps(transform+3, args->context->inst2world[0]+4, 1.e-6f)); + ASSERT(f3_eq_eps(transform+6, args->context->inst2world[0]+8, 1.e-6f)); + ASSERT(f3_eq_eps(transform+9, args->context->inst2world[0]+12,1.e-6f)); + + /* Transform the triangle vertices in world space */ + f3_add(v0, f33_mulf3(v0, transform, v0), transform+9); + f3_add(v1, f33_mulf3(v1, transform, v1), transform+9); + f3_add(v2, f33_mulf3(v2, transform, v2), transform+9); + + flip_surface = inst->flip_surface; + } + + /* Compute the closest point onto the triangle from the submitted point */ + closest_point_triangle(query_pos, v0, v1, v2, closest_point, uv); + + f3_sub(vec, closest_point, query_pos); + dst = f3_len(vec); + 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); + + /* Flip the geometric normal wrt the flip surface flag */ + flip_surface ^= geom->flip_surface; + if(flip_surface) f3_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.prim.prim_id = args->primID; + hit.prim.geom_id = geom->name; + hit.prim.inst_id = inst ? inst->name : S3D_INVALID_ID; + hit.prim.scene_prim_id = + hit.prim.prim_id + + 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 + * submit it to the filter function as the direction corresponding to the + * computed hit. */ + filter = &geom->data.mesh->filter; + if(filter->func + && filter->func(&hit, query_pos, vec, query_data, filter->data)) { + return 0; /* This point is filtered. Discard it! */ + } + + /* Update output data */ + out_hit = args->userPtr; + *out_hit = hit; + + /* Shrink the query radius */ + args->query->radius = dst; + + return 1; /* Notify that the query radius was updated */ +} + +static bool +closest_point_sphere + (struct RTCPointQueryFunctionArguments* args, + struct geometry* geom, + struct geometry* inst, + void* query_data) +{ + struct s3d_hit hit = S3D_HIT_NULL; + struct s3d_hit* out_hit = NULL; + struct hit_filter* filter = NULL; + float query_pos[3]; + float sphere_pos[3]; + float Ng[3]; + float dir[3]; + float uv[2]; + float dst; + float len; + int flip_surface = 0; + ASSERT(args && geom && geom->type == GEOM_SPHERE && args->primID == 1); + + /* Local copy of the query position */ + query_pos[0] = args->query->x; + query_pos[1] = args->query->y; + query_pos[2] = args->query->z; + + f3_set(sphere_pos, geom->data.sphere->pos); + if(args->context->instStackSize) { /* The sphere is instantiated */ + const float* transform; + transform = inst->data.instance->transform; + ASSERT(args->context->instStackSize == 1); + ASSERT(args->similarityScale == 1); + ASSERT(inst && inst->type == GEOM_INSTANCE); + ASSERT(f3_eq_eps(transform+0, args->context->inst2world[0]+0, 1.e-6f)); + ASSERT(f3_eq_eps(transform+3, args->context->inst2world[0]+4, 1.e-6f)); + ASSERT(f3_eq_eps(transform+6, args->context->inst2world[0]+8, 1.e-6f)); + ASSERT(f3_eq_eps(transform+9, args->context->inst2world[0]+12,1.e-6f)); + + /* Transform the sphere position in world space */ + f33_mulf3(sphere_pos, transform, sphere_pos); + f3_add(sphere_pos, transform+9, sphere_pos); + + flip_surface = inst->flip_surface; + } + + /* Compute the distance from the the query pos to the sphere center */ + f3_sub(Ng, query_pos, sphere_pos); + len = f3_len(Ng); + + /* Evaluate the distance from the query pos to the sphere surface */ + dst = fabsf(len - geom->data.sphere->radius); + + /* The closest point onto the sphere is outside the query radius */ + if(dst >= args->query->radius) + return 0; + + /* Compute the uv the of found point */ + f3_divf(Ng, Ng, len); /* Normalize the hit normal */ + sphere_normal_to_uv(Ng, uv); + + /* Flip the geometric normal wrt the flip surface flag */ + flip_surface ^= geom->flip_surface; + if(flip_surface) f3_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.prim.prim_id = args->primID; + hit.prim.geom_id = geom->name; + hit.prim.inst_id = inst ? inst->name : S3D_INVALID_ID; + hit.prim.scene_prim_id = + hit.prim.prim_id + + geom->scene_prim_id_offset + + (inst ? inst->scene_prim_id_offset : 0); + + + /* Use the reversed geometric normal as the hit direction since it is along + * this vector that the closest point was effectively computed */ + f3_minus(dir, Ng); + filter = &geom->data.sphere->filter; + if(filter->func + && filter->func(&hit, query_pos, dir, query_data, filter->data)) { + return 0; + } + + /* Update output data */ + out_hit = args->userPtr; + *out_hit = hit; + + /* Shrink the query radius */ + args->query->radius = dst; + + return 1; /* Notify that the query radius was updated */ +} + +static bool +closest_point(struct RTCPointQueryFunctionArguments* args) +{ + struct point_query_context* ctx = NULL; + struct geometry* geom = NULL; + struct geometry* inst = NULL; + bool query_radius_is_upd = false; + ASSERT(args); + + ctx = CONTAINER_OF(args->context, struct point_query_context, rtc); + if(args->context->instStackSize == 0) { + geom = scene_view_geometry_from_embree_id + (ctx->scnview, args->geomID); + } else { + ASSERT(args->context->instStackSize == 1); + ASSERT(args->similarityScale == 1); + inst = scene_view_geometry_from_embree_id + (ctx->scnview, args->context->instID[0]); + geom = scene_view_geometry_from_embree_id + (inst->data.instance->scnview, args->geomID); + } + + switch(geom->type) { + case GEOM_MESH: + query_radius_is_upd = closest_point_mesh(args, geom, inst, ctx->data); + break; + case GEOM_SPHERE: + query_radius_is_upd = closest_point_sphere(args, geom, inst, ctx->data); + break; + default: FATAL("Unreachable code\n"); break; + } + return query_radius_is_upd; +} + +/******************************************************************************* + * Exported functions + ******************************************************************************/ +res_T +s3d_scene_view_point_query + (struct s3d_scene_view* scnview, + const float pos[3], + const float radius, + void* query_data, + struct s3d_hit* hit) +{ + struct RTCPointQuery query; + struct point_query_context query_ctx; + + if(!scnview || !pos || radius <= 0 || hit) + return RES_BAD_ARG; + if((scnview->mask & S3D_TRACE) == 0) { + log_error(scnview->scn->dev, + "%s: the S3D_TRACE flag is not active onto the submitted scene view.\n", + FUNC_NAME); + return RES_BAD_OP; + } + + *hit = S3D_HIT_NULL; + + /* Initialise the point query */ + query.x = pos[0]; + query.y = pos[1]; + query.z = pos[2]; + query.radius = radius; + query.time = FLT_MAX; /* Invalid fields */ + + /* Initialise the point query context */ + rtcInitPointQueryContext(&query_ctx.rtc); + query_ctx.scnview = scnview; + query_ctx.data = query_data; + + /* Here we go! */ + rtcPointQuery(scnview->rtc_scn, &query, &query_ctx.rtc, closest_point, hit); + + return RES_OK; +} + diff --git a/src/s3d_sphere.h b/src/s3d_sphere.h @@ -100,4 +100,26 @@ sphere_compute_volume(const struct sphere* sphere) return (float)(4*PI*r*r*r/3); } +static FINLINE void +sphere_normal_to_uv(const float normal[3], float uv[2]) +{ + float u, v, cos_theta; + ASSERT(normal && uv && f3_is_normalized(normal)); + + cos_theta = normal[2]; + + v = (1.f - cos_theta) * 0.5f; + if(absf(cos_theta) == 1) { + u = 0; + } else if(eq_epsf(normal[0], 0.f, 1.e-6f)) { + u = normal[1] > 0 ? 0.25f : 0.75f; + } else { + double phi = atan2f(normal[1], normal[0]); /* phi in [-PI, PI] */ + if(phi < 0) phi = 2*PI + phi; /* phi in [0, 2PI] */ + u = (float)(phi / (2*PI)); + } + uv[0] = u; + uv[1] = v; +} + #endif /* S3D_SPHERE_H */