git.mutantstargoat.com Git - summerhack/blob - src/3dengfx/src/3dengfx/ggen.cpp

   1 /*
   2 This file is part of the 3dengfx, 3d visualization system.
   3
   4 Copyright (c) 2004, 2005, 2006 John Tsiombikas <nuclear@siggraph.org>
   5
   6 3dengfx is free software; you can redistribute it and/or modify
   7 it under the terms of the GNU General Public License as published by
   8 the Free Software Foundation; either version 2 of the License, or
   9 (at your option) any later version.
  10
  11 3dengfx is distributed in the hope that it will be useful,
  12 but WITHOUT ANY WARRANTY; without even the implied warranty of
  13 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  14 GNU General Public License for more details.
  15
  16 You should have received a copy of the GNU General Public License
  17 along with 3dengfx; if not, write to the Free Software
  18 Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
  19 */
  20
  21 /* geometry generation
  22  *
  23  * author: John Tsiombikas 2004
  24  * modified:
  25  *              Mihalis Georgoulopoulos 2004, 2005
  26  *              John Tsiombikas 2005
  27  */
  28
  29 #include <time.h>
  30 #include <float.h>
  31 #include "3dengfx_config.h"
  32 #include "gfx/curves.hpp"
  33 #include "ggen.hpp"
  34
  35 #define GGEN_SOURCE
  36 #include "teapot.h"
  37
  38 /* create_cube - (JT)
  39  * creates a subdivided cube
  40  */
  41 void create_cube(TriMesh *mesh, scalar_t size, int subdiv) {
  42         TriMesh tmp;
  43         Matrix4x4 mat;
  44
  45         Vector3 face_vec[] = {
  46                 Vector3(0, 0, -1),
  47                 Vector3(1, 0, 0),
  48                 Vector3(0, 0, 1),
  49                 Vector3(-1, 0, 0),
  50                 Vector3(0, 1, 0),
  51                 Vector3(0, -1, 0)
  52         };
  53
  54         for(int i=0; i<6; i++) {
  55                 create_plane(i ? &tmp : mesh, face_vec[i], Vector2(size, size), subdiv);
  56                 mat.set_translation(face_vec[i] * (size / 2.0));
  57
  58                 if(i) {
  59                         tmp.apply_xform(mat);
  60                         join_tri_mesh(mesh, mesh, &tmp);
  61                 } else {
  62                         mesh->apply_xform(mat);
  63                 }
  64         }
  65 }
  66
  67 /* CreatePlane - (JT)
  68  * creates a planar mesh of arbitrary subdivision
  69  */
  70 void create_plane(TriMesh *mesh, const Vector3 &normal, const Vector2 &size, int subdiv) {
  71         unsigned long vcount = (subdiv + 2) * (subdiv + 2);
  72         unsigned long tcount = (subdiv + 1) * (subdiv + 1) * 2;
  73         int quad_num = tcount / 2;
  74         int quads_row = subdiv + 1;             // quads per row
  75         int vrow = subdiv + 2;          // vertices per row
  76         int vcol = vrow;
  77
  78     Vertex *varray = new Vertex[vcount];
  79         Triangle *tarray = new Triangle[tcount];
  80
  81         for(int j=0; j<vcol; j++) {
  82                 for(int i=0; i<vrow; i++) {
  83
  84                         scalar_t u = (scalar_t)i/(scalar_t)quads_row;
  85                         scalar_t v = (scalar_t)j/(scalar_t)quads_row;
  86                         varray[j * vrow + i] = Vertex(Vector3(u - 0.5f, 1.0f - v - 0.5f, 0), u, 1.0 - v, Color(1.0f));
  87                         varray[j * vrow + i].pos.x *= size.x;
  88                         varray[j * vrow + i].pos.y *= size.y;
  89                 }
  90         }
  91
  92         // first seperate the quads and then triangulate
  93         Quad *quads = new Quad[quad_num];
  94
  95         for(int i=0; i<quad_num; i++) {
  96                 quads[i].vertices[0] = i + i / quads_row;
  97                 quads[i].vertices[1] = quads[i].vertices[0] + 1;
  98                 quads[i].vertices[2] = quads[i].vertices[0] + vrow;
  99                 quads[i].vertices[3] = quads[i].vertices[1] + vrow;
 100         }
 101
 102         for(int i=0; i<quad_num; i++) {
 103                 tarray[i * 2] = Triangle(quads[i].vertices[0], quads[i].vertices[1], quads[i].vertices[3]);
 104                 tarray[i * 2 + 1] = Triangle(quads[i].vertices[0], quads[i].vertices[3], quads[i].vertices[2]);
 105         }
 106
 107         for(unsigned long i=0; i<vcount; i++) {
 108                 varray[i].normal = Vector3(0, 0, -1);
 109         }
 110
 111         // reorient the plane to match the specification
 112         Basis b;
 113         Vector3 n = normal.normalized();
 114         b.k = -n;
 115
 116         if(fabs(dot_product(n, Vector3(1, 0, 0))) < 1.0 - small_number) {
 117                 b.i = Vector3(1, 0, 0);
 118                 b.j = cross_product(b.k, b.i);
 119                 b.i = cross_product(b.j, b.k);
 120         } else {
 121                 b.j = Vector3(0, 1, 0);
 122                 b.i = cross_product(b.j, b.k);
 123                 b.j = cross_product(b.k, b.i);
 124         }
 125
 126         Matrix3x3 rot = b.create_rotation_matrix();
 127
 128         for(unsigned long i=0; i<vcount; i++) {
 129                 varray[i].pos.transform(rot);
 130                 varray[i].normal.transform(rot);
 131         }
 132
 133         mesh->set_data(varray, vcount, tarray, tcount);
 134
 135         delete [] quads;
 136         delete [] varray;
 137         delete [] tarray;
 138 }
 139
 140 /* CreateCylinder - (JT)
 141  * creates a cylinder by extruding a circle along the y axis, with optional
 142  * caps at each end.
 143  */
 144 void create_cylinder(TriMesh *mesh, scalar_t rad, scalar_t len, bool caps, int udiv, int vdiv) {
 145         if(udiv < 3) udiv = 3;
 146         Vector3 *circle = new Vector3[udiv + 1];
 147
 148         // generate the circle
 149         Vector3 cgen(0.0, -len / 2.0, rad);
 150
 151         for(int i=0; i<udiv; i++) {
 152                 Matrix3x3 mat;
 153                 mat.set_rotation(Vector3(0.0, two_pi * (scalar_t)i / (scalar_t)udiv, 0.0));
 154                 circle[i] = cgen.transformed(mat);
 155         }
 156         circle[udiv] = cgen;
 157
 158         // extrude along y
 159         int slices = vdiv + 2;
 160         int vcount = (udiv + 1) * slices + (caps ? 2 * (udiv + 2) : 0);
 161         Vertex *verts = new Vertex[vcount];
 162         scalar_t yoffs = len / (vdiv + 1);
 163
 164         Vertex *vptr = verts;
 165         for(int i=0; i<slices; i++) {
 166                 for(int j=0; j<=udiv; j++) {
 167                         vptr->pos = circle[j];
 168                         vptr->pos.y += yoffs * i;
 169                         vptr->normal = circle[j];
 170                         vptr->normal.y = 0.0;
 171                         vptr->normal /= rad;
 172                         scalar_t u = (scalar_t)j / (scalar_t)udiv;
 173                         scalar_t v = (scalar_t)i / (scalar_t)(vdiv + 1);
 174                         vptr->tex[0] = vptr->tex[1] = TexCoord(u, v);
 175                         vptr++;
 176                 }
 177         }
 178
 179         if(caps) {
 180                 Vertex *cap1 = vptr;
 181                 Vertex *cap2 = vptr + udiv + 2;
 182
 183                 for(int i=0; i<=udiv; i++) {
 184                         circle[i].y = -len/2.0;
 185                         cap1->pos = circle[i];
 186                         cap1->normal = Vector3(0.0, -1.0, 0.0);
 187                         //cap1->tex[0] = vptr->tex[1] = TexCoord((scalar_t)i / (scalar_t)udiv, 0.0);
 188                         cap1->tex[0].u = cap1->tex[1].u = (cap1->pos.x / rad) * 0.5 + 0.5;
 189                         cap1->tex[0].v = cap1->tex[1].v = (cap1->pos.z / rad) * 0.5 + 0.5;
 190
 191                         *cap2 = *cap1;
 192                         cap2->pos.y = len/2.0;
 193                         cap2->normal.y *= -1.0;
 194                         //cap2->tex[0].v = cap2->tex[1].v = 1.0;
 195
 196                         cap1++;
 197                         cap2++;
 198                 }
 199
 200                 *cap1 = Vertex(Vector3(0.0, -len/2.0, 0.0), 0.5, 0.5);
 201                 cap1->normal = Vector3(0.0, -1.0, 0.0);
 202
 203                 *cap2 = *cap1;
 204                 cap2->pos.y = len/2.0;
 205                 cap2->normal.y *= -1.0;
 206                 //cap2->tex[0].v = cap2->tex[1].v = 1.0;
 207         }
 208         delete [] circle;
 209
 210         // triangulate
 211         int tcount = 2 * udiv * (slices - 1) + (caps ? udiv * 2 : 0);
 212         Triangle *triangles = new Triangle[tcount];
 213         Triangle *tptr = triangles;
 214         int v = 0;
 215         for(int i=0; i<slices-1; i++) {
 216                 for(int j=0; j<udiv; j++) {
 217                         *tptr++ = Triangle(v + udiv + 1, v + 1, v + udiv + 2);
 218                         *tptr++ = Triangle(v + udiv + 1, v, v + 1);
 219                         v++;
 220                 }
 221                 v++;
 222         }
 223
 224         if(caps) {
 225                 v += udiv + 1;
 226                 int v2 = v + udiv + 2;
 227                 Triangle *tcap2 = tptr + udiv;
 228                 for(int i=0; i<udiv; i++) {
 229                         *tptr++ = Triangle(v + udiv + 1, v + i + 1, v + i);
 230                         *tcap2++ = Triangle(v2 + udiv + 1, v2 + i, v2 + i + 1);
 231                 }
 232         }
 233
 234         mesh->set_data(verts, vcount, triangles, tcount);
 235
 236         delete [] verts;
 237         delete [] triangles;
 238 }
 239
 240 /* CreateSphere - (MG)
 241  * creates a sphere as a solid of revolution
 242  */
 243 void create_sphere(TriMesh *mesh, scalar_t radius, int subdiv) {
 244         // Solid of revolution. A slice of pi rads is rotated
 245         // for 2pi rads. Subdiv in this revolution should be
 246         // double than subdiv of the slice, because the angle
 247         // is double.
 248
 249         unsigned long edges_pi  = 2 * subdiv;
 250         unsigned long edges_2pi = 4 * subdiv;
 251
 252         unsigned long vcount_pi  = edges_pi  + 1;
 253         unsigned long vcount_2pi = edges_2pi + 1;
 254
 255         unsigned long vcount = vcount_pi * vcount_2pi;
 256
 257         unsigned long quad_count = edges_pi * edges_2pi;
 258         unsigned long tcount = quad_count * 2;
 259
 260         Vertex *varray = new Vertex[vcount];
 261         Triangle *tarray = new Triangle[tcount];
 262
 263         for(unsigned long j = 0; j < vcount_pi; j++) {
 264                 for(unsigned long i = 0; i < vcount_2pi; i++) {
 265
 266                         Vector3 up_vec(0,1,0);
 267
 268                         scalar_t rotx,roty;
 269                         rotx = -(pi * j) / (vcount_pi - 1);
 270                         roty = -(two_pi * i) / (vcount_2pi - 1);
 271
 272                         Matrix4x4 rot_mat;
 273                         rot_mat.set_rotation(Vector3(rotx, 0, 0));
 274                         up_vec.transform(rot_mat);
 275                         rot_mat.set_rotation(Vector3(0, roty, 0));
 276                         up_vec.transform(rot_mat);
 277
 278                         scalar_t u = (scalar_t)i / (scalar_t)(vcount_2pi - 1);
 279                         scalar_t v = 1.0 - (scalar_t)j / (scalar_t)(vcount_pi - 1);
 280                         varray[j * vcount_2pi + i] = Vertex(up_vec * radius, u, v, Color(1.0f));
 281                         varray[j * vcount_2pi + i].normal = up_vec;
 282                 }
 283         }
 284
 285         // first seperate the quads and then triangulate
 286         Quad *quads = new Quad[quad_count];
 287
 288         for(unsigned long i=0; i<quad_count; i++) {
 289
 290                 unsigned long hor_edge,vert_edge;
 291                 hor_edge  = i % edges_2pi;
 292                 vert_edge = i / edges_2pi;
 293
 294                 quads[i].vertices[0] = hor_edge + vert_edge * vcount_2pi;
 295                 quads[i].vertices[1] = quads[i].vertices[0] + 1;
 296                 quads[i].vertices[2] = quads[i].vertices[0] + vcount_2pi;
 297                 quads[i].vertices[3] = quads[i].vertices[1] + vcount_2pi;
 298         }
 299
 300         for(unsigned long i=0; i<quad_count; i++) {
 301                 tarray[i * 2] = Triangle(quads[i].vertices[0], quads[i].vertices[1], quads[i].vertices[3]);
 302                 tarray[i * 2 + 1] = Triangle(quads[i].vertices[0], quads[i].vertices[3], quads[i].vertices[2]);
 303         }
 304
 305         mesh->set_data(varray, vcount, tarray, tcount);
 306
 307         delete [] quads;
 308         delete [] varray;
 309         delete [] tarray;
 310 }
 311
 312 /* CreateTorus - (MG)
 313  * Creates a toroid mesh
 314  */
 315 void create_torus(TriMesh *mesh, scalar_t circle_rad, scalar_t revolv_rad, int subdiv) {
 316         unsigned long edges_2pi  = 4 * subdiv;
 317         unsigned long vcount_2pi = edges_2pi + 1;
 318
 319         unsigned long vcount = vcount_2pi * vcount_2pi;
 320         unsigned long qcount = edges_2pi * edges_2pi;
 321         unsigned long tcount = qcount * 2;
 322
 323         // alloc mamory
 324         Vertex   *varray = new Vertex[vcount];
 325         Quad     *qarray = new Quad[qcount];
 326         Triangle *tarray = new Triangle[tcount];
 327         Vertex   *circle = new Vertex[vcount_2pi];
 328
 329         // first create a circle
 330         // rotation of this circle will produce the torus
 331         for (unsigned long i = 0; i < vcount_2pi; i++)
 332         {
 333                 scalar_t t = (scalar_t)i / (scalar_t)(vcount_2pi - 1);
 334
 335                 Vector3 up_vec = Vector3(0, 1, 0);
 336                 Matrix4x4 rot_mat;
 337                 rot_mat.set_rotation(Vector3(0, 0, two_pi * t));
 338                 up_vec.transform(rot_mat);
 339
 340                 Vector3 pos_vec = up_vec * circle_rad;
 341                 pos_vec += Vector3(revolv_rad, 0, 0);
 342
 343                 circle[i] = Vertex(pos_vec, 0, t, Color(1.0f));
 344                 circle[i].normal = up_vec;
 345         }
 346
 347         // vertex loop
 348         for (unsigned long j = 0; j < vcount_2pi; j++)
 349         {
 350                 for (unsigned long i = 0; i < vcount_2pi; i++)
 351                 {
 352                         scalar_t t = (scalar_t)i / (scalar_t)(vcount_2pi - 1);
 353
 354                         Vector3 pos,nor;
 355                         pos = circle[j].pos;
 356                         nor = circle[j].normal;
 357
 358                         Matrix4x4 rot_mat;
 359                         rot_mat.set_rotation(Vector3(0, two_pi * t, 0));
 360                         pos.transform(rot_mat);
 361                         nor.transform(rot_mat);
 362
 363                         unsigned long index = i + vcount_2pi * j;
 364
 365                         varray[index] = Vertex(pos, t, 1.0 - circle[j].tex[0].v, Color(1.0f));
 366                         varray[index].normal = nor;
 367                 }
 368         } // End vertex loop
 369
 370
 371         for(unsigned long i=0; i<qcount; i++) {
 372
 373                 unsigned long hor_edge,vert_edge;
 374                 hor_edge  = i % edges_2pi;
 375                 vert_edge = i / edges_2pi;
 376
 377                 qarray[i].vertices[0] = hor_edge + vert_edge * vcount_2pi;
 378                 qarray[i].vertices[1] = qarray[i].vertices[0] + 1;
 379                 qarray[i].vertices[2] = qarray[i].vertices[0] + vcount_2pi;
 380                 qarray[i].vertices[3] = qarray[i].vertices[1] + vcount_2pi;
 381         }
 382
 383         for(unsigned long i=0; i<qcount; i++) {
 384                 tarray[i * 2] = Triangle(qarray[i].vertices[0], qarray[i].vertices[1], qarray[i].vertices[3]);
 385                 tarray[i * 2 + 1] = Triangle(qarray[i].vertices[0], qarray[i].vertices[3], qarray[i].vertices[2]);
 386         }
 387
 388         mesh->set_data(varray, vcount, tarray, tcount);
 389
 390         // cleanup
 391         delete [] varray;
 392         delete [] qarray;
 393         delete [] tarray;
 394         delete [] circle;
 395 }
 396
 397
 398 /* create_revolution - (JT)
 399  * Creates a surface of revolution by rotating a curve around the Y axis.
 400  */
 401 void create_revolution(TriMesh *mesh, const Curve &curve, int udiv, int vdiv) {
 402         if(udiv < 3) udiv = 3;
 403         if(vdiv < 1) vdiv = 1;
 404
 405         int slices = udiv;
 406         int stacks = vdiv + 1;
 407
 408         // create the slice that will be revolved to create the mesh
 409         Vector3 *slice = new Vector3[stacks];
 410         Vector3 *slice_normal = new Vector3[stacks];
 411
 412         for(int i=0; i<stacks; i++) {
 413                 scalar_t t = (scalar_t)i / (scalar_t)vdiv;
 414                 slice[i] = curve(t);
 415
 416                 // calculate normal
 417                 Vector3 bitangent = (curve(t + 0.0001) - curve(t - 0.0001)).normalized();
 418                 Vector3 tp1 = slice[i].rotated(Vector3(0.0, DEG_TO_RAD(3), 0.0));
 419                 Vector3 tp2 = slice[i].rotated(Vector3(0.0, -DEG_TO_RAD(3), 0.0));
 420                 Vector3 tangent = (tp1 - tp2).normalized();
 421
 422                 slice_normal[i] = cross_product(tangent, bitangent);
 423         }
 424
 425         int vcount = stacks * slices;
 426         int quad_count = udiv * vdiv;
 427         int tcount = quad_count * 2;
 428
 429         Vertex *varray = new Vertex[vcount];
 430         Triangle *tarray = new Triangle[tcount];
 431
 432         Vertex *vptr = varray;
 433         Triangle *tptr = tarray;
 434
 435         for(int i=0; i<slices; i++) {
 436                 Matrix4x4 rot;
 437                 rot.set_rotation(Vector3(0.0, two_pi * (scalar_t)i / (scalar_t)udiv, 0.0));
 438
 439                 for(int j=0; j<stacks; j++) {
 440                         // create the vertex
 441                         vptr->pos = slice[j].transformed(rot);
 442                         vptr->tex[0].u = vptr->tex[1].u = (scalar_t)i / (scalar_t)udiv;
 443                         vptr->tex[0].v = vptr->tex[1].v = (scalar_t)j / (scalar_t)vdiv;
 444                         vptr->normal = slice_normal[j].transformed(rot);
 445                         vptr++;
 446
 447                         if(j < vdiv) {
 448                                 // create the quad
 449                                 Quad q;
 450                                 q.vertices[0] = j + (i % slices) * stacks;
 451                                 q.vertices[1] = j + ((i + 1) % slices) * stacks;
 452                                 q.vertices[2] = (j + 1) + ((i + 1) % slices) * stacks;
 453                                 q.vertices[3] = (j + 1) + (i % slices) * stacks;
 454
 455                                 // triangulate
 456                                 tptr->vertices[0] = q.vertices[0];
 457                                 tptr->vertices[1] = q.vertices[1];
 458                                 tptr->vertices[2] = q.vertices[2];
 459                                 tptr++;
 460                                 tptr->vertices[0] = q.vertices[0];
 461                                 tptr->vertices[1] = q.vertices[2];
 462                                 tptr->vertices[2] = q.vertices[3];
 463                                 tptr++;
 464                         }
 465                 }
 466         }
 467
 468         mesh->set_data(varray, vcount, tarray, tcount);
 469         delete [] varray;
 470         delete [] tarray;
 471         delete [] slice;
 472         delete [] slice_normal;
 473 }
 474
 475 /* create_revolution - helper function - (JT)
 476  * accepts an array of data points, fits a spline through them and passes the rest
 477  * to the real create_revolution defined above.
 478  */
 479 void create_revolution(TriMesh *mesh, const Vector3 *data, int count, int udiv, int vdiv) {
 480         CatmullRomSplineCurve spline;
 481         for(int i=0; i<count; i++) {
 482                 spline.add_control_point(data[i]);
 483         }
 484         spline.set_arc_parametrization(true);
 485
 486         create_revolution(mesh, spline, udiv, vdiv);
 487 }
 488
 489 /* create_extrusion - (JT)
 490  * Takes a shape and extrudes it along a path.
 491  */
 492 void create_extrusion(TriMesh *mesh, const Curve &shape, const Curve &path, int udiv, int vdiv, scalar_t start_scale, scalar_t end_scale) {
 493         if(udiv < 3) udiv = 3;
 494         Vector3 *shape_pt = new Vector3[udiv + 1];
 495
 496         for(int i=0; i<udiv; i++) {
 497                 shape_pt[i] = shape((scalar_t)i / (scalar_t)udiv);
 498         }
 499         shape_pt[udiv] = shape_pt[0];
 500
 501         // extrude along the spline
 502         int slices = vdiv + 2;
 503         int vcount = (udiv + 1) * slices;
 504         Vertex *verts = new Vertex[vcount];
 505         scalar_t dt = 1.0 / (vdiv + 1);
 506
 507         Vertex *vptr = verts;
 508         for(int i=0; i<slices; i++) {
 509                 for(int j=0; j<=udiv; j++) {
 510                         // XXX FIX THIS
 511                         vptr->pos = shape_pt[j];
 512                         vptr->pos.y += dt * i;
 513                         scalar_t u = (scalar_t)j / (scalar_t)udiv;
 514                         scalar_t v = (scalar_t)i / (scalar_t)(vdiv + 1);
 515                         vptr->tex[0] = vptr->tex[1] = TexCoord(u, v);
 516                         vptr++;
 517                 }
 518         }
 519
 520         delete [] shape_pt;
 521
 522         // triangulate
 523         int tcount = 2 * udiv * (slices - 1);
 524         Triangle *triangles = new Triangle[tcount];
 525         Triangle *tptr = triangles;
 526         int v = 0;
 527         for(int i=0; i<slices-1; i++) {
 528                 for(int j=0; j<udiv; j++) {
 529                         *tptr++ = Triangle(v + udiv + 1, v + 1, v + udiv + 2);
 530                         *tptr++ = Triangle(v + udiv + 1, v, v + 1);
 531                         v++;
 532                 }
 533                 v++;
 534         }
 535
 536         mesh->set_data(verts, vcount, triangles, tcount);
 537         mesh->calculate_normals();
 538
 539         delete [] verts;
 540         delete [] triangles;
 541 }
 542
 543 /* CreateBezierPatch - (MG)
 544  * overloaded function that gets a vector3 array
 545  * and makes a single Bezier patch
 546  */
 547 void create_bezier_patch(TriMesh *mesh, const Vector3 *cp, int subdiv)
 548 {
 549
 550         // make 8 BezierSpline's
 551         BezierSpline u[4], v[4];
 552
 553         u[0].add_control_point(cp[0]);
 554         u[0].add_control_point(cp[1]);
 555         u[0].add_control_point(cp[2]);
 556         u[0].add_control_point(cp[3]);
 557
 558         u[1].add_control_point(cp[4]);
 559         u[1].add_control_point(cp[5]);
 560         u[1].add_control_point(cp[6]);
 561         u[1].add_control_point(cp[7]);
 562
 563         u[2].add_control_point(cp[8]);
 564         u[2].add_control_point(cp[9]);
 565         u[2].add_control_point(cp[10]);
 566         u[2].add_control_point(cp[11]);
 567
 568         u[3].add_control_point(cp[12]);
 569         u[3].add_control_point(cp[13]);
 570         u[3].add_control_point(cp[14]);
 571         u[3].add_control_point(cp[15]);
 572
 573         v[0].add_control_point(cp[0]);
 574         v[0].add_control_point(cp[4]);
 575         v[0].add_control_point(cp[8]);
 576         v[0].add_control_point(cp[12]);
 577
 578         v[1].add_control_point(cp[1]);
 579         v[1].add_control_point(cp[5]);
 580         v[1].add_control_point(cp[9]);
 581         v[1].add_control_point(cp[13]);
 582
 583         v[2].add_control_point(cp[2]);
 584         v[2].add_control_point(cp[6]);
 585         v[2].add_control_point(cp[10]);
 586         v[2].add_control_point(cp[14]);
 587
 588         v[3].add_control_point(cp[3]);
 589         v[3].add_control_point(cp[7]);
 590         v[3].add_control_point(cp[11]);
 591         v[3].add_control_point(cp[15]);
 592
 593         unsigned long edges = subdiv * 2;
 594         unsigned long vrow = edges + 1;
 595         unsigned long qcount = edges * edges;
 596         unsigned long vcount = vrow * vrow;
 597         unsigned long tcount = qcount * 2;
 598
 599         // allocate memory
 600         Vertex *varray = new Vertex[vcount];
 601         Quad *qarray = new Quad[qcount];
 602         Triangle *tarray = new Triangle[tcount];
 603
 604         // Vertex loop
 605         for (unsigned long j=0; j<vrow; j++)
 606         {
 607                 scalar_t tv = (scalar_t)j / (scalar_t)(vrow - 1);
 608                 BezierSpline uc;
 609                 uc.add_control_point(v[0](tv));
 610                 uc.add_control_point(v[1](tv));
 611                 uc.add_control_point(v[2](tv));
 612                 uc.add_control_point(v[3](tv));
 613
 614                 for (unsigned long i=0; i<vrow; i++)
 615                 {
 616                         scalar_t tu = (scalar_t)i / (scalar_t)(vrow - 1);
 617                         BezierSpline vc;
 618                         vc.add_control_point(u[0](tu));
 619                         vc.add_control_point(u[1](tu));
 620                         vc.add_control_point(u[2](tu));
 621                         vc.add_control_point(u[3](tu));
 622
 623                         // get the position
 624                         Vector3 pos = uc(tu);
 625
 626                         // get normal
 627                         Vector3 tan_u,tan_v;
 628                         tan_u = uc.get_tangent(tu);
 629                         tan_v = vc.get_tangent(tv);
 630                         Vector3 normal;
 631                         normal = cross_product(tan_u, tan_v);
 632                         normal.normalize();
 633
 634                         // store vertex
 635                         varray[i + j * vrow] = Vertex(pos, tu, 1.0 - tv, Color(1.0f));
 636                         varray[i + j * vrow].normal = normal;
 637                 }
 638         } // end vertex loop
 639
 640
 641         // first seperate the quads and then triangulate
 642         for(unsigned long i=0; i<qcount; i++) {
 643                 qarray[i].vertices[0] = i + i / edges;
 644                 qarray[i].vertices[1] = qarray[i].vertices[0] + 1;
 645                 qarray[i].vertices[2] = qarray[i].vertices[0] + vrow;
 646                 qarray[i].vertices[3] = qarray[i].vertices[1] + vrow;
 647         }
 648
 649         for(unsigned long i=0; i<qcount; i++) {
 650                 tarray[i * 2] = Triangle(qarray[i].vertices[0], qarray[i].vertices[1], qarray[i].vertices[3]);
 651                 tarray[i * 2 + 1] = Triangle(qarray[i].vertices[0], qarray[i].vertices[3], qarray[i].vertices[2]);
 652         }
 653
 654         mesh->set_data(varray, vcount, tarray, tcount);
 655
 656         // cleanup
 657         delete [] varray;
 658         delete [] qarray;
 659         delete [] tarray;
 660 }
 661
 662 /* CreateBezierPatch - (MG)
 663  * creates a bezier patch (!)
 664  * if the control curves contain more than one
 665  * segments, multiple patches will be included
 666  * in the output TriMesh
 667  */
 668 void create_bezier_patch(TriMesh *mesh, const BezierSpline &u0, const BezierSpline &u1, const BezierSpline &u2, const BezierSpline &u3, int subdiv)
 669 {
 670         // get minimum number of segments
 671         unsigned long min_seg , tmp;
 672         min_seg = u0.get_segment_count();
 673         tmp = u1.get_segment_count(); if (min_seg > tmp) min_seg = tmp;
 674         tmp = u2.get_segment_count(); if (min_seg > tmp) min_seg = tmp;
 675         tmp = u3.get_segment_count(); if (min_seg > tmp) min_seg = tmp;
 676
 677         TriMesh tmp_mesh;
 678         Vector3 *cp = new Vector3[16];
 679         for (unsigned long i=0; i<min_seg; i++)
 680         {
 681                 // fill control point array
 682                 cp[0] = u0.get_control_point(4 * i + 0);
 683                 cp[1] = u0.get_control_point(4 * i + 1);
 684                 cp[2] = u0.get_control_point(4 * i + 2);
 685                 cp[3] = u0.get_control_point(4 * i + 3);
 686
 687                 cp[4] = u1.get_control_point(4 * i + 0);
 688                 cp[5] = u1.get_control_point(4 * i + 1);
 689                 cp[6] = u1.get_control_point(4 * i + 2);
 690                 cp[7] = u1.get_control_point(4 * i + 3);
 691
 692                 cp[8] = u2.get_control_point(4 * i + 0);
 693                 cp[9] = u2.get_control_point(4 * i + 1);
 694                 cp[10] = u2.get_control_point(4 * i + 2);
 695                 cp[11] = u2.get_control_point(4 * i + 3);
 696
 697                 cp[12] = u3.get_control_point(4 * i + 0);
 698                 cp[13] = u3.get_control_point(4 * i + 1);
 699                 cp[14] = u3.get_control_point(4 * i + 2);
 700                 cp[15] = u3.get_control_point(4 * i + 3);
 701
 702                 // Make a single patch and put all patches together
 703                 create_bezier_patch(&tmp_mesh, cp, subdiv);
 704
 705                 join_tri_mesh(mesh, mesh, &tmp_mesh);
 706         }
 707
 708         // cleanup
 709         delete [] cp;
 710 }
 711
 712 /* CreateBezierMesh - (MG)
 713  * tesselates a whole mesh of bezier patches.
 714  * usefull when some patches share vertices
 715  * TODO : Make a bezier patch class , like Triangle
 716  * and Quad. Store indices
 717  */
 718 void create_bezier_mesh(TriMesh *mesh, const Vector3 *cp, unsigned int *patches, int patch_count, int subdiv)
 719 {
 720         TriMesh tmp_mesh;
 721         Vector3 control_pts[16];
 722         for(int i=0; i<patch_count; i++)
 723         {
 724                 for(int j=0; j<16; j++)
 725                 {
 726                         control_pts[j] = cp[ patches[16 * i + j] ];
 727                 }
 728
 729                 create_bezier_patch(&tmp_mesh, control_pts, subdiv);
 730
 731                 join_tri_mesh(mesh, mesh, &tmp_mesh);
 732         }
 733 }
 734
 735 /* CreateTeapot - (MG)
 736  * Creates a teapot TriMesh, using the original
 737  * data file from Newell
 738  */
 739 void create_teapot(TriMesh *mesh, scalar_t size, int subdiv)
 740 {
 741         unsigned int *patches = new unsigned int[teapot_num_patches * 16];
 742
 743         // fix indices to start from zero
 744         for(int i=0; i<teapot_num_patches * 16; i++)
 745         {
 746                 patches[i] = teapot_patches[i] - 1;
 747         }
 748
 749
 750         // rearrange patches to clockwise order
 751         for(int p=0; p<teapot_num_patches; p++)
 752         {
 753                 unsigned int new_order[16];
 754                 for(int j=0; j<4; j++)
 755                 {
 756                         for(int i=0; i<4; i++)
 757                         {
 758                                 new_order[j * 4 + (3 - i)] = patches[p * 16 + j * 4 + i];
 759                         }
 760                 }
 761
 762                 for(int k=0; k<16; k++)
 763                 {
 764                         patches[16 * p + k] = new_order[k];
 765                 }
 766         }
 767
 768         // rearrange vertices to correct axes
 769         Vector3 *vertices = new Vector3[teapot_num_vertices];
 770         for(int i = 0; i < teapot_num_vertices; i++)
 771         {
 772                 vertices[i].x = teapot_vertices[i * 3 + 0] * size;
 773                 vertices[i].z = teapot_vertices[i * 3 + 1] * size;
 774                 vertices[i].y = teapot_vertices[i * 3 + 2] * size;
 775         }
 776
 777         create_bezier_mesh(mesh, vertices, patches, teapot_num_patches, subdiv);
 778
 779         mesh->calculate_normals();
 780
 781         // cleanup
 782         delete [] patches;
 783         delete [] vertices;
 784 }
 785
 786
 787 // fractal stuff ...
 788
 789 /* create_landscape (JT)
 790  * Creates a fractal landscape... (TODO: add algorithm description or something)
 791  */
 792 void create_landscape(TriMesh *mesh, const Vector2 &size, int mesh_detail, scalar_t max_height, int iter, scalar_t roughness, int seed) {
 793         create_plane(mesh, Vector3(0, 1, 0), size, mesh_detail);
 794         roughness *= 0.25;
 795
 796         if(seed == GGEN_RANDOM_SEED) {
 797                 srand(time(0));
 798         } else if(seed != GGEN_NO_RESEED) {
 799                 srand(seed);
 800         }
 801
 802         scalar_t offs = max_height / (scalar_t)iter;
 803
 804         unsigned long vcount = mesh->get_vertex_array()->get_count();
 805         Vertex *varray = mesh->get_mod_vertex_array()->get_mod_data();
 806
 807         for(int i=0; i<iter; i++) {
 808                 // pick a partitioning line (2d)
 809                 Vector2 pt1(frand(size.x) - size.x / 2.0, frand(size.y) - size.y / 2.0);
 810                 Vector2 pt2(frand(size.x) - size.x / 2.0, frand(size.y) - size.y / 2.0);
 811
 812                 // find its normal
 813                 Vector2 normal(pt2.y - pt1.y, pt1.x - pt2.x);
 814
 815                 // classify all points wrt. this line and raise them accordingly.
 816                 for(unsigned long j=0; j<vcount; j++) {
 817                         Vector3 *vpos = &varray[j].pos;
 818                         Vector2 vpos2d(vpos->x, vpos->z);
 819
 820                         scalar_t dist = dist_line(pt1, pt2, vpos2d);
 821
 822                         /* this was considered but it was producing extremely smooth
 823                          * results, which looked unnatural.
 824                          */
 825                         /*
 826                         if(dot_product(normal, vpos2d - pt1) < 0.0) {
 827                                 dist = -dist;
 828                         }
 829                         scalar_t sigmoid = (tanh(dist * size.x * size.y * roughness) + 1.0) * 0.5;
 830                         vpos->y += offs * sigmoid;
 831                         */
 832
 833                         if(dot_product(normal, vpos2d - pt1) > 0.0) {
 834                                 scalar_t sigmoid = tanh(dist * size.x * size.y * roughness);
 835                                 vpos->y += offs * sigmoid;
 836                         }
 837                 }
 838         }
 839
 840         // normalize the landscape in the range [0, max_height)
 841         scalar_t hmin = FLT_MAX, hmax = 0.0;
 842         Vertex *vptr = varray;
 843
 844         for(unsigned long i=0; i<vcount; i++) {
 845                 if(vptr->pos.y > hmax) hmax = vptr->pos.y;
 846                 if(vptr->pos.y < hmin) hmin = vptr->pos.y;
 847                 vptr++;
 848         }
 849
 850         vptr = varray;
 851         for(unsigned long i=0; i<vcount; i++) {
 852                 vptr->pos.y = max_height * (vptr->pos.y - hmin) / (hmax - hmin);
 853                 vptr++;
 854         }
 855
 856         mesh->calculate_normals();
 857 }