rngrd_mesh.c (10523B)
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 |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 /******************************************************************************* 35 * Helper functions 36 ******************************************************************************/ 37 static INLINE int 38 hit_on_edge 39 (const struct s3d_hit* hit, 40 const float org[3], 41 const float dir[3]) 42 { 43 const float on_edge_epsilon = 1.e-4f; 44 45 struct s3d_attrib v0, v1, v2; 46 float E0[3], E1[3], N[3]; 47 float tri_2area; 48 float hit_2area0; 49 float hit_2area1; 50 float hit_2area2; 51 float hit_pos[3]; 52 ASSERT(hit && !S3D_HIT_NONE(hit) && org && dir); 53 54 /* Retrieve the triangle vertices */ 55 S3D(triangle_get_vertex_attrib(&hit->prim, 0, S3D_POSITION, &v0)); 56 S3D(triangle_get_vertex_attrib(&hit->prim, 1, S3D_POSITION, &v1)); 57 S3D(triangle_get_vertex_attrib(&hit->prim, 2, S3D_POSITION, &v2)); 58 59 /* Compute the intersection position */ 60 f3_add(hit_pos, org, f3_mulf(hit_pos, dir, hit->distance)); 61 62 /* Compute the area of the intersected triangle. Actually we compute the 63 * area*2 to save computation time */ 64 f3_sub(E0, v1.value, v0.value); 65 f3_sub(E1, v2.value, v0.value); 66 tri_2area = f3_len(f3_cross(N, E0, E1)); 67 68 /* Calculate the areas of the 3 triangles formed by an edge of the 69 * intersecting triangle and the position of intersection. Actually we 70 * compute the areas*2 to save computation time. */ 71 f3_sub(E0, v0.value, hit_pos); 72 f3_sub(E1, v1.value, hit_pos); 73 hit_2area0 = f3_len(f3_cross(N, E0, E1)); 74 f3_sub(E0, v1.value, hit_pos); 75 f3_sub(E1, v2.value, hit_pos); 76 hit_2area1 = f3_len(f3_cross(N, E0, E1)); 77 f3_sub(E0, v2.value, hit_pos); 78 f3_sub(E1, v0.value, hit_pos); 79 hit_2area2 = f3_len(f3_cross(N, E0, E1)); 80 81 if(hit_2area0/tri_2area < on_edge_epsilon 82 || hit_2area1/tri_2area < on_edge_epsilon 83 || hit_2area2/tri_2area < on_edge_epsilon) 84 return 1; 85 86 return 0; 87 } 88 89 /* Returns 1 if the intersection found is a self-intersection, i.e. if the 90 * intersection triangle is the triangle from which the ray starts. */ 91 static INLINE int 92 self_hit 93 (const struct s3d_hit* hit, 94 const float ray_org[3], 95 const float ray_dir[3], 96 const float ray_range[2], 97 const struct s3d_hit* hit_from) 98 { 99 ASSERT(hit && hit_from); 100 (void)ray_org, (void)ray_dir; 101 102 if(S3D_HIT_NONE(hit_from)) 103 return 0; 104 105 /* The intersected triangle is the one from which the ray starts. We ignore 106 * this intersection */ 107 if(S3D_PRIMITIVE_EQ(&hit->prim, &hit_from->prim)) 108 return 1; 109 110 /* If the intersection is close to the origin of the ray, we check if it is 111 * on an edge/vertex shared by the triangle from which the ray originates. If 112 * yes, we assume self-intersection and ignore it. */ 113 if(hit->distance/ray_range[1] < 1.e-4f 114 && hit_on_edge(hit_from, ray_org, ray_dir) 115 && hit_on_edge(hit, ray_org, ray_dir)) { 116 return 1; 117 } 118 119 /* No self hit */ 120 return 0; 121 } 122 123 static int 124 mesh_filter 125 (const struct s3d_hit* hit, 126 const float ray_org[3], 127 const float ray_dir[3], 128 const float ray_range[2], 129 void* ray_data, 130 void* filter_data) 131 { 132 const struct rngrd_trace_ray_args* ray_args = ray_data; 133 (void)filter_data; 134 135 /* Internally, Star-3D relies on Embree which, due to numerical imprecision, 136 * can find intersections whose distances are not strictly within the ray 137 * range. We reject these intersections. */ 138 if(hit->distance <= ray_range[0] || hit->distance >= ray_range[1]) 139 return 1; 140 141 if(!ray_args) /* Nothing more to do */ 142 return 0; 143 144 /* Discard this intersection */ 145 if(self_hit(hit, ray_org, ray_dir, ray_range, &ray_args->hit_from)) 146 return 1; 147 148 if(!ray_args->filter) /* No user-defined filter functions */ 149 return 0; 150 151 return ray_args->filter /* Invoke user-defined filtering */ 152 (hit, ray_org, ray_dir, ray_range, ray_args->filter_data, filter_data); 153 } 154 155 static void 156 get_indices(const unsigned itri, unsigned ids[3], void* ctx) 157 { 158 struct smsh_desc* desc = ctx; 159 const uint64_t* indices = NULL; 160 ASSERT(itri < desc->ncells && desc->dcell == 3); 161 indices = smsh_desc_get_cell(desc, itri); 162 ids[0] = (unsigned)indices[0]; 163 ids[1] = (unsigned)indices[1]; 164 ids[2] = (unsigned)indices[2]; 165 } 166 167 static void 168 get_position(const unsigned ivert, float pos[3], void* ctx) 169 { 170 struct smsh_desc* desc = ctx; 171 const double* position = NULL; 172 ASSERT(ivert < desc->nnodes && desc->dnode == 3); 173 position = smsh_desc_get_node(desc, ivert); 174 pos[0] = (float)position[0]; 175 pos[1] = (float)position[1]; 176 pos[2] = (float)position[2]; 177 } 178 179 static res_T 180 check_smsh_desc(const struct rngrd* ground, const struct smsh_desc* desc) 181 { 182 res_T res = RES_OK; 183 ASSERT(ground && desc); 184 185 if(desc->dnode != 3 && desc->dcell != 3) { 186 log_err(ground, 187 "The ground mesh must be a 3D triangular mesh " 188 "(dimension of the mesh: %u; dimension of the vertices: %u)\n", 189 desc->dnode, desc->dcell); 190 res = RES_BAD_ARG; 191 goto error; 192 } 193 194 /* Check number of triangles */ 195 if(desc->ncells > UINT_MAX) { 196 log_err(ground, 197 "The number of triangles cannot exceed %lu whereas it is %lu\n", 198 (unsigned long)UINT_MAX, (unsigned long)desc->ncells); 199 res = RES_BAD_ARG; 200 goto error; 201 } 202 203 /* Check number of vertices */ 204 if(desc->nnodes > UINT_MAX) { 205 log_err(ground, 206 "The number of veritces cannot exceed %lu whereas it is %lu\n", 207 (unsigned long)UINT_MAX, (unsigned long)desc->nnodes); 208 res = RES_BAD_ARG; 209 goto error; 210 } 211 212 exit: 213 return res; 214 error: 215 goto exit; 216 } 217 218 static res_T 219 setup_s3d(struct rngrd* ground, struct smsh_desc* smsh_desc) 220 { 221 struct s3d_vertex_data vdata; 222 struct s3d_shape* mesh = NULL; 223 struct s3d_scene* scene = NULL; 224 res_T res = RES_OK; 225 226 res = s3d_device_create 227 (ground->logger, ground->allocator, 0/*Make Star3D quiet*/, &ground->s3d); 228 if(res != RES_OK) goto error; 229 230 res = s3d_shape_create_mesh(ground->s3d, &mesh); 231 if(res != RES_OK) goto error; 232 res = s3d_mesh_set_hit_filter_function(mesh, mesh_filter, NULL); 233 if(res != RES_OK) goto error; 234 res = s3d_scene_create(ground->s3d, &scene); 235 if(res != RES_OK) goto error; 236 res = s3d_scene_attach_shape(scene, mesh); 237 if(res != RES_OK) goto error; 238 239 vdata.usage = S3D_POSITION; 240 vdata.type = S3D_FLOAT3; 241 vdata.get = get_position; 242 res = s3d_mesh_setup_indexed_vertices(mesh, (unsigned)smsh_desc->ncells, 243 get_indices, (unsigned)smsh_desc->nnodes, &vdata, 1, smsh_desc); 244 if(res != RES_OK) goto error; 245 246 res = s3d_scene_view_create(scene, S3D_TRACE, &ground->s3d_view); 247 if(res != RES_OK) goto error; 248 249 exit: 250 if(mesh) S3D(shape_ref_put(mesh)); 251 if(scene) S3D(scene_ref_put(scene)); 252 return res; 253 error: 254 log_err(ground, "Could not setup the Star-3D data structures -- %s\n", 255 res_to_cstr(res)); 256 if(ground->s3d) S3D(device_ref_put(ground->s3d)); 257 if(ground->s3d_view) S3D(scene_view_ref_put(ground->s3d_view)); 258 ground->s3d = NULL; 259 ground->s3d_view = NULL; 260 goto exit; 261 } 262 263 /******************************************************************************* 264 * Exported functions 265 ******************************************************************************/ 266 res_T 267 rngrd_trace_ray 268 (const struct rngrd* ground, 269 struct rngrd_trace_ray_args* args, 270 struct s3d_hit* hit) 271 { 272 float org[3]; 273 float dir[3]; 274 float range[2]; 275 res_T res = RES_OK; 276 277 if(!ground || !args || !hit) { 278 res = RES_BAD_ARG; 279 goto error; 280 } 281 282 f3_set_d3(org, args->ray_org); 283 f3_set_d3(dir, args->ray_dir); 284 f2_set_d2(range, args->ray_range); 285 286 *hit = S3D_HIT_NULL; 287 288 res = s3d_scene_view_trace_ray(ground->s3d_view, org, dir, range, args, hit); 289 if(res != RES_OK) { 290 log_err(ground, 291 "%s: error tracing ray " 292 "(origin = %g, %g, %g; direction = %g, %g ,%g; range = %g, %g)\n", 293 FUNC_NAME, SPLIT3(org), SPLIT3(dir), SPLIT2(range)); 294 goto error; 295 } 296 297 exit: 298 return res; 299 error: 300 goto exit; 301 } 302 303 /******************************************************************************* 304 * Local function 305 ******************************************************************************/ 306 res_T 307 setup_mesh(struct rngrd* ground, const struct rngrd_create_args* args) 308 { 309 struct smsh_create_args smsh_args = SMSH_CREATE_ARGS_DEFAULT; 310 struct smsh_desc smsh_desc = SMSH_DESC_NULL; 311 struct smsh_load_args smsh_load_args = SMSH_LOAD_ARGS_NULL; 312 struct smsh* smsh = NULL; 313 res_T res = RES_OK; 314 ASSERT(ground && args); 315 316 /* Create the Star-Mesh loader */ 317 smsh_args.logger = ground->logger; 318 smsh_args.allocator = ground->allocator; 319 smsh_args.verbose = ground->verbose; 320 res = smsh_create(&smsh_args, &smsh); 321 if(res != RES_OK) goto error; 322 323 /* Load and retrieve the Star-Mesh data */ 324 smsh_load_args.path = args->smsh_filename; 325 smsh_load_args.memory_mapping = 1; 326 res = smsh_load(smsh, &smsh_load_args); 327 if(res != RES_OK) goto error; 328 res = smsh_get_desc(smsh, &smsh_desc); 329 if(res != RES_OK) goto error; 330 res = check_smsh_desc(ground, &smsh_desc); 331 if(res != RES_OK) goto error; 332 333 /* Setup the Star-3D data structures */ 334 res = setup_s3d(ground, &smsh_desc); 335 if(res != RES_OK) goto error; 336 337 ground->ntriangles = smsh_desc.ncells; 338 339 exit: 340 if(smsh) SMSH(ref_put(smsh)); 341 return res; 342 error: 343 if(ground->s3d) { 344 S3D(device_ref_put(ground->s3d)); 345 ground->s3d = NULL; 346 } 347 goto exit; 348 }