schiff_mesh.c (18401B)
1 /* Copyright (C) 2015-2016 CNRS 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 #include "schiff_mesh.h" 17 18 #include <rsys/float2.h> 19 #include <rsys/float3.h> 20 #include <rsys/stretchy_array.h> 21 22 /******************************************************************************* 23 * Helper functions 24 ******************************************************************************/ 25 /* Assume that the vertices are arranged in "phi major" order and that the 26 * lower and upper polar vertices are the two last ones. */ 27 static void 28 setup_sphere_indices 29 (struct darray_uint* indices, 30 const unsigned nthetas, 31 const unsigned nphis) 32 { 33 unsigned i, itheta, iphi; 34 ASSERT(indices && nthetas && nphis); 35 36 /* Define the indices of the contour primitives */ 37 i = 0; 38 FOR_EACH(itheta, 0, nthetas) { 39 const unsigned itheta0 = itheta * (nphis - 1); 40 const unsigned itheta1 = ((itheta + 1) % nthetas) * (nphis - 1); 41 FOR_EACH(iphi, 0, nphis-2) { 42 const unsigned iphi0 = iphi + 0; 43 const unsigned iphi1 = iphi + 1; 44 45 darray_uint_data_get(indices)[i++] = itheta0 + iphi0; 46 darray_uint_data_get(indices)[i++] = itheta0 + iphi1; 47 darray_uint_data_get(indices)[i++] = itheta1 + iphi0; 48 49 darray_uint_data_get(indices)[i++] = itheta1 + iphi0; 50 darray_uint_data_get(indices)[i++] = itheta0 + iphi1; 51 darray_uint_data_get(indices)[i++] = itheta1 + iphi1; 52 } 53 } 54 55 /* Define the indices of the polar primitives */ 56 FOR_EACH(itheta, 0, nthetas) { 57 const unsigned itheta0 = itheta * (nphis - 1); 58 const unsigned itheta1 = ((itheta + 1) % nthetas) * (nphis - 1); 59 60 darray_uint_data_get(indices)[i++] = nthetas * (nphis - 1); 61 darray_uint_data_get(indices)[i++] = itheta0; 62 darray_uint_data_get(indices)[i++] = itheta1; 63 64 darray_uint_data_get(indices)[i++] = nthetas * (nphis - 1) + 1; 65 darray_uint_data_get(indices)[i++] = itheta1 + (nphis - 2); 66 darray_uint_data_get(indices)[i++] = itheta0 + (nphis - 2); 67 } 68 } 69 70 /******************************************************************************* 71 * Local functions 72 ******************************************************************************/ 73 res_T 74 schiff_mesh_init_sphere 75 (struct mem_allocator* allocator, 76 struct schiff_mesh* sphere, 77 const unsigned nthetas) 78 { 79 /* Theta in [0, 2PI[ and phi in [0, PI[ */ 80 const unsigned nphis = (unsigned)(((double)nthetas + 0.5) / 2.0); 81 const double step_theta = 2*PI / (double)nthetas; 82 const double step_phi = PI / (double)nphis; 83 size_t nverts; 84 size_t ntris; 85 size_t i; 86 unsigned itheta, iphi; 87 res_T res = RES_OK; 88 ASSERT(allocator && sphere && nthetas); 89 90 sphere->coordinates = SCHIFF_CARTESIAN; 91 darray_float_init(allocator, &sphere->vertices.cartesian); 92 darray_uint_init(allocator, &sphere->indices); 93 darray_sincos_init(allocator, &sphere->thetas); 94 darray_sincos_init(allocator, &sphere->phis); 95 96 nverts = nthetas*(nphis-1)/* #contour verts */ + 2/* polar verts */; 97 ntris = 2*nthetas*(nphis-2)/* #contour tris */ + 2*nthetas/* #polar tris */; 98 99 res = darray_float_resize(&sphere->vertices.cartesian, nverts * 3/*#coords*/); 100 if(res != RES_OK) goto error; 101 res = darray_uint_resize(&sphere->indices, ntris * 3/*#indices per tri*/); 102 if(res != RES_OK) goto error; 103 res = darray_sincos_resize(&sphere->thetas, nthetas); 104 if(res != RES_OK) goto error; 105 res = darray_sincos_resize(&sphere->phis, nphis); 106 if(res != RES_OK) goto error; 107 108 /* Precompute the cosine/sinus of the theta/phi angles */ 109 FOR_EACH(itheta, 0, nthetas) { 110 const double theta = -PI + (double)itheta * step_theta; 111 darray_sincos_data_get(&sphere->thetas)[itheta].angle = theta; 112 darray_sincos_data_get(&sphere->thetas)[itheta].sinus = (float)sin(theta); 113 darray_sincos_data_get(&sphere->thetas)[itheta].cosine = (float)cos(theta); 114 } 115 FOR_EACH(iphi, 0, nphis-1) { 116 const double phi = -PI/2 + (double)(iphi + 1) * step_phi; 117 darray_sincos_data_get(&sphere->phis)[iphi].angle = (float)phi; 118 darray_sincos_data_get(&sphere->phis)[iphi].sinus = (float)sin(phi); 119 darray_sincos_data_get(&sphere->phis)[iphi].cosine = (float)cos(phi); 120 } 121 122 /* Build the contour vertices */ 123 i = 0; 124 FOR_EACH(itheta, 0, nthetas) { 125 const struct sin_cos* theta = darray_sincos_data_get(&sphere->thetas) + itheta; 126 FOR_EACH(iphi, 0, nphis-1) { 127 const struct sin_cos* phi = darray_sincos_data_get(&sphere->phis) + iphi; 128 darray_float_data_get(&sphere->vertices.cartesian)[i++] = 129 (float)(theta->cosine * phi->cosine); 130 darray_float_data_get(&sphere->vertices.cartesian)[i++] = 131 (float)(theta->sinus * phi->cosine); 132 darray_float_data_get(&sphere->vertices.cartesian)[i++] = 133 (float)phi->sinus; 134 } 135 } 136 /* Setup polar vertices */ 137 f3(darray_float_data_get(&sphere->vertices.cartesian) + i+0, 0.f, 0.f,-1.f); 138 f3(darray_float_data_get(&sphere->vertices.cartesian) + i+3, 0.f, 0.f, 1.f); 139 140 /* Define the indices of the sphere */ 141 setup_sphere_indices(&sphere->indices, nthetas, nphis); 142 143 sphere->is_init = 1; 144 exit: 145 darray_sincos_release(&sphere->thetas); 146 darray_sincos_release(&sphere->phis); 147 return res; 148 error: 149 schiff_mesh_release(sphere); 150 goto exit; 151 } 152 153 res_T 154 schiff_mesh_init_sphere_polar 155 (struct mem_allocator* allocator, 156 struct schiff_mesh* sphere, 157 const unsigned nthetas) 158 { 159 /* Theta in [0, 2PI[ and phi in [0, PI[ */ 160 const unsigned nphis = (unsigned)(((double)nthetas + 0.5) / 2.0); 161 const double step_theta = 2*PI / (double)nthetas; 162 const double step_phi = PI / (double)nphis; 163 size_t nverts; 164 size_t ntris; 165 size_t i; 166 unsigned itheta, iphi; 167 res_T res = RES_OK; 168 ASSERT(allocator && sphere && nthetas); 169 170 sphere->coordinates = SCHIFF_POLAR; 171 darray_uint_init(allocator, &sphere->vertices.polar); 172 darray_uint_init(allocator, &sphere->indices); 173 darray_sincos_init(allocator, &sphere->thetas); 174 darray_sincos_init(allocator, &sphere->phis); 175 176 nverts = nthetas*(nphis-1)/* #contour verts */ + 2/* polar verts */; 177 ntris = 2*nthetas*(nphis-2)/* #contour tris */ + 2*nthetas/* #polar tris */; 178 179 res = darray_uint_resize(&sphere->vertices.polar, nverts * 2/* theta phi*/); 180 if(res != RES_OK) goto error; 181 res = darray_uint_resize(&sphere->indices, ntris * 3/*#indices per tri*/); 182 if(res != RES_OK) goto error; 183 res = darray_sincos_resize(&sphere->thetas, nthetas); 184 if(res != RES_OK) goto error; 185 res = darray_sincos_resize(&sphere->phis, nphis + 1/*Polar*/); 186 if(res != RES_OK) goto error; 187 188 /* Cosine/Sinus of the theta/phi contour angles */ 189 FOR_EACH(itheta, 0, nthetas) { 190 const double theta = -PI + (double)itheta * step_theta; 191 darray_sincos_data_get(&sphere->thetas)[itheta].angle = theta; 192 darray_sincos_data_get(&sphere->thetas)[itheta].sinus = (float)sin(theta); 193 darray_sincos_data_get(&sphere->thetas)[itheta].cosine = (float)cos(theta); 194 } 195 FOR_EACH(iphi, 0, nphis-1) { 196 const double phi = -PI/2 + (double)(iphi + 1) * step_phi; 197 darray_sincos_data_get(&sphere->phis)[iphi].angle = (float)phi; 198 darray_sincos_data_get(&sphere->phis)[iphi].sinus = (float)sin(phi); 199 darray_sincos_data_get(&sphere->phis)[iphi].cosine = (float)cos(phi); 200 } 201 202 /* Cosine/Sinus of the theta/phi polar angles */ 203 darray_sincos_data_get(&sphere->phis)[nphis-1].angle =-(float)PI/2; 204 darray_sincos_data_get(&sphere->phis)[nphis-1].sinus =-1.f; 205 darray_sincos_data_get(&sphere->phis)[nphis-1].cosine = 0.f; 206 darray_sincos_data_get(&sphere->phis)[nphis-0].angle = (float)PI/2; 207 darray_sincos_data_get(&sphere->phis)[nphis-0].sinus = 1.f; 208 darray_sincos_data_get(&sphere->phis)[nphis-0].cosine = 0.f; 209 210 /* Build the contour vertices */ 211 i = 0; 212 FOR_EACH(itheta, 0, nthetas) { 213 FOR_EACH(iphi, 0, nphis-1) { 214 darray_uint_data_get(&sphere->vertices.polar)[i++] = itheta; 215 darray_uint_data_get(&sphere->vertices.polar)[i++] = iphi; 216 } 217 } 218 /* Setup polar vertices */ 219 darray_uint_data_get(&sphere->vertices.polar)[i++] = 0; 220 darray_uint_data_get(&sphere->vertices.polar)[i++] = nphis-1; 221 darray_uint_data_get(&sphere->vertices.polar)[i++] = 0; 222 darray_uint_data_get(&sphere->vertices.polar)[i++] = nphis-0; 223 224 /* Define the indices of the sphere */ 225 setup_sphere_indices(&sphere->indices, nthetas, nphis); 226 227 sphere->is_init = 1; 228 exit: 229 return res; 230 error: 231 schiff_mesh_release(sphere); 232 goto exit; 233 } 234 235 res_T 236 schiff_mesh_init_cylinder 237 (struct mem_allocator* allocator, 238 struct schiff_mesh* cylinder, 239 const unsigned nsteps) 240 { 241 const double step = 2*PI / (double)nsteps; 242 size_t nverts; 243 size_t ntris; 244 unsigned i; 245 res_T res = RES_OK; 246 ASSERT(allocator && cylinder && nsteps); 247 248 cylinder->coordinates = SCHIFF_CARTESIAN; 249 darray_float_init(allocator, &cylinder->vertices.cartesian); 250 darray_uint_init(allocator, &cylinder->indices); 251 252 nverts = nsteps*2/* #contour verts */ + 2/* #polar verts */; 253 ntris = nsteps*2/* #contour tris */ + 2*nsteps/* #cap tris */; 254 255 res = darray_float_resize(&cylinder->vertices.cartesian, nverts*3/*#coords*/); 256 if(res != RES_OK) goto error; 257 res = darray_uint_resize(&cylinder->indices, ntris*3/*#indices per tri*/); 258 if(res != RES_OK) goto error; 259 260 /* Generate the vertex coordinates */ 261 FOR_EACH(i, 0, nsteps) { 262 const float theta = (float)(i* step); 263 const float x = (float)cos(theta); 264 const float y = (float)sin(theta); 265 f3(darray_float_data_get(&cylinder->vertices.cartesian)+(i*2+0)*3, x,y,0); 266 f3(darray_float_data_get(&cylinder->vertices.cartesian)+(i*2+1)*3, x,y,1); 267 } 268 269 /* "Polar" vertices */ 270 f3(darray_float_data_get(&cylinder->vertices.cartesian)+(i*2+0)*3, 0,0,0); 271 f3(darray_float_data_get(&cylinder->vertices.cartesian)+(i*2+1)*3, 0,0,1); 272 273 /* Contour primitives */ 274 FOR_EACH(i, 0, nsteps) { 275 const unsigned id = i * 2; 276 unsigned* iprim = darray_uint_data_get(&cylinder->indices) + id*3; 277 278 iprim[0] = (id + 0); 279 iprim[1] = (id + 1); 280 iprim[2] = (id + 2) % (nsteps*2); 281 282 iprim += 3; 283 284 iprim[0] = (id + 2) % (nsteps*2); 285 iprim[1] = (id + 1); 286 iprim[2] = (id + 3) % (nsteps*2); 287 } 288 289 /* Cap primitives */ 290 FOR_EACH(i, 0, nsteps) { 291 const unsigned id = i* 2; 292 unsigned* iprim = darray_uint_data_get(&cylinder->indices) + (nsteps + i)*6; 293 294 iprim[0] = (nsteps * 2); 295 iprim[1] = (id + 0); 296 iprim[2] = (id + 2) % (nsteps*2); 297 298 iprim += 3; 299 300 iprim[0] = (nsteps * 2) + 1; 301 iprim[1] = (id + 3) % (nsteps*2); 302 iprim[2] = (id + 1); 303 } 304 305 cylinder->is_init = 1; 306 exit: 307 return res; 308 error: 309 schiff_mesh_release(cylinder); 310 goto exit; 311 } 312 313 res_T 314 schiff_mesh_init_helical_pipe 315 (struct mem_allocator* allocator, 316 struct schiff_mesh* helical_pipe, 317 const unsigned nsteps_helicoid, 318 const unsigned nsteps_circle) 319 320 { 321 unsigned* indices; 322 unsigned* bottom_cap; 323 unsigned* top_cap; 324 size_t ihelicoid, icircle; 325 size_t nverts, ntris; 326 size_t ibottom_cap_index, itop_cap_index; 327 size_t ibottom_cap_vertex, itop_cap_vertex; 328 size_t ilast_meridian_vertex; 329 res_T res; 330 ASSERT(allocator && helical_pipe && nsteps_helicoid>=2 && nsteps_circle>=4); 331 332 helical_pipe->coordinates = SCHIFF_NO_COORDINATE; 333 darray_uint_init(allocator, &helical_pipe->indices); 334 335 nverts = (nsteps_helicoid+1)*nsteps_circle/*#contour*/ + 2/*#cap*/; 336 ntris = (nsteps_helicoid)*nsteps_circle*2/*#contour*/ + 2*nsteps_circle/*#cap*/; 337 338 res = darray_uint_resize(&helical_pipe->indices, ntris*3/*#indices per tri*/); 339 if(res != RES_OK) goto error; 340 341 indices = darray_uint_data_get(&helical_pipe->indices); 342 343 /* Setup the indices toward the contour primitives */ 344 FOR_EACH(ihelicoid, 0, nsteps_helicoid) { 345 /* Offset toward the first index of the current meridian */ 346 const size_t offset = ihelicoid * nsteps_circle; 347 348 FOR_EACH(icircle, 0, nsteps_circle) { 349 size_t prim_ids[4]; 350 const size_t id = icircle + offset; 351 unsigned* iprim = indices + id*3/*#ids per prim*/*2/*#prims per step*/; 352 353 /* Indices of the quad */ 354 prim_ids[0] = icircle + offset; 355 prim_ids[1] = icircle + nsteps_circle + offset; 356 prim_ids[2] = (icircle + 1) % nsteps_circle + nsteps_circle + offset; 357 prim_ids[3] = (icircle + 1) % nsteps_circle + offset; 358 /* First triangle of the quad */ 359 iprim[0] = (unsigned)prim_ids[0]; 360 iprim[1] = (unsigned)prim_ids[3]; 361 iprim[2] = (unsigned)prim_ids[1]; 362 /* Second triangule of the quad */ 363 iprim[3] = (unsigned)prim_ids[1]; 364 iprim[4] = (unsigned)prim_ids[3]; 365 iprim[5] = (unsigned)prim_ids[2]; 366 } 367 } 368 369 /* Define the index of the cap vertices */ 370 ibottom_cap_vertex = nverts - 2; /* Index toward the bottom cap vertex */ 371 itop_cap_vertex = nverts - 1; /* Index toward the top cap vertex */ 372 373 /* Index of the first vertex of the last meridian circle */ 374 ilast_meridian_vertex = nsteps_helicoid * nsteps_circle; 375 376 ibottom_cap_index = /* Index of the 1st bottom cap index */ 377 nsteps_helicoid 378 * nsteps_circle 379 * 2 /*#prims per step*/ 380 * 3 /*#ids per prim*/; 381 itop_cap_index = /* Index of the 1st top cap index */ 382 ibottom_cap_index 383 + nsteps_circle * 3/*# ids per primitive */; 384 385 bottom_cap = indices + ibottom_cap_index ; 386 top_cap = indices + itop_cap_index; 387 388 FOR_EACH(icircle, 0, nsteps_circle) { 389 const size_t id = icircle * 3/*#ids per prim*/; 390 391 bottom_cap[id + 0] = (unsigned)icircle; 392 bottom_cap[id + 1] = (unsigned)ibottom_cap_vertex; 393 bottom_cap[id + 2] = (unsigned)((icircle + 1) % nsteps_circle); 394 395 top_cap[id + 0] = (unsigned)((icircle+1)%nsteps_circle + ilast_meridian_vertex); 396 top_cap[id + 1] = (unsigned)itop_cap_vertex; 397 top_cap[id + 2] = (unsigned)(icircle + ilast_meridian_vertex); 398 } 399 400 helical_pipe->is_init = 1; 401 exit: 402 return res; 403 error: 404 schiff_mesh_release(helical_pipe); 405 goto exit; 406 } 407 408 res_T 409 schiff_mesh_helical_pipe_create_vertices 410 (const struct schiff_mesh* mesh, 411 const double pitch, 412 const double height, 413 const double hradius, /* Helicoid radius */ 414 const double cradius, /* Circle radius */ 415 const unsigned nsteps_helicoid, 416 const unsigned nsteps_circle, 417 size_t* nvertices, 418 float** out_vertices) 419 { 420 double* meridian = NULL; 421 float* vertices = NULL; 422 double c; 423 double phi_max; 424 double step_helicoid, step_circle; 425 double rcp_sqrt_hradius2_add_c2; /* Precomputed value */ 426 size_t ihelicoid, icircle; 427 size_t icoord, icoord_top, icoord_bottom; 428 size_t nverts; 429 res_T res = RES_OK; 430 431 ASSERT(mesh && pitch > 0 && height > 0 && hradius > 0 && cradius > 0); 432 ASSERT(nvertices && out_vertices); 433 ASSERT(mesh->coordinates == SCHIFF_NO_COORDINATE); 434 ASSERT(nsteps_helicoid * nsteps_circle * 2 * 3 + nsteps_circle*3*2/*cap*/ 435 == darray_uint_size_get(&mesh->indices)); 436 (void)mesh; 437 438 nverts = (nsteps_helicoid+1) * nsteps_circle/*#contour*/ + 2/*#cap*/; 439 if(!sa_add(vertices, nverts * 3/*#coords per vertex*/)) { 440 res = RES_MEM_ERR; 441 goto error; 442 } 443 if(!sa_add(meridian, nsteps_circle * 3/*#coords per vertex*/)) { 444 res = RES_MEM_ERR; 445 goto error; 446 } 447 448 c = pitch / 2*PI; 449 phi_max = height * 2*PI / pitch; 450 step_helicoid = phi_max / (double)nsteps_helicoid; 451 step_circle = 2*PI / (double)nsteps_circle; 452 rcp_sqrt_hradius2_add_c2 = 1.0 / sqrt(hradius*hradius + c*c); 453 454 /* Compute the meridian positions */ 455 FOR_EACH(icircle, 0, nsteps_circle) { 456 const double theta = (double)icircle * step_circle; 457 const double cos_theta = cos(theta); 458 const double sin_theta = sin(theta); 459 const size_t id = icircle *3; 460 461 meridian[id + 0] = cradius * cos_theta + hradius; 462 meridian[id + 1] = -cradius * c * rcp_sqrt_hradius2_add_c2 * sin_theta; 463 meridian[id + 2] = cradius*hradius * rcp_sqrt_hradius2_add_c2 * sin_theta; 464 } 465 466 icoord = 0; 467 icoord_bottom = 468 (nsteps_helicoid + 1) 469 * nsteps_circle 470 * 3/*#coords per vertex*/; 471 icoord_top = icoord_bottom + 3/*#coords per vertex */; 472 473 FOR_EACH(ihelicoid, 0, nsteps_helicoid + 1) { 474 const double phi = (double)ihelicoid * step_helicoid; 475 const double cos_phi = cos(phi); 476 const double sin_phi = sin(phi); 477 478 FOR_EACH(icircle, 0, nsteps_circle) { 479 const double* pos = meridian + icircle*3; 480 vertices[icoord + 0] = (float)(pos[0]*cos_phi - pos[1]*sin_phi); 481 vertices[icoord + 1] = (float)(pos[0]*sin_phi + pos[1]*cos_phi); 482 vertices[icoord + 2] = (float)(pos[2] + c * phi); 483 icoord += 3; 484 } 485 486 /* Top cap vertex */ 487 if(ihelicoid == 0) { 488 vertices[icoord_bottom + 0] = (float)(hradius * cos_phi); 489 vertices[icoord_bottom + 1] = (float)(hradius * sin_phi); 490 vertices[icoord_bottom + 2] = (float)(c*phi); 491 492 /* Bottom cap vertex */ 493 } else if(ihelicoid == nsteps_helicoid) { 494 vertices[icoord_top + 0] = (float)(hradius * cos_phi); 495 vertices[icoord_top + 1] = (float)(hradius * sin_phi); 496 vertices[icoord_top + 2] = (float)(c * phi); 497 } 498 } 499 500 exit: 501 if(meridian) sa_release(meridian); 502 *nvertices = nverts; 503 *out_vertices = vertices; 504 return res; 505 error: 506 if(vertices) { 507 sa_release(vertices); 508 vertices = NULL; 509 nverts = 0; 510 } 511 goto exit; 512 } 513 514 void 515 schiff_mesh_helical_pipe_destroy_vertices 516 (const struct schiff_mesh* mesh, float* out_vertices) 517 { 518 ASSERT(mesh && out_vertices); 519 (void)mesh; 520 sa_release(out_vertices); 521 } 522 523 void 524 schiff_mesh_release(struct schiff_mesh* mesh) 525 { 526 ASSERT(mesh); 527 if(!mesh->is_init) 528 return; /* The mesh is not initialised */ 529 530 switch(mesh->coordinates) { 531 case SCHIFF_CARTESIAN: 532 darray_float_release(&mesh->vertices.cartesian); 533 break; 534 case SCHIFF_POLAR: 535 darray_uint_release(&mesh->vertices.polar); 536 darray_sincos_release(&mesh->thetas); 537 darray_sincos_release(&mesh->phis); 538 break; 539 case SCHIFF_NO_COORDINATE: /* Do nothing */ break; 540 default: FATAL("Unreachable code\n"); break; 541 } 542 darray_uint_release(&mesh->indices); 543 mesh->coordinates = SCHIFF_NO_COORDINATE; 544 mesh->is_init = 0; 545 } 546