| 688 |
khor |
1 |
/* ----------------------------------------------------------------------- *
|
|
|
2 |
* This file is part of GEL, http://www.imm.dtu.dk/GEL
|
|
|
3 |
* Copyright (C) the authors and DTU Informatics
|
|
|
4 |
* For license and list of authors, see ../../doc/intro.pdf
|
|
|
5 |
* ----------------------------------------------------------------------- */
|
|
|
6 |
|
|
|
7 |
#include "cleanup.h"
|
|
|
8 |
|
|
|
9 |
#include "../CGLA/Vec3f.h"
|
|
|
10 |
#include "../CGLA/Vec3d.h"
|
|
|
11 |
#include "../Geometry/QEM.h"
|
|
|
12 |
#include "../Geometry/KDTree.h"
|
|
|
13 |
|
|
|
14 |
#include "Manifold.h"
|
|
|
15 |
|
|
|
16 |
namespace HMesh
|
|
|
17 |
{
|
|
|
18 |
using namespace std;
|
|
|
19 |
using namespace CGLA;
|
|
|
20 |
using namespace Geometry;
|
|
|
21 |
|
|
|
22 |
namespace
|
|
|
23 |
{
|
|
|
24 |
/// get the minimum angle between 3 vertices
|
|
|
25 |
float min_angle(const Vec3d& v0, const Vec3d& v1, const Vec3d& v2)
|
|
|
26 |
{
|
|
|
27 |
Vec3d a = normalize(v1-v0);
|
|
|
28 |
Vec3d b = normalize(v2-v1);
|
|
|
29 |
Vec3d c = normalize(v0-v2);
|
|
|
30 |
|
|
|
31 |
return min(dot(a, -c), min(dot(b, -a), dot(c, -b)));
|
|
|
32 |
}
|
|
|
33 |
|
|
|
34 |
/// need description
|
|
|
35 |
float edge_error(const Manifold& m, HalfEdgeID h, const Vec3d& pa, const Vec3d& pb)
|
|
|
36 |
{
|
|
|
37 |
QEM q;
|
|
|
38 |
Walker j = m.walker(h);
|
|
|
39 |
|
|
|
40 |
FaceID f = j.face();
|
|
|
41 |
if(f != InvalidFaceID)
|
|
|
42 |
q += QEM(Vec3d(0), Vec3d(normal(m, f)));
|
|
|
43 |
|
|
|
44 |
f = j.opp().face();
|
|
|
45 |
if(f != InvalidFaceID)
|
|
|
46 |
q += QEM(Vec3d(0), Vec3d(normal(m, f)));
|
|
|
47 |
|
|
|
48 |
return q.error(pb - pa);
|
|
|
49 |
}
|
|
|
50 |
|
|
|
51 |
/// need description
|
|
|
52 |
float vertex_error(const Manifold& m, VertexID v, const Vec3d& pb)
|
|
|
53 |
{
|
|
|
54 |
QEM q;
|
|
|
55 |
Vec3d pa(m.pos(v));
|
|
|
56 |
|
|
|
57 |
for(Walker vj = m.walker(v); !vj.full_circle(); vj = vj.circulate_vertex_cw()){
|
|
|
58 |
FaceID f = vj.face();
|
|
|
59 |
if(f != InvalidFaceID)
|
|
|
60 |
q += QEM(Vec3d(0), Vec3d(normal(m, f)));
|
|
|
61 |
}
|
|
|
62 |
return q.error(pb - pa);
|
|
|
63 |
}
|
|
|
64 |
}
|
|
|
65 |
|
|
|
66 |
void remove_caps(Manifold& m, float thresh)
|
|
|
67 |
{
|
|
|
68 |
for(FaceIDIterator f = m.faces_begin(); f != m.faces_end(); ++f){
|
|
|
69 |
Vec3d p[3];
|
|
|
70 |
HalfEdgeID he[3];
|
|
|
71 |
VertexID vh[3];
|
|
|
72 |
|
|
|
73 |
// store ids of vertices and halfedges and vertex positions of face
|
|
|
74 |
size_t n = 0;
|
|
|
75 |
for(Walker fj = m.walker(*f); !fj.full_circle(); fj = fj.circulate_face_cw(), ++n){
|
|
|
76 |
vh[n] = fj.vertex();
|
|
|
77 |
p[n] = Vec3d(m.pos(vh[n]));
|
|
|
78 |
|
|
|
79 |
// original face circulator implementation returned next halfedge, jumper doesn't. Can this be optimized?
|
|
|
80 |
he[n] = fj.halfedge();
|
|
|
81 |
}
|
|
|
82 |
assert(n == 3);
|
|
|
83 |
|
|
|
84 |
// calculate the edge lengths of face
|
|
|
85 |
bool is_collapsed = false;
|
|
|
86 |
Vec3d edges[3];
|
|
|
87 |
for(size_t i = 0; i < 3; ++i){
|
|
|
88 |
edges[i] = p[(i+1)%3] - p[i];
|
|
|
89 |
float l = length(edges[i]);
|
|
|
90 |
if(l < 1e-20)
|
|
|
91 |
is_collapsed = true;
|
|
|
92 |
else
|
|
|
93 |
edges[i] /= l;
|
|
|
94 |
|
|
|
95 |
}
|
|
|
96 |
// an edge length was close to 1e-20, thus collapsed
|
|
|
97 |
if(is_collapsed)
|
|
|
98 |
continue;
|
|
|
99 |
|
|
|
100 |
for(size_t i = 0; i < 3; ++i){
|
|
|
101 |
float ang = acos(max(-1.0, min(1.0, dot(-edges[(i+2)%3], edges[i]))));
|
|
|
102 |
|
|
|
103 |
// flip long edge of current face if angle exceeds cap threshold and result is better than current cap
|
|
|
104 |
if(ang > thresh){
|
|
|
105 |
size_t iplus1 = (i+1)%3;
|
|
|
106 |
Vec3d edge_dir = edges[iplus1];
|
|
|
107 |
|
|
|
108 |
Walker j = m.walker(he[iplus1]);
|
|
|
109 |
Vec3d v0(m.pos(j.vertex()));
|
|
|
110 |
Vec3d v1(m.pos(j.next().vertex()));
|
|
|
111 |
Vec3d v2(m.pos(j.opp().vertex()));
|
|
|
112 |
Vec3d v3(m.pos(j.opp().next().vertex()));
|
|
|
113 |
|
|
|
114 |
float m1 = min(min_angle(v0, v1, v2), min_angle(v0, v2, v3));
|
|
|
115 |
float m2 = min(min_angle(v0, v1, v3), min_angle(v1, v2, v3));
|
|
|
116 |
|
|
|
117 |
if(m1 < m2){
|
|
|
118 |
// If the "cap vertex" projected onto the long edge is better in the
|
|
|
119 |
// sense that there is less geometric error after the flip, then we
|
|
|
120 |
// use the projected vertex. In other words, we see if it pays to snap
|
|
|
121 |
// to the edge.
|
|
|
122 |
Vec3d pprj = edge_dir * dot(edge_dir, p[i]-p[iplus1])+p[iplus1];
|
|
|
123 |
if(edge_error(m, he[iplus1], pprj, Vec3d(m.pos(vh[i]))) > vertex_error(m, vh[i], pprj))
|
|
|
124 |
m.pos(vh[i]) = pprj;
|
|
|
125 |
// flip if legal
|
|
|
126 |
if(precond_flip_edge(m, he[iplus1]))
|
|
|
127 |
m.flip_edge(he[iplus1]);
|
|
|
128 |
break;
|
|
|
129 |
}
|
|
|
130 |
}
|
|
|
131 |
}
|
|
|
132 |
}
|
|
|
133 |
|
|
|
134 |
}
|
|
|
135 |
|
|
|
136 |
void remove_needles(Manifold& m, float thresh)
|
|
|
137 |
{
|
|
|
138 |
bool did_work = false;
|
|
|
139 |
|
|
|
140 |
// remove needles until no more can be removed
|
|
|
141 |
do{
|
|
|
142 |
did_work = false;
|
|
|
143 |
for(VertexIDIterator v = m.vertices_begin(); v != m.vertices_end(); ++v){
|
|
|
144 |
// don't attempt to remove needles if vertex is boundary
|
|
|
145 |
if(boundary(m, *v))
|
|
|
146 |
continue;
|
|
|
147 |
|
|
|
148 |
for(Walker vj = m.walker(*v); !vj.full_circle(); vj = vj.circulate_vertex_cw()){
|
|
|
149 |
// don't attempt to remove needles if vertex of jumper halfedge is boundary
|
|
|
150 |
// if(boundary(m, vj.vertex()))
|
|
|
151 |
// continue;
|
|
|
152 |
|
|
|
153 |
HalfEdgeID h = vj.opp().halfedge();
|
|
|
154 |
// VertexID n = vj.vertex();
|
|
|
155 |
float dist = length(m, h);
|
|
|
156 |
|
|
|
157 |
// collapse edge if allowed and needle is present
|
|
|
158 |
if(dist < thresh && precond_collapse_edge(m, h)){
|
|
|
159 |
// if(vertex_error(m, *v, Vec3d(m.pos(n))) < vertex_error(m, n, Vec3d(m.pos(*v))))
|
|
|
160 |
// m.pos(*v) = m.pos(n);
|
|
|
161 |
m.collapse_edge(h);
|
|
|
162 |
did_work = true;
|
|
|
163 |
break;
|
|
|
164 |
}
|
|
|
165 |
}
|
|
|
166 |
}
|
|
|
167 |
}
|
|
|
168 |
while(did_work);
|
|
|
169 |
}
|
|
|
170 |
|
|
|
171 |
int remove_duplicates(Manifold& m, double rad)
|
|
|
172 |
{
|
|
|
173 |
KDTree<Vec3d, VertexID> vtree;
|
|
|
174 |
|
|
|
175 |
for(VertexID v: m.vertices())
|
|
|
176 |
vtree.insert(m.pos(v), v);
|
|
|
177 |
vtree.build();
|
|
|
178 |
|
|
|
179 |
VertexAttributeVector<int> vertex_dup(m.allocated_vertices(),0);
|
|
|
180 |
vector<vector<HalfEdgeID> > clustered_halfedges;
|
|
|
181 |
|
|
|
182 |
int cnt = 0;
|
|
|
183 |
for(VertexID v: m.vertices())
|
|
|
184 |
{
|
|
|
185 |
vector<Vec3d> keys;
|
|
|
186 |
vector<VertexID> vals;
|
|
|
187 |
int n = vtree.in_sphere(m.pos(v), rad, keys, vals);
|
|
|
188 |
|
|
|
189 |
if(n>1) {
|
|
|
190 |
sort(vals.begin(), vals.end());
|
|
|
191 |
|
|
|
192 |
auto v2remove = vals.begin();
|
|
|
193 |
for(v2remove++; v2remove != vals.end(); v2remove++) {
|
|
|
194 |
m.remove_vertex(*v2remove);
|
|
|
195 |
++cnt;
|
|
|
196 |
}
|
|
|
197 |
}
|
|
|
198 |
}
|
|
|
199 |
return cnt;
|
|
|
200 |
}
|
|
|
201 |
|
|
|
202 |
|
|
|
203 |
int stitch_mesh(Manifold& m, double rad)
|
|
|
204 |
{
|
|
|
205 |
KDTree<Vec3d, VertexID> vtree;
|
|
|
206 |
|
|
|
207 |
for(VertexIDIterator vid = m.vertices_begin(); vid != m.vertices_end(); ++vid)
|
|
|
208 |
if(boundary(m, *vid))
|
|
|
209 |
vtree.insert(m.pos(*vid), *vid);
|
|
|
210 |
vtree.build();
|
|
|
211 |
|
|
|
212 |
VertexAttributeVector<int> cluster_id(m.allocated_vertices(),-1);
|
|
|
213 |
vector<vector<HalfEdgeID> > clustered_halfedges;
|
|
|
214 |
|
|
|
215 |
int cluster_ctr=0;
|
|
|
216 |
for(VertexIDIterator vid = m.vertices_begin(); vid != m.vertices_end(); ++vid)
|
|
|
217 |
if(boundary(m, *vid) && cluster_id[*vid] == -1)
|
|
|
218 |
{
|
|
|
219 |
vector<Vec3d> keys;
|
|
|
220 |
vector<VertexID> vals;
|
|
|
221 |
int n = vtree.in_sphere(m.pos(*vid), rad, keys, vals);
|
|
|
222 |
|
|
|
223 |
vector<HalfEdgeID> boundary_edges;
|
|
|
224 |
for(int i=0;i<n;++i)
|
|
|
225 |
{
|
|
|
226 |
cluster_id[vals[i]] = cluster_ctr;
|
|
|
227 |
Walker w = m.walker(vals[i]);
|
|
|
228 |
boundary_edges.push_back(w.halfedge());
|
|
|
229 |
}
|
|
|
230 |
clustered_halfedges.push_back(boundary_edges);
|
|
|
231 |
++cluster_ctr;
|
|
|
232 |
}
|
|
|
233 |
|
|
|
234 |
int unstitched=0;
|
|
|
235 |
for(HalfEdgeIDIterator hid = m.halfedges_begin(); hid != m.halfedges_end(); ++hid)
|
|
|
236 |
{
|
|
|
237 |
HalfEdgeID h0 = *hid;
|
|
|
238 |
Walker w = m.walker(h0);
|
|
|
239 |
if(w.face() == InvalidFaceID)
|
|
|
240 |
{
|
|
|
241 |
VertexID v0 = w.opp().vertex();
|
|
|
242 |
VertexID v1 = w.vertex();
|
|
|
243 |
|
|
|
244 |
int cid = cluster_id[v1];
|
|
|
245 |
vector<HalfEdgeID>& stitch_candidates = clustered_halfedges[cid];
|
|
|
246 |
size_t i=0;
|
|
|
247 |
for(;i<stitch_candidates.size(); ++i)
|
|
|
248 |
{
|
|
|
249 |
HalfEdgeID h1 = stitch_candidates[i];
|
|
|
250 |
if(m.in_use(h1))
|
|
|
251 |
{
|
|
|
252 |
Walker w = m.walker(h1);
|
|
|
253 |
if(cluster_id[w.vertex()] == cluster_id[v0])
|
|
|
254 |
if(m.stitch_boundary_edges(h0,h1))
|
|
|
255 |
break;
|
|
|
256 |
}
|
|
|
257 |
|
|
|
258 |
}
|
|
|
259 |
if(i == stitch_candidates.size())
|
|
|
260 |
++unstitched;
|
|
|
261 |
}
|
|
|
262 |
}
|
|
|
263 |
return unstitched;
|
|
|
264 |
}
|
|
|
265 |
|
|
|
266 |
void remove_valence_two_vertices(Manifold & m)
|
|
|
267 |
{
|
|
|
268 |
for(VertexID v: m.vertices())
|
|
|
269 |
if(valency(m,v)==2)
|
|
|
270 |
{
|
|
|
271 |
Walker w = m.walker(v);
|
|
|
272 |
if(precond_collapse_edge(m, w.halfedge()))
|
|
|
273 |
m.collapse_edge(w.halfedge());
|
|
|
274 |
}
|
|
|
275 |
}
|
|
|
276 |
|
|
|
277 |
|
|
|
278 |
void stitch_more(Manifold& mani, double rad)
|
|
|
279 |
{
|
|
|
280 |
int unstitched, iter;
|
|
|
281 |
for(iter=0;iter<2;++iter)
|
|
|
282 |
{
|
|
|
283 |
unstitched = stitch_mesh(mani, rad);
|
|
|
284 |
if(unstitched == 0)
|
|
|
285 |
break;
|
|
|
286 |
|
|
|
287 |
vector<HalfEdgeID> hvec;
|
|
|
288 |
for(HalfEdgeIDIterator h = mani.halfedges_begin(); h != mani.halfedges_end();++h)
|
|
|
289 |
if(mani.walker(*h).face() == InvalidFaceID)
|
|
|
290 |
hvec.push_back(*h);
|
|
|
291 |
for(size_t i=0;i<hvec.size(); ++i)
|
|
|
292 |
mani.split_edge(hvec[i]);
|
|
|
293 |
|
|
|
294 |
}
|
|
|
295 |
}
|
|
|
296 |
|
|
|
297 |
|
|
|
298 |
void close_holes(Manifold& m)
|
|
|
299 |
{
|
|
|
300 |
for(HalfEdgeIDIterator h = m.halfedges_begin(); h != m.halfedges_end(); ++h){
|
|
|
301 |
m.close_hole(*h);
|
|
|
302 |
}
|
|
|
303 |
}
|
|
|
304 |
|
|
|
305 |
void flip_orientation(Manifold& m)
|
|
|
306 |
{
|
|
|
307 |
vector<Vec3d> vertices;
|
|
|
308 |
VertexAttributeVector<int> idvec;
|
|
|
309 |
int k=0;
|
|
|
310 |
for(VertexID v: m.vertices()) {
|
|
|
311 |
idvec[v] = k++;
|
|
|
312 |
vertices.push_back(m.pos(v));
|
|
|
313 |
}
|
|
|
314 |
vector<int> indices;
|
|
|
315 |
vector<int> faces;
|
|
|
316 |
for(FaceID f: m.faces()) {
|
|
|
317 |
faces.push_back(no_edges(m, f));
|
|
|
318 |
circulate_face_cw(m, f, (std::function<void(VertexID)>)[&](VertexID v){
|
|
|
319 |
indices.push_back(idvec[v]);
|
|
|
320 |
});
|
|
|
321 |
}
|
|
|
322 |
m.clear();
|
|
|
323 |
m.build(vertices.size(), vertices[0].get(), faces.size(), &faces[0], &indices[0]);
|
|
|
324 |
}
|
|
|
325 |
|
|
|
326 |
}
|