sdis_scene_Xd.h (44696B)
1 /* Copyright (C) 2016-2025 |Méso|Star> (contact@meso-star.com) 2 * 3 * This program is free software: you can redistribute it and/or modify 4 * it under the terms of the GNU General Public License as published by 5 * the Free Software Foundation, either version 3 of the License, or 6 * (at your option) any later version. 7 * 8 * This program is distributed in the hope that it will be useful, 9 * but WITHOUT ANY WARRANTY; without even the implied warranty of 10 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 11 * GNU General Public License for more details. 12 * 13 * You should have received a copy of the GNU General Public License 14 * along with this program. If not, see <http://www.gnu.org/licenses/>. */ 15 16 #ifndef SDIS_SCENE_XD_H 17 #define SDIS_SCENE_XD_H 18 19 #include "sdis_interface_c.h" 20 #include "sdis_log.h" 21 #include "sdis_medium_c.h" 22 #include "sdis_scene_c.h" 23 24 #include <star/ssp.h> 25 #include <star/senc2d.h> 26 #include <star/senc3d.h> 27 28 #include <rsys/cstr.h> 29 #include <rsys/float22.h> 30 #include <rsys/float33.h> 31 #include <rsys/rsys.h> 32 33 /* Emperical cos threshold defining if an angle is sharp */ 34 #define SHARP_ANGLE_COS_THRESOLD -0.70710678 /* ~ cos(3*PI/4) */ 35 36 /******************************************************************************* 37 * Define the helper functions and the data types used by the scene 38 * independently of its dimension, i.e. 2D or 3D. 39 ******************************************************************************/ 40 /* Context used to wrap the user geometry and interfaces to Star-Enc */ 41 struct geometry { 42 void (*indices)(const size_t iprim, size_t ids[], void*); 43 void (*interf)(const size_t iprim, struct sdis_interface**, void*); 44 void (*position)(const size_t ivert, double pos[], void*); 45 void* data; 46 }; 47 48 /* Fetch the media split by the primitive `iprim'. This first and second media 49 * are the media from the front face side and back face side of the primitive, 50 * respectively. */ 51 static void 52 geometry_media(const unsigned iprim, unsigned media[2], void* data) 53 { 54 struct geometry* ctx = data; 55 struct sdis_interface* interf; 56 ASSERT(ctx && media); 57 ctx->interf(iprim, &interf, ctx->data); 58 media[0] = medium_get_id(interf->medium_front); 59 media[1] = medium_get_id(interf->medium_back); 60 } 61 62 /* Register the submitted medium against the scene if it is not already 63 * registered. On registration, no reference is taken onto the medium; the 64 * scene references its media through its interfaces and it is thus useless to 65 * take another reference onto them. */ 66 static res_T 67 register_medium(struct sdis_scene* scn, struct sdis_medium* mdm) 68 { 69 unsigned id; 70 size_t nmedia; 71 res_T res = RES_OK; 72 ASSERT(scn && mdm); 73 74 /* Check that the medium is already registered against the scene */ 75 id = medium_get_id(mdm); 76 nmedia = darray_medium_size_get(&scn->media); 77 if(id >= nmedia) { 78 res = darray_medium_resize(&scn->media, id + 1); 79 if(res != RES_OK) return res; 80 } 81 if(darray_medium_cdata_get(&scn->media)[id]) { 82 ASSERT(darray_medium_cdata_get(&scn->media)[id] == mdm); 83 } else { 84 /* Do not take a reference onto the medium since we already take a 85 * reference onto at least one interface that uses it, and thus that has a 86 * reference onto it */ 87 darray_medium_data_get(&scn->media)[id] = mdm; 88 } 89 return RES_OK; 90 } 91 92 /* Release the reference toward the interfaces and thus clear the list of scene 93 * interfaces, the list of scene media, and the list of per-primitive 94 * interface. */ 95 static void 96 clear_properties(struct sdis_scene* scn) 97 { 98 size_t i; 99 ASSERT(scn); 100 FOR_EACH(i, 0, darray_interf_size_get(&scn->interfaces)) { 101 if(darray_interf_cdata_get(&scn->interfaces)[i]) { 102 SDIS(interface_ref_put(darray_interf_data_get(&scn->interfaces)[i])); 103 } 104 } 105 darray_interf_clear(&scn->interfaces); 106 darray_medium_clear(&scn->media); 107 darray_prim_prop_clear(&scn->prim_props); 108 } 109 110 static INLINE int 111 check_sdis_scene_create_args(const struct sdis_scene_create_args* args) 112 { 113 return args 114 && args->get_indices 115 && args->get_interface 116 && args->get_position 117 && args->nprimitives 118 && args->nprimitives < UINT_MAX 119 && args->nvertices 120 && args->nvertices < UINT_MAX 121 && args->fp_to_meter > 0; 122 } 123 124 static INLINE res_T 125 check_sdis_scene_find_closest_point_args 126 (const struct sdis_scene_find_closest_point_args* args) 127 { 128 /* Undefined input arguments */ 129 if(!args) return RES_BAD_ARG; 130 131 /* Invalid radius */ 132 if(args->radius <= 0) return RES_BAD_ARG; 133 134 return RES_OK; 135 } 136 137 #endif /* SDIS_SCENE_XD_H */ 138 139 #include "sdis_device_c.h" 140 141 #include <rsys/float2.h> 142 #include <rsys/float3.h> 143 #include <rsys/double2.h> 144 #include <rsys/double3.h> 145 #include <rsys/mem_allocator.h> 146 147 #include <limits.h> 148 149 /* Check the submitted dimension and include its specific headers */ 150 #define SENCXD_DIM SDIS_XD_DIMENSION 151 #if (SDIS_XD_DIMENSION == 2) 152 #include <star/sencX2d.h> 153 #include <star/s2d.h> 154 #elif (SDIS_XD_DIMENSION == 3) 155 #include <star/sencX3d.h> 156 #include <star/s3d.h> 157 #else 158 #error "Invalid SDIS_XD_DIMENSION value." 159 #endif 160 161 #include "sdis_Xd_begin.h" 162 163 #if DIM == 2 164 #define HIT_ON_BOUNDARY hit_on_vertex 165 #else 166 #define HIT_ON_BOUNDARY hit_on_edge 167 #endif 168 169 /******************************************************************************* 170 * Helper functions 171 ******************************************************************************/ 172 #if DIM == 2 173 #define ON_VERTEX_EPSILON 1.e-4f 174 /* Check that `hit' roughly lies on a vertex. */ 175 static INLINE int 176 hit_on_vertex 177 (const struct s2d_hit* hit, 178 const float org[2], 179 const float dir[2]) 180 { 181 struct s2d_attrib v0, v1; 182 float E[2]; 183 float hit_pos[2]; 184 float segment_len; 185 float hit_len0; 186 float hit_len1; 187 ASSERT(hit && !S2D_HIT_NONE(hit) && org && dir); 188 189 /* Rertieve the segment vertices */ 190 S2D(segment_get_vertex_attrib(&hit->prim, 0, S2D_POSITION, &v0)); 191 S2D(segment_get_vertex_attrib(&hit->prim, 1, S2D_POSITION, &v1)); 192 193 /* Compute the length of the segment */ 194 segment_len = f2_len(f2_sub(E, v1.value, v0.value)); 195 196 /* Compute the hit position onto the segment */ 197 f2_add(hit_pos, org, f2_mulf(hit_pos, dir, hit->distance)); 198 199 /* Compute the length from hit position to segment vertices */ 200 hit_len0 = f2_len(f2_sub(E, v0.value, hit_pos)); 201 hit_len1 = f2_len(f2_sub(E, v1.value, hit_pos)); 202 203 if(hit_len0 / segment_len < ON_VERTEX_EPSILON 204 || hit_len1 / segment_len < ON_VERTEX_EPSILON) 205 return 1; 206 return 0; 207 } 208 209 static int 210 hit_shared_vertex 211 (const struct s2d_primitive* seg0, 212 const struct s2d_primitive* seg1, 213 const float pos0[2], /* Tested position onto the segment 0 */ 214 const float pos1[2]) /* Tested Position onto the segment 1 */ 215 { 216 struct s2d_attrib seg0_vertices[2]; /* Vertex positions of the segment 0 */ 217 struct s2d_attrib seg1_vertices[2]; /* Vertex positions of the segment 1 */ 218 float d0[2], d1[2]; /* temporary vector */ 219 float seg0_len, seg1_len; /* Length of the segments */ 220 float tmp0_len, tmp1_len; 221 float cos_normals; 222 int seg0_vert = -1; /* Id of the shared vertex for the segment 0 */ 223 int seg1_vert = -1; /* Id of the shared vertex for the segment 1 */ 224 int seg0_ivertex, seg1_ivertex; 225 ASSERT(seg0 && seg1 && pos0 && pos1); 226 227 /* Fetch the vertices of the segment 0 */ 228 S2D(segment_get_vertex_attrib(seg0, 0, S2D_POSITION, &seg0_vertices[0])); 229 S2D(segment_get_vertex_attrib(seg0, 1, S2D_POSITION, &seg0_vertices[1])); 230 231 /* Fetch the vertices of the segment 1 */ 232 S2D(segment_get_vertex_attrib(seg1, 0, S2D_POSITION, &seg1_vertices[0])); 233 S2D(segment_get_vertex_attrib(seg1, 1, S2D_POSITION, &seg1_vertices[1])); 234 235 /* Look for the vertex shared by the 2 segments */ 236 for(seg0_ivertex = 0; seg0_ivertex < 2 && seg0_vert < 0; ++seg0_ivertex) { 237 for(seg1_ivertex = 0; seg1_ivertex < 2 && seg1_vert < 0; ++seg1_ivertex) { 238 const int vertex_eq = f2_eq_eps 239 (seg0_vertices[seg0_ivertex].value, 240 seg1_vertices[seg1_ivertex].value, 241 1.e-6f); 242 if(vertex_eq) { 243 seg0_vert = seg0_ivertex; 244 seg1_vert = seg1_ivertex; 245 /* We assume that the segments are not degenerated. As a consequence we 246 * can break here since a vertex of the segment 0 can be equal to at most 247 * one vertex of the segment 1 */ 248 break; 249 } 250 }} 251 252 /* The segments do not have a common vertex */ 253 if(seg0_vert < 0) return 0; 254 255 /* Compute the dirctions from shared vertex to the opposite segment vertex */ 256 f2_sub(d0, seg0_vertices[(seg0_vert+1)%2].value, seg0_vertices[seg0_vert].value); 257 f2_sub(d1, seg1_vertices[(seg1_vert+1)%2].value, seg1_vertices[seg1_vert].value); 258 259 /* Compute the cosine between the segments */ 260 seg0_len = f2_normalize(d0, d0); 261 seg1_len = f2_normalize(d1, d1); 262 cos_normals = f2_dot(d0, d1); 263 264 /* The angle formed by the 2 segments is sharp. Do not filter the hit */ 265 if(cos_normals > SHARP_ANGLE_COS_THRESOLD) return 0; 266 267 /* Compute the length from pos<0|1> to shared vertex */ 268 f2_sub(d0, seg0_vertices[seg0_vert].value, pos0); 269 f2_sub(d1, seg1_vertices[seg1_vert].value, pos1); 270 tmp0_len = f2_len(d0); 271 tmp1_len = f2_len(d1); 272 273 return (eq_epsf(seg0_len, 0, 1.e-6f) || tmp0_len/seg0_len < ON_VERTEX_EPSILON) 274 && (eq_epsf(seg1_len, 0, 1.e-6f) || tmp1_len/seg1_len < ON_VERTEX_EPSILON); 275 } 276 277 #else /* DIM == 3 */ 278 #define ON_EDGE_EPSILON 1.e-4f 279 /* Check that `hit' roughly lies on an edge. */ 280 static INLINE int 281 hit_on_edge 282 (const struct s3d_hit* hit, 283 const float org[3], 284 const float dir[3]) 285 { 286 struct s3d_attrib v0, v1, v2; 287 float E0[3], E1[3], N[3]; 288 float tri_2area; 289 float hit_2area0; 290 float hit_2area1; 291 float hit_2area2; 292 float hit_pos[3]; 293 ASSERT(hit && !S3D_HIT_NONE(hit) && org && dir); 294 295 /* Retrieve the triangle vertices */ 296 S3D(triangle_get_vertex_attrib(&hit->prim, 0, S3D_POSITION, &v0)); 297 S3D(triangle_get_vertex_attrib(&hit->prim, 1, S3D_POSITION, &v1)); 298 S3D(triangle_get_vertex_attrib(&hit->prim, 2, S3D_POSITION, &v2)); 299 300 /* Compute the triangle area * 2 */ 301 f3_sub(E0, v1.value, v0.value); 302 f3_sub(E1, v2.value, v0.value); 303 tri_2area = f3_len(f3_cross(N, E0, E1)); 304 305 /* Compute the hit position */ 306 f3_add(hit_pos, org, f3_mulf(hit_pos, dir, hit->distance)); 307 308 /* Compute areas */ 309 f3_sub(E0, v0.value, hit_pos); 310 f3_sub(E1, v1.value, hit_pos); 311 hit_2area0 = f3_len(f3_cross(N, E0, E1)); 312 f3_sub(E0, v1.value, hit_pos); 313 f3_sub(E1, v2.value, hit_pos); 314 hit_2area1 = f3_len(f3_cross(N, E0, E1)); 315 f3_sub(E0, v2.value, hit_pos); 316 f3_sub(E1, v0.value, hit_pos); 317 hit_2area2 = f3_len(f3_cross(N, E0, E1)); 318 319 if(hit_2area0 / tri_2area < ON_EDGE_EPSILON 320 || hit_2area1 / tri_2area < ON_EDGE_EPSILON 321 || hit_2area2 / tri_2area < ON_EDGE_EPSILON) 322 return 1; 323 324 return 0; 325 } 326 327 static int 328 hit_shared_edge 329 (const struct s3d_primitive* tri0, 330 const struct s3d_primitive* tri1, 331 const float uv0[2], /* Barycentric coordinates of tested position on tri0 */ 332 const float uv1[2], /* Barycentric coordinates of tested position on tri1 */ 333 const float pos0[3], /* Tested position onto the triangle 0 */ 334 const float pos1[3]) /* Tested Position onto the triangle 1 */ 335 { 336 struct s3d_attrib tri0_vertices[3]; /* Vertex positions of the triangle 0 */ 337 struct s3d_attrib tri1_vertices[3]; /* Vertex positions of the triangle 1 */ 338 int tri0_edge[2] = {-1, -1}; /* Shared edge vertex ids for the triangle 0 */ 339 int tri1_edge[2] = {-1, -1}; /* Shared edge vertex ids for the triangle 1 */ 340 int edge_ivertex = 0; /* Temporary variable */ 341 int tri0_ivertex, tri1_ivertex; 342 ASSERT(tri0 && tri1 && pos0 && pos1); 343 344 /* Fetch the vertices of the triangle 0 */ 345 S3D(triangle_get_vertex_attrib(tri0, 0, S3D_POSITION, &tri0_vertices[0])); 346 S3D(triangle_get_vertex_attrib(tri0, 1, S3D_POSITION, &tri0_vertices[1])); 347 S3D(triangle_get_vertex_attrib(tri0, 2, S3D_POSITION, &tri0_vertices[2])); 348 349 /* Fetch the vertices of the triangle 1 */ 350 S3D(triangle_get_vertex_attrib(tri1, 0, S3D_POSITION, &tri1_vertices[0])); 351 S3D(triangle_get_vertex_attrib(tri1, 1, S3D_POSITION, &tri1_vertices[1])); 352 S3D(triangle_get_vertex_attrib(tri1, 2, S3D_POSITION, &tri1_vertices[2])); 353 354 /* Look for the vertices shared by the 2 triangles */ 355 for(tri0_ivertex=0; tri0_ivertex < 3 && edge_ivertex < 2; ++tri0_ivertex) { 356 for(tri1_ivertex=0; tri1_ivertex < 3 && edge_ivertex < 2; ++tri1_ivertex) { 357 const int vertex_eq = f3_eq_eps 358 (tri0_vertices[tri0_ivertex].value, 359 tri1_vertices[tri1_ivertex].value, 360 1.e-6f); 361 if(vertex_eq) { 362 tri0_edge[edge_ivertex] = tri0_ivertex; 363 tri1_edge[edge_ivertex] = tri1_ivertex; 364 ++edge_ivertex; 365 /* We assume that the triangles are not degenerated. As a consequence we 366 * can break here since a vertex of the triangle 0 can be equal to at 367 * most one vertex of the triangle 1 */ 368 break; 369 } 370 }} 371 372 /* The triangles do not have a common edge */ 373 if(edge_ivertex == 0) { 374 return 0; 375 376 /* The triangles have a common vertex */ 377 } else if(edge_ivertex == 1) { 378 float bcoord0, bcoord1; 379 int hit_vertex; 380 381 /* Retrieve the barycentric coordinate of the position on triangle 0 382 * corresponding to the vertex shared between the 2 triangles. */ 383 switch(tri0_edge[0]) { 384 case 0: bcoord0 = uv0[0]; break; 385 case 1: bcoord0 = uv0[1]; break; 386 case 2: bcoord0 = CLAMP(1.f - uv0[0] - uv0[1], 0.f, 1.f); break; 387 default: FATAL("Unreachable code\n"); break; 388 } 389 390 /* Retrieve the barycentric coordinate of the position on triangle 1 391 * corresponding to the vertex shared between the 2 triangles. */ 392 switch(tri1_edge[0]) { 393 case 0: bcoord1 = uv1[0]; break; 394 case 1: bcoord1 = uv1[1]; break; 395 case 2: bcoord1 = CLAMP(1.f - uv0[0] - uv0[1], 0.f, 1.f); break; 396 default: FATAL("Unreachable code\n"); break; 397 } 398 399 /* Check that the both positions lie on the shared vertex */ 400 hit_vertex = eq_epsf(1.f, bcoord0, ON_VERTEX_EPSILON) 401 && eq_epsf(1.f, bcoord1, ON_VERTEX_EPSILON); 402 return hit_vertex; 403 404 /* The triangles have a common edge */ 405 } else { 406 float E0[3], E1[3]; /* Temporary variables storing triangle edges */ 407 float N0[3], N1[3]; /* Temporary Normals */ 408 float tri0_2area, tri1_2area; /* 2*area of the submitted triangles */ 409 float tmp0_2area, tmp1_2area; 410 float cos_normals; 411 int iv0, iv1, iv2; 412 int hit_edge; 413 414 /* Ensure that the vertices of the shared edge are registered in the right 415 * order regarding the triangle vertices, i.e. (0,1), (1,2) or (2,0) */ 416 if((tri0_edge[0]+1)%3 != tri0_edge[1]) SWAP(int, tri0_edge[0], tri0_edge[1]); 417 if((tri1_edge[0]+1)%3 != tri1_edge[1]) SWAP(int, tri1_edge[0], tri1_edge[1]); 418 419 /* Compute the shared edge normal lying in the triangle 0 plane */ 420 iv0 = tri0_edge[0]; 421 iv1 = tri0_edge[1]; 422 iv2 = (tri0_edge[1]+1) % 3; 423 f3_sub(E0, tri0_vertices[iv1].value, tri0_vertices[iv0].value); 424 f3_sub(E1, tri0_vertices[iv2].value, tri0_vertices[iv0].value); 425 f3_cross(N0, E0, E1); /* Triangle 0 normal */ 426 tri0_2area = f3_len(N0); 427 f3_cross(N0, N0, E0); 428 429 /* Compute the shared edge normal lying in the triangle 1 plane */ 430 iv0 = tri1_edge[0]; 431 iv1 = tri1_edge[1]; 432 iv2 = (tri1_edge[1]+1) % 3; 433 f3_sub(E0, tri1_vertices[iv1].value, tri1_vertices[iv0].value); 434 f3_sub(E1, tri1_vertices[iv2].value, tri1_vertices[iv0].value); 435 f3_cross(N1, E0, E1); 436 tri1_2area = f3_len(N1); 437 f3_cross(N1, N1, E0); 438 439 /* Compute the cosine between the 2 edge normals */ 440 f3_normalize(N0, N0); 441 f3_normalize(N1, N1); 442 cos_normals = f3_dot(N0, N1); 443 444 /* The angle formed by the 2 triangles is sharp */ 445 if(cos_normals > SHARP_ANGLE_COS_THRESOLD) return 0; 446 447 /* Compute the 2 times the area of the (pos0, shared_edge.vertex0, 448 * shared_edge.vertex1) triangles */ 449 f3_sub(E0, tri0_vertices[tri0_edge[0]].value, pos0); 450 f3_sub(E1, tri0_vertices[tri0_edge[1]].value, pos0); 451 tmp0_2area = f3_len(f3_cross(N0, E0, E1)); 452 453 /* Compute the 2 times the area of the (pos1, shared_edge.vertex0, 454 * shared_edge.vertex1) triangles */ 455 f3_sub(E0, tri1_vertices[tri1_edge[0]].value, pos1); 456 f3_sub(E1, tri1_vertices[tri1_edge[1]].value, pos1); 457 tmp1_2area = f3_len(f3_cross(N1, E0, E1)); 458 459 hit_edge = 460 ( eq_epsf(tri0_2area, 0, 1.e-6f) 461 || eq_epsf(tmp0_2area, 0, 1.e-6f) 462 || tmp0_2area/tri0_2area < ON_EDGE_EPSILON); 463 hit_edge = hit_edge && 464 ( eq_epsf(tri1_2area, 0, 1.e-6f) 465 || eq_epsf(tmp1_2area, 0, 1.e-6f) 466 || tmp1_2area/tri1_2area < ON_EDGE_EPSILON); 467 return hit_edge; 468 } 469 } 470 #undef ON_EDGE_EPSILON 471 #endif /* DIM == 2 */ 472 473 /* Avoid self-intersection for a ray starting from a planar primitive, i.e. a 474 * triangle or a line segment */ 475 static int 476 XD(hit_filter_function) 477 (const struct sXd(hit)* hit, 478 const float org[DIM], 479 const float dir[DIM], 480 const float range[2], 481 void* query_data, 482 void* global_data) 483 { 484 const struct hit_filter_data* filter_data = query_data; 485 const struct sXd(hit)* hit_from = NULL; 486 (void)org, (void)dir, (void)global_data, (void)range; 487 488 /* No user defined data. Do not filter */ 489 if(!filter_data) return 0; 490 491 /* Call the custom filter function if it exists 492 * or perform regular filtering otherwise */ 493 if(filter_data->XD(custom_filter)) { 494 return filter_data->XD(custom_filter) 495 (hit, org, dir, range, filter_data->custom_filter_data, global_data); 496 } 497 498 /* There is no intersection to discard */ 499 hit_from = &filter_data->XD(hit); 500 if(SXD_HIT_NONE(hit_from)) return 0; 501 502 if(SXD_PRIMITIVE_EQ(&hit_from->prim, &hit->prim)) return 1; 503 504 /* No displacement => assume self intersection in all situations */ 505 if(hit->distance <= 0) return 1; 506 507 if(eq_epsf(hit->distance, 0, (float)filter_data->epsilon)) { 508 float pos[DIM]; 509 int reject_hit = 0; 510 fX(add)(pos, org, fX(mulf)(pos, dir, hit->distance)); 511 /* If the targeted point is near of the origin, check that it lies on an 512 * edge/vertex shared by the 2 primitives */ 513 #if DIM == 2 514 reject_hit = hit_shared_vertex(&hit_from->prim, &hit->prim, org, pos); 515 #else 516 reject_hit = hit_shared_edge 517 (&hit_from->prim, &hit->prim, hit_from->uv, hit->uv, org, pos); 518 #endif 519 if(reject_hit) return 1; 520 } 521 522 /* If the hit to be considered is (approximately) on a boundary between 2 523 * primitives, it may belong to the wrong primitive, i.e. the one that doesn't 524 * "face" the direction of the ray. We therefore check the enclosure towards 525 * which it is directed, and reject it if it is not the same as the one from 526 * which the ray originates */ 527 if(filter_data->scn && HIT_ON_BOUNDARY(hit, org, dir)) { 528 unsigned enc_ids[2] = {ENCLOSURE_ID_NULL, ENCLOSURE_ID_NULL}; 529 unsigned chk_enc_id = ENCLOSURE_ID_NULL; 530 531 scene_get_enclosure_ids(filter_data->scn, hit->prim.prim_id, enc_ids); 532 chk_enc_id = fX(dot)(dir, hit->normal) < 0 533 ? enc_ids[0] /* Front */ 534 : enc_ids[1]; /* Back */ 535 536 if(chk_enc_id != filter_data->enc_id) return 1; 537 } 538 return 0; 539 } 540 541 /* Retrieve the indices of `struct geometry' primitive */ 542 static void 543 XD(geometry_indices)(const unsigned iprim, unsigned out_ids[DIM], void* data) 544 { 545 struct geometry* ctx = data; 546 size_t ids[DIM]; 547 int i; 548 ASSERT(ctx && out_ids); 549 ctx->indices(iprim, ids, ctx->data); 550 FOR_EACH(i, 0, DIM) out_ids[i] = (unsigned)ids[i]; 551 } 552 553 /* Retrieve the coordinates of `struct geometry' vertex */ 554 static void 555 XD(geometry_position)(const unsigned ivert, double out_pos[DIM], void* data) 556 { 557 struct geometry* ctx = data; 558 double pos[DIM]; 559 int i; 560 ASSERT(ctx && out_pos); 561 ctx->position(ivert, pos, ctx->data); 562 FOR_EACH(i, 0, DIM) out_pos[i] = pos[i]; 563 } 564 565 /* Retrieve the indices of a primitive of a Star-EncXD descriptor */ 566 static void 567 XD(scene_indices)(const unsigned iprim, unsigned ids[DIM], void* data) 568 { 569 struct sencXd(scene)* scn = data; 570 SENCXD(scene_get_primitive(scn, iprim, ids)); 571 } 572 573 /* Retrieve the coordinates of a vertex of a Star-EncXD descriptor */ 574 static void 575 XD(scene_position)(const unsigned ivert, float out_pos[DIM], void* data) 576 { 577 struct sencXd(scene)* scn = data; 578 double pos[3]; 579 int i; 580 SENCXD(scene_get_vertex(scn, ivert, pos)); 581 FOR_EACH(i, 0, DIM) out_pos[i] = (float)pos[i]; 582 } 583 584 /* Retrieve the indices of a primitive of a Star-EncXD enclosure */ 585 static void 586 XD(enclosure_indices)(const unsigned iprim, unsigned ids[DIM], void* data) 587 { 588 struct sencXd(enclosure)* enc = data; 589 SENCXD(enclosure_get_primitive(enc, iprim, ids)); 590 } 591 592 /* Retrieve the coordinates of a vertex of a Star-EncXD encolsure */ 593 static void 594 XD(enclosure_position)(const unsigned ivert, float out_pos[DIM], void* data) 595 { 596 struct sencXd(enclosure)* enc = data; 597 double pos[DIM]; 598 int i; 599 ASSERT(out_pos); 600 SENCXD(enclosure_get_vertex(enc, ivert, pos)); 601 FOR_EACH(i, 0, DIM) out_pos[i] = (float)pos[i]; 602 } 603 604 /* Use Star-EncXD to analyze the user defined data. It essentially cleans-up 605 * the geometry and extracts the enclosures wrt to the submitted media. Note 606 * that data inconsistencies are also detected by this function. In this case 607 * an error is returned. */ 608 static res_T 609 XD(run_analyze) 610 (struct sdis_scene* scn, 611 const size_t nprims, /* #primitives */ 612 sdis_get_primitive_indices_T indices, 613 sdis_get_primitive_interface_T interf, 614 const size_t nverts, /* #vertices */ 615 sdis_get_vertex_position_T position, 616 void* ctx, 617 struct sencXd(scene)** out_scn) 618 { 619 struct geometry geom; 620 struct sencXd(device)* senc = NULL; 621 struct sencXd(scene)* senc_scn = NULL; 622 unsigned count; 623 res_T res = RES_OK; 624 ASSERT(scn && nprims && indices && interf && nverts && position && out_scn); 625 626 res = sencXd(device_create)(scn->dev->logger, scn->dev->allocator, 627 scn->dev->nthreads, scn->dev->verbose, &senc); 628 if(res != RES_OK) goto error; 629 630 /* Setup the geometry data */ 631 geom.indices = indices; 632 geom.interf = interf; 633 geom.position = position; 634 geom.data = ctx; 635 res = sencXd(scene_create)(senc, 636 SENCXD_(CONVENTION_NORMAL_BACK) | SENCXD_(CONVENTION_NORMAL_OUTSIDE), 637 (unsigned)nprims, XD(geometry_indices), geometry_media, 638 (unsigned)nverts, XD(geometry_position), &geom, &senc_scn); 639 if(res != RES_OK) goto error; 640 /* With il-formed scenes, scene creation can success without being able 641 * to extract enclosures; in this case just fail */ 642 res = sencXd(scene_get_enclosure_count(senc_scn, &count)); 643 if(res != RES_OK) { 644 count = 0; 645 #if DIM==2 646 senc2d_scene_get_overlapping_segments_count(senc_scn, &count); 647 if(count > 0) 648 log_err(scn->dev, 649 "%s: the scene includes overlapping segments.\n", 650 FUNC_NAME); 651 #else 652 senc3d_scene_get_overlapping_triangles_count(senc_scn, &count); 653 if(count > 0) 654 log_err(scn->dev, 655 "%s: the scene includes overlapping triangles.\n", 656 FUNC_NAME); 657 #endif 658 goto error; 659 } 660 661 exit: 662 if(senc) SENCXD(device_ref_put(senc)); 663 if(out_scn) *out_scn = senc_scn; 664 return res; 665 error: 666 if(senc_scn) { 667 SENCXD(scene_ref_put(senc_scn)); 668 senc_scn = NULL; 669 } 670 goto exit; 671 } 672 673 /* Register the media and the interfaces, map each primitive to its interface 674 * and associated to each primitive the identifier of the enclosures that it 675 * splits. */ 676 static res_T 677 XD(setup_properties) 678 (struct sdis_scene* scn, 679 struct sencXd(scene)* senc_scn, 680 void (*interf)(const size_t itri, struct sdis_interface**, void*), 681 void* ctx) 682 { 683 unsigned iprim, nprims; 684 res_T res = RES_OK; 685 ASSERT(scn && senc_scn && interf); 686 687 clear_properties(scn); 688 689 SENCXD(scene_get_primitives_count(senc_scn, &nprims)); 690 FOR_EACH(iprim, 0, nprims) { 691 struct prim_prop* prim_prop; 692 struct sdis_interface* itface; 693 unsigned enclosures[2]; 694 unsigned id; 695 int i; 696 double* enc_upper_bound; 697 size_t ninterfaces; 698 699 /* Fetch the enclosures that the segment/triangle splits */ 700 SENCXD(scene_get_primitive_enclosures(senc_scn, iprim, enclosures)); 701 702 /* Fetch the interface of the primitive */ 703 interf(iprim, &itface, ctx); 704 705 /* Check that the interface is already registered against the scene */ 706 id = interface_get_id(itface); 707 ninterfaces = darray_interf_size_get(&scn->interfaces); 708 if(id >= ninterfaces) { 709 res = darray_interf_resize(&scn->interfaces, id + 1); 710 if(res != RES_OK) goto error; 711 } 712 if(darray_interf_cdata_get(&scn->interfaces)[id]) { 713 ASSERT(darray_interf_cdata_get(&scn->interfaces)[id] == itface); 714 } else { 715 SDIS(interface_ref_get(itface)); 716 darray_interf_data_get(&scn->interfaces)[id] = itface; 717 } 718 719 /* Register the interface media against the scene */ 720 res = register_medium(scn, itface->medium_front); 721 if(res != RES_OK) goto error; 722 res = register_medium(scn, itface->medium_back); 723 if(res != RES_OK) goto error; 724 725 /* Allocate primitive properties */ 726 res = darray_prim_prop_resize(&scn->prim_props, iprim+1); 727 if(res != RES_OK) goto error; 728 729 /* Setup primitive properties */ 730 prim_prop = darray_prim_prop_data_get(&scn->prim_props) + iprim; 731 prim_prop->interf = itface; 732 prim_prop->front_enclosure = enclosures[0]; 733 prim_prop->back_enclosure = enclosures[1]; 734 735 /* Build per-interface hc upper bounds in a tmp table */ 736 FOR_EACH(i, 0, 2) { 737 double hc_ub = interface_get_convection_coef_upper_bound(itface); 738 enc_upper_bound = htable_d_find(&scn->tmp_hc_ub, enclosures+i); 739 if(!enc_upper_bound) { 740 res = htable_d_set(&scn->tmp_hc_ub, enclosures+i, &hc_ub); 741 } else { 742 *enc_upper_bound = MMAX(*enc_upper_bound, hc_ub); 743 } 744 } 745 } 746 747 exit: 748 return res; 749 error: 750 clear_properties(scn); 751 goto exit; 752 } 753 754 /* Build the Star-XD scene view of the whole scene */ 755 static res_T 756 XD(setup_scene_geometry)(struct sdis_scene* scn, struct sencXd(scene)* senc_scn) 757 { 758 struct sXd(device)* sXd_dev = NULL; 759 struct sXd(shape)* sXd_shape = NULL; 760 struct sXd(scene)* sXd_scn = NULL; 761 struct sXd(vertex_data) vdata = SXD_VERTEX_DATA_NULL; 762 unsigned nprims, nverts; 763 res_T res = RES_OK; 764 ASSERT(scn && senc_scn); 765 766 SENCXD(scene_get_vertices_count(senc_scn, &nverts)); 767 768 /* Setup the vertex data */ 769 vdata.usage = SXD_POSITION; 770 vdata.type = SXD_FLOATX; 771 vdata.get = XD(scene_position); 772 773 /* Create the Star-XD geometry of the whole scene */ 774 #define CALL(Func) { if(RES_OK != (res = Func)) goto error; } (void)0 775 sXd_dev = scn->dev->sXd(dev); 776 SENCXD(scene_get_primitives_count(senc_scn, &nprims)); 777 #if DIM == 2 778 CALL(sXd(shape_create_line_segments)(sXd_dev, &sXd_shape)); 779 CALL(sXd(line_segments_set_hit_filter_function)(sXd_shape, 780 XD(hit_filter_function), NULL)); 781 CALL(sXd(line_segments_setup_indexed_vertices)(sXd_shape, nprims, 782 XD(scene_indices), nverts, &vdata, 1, senc_scn)); 783 #else 784 CALL(sXd(shape_create_mesh)(sXd_dev, &sXd_shape)); 785 CALL(sXd(mesh_set_hit_filter_function)(sXd_shape, XD(hit_filter_function), NULL)); 786 CALL(sXd(mesh_setup_indexed_vertices)(sXd_shape, nprims, XD(scene_indices), 787 nverts, &vdata, 1, senc_scn)); 788 #endif 789 CALL(sXd(scene_create)(sXd_dev, &sXd_scn)); 790 CALL(sXd(scene_attach_shape)(sXd_scn, sXd_shape)); 791 CALL(sXd(scene_view_create)(sXd_scn, SXD_TRACE|SXD_GET_PRIMITIVE, 792 &scn->sXd(view))); 793 #undef CALL 794 795 exit: 796 if(sXd_shape) SXD(shape_ref_put(sXd_shape)); 797 if(sXd_scn) SXD(scene_ref_put(sXd_scn)); 798 return res; 799 error: 800 if(scn->sXd(view)) SXD(scene_view_ref_put(scn->sXd(view))); 801 goto exit; 802 } 803 804 /* Build the Star-XD scene view of a specific enclosure and map their local 805 * primitive id to their primitive id in the whole scene */ 806 static res_T 807 XD(register_enclosure)(struct sdis_scene* scn, struct sencXd(enclosure)* enc) 808 { 809 struct sXd(device)* sXd_dev = NULL; 810 struct sXd(scene)* sXd_scn = NULL; 811 struct sXd(shape)* sXd_shape = NULL; 812 struct sXd(vertex_data) vdata = SXD_VERTEX_DATA_NULL; 813 struct enclosure enc_dummy; 814 struct enclosure* enc_data; 815 float S, V; 816 double* p_ub; 817 unsigned iprim, nprims, nverts; 818 struct sencXd(enclosure_header) header; 819 res_T res = RES_OK; 820 ASSERT(scn && enc); 821 822 enclosure_init(scn->dev->allocator, &enc_dummy); 823 824 SENCXD(enclosure_get_header(enc, &header)); 825 sXd_dev = scn->dev->sXd(dev); 826 nprims = header.primitives_count; 827 nverts = header.vertices_count; 828 829 /* Register the enclosure into the scene. Use a dummy data on their 830 * registration; in order to avoid a costly copy, we are going to setup the 831 * registered data rather than a local data that would be then registered. In 832 * other words, the following hash table registration can be seen as an 833 * allocation of the enclosure data to setup. */ 834 res = htable_enclosure_set(&scn->enclosures, &header.enclosure_id, &enc_dummy); 835 if(res != RES_OK) goto error; 836 837 /* Fetch the data of the registered enclosure */ 838 enc_data = htable_enclosure_find(&scn->enclosures, &header.enclosure_id); 839 ASSERT(enc_data != NULL); 840 841 /* Setup the medium id of the enclosure */ 842 if(header.enclosed_media_count > 1) { 843 enc_data->medium_id = MEDIUM_ID_MULTI; 844 } else { 845 SENCXD(enclosure_get_medium(enc, 0, &enc_data->medium_id)); 846 } 847 848 /* Do not configure the enclosure geometry for enclosures that are infinite 849 * or composed of several media, i.e. that define boundary conditions */ 850 if(header.is_infinite || header.enclosed_media_count > 1) 851 goto exit; 852 853 /* Setup the vertex data */ 854 vdata.usage = SXD_POSITION; 855 vdata.type = SXD_FLOATX; 856 vdata.get = XD(enclosure_position); 857 858 /* Create the Star-XD geometry */ 859 #define CALL(Func) { if(RES_OK != (res = Func)) goto error; } (void)0 860 #if DIM == 2 861 CALL(sXd(shape_create_line_segments)(sXd_dev, &sXd_shape)); 862 CALL(sXd(line_segments_setup_indexed_vertices)(sXd_shape, nprims, 863 XD(enclosure_indices), nverts, &vdata, 1, enc)); 864 #else 865 CALL(sXd(shape_create_mesh)(sXd_dev, &sXd_shape)); 866 CALL(sXd(mesh_setup_indexed_vertices)(sXd_shape, nprims, XD(enclosure_indices), 867 nverts, &vdata, 1, enc)); 868 #endif 869 CALL(sXd(scene_create)(sXd_dev, &sXd_scn)); 870 CALL(sXd(scene_attach_shape)(sXd_scn, sXd_shape)); 871 CALL(sXd(scene_view_create)(sXd_scn, SXD_SAMPLE|SXD_TRACE, &enc_data->sXd(view))); 872 873 /* Compute the S/V ratio */ 874 #if DIM == 2 875 CALL(s2d_scene_view_compute_contour_length(enc_data->s2d_view, &S)); 876 CALL(s2d_scene_view_compute_area(enc_data->s2d_view, &V)); 877 #else 878 CALL(s3d_scene_view_compute_area(enc_data->s3d_view, &S)); 879 CALL(s3d_scene_view_compute_volume(enc_data->s3d_view, &V)); 880 #endif 881 enc_data->V = V; 882 enc_data->S_over_V = S / V; 883 ASSERT(enc_data->S_over_V >= 0); 884 #undef CALL 885 886 /* Set enclosure hc upper bound regardless of its media being a fluid */ 887 p_ub = htable_d_find(&scn->tmp_hc_ub, &header.enclosure_id); 888 ASSERT(p_ub); 889 enc_data->hc_upper_bound = *p_ub; 890 891 /* Define the identifier of the enclosure primitives in the whole scene */ 892 res = darray_uint_resize(&enc_data->local2global, nprims); 893 if(res != RES_OK) goto error; 894 FOR_EACH(iprim, 0, nprims) { 895 enum sencXd(side) side; 896 SENCXD(enclosure_get_primitive_id 897 (enc, iprim, darray_uint_data_get(&enc_data->local2global)+iprim, &side)); 898 } 899 900 exit: 901 enclosure_release(&enc_dummy); 902 if(sXd_shape) SXD(shape_ref_put(sXd_shape)); 903 if(sXd_scn) SXD(scene_ref_put(sXd_scn)); 904 return res; 905 error: 906 htable_enclosure_erase(&scn->enclosures, &header.enclosure_id); 907 goto exit; 908 } 909 910 /* Build the Star-XD scene view and define its associated data of the finite 911 * fluid enclosures */ 912 static res_T 913 XD(setup_enclosures)(struct sdis_scene* scn, struct sencXd(scene)* senc_scn) 914 { 915 struct sencXd(enclosure)* enc = NULL; 916 unsigned ienc, nencs; 917 int inner_multi = 0; 918 res_T res = RES_OK; 919 ASSERT(scn && senc_scn); 920 921 SENCXD(scene_get_enclosure_count(senc_scn, &nencs)); 922 FOR_EACH(ienc, 0, nencs) { 923 struct sencXd(enclosure_header) header; 924 925 SENCXD(scene_get_enclosure(senc_scn, ienc, &enc)); 926 SENCXD(enclosure_get_header(enc, &header)); 927 928 if(header.is_infinite) { 929 ASSERT(scn->outer_enclosure_id == UINT_MAX); /* Not set yet */ 930 scn->outer_enclosure_id = ienc; 931 } 932 933 if(header.enclosed_media_count != 1 && !header.is_infinite) 934 inner_multi++; 935 936 res = XD(register_enclosure)(scn, enc); 937 if(res != RES_OK) goto error; 938 939 SENCXD(enclosure_ref_put(enc)); 940 enc = NULL; 941 } 942 943 if(inner_multi) { 944 log_info(scn->dev, 945 "# Found %d internal enclosure(s) with more than 1 medium.\n", 946 inner_multi); 947 } 948 949 /* tmp table no more useful */ 950 htable_d_purge(&scn->tmp_hc_ub); 951 exit: 952 if(enc) SENCXD(enclosure_ref_put(enc)); 953 return res; 954 error: 955 goto exit; 956 } 957 958 #if DIM == 2 959 static res_T 960 setup_primitive_keys_2d(struct sdis_scene* scn, struct senc2d_scene* senc_scn) 961 { 962 unsigned iprim = 0; 963 unsigned nprims = 0; 964 res_T res = RES_OK; 965 ASSERT(scn && senc_scn); 966 967 SENC2D(scene_get_primitives_count(senc_scn, &nprims)); 968 969 FOR_EACH(iprim, 0, nprims) { 970 struct s2d_primitive prim = S2D_PRIMITIVE_NULL; 971 struct sdis_primkey key = SDIS_PRIMKEY_NULL; 972 unsigned ids[2] = {0,0}; 973 double v0[2] = {0,0}; 974 double v1[2] = {0,0}; 975 976 /* Retrieve positions from Star-Enclosre, not Star-2D. Star-Enclosure keeps 977 * the positions submitted by the user as they are, without any 978 * transformation or conversion (Star-2D converts them to float). This 979 * ensures that the caller can construct the same key from his data */ 980 SENC2D(scene_get_primitive(senc_scn, iprim, ids)); 981 SENC2D(scene_get_vertex(senc_scn, ids[0], v0)); 982 SENC2D(scene_get_vertex(senc_scn, ids[1], v1)); 983 S2D(scene_view_get_primitive(scn->s2d_view, iprim, &prim)); 984 985 sdis_primkey_2d_setup(&key, v0, v1); 986 987 res = htable_key2prim2d_set(&scn->key2prim2d, &key, &prim); 988 if(res != RES_OK) goto error; 989 } 990 991 exit: 992 return res; 993 error: 994 htable_key2prim2d_purge(&scn->key2prim2d); 995 goto exit; 996 } 997 998 #elif DIM == 3 999 static res_T 1000 setup_primitive_keys_3d(struct sdis_scene* scn, struct senc3d_scene* senc_scn) 1001 { 1002 unsigned iprim = 0; 1003 unsigned nprims = 0; 1004 res_T res = RES_OK; 1005 ASSERT(scn && senc_scn); 1006 1007 SENC3D(scene_get_primitives_count(senc_scn, &nprims)); 1008 1009 FOR_EACH(iprim, 0, nprims) { 1010 struct s3d_primitive prim = S3D_PRIMITIVE_NULL; 1011 struct sdis_primkey key = SDIS_PRIMKEY_NULL; 1012 unsigned ids[3] = {0}; 1013 double v0[3] = {0}; 1014 double v1[3] = {0}; 1015 double v2[3] = {0}; 1016 1017 /* Retrieve positions from Star-Enclosre, not Star-3D. Star-Enclosure keeps 1018 * the positions submitted by the user as they are, without any 1019 * transformation or conversion (Star-3D converts them to float). This 1020 * ensures that the caller can construct the same key from his data */ 1021 SENC3D(scene_get_primitive(senc_scn, iprim, ids)); 1022 SENC3D(scene_get_vertex(senc_scn, ids[0], v0)); 1023 SENC3D(scene_get_vertex(senc_scn, ids[1], v1)); 1024 SENC3D(scene_get_vertex(senc_scn, ids[2], v2)); 1025 S3D(scene_view_get_primitive(scn->s3d_view, iprim, &prim)); 1026 1027 sdis_primkey_setup(&key, v0, v1, v2); 1028 1029 res = htable_key2prim3d_set(&scn->key2prim3d, &key, &prim); 1030 if(res != RES_OK) goto error; 1031 } 1032 1033 exit: 1034 return res; 1035 error: 1036 htable_key2prim3d_purge(&scn->key2prim3d); 1037 goto exit; 1038 } 1039 #endif 1040 1041 /* Create a Stardis scene */ 1042 static res_T 1043 XD(scene_create) 1044 (struct sdis_device* dev, 1045 const struct sdis_scene_create_args* args, 1046 struct sdis_scene** out_scn) 1047 { 1048 struct sencXd(scene)* senc_scn = NULL; 1049 struct sdis_scene* scn = NULL; 1050 res_T res = RES_OK; 1051 1052 if(!dev || !check_sdis_scene_create_args(args) || !out_scn) { 1053 res = RES_BAD_ARG; 1054 goto error; 1055 } 1056 1057 scn = MEM_CALLOC(dev->allocator, 1, sizeof(struct sdis_scene)); 1058 if(!scn) { 1059 res = RES_MEM_ERR; 1060 log_err(dev, "%s: unabale to allocate the scene -- %s\n", 1061 FUNC_NAME, res_to_cstr(res)); 1062 goto error; 1063 } 1064 1065 ref_init(&scn->ref); 1066 SDIS(device_ref_get(dev)); 1067 scn->dev = dev; 1068 scn->fp_to_meter = args->fp_to_meter; 1069 scn->tmin = args->t_range[0]; 1070 scn->tmax = args->t_range[1]; 1071 scn->outer_enclosure_id = UINT_MAX; 1072 darray_interf_init(dev->allocator, &scn->interfaces); 1073 darray_medium_init(dev->allocator, &scn->media); 1074 darray_prim_prop_init(dev->allocator, &scn->prim_props); 1075 htable_enclosure_init(dev->allocator, &scn->enclosures); 1076 htable_d_init(dev->allocator, &scn->tmp_hc_ub); 1077 htable_key2prim2d_init(dev->allocator, &scn->key2prim2d); 1078 htable_key2prim3d_init(dev->allocator, &scn->key2prim3d); 1079 1080 if(args->source) { 1081 SDIS(source_ref_get(args->source)); 1082 scn->source = args->source; 1083 } 1084 1085 if(args->radenv) { 1086 SDIS(radiative_env_ref_get(args->radenv)); 1087 scn->radenv = args->radenv; 1088 } 1089 1090 res = XD(run_analyze) 1091 (scn, 1092 args->nprimitives, 1093 args->get_indices, 1094 args->get_interface, 1095 args->nvertices, 1096 args->get_position, 1097 args->context, 1098 &senc_scn); 1099 if(res != RES_OK) { 1100 log_err(dev, "%s: unable to analyze the scene -- %s\n", 1101 FUNC_NAME, res_to_cstr(res)); 1102 goto error; 1103 } 1104 res = XD(setup_properties)(scn, senc_scn, args->get_interface, args->context); 1105 if(res != RES_OK) { 1106 log_err(dev, "%s: unable to configure interfaces and media -- %s\n", 1107 FUNC_NAME, res_to_cstr(res)); 1108 goto error; 1109 } 1110 res = XD(setup_scene_geometry)(scn, senc_scn); 1111 if(res != RES_OK) { 1112 log_err(dev, "%s: unable to configure scene geometry -- %s\n", 1113 FUNC_NAME, res_to_cstr(res)); 1114 goto error; 1115 } 1116 res = XD(setup_enclosures)(scn, senc_scn); 1117 if(res != RES_OK) { 1118 log_err(dev, "%s: unable to configure enclosures -- %s\n", 1119 FUNC_NAME, res_to_cstr(res)); 1120 goto error; 1121 } 1122 res = XD(setup_primitive_keys)(scn, senc_scn); 1123 if(res != RES_OK) { 1124 log_err(dev, "%s: unable to configure primitive keys -- %s\n", 1125 FUNC_NAME, res_to_cstr(res)); 1126 goto error; 1127 } 1128 scn->sencXd(scn) = senc_scn; 1129 1130 exit: 1131 if(out_scn) *out_scn = scn; 1132 return res; 1133 error: 1134 if(senc_scn) SENCXD(scene_ref_put(senc_scn)); 1135 if(scn) { 1136 SDIS(scene_ref_put(scn)); 1137 scn = NULL; 1138 } 1139 goto exit; 1140 } 1141 1142 static res_T 1143 XD(scene_find_closest_point) 1144 (const struct sdis_scene* scn, 1145 const struct sdis_scene_find_closest_point_args* args, 1146 size_t* iprim, 1147 double uv[2]) 1148 { 1149 struct sXd(hit) hit; 1150 float query_pos[DIM]; 1151 float query_radius; 1152 res_T res = RES_OK; 1153 1154 if(!scn || !iprim || !uv || scene_is_2d(scn) != (DIM == 2)) { 1155 res = RES_BAD_ARG; 1156 goto error; 1157 } 1158 res = check_sdis_scene_find_closest_point_args(args); 1159 if(res != RES_OK) goto error; 1160 1161 /* Avoid a null query radius due to casting in single-precision */ 1162 query_radius = MMAX((float)args->radius, FLT_MIN); 1163 1164 fX_set_dX(query_pos, args->position); 1165 1166 /* Do not filter anything */ 1167 if(!args->XD(filter)) { 1168 res = sXd(scene_view_closest_point) 1169 (scn->sXd(view), query_pos, query_radius, NULL, &hit); 1170 1171 /* Filter points according to user-defined filter function */ 1172 } else { 1173 struct hit_filter_data filter_data = HIT_FILTER_DATA_NULL; 1174 filter_data.XD(custom_filter) = args->XD(filter); 1175 filter_data.custom_filter_data = args->filter_data; 1176 res = sXd(scene_view_closest_point) 1177 (scn->sXd(view), query_pos, query_radius, &filter_data, &hit); 1178 } 1179 1180 if(res != RES_OK) { 1181 log_err(scn->dev, 1182 "%s: error querying the closest position at `"FORMAT_VECX"' " 1183 "for a radius of %g -- %s.\n", 1184 FUNC_NAME, SPLITX(query_pos), query_radius, res_to_cstr(res)); 1185 goto error; 1186 } 1187 1188 if(SXD_HIT_NONE(&hit)) { 1189 *iprim = SDIS_PRIMITIVE_NONE; 1190 } else { 1191 *iprim = hit.prim.scene_prim_id; 1192 #if DIM == 2 1193 uv[0] = hit.u; 1194 #else 1195 uv[0] = hit.uv[0]; 1196 uv[1] = hit.uv[1]; 1197 #endif 1198 } 1199 1200 exit: 1201 return res; 1202 error: 1203 goto exit; 1204 } 1205 1206 /******************************************************************************* 1207 * Local functions 1208 ******************************************************************************/ 1209 static res_T 1210 XD(scene_get_enclosure_id) 1211 (struct sdis_scene* scn, 1212 const double pos[DIM], 1213 unsigned* out_enc_id) 1214 { 1215 size_t iprim, nprims; 1216 float P[DIM]; 1217 /* Range of the parametric coordinate into which positions are challenged */ 1218 #if DIM == 2 1219 float st[3]; 1220 #else 1221 float st[3][2]; 1222 #endif 1223 size_t nsteps = 3; 1224 unsigned enc_id = ENCLOSURE_ID_NULL; 1225 res_T res = RES_OK; 1226 ASSERT(scn && pos); 1227 1228 #if DIM == 2 1229 st[0] = 0.25f; 1230 st[1] = 0.50f; 1231 st[2] = 0.75f; 1232 #else 1233 f2(st[0], 1.f/6.f, 5.f/12.f); 1234 f2(st[1], 5.f/12.f, 1.f/6.f); 1235 f2(st[2], 5.f/12.f, 5.f/12.f); 1236 #endif 1237 1238 fX_set_dX(P, pos); 1239 1240 SXD(scene_view_primitives_count(scn->sXd(view), &nprims)); 1241 FOR_EACH(iprim, 0, nprims) { 1242 struct sXd(hit) hit; 1243 struct sXd(attrib) attr; 1244 struct sXd(primitive) prim; 1245 size_t iprim2; 1246 const float range[2] = {FLT_MIN, FLT_MAX}; 1247 float N[DIM] = {0}; 1248 float dir[DIM], cos_N_dir; 1249 size_t istep = 0; 1250 1251 /* 1 primitive over 2, take a primitive from the end of the primitive list. 1252 * When primitives are sorted in a coherent manner regarding their 1253 * position, this strategy avoids to test primitives that are going to be 1254 * rejected of the same manner due to possible numerical issues of the 1255 * resulting intersection. */ 1256 if((iprim % 2) == 0) { 1257 iprim2 = iprim / 2; 1258 } else { 1259 iprim2 = nprims - 1 - (iprim / 2); 1260 } 1261 1262 do { 1263 /* Retrieve a position onto the primitive */ 1264 SXD(scene_view_get_primitive(scn->sXd(view), (unsigned)iprim2, &prim)); 1265 SXD(primitive_get_attrib(&prim, SXD_POSITION, st[istep], &attr)); 1266 1267 /* Trace a ray from the random walk vertex toward the retrieved primitive 1268 * position */ 1269 fX(normalize)(dir, fX(sub)(dir, attr.value, P)); 1270 SXD(scene_view_trace_ray(scn->sXd(view), P, dir, range, NULL, &hit)); 1271 1272 /* Try another position onto the current primitive if there is no 1273 * intersection or if it is on a vertex/edge */ 1274 } while((SXD_HIT_NONE(&hit) || HIT_ON_BOUNDARY(&hit, P, dir)) 1275 && ++istep < nsteps); 1276 1277 /* No valid intersection is found on the current primitive. 1278 * Challenge another. */ 1279 if(istep >= nsteps) continue; 1280 1281 fX(normalize)(N, hit.normal); 1282 cos_N_dir = fX(dot)(N, dir); 1283 1284 /* Not too close and not roughly orthognonal */ 1285 if(hit.distance > 1.e-6 && absf(cos_N_dir) > 1.e-2f) { 1286 unsigned enc_ids[2]; 1287 scene_get_enclosure_ids(scn, hit.prim.prim_id, enc_ids); 1288 enc_id = cos_N_dir < 0 ? enc_ids[0] : enc_ids[1]; 1289 1290 break; /* That's all folks */ 1291 } 1292 } 1293 1294 if(iprim >= nprims) { 1295 log_warn(scn->dev, 1296 "%s: cannot retrieve current enclosure at {%g, %g, %g}.\n", 1297 FUNC_NAME, P[0], P[1], DIM == 3 ? P[2] : 0); 1298 res = RES_BAD_OP; 1299 goto error; 1300 } 1301 1302 if(iprim > 10 && iprim > (size_t)((double)nprims * 0.05)) { 1303 log_warn(scn->dev, 1304 "%s: performance issue. Up to %lu primitives were tested to find " 1305 "current enclosure at {%g, %g, %g}.\n", 1306 FUNC_NAME, (unsigned long)iprim, P[0], P[1], DIM == 3 ? P[2] : 0); 1307 } 1308 1309 exit: 1310 *out_enc_id = enc_id; 1311 return res; 1312 error: 1313 enc_id = ENCLOSURE_ID_NULL; 1314 goto exit; 1315 } 1316 1317 static res_T 1318 XD(scene_get_enclosure_id_in_closed_boundaries) 1319 (struct sdis_scene* scn, 1320 const double pos[DIM], 1321 unsigned* out_enc_id) 1322 { 1323 float dirs[6][3] = {{1,0,0},{-1,0,0},{0,1,0},{0,-1,0},{0,0,1},{0,0,-1}}; 1324 float frame[DIM*DIM]; 1325 float P[DIM]; 1326 unsigned enc_id = ENCLOSURE_ID_NULL; 1327 int idir; 1328 res_T res = RES_OK; 1329 ASSERT(scn && pos && out_enc_id); 1330 1331 /* Build a frame that will be used to rotate the main axis by PI/4 around 1332 * each axis. This can avoid numerical issues when geometry is discretized 1333 * along the main axis */ 1334 #if DIM == 2 1335 f22_rotation(frame, (float)PI/4); 1336 #else 1337 f33_rotation(frame, (float)PI/4, (float)PI/4, (float)PI/4); 1338 #endif 1339 1340 fX_set_dX(P, pos); 1341 FOR_EACH(idir, 0, 2*DIM) { 1342 struct sXd(hit) hit; 1343 float N[DIM] = {0}; 1344 const float range[2] = {FLT_MIN, FLT_MAX}; 1345 float cos_N_dir; 1346 1347 /* Transform the directions to avoid to be aligned with the axis */ 1348 fXX_mulfX(dirs[idir], frame, dirs[idir]); 1349 1350 /* Trace a ray from the random walk vertex toward the retrieved primitive 1351 * position */ 1352 SXD(scene_view_trace_ray(scn->sXd(view), P, dirs[idir], range, NULL, &hit)); 1353 1354 /* Unforeseen error. One has to intersect a primitive ! */ 1355 if(SXD_HIT_NONE(&hit)) continue; 1356 1357 /* Discard a hits if it lies on an edge/point */ 1358 if(HIT_ON_BOUNDARY(&hit, P, dirs[idir])) continue; 1359 1360 fX(normalize)(N, hit.normal); 1361 cos_N_dir = fX(dot)(N, dirs[idir]); 1362 1363 /* Not too close and not roughly orthogonal */ 1364 if(hit.distance > 1.e-6 && absf(cos_N_dir) > 1.e-2f) { 1365 unsigned enc_ids[2]; 1366 scene_get_enclosure_ids(scn, hit.prim.prim_id, enc_ids); 1367 enc_id = cos_N_dir < 0 ? enc_ids[0] : enc_ids[1]; 1368 1369 break; /* That's all folks */ 1370 } 1371 } 1372 if(idir >= 2*DIM) { /* Fallback to scene_get_enclosure_id function */ 1373 res = XD(scene_get_enclosure_id)(scn, pos, &enc_id); 1374 if(res != RES_OK) goto error; 1375 } 1376 1377 exit: 1378 *out_enc_id = enc_id; 1379 return res; 1380 error: 1381 enc_id = ENCLOSURE_ID_NULL; 1382 goto exit; 1383 } 1384 1385 #if (SDIS_XD_DIMENSION == 2) 1386 #include <star/sencX2d_undefs.h> 1387 #else /* SDIS_XD_DIMENSION == 3 */ 1388 #include <star/sencX3d_undefs.h> 1389 #endif 1390 1391 #undef HIT_ON_BOUNDARY 1392 1393 #include "sdis_Xd_end.h"