////////////////////////////////////////////////////////////////////////////// // // Author : Josh Grant // Date : November 26th, 2001 // File : MarchingCubes.cpp // Description : Header file defining the MarchingCubes class // ////////////////////////////////////////////////////////////////////////////// #include #include "MarchingCubes.h" ////////////////////////////////////////////////////////////////////////////// // // MarchingCubes class // // All tables were based off of code written by Cory Bloyd // (corysama@yahoo.com), from a nice description of the Marching // Cubes algorithm located at // http://www.swin.edu.au/astronomy/pbourke/modelling/polygonise/ // ////////////////////////////////////////////////////////////////////////////// // the blending percentage must be less than 0.5 #define MC_MAX_BLEND_PERCENT 0.499999f // blending technique added by Gokhan Kisacikoglu namespace blend { // The blending is nothing but a simple test to the closest voxel // corner. If we are too close, we accumulate the point and // average it later. The accumulated point will be kept in the // 'verts' for conservative memory management, but we have to // save at least its index and the number of accumulated points. // class point { public: point() : m_location(0, 0, 0), m_normal(0, 0, 0), m_vertexIndex(-1), m_hits(0) {} void set( int _vertexIndex, const SbVec3f &_location, const SbVec3f &_normal ) { m_vertexIndex = _vertexIndex; m_location = _location; m_normal = _normal; m_hits = 1; } int add( const SbVec3f &_location, const SbVec3f &_normal ) { m_location += _location; m_normal += _normal; ++ m_hits; return m_vertexIndex; } int hits() const { return m_hits; } SbVec3f location() const { return m_location / float(m_hits); } SbVec3f normal() const { return m_normal / float(m_hits); } int vertexIndex() const { return m_vertexIndex; } private: SbVec3f m_location, m_normal; int m_vertexIndex, m_hits; }; } // For each of the possible vertex states listed in this table there is a // specific triangulation of the edge intersection points. The table lists // all of them in the form of 0-5 edge triples with the list terminated by // the invalid value -1. For example: triTable[3] list the 2 triangles // formed when cube[0] and cube[1] are inside of the surface, but the rest of // the cube is not. static const short triTable[256][16] = { {-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, { 0, 8, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, { 0, 1, 9, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, { 1, 8, 3, 9, 8, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, { 1, 2, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, { 0, 8, 3, 1, 2, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, { 9, 2, 10, 0, 2, 9, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, { 2, 8, 3, 2, 10, 8, 10, 9, 8, -1, -1, -1, -1, -1, -1, -1}, { 3, 11, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, { 0, 11, 2, 8, 11, 0, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, { 1, 9, 0, 2, 3, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, { 1, 11, 2, 1, 9, 11, 9, 8, 11, -1, -1, -1, -1, -1, -1, -1}, { 3, 10, 1, 11, 10, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, { 0, 10, 1, 0, 8, 10, 8, 11, 10, -1, -1, -1, -1, -1, -1, -1}, { 3, 9, 0, 3, 11, 9, 11, 10, 9, -1, -1, -1, -1, -1, -1, -1}, { 9, 8, 10, 10, 8, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, { 4, 7, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, { 4, 3, 0, 7, 3, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, { 0, 1, 9, 8, 4, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, { 4, 1, 9, 4, 7, 1, 7, 3, 1, -1, -1, -1, -1, -1, -1, -1}, { 1, 2, 10, 8, 4, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, { 3, 4, 7, 3, 0, 4, 1, 2, 10, -1, -1, -1, -1, -1, -1, -1}, { 9, 2, 10, 9, 0, 2, 8, 4, 7, -1, -1, -1, -1, -1, -1, -1}, { 2, 10, 9, 2, 9, 7, 2, 7, 3, 7, 9, 4, -1, -1, -1, -1}, { 8, 4, 7, 3, 11, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, {11, 4, 7, 11, 2, 4, 2, 0, 4, -1, -1, -1, -1, -1, -1, -1}, { 9, 0, 1, 8, 4, 7, 2, 3, 11, -1, -1, -1, -1, -1, -1, -1}, { 4, 7, 11, 9, 4, 11, 9, 11, 2, 9, 2, 1, -1, -1, -1, -1}, { 3, 10, 1, 3, 11, 10, 7, 8, 4, -1, -1, -1, -1, -1, -1, -1}, { 1, 11, 10, 1, 4, 11, 1, 0, 4, 7, 11, 4, -1, -1, -1, -1}, { 4, 7, 8, 9, 0, 11, 9, 11, 10, 11, 0, 3, -1, -1, -1, -1}, { 4, 7, 11, 4, 11, 9, 9, 11, 10, -1, -1, -1, -1, -1, -1, -1}, { 9, 5, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, { 9, 5, 4, 0, 8, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, { 0, 5, 4, 1, 5, 0, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, { 8, 5, 4, 8, 3, 5, 3, 1, 5, -1, -1, -1, -1, -1, -1, -1}, { 1, 2, 10, 9, 5, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, { 3, 0, 8, 1, 2, 10, 4, 9, 5, -1, -1, -1, -1, -1, -1, -1}, { 5, 2, 10, 5, 4, 2, 4, 0, 2, -1, -1, -1, -1, -1, -1, -1}, { 2, 10, 5, 3, 2, 5, 3, 5, 4, 3, 4, 8, -1, -1, -1, -1}, { 9, 5, 4, 2, 3, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, { 0, 11, 2, 0, 8, 11, 4, 9, 5, -1, -1, -1, -1, -1, -1, -1}, { 0, 5, 4, 0, 1, 5, 2, 3, 11, -1, -1, -1, -1, -1, -1, -1}, { 2, 1, 5, 2, 5, 8, 2, 8, 11, 4, 8, 5, -1, -1, -1, -1}, {10, 3, 11, 10, 1, 3, 9, 5, 4, -1, -1, -1, -1, -1, -1, -1}, { 4, 9, 5, 0, 8, 1, 8, 10, 1, 8, 11, 10, -1, -1, -1, -1}, { 5, 4, 0, 5, 0, 11, 5, 11, 10, 11, 0, 3, -1, -1, -1, -1}, { 5, 4, 8, 5, 8, 10, 10, 8, 11, -1, -1, -1, -1, -1, -1, -1}, { 9, 7, 8, 5, 7, 9, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, { 9, 3, 0, 9, 5, 3, 5, 7, 3, -1, -1, -1, -1, -1, -1, -1}, { 0, 7, 8, 0, 1, 7, 1, 5, 7, -1, -1, -1, -1, -1, -1, -1}, { 1, 5, 3, 3, 5, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, { 9, 7, 8, 9, 5, 7, 10, 1, 2, -1, -1, -1, -1, -1, -1, -1}, {10, 1, 2, 9, 5, 0, 5, 3, 0, 5, 7, 3, -1, -1, -1, -1}, { 8, 0, 2, 8, 2, 5, 8, 5, 7, 10, 5, 2, -1, -1, -1, -1}, { 2, 10, 5, 2, 5, 3, 3, 5, 7, -1, -1, -1, -1, -1, -1, -1}, { 7, 9, 5, 7, 8, 9, 3, 11, 2, -1, -1, -1, -1, -1, -1, -1}, { 9, 5, 7, 9, 7, 2, 9, 2, 0, 2, 7, 11, -1, -1, -1, -1}, { 2, 3, 11, 0, 1, 8, 1, 7, 8, 1, 5, 7, -1, -1, -1, -1}, {11, 2, 1, 11, 1, 7, 7, 1, 5, -1, -1, -1, -1, -1, -1, -1}, { 9, 5, 8, 8, 5, 7, 10, 1, 3, 10, 3, 11, -1, -1, -1, -1}, { 5, 7, 0, 5, 0, 9, 7, 11, 0, 1, 0, 10, 11, 10, 0, -1}, {11, 10, 0, 11, 0, 3, 10, 5, 0, 8, 0, 7, 5, 7, 0, -1}, {11, 10, 5, 7, 11, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, {10, 6, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, { 0, 8, 3, 5, 10, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, { 9, 0, 1, 5, 10, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, { 1, 8, 3, 1, 9, 8, 5, 10, 6, -1, -1, -1, -1, -1, -1, -1}, { 1, 6, 5, 2, 6, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, { 1, 6, 5, 1, 2, 6, 3, 0, 8, -1, -1, -1, -1, -1, -1, -1}, { 9, 6, 5, 9, 0, 6, 0, 2, 6, -1, -1, -1, -1, -1, -1, -1}, { 5, 9, 8, 5, 8, 2, 5, 2, 6, 3, 2, 8, -1, -1, -1, -1}, { 2, 3, 11, 10, 6, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, {11, 0, 8, 11, 2, 0, 10, 6, 5, -1, -1, -1, -1, -1, -1, -1}, { 0, 1, 9, 2, 3, 11, 5, 10, 6, -1, -1, -1, -1, -1, -1, -1}, { 5, 10, 6, 1, 9, 2, 9, 11, 2, 9, 8, 11, -1, -1, -1, -1}, { 6, 3, 11, 6, 5, 3, 5, 1, 3, -1, -1, -1, -1, -1, -1, -1}, { 0, 8, 11, 0, 11, 5, 0, 5, 1, 5, 11, 6, -1, -1, -1, -1}, { 3, 11, 6, 0, 3, 6, 0, 6, 5, 0, 5, 9, -1, -1, -1, -1}, { 6, 5, 9, 6, 9, 11, 11, 9, 8, -1, -1, -1, -1, -1, -1, -1}, { 5, 10, 6, 4, 7, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, { 4, 3, 0, 4, 7, 3, 6, 5, 10, -1, -1, -1, -1, -1, -1, -1}, { 1, 9, 0, 5, 10, 6, 8, 4, 7, -1, -1, -1, -1, -1, -1, -1}, {10, 6, 5, 1, 9, 7, 1, 7, 3, 7, 9, 4, -1, -1, -1, -1}, { 6, 1, 2, 6, 5, 1, 4, 7, 8, -1, -1, -1, -1, -1, -1, -1}, { 1, 2, 5, 5, 2, 6, 3, 0, 4, 3, 4, 7, -1, -1, -1, -1}, { 8, 4, 7, 9, 0, 5, 0, 6, 5, 0, 2, 6, -1, -1, -1, -1}, { 7, 3, 9, 7, 9, 4, 3, 2, 9, 5, 9, 6, 2, 6, 9, -1}, { 3, 11, 2, 7, 8, 4, 10, 6, 5, -1, -1, -1, -1, -1, -1, -1}, { 5, 10, 6, 4, 7, 2, 4, 2, 0, 2, 7, 11, -1, -1, -1, -1}, { 0, 1, 9, 4, 7, 8, 2, 3, 11, 5, 10, 6, -1, -1, -1, -1}, { 9, 2, 1, 9, 11, 2, 9, 4, 11, 7, 11, 4, 5, 10, 6, -1}, { 8, 4, 7, 3, 11, 5, 3, 5, 1, 5, 11, 6, -1, -1, -1, -1}, { 5, 1, 11, 5, 11, 6, 1, 0, 11, 7, 11, 4, 0, 4, 11, -1}, { 0, 5, 9, 0, 6, 5, 0, 3, 6, 11, 6, 3, 8, 4, 7, -1}, { 6, 5, 9, 6, 9, 11, 4, 7, 9, 7, 11, 9, -1, -1, -1, -1}, {10, 4, 9, 6, 4, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, { 4, 10, 6, 4, 9, 10, 0, 8, 3, -1, -1, -1, -1, -1, -1, -1}, {10, 0, 1, 10, 6, 0, 6, 4, 0, -1, -1, -1, -1, -1, -1, -1}, { 8, 3, 1, 8, 1, 6, 8, 6, 4, 6, 1, 10, -1, -1, -1, -1}, { 1, 4, 9, 1, 2, 4, 2, 6, 4, -1, -1, -1, -1, -1, -1, -1}, { 3, 0, 8, 1, 2, 9, 2, 4, 9, 2, 6, 4, -1, -1, -1, -1}, { 0, 2, 4, 4, 2, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, { 8, 3, 2, 8, 2, 4, 4, 2, 6, -1, -1, -1, -1, -1, -1, -1}, {10, 4, 9, 10, 6, 4, 11, 2, 3, -1, -1, -1, -1, -1, -1, -1}, { 0, 8, 2, 2, 8, 11, 4, 9, 10, 4, 10, 6, -1, -1, -1, -1}, { 3, 11, 2, 0, 1, 6, 0, 6, 4, 6, 1, 10, -1, -1, -1, -1}, { 6, 4, 1, 6, 1, 10, 4, 8, 1, 2, 1, 11, 8, 11, 1, -1}, { 9, 6, 4, 9, 3, 6, 9, 1, 3, 11, 6, 3, -1, -1, -1, -1}, { 8, 11, 1, 8, 1, 0, 11, 6, 1, 9, 1, 4, 6, 4, 1, -1}, { 3, 11, 6, 3, 6, 0, 0, 6, 4, -1, -1, -1, -1, -1, -1, -1}, { 6, 4, 8, 11, 6, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, { 7, 10, 6, 7, 8, 10, 8, 9, 10, -1, -1, -1, -1, -1, -1, -1}, { 0, 7, 3, 0, 10, 7, 0, 9, 10, 6, 7, 10, -1, -1, -1, -1}, {10, 6, 7, 1, 10, 7, 1, 7, 8, 1, 8, 0, -1, -1, -1, -1}, {10, 6, 7, 10, 7, 1, 1, 7, 3, -1, -1, -1, -1, -1, -1, -1}, { 1, 2, 6, 1, 6, 8, 1, 8, 9, 8, 6, 7, -1, -1, -1, -1}, { 2, 6, 9, 2, 9, 1, 6, 7, 9, 0, 9, 3, 7, 3, 9, -1}, { 7, 8, 0, 7, 0, 6, 6, 0, 2, -1, -1, -1, -1, -1, -1, -1}, { 7, 3, 2, 6, 7, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, { 2, 3, 11, 10, 6, 8, 10, 8, 9, 8, 6, 7, -1, -1, -1, -1}, { 2, 0, 7, 2, 7, 11, 0, 9, 7, 6, 7, 10, 9, 10, 7, -1}, { 1, 8, 0, 1, 7, 8, 1, 10, 7, 6, 7, 10, 2, 3, 11, -1}, {11, 2, 1, 11, 1, 7, 10, 6, 1, 6, 7, 1, -1, -1, -1, -1}, { 8, 9, 6, 8, 6, 7, 9, 1, 6, 11, 6, 3, 1, 3, 6, -1}, { 0, 9, 1, 11, 6, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, { 7, 8, 0, 7, 0, 6, 3, 11, 0, 11, 6, 0, -1, -1, -1, -1}, { 7, 11, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, { 7, 6, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, { 3, 0, 8, 11, 7, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, { 0, 1, 9, 11, 7, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, { 8, 1, 9, 8, 3, 1, 11, 7, 6, -1, -1, -1, -1, -1, -1, -1}, {10, 1, 2, 6, 11, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, { 1, 2, 10, 3, 0, 8, 6, 11, 7, -1, -1, -1, -1, -1, -1, -1}, { 2, 9, 0, 2, 10, 9, 6, 11, 7, -1, -1, -1, -1, -1, -1, -1}, { 6, 11, 7, 2, 10, 3, 10, 8, 3, 10, 9, 8, -1, -1, -1, -1}, { 7, 2, 3, 6, 2, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, { 7, 0, 8, 7, 6, 0, 6, 2, 0, -1, -1, -1, -1, -1, -1, -1}, { 2, 7, 6, 2, 3, 7, 0, 1, 9, -1, -1, -1, -1, -1, -1, -1}, { 1, 6, 2, 1, 8, 6, 1, 9, 8, 8, 7, 6, -1, -1, -1, -1}, {10, 7, 6, 10, 1, 7, 1, 3, 7, -1, -1, -1, -1, -1, -1, -1}, {10, 7, 6, 1, 7, 10, 1, 8, 7, 1, 0, 8, -1, -1, -1, -1}, { 0, 3, 7, 0, 7, 10, 0, 10, 9, 6, 10, 7, -1, -1, -1, -1}, { 7, 6, 10, 7, 10, 8, 8, 10, 9, -1, -1, -1, -1, -1, -1, -1}, { 6, 8, 4, 11, 8, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, { 3, 6, 11, 3, 0, 6, 0, 4, 6, -1, -1, -1, -1, -1, -1, -1}, { 8, 6, 11, 8, 4, 6, 9, 0, 1, -1, -1, -1, -1, -1, -1, -1}, { 9, 4, 6, 9, 6, 3, 9, 3, 1, 11, 3, 6, -1, -1, -1, -1}, { 6, 8, 4, 6, 11, 8, 2, 10, 1, -1, -1, -1, -1, -1, -1, -1}, { 1, 2, 10, 3, 0, 11, 0, 6, 11, 0, 4, 6, -1, -1, -1, -1}, { 4, 11, 8, 4, 6, 11, 0, 2, 9, 2, 10, 9, -1, -1, -1, -1}, {10, 9, 3, 10, 3, 2, 9, 4, 3, 11, 3, 6, 4, 6, 3, -1}, { 8, 2, 3, 8, 4, 2, 4, 6, 2, -1, -1, -1, -1, -1, -1, -1}, { 0, 4, 2, 4, 6, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, { 1, 9, 0, 2, 3, 4, 2, 4, 6, 4, 3, 8, -1, -1, -1, -1}, { 1, 9, 4, 1, 4, 2, 2, 4, 6, -1, -1, -1, -1, -1, -1, -1}, { 8, 1, 3, 8, 6, 1, 8, 4, 6, 6, 10, 1, -1, -1, -1, -1}, {10, 1, 0, 10, 0, 6, 6, 0, 4, -1, -1, -1, -1, -1, -1, -1}, { 4, 6, 3, 4, 3, 8, 6, 10, 3, 0, 3, 9, 10, 9, 3, -1}, {10, 9, 4, 6, 10, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, { 4, 9, 5, 7, 6, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, { 0, 8, 3, 4, 9, 5, 11, 7, 6, -1, -1, -1, -1, -1, -1, -1}, { 5, 0, 1, 5, 4, 0, 7, 6, 11, -1, -1, -1, -1, -1, -1, -1}, {11, 7, 6, 8, 3, 4, 3, 5, 4, 3, 1, 5, -1, -1, -1, -1}, { 9, 5, 4, 10, 1, 2, 7, 6, 11, -1, -1, -1, -1, -1, -1, -1}, { 6, 11, 7, 1, 2, 10, 0, 8, 3, 4, 9, 5, -1, -1, -1, -1}, { 7, 6, 11, 5, 4, 10, 4, 2, 10, 4, 0, 2, -1, -1, -1, -1}, { 3, 4, 8, 3, 5, 4, 3, 2, 5, 10, 5, 2, 11, 7, 6, -1}, { 7, 2, 3, 7, 6, 2, 5, 4, 9, -1, -1, -1, -1, -1, -1, -1}, { 9, 5, 4, 0, 8, 6, 0, 6, 2, 6, 8, 7, -1, -1, -1, -1}, { 3, 6, 2, 3, 7, 6, 1, 5, 0, 5, 4, 0, -1, -1, -1, -1}, { 6, 2, 8, 6, 8, 7, 2, 1, 8, 4, 8, 5, 1, 5, 8, -1}, { 9, 5, 4, 10, 1, 6, 1, 7, 6, 1, 3, 7, -1, -1, -1, -1}, { 1, 6, 10, 1, 7, 6, 1, 0, 7, 8, 7, 0, 9, 5, 4, -1}, { 4, 0, 10, 4, 10, 5, 0, 3, 10, 6, 10, 7, 3, 7, 10, -1}, { 7, 6, 10, 7, 10, 8, 5, 4, 10, 4, 8, 10, -1, -1, -1, -1}, { 6, 9, 5, 6, 11, 9, 11, 8, 9, -1, -1, -1, -1, -1, -1, -1}, { 3, 6, 11, 0, 6, 3, 0, 5, 6, 0, 9, 5, -1, -1, -1, -1}, { 0, 11, 8, 0, 5, 11, 0, 1, 5, 5, 6, 11, -1, -1, -1, -1}, { 6, 11, 3, 6, 3, 5, 5, 3, 1, -1, -1, -1, -1, -1, -1, -1}, { 1, 2, 10, 9, 5, 11, 9, 11, 8, 11, 5, 6, -1, -1, -1, -1}, { 0, 11, 3, 0, 6, 11, 0, 9, 6, 5, 6, 9, 1, 2, 10, -1}, {11, 8, 5, 11, 5, 6, 8, 0, 5, 10, 5, 2, 0, 2, 5, -1}, { 6, 11, 3, 6, 3, 5, 2, 10, 3, 10, 5, 3, -1, -1, -1, -1}, { 5, 8, 9, 5, 2, 8, 5, 6, 2, 3, 8, 2, -1, -1, -1, -1}, { 9, 5, 6, 9, 6, 0, 0, 6, 2, -1, -1, -1, -1, -1, -1, -1}, { 1, 5, 8, 1, 8, 0, 5, 6, 8, 3, 8, 2, 6, 2, 8, -1}, { 1, 5, 6, 2, 1, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, { 1, 3, 6, 1, 6, 10, 3, 8, 6, 5, 6, 9, 8, 9, 6, -1}, {10, 1, 0, 10, 0, 6, 9, 5, 0, 5, 6, 0, -1, -1, -1, -1}, { 0, 3, 8, 5, 6, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, {10, 5, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, {11, 5, 10, 7, 5, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, {11, 5, 10, 11, 7, 5, 8, 3, 0, -1, -1, -1, -1, -1, -1, -1}, { 5, 11, 7, 5, 10, 11, 1, 9, 0, -1, -1, -1, -1, -1, -1, -1}, {10, 7, 5, 10, 11, 7, 9, 8, 1, 8, 3, 1, -1, -1, -1, -1}, {11, 1, 2, 11, 7, 1, 7, 5, 1, -1, -1, -1, -1, -1, -1, -1}, { 0, 8, 3, 1, 2, 7, 1, 7, 5, 7, 2, 11, -1, -1, -1, -1}, { 9, 7, 5, 9, 2, 7, 9, 0, 2, 2, 11, 7, -1, -1, -1, -1}, { 7, 5, 2, 7, 2, 11, 5, 9, 2, 3, 2, 8, 9, 8, 2, -1}, { 2, 5, 10, 2, 3, 5, 3, 7, 5, -1, -1, -1, -1, -1, -1, -1}, { 8, 2, 0, 8, 5, 2, 8, 7, 5, 10, 2, 5, -1, -1, -1, -1}, { 9, 0, 1, 5, 10, 3, 5, 3, 7, 3, 10, 2, -1, -1, -1, -1}, { 9, 8, 2, 9, 2, 1, 8, 7, 2, 10, 2, 5, 7, 5, 2, -1}, { 1, 3, 5, 3, 7, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, { 0, 8, 7, 0, 7, 1, 1, 7, 5, -1, -1, -1, -1, -1, -1, -1}, { 9, 0, 3, 9, 3, 5, 5, 3, 7, -1, -1, -1, -1, -1, -1, -1}, { 9, 8, 7, 5, 9, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, { 5, 8, 4, 5, 10, 8, 10, 11, 8, -1, -1, -1, -1, -1, -1, -1}, { 5, 0, 4, 5, 11, 0, 5, 10, 11, 11, 3, 0, -1, -1, -1, -1}, { 0, 1, 9, 8, 4, 10, 8, 10, 11, 10, 4, 5, -1, -1, -1, -1}, {10, 11, 4, 10, 4, 5, 11, 3, 4, 9, 4, 1, 3, 1, 4, -1}, { 2, 5, 1, 2, 8, 5, 2, 11, 8, 4, 5, 8, -1, -1, -1, -1}, { 0, 4, 11, 0, 11, 3, 4, 5, 11, 2, 11, 1, 5, 1, 11, -1}, { 0, 2, 5, 0, 5, 9, 2, 11, 5, 4, 5, 8, 11, 8, 5, -1}, { 9, 4, 5, 2, 11, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, { 2, 5, 10, 3, 5, 2, 3, 4, 5, 3, 8, 4, -1, -1, -1, -1}, { 5, 10, 2, 5, 2, 4, 4, 2, 0, -1, -1, -1, -1, -1, -1, -1}, { 3, 10, 2, 3, 5, 10, 3, 8, 5, 4, 5, 8, 0, 1, 9, -1}, { 5, 10, 2, 5, 2, 4, 1, 9, 2, 9, 4, 2, -1, -1, -1, -1}, { 8, 4, 5, 8, 5, 3, 3, 5, 1, -1, -1, -1, -1, -1, -1, -1}, { 0, 4, 5, 1, 0, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, { 8, 4, 5, 8, 5, 3, 9, 0, 5, 0, 3, 5, -1, -1, -1, -1}, { 9, 4, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, { 4, 11, 7, 4, 9, 11, 9, 10, 11, -1, -1, -1, -1, -1, -1, -1}, { 0, 8, 3, 4, 9, 7, 9, 11, 7, 9, 10, 11, -1, -1, -1, -1}, { 1, 10, 11, 1, 11, 4, 1, 4, 0, 7, 4, 11, -1, -1, -1, -1}, { 3, 1, 4, 3, 4, 8, 1, 10, 4, 7, 4, 11, 10, 11, 4, -1}, { 4, 11, 7, 9, 11, 4, 9, 2, 11, 9, 1, 2, -1, -1, -1, -1}, { 9, 7, 4, 9, 11, 7, 9, 1, 11, 2, 11, 1, 0, 8, 3, -1}, {11, 7, 4, 11, 4, 2, 2, 4, 0, -1, -1, -1, -1, -1, -1, -1}, {11, 7, 4, 11, 4, 2, 8, 3, 4, 3, 2, 4, -1, -1, -1, -1}, { 2, 9, 10, 2, 7, 9, 2, 3, 7, 7, 4, 9, -1, -1, -1, -1}, { 9, 10, 7, 9, 7, 4, 10, 2, 7, 8, 7, 0, 2, 0, 7, -1}, { 3, 7, 10, 3, 10, 2, 7, 4, 10, 1, 10, 0, 4, 0, 10, -1}, { 1, 10, 2, 8, 7, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, { 4, 9, 1, 4, 1, 7, 7, 1, 3, -1, -1, -1, -1, -1, -1, -1}, { 4, 9, 1, 4, 1, 7, 0, 8, 1, 8, 7, 1, -1, -1, -1, -1}, { 4, 0, 3, 7, 4, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, { 4, 8, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, { 9, 10, 8, 10, 11, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, { 3, 0, 9, 3, 9, 11, 11, 9, 10, -1, -1, -1, -1, -1, -1, -1}, { 0, 1, 10, 0, 10, 8, 8, 10, 11, -1, -1, -1, -1, -1, -1, -1}, { 3, 1, 10, 11, 3, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, { 1, 2, 11, 1, 11, 9, 9, 11, 8, -1, -1, -1, -1, -1, -1, -1}, { 3, 0, 9, 3, 9, 11, 1, 2, 9, 2, 11, 9, -1, -1, -1, -1}, { 0, 2, 11, 8, 0, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, { 3, 2, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, { 2, 3, 8, 2, 8, 10, 10, 8, 9, -1, -1, -1, -1, -1, -1, -1}, { 9, 10, 2, 0, 9, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, { 2, 3, 8, 2, 8, 10, 0, 1, 8, 1, 10, 8, -1, -1, -1, -1}, { 1, 10, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, { 1, 3, 8, 9, 1, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, { 0, 9, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, { 0, 3, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, {-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1} }; // offsets used to determine which vertex of a cube has a value >= to the // isoval. static const int cubeVertOffset[8][3] = { { 0, 0, 0}, {-1, 0, 0}, {-1, -1, 0}, { 0, -1, 0}, { 0, 0, -1}, {-1, 0, -1}, {-1, -1, -1}, { 0, -1, -1} }; // Inventor stuff SO_ENGINE_SOURCE(MarchingCubes); // Inventor stuff void MarchingCubes::initClass () { SO_ENGINE_INIT_CLASS(MarchingCubes, SoEngine, "Engine"); } ////////////////////////////////////////////////////////////////////////////// // // Description: // Default constructor initializes all class fields // ////////////////////////////////////////////////////////////////////////////// MarchingCubes::MarchingCubes () { SO_ENGINE_CONSTRUCTOR(MarchingCubes); dims[0] = dims[1] = dims[2] = 0; SO_ENGINE_ADD_INPUT(data, (dims, 0, false)); SO_ENGINE_ADD_INPUT(isoValue, (1.0)); SO_ENGINE_ADD_INPUT(blendPercent, (MC_MAX_BLEND_PERCENT)); SO_ENGINE_ADD_OUTPUT(points, SoMFVec3f); SO_ENGINE_ADD_OUTPUT(normals, SoMFVec3f); SO_ENGINE_ADD_OUTPUT(indexes, SoMFInt32); SO_ENGINE_ADD_OUTPUT(numTriangles, SoSFFloat); computeGradients = true; // surface will be created with coordinates between 0 and 1, making it // easier to translate and scale min = SbVec3f(0.0, 0.0, 0.0); max = SbVec3f(1.0, 1.0, 1.0); } MarchingCubes::~MarchingCubes () { } ////////////////////////////////////////////////////////////////////////////// // // Description: // Whenever the data field changes set the flag to recompute the // gradients. // ////////////////////////////////////////////////////////////////////////////// void MarchingCubes::inputChanged (SoField *which) { if (which == &data) computeGradients = true; } ////////////////////////////////////////////////////////////////////////////// // // Description: // Called by Inventor each time an output has been called and an input has // been modified. The points and the indexes will be computed here. // ////////////////////////////////////////////////////////////////////////////// void MarchingCubes::evaluate () { //currVert = 0; // reset the arrays holding the vertices and indexes before they are // modified. verts.setNum(0); norms.setNum(0); inds.setNum(0); // grab the current value of isoValue so that only one value will be used // during the computation. isoval = isoValue.getValue(); bperc = fabs(blendPercent.getValue()); bperc = (bperc < 0.5) ? bperc : MC_MAX_BLEND_PERCENT; // calculate the gradient at each grid point, this is called whenever the // data has been modified. if (computeGradients) { // get the latest values values = data.getValue(dims); // compute the size of each voxel (dependent on the dimensions of the // data) SbVec3f size = max - min; voxel[0] = (dims[0] > 1) ? size[0] / (float)(dims[0] - 1) : 0.0; voxel[1] = (dims[1] > 1) ? size[1] / (float)(dims[1] - 1) : 0.0; voxel[2] = (dims[2] > 1) ? size[2] / (float)(dims[2] - 1) : 0.0; // there are actually more than the standard (nx - 1)*(ny - 1)*(nz - 1) // number of cubes for the dataset. To speed computation up by // eliminating constant index bounds checking, a cube buffer is placed // around the actual cubes being used. This way if the indexing goes one // under zero or over the max it won't core dump :) cubeDims[0] = dims[0] + 1; cubeDims[1] = dims[1] + 1; cubeDims[2] = dims[2] + 1; // set the total number of cubes cubes.setNum(cubeDims[0] * cubeDims[1] * cubeDims[2]); // set the total number of edges being looked at. edgeIndex.setNum(3 * dims[0] * dims[1] * dims[2]); // calculate the gradients at each grid point. these gradients will be // used to interpolate the values at the surface points. calcGradients(); computeGradients = false; } // calculate the points on the surface calcPoints(); // create the triangle triples defining the triangles of the surface calcIndexes(); // send the results to the connected fields SO_ENGINE_OUTPUT(points, SoMFVec3f, setValues(0, verts.getNum(), verts.getValues(0))); SO_ENGINE_OUTPUT(normals, SoMFVec3f, setValues(0, norms.getNum(), norms.getValues(0))); int i = inds.getNum(); SO_ENGINE_OUTPUT(indexes, SoMFInt32, setNum(i)); SO_ENGINE_OUTPUT(indexes, SoMFInt32, setValues(0, i, inds.getValues(0))); SO_ENGINE_OUTPUT(numTriangles, SoSFFloat, setValue(triNum)); } ////////////////////////////////////////////////////////////////////////////// // // Description: // To avoid calculating points multiple times a "Marching Edges" // technique is used (this may as well be called something else by someone // else, I do not know. I did not research on this technique, just // thought of it one day while waiting in an airport.) Instead of // searching each cube of data one by one and accessing 8 different points // at once only 2 are looked at each time, and less assigments are made. // The two pulled out make up an edge on a cube. If the isoValue is // inbetween the two edge values then a point is computed for the // surface. To eliminate additional calls to the same data point the // order of operation is done in the following manner... // // V2 // | V1 // y| / // | / z // |/ // V0-------V3 // x // // V0 and V1 are checked, // V0 and V2 are checked, // V0 and V3 are checked // // Every time a point is added to the 'verts' array its index is stored in // the 'edgeIndex' array. This way when building the triangles with the // Marching Cubes tables the proper edge can point to the proper point. // ////////////////////////////////////////////////////////////////////////////// void MarchingCubes::calcPoints () { int i, j, k, ci, cj, ck, cv; int currCube, currEdge, currVert, lastVert; int offset[3], checkDims[3]; int cyOffset, czOffset; int32_t *cptr, *eptr; int v1, v2, index; float delta, percent; SbVec3f basePt, pt, norm; const SbVec3f *gptr; // the point blending: we are allocating the voxel structure // in order to process as locate the blended points as possible later // blend :: point *bTableP, *bPntP; bool nearCorner; bTableP = new blend :: point[ dims[0] * dims[1] * dims[2] ]; // add a point to the vertice list. A point is added when it is determined // that the isosurface goes through an edge. Sets the normal as well. // (accumulate if too close to the voxel corner except for the boundaries) // #define MC_ADD_POINT(comp, v1, v2) \ { \ pt = basePt; \ delta = values[v2] - values[v1]; \ percent = (delta == 0.0) ? 0.5 : (isoval - values[v1]) / delta; \ pt[comp] += voxel[comp] * percent; \ norm = gptr[v1] + (gptr[v2] - gptr[v1]) * percent; \ norm.normalize(); \ \ nearCorner = ( fabs( percent ) < bperc); \ if ( i && j && k && \ i < checkDims[0] && j < checkDims[1] && k < checkDims[2] && \ ( nearCorner || ( fabs( 1.0f - percent ) < bperc ) ) ) \ { \ bPntP = bTableP + ((nearCorner) ? v1 : v2); \ if ( bPntP->hits() ) \ { \ lastVert = bPntP->add( pt, norm ); \ } \ else \ { \ bPntP->set( currVert, pt, norm ); \ lastVert = currVert; \ ++ currVert; \ } \ } \ else \ { \ verts.set1Value(currVert, pt); \ norms.set1Value(currVert, norm); \ lastVert = currVert; \ ++ currVert; \ } \ eptr[currEdge] = lastVert; \ } // check an edge to see if v2 is > the isoval. This is called // when it is already known that the v1 is <= to the isoval #define MC_CHECK_EDGE_1(comp) \ { \ v2 = v1 + offset[comp]; \ if (values[v2] > isoval) { \ MC_ADD_POINT(comp, v1, v2); \ } \ } \ currEdge++; // check an edge to see if v2 is <= to the isoval. This is called when it // is already know that v1 is > than the isoval. #define MC_CHECK_EDGE_2(comp) \ { \ v2 = v1 + offset[comp]; \ if (values[v2] <= isoval) { \ MC_ADD_POINT(comp, v1, v2); \ } \ } \ currEdge++; // set offsets for each direction offset[0] = 1; offset[1] = dims[0]; offset[2] = dims[0]*dims[1]; // just a temporary variable to avoid checking (dims[x] - 1) against a // number, the less CPU cycles the better. checkDims[0] = dims[0] - 1; checkDims[1] = dims[1] - 1; checkDims[2] = dims[2] - 1; // set offsets for each direction when computing the index of surrounding // cubes cyOffset = cubeDims[0]; czOffset = cubeDims[0] * cubeDims[1]; // initialize all variables gptr = gradients.getValues(0); eptr = edgeIndex.startEditing(); cptr = cubes.startEditing(); memset(cptr, 0, sizeof(int32_t)*cubeDims[0]*cubeDims[1]*cubeDims[2]); for (i = 0; i < dims[0]*dims[1]*dims[2]; i++) eptr[i] = -1; currCube = 0; currVert = 0; currEdge = 0; v1 = 0; // for each vertice, check the edges it is connected to in the positive // x,y,z direction. for (k = 0, ck = 1; k < dims[2]; k++, ck++) { basePt[2] = min[2] + k*voxel[2]; for (j = 0, cj = 1; j < dims[1]; j++, cj++) { basePt[1] = min[1] + j*voxel[1]; for (i = 0, ci = 1; i < dims[0]; i++, ci++) { basePt[0] = min[0] + i*voxel[0]; // most of the time all edges in each direction will be checked if (k < checkDims[2] && j < checkDims[1] && i < checkDims[0]) { // if the value at data point 'v1' is <= isoval then the cubes that // contain the point as a vertice need to have the bits toggled. if (values[v1] <= isoval) { // compute index of the current cube, this is always the cube // where 'v1' is vertex 0 currCube = ck*czOffset + cj*cyOffset + ci; for (cv = 0; cv < 8; cv++) { index = currCube + cubeVertOffset[cv][2]*czOffset + cubeVertOffset[cv][1]*cyOffset + cubeVertOffset[cv][0]; cptr[index] |= (1 << cv); } MC_CHECK_EDGE_1(2); MC_CHECK_EDGE_1(1); MC_CHECK_EDGE_1(0); } else { MC_CHECK_EDGE_2(2); MC_CHECK_EDGE_2(1); MC_CHECK_EDGE_2(0); } } else { if (values[v1] <= isoval) { currCube = ck*czOffset + cj*cyOffset + ci; for (cv = 0; cv < 8; cv++) { index = currCube + cubeVertOffset[cv][2]*czOffset + cubeVertOffset[cv][1]*cyOffset + cubeVertOffset[cv][0]; cptr[index] |= (1 << cv); } if (k < checkDims[2]) MC_CHECK_EDGE_1(2); if (j < checkDims[1]) MC_CHECK_EDGE_1(1); if (i < checkDims[0]) MC_CHECK_EDGE_1(0); } else { if (k < checkDims[2]) MC_CHECK_EDGE_2(2); if (j < checkDims[1]) MC_CHECK_EDGE_2(1); if (i < checkDims[0]) MC_CHECK_EDGE_2(0); } } v1++; } } } bPntP = bTableP; for (k = 0; k < dims[2]; ++k) { for (j = 0; j < dims[1]; ++j) { for (i = 0; i < dims[0]; ++i, ++ bPntP) { if ( bPntP->hits() ) { verts.set1Value(bPntP->vertexIndex(), bPntP->location()); norms.set1Value(bPntP->vertexIndex(), bPntP->normal()); } } } } delete [] bTableP; } ////////////////////////////////////////////////////////////////////////////// // // Description: // Create the index triples which will define the triangles used to create // the surface. This is done by iterating through each cube in 'cubes' (not // including the buffer cubes) and using the value as an index into the // 'triTable' array. The 'triTable' array is used to determine which // edges the surface is intersecting. Based on this information the edge // index is computed and then the index of the point is found. // ////////////////////////////////////////////////////////////////////////////// void MarchingCubes::calcIndexes () { int c, e, i, j, k, t; int currIndex; int v1, v2, v3; const int32_t *cptr, *eptr; // these offsets are used to index into the 'edgeIndex' array. The edges // of a cube according to the Marching Cubes table 'triTable' are different // then how the edges are stored in the 'edgeIndex' array. The offset of // an edge is found by indexing into the array by the edge number (0-11). // Notice these offsets are a function of the dimensions of the data. const int32_t offset[12] = { 2, 4, 3*dims[0] + 2, 1, 3*dims[1]*dims[0] + 2, 3*dims[1]*dims[0] + 4, 3*dims[1]*dims[0] + 3*dims[0] + 2, 3*dims[1]*dims[0] + 1, 0, 3, 3*dims[0] + 3, 3*dims[0] }; // get pointers to the cubes and edges cptr = cubes.getValues(0); eptr = edgeIndex.getValues(0); e = currIndex = triNum = 0; // for each cube use the value as an index into 'triTable' and determine if // the cube should have any triangles in it. for (k = 1; k < cubeDims[2] - 1; k++) { for (j = 1; j < cubeDims[1] - 1; j++) { for (i = 1; i < cubeDims[0] - 1; i++) { // create the cube index c = cubeDims[0] * (k*cubeDims[1] + j) + i; // at most there will be 5 triangles in a cube, cycle through until a // -1 is found (signifying no more triangles) for (t = 0; t < 15; t += 3) { if (triTable[cptr[c]][t] < 0) break; // extract the indexes by... // - look up triangle combination in 'triTable' // - use 't' to index into 'triTable' and determine which edges // make up the triangle // - use edge value as an index into the 'offset' array // - add offset to the current edge (always the lower left edge // of the current cube pointing in the z-direction) // - index into the 'edgeIndex' array and get the index of the // point corresponding to that edge. v1 = eptr[offset[triTable[cptr[c]][t + 0]] + e]; v2 = eptr[offset[triTable[cptr[c]][t + 1]] + e]; v3 = eptr[offset[triTable[cptr[c]][t + 2]] + e]; // the blending of vertices may have edges pointing to the same // vertex, this would be a waste of a triangle so it is ignored. if ( v1 == v2 || v2 == v3 || v1 == v3 ) continue; // just incase something went wrong a check can be made //if (v1 < 0 || v2 < 0 || v3 < 0) { // cout << "no point " << v1 << " " << v2 << " " << v3 << endl; // continue; //} // set the values to form a triangle inds.set1Value(currIndex++, v1); inds.set1Value(currIndex++, v2); inds.set1Value(currIndex++, v3); inds.set1Value(currIndex++, -1); triNum++; } c++; e += 3; } e += 3; } e += 3*dims[0]; } } ////////////////////////////////////////////////////////////////////////////// // // Description: // Computes the gradients at each grid point. This is called only when // the data has been changed. This way when computing normals at each // triangle vertex, tri-linear interpolation will be used to efficiently // calculate the values. // ////////////////////////////////////////////////////////////////////////////// void MarchingCubes::calcGradients () { int i, j, k; int index, xOffset, yOffset, zOffset; SbVec3f *grad; // set the total number of gradient values gradients.setNum(dims[0] * dims[1] * dims[2]); // macro used to determine normal for a given component #define MC_COMPONENT(c, i, offset) \ if (i == 0) \ grad[index][c] = (values[index + offset] - values[index]) / voxel[c]; \ else if (i == dims[c] - 1) \ grad[index][c] = (values[index] - values[index - offset]) / voxel[c]; \ else \ grad[index][c] = (values[index + offset] - values[index - offset]) / \ (2.0 * voxel[c]); // get pointer to start editing gradients grad = gradients.startEditing(); index = 0; // set offsets for each direction xOffset = 1; yOffset = dims[0]; zOffset = dims[0]*dims[1]; for (k = 0; k < dims[2]; k++) { for (j = 0; j < dims[1]; j++) { for (i = 0; i < dims[0]; i++) { MC_COMPONENT(0, i, xOffset); MC_COMPONENT(1, j, yOffset); MC_COMPONENT(2, k, zOffset); grad[index].normalize(); index++; } } } gradients.finishEditing(); }