iDevGames Forums
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 Smile

Code:
template < class T >
struct Vec2Comparator_
{
    bool operator()( const Vec2<T> &a, const Vec2<T> &b ) const
    {
        if ( a.x != b.x ) return a.x < b.x;
        return a.y < b.y;              
    }
};      

typedef Vec2Comparator_<int> Vec2iComparator;


class IslandVoxelSpaceAdapter
{
    private:
    
        const std::map<Vec2i,Vertex*,Vec2iComparator> &_voxels;
        Vec2i _min, _max;

    public:
    
        IslandVoxelSpaceAdapter( const std::map<Vec2i,Vertex*,Vec2iComparator> &voxels ):
            _voxels(voxels),
            _min(INT_MAX, INT_MAX ),
            _max( INT_MIN, INT_MIN ),
        {
            for( std::map<Vec2i,Vertex*,Vec2iComparator>::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<Vec2i,Vertex*,Vec2iComparator>::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<Vec2i,Vertex*,Vec2iComparator> _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.