2 This file is part of the 3dengfx, realtime visualization system.
4 Copyright (c) 2004, 2005 John Tsiombikas <nuclear@siggraph.org>
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.
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.
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
21 /* Scalar fields and polygonization
23 * Author: Mihalis Georgoulopoulos 2005
26 #define SCFIELD_SOURCE
27 #include "mcube_tables.h"
28 #include "scfield.hpp"
29 #include "3dengfx/3denginefx.hpp"
32 #define EDGE_NOT_ASSOCIATED 0xFFFFFFFF
41 * adds a vertex and returns its index
43 unsigned int ScalarField::add_vertex(const Vertex &vert)
45 verts.push_back(vert);
46 return verts.size() - 1;
51 * clears std::vector's that hold mesh data and resets edges table
53 void ScalarField::clear()
55 verts.erase(verts.begin(), verts.end());
56 tris.erase(tris.begin(), tris.end());
58 unsigned int num_bytes = dimensions * dimensions * dimensions * sizeof(unsigned int);
60 memset(edges_x, 0xFF, num_bytes);
61 memset(edges_y, 0xFF, num_bytes);
62 memset(edges_z, 0xFF, num_bytes);
67 * Evaluates all values with the external Evaluate function (if specified)
69 void ScalarField::evaluate_all(scalar_t t)
76 for (unsigned int z=0; z<dimensions; z++)
78 for (unsigned int y=0; y<dimensions; y++)
80 for (unsigned int x=0; x<dimensions; x++)
82 set_value(x, y, z, evaluate(get_position(x, y, z), t));
92 void ScalarField::process_cell(int x, int y, int z, scalar_t isolevel)
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;
104 int edge_index = cube_edge_flags[cube_index];
106 scalar_t p , val1, val2;
109 if ( (edge_index & 1) && (get_edge(x, y, z, 0) == EDGE_NOT_ASSOCIATED) )
111 val1 = get_value(x, y, z, 0);
112 val2 = get_value(x, y, z, 1);
113 p = (isolevel - val1) / (val2 - val1);
115 vec1 = get_position(x, y, z, 0);
116 vec2 = get_position(x, y, z, 1);
118 set_edge(x, y, z, 0, add_vertex( Vertex( vec1 + p * (vec2 - vec1) ) ) );
121 if ( (edge_index & 2) && (get_edge(x, y, z, 1) == EDGE_NOT_ASSOCIATED) )
123 val1 = get_value(x, y, z, 1);
124 val2 = get_value(x, y, z, 2);
125 p = (isolevel - val1) / (val2 - val1);
127 vec1 = get_position(x, y, z, 1);
128 vec2 = get_position(x, y, z, 2);
130 set_edge(x, y, z, 1, add_vertex( Vertex( vec1 + p * (vec2 - vec1) ) ) );
133 if ( (edge_index & 4) && (get_edge(x, y, z, 2) == EDGE_NOT_ASSOCIATED) )
135 val1 = get_value(x, y, z, 2);
136 val2 = get_value(x, y, z, 3);
137 p = (isolevel - val1) / (val2 - val1);
139 vec1 = get_position(x, y, z, 2);
140 vec2 = get_position(x, y, z, 3);
142 set_edge(x, y, z, 2, add_vertex( Vertex( vec1 + p * (vec2 - vec1) ) ) );
145 if ( (edge_index & 8) && (get_edge(x, y, z, 3) == EDGE_NOT_ASSOCIATED) )
147 val1 = get_value(x, y, z, 3);
148 val2 = get_value(x, y, z, 0);
149 p = (isolevel - val1) / (val2 - val1);
151 vec1 = get_position(x, y, z, 0);
152 vec2 = get_position(x, y, z, 1);
154 set_edge(x, y, z, 3, add_vertex( Vertex( vec1 + p * (vec2 - vec1) ) ) );
157 if ( (edge_index & 16) && (get_edge(x, y, z, 4) == EDGE_NOT_ASSOCIATED) )
159 val1 = get_value(x, y, z, 4);
160 val2 = get_value(x, y, z, 5);
161 p = (isolevel - val1) / (val2 - val1);
163 vec1 = get_position(x, y, z, 4);
164 vec2 = get_position(x, y, z, 5);
166 set_edge(x, y, z, 4, add_vertex( Vertex( vec1 + p * (vec2 - vec1) ) ) );
169 if ( (edge_index & 32) && (get_edge(x, y, z, 5) == EDGE_NOT_ASSOCIATED) )
171 val1 = get_value(x, y, z, 5);
172 val2 = get_value(x, y, z, 6);
173 p = (isolevel - val1) / (val2 - val1);
175 vec1 = get_position(x, y, z, 5);
176 vec2 = get_position(x, y, z, 6);
178 set_edge(x, y, z, 5, add_vertex( Vertex( vec1 + p * (vec2 - vec1) ) ) );
181 if ( (edge_index & 64) && (get_edge(x, y, z, 6) == EDGE_NOT_ASSOCIATED) )
183 val1 = get_value(x, y, z, 6);
184 val2 = get_value(x, y, z, 7);
185 p = (isolevel - val1) / (val2 - val1);
187 vec1 = get_position(x, y, z, 6);
188 vec2 = get_position(x, y, z, 7);
190 set_edge(x, y, z, 6, add_vertex( Vertex( vec1 + p * (vec2 - vec1) ) ) );
193 if ( (edge_index & 128) && (get_edge(x, y, z, 7) == EDGE_NOT_ASSOCIATED) )
195 val1 = get_value(x, y, z, 7);
196 val2 = get_value(x, y, z, 4);
197 p = (isolevel - val1) / (val2 - val1);
199 vec1 = get_position(x, y, z, 7);
200 vec2 = get_position(x, y, z, 4);
202 set_edge(x, y, z, 7, add_vertex( Vertex( vec1 + p * (vec2 - vec1) ) ) );
205 if ( (edge_index & 256) && (get_edge(x, y, z, 8) == EDGE_NOT_ASSOCIATED) )
207 val1 = get_value(x, y, z, 0);
208 val2 = get_value(x, y, z, 4);
209 p = (isolevel - val1) / (val2 - val1);
211 vec1 = get_position(x, y, z, 0);
212 vec2 = get_position(x, y, z, 4);
214 set_edge(x, y, z, 8, add_vertex( Vertex( vec1 + p * (vec2 - vec1) ) ) );
217 if ( (edge_index & 512) && (get_edge(x, y, z, 9) == EDGE_NOT_ASSOCIATED) )
219 val1 = get_value(x, y, z, 1);
220 val2 = get_value(x, y, z, 5);
221 p = (isolevel - val1) / (val2 - val1);
223 vec1 = get_position(x, y, z, 1);
224 vec2 = get_position(x, y, z, 5);
226 set_edge(x, y, z, 9, add_vertex( Vertex( vec1 + p * (vec2 - vec1) ) ) );
229 if ( (edge_index & 1024) && (get_edge(x, y, z, 10) == EDGE_NOT_ASSOCIATED) )
231 val1 = get_value(x, y, z, 2);
232 val2 = get_value(x, y, z, 6);
233 p = (isolevel - val1) / (val2 - val1);
235 vec1 = get_position(x, y, z, 2);
236 vec2 = get_position(x, y, z, 6);
238 set_edge(x, y, z, 10, add_vertex( Vertex( vec1 + p * (vec2 - vec1) ) ) );
241 if ( (edge_index & 2048) && (get_edge(x, y, z, 11) == EDGE_NOT_ASSOCIATED) )
243 val1 = get_value(x, y, z, 3);
244 val2 = get_value(x, y, z, 7);
245 p = (isolevel - val1) / (val2 - val1);
247 vec1 = get_position(x, y, z, 3);
248 vec2 = get_position(x, y, z, 7);
250 set_edge(x, y, z, 11, add_vertex( Vertex( vec1 + p * (vec2 - vec1) ) ) );
254 unsigned int p1, p2, p3;
256 if (tri_table[cube_index][0] != -1)
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));
264 if (tri_table[cube_index][3] != -1)
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));
272 if (tri_table[cube_index][6] != -1)
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));
280 if (tri_table[cube_index][9] != -1)
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));
288 if (tri_table[cube_index][12] != -1)
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));
299 * returns the index to the values array for the specified coords
301 unsigned int ScalarField::get_value_index(int x, int y, int z)
303 return x + y * dimensions + z * dimensions * dimensions;
307 Vector3 ScalarField::def_eval_normals(const Vector3 &vec, scalar_t t) {
308 if(!evaluate) return Vector3(0, 0, 0);
310 Vector3 diff = cell_size * 0.25;
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);
317 return grad.normalized();
326 ScalarField::ScalarField()
329 edges_x = edges_y = edges_z = 0;
331 from = to = cell_size = Vector3(0, 0, 0);
336 ScalarField::ScalarField(unsigned int dimensions, const Vector3 &from, const Vector3 &to)
338 this->dimensions = dimensions;
341 this->cell_size = (to - from) / (dimensions - 1);
344 edges_x = edges_y = edges_z = 0;
349 set_dimensions(dimensions);
352 ScalarField::~ScalarField()
364 void ScalarField::set_dimensions(unsigned int dimensions)
366 this->dimensions = dimensions;
376 values = new scalar_t [dimensions * dimensions * dimensions];
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];
387 void ScalarField::set_value(int x, int y, int z, scalar_t value)
389 values[get_value_index(x, y, z)] = value;
392 scalar_t ScalarField::get_value(int x, int y, int z)
394 return values[get_value_index(x, y, z)];
398 // get / set relative to cell
399 void ScalarField::set_value(int cx, int cy, int cz, int vert_index, scalar_t value)
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);
412 scalar_t ScalarField::get_value(int cx, int cy, int cz, int vert_index)
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);
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)
430 unsigned int d = dimensions;
431 unsigned int d2 = dimensions * dimensions;
434 edges_x[cx + 0 + (cy + 0) * d + (cz + 1) * d2] = index;
436 edges_z[cx + 1 + (cy + 0) * d + (cz + 0) * d2] = index;
438 edges_x[cx + 0 + (cy + 0) * d + (cz + 0) * d2] = index;
440 edges_z[cx + 0 + (cy + 0) * d + (cz + 0) * d2] = index;
442 edges_x[cx + 0 + (cy + 1) * d + (cz + 1) * d2] = index;
444 edges_z[cx + 1 + (cy + 1) * d + (cz + 0) * d2] = index;
446 edges_x[cx + 0 + (cy + 1) * d + (cz + 0) * d2] = index;
448 edges_z[cx + 0 + (cy + 1) * d + (cz + 0) * d2] = index;
450 edges_y[cx + 0 + (cy + 0) * d + (cz + 1) * d2] = index;
452 edges_y[cx + 1 + (cy + 0) * d + (cz + 1) * d2] = index;
454 edges_y[cx + 1 + (cy + 0) * d + (cz + 0) * d2] = index;
456 edges_y[cx + 0 + (cy + 0) * d + (cz + 0) * d2] = index;
459 unsigned int ScalarField::get_edge(int cx, int cy, int cz, int edge)
461 unsigned int d = dimensions;
462 unsigned int d2 = dimensions * dimensions;
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];
481 Vector3 ScalarField::get_position(int x, int y, int z)
484 vx = from.x + cell_size.x * x;
485 vy = from.y + cell_size.y * y;
486 vz = from.z + cell_size.z * z;
488 return Vector3(vx, vy, vz);
491 Vector3 ScalarField::get_position(int cx, int cy, int cz, int vert_index)
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);
502 return Vector3(0, 0, 0);
505 void ScalarField::set_from_to(const Vector3 &from, const Vector3 &to)
509 this->cell_size = (to - from) / (dimensions - 1);
512 void ScalarField::draw_field(bool full)
520 glColor4f(1.0f, 1.0f, 1.0f, 1.0f);
522 for (unsigned int j=0; j<dimensions; j++)
524 for (unsigned int i=0; i<dimensions; i++)
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);
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);
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);
547 glColor4f(1.0f, 1.0f, 1.0f, 1.0f);
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);
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);
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);
573 Vector3 ScalarField::get_from()
578 Vector3 ScalarField::get_to()
584 void ScalarField::set_evaluator(scalar_t (*evaluate) (const Vector3 &vec, scalar_t t))
586 this->evaluate = evaluate;
589 void ScalarField::set_normal_evaluator(Vector3 (*get_normal) (const Vector3 &vec, scalar_t t))
591 this->get_normal = get_normal;
594 // last but not least
595 void ScalarField::triangulate(TriMesh *mesh, scalar_t isolevel, scalar_t t, bool calc_normals)
597 // Reset mesh and edges table
604 for (unsigned int z=0; z<dimensions-1; z++)
606 for (unsigned int y=0; y<dimensions-1; y++)
608 for (unsigned int x=0; x<dimensions-1; x++)
610 process_cell(x, y, z, isolevel);
616 Vertex *varray = new Vertex[verts.size()];
617 Triangle *tarray = new Triangle[tris.size()];
619 for (unsigned int i=0; i<verts.size(); i++)
621 varray[i] = verts[i];
624 for (unsigned int i=0; i<tris.size(); i++)
629 bool need_normals = calc_normals;
631 // calculate normals if needed
633 unsigned int vsz = verts.size();
635 for(unsigned int i=0; i<vsz; i++) {
636 varray[i].normal = get_normal(varray[i].pos, t);
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);
643 need_normals = false;
647 mesh->set_data(varray, verts.size(), tarray, tris.size());
652 // as a final resort, if we could not calculate normals any other way
653 // use the regular mesh normal calculation function.
655 mesh->calculate_normals_by_index();