include/shape/mesher3d.hpp

Go to the documentation of this file.

/*****************************************************************************
 * M a t h e m a g i x
 *****************************************************************************
 * 2011-04-26
 * Bernard Mourrain
 *****************************************************************************
 *               Copyright (C) 2011 INRIA Sophia-Antipolis
 *****************************************************************************
 * Comments :
 ****************************************************************************/
# ifndef shape_mesher3d_hpp
# define shape_mesher3d_hpp

# include <shape/subdivision.hpp>
# include <shape/topology3d.hpp>
# include <shape/dualize3d.hpp>
# include <shape/polygonizer.hpp>

# define TMPL   template<class C,class V,class Cell>
# define TMPL1  template<class V>
# define SELF   mesher3d<C,V,Cell>
//====================================================================
namespace mmx { namespace shape {
//====================================================================
class mesher3d_def {};

template<class C,class V>
struct use< mesher3d_def,C,V> {

    template<class POINT, class CELL> static void
    edge_points(POINT* edges, const POINT* points, const CELL* cl, int cubeindex) {

        //Point* edges[12];
        //int cubeindex=marching_cube::index(values);

        if (marching_cube_edge_table[cubeindex] == 0) return ;

        if (marching_cube_edge_table[cubeindex] & 1)
            cl->point(edges[0], points[0],points[1]);
        if (marching_cube_edge_table[cubeindex] & 2)
            cl->point(edges[1], points[1],points[2]);
        if (marching_cube_edge_table[cubeindex] & 4)
            cl->point(edges[2],  points[2],points[3]);
        if (marching_cube_edge_table[cubeindex] & 8)
            cl->point(edges[3],  points[3],points[0]);
        if (marching_cube_edge_table[cubeindex] & 16)
            cl->point(edges[4],  points[4],points[5]);
        if (marching_cube_edge_table[cubeindex] & 32)
            cl->point(edges[5],  points[5],points[6]);
        if (marching_cube_edge_table[cubeindex] & 64)
            cl->point(edges[6],  points[6],points[7]);
        if (marching_cube_edge_table[cubeindex] & 128)
            cl->point(edges[7],  points[7],points[4]);
        if (marching_cube_edge_table[cubeindex] & 256)
            cl->point(edges[8],  points[0],points[4]);
        if (marching_cube_edge_table[cubeindex] & 512)
            cl->point(edges[9],  points[1],points[5]);
        if (marching_cube_edge_table[cubeindex] & 1024)
            cl->point(edges[10],  points[2],points[6]);
        if (marching_cube_edge_table[cubeindex] & 2048)
            cl->point(edges[11],  points[3],points[7]);
    }
};
//--------------------------------------------------------------------

template<class C,class V=default_env, class Cell=cell3d<C,V> >
struct mesher3d : public PROCESS_OF(V)
{
    typedef subdivision<C,V,Cell> Subdivisor;

    typedef Cell                                 Input;
    typedef typename use<mesh3d_def,C,V>::Mesh   Output;
    typedef use<topology3d_def,C,V>              Topology;

    typedef typename Output::Point      Point;
    typedef typename Output::Edge       Edge;

    mesher3d(double e1= 0.1, double e2=0.01);
    ~mesher3d(void) ;

    void set_input (Cell* cl) { m_input= cl; }

    void set_smoothness(double e) { m_smooth = e; }
    void set_precision (double e) { m_prec   = e; }

    double get_smoothness(void) { return m_smooth; }
    double get_precision (void) { return m_prec;   }


    Input*  get_input  (void)     { return m_input; }
    Output* get_output (void)     { return m_output; }

    void run(void);
    void clear(void);

    private:
    double          m_smooth ;
    double          m_prec ;

