diff --git a/Polyhedron/demo/Polyhedron/include/CGAL/enriched_polyhedron.h b/Polyhedron/demo/Polyhedron/include/CGAL/enriched_polyhedron.h new file mode 100644 index 00000000000..6c2eec59103 --- /dev/null +++ b/Polyhedron/demo/Polyhedron/include/CGAL/enriched_polyhedron.h @@ -0,0 +1,923 @@ +// -*- tab-width: 2 -*- + +/////////////////////////////////////////////////////////////////////////// +// // +// Class: Enriched_polyhedron // +// // +/////////////////////////////////////////////////////////////////////////// + +#ifndef _POLYGON_MESH_ +#define _POLYGON_MESH_ + +// CGAL stuff +#include +#include +#include +#include + +#include + +// a refined facet with a normal and a tag +template +class Enriched_facet : public CGAL::HalfedgeDS_face_base +{ + // tag + int m_tag; + + // normal + Norm m_normal; + +public: + + // life cycle + // no constructors to repeat, since only + // default constructor mandatory + + Enriched_facet() + { + } + + // tag + const int& tag() const { return m_tag; } + int& tag() { return m_tag; } + void tag(const int& i) { m_tag = i; } + + // normal + typedef Norm Normal_3; + Normal_3& normal() { return m_normal; } + const Normal_3& normal() const { return m_normal; } +}; + +// a refined halfedge with a general tag and +// a binary tag to indicate wether it belongs +// to the control mesh or not +template +class Enriched_halfedge : public CGAL::HalfedgeDS_halfedge_base +{ +private: + + // general purpose tag + int m_tag; + + // option for edge superimposing + bool m_control_edge; + bool m_sharp; + +public: + + // life cycle + Enriched_halfedge() + { + m_control_edge = true; + m_sharp = false; + } + + // tag + const int& tag() const { return m_tag; } + int& tag() { return m_tag; } + void tag(const int& t) { m_tag = t; } + + // control edge + bool& control_edge() { return m_control_edge; } + const bool& control_edge() const { return m_control_edge; } + + // sharp + bool& sharp() { return m_sharp; } + const bool& sharp() const { return m_sharp; } +}; + + + +// a refined vertex with a normal and a tag +template +class Enriched_vertex : public CGAL::HalfedgeDS_vertex_base +{ + // tag + int m_tag; + + // normal + Norm m_normal; + +public: + // life cycle + Enriched_vertex() {} + // repeat mandatory constructors + Enriched_vertex(const P& pt) + : CGAL::HalfedgeDS_vertex_base(pt) + { + } + + // normal + typedef Norm Normal_3; + Normal_3& normal() { return m_normal; } + const Normal_3& normal() const { return m_normal; } + + // tag + int& tag() { return m_tag; } + const int& tag() const { return m_tag; } + void tag(const int& t) { m_tag = t; } +}; + +// A redefined items class for the Polyhedron_3 +// with a refined vertex class that contains a +// member for the normal vector and a refined +// facet with a normal vector instead of the +// plane equation (this is an alternative +// solution instead of using +// Polyhedron_traits_with_normals_3). + +struct Enriched_items : public CGAL::Polyhedron_items_3 +{ + // wrap vertex + template struct Vertex_wrapper + { + typedef typename Traits::Point_3 Point; + typedef typename Traits::Vector_3 Normal; + typedef Enriched_vertex Vertex; + }; + + // wrap face + template struct Face_wrapper + { + typedef typename Traits::Point_3 Point; + typedef typename Traits::Vector_3 Normal; + typedef Enriched_facet Face; + }; + + // wrap halfedge + template struct Halfedge_wrapper + { + typedef typename Traits::Vector_3 Normal; + typedef Enriched_halfedge Halfedge; + }; +}; + +static const double PI = 3.1415926535897932384626; + +// compute facet normal +struct Facet_normal // (functor) +{ + template void operator()(Facet& f) + { + typename Facet::Normal_3 sum = CGAL::NULL_VECTOR; + typename Facet::Halfedge_around_facet_circulator h = f.facet_begin(); + do + { + typename Facet::Normal_3 normal = CGAL::cross_product(h->next()->vertex()->point() - h->vertex()->point(), h->next()->next()->vertex()->point() - h->next()->vertex()->point()); + double sqnorm = normal * normal; + if (sqnorm != 0) + normal = normal / (float)std::sqrt(sqnorm); + sum = sum + normal; + } while (++h != f.facet_begin()); + float sqnorm = sum * sum; + if (sqnorm != 0.0) + f.normal() = sum / std::sqrt(sqnorm); + else + { + f.normal() = CGAL::NULL_VECTOR; + // TRACE("degenerate face\n"); + } + } +}; + +// compute vertex normal +struct Vertex_normal // (functor) +{ + template void operator()(Vertex& v) + { + typename Vertex::Normal_3 normal = CGAL::NULL_VECTOR; + typename Vertex::Halfedge_around_vertex_const_circulator pHalfedge = + v.vertex_begin(); + typename Vertex::Halfedge_around_vertex_const_circulator begin = + pHalfedge; + CGAL_For_all(pHalfedge,begin) + if(!pHalfedge->is_border()) + normal = normal + pHalfedge->facet()->normal(); + float sqnorm = normal * normal; + if (sqnorm != 0.0f) + v.normal() = normal / (float)std::sqrt(sqnorm); + else + v.normal() = CGAL::NULL_VECTOR; + } +}; + +//********************************************************* +template +#endif + class HDS = CGAL::HalfedgeDS_default > +class Enriched_polyhedron : public CGAL::Polyhedron_3 +{ +public : + typedef typename kernel::FT FT; + typedef typename kernel::Point_3 Point; + typedef typename kernel::Vector_3 Vector; + typedef typename kernel::Iso_cuboid_3 Iso_cuboid; + typedef CGAL::Polyhedron_3 Base; + + typedef typename Base::Vertex_handle Vertex_handle; + typedef typename Base::Vertex_iterator Vertex_iterator; + typedef typename Base::Halfedge_handle Halfedge_handle; + typedef typename Base::Halfedge_iterator Halfedge_iterator; + typedef typename Base::Halfedge_around_facet_circulator Halfedge_around_facet_circulator; + typedef typename Base::Halfedge_around_vertex_circulator Halfedge_around_vertex_circulator; + typedef typename Base::Edge_iterator Edge_iterator; + typedef typename Base::Facet_iterator Facet_iterator; + typedef typename Base::Facet_handle Facet_handle; + typedef typename Base::Facet Facet; + + using Base::vertices_begin; + using Base::vertices_end; + using Base::edges_begin; + using Base::edges_end; + using Base::halfedges_begin; + using Base::halfedges_end; + using Base::facets_begin; + using Base::facets_end; + using Base::size_of_halfedges; + using Base::size_of_facets; + using Base::size_of_vertices; +private : + + // bounding box + Iso_cuboid m_bbox; + + // type + bool m_pure_quad; + bool m_pure_triangle; + +public : + enum Vertex_type { SMOOTH, // 0 sharp edge + DART, // 1 + CREASE_REGULAR, // 2 - with two non-sharp edges on each side + CREASE_IRREGULAR, // 2 + CORNER }; // 3 and more + +public : + + // life cycle + Enriched_polyhedron() + { + m_pure_quad = false; + m_pure_triangle = false; + } + virtual ~Enriched_polyhedron() + { + } + + // type + bool is_pure_triangle() { return m_pure_triangle; } + bool is_pure_quad() { return m_pure_quad; } + + // normals (per facet, then per vertex) + void compute_normals_per_facet() + { + std::for_each(facets_begin(),facets_end(),::Facet_normal()); + } + void compute_normals_per_vertex() + { + std::for_each(vertices_begin(),vertices_end(),::Vertex_normal()); + } + void compute_normals() + { + compute_normals_per_facet(); + compute_normals_per_vertex(); + } + + // bounding box + Iso_cuboid& bbox() { return m_bbox; } + const Iso_cuboid bbox() const { return m_bbox; } + + // compute bounding box + void compute_bounding_box() + { + CGAL_assertion(size_of_vertices() != 0); + + FT xmin,xmax,ymin,ymax,zmin,zmax; + Vertex_iterator pVertex = vertices_begin(); + xmin = xmax = pVertex->point().x(); + ymin = ymax = pVertex->point().y(); + zmin = zmax = pVertex->point().z(); + for(; + pVertex != vertices_end(); + pVertex++) + { + const Point& p = pVertex->point(); + + xmin = std::min(xmin,p.x()); + ymin = std::min(ymin,p.y()); + zmin = std::min(zmin,p.z()); + + xmax = std::max(xmax,p.x()); + ymax = std::max(ymax,p.y()); + zmax = std::max(zmax,p.z()); + } + m_bbox = Iso_cuboid(xmin,ymin,zmin, + xmax,ymax,zmax); + } + + // bounding box + FT xmin() { return m_bbox.xmin(); } + FT xmax() { return m_bbox.xmax(); } + FT ymin() { return m_bbox.ymin(); } + FT ymax() { return m_bbox.ymax(); } + FT zmin() { return m_bbox.zmin(); } + FT zmax() { return m_bbox.zmax(); } + + Point center() + { + FT cx = (FT)0.5 * (xmin() + xmax()); + FT cy = (FT)0.5 * (ymin() + ymax()); + FT cz = (FT)0.5 * (zmin() + zmax()); + return Point(cx,cy,cz); + } + + FT size() + { + FT dx = xmax() - xmin(); + FT dy = ymax() - ymin(); + FT dz = zmax() - zmin(); + return std::max(dx,std::max(dy,dz)); + } + + // degree of a face + static unsigned int degree(Facet_handle pFace) + { + return CGAL::circulator_size(pFace->facet_begin()); + } + + // valence of a vertex + static unsigned int valence(Vertex_handle pVertex) + { + return CGAL::circulator_size(pVertex->vertex_begin()); + } + + // check wether a vertex is on a boundary or not + static bool is_border(Vertex_handle pVertex) + { + Halfedge_around_vertex_circulator pHalfEdge = pVertex->vertex_begin(); + if(pHalfEdge == NULL) // isolated vertex + return true; + Halfedge_around_vertex_circulator d = pHalfEdge; + CGAL_For_all(pHalfEdge,d) + if(pHalfEdge->is_border()) + return true; + return false; + } + + // get any border halfedge attached to a vertex + Halfedge_handle get_border_halfedge(Vertex_handle pVertex) + { + Halfedge_around_vertex_circulator pHalfEdge = pVertex->vertex_begin(); + Halfedge_around_vertex_circulator d = pHalfEdge; + CGAL_For_all(pHalfEdge,d) + if(pHalfEdge->is_border()) + return pHalfEdge; + return NULL; + } + + // tag all vertices + void tag_vertices(const int tag) + { + for(Vertex_iterator vit = vertices_begin(), + end = vertices_end(); + vit != end; + ++vit) + vit->tag(tag); + } + + // tag all halfedges + void tag_halfedges(const int tag) + { + for(Halfedge_iterator pHalfedge = halfedges_begin(); + pHalfedge != halfedges_end(); + pHalfedge++) + pHalfedge->tag(tag); + } + + // tag all facets + void tag_facets(const int tag) + { + for(Facet_iterator pFacet = facets_begin(); + pFacet != facets_end(); + pFacet++) + pFacet->tag() = tag; + } + + // set index for all vertices + void set_index_vertices() + { + int index = 0; + for(Vertex_iterator pVertex = vertices_begin(); + pVertex != vertices_end(); + pVertex++) + pVertex->tag(index++); + } + + // set index for all edges + void set_index_edges() + { + int index = 0; + for(Edge_iterator he = edges_begin(); + he != edges_end(); + he++) + { + he->tag(index); + he->opposite()->tag(index); + index++; + } + } + + // is pure degree ? + bool is_pure_degree(unsigned int d) + { + for(Facet_iterator pFace = facets_begin(); + pFace != facets_end(); + pFace++) + if(degree(pFace) != d) + return false; + return true; + } + + // compute type + void compute_type() + { + m_pure_quad = is_pure_degree(4); + m_pure_triangle = is_pure_degree(3); + } + + // compute facet center + void compute_facet_center(Facet_handle pFace, + Point& center) + { + Halfedge_around_facet_circulator pHalfEdge = pFace->facet_begin(); + Halfedge_around_facet_circulator end = pHalfEdge; + Vector vec(0.0,0.0,0.0); + int degree = 0; + CGAL_For_all(pHalfEdge,end) + { + vec = vec + (pHalfEdge->vertex()->point()-CGAL::ORIGIN); + degree++; + } + center = CGAL::ORIGIN + (vec/FT(degree)); + } + + + void gl_draw_direct_triangles(bool smooth_shading, + bool use_normals, + bool inverse_normals = false) + { + // draw triangles + ::glBegin(GL_TRIANGLES); + Facet_iterator pFacet = facets_begin(); + for(;pFacet != facets_end();pFacet++) + gl_draw_facet(pFacet,smooth_shading,use_normals,inverse_normals); + ::glEnd(); // end polygon assembly + } + + + void gl_draw_direct(bool smooth_shading, + bool use_normals, + bool inverse_normals = false) + { + // draw polygons + Facet_iterator pFacet = facets_begin(); + for(;pFacet != facets_end();pFacet++) + { + // begin polygon assembly + ::glBegin(GL_POLYGON); + gl_draw_facet(pFacet,smooth_shading,use_normals,inverse_normals); + ::glEnd(); // end polygon assembly + } + } + + void gl_draw_facet(Facet_handle pFacet, + bool smooth_shading, + bool use_normals, + bool inverse_normals = false) + { + // one normal per face + if(use_normals && !smooth_shading) + { + const typename Facet::Normal_3& normal = pFacet->normal(); + if(inverse_normals) + ::glNormal3f(-normal[0],-normal[1],-normal[2]); + else + ::glNormal3f(normal[0],normal[1],normal[2]); + } + + // revolve around current face to get vertices + Halfedge_around_facet_circulator pHalfedge = pFacet->facet_begin(); + do + { + // one normal per vertex + if(use_normals && smooth_shading) + { + const typename Facet::Normal_3& normal = pHalfedge->vertex()->normal(); + if(inverse_normals) + ::glNormal3f(-normal[0],-normal[1],-normal[2]); + else + ::glNormal3f(normal[0],normal[1],normal[2]); } + + // polygon assembly is performed per vertex + const Point& point = pHalfedge->vertex()->point(); + ::glVertex3d(point[0],point[1],point[2]); + } + while(++pHalfedge != pFacet->facet_begin()); + } + + // superimpose edges + void superimpose_edges(bool skip_ordinary_edges = true, + bool skip_control_edges = false) + { + ::glBegin(GL_LINES); + for(Edge_iterator h = edges_begin(); + h != edges_end(); + h++) + { + if(h->sharp()) + continue; + + // ignore this edges + if(skip_ordinary_edges && !h->control_edge()) + continue; + + // ignore control edges + if(skip_control_edges && h->control_edge()) + continue; + + // assembly and draw line segment + const Point& p1 = h->prev()->vertex()->point(); + const Point& p2 = h->vertex()->point(); + ::glVertex3f(p1[0],p1[1],p1[2]); + ::glVertex3f(p2[0],p2[1],p2[2]); + } + ::glEnd(); + } + + bool is_sharp(Halfedge_handle he, + const double angle_sharp) + { + Facet_handle f1 = he->facet(); + Facet_handle f2 = he->opposite()->facet(); + if(f1 == Facet_handle() || f2 == Facet_handle() ) + return false; + const Vector& n1 = f1->normal(); + const Vector& n2 = f2->normal(); + if(angle_deg(n1,n2) > angle_sharp) + return true; + else + return false; + } + + Vertex_type type(Vertex_handle v) + { + unsigned int nb = nb_sharp_edges(v); + switch(nb) + { + case 0: + return SMOOTH; + case 1: + return DART; + case 2: // crease vertex - may be regular or not + return crease_type(v); + default: // 3 and more + return CORNER; + } + } + + // regular crease vertex must have valence 6, with + // exactly two smooth edges on each side. + Vertex_type crease_type(Vertex_handle v) + { + if(valence(v) != 6) + return CREASE_IRREGULAR; + + // valence = 6 - let us check regularity + + // pick first sharp edge + Halfedge_around_vertex_circulator he = v->vertex_begin(); + Halfedge_around_vertex_circulator end = he; + CGAL_For_all(he,end) + if(he->sharp()) + break; + + // next two must be smooth + for(int i=0;i<2;i++) + if(++he->sharp()) + return CREASE_IRREGULAR; + + // next one must be sharp + if(!++he->sharp()) + return CREASE_IRREGULAR; + + // next two must be smooth + for(int i=0;i<2;i++) + if(++he->sharp()) + return CREASE_IRREGULAR; + + return CREASE_REGULAR; + } + + // return true if succeds + void incident_points_on_crease(Vertex_handle v, + Point& a, + Point& b) + { +#ifdef _DEBUG + Vertex_type vertex_type = type(v,angle_sharp); + ASSERT(vertex_type == CREASE_IRREGULAR || vertex_type == CREASE_REGULAR); +#endif // _DEBUG + + // pick first sharp edge + Halfedge_around_vertex_circulator he = v->vertex_begin(); + Halfedge_around_vertex_circulator end = he; + CGAL_For_all(he,end) + if(he->sharp()) + break; + a = he->opposite()->vertex()->point(); + + // pick next sharp edge + he++; + CGAL_For_all(he,end) + if(he->sharp()) + break; + b = he->opposite()->vertex()->point(); + } + + // return nb of sharp edges incident to v + unsigned int nb_sharp_edges(Vertex_handle v) const + { + Halfedge_around_vertex_circulator he = v->vertex_begin(); + Halfedge_around_vertex_circulator end = he; + unsigned int nb_sharp_edges = 0; + CGAL_For_all(he,end) + if(he->sharp()) + nb_sharp_edges++; + return nb_sharp_edges; + } + + // Angle between two vectors (in degrees) + // we use this formula + // uv = |u||v| cos(u,v) + // u ^ v = w + // |w| = |u||v| |sin(u,v)| + //************************************************** + static double angle_deg(const Vector &u, + const Vector &v) + { + static const double conv = 1.0/PI*180.0; + return conv * angle_rad(u,v); + } + + static FT len(const Vector &v) + { + return (FT)std::sqrt(CGAL_NTS to_double(v*v)); + } + + // Angle between two vectors (in rad) + // uv = |u||v| cos(u,v) + // u ^ v = w + // |w| = |u||v| |sin(u,v)| + //************************************************** + static double angle_rad(const Vector &u, + const Vector &v) + { + // check + double product = len(u)*len(v); + if(product == 0) + return 0.0; + + // cosine + double dot = (u*v); + double cosine = dot / product; + + // sine + Vector w = CGAL::cross_product(u,v); + double AbsSine = len(w) / product; + + if(cosine >= 0) + return std::asin(fix_sine(AbsSine)); + else + return PI-std::asin(fix_sine(AbsSine)); + } + + //********************************************** + // fix sine + //********************************************** + static double fix_sine(double sine) + { + if(sine >= 1) + return 1; + else + if(sine <= -1) + return -1; + else + return sine; + } + + unsigned int tag_sharp_edges(const double angle_sharp) + { + unsigned int nb = 0; + for(Halfedge_iterator he = edges_begin(); + he != edges_end(); + he++) + { + const bool tag = is_sharp(he,angle_sharp); + he->sharp() = tag; + he->opposite()->sharp() = tag; + nb += tag ? 1 : 0; + } + return nb; + } + + // draw edges + void gl_draw_sharp_edges(const float line_width, + unsigned char r, + unsigned char g, + unsigned char b) + { + ::glLineWidth(line_width); + ::glColor3ub(r,g,b); + + ::glBegin(GL_LINES); + for(Halfedge_iterator he = edges_begin(); + he != edges_end(); + he++) + { + if(he->sharp()) + { + const Point& a = he->opposite()->vertex()->point(); + const Point& b = he->vertex()->point(); + ::glVertex3d(a[0],a[1],a[2]); + ::glVertex3d(b[0],b[1],b[2]); + } + } + ::glEnd(); + } + + + void gl_draw_boundaries() + { + ::glBegin(GL_LINES); + for(Halfedge_iterator he = halfedges_begin(); + he != halfedges_end(); + he++) + { + if(he->is_border()) + { + const Point& a = he->vertex()->point(); + const Point& b = he->opposite()->vertex()->point(); + ::glVertex3d(a.x(),a.y(),a.z()); + ::glVertex3d(b.x(),b.y(),b.z()); + } + } + ::glEnd(); + } + + // draw bounding box + void gl_draw_bounding_box() + { + ::glBegin(GL_LINES); + + // along x axis + ::glVertex3f(m_bbox.xmin(),m_bbox.ymin(),m_bbox.zmin()); + ::glVertex3f(m_bbox.xmax(),m_bbox.ymin(),m_bbox.zmin()); + ::glVertex3f(m_bbox.xmin(),m_bbox.ymin(),m_bbox.zmax()); + ::glVertex3f(m_bbox.xmax(),m_bbox.ymin(),m_bbox.zmax()); + ::glVertex3f(m_bbox.xmin(),m_bbox.ymax(),m_bbox.zmin()); + ::glVertex3f(m_bbox.xmax(),m_bbox.ymax(),m_bbox.zmin()); + ::glVertex3f(m_bbox.xmin(),m_bbox.ymax(),m_bbox.zmax()); + ::glVertex3f(m_bbox.xmax(),m_bbox.ymax(),m_bbox.zmax()); + + // along y axis + ::glVertex3f(m_bbox.xmin(),m_bbox.ymin(),m_bbox.zmin()); + ::glVertex3f(m_bbox.xmin(),m_bbox.ymax(),m_bbox.zmin()); + ::glVertex3f(m_bbox.xmin(),m_bbox.ymin(),m_bbox.zmax()); + ::glVertex3f(m_bbox.xmin(),m_bbox.ymax(),m_bbox.zmax()); + ::glVertex3f(m_bbox.xmax(),m_bbox.ymin(),m_bbox.zmin()); + ::glVertex3f(m_bbox.xmax(),m_bbox.ymax(),m_bbox.zmin()); + ::glVertex3f(m_bbox.xmax(),m_bbox.ymin(),m_bbox.zmax()); + ::glVertex3f(m_bbox.xmax(),m_bbox.ymax(),m_bbox.zmax()); + + // along z axis + ::glVertex3f(m_bbox.xmin(),m_bbox.ymin(),m_bbox.zmin()); + ::glVertex3f(m_bbox.xmin(),m_bbox.ymin(),m_bbox.zmax()); + ::glVertex3f(m_bbox.xmin(),m_bbox.ymax(),m_bbox.zmin()); + ::glVertex3f(m_bbox.xmin(),m_bbox.ymax(),m_bbox.zmax()); + ::glVertex3f(m_bbox.xmax(),m_bbox.ymin(),m_bbox.zmin()); + ::glVertex3f(m_bbox.xmax(),m_bbox.ymin(),m_bbox.zmax()); + ::glVertex3f(m_bbox.xmax(),m_bbox.ymax(),m_bbox.zmin()); + ::glVertex3f(m_bbox.xmax(),m_bbox.ymax(),m_bbox.zmax()); + + ::glEnd(); + } + + // count #boundaries + unsigned int nb_boundaries() + { + unsigned int nb = 0; + tag_halfedges(0); + for(Halfedge_iterator he = halfedges_begin(); + he != halfedges_end(); + he++) + { + if(he->is_border() && he->tag() == 0) + { + nb++; + Halfedge_handle curr = he; + do + { + curr = curr->next(); + curr->tag(1); + } + while(curr != he); + } + } + return nb; + } + + // tag component + void tag_component(Facet_handle pSeedFacet, + const int tag_free, + const int tag_done) +{ + pSeedFacet->tag() = tag_done; + std::list facets; + facets.push_front(pSeedFacet); + while(!facets.empty()) + { + Facet_handle pFacet = facets.front(); + facets.pop_front(); + pFacet->tag() = tag_done; + Halfedge_around_facet_circulator pHalfedge = pFacet->facet_begin(); + Halfedge_around_facet_circulator end = pHalfedge; + CGAL_For_all(pHalfedge,end) + { + Facet_handle pNFacet = pHalfedge->opposite()->facet(); + if(pNFacet != NULL && pNFacet->tag() == tag_free) + { + facets.push_front(pNFacet); + pNFacet->tag() = tag_done; + } + } + } +} + + // count #components + unsigned int nb_components() + { + unsigned int nb = 0; + tag_facets(0); + for(Facet_iterator pFacet = facets_begin(); + pFacet != facets_end(); + pFacet++) + { + if(pFacet->tag() == 0) + { + nb++; + tag_component(pFacet,0,1); + } + } + return nb; + } + + // compute the genus + // V - E + F + B = 2 (C - G) + // C -> #connected components + // G : genus + // B : #boundaries + int genus() + { + int c = nb_components(); + int b = nb_boundaries(); + int v = size_of_vertices(); + int e = size_of_halfedges()/2; + int f = size_of_facets(); + return genus(c,v,f,e,b); + } + int genus(int c, + int v, + int f, + int e, + int b) + { + return (2*c+e-b-f-v)/2; + } + +}; + + +#endif // _POLYGON_MESH_