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

   1 /*
   2 This file is part of the 3dengfx, realtime visualization system.
   3
   4 Copyright (c) 2004, 2005 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 /* Scalar fields and polygonization
  22  *
  23  * Author: Mihalis Georgoulopoulos 2005
  24  */
  25
  26 #define SCFIELD_SOURCE
  27 #include "mcube_tables.h"
  28 #include "scfield.hpp"
  29 #include "3dengfx/3denginefx.hpp"
  30
  31 // don't change this
  32 #define EDGE_NOT_ASSOCIATED             0xFFFFFFFF
  33
  34 /* -----------------
  35  * private functions
  36  * -----------------
  37  */
  38
  39 /*
  40  * AddVertex
  41  * adds a vertex and returns its index
  42  */
  43 unsigned int ScalarField::add_vertex(const Vertex &vert)
  44 {
  45         verts.push_back(vert);
  46         return verts.size() - 1;
  47 }
  48
  49 /*
  50  * Clear
  51  * clears std::vector's that hold mesh data and resets edges table
  52  */
  53 void ScalarField::clear()
  54 {
  55         verts.erase(verts.begin(), verts.end());
  56         tris.erase(tris.begin(), tris.end());
  57
  58         unsigned int num_bytes = dimensions * dimensions * dimensions * sizeof(unsigned int);
  59
  60         memset(edges_x, 0xFF, num_bytes);
  61         memset(edges_y, 0xFF, num_bytes);
  62         memset(edges_z, 0xFF, num_bytes);
  63 }
  64
  65 /*
  66  * EvaluateAll
  67  * Evaluates all values with the external Evaluate function (if specified)
  68  */
  69 void ScalarField::evaluate_all(scalar_t t)
  70 {
  71         if (!evaluate)
  72         {
  73                 return;
  74         }
  75
  76         for (unsigned int z=0; z<dimensions; z++)
  77         {
  78                 for (unsigned int y=0; y<dimensions; y++)
  79                 {
  80                         for (unsigned int x=0; x<dimensions; x++)
  81                         {
  82                                 set_value(x, y, z, evaluate(get_position(x, y, z), t));
  83                         }
  84                 }
  85         }
  86 }
  87
  88
  89 /*
  90  * ProcesssCell
  91  */
  92 void ScalarField::process_cell(int x, int y, int z, scalar_t isolevel)
  93 {
  94         unsigned char cube_index = 0;
  95         if(get_value(x, y, z, 0) < isolevel) cube_index |= 1;
  96         if(get_value(x, y, z, 1) < isolevel) cube_index |= 2;
  97         if(get_value(x, y, z, 2) < isolevel) cube_index |= 4;
  98         if(get_value(x, y, z, 3) < isolevel) cube_index |= 8;
  99         if(get_value(x, y, z, 4) < isolevel) cube_index |= 16;
 100         if(get_value(x, y, z, 5) < isolevel) cube_index |= 32;
 101         if(get_value(x, y, z, 6) < isolevel) cube_index |= 64;
 102         if(get_value(x, y, z, 7) < isolevel) cube_index |= 128;
 103
 104         int edge_index = cube_edge_flags[cube_index];
 105
 106         scalar_t p , val1, val2;
 107         Vector3 vec1, vec2;
 108
 109         if ( (edge_index & 1) && (get_edge(x, y, z, 0) == EDGE_NOT_ASSOCIATED) )
 110         {
 111                 val1 = get_value(x, y, z, 0);
 112                 val2 = get_value(x, y, z, 1);
 113                 p = (isolevel - val1) / (val2 - val1);
 114
 115                 vec1 = get_position(x, y, z, 0);
 116                 vec2 = get_position(x, y, z, 1);
 117
 118                 set_edge(x, y, z, 0, add_vertex( Vertex( vec1 + p * (vec2 - vec1) ) ) );
 119         }
 120
 121         if ( (edge_index & 2) && (get_edge(x, y, z, 1) == EDGE_NOT_ASSOCIATED) )
 122         {
 123                 val1 = get_value(x, y, z, 1);
 124                 val2 = get_value(x, y, z, 2);
 125                 p = (isolevel - val1) / (val2 - val1);
 126
 127                 vec1 = get_position(x, y, z, 1);
 128                 vec2 = get_position(x, y, z, 2);
 129
 130                 set_edge(x, y, z, 1, add_vertex( Vertex( vec1 + p * (vec2 - vec1) ) ) );
 131         }
 132
 133         if ( (edge_index & 4) && (get_edge(x, y, z, 2) == EDGE_NOT_ASSOCIATED) )
 134         {
 135                 val1 = get_value(x, y, z, 2);
 136                 val2 = get_value(x, y, z, 3);
 137                 p = (isolevel - val1) / (val2 - val1);
 138
 139                 vec1 = get_position(x, y, z, 2);
 140                 vec2 = get_position(x, y, z, 3);
 141
 142                 set_edge(x, y, z, 2, add_vertex( Vertex( vec1 + p * (vec2 - vec1) ) ) );
 143         }
 144
 145         if ( (edge_index & 8) && (get_edge(x, y, z, 3) == EDGE_NOT_ASSOCIATED) )
 146         {
 147                 val1 = get_value(x, y, z, 3);
 148                 val2 = get_value(x, y, z, 0);
 149                 p = (isolevel - val1) / (val2 - val1);
 150
 151                 vec1 = get_position(x, y, z, 0);
 152                 vec2 = get_position(x, y, z, 1);
 153
 154                 set_edge(x, y, z, 3, add_vertex( Vertex( vec1 + p * (vec2 - vec1) ) ) );
 155         }
 156
 157         if ( (edge_index & 16) && (get_edge(x, y, z, 4) == EDGE_NOT_ASSOCIATED) )
 158         {
 159                 val1 = get_value(x, y, z, 4);
 160                 val2 = get_value(x, y, z, 5);
 161                 p = (isolevel - val1) / (val2 - val1);
 162
 163                 vec1 = get_position(x, y, z, 4);
 164                 vec2 = get_position(x, y, z, 5);
 165
 166                 set_edge(x, y, z, 4, add_vertex( Vertex( vec1 + p * (vec2 - vec1) ) ) );
 167         }
 168
 169         if ( (edge_index & 32) && (get_edge(x, y, z, 5) == EDGE_NOT_ASSOCIATED) )
 170         {
 171                 val1 = get_value(x, y, z, 5);
 172                 val2 = get_value(x, y, z, 6);
 173                 p = (isolevel - val1) / (val2 - val1);
 174
 175                 vec1 = get_position(x, y, z, 5);
 176                 vec2 = get_position(x, y, z, 6);
 177
 178                 set_edge(x, y, z, 5, add_vertex( Vertex( vec1 + p * (vec2 - vec1) ) ) );
 179         }
 180
 181         if ( (edge_index & 64) && (get_edge(x, y, z, 6) == EDGE_NOT_ASSOCIATED) )
 182         {
 183                 val1 = get_value(x, y, z, 6);
 184                 val2 = get_value(x, y, z, 7);
 185                 p = (isolevel - val1) / (val2 - val1);
 186
 187                 vec1 = get_position(x, y, z, 6);
 188                 vec2 = get_position(x, y, z, 7);
 189
 190                 set_edge(x, y, z, 6, add_vertex( Vertex( vec1 + p * (vec2 - vec1) ) ) );
 191         }
 192
 193         if ( (edge_index & 128) && (get_edge(x, y, z, 7) == EDGE_NOT_ASSOCIATED) )
 194         {
 195                 val1 = get_value(x, y, z, 7);
 196                 val2 = get_value(x, y, z, 4);
 197                 p = (isolevel - val1) / (val2 - val1);
 198
 199                 vec1 = get_position(x, y, z, 7);
 200                 vec2 = get_position(x, y, z, 4);
 201
 202                 set_edge(x, y, z, 7, add_vertex( Vertex( vec1 + p * (vec2 - vec1) ) ) );
 203         }
 204
 205         if ( (edge_index & 256) && (get_edge(x, y, z, 8) == EDGE_NOT_ASSOCIATED) )
 206         {
 207                 val1 = get_value(x, y, z, 0);
 208                 val2 = get_value(x, y, z, 4);
 209                 p = (isolevel - val1) / (val2 - val1);
 210
 211                 vec1 = get_position(x, y, z, 0);
 212                 vec2 = get_position(x, y, z, 4);
 213
 214                 set_edge(x, y, z, 8, add_vertex( Vertex( vec1 + p * (vec2 - vec1) ) ) );
 215         }
 216
 217         if ( (edge_index & 512) && (get_edge(x, y, z, 9) == EDGE_NOT_ASSOCIATED) )
 218         {
 219                 val1 = get_value(x, y, z, 1);
 220                 val2 = get_value(x, y, z, 5);
 221                 p = (isolevel - val1) / (val2 - val1);
 222
 223                 vec1 = get_position(x, y, z, 1);
 224                 vec2 = get_position(x, y, z, 5);
 225
 226                 set_edge(x, y, z, 9, add_vertex( Vertex( vec1 + p * (vec2 - vec1) ) ) );
 227         }
 228
 229         if ( (edge_index & 1024) && (get_edge(x, y, z, 10) == EDGE_NOT_ASSOCIATED) )
 230         {
 231                 val1 = get_value(x, y, z, 2);
 232                 val2 = get_value(x, y, z, 6);
 233                 p = (isolevel - val1) / (val2 - val1);
 234
 235                 vec1 = get_position(x, y, z, 2);
 236                 vec2 = get_position(x, y, z, 6);
 237
 238                 set_edge(x, y, z, 10, add_vertex( Vertex( vec1 + p * (vec2 - vec1) ) ) );
 239         }
 240
 241         if ( (edge_index & 2048) && (get_edge(x, y, z, 11) == EDGE_NOT_ASSOCIATED) )
 242         {
 243                 val1 = get_value(x, y, z, 3);
 244                 val2 = get_value(x, y, z, 7);
 245                 p = (isolevel - val1) / (val2 - val1);
 246
 247                 vec1 = get_position(x, y, z, 3);
 248                 vec2 = get_position(x, y, z, 7);
 249
 250                 set_edge(x, y, z, 11, add_vertex( Vertex( vec1 + p * (vec2 - vec1) ) ) );
 251         }
 252
 253         // Add triangles
 254         unsigned int p1, p2, p3;
 255
 256         if (tri_table[cube_index][0] != -1)
 257         {
 258                 p1 = get_edge(x, y, z, tri_table[cube_index][0]);
 259                 p2 = get_edge(x, y, z, tri_table[cube_index][1]);
 260                 p3 = get_edge(x, y, z, tri_table[cube_index][2]);
 261                 tris.push_back(Triangle(p2, p1, p3));
 262         }
 263
 264         if (tri_table[cube_index][3] != -1)
 265         {
 266                 p1 = get_edge(x, y, z, tri_table[cube_index][3]);
 267                 p2 = get_edge(x, y, z, tri_table[cube_index][4]);
 268                 p3 = get_edge(x, y, z, tri_table[cube_index][5]);
 269                 tris.push_back(Triangle(p2, p1, p3));
 270         }
 271
 272         if (tri_table[cube_index][6] != -1)
 273         {
 274                 p1 = get_edge(x, y, z, tri_table[cube_index][6]);
 275                 p2 = get_edge(x, y, z, tri_table[cube_index][7]);
 276                 p3 = get_edge(x, y, z, tri_table[cube_index][8]);
 277                 tris.push_back(Triangle(p2, p1, p3));
 278         }
 279
 280         if (tri_table[cube_index][9] != -1)
 281         {
 282                 p1 = get_edge(x, y, z, tri_table[cube_index][9]);
 283                 p2 = get_edge(x, y, z, tri_table[cube_index][10]);
 284                 p3 = get_edge(x, y, z, tri_table[cube_index][11]);
 285                 tris.push_back(Triangle(p2, p1, p3));
 286         }
 287
 288         if (tri_table[cube_index][12] != -1)
 289         {
 290                 p1 = get_edge(x, y, z, tri_table[cube_index][12]);
 291                 p2 = get_edge(x, y, z, tri_table[cube_index][13]);
 292                 p3 = get_edge(x, y, z, tri_table[cube_index][14]);
 293                 tris.push_back(Triangle(p2, p1, p3));
 294         }
 295 }
 296
 297 /*
 298  * GetValueIndex
 299  * returns the index to the values array for the specified coords
 300  */
 301 unsigned int ScalarField::get_value_index(int x, int y, int z)
 302 {
 303         return x + y * dimensions + z * dimensions * dimensions;
 304 }
 305
 306
 307 Vector3 ScalarField::def_eval_normals(const Vector3 &vec, scalar_t t) {
 308         if(!evaluate) return Vector3(0, 0, 0);
 309
 310         Vector3 diff = cell_size * 0.25;
 311
 312         Vector3 grad;
 313         grad.x = evaluate(vec + Vector3(diff.x, 0, 0), t) - evaluate(vec + Vector3(-diff.x, 0, 0), t);
 314         grad.y = evaluate(vec + Vector3(0, diff.y, 0), t) - evaluate(vec + Vector3(0, -diff.y, 0), t);
 315         grad.z = evaluate(vec + Vector3(0, 0, diff.z), t) - evaluate(vec + Vector3(0, 0, -diff.z), t);
 316
 317         return grad.normalized();
 318 }
 319
 320 /* --------------
 321  * public methods
 322  * --------------
 323  */
 324
 325 // constructor
 326 ScalarField::ScalarField()
 327 {
 328         values = 0;
 329         edges_x = edges_y = edges_z = 0;
 330         dimensions = 0;
 331         from = to = cell_size =  Vector3(0, 0, 0);
 332         evaluate = 0;
 333         get_normal = 0;
 334 }
 335
 336 ScalarField::ScalarField(unsigned int dimensions, const Vector3 &from, const Vector3 &to)
 337 {
 338         this->dimensions = dimensions;
 339         this->from = from;
 340         this->to = to;
 341         this->cell_size =  (to - from) / (dimensions - 1);
 342
 343         values = 0;
 344         edges_x = edges_y = edges_z = 0;
 345
 346         evaluate = 0;
 347         get_normal = 0;
 348
 349         set_dimensions(dimensions);
 350 }
 351
 352 ScalarField::~ScalarField()
 353 {
 354         if (values)
 355                 delete [] values;
 356         if (edges_x)
 357                 delete [] edges_x;
 358         if (edges_y)
 359                 delete [] edges_y;
 360         if (edges_z)
 361                 delete [] edges_z;
 362 }
 363
 364 void ScalarField::set_dimensions(unsigned int dimensions)
 365 {
 366         this->dimensions = dimensions;
 367         if (values)
 368                 delete [] values;
 369         if (edges_x)
 370                 delete [] edges_x;
 371         if (edges_y)
 372                 delete [] edges_y;
 373         if (edges_z)
 374                 delete [] edges_z;
 375
 376         values = new scalar_t [dimensions * dimensions * dimensions];
 377
 378         unsigned int edges_per_dim = dimensions * dimensions * dimensions;
 379         edges_x = new unsigned int[edges_per_dim];
 380         edges_y = new unsigned int[edges_per_dim];
 381         edges_z = new unsigned int[edges_per_dim];
 382
 383         clear();
 384 }
 385
 386 // Get / Set
 387 void ScalarField::set_value(int x, int y, int z, scalar_t value)
 388 {
 389         values[get_value_index(x, y, z)] = value;
 390 }
 391
 392 scalar_t ScalarField::get_value(int x, int y, int z)
 393 {
 394         return values[get_value_index(x, y, z)];
 395 }
 396
 397
 398 // get / set relative to cell
 399 void ScalarField::set_value(int cx, int cy, int cz, int vert_index, scalar_t value)
 400 {
 401         if (vert_index == 0) set_value(cx + 0, cy + 0, cz + 1, value);
 402         else if (vert_index == 1) set_value(cx + 1, cy + 0, cz + 1, value);
 403         else if (vert_index == 2) set_value(cx + 1, cy + 0, cz + 0, value);
 404         else if (vert_index == 3) set_value(cx + 0, cy + 0, cz + 0, value);
 405         else if (vert_index == 4) set_value(cx + 0, cy + 1, cz + 1, value);
 406         else if (vert_index == 5) set_value(cx + 1, cy + 1, cz + 1, value);
 407         else if (vert_index == 6) set_value(cx + 1, cy + 1, cz + 0, value);
 408         else if (vert_index == 7) set_value(cx + 0, cy + 1, cz + 0, value);
 409 }
 410
 411
 412 scalar_t ScalarField::get_value(int cx, int cy, int cz, int vert_index)
 413 {
 414         if (vert_index == 0) return get_value(cx + 0, cy + 0, cz + 1);
 415         if (vert_index == 1) return get_value(cx + 1, cy + 0, cz + 1);
 416         if (vert_index == 2) return get_value(cx + 1, cy + 0, cz + 0);
 417         if (vert_index == 3) return get_value(cx + 0, cy + 0, cz + 0);
 418         if (vert_index == 4) return get_value(cx + 0, cy + 1, cz + 1);
 419         if (vert_index == 5) return get_value(cx + 1, cy + 1, cz + 1);
 420         if (vert_index == 6) return get_value(cx + 1, cy + 1, cz + 0);
 421         if (vert_index == 7) return get_value(cx + 0, cy + 1, cz + 0);
 422
 423         return 0;
 424 }
 425
 426 // edges are addressed relatively to a cell (cx, cy, cz)
 427 // and the cell's edge number
 428 void ScalarField::set_edge(int cx, int cy, int cz, int edge, unsigned int index)
 429 {
 430         unsigned int d = dimensions;
 431         unsigned int d2 = dimensions * dimensions;
 432
 433         if (edge == 0)
 434                 edges_x[cx + 0 + (cy + 0) * d + (cz + 1) * d2] =  index;
 435         else if (edge == 1)
 436                 edges_z[cx + 1 + (cy + 0) * d + (cz + 0) * d2] =  index;
 437         else if (edge == 2)
 438                 edges_x[cx + 0 + (cy + 0) * d + (cz + 0) * d2] =  index;
 439         else if (edge == 3)
 440                 edges_z[cx + 0 + (cy + 0) * d + (cz + 0) * d2] =  index;
 441         else if (edge == 4)
 442                 edges_x[cx + 0 + (cy + 1) * d + (cz + 1) * d2] =  index;
 443         else if (edge == 5)
 444                 edges_z[cx + 1 + (cy + 1) * d + (cz + 0) * d2] =  index;
 445         else if (edge == 6)
 446                 edges_x[cx + 0 + (cy + 1) * d + (cz + 0) * d2] =  index;
 447         else if (edge == 7)
 448                 edges_z[cx + 0 + (cy + 1) * d + (cz + 0) * d2] =  index;
 449         else if (edge == 8)
 450                 edges_y[cx + 0 + (cy + 0) * d + (cz + 1) * d2] =  index;
 451         else if (edge == 9)
 452                 edges_y[cx + 1 + (cy + 0) * d + (cz + 1) * d2] =  index;
 453         else if (edge == 10)
 454                 edges_y[cx + 1 + (cy + 0) * d + (cz + 0) * d2] =  index;
 455         else if (edge == 11)
 456                 edges_y[cx + 0 + (cy + 0) * d + (cz + 0) * d2] =  index;
 457 }
 458
 459 unsigned int ScalarField::get_edge(int cx, int cy, int cz, int edge)
 460 {
 461         unsigned int d = dimensions;
 462         unsigned int d2 = dimensions * dimensions;
 463
 464         if (edge == 0)  return edges_x[cx + 0 + (cy + 0) * d + (cz + 1) * d2];
 465         if (edge == 1)  return edges_z[cx + 1 + (cy + 0) * d + (cz + 0) * d2];
 466         if (edge == 2)  return edges_x[cx + 0 + (cy + 0) * d + (cz + 0) * d2];
 467         if (edge == 3)  return edges_z[cx + 0 + (cy + 0) * d + (cz + 0) * d2];
 468         if (edge == 4)  return edges_x[cx + 0 + (cy + 1) * d + (cz + 1) * d2];
 469         if (edge == 5)  return edges_z[cx + 1 + (cy + 1) * d + (cz + 0) * d2];
 470         if (edge == 6)  return edges_x[cx + 0 + (cy + 1) * d + (cz + 0) * d2];
 471         if (edge == 7)  return edges_z[cx + 0 + (cy + 1) * d + (cz + 0) * d2];
 472         if (edge == 8)  return edges_y[cx + 0 + (cy + 0) * d + (cz + 1) * d2];
 473         if (edge == 9)  return edges_y[cx + 1 + (cy + 0) * d + (cz + 1) * d2];
 474         if (edge == 10) return edges_y[cx + 1 + (cy + 0) * d + (cz + 0) * d2];
 475         if (edge == 11) return edges_y[cx + 0 + (cy + 0) * d + (cz + 0) * d2];
 476
 477         return 0;
 478 }
 479
 480 // Position in space
 481 Vector3 ScalarField::get_position(int x, int y, int z)
 482 {
 483         scalar_t vx, vy, vz;
 484         vx = from.x + cell_size.x * x;
 485         vy = from.y + cell_size.y * y;
 486         vz = from.z + cell_size.z * z;
 487
 488         return Vector3(vx, vy, vz);
 489 }
 490
 491 Vector3 ScalarField::get_position(int cx, int cy, int cz, int vert_index)
 492 {
 493         if (vert_index == 0) return get_position(cx + 0, cy + 0, cz + 1);
 494         if (vert_index == 1) return get_position(cx + 1, cy + 0, cz + 1);
 495         if (vert_index == 2) return get_position(cx + 1, cy + 0, cz + 0);
 496         if (vert_index == 3) return get_position(cx + 0, cy + 0, cz + 0);
 497         if (vert_index == 4) return get_position(cx + 0, cy + 1, cz + 1);
 498         if (vert_index == 5) return get_position(cx + 1, cy + 1, cz + 1);
 499         if (vert_index == 6) return get_position(cx + 1, cy + 1, cz + 0);
 500         if (vert_index == 7) return get_position(cx + 0, cy + 1, cz + 0);
 501
 502         return Vector3(0, 0, 0);
 503 }
 504
 505 void ScalarField::set_from_to(const Vector3 &from, const Vector3 &to)
 506 {
 507         this->from = from;
 508         this->to = to;
 509         this->cell_size = (to - from) / (dimensions - 1);
 510 }
 511
 512 void ScalarField::draw_field(bool full)
 513 {
 514         set_lighting(false);
 515
 516         if (full)
 517         {
 518                 glBegin(GL_LINES);
 519                 {
 520                         glColor4f(1.0f, 1.0f, 1.0f, 1.0f);
 521
 522                         for (unsigned int j=0; j<dimensions; j++)
 523                         {
 524                                 for (unsigned int i=0; i<dimensions; i++)
 525                                 {
 526                                         // x lines
 527                                         glVertex3f(from.x, from.y + i * cell_size.y, from.z + j * cell_size.z);
 528                                         glVertex3f(  to.x, from.y + i * cell_size.y, from.z + j * cell_size.z);
 529
 530                                         // y lines
 531                                         glVertex3f(from.x + i * cell_size.x, from.y, from.z + j * cell_size.z);
 532                                         glVertex3f(from.x + i * cell_size.x,   to.y, from.z + j * cell_size.z);
 533
 534                                         // z lines
 535                                         glVertex3f(from.x + i * cell_size.x, from.y + j * cell_size.y, from.z);
 536                                         glVertex3f(from.x + i * cell_size.x, from.y + j * cell_size.y,   to.z);
 537
 538                                 }
 539                         }
 540                 }
 541                 glEnd();
 542         }
 543         else
 544         {
 545                 glBegin(GL_LINES);
 546                 {
 547                         glColor4f(1.0f, 1.0f, 1.0f, 1.0f);
 548
 549                         // x lines
 550                         glVertex3f(from.x, from.y, from.z); glVertex3f(  to.x, from.y, from.z);
 551                         glVertex3f(from.x,   to.y, from.z); glVertex3f(  to.x,   to.y, from.z);
 552                         glVertex3f(from.x,   to.y,   to.z); glVertex3f(  to.x,   to.y,   to.z);
 553                         glVertex3f(from.x, from.y,   to.z); glVertex3f(  to.x, from.y,   to.z);
 554
 555                         // y lines
 556                         glVertex3f(from.x, from.y, from.z); glVertex3f(from.x,   to.y, from.z);
 557                         glVertex3f(  to.x, from.y, from.z); glVertex3f(  to.x,   to.y, from.z);
 558                         glVertex3f(  to.x, from.y,   to.z); glVertex3f(  to.x,   to.y,   to.z);
 559                         glVertex3f(from.x, from.y,   to.z); glVertex3f(from.x,   to.y,   to.z);
 560
 561                         // z lines
 562                         glVertex3f(from.x, from.y, from.z); glVertex3f(from.x, from.y,   to.z);
 563                         glVertex3f(  to.x, from.y, from.z); glVertex3f(  to.x, from.y,   to.z);
 564                         glVertex3f(  to.x,   to.y, from.z); glVertex3f(  to.x,   to.y,   to.z);
 565                         glVertex3f(from.x,   to.y, from.z); glVertex3f(from.x,   to.y,   to.z);
 566
 567                 }
 568                 glEnd();
 569         }
 570         set_lighting(true);
 571 }
 572
 573 Vector3 ScalarField::get_from()
 574 {
 575         return this->from;
 576 }
 577
 578 Vector3 ScalarField::get_to()
 579 {
 580         return this->to;
 581 }
 582
 583 // Evaluators
 584 void ScalarField::set_evaluator(scalar_t (*evaluate) (const Vector3 &vec, scalar_t t))
 585 {
 586         this->evaluate = evaluate;
 587 }
 588
 589 void ScalarField::set_normal_evaluator(Vector3 (*get_normal) (const Vector3 &vec, scalar_t t))
 590 {
 591         this->get_normal = get_normal;
 592 }
 593
 594 // last but not least
 595 void ScalarField::triangulate(TriMesh *mesh, scalar_t isolevel, scalar_t t, bool calc_normals)
 596 {
 597         // Reset mesh and edges table
 598         clear();
 599
 600         // Evaluate
 601         evaluate_all(t);
 602
 603         // triangulate
 604         for (unsigned int z=0; z<dimensions-1; z++)
 605         {
 606                 for (unsigned int y=0; y<dimensions-1; y++)
 607                 {
 608                         for (unsigned int x=0; x<dimensions-1; x++)
 609                         {
 610                                 process_cell(x, y, z, isolevel);
 611                         }
 612                 }
 613         }
 614
 615         // Generate TriMesh
 616         Vertex *varray = new Vertex[verts.size()];
 617         Triangle *tarray = new Triangle[tris.size()];
 618
 619         for (unsigned int i=0; i<verts.size(); i++)
 620         {
 621                 varray[i] = verts[i];
 622         }
 623
 624         for (unsigned int i=0; i<tris.size(); i++)
 625         {
 626                 tarray[i] = tris[i];
 627         }
 628
 629         bool need_normals = calc_normals;
 630
 631         // calculate normals if needed
 632         if(need_normals) {
 633                 unsigned int vsz = verts.size();
 634                 if(get_normal) {
 635                         for(unsigned int i=0; i<vsz; i++) {
 636                                 varray[i].normal = get_normal(varray[i].pos, t);
 637                         }
 638                         need_normals = false;
 639                 } else if(evaluate) {
 640                         for(unsigned int i=0; i<vsz; i++) {
 641                                 varray[i].normal = def_eval_normals(varray[i].pos, t);
 642                         }
 643                         need_normals = false;
 644                 }
 645         }
 646
 647         mesh->set_data(varray, verts.size(), tarray, tris.size());
 648
 649         delete [] varray;
 650         delete [] tarray;
 651
 652         // as a final resort, if we could not calculate normals any other way
 653         // use the regular mesh normal calculation function.
 654         if(need_normals) {
 655                 mesh->calculate_normals_by_index();
 656         }
 657 }
 658
 659