MarchingSquares - Printable Version +- iDevGames Forums (http://www.idevgames.com/forums) +-- Forum: Development Zone (/forum-3.html) +--- Forum: Game Programming Fundamentals (/forum-7.html) +--- Thread: MarchingSquares (/thread-8761.html) MarchingSquares - TomorrowPlusX - Mar 27, 2011 07:19 AM Hey everybody, I've been working on a 2d voxel based game with dynamic terrain, blah blah. In the interest of proof-of-concepting I used 3d marching cubes to compute the shape of the voxel space, which I basically projected to 2D, and then triangulated. It worked, but was a hotspot in my profiling. I did some googling to see if there were any 2D versions... but as it happens, "MarchingSquares" is a totally different algorithm. So I put pen to paper and worked out the edge tables and whatnot for a 2D variant of MarchingCubes -- much easier given there are only 16 intersection types instead of 256 -- and got it working yesterday. Since I've gotten so much help in this forum over the years I'm contributing it here, for anybody who may need it in the future. It's C++ ( I know, yuck, whatever ) and I thin kit's fairly clean. It does depend on you providing some kind of Vec2f and Vec2i vector type with basic math operators. I'm using libcinder but you can provide your own easily enough. First, here's my marching squares implementation: Code: ```    #pragma once     ///////////////////////////////////////////////////////////////////////     // Adapted from "Polygonising A Scalar Field" by Paul Bourke     // http://paulbourke.net/geometry/polygonise/     namespace marching_squares     {         using namespace ci;         const int EdgeTable[16] = {             0x0,     //0000,             0x9,     //1001,             0x3,     //0011             0xa,     //1010             0x6,     //0110,             0xf,     //1111,             0x5,     //0101             0xc,     //1100             0xc,     //1100             0x5,     //0101             0xf,     //1111,             0x6,     //0110,             0xa,     //1010             0x3,     //0011             0x9,     //1001,             0x0,     //0000         };         // note the sub array is length 5, with a -1 at the end to mark end of array         const int SegmentTable[16][5] = {             {-1,-1,-1,-1,-1},             {0,3,-1,-1,-1},             {1,0,-1,-1,-1},             {1,3,-1,-1,-1},             {2,1,-1,-1,-1},             {2,1,0,3,-1},             {2,0,-1,-1,-1},             {2,3,-1,-1,-1},             {3,2,-1,-1,-1},             {0,2,-1,-1,-1},             {1,0,3,2,-1},             {1,2,-1,-1,-1},             {3,1,-1,-1,-1},             {0,1,-1,-1,-1},             {3,0,-1,-1,-1},             {-1,-1,-1,-1,-1}         };         /**             @class Segment             Segments are submitted to user via segment callback param of march()         */         struct Segment         {             Vec2f a,b;             Segment( void ) {}                     Segment( const Vec2f &A, const Vec2f &B    ):a(A),b(B){}                 };         #pragma mark -         #pragma mark MarchingSquares         struct GridCell         {             Vec2f v[4];             float val[4];             GridCell(){}         };         /*             Linearly interpolate the position where an isosurface cuts             an edge between two vertices, each with their own scalar value         */         inline ci::Vec2f VertexInterp( float isolevel, const Vec2f &p1, const Vec2f &p2, float valp1, float valp2 )         {             float mu;             Vec2f p;             if ( std::abs(isolevel-valp1) < EPSILON )                 return p1;             if ( std::abs(isolevel-valp2) < EPSILON )                 return p2;             if ( std::abs(valp1-valp2) < EPSILON )                 return p1;             mu = (isolevel - valp1) / (valp2 - valp1);             p.x = p1.x + mu * (p2.x - p1.x);             p.y = p1.y + mu * (p2.y - p1.y);             return p;         }         template< class VOXELSTORE >         bool GetGridCell( int x, int y, const VOXELSTORE &voxels, GridCell &cell )         {             // store the value in the voxel array -- clockwise             cell.val[0] = voxels.valueAt( x  , y   );             cell.val[1] = voxels.valueAt( x+1, y   );             cell.val[2] = voxels.valueAt( x+1, y+1 );             cell.val[3] = voxels.valueAt( x  , y+1 );             if (cell.val[0] > 0.0f ||                 cell.val[1] > 0.0f ||                 cell.val[2] > 0.0f ||                 cell.val[3] > 0.0f )             {                 // store the location in the voxel array                 cell.v[0] = Vec2f( x  , y   );                 cell.v[1] = Vec2f( x+1, y   );                 cell.v[2] = Vec2f( x+1, y+1 );                 cell.v[3] = Vec2f( x  , y+1 );                 return true;             }             return false;         }         /*             Given a GridCell and isolevel, compute the Segments required to represent             the isosurface through the cell.             Return the number of segments, the array 'segments' will be populated             with up to 2 Segments. Returns the number of Segments computed.         */         inline int Polygonise( const GridCell &grid, float isolevel, Segment *segments )         {             //             //    Determine the index into the edge table which             //    tells us which vertices are inside of the surface             //             int squareIndex = 0;             if (grid.val[0] < isolevel) squareIndex |= 1;             if (grid.val[1] < isolevel) squareIndex |= 2;             if (grid.val[2] < isolevel) squareIndex |= 4;             if (grid.val[3] < isolevel) squareIndex |= 8;             //             //    Square is entirely in/out of the surface             //             if (EdgeTable[squareIndex] == 0)                 return 0;             //             //    Find the vertices where the surface intersects the cube             //             Vec2f vertlist[4];             for ( int i = 0; i < 4; i++ ) vertlist[i] = Vec2f(0,0);             if (EdgeTable[squareIndex] & 1)             {                 vertlist[0] = VertexInterp(isolevel,grid.v[0],grid.v[1],grid.val[0],grid.val[1]);             }             if (EdgeTable[squareIndex] & 2)             {                 vertlist[1] = VertexInterp(isolevel,grid.v[1],grid.v[2],grid.val[1],grid.val[2]);             }             if (EdgeTable[squareIndex] & 4)             {                 vertlist[2] = VertexInterp(isolevel,grid.v[2],grid.v[3],grid.val[2],grid.val[3]);             }             if (EdgeTable[squareIndex] & 8)             {                 vertlist[3] = VertexInterp(isolevel,grid.v[3],grid.v[0],grid.val[3],grid.val[0]);             }             //             //    Finally, generate line segments             //             int nSegments = 0;             for ( int i = 0; SegmentTable[squareIndex][i] != -1; i += 2 )             {                 segments[nSegments].a = vertlist[SegmentTable[squareIndex][i  ]];                 segments[nSegments].b = vertlist[SegmentTable[squareIndex][i+1]];                 nSegments++;             }             return nSegments;         }             #pragma mark -         #pragma mark Public API         /**             March a 2D voxel space, invoking the Segment callback on each generated Segment             VOXELSPACE is an object with the following interface:                 - Vec2i min() const;                 - Vec2i max() const;                 - float valueAt( int x, int y ) const;             SEGCALLBACK is an object with the following interface:                 - void operator()( int x, int y, const marching_squares::Segment &seg )         */         template< class VOXELSTORE, class SEGCALLBACK >         void march( const VOXELSTORE &voxels, SEGCALLBACK &sc, float isolevel = 0.5f )         {             Segment segments[2];             GridCell cell;                     Vec2i min = voxels.min(), max = voxels.max();             for ( int y = min.y; y < max.y; y++ )             {                 for ( int x = min.x; x < max.x; x++ )                 {                     if ( GetGridCell( x,y, voxels, cell ) )                     {                                                 for ( int s = 0, nSegments = Polygonise( cell, isolevel, segments ); s < nSegments; s++ )                         {                             sc( x,y, segments[s] );                         }                     }                 }             }             }     }``` Now, here's demo code that calls marching squares and generates an array of continous clockwise windings from the unconnected segments marching squares generates. Holes in your isosurface will be wound in the counter-clockwise. Note, this is demo code -- for dumb reasons my pretend voxel store is using a std::map instead of a 2d array so it's slow. The actual game uses faster lookups, but that's an irrelevant detail here Code: ```template < class T > struct Vec2Comparator_ {     bool operator()( const Vec2 &a, const Vec2 &b ) const     {         if ( a.x != b.x ) return a.x < b.x;         return a.y < b.y;                   } };       typedef Vec2Comparator_ Vec2iComparator; class IslandVoxelSpaceAdapter {     private:              const std::map &_voxels;         Vec2i _min, _max;     public:              IslandVoxelSpaceAdapter( const std::map &voxels ):             _voxels(voxels),             _min(INT_MAX, INT_MAX ),             _max( INT_MIN, INT_MIN ),         {             for( std::map::const_iterator it(voxels.begin()),end(voxels.end()); it != end; ++it )             {                 _min = ::min( _min, it->second->ordinalPosition );                 _max = ::max( _max, it->second->ordinalPosition );             }         }         ~IslandVoxelSpaceAdapter(){}         /**             To march edges correctly, we need to outset by one         */         Vec2i min() const { return Vec2i( _min.x - 1, _min.y - 1 ); }         /**             To march edges correctly, we need to outset by one         */         Vec2i max() const { return Vec2i( _max.x + 1, _max.y + 1 ); }                  /**             Get the occupation value at a voxel coord. We return 0 for values outside the island.         */         float valueAt( int x, int y ) const         {             if ( x < _min.x || x > _max.x || y < _min.y || y > _max.y ) return 0.0f;             std::map::const_iterator pos = _voxels.find( Vec2i( x, y ));             if ( pos != _voxels.end() )             {                 return pos->second->occupation;             }                          return 0.0f;         } }; // minimum float delta we care about const float V_EPSILON = 1.0f / 256.0f; const float V_SCALE = 256.0f; struct PerimeterGenerator {     public:              enum winding {             CLOCKWISE,             COUNTER_CLOCKWISE         };     private:              typedef std::pair< Vec2i, Vec2i > Edge;             std::map< Vec2i, Edge, Vec2iComparator > _edgesByFirstVertex;         winding _winding;              public:              PerimeterGenerator( winding w ):             _winding( w )         {}                  inline void operator()( int x, int y, const ms::Segment &seg )         {             switch( _winding )             {                 case CLOCKWISE:                 {                     Edge e( scaleUp(seg.a), scaleUp(seg.b));                     if ( e.first != e.second ) _edgesByFirstVertex[e.first] = e;                     break;                 }                                  case COUNTER_CLOCKWISE:                 {                     Edge e( scaleUp(seg.b), scaleUp(seg.a));                     if ( e.first != e.second ) _edgesByFirstVertex[e.first] = e;                     break;                 }             }         }                  void generate( std::vector< Vec2fVec > &perimeters )         {               perimeters.clear();                          while( !_edgesByFirstVertex.empty() )             {                 std::map< Vec2i, Edge, Vec2iComparator >::iterator                     begin = _edgesByFirstVertex.begin(),                     end = _edgesByFirstVertex.end(),                     it = begin;                                      int count = _edgesByFirstVertex.size();                                 perimeters.push_back( Vec2fVec() );                 do                 {                     perimeters.back().push_back( scaleDown(it->first));                     _edgesByFirstVertex.erase(it);                     it = _edgesByFirstVertex.find(it->second.second);                     count--;                                  } while( it != begin && it != end && count > 0 );             }         }              private:              inline Vec2i scaleUp( const Vec2f &v ) const         {             return Vec2i( lrintf( V_SCALE * v.x ), lrintf( V_SCALE * v.y ));         }                  inline Vec2f scaleDown( const Vec2i &v ) const         {             return Vec2f( float(v.x) / V_SCALE, float(v.y) / V_SCALE );         }          };``` Finally, here's code that puts it all together Code: ```//std::map _verticesByOrdinalPosition = ... IslandVoxelSpaceAdapter voxels( _verticesByOrdinalPosition ); PerimeterGenerator pgen( PerimeterGenerator::CLOCKWISE ); // march IslandVoxelSpaceAdapter and collect segments in PerimeterGenerator ms::march( voxels, pgen ); // stitch segments into an array of contiguous loops std::vector< Vec2fVec > perimeters; pgen.generate( perimeters ); for ( std::vector< Vec2fVec >::iterator it( perimeters.begin()), end( perimeters.end()); it != end; ++it ) {     //triangulate, or whatever     }``` EDIT Memory-smasher bugfix -- SegmentTable's subarrays need to be length 5 with a terminator -1 value. SADFACE.