star-gs

Literate program for a geometric sensitivity calculation
git clone git://git.meso-star.fr/star-gs.git
Log | Files | Refs | README | LICENSE

sgs_geometry.c (15698B)


      1 /* Copyright (C) 2021-2023 Centre National de la Recherche Scientifique
      2  * Copyright (C) 2021-2023 INSA Lyon
      3  * Copyright (C) 2021-2023 Institut Mines Télécom Albi-Carmaux
      4  * Copyright (C) 2021-2023 |Méso|Star> (contact@meso-star.com)
      5  * Copyright (C) 2021-2023 Institut Pascal
      6  * Copyright (C) 2021-2023 PhotonLyX (info@photonlyx.com)
      7  * Copyright (C) 2021-2023 Université de Lorraine
      8  * Copyright (C) 2021-2023 Université Paul Sabatier
      9  * Copyright (C) 2021-2023 Université Toulouse - Jean Jaurès
     10  *
     11  * This program is free software: you can redistribute it and/or modify
     12  * it under the terms of the GNU General Public License as published by
     13  * the Free Software Foundation, either version 3 of the License, or
     14  * (at your option) any later version.
     15  *
     16  * This program is distributed in the hope that it will be useful,
     17  * but WITHOUT ANY WARRANTY; without even the implied warranty of
     18  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
     19  * GNU General Public License for more details.
     20  *
     21  * You should have received a copy of the GNU General Public License
     22  * along with this program. If not, see <http://www.gnu.org/licenses/>. */
     23 
     24 #include "sgs.h"
     25 #include "sgs_geometry.h"
     26 #include "sgs_geometry_c.h"
     27 #include "sgs_log.h"
     28 
     29 #include <star/s3d.h>
     30 #include <star/ssp.h>
     31 
     32 #include <rsys/cstr.h>
     33 #include <rsys/float2.h>
     34 #include <rsys/float3.h>
     35 
     36 struct hit_filter_context {
     37   float range[2];
     38   enum sgs_surface_type surface_from;
     39   const int* prim2surface; /* Map a primitive id to its surface type */
     40 };
     41 #define HIT_FILTER_CONTEXT_NULL__ {{0,FLT_MAX}, SGS_SURFACE_NONE, NULL}
     42 static const struct hit_filter_context HIT_FILTER_CONTEXT_NULL =
     43   HIT_FILTER_CONTEXT_NULL__;
     44 
     45 struct mesh {
     46   const double* vertices;
     47   const size_t* triangles;
     48   size_t nvertices;
     49   size_t ntriangles;
     50 };
     51 static const struct mesh MESH_NULL = {NULL, NULL, 0, 0};
     52 
     53 /*******************************************************************************
     54  * Helper functions
     55  ******************************************************************************/
     56 static int
     57 hit_filter
     58   (const struct s3d_hit* hit,
     59    const float ray_org[3],
     60    const float ray_dir[3],
     61    const float ray_range[2],
     62    void* ray_data,
     63    void* filter_data)
     64 {
     65   const struct hit_filter_context* ctx = ray_data;
     66   enum sgs_surface_type surface;
     67   (void)ray_org, (void)ray_dir, (void)ray_range, (void)filter_data;
     68 
     69   if(!ctx) /* Nothing to do */
     70     return 0;
     71 
     72   /* Internally, Star-3D relies on Embree that, due to numerically inaccuracy,
     73    * can find hits whose distances are not strictly included in the submitted
     74    * ray range. Discard these hits. */
     75   if(hit->distance <= ctx->range[0] || hit->distance >= ctx->range[1])
     76     return 1;
     77 
     78   /* The hit surface is the one from which the the ray starts. Discard it */
     79   surface = ctx->prim2surface[hit->prim.prim_id];
     80   if(surface == ctx->surface_from)
     81     return 1;
     82 
     83   return 0; /* No filtering */
     84 }
     85 
     86 static void
     87 mesh_get_vertex(const unsigned ivert, float pos[3], void* ctx)
     88 {
     89   struct mesh* mesh = ctx;
     90   ASSERT(pos && ctx && ivert < mesh->nvertices);
     91   pos[0] = (float)mesh->vertices[ivert*3+0];
     92   pos[1] = (float)mesh->vertices[ivert*3+1];
     93   pos[2] = (float)mesh->vertices[ivert*3+2];
     94 }
     95 static void
     96 mesh_get_triangle(const unsigned itri, unsigned ids[3], void* ctx)
     97 {
     98   struct mesh* mesh = ctx;
     99   ASSERT(ids && ctx && itri < mesh->ntriangles);
    100   ids[0] = (unsigned)mesh->triangles[itri*3+0];
    101   ids[1] = (unsigned)mesh->triangles[itri*3+1];
    102   ids[2] = (unsigned)mesh->triangles[itri*3+2];
    103 }
    104 
    105 static res_T
    106 mesh_create_view
    107   (struct sgs* sgs,
    108    struct mesh* mesh,
    109    const int view_mask,
    110    struct s3d_scene_view** out_view)
    111 {
    112   struct s3d_vertex_data vdata = S3D_VERTEX_DATA_NULL;
    113   struct s3d_shape* shape = NULL;
    114   struct s3d_scene* scene = NULL;
    115   struct s3d_scene_view* view = NULL;
    116   res_T res = RES_OK;
    117   ASSERT(sgs && mesh && out_view);
    118 
    119   /* Setup the Star-3D vertex data */
    120   vdata.usage = S3D_POSITION;
    121   vdata.type = S3D_FLOAT3;
    122   vdata.get = mesh_get_vertex;
    123 
    124   /* Setup the Star-3D ray-tracing scene */
    125   res = s3d_shape_create_mesh(sgs_get_s3d_device(sgs), &shape);
    126   if(res != RES_OK) goto error;
    127   res = s3d_mesh_setup_indexed_vertices(shape, (unsigned)mesh->ntriangles,
    128     mesh_get_triangle, (unsigned)mesh->nvertices, &vdata, 1, mesh);
    129   if(res != RES_OK) goto error;
    130   res = s3d_mesh_set_hit_filter_function(shape, hit_filter, NULL);
    131   if(res != RES_OK) goto error;
    132   res = s3d_scene_create(sgs_get_s3d_device(sgs), &scene);
    133   if(res != RES_OK) goto error;
    134   res = s3d_scene_attach_shape(scene, shape);
    135   if(res != RES_OK) goto error;
    136 
    137   /* Create the scene view */
    138   res = s3d_scene_view_create(scene, view_mask, &view);
    139   if(res != RES_OK) goto error;
    140 
    141 exit:
    142   if(shape) S3D(shape_ref_put(shape));
    143   if(scene) S3D(scene_ref_put(scene));
    144   *out_view = view;
    145   return res;
    146 error:
    147   if(view) {
    148     S3D(scene_view_ref_put(view));
    149     view = NULL;
    150   }
    151   goto exit;
    152 }
    153 
    154 static void
    155 release_geometry(ref_T* ref)
    156 {
    157   struct sgs_geometry* geom = CONTAINER_OF(ref, struct sgs_geometry, ref);
    158   struct sgs* sgs = NULL;
    159   ASSERT(ref);
    160 
    161   if(geom->view_sp) S3D(scene_view_ref_put(geom->view_sp));
    162   if(geom->view_rt) S3D(scene_view_ref_put(geom->view_rt));
    163   darray_double_release(&geom->verts);
    164   darray_size_t_release(&geom->tris);
    165   darray_int_release(&geom->tris_rt_type);
    166   darray_int_release(&geom->tris_sp_type);
    167 
    168   sgs = geom->sgs;
    169   MEM_RM(sgs_get_allocator(sgs), geom);
    170 }
    171 
    172 /*******************************************************************************
    173  * API functions
    174  ******************************************************************************/
    175 void
    176 sgs_geometry_ref_get(struct sgs_geometry* geom)
    177 {
    178   ASSERT(geom);
    179   ref_get(&geom->ref);
    180 }
    181 
    182 void
    183 sgs_geometry_ref_put(struct sgs_geometry* geom)
    184 {
    185   ASSERT(geom);
    186   ref_put(&geom->ref, release_geometry);
    187 }
    188 
    189 void
    190 sgs_geometry_get_aabb
    191   (const struct sgs_geometry* geom,
    192    double low[3],
    193    double upp[3])
    194 {
    195   ASSERT(geom && low && upp);
    196   d3_set(low, geom->low);
    197   d3_set(upp, geom->upp);
    198 }
    199 
    200 int
    201 sgs_geometry_get_sampling_mask(const struct sgs_geometry* geom)
    202 {
    203   ASSERT(geom);
    204   return geom->sampling_mask;
    205 }
    206 
    207 double
    208 sgs_geometry_compute_sampling_area(const struct sgs_geometry* geom)
    209 {
    210   float area;
    211   ASSERT(geom);
    212   S3D(scene_view_compute_area(geom->view_sp, &area));
    213   return area;
    214 }
    215 
    216 void
    217 sgs_geometry_trace_ray
    218   (const struct sgs_geometry* geom,
    219    const double ray_org[3],
    220    const double ray_dir[3],
    221    const double ray_range[2],
    222    const enum sgs_surface_type surface_from,
    223    struct sgs_hit* out_hit)
    224 {
    225   struct hit_filter_context filter_ctx = HIT_FILTER_CONTEXT_NULL;
    226   struct s3d_hit hit = S3D_HIT_NULL;
    227   float org[3];
    228   float dir[3];
    229   float range[2];
    230   ASSERT(geom && ray_org && ray_dir && ray_range && out_hit);
    231 
    232   /* Convert ray data to single precision */
    233   f3_set_d3(org, ray_org);
    234   f3_set_d3(dir, ray_dir);
    235   f2_set_d2(range, ray_range);
    236 
    237   /* Ensure the normalization of the ray direction after its conversion to
    238    * single precision */
    239   f3_normalize(dir, dir);
    240 
    241   /* Setup the hit filter data */
    242   filter_ctx.range[0] = range[0];
    243   filter_ctx.range[1] = range[1];
    244   filter_ctx.surface_from = surface_from;
    245   filter_ctx.prim2surface = darray_int_cdata_get(&geom->tris_rt_type);
    246 
    247   /* Trace the ray */
    248   S3D(scene_view_trace_ray(geom->view_rt, org, dir, range, &filter_ctx, &hit));
    249 
    250   /* Setup the output hit */
    251   *out_hit = SGS_HIT_NULL;
    252   if(!S3D_HIT_NONE(&hit)) {
    253     d3_set_f3(out_hit->normal, hit.normal);
    254     d3_normalize(out_hit->normal, out_hit->normal);
    255     out_hit->distance = hit.distance;
    256     out_hit->surface = darray_int_cdata_get(&geom->tris_rt_type)[hit.prim.prim_id];
    257   }
    258 }
    259 
    260 extern LOCAL_SYM void
    261 sgs_geometry_sample_sensitivity_source
    262   (const struct sgs_geometry* geom,
    263    struct ssp_rng* rng,
    264    struct sgs_fragment* frag)
    265 {
    266   struct s3d_primitive prim = S3D_PRIMITIVE_NULL;
    267   struct s3d_attrib position;
    268   struct s3d_attrib normal;
    269   float r0, r1, r2;
    270   float st[2];
    271   ASSERT(geom && rng && frag);
    272 
    273   /* Sample the geometry */
    274   r0 = ssp_rng_canonical_float(rng);
    275   r1 = ssp_rng_canonical_float(rng);
    276   r2 = ssp_rng_canonical_float(rng);
    277   S3D(scene_view_sample(geom->view_sp, r0, r1, r2, &prim, st));
    278 
    279   /* Retrieve sampling point attributes */
    280   S3D(primitive_get_attrib(&prim, S3D_POSITION, st, &position));
    281   S3D(primitive_get_attrib(&prim, S3D_GEOMETRY_NORMAL, st, &normal));
    282 
    283   /* Setup the output variable */
    284   d3_set_f3(frag->position, position.value);
    285   d3_set_f3(frag->normal, normal.value);
    286   d3_normalize(frag->normal, frag->normal);
    287   frag->surface = darray_int_cdata_get(&geom->tris_sp_type)[prim.prim_id];
    288 }
    289 
    290 res_T
    291 sgs_geometry_dump_vtk(const struct sgs_geometry* geom, FILE* stream)
    292 {
    293   size_t i = 0;
    294   size_t nverts = 0;
    295   size_t ntris = 0;
    296   res_T res = RES_OK;
    297   int err = 0;
    298   ASSERT(geom);
    299 
    300   nverts = geometry_get_nvertices(geom);
    301   ntris = geometry_get_ntriangles(geom);
    302 
    303   #define FPRINTF(Fmt, Args) {                                                 \
    304     err = fprintf(stream, Fmt COMMA_##Args LIST_##Args);                       \
    305     if(err < 0) {                                                              \
    306       sgs_log_err(geom->sgs, "Error writing data.\n");                         \
    307       res = RES_IO_ERR;                                                        \
    308       goto error;                                                              \
    309     }                                                                          \
    310   } (void)0
    311 
    312   /* Write header */
    313   FPRINTF("# vtk DataFile Version 2.0\n", ARG0());
    314   FPRINTF("Geometry\n", ARG0());
    315   FPRINTF("ASCII\n", ARG0());
    316   FPRINTF("DATASET POLYDATA\n", ARG0());
    317 
    318   /* Write the vertices */
    319   FPRINTF("POINTS %lu float\n", ARG1((unsigned long)nverts));
    320   FOR_EACH(i, 0, nverts) {
    321     const double* pos = darray_double_cdata_get(&geom->verts) + i*3;
    322     FPRINTF("%g %g %g\n", ARG3(pos[0], pos[1], pos[2]));
    323   }
    324 
    325   /* Write the triangles */
    326   FPRINTF("POLYGONS %lu %lu\n",
    327     ARG2((unsigned long)ntris, (unsigned long)ntris*4));
    328   FOR_EACH(i, 0, ntris) {
    329     const size_t* ids = darray_size_t_cdata_get(&geom->tris) + i*3;
    330     FPRINTF("3 %lu %lu %lu\n", ARG3(ids[0], ids[1], ids[2]));
    331   }
    332 
    333   FPRINTF("CELL_DATA %lu\n", ARG1((unsigned long)ntris));
    334 
    335   /* Write the triangle type */
    336   FPRINTF("SCALARS Triangle_Type integer 1\n", ARG0());
    337   FPRINTF("LOOKUP_TABLE default\n", ARG0());
    338   FOR_EACH(i, 0, ntris) {
    339     const enum sgs_surface_type type =
    340       darray_int_cdata_get(&geom->tris_rt_type)[i];
    341     FPRINTF("%i\n", ARG1(type));
    342   }
    343 
    344   /* Write the sampling type */
    345   FPRINTF("SCALARS Sampling_Surface integer 1\n", ARG0());
    346   FPRINTF("LOOKUP_TABLE Sampling_Type\n", ARG0());
    347   FOR_EACH(i, 0, ntris) {
    348     const int type = darray_int_cdata_get(&geom->tris_rt_type)[i];
    349     FPRINTF("%i\n", ARG1((type & geom->sampling_mask) != 0));
    350   }
    351   FPRINTF("LOOKUP_TABLE Sampling_Type 2\n", ARG0());
    352   FPRINTF("0 0 1 1\n", ARG0());
    353   FPRINTF("1 0 0 1\n", ARG0());
    354 
    355   #undef FPRINTF
    356 
    357 exit:
    358   return res;
    359 error:
    360   goto exit;
    361 }
    362 
    363 /*******************************************************************************
    364  * Local functions
    365  ******************************************************************************/
    366 res_T
    367 geometry_create
    368   (struct sgs* sgs,
    369    const int sampling_mask,
    370    struct sgs_geometry** out_geom)
    371 {
    372   struct sgs_geometry* geom = NULL;
    373   struct mem_allocator* allocator = NULL;
    374   res_T res = RES_OK;
    375   ASSERT(sgs && out_geom);
    376   ASSERT(sampling_mask != 0);
    377 
    378   allocator = sgs_get_allocator(sgs);
    379   geom = MEM_CALLOC(allocator, 1, sizeof(struct sgs_geometry));
    380   if(!geom) {
    381     sgs_log_err(sgs, "Could not allocate the geometry.\n");
    382     res = RES_MEM_ERR;
    383     goto error;
    384   }
    385   ref_init(&geom->ref);
    386   geom->sgs = sgs;
    387   darray_double_init(allocator, &geom->verts);
    388   darray_size_t_init(allocator, &geom->tris);
    389   darray_int_init(allocator, &geom->tris_rt_type);
    390   darray_int_init(allocator, &geom->tris_sp_type);
    391   geom->sampling_mask = sampling_mask;
    392 
    393 exit:
    394   *out_geom = geom;
    395   return res;
    396 error:
    397   if(geom) {
    398     sgs_geometry_ref_put(geom);
    399     geom = NULL;
    400   }
    401   goto exit;
    402 }
    403 
    404 void
    405 geometry_compute_aabb(struct sgs_geometry* geom)
    406 {
    407   const double* coords = NULL;
    408   size_t i;
    409   ASSERT(geom);
    410   d3_splat(geom->low, DBL_MAX);
    411   d3_splat(geom->upp,-DBL_MAX);
    412 
    413   coords = darray_double_cdata_get(&geom->verts);
    414   FOR_EACH(i, 0, darray_double_size_get(&geom->verts)/3) {
    415     const double* vert = coords + i*3;
    416     geom->low[0] = MMIN(geom->low[0], vert[0]);
    417     geom->low[1] = MMIN(geom->low[1], vert[1]);
    418     geom->low[2] = MMIN(geom->low[2], vert[2]);
    419     geom->upp[0] = MMAX(geom->upp[0], vert[0]);
    420     geom->upp[1] = MMAX(geom->upp[1], vert[1]);
    421     geom->upp[2] = MMAX(geom->upp[2], vert[2]);
    422   }
    423 }
    424 
    425 res_T
    426 geometry_setup_view_rt(struct sgs_geometry* geom)
    427 {
    428   struct mesh mesh = MESH_NULL;
    429   res_T res = RES_OK;
    430   ASSERT(geom);
    431 
    432   /* Clean up an eventual previous scene view */
    433   if(geom->view_rt) {
    434     S3D(scene_view_ref_put(geom->view_rt));
    435     geom->view_rt = NULL;
    436   }
    437 
    438   /* Setup the mesh used during the ray-tracing */
    439   mesh.vertices = darray_double_cdata_get(&geom->verts);
    440   mesh.triangles = darray_size_t_cdata_get(&geom->tris);
    441   mesh.nvertices = darray_double_size_get(&geom->verts)/3/*#coords per vertex*/;
    442   mesh.ntriangles = darray_size_t_size_get(&geom->tris)/3/*#ids per triangle*/;
    443 
    444   /* Create the scene view */
    445   res = mesh_create_view(geom->sgs, &mesh, S3D_TRACE, &geom->view_rt);
    446   if(res != RES_OK) {
    447     sgs_log_err(geom->sgs,
    448       "Error creating the Star-3D scene view to ray-trace -- %s.\n",
    449       res_to_cstr(res));
    450     goto error;
    451   }
    452 
    453 exit:
    454   return res;
    455 error:
    456   if(geom->view_rt) {
    457     S3D(scene_view_ref_put(geom->view_rt));
    458     geom->view_rt = NULL;
    459   }
    460   goto exit;
    461 }
    462 
    463 res_T
    464 geometry_setup_view_sp(struct sgs_geometry* geom, const int sampling_mask)
    465 {
    466   struct darray_size_t tris_sp;
    467   struct mesh mesh = MESH_NULL;
    468   const int* types = NULL;
    469   const size_t* ids = NULL;
    470   size_t itri = 0;
    471   size_t ntris = 0;
    472   res_T res = RES_OK;
    473   ASSERT(geom && sampling_mask);
    474 
    475   darray_size_t_init(sgs_get_allocator(geom->sgs), &tris_sp);
    476 
    477   /* Clean up an eventual previous scene view */
    478   if(geom->view_sp) {
    479     S3D(scene_view_ref_put(geom->view_sp));
    480     geom->view_sp = NULL;
    481   }
    482 
    483   /* Fetch the geometry data */
    484   ntris = geometry_get_ntriangles(geom);
    485   ids = darray_size_t_cdata_get(&geom->tris);
    486   types = darray_int_cdata_get(&geom->tris_rt_type);
    487 
    488   /* Find the triangles to sample */
    489   FOR_EACH(itri, 0, ntris) {
    490     if((types[itri] & sampling_mask) != 0) {
    491 
    492       /* Register the triangle as a triangle to sample */
    493       res = darray_size_t_push_back(&tris_sp, &ids[itri*3+0]);
    494       if(res != RES_OK) goto error;
    495       res = darray_size_t_push_back(&tris_sp, &ids[itri*3+1]);
    496       if(res != RES_OK) goto error;
    497       res = darray_size_t_push_back(&tris_sp, &ids[itri*3+2]);
    498       if(res != RES_OK) goto error;
    499 
    500       /* Save the surface type of the triangle */
    501       res = darray_int_push_back(&geom->tris_sp_type, &types[itri]);
    502       if(res != RES_OK) goto error;
    503     }
    504   }
    505 
    506   /* Setup the mesh to sample */
    507   mesh.vertices = darray_double_cdata_get(&geom->verts);
    508   mesh.triangles = darray_size_t_cdata_get(&tris_sp);
    509   mesh.nvertices = darray_double_size_get(&geom->verts)/3/*#coords per vertex*/;
    510   mesh.ntriangles = darray_size_t_size_get(&tris_sp)/3/*#ids per triangle*/;
    511 
    512   /* Create the scene view */
    513   res = mesh_create_view(geom->sgs, &mesh, S3D_SAMPLE, &geom->view_sp);
    514   if(res != RES_OK) {
    515     sgs_log_err(geom->sgs,
    516       "Error creating the Star-3D scene view to sample -- %s.\n",
    517       res_to_cstr(res));
    518     goto error;
    519   }
    520 
    521 exit:
    522   darray_size_t_release(&tris_sp);
    523   return res;
    524 error:
    525   if(geom->view_sp) {
    526     S3D(scene_view_ref_put(geom->view_sp));
    527     geom->view_sp = NULL;
    528   }
    529   goto exit;
    530 }