1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
|
#include <CGAL/Simple_cartesian.h>
#include <CGAL/Polyhedron_3.h>
#include <CGAL/IO/Polyhedron_iostream.h>
#include <iostream>
#include <algorithm>
#include <vector>
#include <cmath>
typedef CGAL::Simple_cartesian<double> Kernel;
typedef Kernel::Vector_3 Vector;
typedef Kernel::Point_3 Point;
typedef CGAL::Polyhedron_3<Kernel> Polyhedron;
typedef Polyhedron::Vertex Vertex;
typedef Polyhedron::Vertex_iterator Vertex_iterator;
typedef Polyhedron::Halfedge_handle Halfedge_handle;
typedef Polyhedron::Edge_iterator Edge_iterator;
typedef Polyhedron::Facet_iterator Facet_iterator;
typedef Polyhedron::Halfedge_around_vertex_const_circulator HV_circulator;
typedef Polyhedron::Halfedge_around_facet_circulator HF_circulator;
void create_center_vertex( Polyhedron& P, Facet_iterator f) {
Vector vec( 0.0, 0.0, 0.0);
std::size_t order = 0;
HF_circulator h = f->facet_begin();
do {
vec = vec + ( h->vertex()->point() - CGAL::ORIGIN);
++ order;
} while ( ++h != f->facet_begin());
CGAL_assertion( order >= 3); // guaranteed by definition of polyhedron
Point center = CGAL::ORIGIN + (vec / static_cast<double>(order));
Halfedge_handle new_center = P.create_center_vertex( f->halfedge());
new_center->vertex()->point() = center;
}
struct Smooth_old_vertex {
Point operator()( const Vertex& v) const {
CGAL_precondition((CGAL::circulator_size( v.vertex_begin()) & 1) == 0);
std::size_t degree = CGAL::circulator_size( v.vertex_begin()) / 2;
double alpha = ( 4.0 - 2.0 * std::cos( 2.0 * CGAL_PI / degree)) / 9.0;
Vector vec = (v.point() - CGAL::ORIGIN) * ( 1.0 - alpha);
HV_circulator h = v.vertex_begin();
do {
vec = vec + ( h->opposite()->vertex()->point() - CGAL::ORIGIN)
* alpha / static_cast<double>(degree);
++ h;
CGAL_assertion( h != v.vertex_begin()); // even degree guaranteed
++ h;
} while ( h != v.vertex_begin());
return (CGAL::ORIGIN + vec);
}
};
void flip_edge( Polyhedron& P, Halfedge_handle e) {
Halfedge_handle h = e->next();
P.join_facet( e);
P.split_facet( h, h->next()->next());
}
void subdiv( Polyhedron& P) {
if ( P.size_of_facets() == 0)
return;
// We use that new vertices/halfedges/facets are appended at the end.
std::size_t nv = P.size_of_vertices();
Vertex_iterator last_v = P.vertices_end();
-- last_v; // the last of the old vertices
Edge_iterator last_e = P.edges_end();
-- last_e; // the last of the old edges
Facet_iterator last_f = P.facets_end();
-- last_f; // the last of the old facets
Facet_iterator f = P.facets_begin(); // create new center vertices
do {
create_center_vertex( P, f);
} while ( f++ != last_f);
std::vector<Point> pts; // smooth the old vertices
pts.reserve( nv); // get intermediate space for the new points
++ last_v; // make it the past-the-end position again
std::transform( P.vertices_begin(), last_v, std::back_inserter( pts),
Smooth_old_vertex());
std::copy( pts.begin(), pts.end(), P.points_begin());
Edge_iterator e = P.edges_begin(); // flip the old edges
++ last_e; // make it the past-the-end position again
while ( e != last_e) {
Halfedge_handle h = e;
++e; // careful, incr. before flip since flip destroys current edge
flip_edge( P, h);
};
CGAL_postcondition( P.is_valid());
}
int main() {
Polyhedron P;
std::cin >> P;
P.normalize_border();
if ( P.size_of_border_edges() != 0) {
std::cerr << "The input object has border edges. Cannot subdivide."
<< std::endl;
std::exit(1);
}
subdiv( P);
std::cout << P;
return 0;
}
|
