diff --git a/Packages/Tutorial/tutorial/Polyhedron/doc/makefile b/Packages/Tutorial/tutorial/Polyhedron/doc/makefile index 65b0bca1c06..e3ba8c575dd 100644 --- a/Packages/Tutorial/tutorial/Polyhedron/doc/makefile +++ b/Packages/Tutorial/tutorial/Polyhedron/doc/makefile @@ -14,12 +14,12 @@ all: make ps make pdf make clean - \rm *.ps + \rm *.ps *.tmp clear ls -l clean: - \rm *.aux *.dvi *.log + \rm *.aux *.dvi *.log *.out *.blg *.bbl superclean: \rm *.ps *.pdf *~ diff --git a/Packages/Tutorial/tutorial/Polyhedron/doc/tutorial.pdf b/Packages/Tutorial/tutorial/Polyhedron/doc/tutorial.pdf index 23f591b7d61..4f266a92586 100644 Binary files a/Packages/Tutorial/tutorial/Polyhedron/doc/tutorial.pdf and b/Packages/Tutorial/tutorial/Polyhedron/doc/tutorial.pdf differ diff --git a/Packages/Tutorial/tutorial/Polyhedron/doc/tutorial.tex b/Packages/Tutorial/tutorial/Polyhedron/doc/tutorial.tex index 246a3beeb61..b8e907dc9ca 100644 --- a/Packages/Tutorial/tutorial/Polyhedron/doc/tutorial.tex +++ b/Packages/Tutorial/tutorial/Polyhedron/doc/tutorial.tex @@ -931,6 +931,398 @@ Polyhedron (see Section~\ref{sec:builder}). Subdivision engine +{ \scriptsize +\begin{verbatim} + +#include "enriched_polyhedron.h" +#include "builder.h" + +template +class CModifierQuadTriangle : public CGAL::Modifier_base +{ +private: + typedef typename kernel::FT FT; + typedef typename HDS::Vertex Vertex; + typedef typename Vertex::Point Point; + typedef typename HDS::Face_handle Face_handle; + typedef typename HDS::Halfedge_handle Halfedge_handle; + typedef typename CGAL::Enriched_polyhedron_incremental_builder_3 builder; + typedef typename Polyhedron::Halfedge_around_vertex_circulator Halfedge_around_vertex_circulator; + typedef typename Polyhedron::Facet_handle Facet_handle; + typedef typename Polyhedron::Vertex_handle Vertex_handle; + + Polyhedron *m_pMesh; + +public: + + // life cycle + CModifierQuadTriangle(Polyhedron *pMesh) + { + CGAL_assertion(pMesh != NULL); + m_pMesh = pMesh; + } + ~CModifierQuadTriangle() {} + + // subdivision + void operator()( HDS& hds) + { + builder B(hds,true); + B.begin_surface(3,1,6); + add_vertices(B); + add_facets(B); + B.end_surface(); + } + +private: + + //*************************************************** + // add vertices + //*************************************************** + void add_vertices(builder &B) + { + // put original vertices + int index = 0; + typename Polyhedron::Vertex_iterator pVertex; + for(pVertex = m_pMesh->vertices_begin(); + pVertex != m_pMesh->vertices_end(); + pVertex++) + { + pVertex->tag(index); + B.add_vertex(pVertex->point()); + index++; + } + + // as many as #edges + m_pMesh->tag_halfedges(-1); + typename Polyhedron::Halfedge_iterator pHalfedge; + for(pHalfedge = m_pMesh->halfedges_begin(); + pHalfedge != m_pMesh->halfedges_end(); + pHalfedge++) + { + if(pHalfedge->tag() != -1) + continue; + + // simple edge bissection + const Point& p1 = pHalfedge->vertex()->point(); + const Point& p2 = pHalfedge->opposite()->vertex()->point(); + Point point = Point(0.5f*(p1.x()+p2.x()), + 0.5f*(p1.y()+p2.y()), + 0.5f*(p1.z()+p2.z())); + B.add_vertex(point); + + // put vertex indices on both halfedges + pHalfedge->tag(index); + pHalfedge->opposite()->tag(index); + index++; + } + + // and as many as #facets with degree > 3 + m_pMesh->tag_facets(-1); + typename Polyhedron::Facet_iterator pFacet; + for(pFacet = m_pMesh->facets_begin(); + pFacet != m_pMesh->facets_end(); + pFacet++) + { + unsigned int degree = Polyhedron::degree(pFacet); + CGAL_assertion(degree >= 3); + // forget about triangles, they are + // simply 1-to-4 subdivided by edge bisection + if(degree == 3) + continue; + + // barycentric subdivision + Point barycenter; + m_pMesh->compute_facet_center(pFacet,barycenter); + B.add_vertex(barycenter); + pFacet->tag(index); + index++; + } + } + + //*************************************************** + // add facets + //*************************************************** + void add_facets(builder &B) + { + typename Polyhedron::Facet_iterator pFacet; + for(pFacet = m_pMesh->facets_begin(); + pFacet != m_pMesh->facets_end(); + pFacet++) + { + unsigned int degree = Polyhedron::degree(pFacet); + CGAL_assertion(degree >= 3); + + if(degree == 3) + { + typename Polyhedron::Halfedge_handle pHalfedge = pFacet->halfedge(); + int i0 = pHalfedge->tag(); + int i1 = pHalfedge->vertex()->tag(); + int i2 = pHalfedge->next()->tag(); + int i3 = pHalfedge->next()->vertex()->tag(); + int i4 = pHalfedge->next()->next()->tag(); + int i5 = pHalfedge->next()->next()->vertex()->tag(); + CGAL_assertion(i0 >= 0); + CGAL_assertion(i1 >= 0); + CGAL_assertion(i2 >= 0); + CGAL_assertion(i3 >= 0); + CGAL_assertion(i4 >= 0); + CGAL_assertion(i5 >= 0); + + // create 4 triangles + B.begin_facet(); + B.add_vertex_to_facet(i1); + B.add_vertex_to_facet(i2); + B.add_vertex_to_facet(i0); + const Halfedge_handle& h1 = B.end_facet(); + B.begin_facet(); + B.add_vertex_to_facet(i3); + B.add_vertex_to_facet(i4); + B.add_vertex_to_facet(i2); + const Halfedge_handle& h2 = B.end_facet(); + B.begin_facet(); + B.add_vertex_to_facet(i5); + B.add_vertex_to_facet(i0); + B.add_vertex_to_facet(i4); + const Halfedge_handle& h3 = B.end_facet(); + + // center face + B.begin_facet(); + B.add_vertex_to_facet(i0); + B.add_vertex_to_facet(i2); + B.add_vertex_to_facet(i4); + Halfedge_handle h4 = B.end_facet(); + + h1->control_edge(false); + h1->next()->control_edge(false); + h1->next()->next()->control_edge(false); + + h2->control_edge(false); + h2->next()->control_edge(false); + h2->next()->next()->control_edge(false); + + h3->control_edge(false); + h3->next()->control_edge(false); + h3->next()->next()->control_edge(false); + + h4->control_edge(false); + h4->next()->control_edge(false); + h4->next()->next()->control_edge(false); + + if(pHalfedge->control_edge()) + { + h1->control_edge(true); + h3->next()->control_edge(true); + } + if(pHalfedge->next()->control_edge()) + { + h1->next()->control_edge(true); + h2->control_edge(true); + } + if(pHalfedge->next()->next()->control_edge()) + { + h2->next()->control_edge(true); + h3->control_edge(true); + } + } + else + { + // i1: index of barycenter vertex + int i1 = pFacet->tag(); + CGAL_assertion(i1 >= 0); + + // for each halfedge + typename Polyhedron::Halfedge_around_facet_circulator h; + h = pFacet->facet_begin(); + do + { + // i2,3,4: indices of three consecutive + // vertices on halfedges + int i2 = h->tag(); + int i3 = h->vertex()->tag(); + int i4 = h->next()->tag(); + CGAL_assertion(i2 >= 0); + CGAL_assertion(i3 >= 0); + CGAL_assertion(i4 >= 0); + + // create a quad + B.begin_facet(); + B.add_vertex_to_facet(i3); + B.add_vertex_to_facet(i4); + B.add_vertex_to_facet(i1); + B.add_vertex_to_facet(i2); + const Halfedge_handle& pNewHalfedge = B.end_facet(); + + pNewHalfedge->control_edge(false); + pNewHalfedge->next()->control_edge(false); + pNewHalfedge->next()->next()->control_edge(false); + pNewHalfedge->next()->next()->next()->control_edge(false); + + if(h->control_edge()) + pNewHalfedge->control_edge(true); + if(h->next()->control_edge()) + pNewHalfedge->next()->control_edge(true); + } + while(++h != pFacet->facet_begin()); + } + } + } + + //*************************************************** + // correcting factor + //*************************************************** + static float correcting_factor(unsigned int ne, + unsigned int nq) + { + if(ne == 2 && nq == 1) + return -0.20505f; + if(ne == 3 && nq == 1) + return 0.80597f; + if(ne == 3 && nq == 2) + return 0.61539f; + if(ne == 4 && nq == 1) + return 0.34792f; + if(ne == 4 && nq == 2) + return 0.21380f; + if(ne == 4 && nq == 3) + return 0.10550f; + return 0.0f; + } +public: + //*************************************************** + // smooth vertex positions + //*************************************************** + static void smooth(Polyhedron *pMesh, + bool smooth_boundary = true) + { + CGAL_assertion(pMesh != NULL); + + // alloc position vectors + unsigned int nb_vertices = pMesh->size_of_vertices(); + FT *pPos = new FT[3*nb_vertices]; + CGAL_assertion(pPos != NULL); + + // compute new positions + unsigned int index = 0; + typename Polyhedron::Vertex_iterator pVertex; + for(pVertex = pMesh->vertices_begin(); + pVertex != pMesh->vertices_end(); + pVertex++) + { + // border vertices will not move + if(Polyhedron::is_border(pVertex)) + { + // do not smooth it + const Point& curr = pVertex->point(); + if(!smooth_boundary) + { + pPos[3*index] = curr.x(); + pPos[3*index+1] = curr.y(); + pPos[3*index+2] = curr.z(); + } + // smooth using [1/4 1/2 1/4] cubic B-spline averaging mask + else + { + const typename Polyhedron::Halfedge_handle& pHalfedge = + pMesh->get_border_halfedge(pVertex); + CGAL_assertion(pHalfedge != NULL); + const Point& next = pHalfedge->vertex()->point(); + const Point& prev = pHalfedge->prev()->prev()->vertex()->point(); + pPos[3*index] = 0.25f*prev.x() + 0.5f*curr.x() + 0.25f*next.x(); + pPos[3*index+1] = 0.25f*prev.y() + 0.5f*curr.y() + 0.25f*next.y(); + pPos[3*index+2] = 0.25f*prev.z() + 0.5f*curr.z() + 0.25f*next.z(); + } + } // end is border + else + { + + unsigned int nb_quads = 0; + unsigned int nb_edges = 0; + + // rotate around vertex to count #edges and #quads + Halfedge_around_vertex_circulator + pHalfEdge = pVertex->vertex_begin(); + Halfedge_around_vertex_circulator end = pHalfEdge; + CGAL_For_all(pHalfEdge,end) + { + const Facet_handle& pFacet = pHalfEdge->facet(); + CGAL_assertion(pFacet != NULL); + unsigned int degree = Polyhedron::degree(pFacet); + CGAL_assertion(degree == 4 || degree == 3); + if(degree == 4) + nb_quads++; + nb_edges++; + } + + // compute coefficients + FT ne = (FT)nb_edges; + FT nq = (FT)nb_quads; + FT alpha = 1.0f / (1.0f + ne/2.0f + nq/4.0f); + FT beta = alpha / 2.0f; // edges + FT gamma = alpha / 4.0f; // corners of incident quads + FT eta = correcting_factor(nb_edges,nb_quads); + + // new position + pPos[3*index] = alpha * pVertex->point().x(); + pPos[3*index+1] = alpha * pVertex->point().y(); + pPos[3*index+2] = alpha * pVertex->point().z(); + + // rotate around vertex to compute new position + pHalfEdge = pVertex->vertex_begin(); + end = pHalfEdge; + CGAL_For_all(pHalfEdge,end) + { + const Facet_handle& pFacet = pHalfEdge->facet(); + CGAL_assertion(pFacet != NULL); + unsigned int degree = Polyhedron::degree(pFacet); + CGAL_assertion(degree == 4 || degree == 3); + + // add edge-vertex contribution + const Point& point = pHalfEdge->prev()->vertex()->point(); + pPos[3*index] += beta * point.x(); + pPos[3*index+1] += beta * point.y(); + pPos[3*index+2] += beta * point.z(); + + // add corner vertex contribution + if(degree == 4) + { + const Point& corner = pHalfEdge->next()->next()->vertex()->point(); + pPos[3*index] += gamma * corner.x(); + pPos[3*index+1] += gamma * corner.y(); + pPos[3*index+2] += gamma * corner.z(); + } + } + + // apply correction + pPos[3*index] = pPos[3*index] + + eta*(pPos[3*index]-pVertex->point().x()); + pPos[3*index+1] = pPos[3*index+1] + + eta*(pPos[3*index+1]-pVertex->point().y()); + pPos[3*index+2] = pPos[3*index+2] + + eta*(pPos[3*index+2]-pVertex->point().z()); + + } // end !is border + index++; + } + + // set new positions + index = 0; + for(pVertex = pMesh->vertices_begin(); + pVertex != pMesh->vertices_end(); + pVertex++) + { + Point& point = pVertex->point(); + point = Point(pPos[3*index], + pPos[3*index+1], + pPos[3*index+2]); + index++; + } + + // cleanup + delete [] pPos; + } +}; +\end{verbatim}} + { \scriptsize \begin{verbatim} template @@ -961,42 +1353,6 @@ public: }; \end{verbatim}} -Subdivision using a modified incremental builder - -{ \scriptsize -\begin{verbatim} -template -class CModifierQuadTriangle : public CGAL::Modifier_base -{ -private: - Polyhedron *m_pMesh; - typedef typename CGAL::Enriched_builder builder; - -public: - - // life cycle - CModifierQuadTriangle(Polyhedron *pMesh) - { - m_pMesh = pMesh; - } - ~CModifierQuadTriangle() {} - - // subdivision - void operator()( HDS& hds) - { - builder B(hds,true); - B.begin_surface(3,1,6); - add_vertices(B); - add_facets(B); - B.end_surface(); - } -... -}; -\end{verbatim}} - - - - % SurfLab \subsection{Subdivision using a rule template}