    Cell*           m_input;
    Output*         m_output;

};
//--------------------------------------------------------------------
TMPL SELF::mesher3d(double e1, double e2): m_smooth(e1), m_prec(e2)
{
    m_output = new Output;
}
//--------------------------------------------------------------------
TMPL SELF::~mesher3d(void) {
    delete m_output;
}
//--------------------------------------------------------------------
TMPL void SELF::run(void) {

    std::cout<< ">> Prec "<< m_smooth <<" " << m_prec<<"\n";

    Subdivisor* sbd = new Subdivisor(m_smooth,m_prec);
    sbd->set_input(this->get_input());
    sbd->run();

    std::cout<< ">> leaves    = "<< sbd->get_output()->m_leaves.size()<<"\n";

    typedef dualize<C,V,Cell>  Dualize;
    Dualize* dl = new Dualize;
    dl->set_input(sbd->get_output());
    dl->run();

    std::cout<< "\n>> Cells      = "<< sbd->get_output()->m_leaves.size()<<"\n";
    std::cout<< ">> Dual edges = "<< dl->get_output()->m_edges.size()/2<<"\n";
    std::cout<< ">> Dual faces = "<< dl->get_output()->m_faces.size()/4<<"\n";
    std::cout<< ">> Dual cells = "<< dl->get_output()->m_cells.size()/8<<"\n";


    Point  points[8];
    Point  edges[12];
    C      values[8];
    int    idx;

    for(unsigned i=0; i< dl->get_output()->m_edges.size(); i+=2) {

        Cell* c1= dl->get_output()->m_edges[i];
        Cell* c2= dl->get_output()->m_edges[i+1];

        //std::cout <<"center value "<< c1->center_value()<< " " <<c2->center_value()<< std::endl;

        //Edge* e= new Edge(new Point(c1->center()), new Point(c2->center()));
        //this->get_output()->insert_face(new Point(c1->center()), new Point(c2->center()), new Point(c1->center()));

        if (c1 != c2) {
            C xm=std::max(c1->xmin(),c2->xmin()), xM=std::min(c1->xmax(),c2->xmax()),
                    ym=std::max(c1->ymin(),c2->ymin()), yM=std::min(c1->ymax(),c2->ymax()),
                    zm=std::max(c1->zmin(),c2->zmin()), zM=std::min(c1->zmax(),c2->zmax());

            if (zm == zM) {
                points[0].x()= xm; points[0].y()=ym; points[0].z()=zm; values[0]= c1->value(xm,ym,zm);
                points[1].x()= xM; points[1].y()=ym; points[1].z()=zm; values[1]= c1->value(xM,ym,zm);
                points[2].x()= xM; points[2].y()=yM; points[2].z()=zm; values[2]= c1->value(xM,yM,zm);
                points[3].x()= xm; points[3].y()=yM; points[3].z()=zm; values[3]= c1->value(xm,yM,zm);

                //this->get_output()->insert_face(new Point(points[0]),new Point(points[1]),new Point(points[2]),new Point(points[3]));
                if(c1->is_active()) {
                    for(unsigned j=4; j<8; j++) {
                        points[j] = c1->center();
                        values[j] = c1->center_value();
                    }

                    idx = marching_cube::index(values);
                    //marching_cube::edge_points(edges, points, values, idx);
                    use<mesher3d_def,C,V>::edge_points(edges, points, c1, idx);
                    marching_cube::polygonize(*this->get_output(), edges, idx);
                }

                if (c2->is_active()) {
                    for(unsigned j=4; j<8; j++) {
                        points[j] = c2->center();
                        values[j] = c2->center_value();
                    }

                    idx = marching_cube::index(values);
                    //marching_cube::edge_points(edges, points, values, idx);
                    use<mesher3d_def,C,V>::edge_points(edges, points, c2, idx);
                    marching_cube::polygonize(*this->get_output(), edges, idx);
                }

            } else if (ym == yM) {
                points[0].x()= xm; points[0].y()=ym; points[0].z()=zm; values[0]= c1->value(xm,ym,zm);
                points[1].x()= xM; points[1].y()=ym; points[1].z()=zm; values[1]= c1->value(xM,ym,zm);
                points[2].x()= xM; points[2].y()=ym; points[2].z()=zM; values[2]= c1->value(xM,ym,zM);
                points[3].x()= xm; points[3].y()=ym; points[3].z()=zM; values[3]= c1->value(xm,ym,zM);

                //this->get_output()->insert_face(new Point(points[0]),new Point(points[1]),new Point(points[2]),new Point(points[3]));
                if(c1->is_active()) {
                    for(unsigned j=4; j<8; j++) {
                        points[j] = c1->center();
                        values[j] = c1->center_value();
                    }

                    idx = marching_cube::index(values);
                    //marching_cube::edge_points(edges, points, values, idx);
                    use<mesher3d_def,C,V>::edge_points(edges, points, c1, idx);
                    marching_cube::polygonize(*this->get_output(), edges, idx);
                    //                marching_cube::polygonize(*this->get_output(), points, values);
                }

                if(c2->is_active()) {
                    for(unsigned j=4; j<8; j++) {
                        points[j] = c2->center();
                        values[j] = c2->center_value();
                    }
                    idx = marching_cube::index(values);
                    //marching_cube::edge_points(edges, points, values, idx);
                    use<mesher3d_def,C,V>::edge_points(edges, points, c2, idx);
                    marching_cube::polygonize(*this->get_output(), edges, idx);
                    //                marching_cube::polygonize(*this->get_output(), points, values);
                }
            } else if (xm == xM ) {
                points[0].x()= xm; points[0].y()=ym; points[0].z()=zm; values[0]= c1->value(xm,ym,zm);
                points[1].x()= xm; points[1].y()=yM; points[1].z()=zm; values[1]= c1->value(xm,yM,zm);
                points[2].x()= xm; points[2].y()=yM; points[2].z()=zM; values[2]= c1->value(xm,yM,zM);
                points[3].x()= xm; points[3].y()=ym; points[3].z()=zM; values[3]= c1->value(xm,ym,zM);

                //this->get_output()->insert_face(new Point(points[0]),new Point(points[1]),new Point(points[2]),new Point(points[3]));
                if(c1->is_active()) {
                    for(unsigned j=4; j<8; j++) {
                        points[j] = c1->center();
                        values[j] = c1->center_value();
                    }

                    idx = marching_cube::index(values);
                    //marching_cube::edge_points(edges, points, values, idx);
                    use<mesher3d_def,C,V>::edge_points(edges, points, c1, idx);
                    marching_cube::polygonize(*this->get_output(), edges, idx);
                }
                //                marching_cube::polygonize(*this->get_output(), points, values);

                if(c2->is_active()) {
                    for(unsigned j=4; j<8; j++) {
                        points[j] = c2->center();
                        values[j] = c2->center_value();
                    }
                    idx = marching_cube::index(values);
                    //marching_cube::edge_points(edges, points, values, idx);
                    use<mesher3d_def,C,V>::edge_points(edges, points, c2, idx);
                    marching_cube::polygonize(*this->get_output(), edges, idx);
                }
                //                marching_cube::polygonize(*this->get_output(), points, values);
            }
        }
        //bool mc  = true;  bool ctr = false;
        //bool mc  = false; bool ctr = true;

        //        typedef polygonizer<C,V, Seq<Cell*>, Output >  Polygonizer;
        //    Polygonizer* plg = new Polygonizer;

        //        if(mc) {
        //            plg->set_input(&sbd->get_output()->m_leaves);
        //            plg->set_output(this->get_output());
        //            plg->run();
        //        }

        //        if(ctr) {
        //            plg->set_input(&dl->get_output()->m_cells);
        //            plg->set_output(this->get_output());
        //            plg->run();
        //        }

        // if (ctr) {
        //   foreach(Cell * cl, sbd->output()->m_leaves) {
        //      //marching_cube::centralise(*cl, *cl);
        //     Point* p= new Point(cl->center());
        //      //        std::cout<<"center "<<p->x()<<" "<<p->y()<<" "<<p->z() <<std::endl;
        //     cl->insert(p);
        //     foreach(Point* p, cl->m_points) this->output()->insert(p);

        //      for(unsigned i=0; i< dl->output()->m_faces.size();i+=4){
        //        Face* f1=new Face;
        //        if(dl->output()->m_faces[i]->m_points.size()>0)
        //         f1->insert(dl->output()->m_faces[i]->m_points[0]);
        //        if(dl->output()->m_faces[i+1]->m_points.size()>0)
        //         f1->insert(dl->output()->m_faces[i+1]->m_points[0]);
        //        if(dl->output()->m_faces[i+2]->m_points.size()>0)
        //         f1->insert(dl->output()->m_faces[i+2]->m_points[0]);
        //        if(f1->size()>2)
        //         this->output()->insert(f1);
        //        else
        //         std::cout<<"Face 1 with "<<f1->size()<<" element"<<std::endl;

        //        Face* f2=new Face;
        //        if(dl->output()->m_faces[i]->m_points.size()>0)
        //         f2->insert(dl->output()->m_faces[i]->m_points[0]);
        //        if(dl->output()->m_faces[i+2]->m_points.size()>0)
        //         f2->insert(dl->output()->m_faces[i+2]->m_points[0]);
        //        if(dl->output()->m_faces[i+3]->m_points.size()>0)
        //         f2->insert(dl->output()->m_faces[i+3]->m_points[0]);
        //        if(f2->size()>2)
        //         this->output()->insert(f2);
        //        else
        //         std::cout<<"Face 2 with "<<f2->size()<<" element"<<std::endl;
        //      }
        //   }
        //   std::cout<< "Triangles  = "<< this->output()->m_faces.size()<<"\n";
        // }// ctr

    }
}

//--------------------------------------------------------------------
TMPL void SELF::clear(void) {
    this->output()->clear();
}
//====================================================================
} ; // namespace shape
              } ; // namespace mmx
//====================================================================
# undef TMPL
# undef TMPL1
# undef Shape
# undef Graphic
# undef SELF
# endif

---