include/shape/marching_cube.hpp

Go to the documentation of this file.

# ifndef shape_marching_cube_hpp
# define shape_marching_cube_hpp
/*



*/
# define ULONG unsigned long 
# define TMPL template<class K>

namespace mmx { namespace shape {

extern int marching_cube_edge_table[256];
extern int marching_cube_tri_table[256][16] ;

template <class C, class POINT>
void mc_interpolate(POINT& p, const POINT& p1, const POINT& p2, C valp1, const C& valp2, const C& iso=0)
{
    if (abs(iso-valp1) < 0.00001) p=p1;
    if (abs(iso-valp2) < 0.00001) p=p2;
    if (abs(valp1-valp2) < 0.00001) p=p1;

    C mu = (iso - valp1) / (valp2 - valp1);

    p.x() = p1.x() + mu * (p2.x() - p1.x());
    p.y() = p1.y() + mu * (p2.y() - p1.y());
    p.z() = p1.z() + mu * (p2.z() - p1.z());
}
//----------------------------------------------------------------
struct marching_square {

public:

    marching_square(void) {};
    ~marching_square(void){};

    template<class Mesh, class POINTS, class VALUES, class C>
    static bool polygonize(Mesh& res, const POINTS& points, const VALUES& values, const C& lvl);

};

template<class Mesh, class POINTS, class VALUES, class C> bool
marching_square::polygonize(Mesh& res, const POINTS& P, const VALUES& val, const C& lvl) {

    typedef typename Mesh::Point       Point;
    typedef typename Mesh::Edge        Edge;

    Point pt[4];

    int nrPoint=0;

    if(val[0]>lvl)
    {
        if(val[1]<=lvl)
        {
            mc_interpolate(pt[nrPoint++], P[0], P[1], val[0], val[1], lvl);
        }

        if(val[3]<=lvl)
        {
            mc_interpolate(pt[nrPoint++], P[0], P[3], val[0], val[3], lvl);
        }
    }

    if(val[1]>lvl)
    {
        if(val[0]<=lvl)
        {
            mc_interpolate(pt[nrPoint++],  P[1], P[0], val[1], val[0], lvl);
        }

        if(val[2]<=lvl)
        {
            mc_interpolate(pt[nrPoint++],  P[1], P[2], val[1], val[2], lvl);
        }
    }


    if(val[2]>lvl)
    {
        if(val[1]<=lvl)
        {
            mc_interpolate(pt[nrPoint++],  P[2], P[1], val[2], val[1], lvl);
        }

        if(val[3]<=lvl)
        {
            mc_interpolate(pt[nrPoint++],  P[2], P[3], val[2], val[3],lvl);
        }
    }

    if(val[3]>lvl)
    {
        if(val[2]<=lvl)
        {
            mc_interpolate(pt[nrPoint++],  P[3], P[2], val[3], val[2], lvl);
        }

        if(val[0]<=lvl)
        {
            mc_interpolate(pt[nrPoint++],  P[3], P[0], val[3], val[0], lvl);
        }
    }

    Point* p1, *p2;
    for (int i=0;i < nrPoint;i+=2) {
        p1 = new Point(pt[i]);
        res.insert_vertex(p1);
        p1->set_index(res.vertices().size()-1);
        p2 = new Point(pt[i+1]);
        res.insert_vertex(p2);
        p2->set_index(res.vertices().size()-1);
        Edge* e= new Edge(p1, p2);
        res.insert_edge(e);
    }
    return true;
}

//************************************************************************
// Marching cube for smooth surface
//************************************************************************
struct marching_cube {

public:

    marching_cube(void) {};
    ~marching_cube(void){};

    template<class C>
    static int index (C* t, C iso=0);

    template<class POINT, class VALUES>
    static bool edge_points(POINT* edges, const POINT* points, const VALUES& values, int idx);

    template<class Mesh, class Cell>
    static bool polygonise(Mesh& res, const Cell& grid);

    template<class Mesh, class POINTS>
    static bool polygonize(Mesh& res, const POINTS& edges, int cubeindex);

    template<class Mesh, class POINTS, class VALUES>
    static bool polygonize(Mesh& res, const POINTS& points, const VALUES& values);

