rngrd_mesh.c (11696B)
1 /* Copyright (C) 2022, 2023, 2025 Centre National de la Recherche Scientifique 2 * Copyright (C) 2022, 2023, 2025 Institut Pierre-Simon Laplace 3 * Copyright (C) 2022, 2023, 2025 Institut de Physique du Globe de Paris 4 * Copyright (C) 2022, 2023, 2025, 2026 |Méso|Star> (contact@meso-star.com) 5 * Copyright (C) 2022, 2023, 2025 Observatoire de Paris 6 * Copyright (C) 2022, 2023, 2025 Université de Reims Champagne-Ardenne 7 * Copyright (C) 2022, 2023, 2025 Université de Versaille Saint-Quentin 8 * Copyright (C) 2022, 2023, 2025 Université Paul Sabatier 9 * 10 * This program is free software: you can redistribute it and/or modify 11 * it under the terms of the GNU General Public License as published by 12 * the Free Software Foundation, either version 3 of the License, or 13 * (at your option) any later version. 14 * 15 * This program is distributed in the hope that it will be useful, 16 * but WITHOUT ANY WARRANTY; without even the implied warranty of 17 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 18 * GNU General Public License for more details. 19 * 20 * You should have received a copy of the GNU General Public License 21 * along with this program. If not, see <http://www.gnu.org/licenses/>. */ 22 23 #include "rngrd.h" 24 #include "rngrd_c.h" 25 #include "rngrd_log.h" 26 27 #include <star/s3d.h> 28 #include <star/smsh.h> 29 30 #include <rsys/cstr.h> 31 #include <rsys/float2.h> 32 #include <rsys/float3.h> 33 34 struct filter_args { 35 struct s3d_hit hit; 36 37 s3d_hit_filter_function_T func; /* NULL <=> Stop RT at 1st hit triangle */ 38 void* context; /* User data send to the filter function */ 39 }; 40 static const struct filter_args FILTER_ARGS_NULL = {S3D_HIT_NULL, NULL, NULL}; 41 42 /******************************************************************************* 43 * Helper functions 44 ******************************************************************************/ 45 static INLINE int 46 hit_on_edge 47 (const struct s3d_hit* hit, 48 const float org[3], 49 const float dir[3]) 50 { 51 const float on_edge_epsilon = 1.e-4f; 52 53 struct s3d_attrib v0, v1, v2; 54 float E0[3], E1[3], N[3]; 55 float tri_2area; 56 float hit_2area0; 57 float hit_2area1; 58 float hit_2area2; 59 float hit_pos[3]; 60 ASSERT(hit && !S3D_HIT_NONE(hit) && org && dir); 61 62 /* Retrieve the triangle vertices */ 63 S3D(triangle_get_vertex_attrib(&hit->prim, 0, S3D_POSITION, &v0)); 64 S3D(triangle_get_vertex_attrib(&hit->prim, 1, S3D_POSITION, &v1)); 65 S3D(triangle_get_vertex_attrib(&hit->prim, 2, S3D_POSITION, &v2)); 66 67 /* Compute the intersection position */ 68 f3_add(hit_pos, org, f3_mulf(hit_pos, dir, hit->distance)); 69 70 /* Compute the area of the intersected triangle. Actually we compute the 71 * area*2 to save computation time */ 72 f3_sub(E0, v1.value, v0.value); 73 f3_sub(E1, v2.value, v0.value); 74 tri_2area = f3_len(f3_cross(N, E0, E1)); 75 76 /* Calculate the areas of the 3 triangles formed by an edge of the 77 * intersecting triangle and the position of intersection. Actually we 78 * compute the areas*2 to save computation time. */ 79 f3_sub(E0, v0.value, hit_pos); 80 f3_sub(E1, v1.value, hit_pos); 81 hit_2area0 = f3_len(f3_cross(N, E0, E1)); 82 f3_sub(E0, v1.value, hit_pos); 83 f3_sub(E1, v2.value, hit_pos); 84 hit_2area1 = f3_len(f3_cross(N, E0, E1)); 85 f3_sub(E0, v2.value, hit_pos); 86 f3_sub(E1, v0.value, hit_pos); 87 hit_2area2 = f3_len(f3_cross(N, E0, E1)); 88 89 if(hit_2area0/tri_2area < on_edge_epsilon 90 || hit_2area1/tri_2area < on_edge_epsilon 91 || hit_2area2/tri_2area < on_edge_epsilon) 92 return 1; 93 94 return 0; 95 } 96 97 /* Returns 1 if the intersection found is a self-intersection, i.e. if the 98 * intersection triangle is the triangle from which the ray starts. */ 99 static INLINE int 100 self_hit 101 (const struct s3d_hit* hit, 102 const float ray_org[3], 103 const float ray_dir[3], 104 const float ray_range[2], 105 const struct s3d_hit* hit_from) 106 { 107 ASSERT(hit && hit_from); 108 (void)ray_org, (void)ray_dir; 109 110 if(S3D_HIT_NONE(hit_from)) 111 return 0; 112 113 /* The intersected triangle is the one from which the ray starts. We ignore 114 * this intersection */ 115 if(S3D_PRIMITIVE_EQ(&hit->prim, &hit_from->prim)) 116 return 1; 117 118 /* If the intersection is close to the origin of the ray, we check if it is 119 * on an edge/vertex shared by the triangle from which the ray originates. If 120 * yes, we assume self-intersection and ignore it. */ 121 if(hit->distance/ray_range[1] < 1.e-4f 122 && hit_on_edge(hit_from, ray_org, ray_dir) 123 && hit_on_edge(hit, ray_org, ray_dir)) { 124 return 1; 125 } 126 127 /* No self hit */ 128 return 0; 129 } 130 131 static int 132 mesh_filter 133 (const struct s3d_hit* hit, 134 const float ray_org[3], 135 const float ray_dir[3], 136 const float ray_range[2], 137 void* ray_data, 138 void* filter_data) 139 { 140 const struct filter_args* filter = ray_data; 141 (void)filter_data; 142 143 /* Internally, Star-3D relies on Embree which, due to numerical imprecision, 144 * can find intersections whose distances are not strictly within the ray 145 * range. We reject these intersections. */ 146 if(hit->distance <= ray_range[0] || hit->distance >= ray_range[1]) 147 return 1; 148 149 if(!filter) /* Nothing more to do */ 150 return 0; 151 152 /* Discard this intersection */ 153 if(self_hit(hit, ray_org, ray_dir, ray_range, &filter->hit)) 154 return 1; 155 156 if(!filter->func) /* No user-defined filter functions */ 157 return 0; 158 159 return filter->func /* Invoke user-defined filtering */ 160 (hit, ray_org, ray_dir, ray_range, filter->context, filter_data); 161 } 162 163 static void 164 get_indices(const unsigned itri, unsigned ids[3], void* ctx) 165 { 166 struct smsh_desc* desc = ctx; 167 const uint64_t* indices = NULL; 168 ASSERT(itri < desc->ncells && desc->dcell == 3); 169 indices = smsh_desc_get_cell(desc, itri); 170 ids[0] = (unsigned)indices[0]; 171 ids[1] = (unsigned)indices[1]; 172 ids[2] = (unsigned)indices[2]; 173 } 174 175 static void 176 get_position(const unsigned ivert, float pos[3], void* ctx) 177 { 178 struct smsh_desc* desc = ctx; 179 const double* position = NULL; 180 ASSERT(ivert < desc->nnodes && desc->dnode == 3); 181 position = smsh_desc_get_node(desc, ivert); 182 pos[0] = (float)position[0]; 183 pos[1] = (float)position[1]; 184 pos[2] = (float)position[2]; 185 } 186 187 static res_T 188 check_smsh_desc(const struct rngrd* ground, const struct smsh_desc* desc) 189 { 190 res_T res = RES_OK; 191 ASSERT(ground && desc); 192 193 if(desc->dnode != 3 && desc->dcell != 3) { 194 log_err(ground, 195 "The ground mesh must be a 3D triangular mesh " 196 "(dimension of the mesh: %u; dimension of the vertices: %u)\n", 197 desc->dnode, desc->dcell); 198 res = RES_BAD_ARG; 199 goto error; 200 } 201 202 /* Check number of triangles */ 203 if(desc->ncells > UINT_MAX) { 204 log_err(ground, 205 "The number of triangles cannot exceed %lu whereas it is %lu\n", 206 (unsigned long)UINT_MAX, (unsigned long)desc->ncells); 207 res = RES_BAD_ARG; 208 goto error; 209 } 210 211 /* Check number of vertices */ 212 if(desc->nnodes > UINT_MAX) { 213 log_err(ground, 214 "The number of veritces cannot exceed %lu whereas it is %lu\n", 215 (unsigned long)UINT_MAX, (unsigned long)desc->nnodes); 216 res = RES_BAD_ARG; 217 goto error; 218 } 219 220 exit: 221 return res; 222 error: 223 goto exit; 224 } 225 226 static res_T 227 setup_s3d(struct rngrd* ground, struct smsh_desc* smsh_desc) 228 { 229 struct s3d_vertex_data vdata; 230 struct s3d_shape* mesh = NULL; 231 struct s3d_scene* scene = NULL; 232 res_T res = RES_OK; 233 234 res = s3d_device_create 235 (ground->logger, ground->allocator, 0/*Make Star3D quiet*/, &ground->s3d); 236 if(res != RES_OK) goto error; 237 238 res = s3d_shape_create_mesh(ground->s3d, &mesh); 239 if(res != RES_OK) goto error; 240 res = s3d_mesh_set_hit_filter_function(mesh, mesh_filter, NULL); 241 if(res != RES_OK) goto error; 242 res = s3d_scene_create(ground->s3d, &scene); 243 if(res != RES_OK) goto error; 244 res = s3d_scene_attach_shape(scene, mesh); 245 if(res != RES_OK) goto error; 246 247 vdata.usage = S3D_POSITION; 248 vdata.type = S3D_FLOAT3; 249 vdata.get = get_position; 250 res = s3d_mesh_setup_indexed_vertices(mesh, (unsigned)smsh_desc->ncells, 251 get_indices, (unsigned)smsh_desc->nnodes, &vdata, 1, smsh_desc); 252 if(res != RES_OK) goto error; 253 254 res = s3d_scene_view_create(scene, S3D_TRACE, &ground->s3d_view); 255 if(res != RES_OK) goto error; 256 257 exit: 258 if(mesh) S3D(shape_ref_put(mesh)); 259 if(scene) S3D(scene_ref_put(scene)); 260 return res; 261 error: 262 log_err(ground, "Could not setup the Star-3D data structures -- %s\n", 263 res_to_cstr(res)); 264 if(ground->s3d) S3D(device_ref_put(ground->s3d)); 265 if(ground->s3d_view) S3D(scene_view_ref_put(ground->s3d_view)); 266 ground->s3d = NULL; 267 ground->s3d_view = NULL; 268 goto exit; 269 } 270 271 /******************************************************************************* 272 * Exported functions 273 ******************************************************************************/ 274 res_T 275 rngrd_trace_ray 276 (const struct rngrd* ground, 277 struct rngrd_trace_ray_args* args, 278 struct s3d_hit* hit) 279 { 280 struct filter_args filter = FILTER_ARGS_NULL; 281 float org[3]; 282 float dir[3]; 283 float range[2]; 284 res_T res = RES_OK; 285 286 if(!ground || !args || !hit) { 287 res = RES_BAD_ARG; 288 goto error; 289 } 290 291 f3_set_d3(org, args->ray_org); 292 f3_set_d3(dir, args->ray_dir); 293 f2_set_d2(range, args->ray_range); 294 295 filter.hit = args->hit_from; 296 filter.func = args->filter; 297 filter.context = args->filter_data; 298 299 *hit = S3D_HIT_NULL; 300 301 res = s3d_scene_view_trace_ray(ground->s3d_view, org, dir, range, &filter, hit); 302 if(res != RES_OK) { 303 log_err(ground, 304 "%s: error tracing ray " 305 "(origin = %g, %g, %g; direction = %g, %g ,%g; range = %g, %g)\n", 306 FUNC_NAME, SPLIT3(org), SPLIT3(dir), SPLIT2(range)); 307 goto error; 308 } 309 310 exit: 311 return res; 312 error: 313 goto exit; 314 } 315 316 res_T 317 rngrd_closest_point 318 (const struct rngrd* ground, 319 struct rngrd_closest_point_args* args, 320 struct s3d_hit* hit) 321 { 322 struct filter_args filter = FILTER_ARGS_NULL; 323 float pos[3]; 324 float radius; 325 res_T res = RES_OK; 326 327 if(!ground || !args || !hit) { 328 res = RES_BAD_ARG; 329 goto error; 330 } 331 332 f3_set_d3(pos, args->position); 333 radius = (float)args->radius; 334 filter.func = args->filter; 335 filter.context = args->filter_data; 336 337 *hit = S3D_HIT_NULL; 338 339 res = s3d_scene_view_closest_point(ground->s3d_view, pos, radius, &filter, hit); 340 if(res != RES_OK) { 341 log_err(ground, 342 "%s: error finding the closest point " 343 "(position = %g, %g, %g; radius = %g)\n", 344 FUNC_NAME, SPLIT3(pos), radius); 345 goto error; 346 } 347 348 exit: 349 return res; 350 error: 351 goto exit; 352 } 353 354 /******************************************************************************* 355 * Local function 356 ******************************************************************************/ 357 res_T 358 setup_mesh(struct rngrd* ground, const struct rngrd_create_args* args) 359 { 360 struct smsh_create_args smsh_args = SMSH_CREATE_ARGS_DEFAULT; 361 struct smsh_desc smsh_desc = SMSH_DESC_NULL; 362 struct smsh_load_args smsh_load_args = SMSH_LOAD_ARGS_NULL; 363 struct smsh* smsh = NULL; 364 res_T res = RES_OK; 365 ASSERT(ground && args); 366 367 /* Create the Star-Mesh loader */ 368 smsh_args.logger = ground->logger; 369 smsh_args.allocator = ground->allocator; 370 smsh_args.verbose = ground->verbose; 371 res = smsh_create(&smsh_args, &smsh); 372 if(res != RES_OK) goto error; 373 374 /* Load and retrieve the Star-Mesh data */ 375 smsh_load_args.path = args->smsh_filename; 376 smsh_load_args.memory_mapping = 1; 377 res = smsh_load(smsh, &smsh_load_args); 378 if(res != RES_OK) goto error; 379 res = smsh_get_desc(smsh, &smsh_desc); 380 if(res != RES_OK) goto error; 381 res = check_smsh_desc(ground, &smsh_desc); 382 if(res != RES_OK) goto error; 383 384 /* Setup the Star-3D data structures */ 385 res = setup_s3d(ground, &smsh_desc); 386 if(res != RES_OK) goto error; 387 388 ground->ntriangles = smsh_desc.ncells; 389 390 exit: 391 if(smsh) SMSH(ref_put(smsh)); 392 return res; 393 error: 394 if(ground->s3d) { 395 S3D(device_ref_put(ground->s3d)); 396 ground->s3d = NULL; 397 } 398 goto exit; 399 }