test_sschiff_utils.h (10796B)
1 /* Copyright (C) 2015, 2016, 2026 Centre National de la Recherche Scientifique 2 * Copyright (C) 2026 Clermont Auvergne INP 3 * Copyright (C) 2026 Institut Mines Télécom Albi-Carmaux 4 * Copyright (C) 2020, 2021, 2023, 2026 |Méso|Star> (contact@meso-star.com) 5 * Copyright (C) 2026 Université de Lorraine 6 * Copyright (C) 2026 Université de Toulouse 7 * 8 * This program is free software: you can redistribute it and/or modify 9 * it under the terms of the GNU General Public License as published by 10 * the Free Software Foundation, either version 3 of the License, or 11 * (at your option) any later version. 12 * 13 * This program is distributed in the hope that it will be useful, 14 * but WITHOUT ANY WARRANTY; without even the implied warranty of 15 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 16 * GNU General Public License for more details. 17 * 18 * You should have received a copy of the GNU General Public License 19 * along with this program. If not, see <http://www.gnu.org/licenses/>. */ 20 21 #ifndef TEST_SSCHIFF_UTILS_H 22 #define TEST_SSCHIFF_UTILS_H 23 24 #include <rsys/float3.h> 25 #include <rsys/stretchy_array.h> 26 #include <rsys/mem_allocator.h> 27 28 #include <star/s3d.h> 29 30 /******************************************************************************* 31 * Helper geometry data structure 32 ******************************************************************************/ 33 struct geometry { 34 float* vertices; /* List of float[3] */ 35 unsigned* indices; /* List of unsigned[3] */ 36 }; 37 38 static INLINE void 39 geometry_init_sphere(struct geometry* sphere, const unsigned ntheta) 40 { 41 const unsigned nphi = (unsigned)(((double)ntheta + 0.5) / 2.0); 42 const double step_theta = 2*PI / (double)ntheta; 43 const double step_phi = PI / (double)nphi; 44 double* cos_theta = NULL; 45 double* sin_theta = NULL; 46 double* cos_phi = NULL; 47 double* sin_phi = NULL; 48 unsigned itheta, iphi; 49 50 CHK(sphere != NULL); 51 CHK(ntheta != 0); 52 53 sphere->vertices = NULL; 54 sphere->indices = NULL; 55 56 /* Precompute the cosine/sinus of the theta/phi angles */ 57 FOR_EACH(itheta, 0, ntheta) { 58 const double theta = -PI + (double)itheta * step_theta; 59 sa_push(cos_theta, cos(theta)); 60 sa_push(sin_theta, sin(theta)); 61 } 62 FOR_EACH(iphi, 1, nphi) { 63 const double phi = -PI/2 + (double)iphi * step_phi; 64 sa_push(cos_phi, cos(phi)); 65 sa_push(sin_phi, sin(phi)); 66 } 67 /* Compute the contour vertices */ 68 FOR_EACH(itheta, 0, ntheta) { 69 FOR_EACH(iphi, 0, nphi-1) { 70 sa_push(sphere->vertices, (float)(cos_theta[itheta]*cos_phi[iphi])); 71 sa_push(sphere->vertices, (float)(sin_theta[itheta]*cos_phi[iphi])); 72 sa_push(sphere->vertices, (float)sin_phi[iphi]); 73 } 74 } 75 /* Compute the polar vertices */ 76 f3(sa_add(sphere->vertices, 3), 0.f, 0.f,-1.f); 77 f3(sa_add(sphere->vertices, 3), 0.f, 0.f, 1.f); 78 79 /* Define the indices of the contour primitives */ 80 FOR_EACH(itheta, 0, ntheta) { 81 const unsigned itheta0 = itheta * (nphi - 1); 82 const unsigned itheta1 = ((itheta + 1) % ntheta) * (nphi - 1); 83 FOR_EACH(iphi, 0, nphi-2) { 84 const unsigned iphi0 = iphi + 0; 85 const unsigned iphi1 = iphi + 1; 86 87 sa_push(sphere->indices, itheta0 + iphi0); 88 sa_push(sphere->indices, itheta0 + iphi1); 89 sa_push(sphere->indices, itheta1 + iphi0); 90 91 sa_push(sphere->indices, itheta1 + iphi0); 92 sa_push(sphere->indices, itheta0 + iphi1); 93 sa_push(sphere->indices, itheta1 + iphi1); 94 } 95 } 96 97 /* Define the indices of the polar primitives */ 98 FOR_EACH(itheta, 0, ntheta) { 99 const unsigned itheta0 = itheta * (nphi - 1); 100 const unsigned itheta1 = ((itheta + 1) % ntheta) * (nphi - 1); 101 102 sa_push(sphere->indices, ntheta * (nphi - 1)); 103 sa_push(sphere->indices, itheta0); 104 sa_push(sphere->indices, itheta1); 105 106 sa_push(sphere->indices, ntheta * (nphi - 1) + 1); 107 sa_push(sphere->indices, itheta1 + (nphi - 2)); 108 sa_push(sphere->indices, itheta0 + (nphi - 2)); 109 } 110 111 /* Release the intermediary data structure */ 112 sa_release(cos_theta); 113 sa_release(sin_theta); 114 sa_release(cos_phi); 115 sa_release(sin_phi); 116 } 117 118 static INLINE void 119 geometry_init_cylinder(struct geometry* geometry, const unsigned nsteps) 120 { 121 const double step = 2*PI / (double)nsteps; 122 unsigned istep; 123 124 CHK(geometry != NULL); 125 CHK(nsteps != 0); 126 127 geometry->vertices = NULL; 128 geometry->indices = NULL; 129 130 /* Generate the vertex coordinates */ 131 FOR_EACH(istep, 0, nsteps) { 132 const float theta = (float)(istep * step); 133 const float x = (float)cos(theta); 134 const float y = (float)sin(theta); 135 f3(sa_add(geometry->vertices, 3), x, y, -1.f); 136 f3(sa_add(geometry->vertices, 3), x, y, 0.f); 137 } 138 139 /* "Polar" vertices */ 140 f3(sa_add(geometry->vertices, 3), 0.f, 0.f, -1.f); 141 f3(sa_add(geometry->vertices, 3), 0.f, 0.f, 0.f); 142 143 /* Contour primitives */ 144 FOR_EACH(istep, 0, nsteps) { 145 const unsigned id = istep * 2; 146 unsigned* iprim; 147 148 iprim = sa_add(geometry->indices, 3); 149 iprim[0] = (id + 0); 150 iprim[1] = (id + 1); 151 iprim[2] = (id + 2) % (nsteps*2); 152 153 iprim = sa_add(geometry->indices, 3); 154 iprim[0] = (id + 2) % (nsteps*2); 155 iprim[1] = (id + 1); 156 iprim[2] = (id + 3) % (nsteps*2); 157 } 158 159 /* Cap primitives */ 160 FOR_EACH(istep, 0, nsteps) { 161 const unsigned id = istep * 2; 162 unsigned* iprim; 163 164 iprim = sa_add(geometry->indices, 3); 165 iprim[0] = (nsteps * 2); 166 iprim[1] = (id + 0); 167 iprim[2] = (id + 2) % (nsteps*2); 168 169 iprim = sa_add(geometry->indices, 3); 170 iprim[0] = (nsteps * 2) + 1; 171 iprim[1] = (id + 3) % (nsteps*2); 172 iprim[2] = (id + 1); 173 } 174 175 } 176 177 static INLINE void 178 geometry_release(struct geometry* geometry) 179 { 180 CHK(geometry != NULL); 181 sa_release(geometry->vertices); 182 sa_release(geometry->indices); 183 geometry->vertices = NULL; 184 geometry->indices = NULL; 185 } 186 187 static INLINE void 188 geometry_dump(struct geometry* geom, FILE* file) 189 { 190 size_t i; 191 CHK(geom != NULL); 192 CHK(file != NULL); 193 194 CHK(sa_size(geom->vertices)%3 == 0); /* Ensure 3D position */ 195 CHK(sa_size(geom->indices)%3 == 0); /* Ensure triangular primitives */ 196 197 FOR_EACH(i, 0, sa_size(geom->vertices)/3) { 198 fprintf(file, "v %f %f %f\n", 199 geom->vertices[i*3+0], 200 geom->vertices[i*3+1], 201 geom->vertices[i*3+2]); 202 } 203 FOR_EACH(i, 0, sa_size(geom->indices)/3) { 204 fprintf(file, "f %d %d %d\n", 205 geom->indices[i*3+0] + 1, 206 geom->indices[i*3+1] + 1, 207 geom->indices[i*3+2] + 1); 208 } 209 } 210 211 /******************************************************************************* 212 * Cylinder shape 213 ******************************************************************************/ 214 struct cylinder { 215 struct geometry* geometry; 216 float radius; 217 float height; 218 }; 219 220 static INLINE void 221 cylinder_get_indices(const unsigned itri, unsigned ids[3], void* ctx) 222 { 223 struct cylinder* cylinder = ctx; 224 const size_t i = itri * 3; 225 226 CHK(sa_size(cylinder->geometry->indices) % 3 == 0); 227 CHK(itri < sa_size(cylinder->geometry->indices) / 3); 228 ids[0] = cylinder->geometry->indices[i + 0]; 229 ids[1] = cylinder->geometry->indices[i + 1]; 230 ids[2] = cylinder->geometry->indices[i + 2]; 231 } 232 233 static INLINE void 234 cylinder_get_position(const unsigned ivert, float vertex[3], void* ctx) 235 { 236 struct cylinder* cylinder = ctx; 237 const size_t i = ivert * 3; 238 239 CHK(sa_size(cylinder->geometry->vertices) % 3 == 0); 240 CHK(ivert < sa_size(cylinder->geometry->vertices) / 3); 241 vertex[0] = cylinder->geometry->vertices[i + 0] * cylinder->radius; 242 vertex[1] = cylinder->geometry->vertices[i + 1] * cylinder->radius; 243 vertex[2] = cylinder->geometry->vertices[i + 2] * cylinder->height; 244 } 245 246 static INLINE res_T 247 cylinder_setup_s3d_shape(struct cylinder* cylinder, struct s3d_shape* shape) 248 { 249 struct s3d_vertex_data attrib; 250 size_t nverts, nprims; 251 252 CHK(cylinder != NULL); 253 CHK(shape != NULL); 254 255 attrib.usage = S3D_POSITION; 256 attrib.type = S3D_FLOAT3; 257 attrib.get = cylinder_get_position; 258 259 nverts = sa_size(cylinder->geometry->vertices) / 3/*#coords*/; 260 nprims = sa_size(cylinder->geometry->indices) / 3/*#indices per prim*/; 261 262 return s3d_mesh_setup_indexed_vertices(shape, (unsigned)nprims, 263 cylinder_get_indices, (unsigned)nverts, &attrib, 1, cylinder); 264 } 265 266 /******************************************************************************* 267 * Spherical shape 268 ******************************************************************************/ 269 struct sphere { 270 struct geometry* geometry; 271 float radius; 272 }; 273 274 static INLINE void 275 sphere_get_indices(const unsigned itri, unsigned ids[3], void* ctx) 276 { 277 struct sphere* sphere = ctx; 278 const size_t i = itri * 3; 279 280 CHK(sa_size(sphere->geometry->indices) % 3 == 0); 281 CHK(itri < sa_size(sphere->geometry->indices) / 3); 282 ids[0] = sphere->geometry->indices[i + 0]; 283 ids[1] = sphere->geometry->indices[i + 1]; 284 ids[2] = sphere->geometry->indices[i + 2]; 285 } 286 287 static INLINE void 288 sphere_get_position(const unsigned ivert, float vertex[3], void* ctx) 289 { 290 struct sphere* sphere = ctx; 291 const size_t i = ivert * 3; 292 293 CHK(sa_size(sphere->geometry->vertices) % 3 == 0); 294 CHK(ivert < sa_size(sphere->geometry->vertices) / 3); 295 vertex[0] = sphere->geometry->vertices[i + 0] * sphere->radius; 296 vertex[1] = sphere->geometry->vertices[i + 1] * sphere->radius; 297 vertex[2] = sphere->geometry->vertices[i + 2] * sphere->radius; 298 } 299 300 static INLINE res_T 301 sphere_setup_s3d_shape(struct sphere* sphere, struct s3d_shape* shape) 302 { 303 struct s3d_vertex_data attrib; 304 size_t nverts, nprims; 305 306 CHK(sphere != NULL); 307 CHK(shape != NULL); 308 309 attrib.usage = S3D_POSITION; 310 attrib.type = S3D_FLOAT3; 311 attrib.get = sphere_get_position; 312 313 nverts = sa_size(sphere->geometry->vertices) / 3/*#coords*/; 314 nprims = sa_size(sphere->geometry->indices) / 3/*#indices per prim*/; 315 316 return s3d_mesh_setup_indexed_vertices(shape, (unsigned)nprims, 317 sphere_get_indices, (unsigned)nverts, &attrib, 1, sphere); 318 } 319 320 /******************************************************************************* 321 * Miscellaneous functions 322 ******************************************************************************/ 323 static INLINE void 324 compute_estimation_intersection 325 (double intersection[2], 326 const double scale, 327 const struct sschiff_state* state0, 328 const struct sschiff_state* state1) 329 { 330 double interval0[2], interval1[2]; 331 CHK(scale > 0); 332 interval0[0] = state0->E - scale*state0->SE; 333 interval0[1] = state0->E + scale*state0->SE; 334 interval1[0] = state1->E - scale*state1->SE; 335 interval1[1] = state1->E + scale*state1->SE; 336 intersection[0] = MMAX(interval0[0], interval1[0]); 337 intersection[1] = MMIN(interval0[1], interval1[1]); 338 } 339 340 static void 341 check_memory_allocator(struct mem_allocator* allocator) 342 { 343 if(MEM_ALLOCATED_SIZE(allocator)) { 344 char dump[512]; 345 MEM_DUMP(allocator, dump, sizeof(dump)/sizeof(char)); 346 fprintf(stderr, "%s\n", dump); 347 FATAL("Memory leaks\n"); 348 } 349 } 350 351 #endif /* TEST_SSCHIFF_UTILS_H */ 352