    template<class Mesh, class Cell>
    static bool centralise(Mesh& res, const Cell& grid);

};

//--------------------------------------------------------------------
template<class C>
int marching_cube::index (C* t, C iso) {
    int CubeIndex=0;
    int s=1;
    for (unsigned i=0; i< 8; i++) {
        if (t[i] < iso) CubeIndex |= s;
        s*=2;
    }
    return CubeIndex;
}
//************************************************************************
// Polygonise the cell in which isosurface pass through
//************************************************************************
template<class Mesh, class Cell> bool 
marching_cube::polygonise(Mesh& res, const Cell& bcell) {

    typedef typename Mesh::Point       Point;
    typedef typename Mesh::Face        Face;

    int CubeIndex=bcell.mc_index();

    if (marching_cube_edge_table[CubeIndex] == 0) return false;

    Point* CELL[12];
    if(bcell.edge_point(CELL, marching_cube_edge_table[CubeIndex])) {
        for (int i=0;marching_cube_tri_table[CubeIndex][i]!=-1;i+=3) {
            Face* F= new Face;
            for(int j=0;j<3;j++) {
                Point* p=CELL[marching_cube_tri_table[CubeIndex][i+j]];
                res.insert(p);
                F->insert(p);
            }
            res.insert(F);
        }
        return true;
    } else
        return false;
}

//************************************************************************
// Polygonise the cell in which isosurface pass through
//************************************************************************
template<class POINT, class VALUES> bool
marching_cube::edge_points(POINT* edges, const POINT* points, const VALUES& values, int cubeindex) {

    //Point* edges[12];

    //int cubeindex=marching_cube::index(values);
    //  std::cout<< cubeindex <<std::endl;
    if (marching_cube_edge_table[cubeindex] == 0) return false;

    if (marching_cube_edge_table[cubeindex] & 1)
        mc_interpolate(edges[0], points[0],points[1],values[0],values[1]);
    if (marching_cube_edge_table[cubeindex] & 2)
        mc_interpolate(edges[1], points[1],points[2],values[1],values[2]);
    if (marching_cube_edge_table[cubeindex] & 4)
        mc_interpolate(edges[2],  points[2],points[3],values[2],values[3]);
    if (marching_cube_edge_table[cubeindex] & 8)
        mc_interpolate(edges[3],  points[3],points[0],values[3],values[0]);
    if (marching_cube_edge_table[cubeindex] & 16)
        mc_interpolate(edges[4],  points[4],points[5],values[4],values[5]);
    if (marching_cube_edge_table[cubeindex] & 32)
        mc_interpolate(edges[5],  points[5],points[6],values[5],values[6]);
    if (marching_cube_edge_table[cubeindex] & 64)
        mc_interpolate(edges[6],  points[6],points[7],values[6],values[7]);
    if (marching_cube_edge_table[cubeindex] & 128)
        mc_interpolate(edges[7],  points[7],points[4],values[7],values[4]);
    if (marching_cube_edge_table[cubeindex] & 256)
        mc_interpolate(edges[8],  points[0],points[4],values[0],values[4]);
    if (marching_cube_edge_table[cubeindex] & 512)
        mc_interpolate(edges[9],  points[1],points[5],values[1],values[5]);
    if (marching_cube_edge_table[cubeindex] & 1024)
        mc_interpolate(edges[10],  points[2],points[6],values[2],values[6]);
    if (marching_cube_edge_table[cubeindex] & 2048)
        mc_interpolate(edges[11],  points[3],points[7],values[3],values[7]);
    return true;
}

//************************************************************************
// Polygonise the cell in which isosurface pass through
//************************************************************************
template<class Mesh, class POINTS> bool
marching_cube::polygonize(Mesh& res, const POINTS& edges, int cubeindex) {

    typedef typename Mesh::Point Point;
    for (int i=0;marching_cube_tri_table[cubeindex][i]!=-1;i+=3) {
//        Face* F= new Face;
//        for(int j=0;j<3;j++) {
//            Point* p=edges[marching_cube_tri_table[cubeindex][i+j]];
//            res.insert_vertex(edges[marching_cube_tri_table[cubeindex][i+j]]);
//        }
        res.insert_face(new Point(edges[marching_cube_tri_table[cubeindex][i]]),
                        new Point(edges[marching_cube_tri_table[cubeindex][i+1]]),
                        new Point(edges[marching_cube_tri_table[cubeindex][i+2]]));
    }
    return true;
}

template<class Mesh, class POINTS, class VALUES> bool 
marching_cube::polygonize(Mesh& res, const POINTS& points, const VALUES& values) {

    typedef typename Mesh::Point       Point;
    typedef typename Mesh::Face        Face;

    int cubeindex=marching_cube::index(values);
    Point edges[12];
    marching_cube::edge_points(edges, points, values, cubeindex);
    marching_cube::polygonize(res, edges, cubeindex);
    return true;
}

//--------------------------------------------------------------------
// Compute a "center" for the bcell
//--------------------------------------------------------------------
template<class Mesh, class Cell> bool 
marching_cube::centralise(Mesh& res, const Cell& bcell) {

    typedef typename Mesh::Point       Point;
    typedef typename Mesh::Face        Face;

    int CubeIndex=bcell.mc_index();
    //  std::cout<<CubeIndex<<std::endl;
    if (marching_cube_edge_table[CubeIndex] == 0) return false;

    Point* CELL[12];

    if(bcell.edge_point(CELL, marching_cube_edge_table[CubeIndex])) {
        Point* c=new Point(0,0,0);
        unsigned ntr=0;
        for (int i=0;marching_cube_tri_table[CubeIndex][i]!=-1;i+=3)
            ntr++;
        ntr/=2;

        for(int j=0;j<3;j++) {
            Point* p=CELL[marching_cube_tri_table[CubeIndex][3*ntr+j]];
            c->x()+= p->x();
            c->y()+= p->y();
            c->z()+= p->z();
        }
        c->x()/=3.;
        c->y()/=3.;
        c->z()/=3.;
        //   c->x()=(bcell.xmin()+bcell.xmax())/2;
        //   c->y()=(bcell.ymin()+bcell.ymax())/2;
        //   c->z()=(bcell.zmin()+bcell.zmax())/2;
        //std::cout<<*c<<std::endl;
        res.insert(c);
        return true;
    } else
        return false;
}

//====================================================================
} //namespace shape
              }//namespace mmx
# undef TMPL
# undef ULONG
# endif //shape_marching_cube_hpp

---