1 | /************************************************************************
|
---|
2 |
|
---|
3 | Surface simplification using quadric error metrics
|
---|
4 |
|
---|
5 | Copyright (C) 1998 Michael Garland. See "COPYING.txt" for details.
|
---|
6 |
|
---|
7 | $Id: MxQSlim.cxx,v 1.1 2002/09/24 16:53:54 wimmer Exp $
|
---|
8 |
|
---|
9 | ************************************************************************/
|
---|
10 |
|
---|
11 | #include "stdmix.h"
|
---|
12 | #include "MxQSlim.h"
|
---|
13 | #include "MxGeom3D.h"
|
---|
14 | #include "MxVector.h"
|
---|
15 |
|
---|
16 | typedef MxQuadric3 Quadric;
|
---|
17 |
|
---|
18 | MxQSlim::MxQSlim(MxStdModel& _m)
|
---|
19 | : MxStdSlim(&_m),
|
---|
20 | quadrics(_m.vert_count())
|
---|
21 | {
|
---|
22 | // Externally visible variables
|
---|
23 | object_transform = NULL;
|
---|
24 | }
|
---|
25 |
|
---|
26 | void MxQSlim::initialize()
|
---|
27 | {
|
---|
28 | collect_quadrics();
|
---|
29 | if( boundary_weight > 0.0 )
|
---|
30 | constrain_boundaries();
|
---|
31 | if( object_transform )
|
---|
32 | transform_quadrics(*object_transform);
|
---|
33 |
|
---|
34 | is_initialized = true;
|
---|
35 | }
|
---|
36 |
|
---|
37 | void MxQSlim::collect_quadrics()
|
---|
38 | {
|
---|
39 | uint j;
|
---|
40 |
|
---|
41 | for(j=0; j<quadrics.length(); j++)
|
---|
42 | quadrics(j).clear();
|
---|
43 |
|
---|
44 | for(MxFaceID i=0; i<m->face_count(); i++)
|
---|
45 | {
|
---|
46 | MxFace& f = m->face(i);
|
---|
47 |
|
---|
48 | Vec3 v1(m->vertex(f(0)));
|
---|
49 | Vec3 v2(m->vertex(f(1)));
|
---|
50 | Vec3 v3(m->vertex(f(2)));
|
---|
51 |
|
---|
52 | Vec4 p = (weighting_policy==MX_WEIGHT_RAWNORMALS) ?
|
---|
53 | triangle_raw_plane(v1, v2, v3):
|
---|
54 | triangle_plane(v1, v2, v3);
|
---|
55 | Quadric Q(p[X], p[Y], p[Z], p[W], m->compute_face_area(i));
|
---|
56 |
|
---|
57 | switch( weighting_policy )
|
---|
58 | {
|
---|
59 | case MX_WEIGHT_ANGLE:
|
---|
60 | for(j=0; j<3; j++)
|
---|
61 | {
|
---|
62 | Quadric Q_j = Q;
|
---|
63 | Q_j *= m->compute_corner_angle(i, j);
|
---|
64 | quadrics(f[j]) += Q_j;
|
---|
65 | }
|
---|
66 | break;
|
---|
67 | case MX_WEIGHT_AREA:
|
---|
68 | case MX_WEIGHT_AREA_AVG:
|
---|
69 | Q *= Q.area();
|
---|
70 | // no break: fallthrough
|
---|
71 | default:
|
---|
72 | quadrics(f[0]) += Q;
|
---|
73 | quadrics(f[1]) += Q;
|
---|
74 | quadrics(f[2]) += Q;
|
---|
75 | break;
|
---|
76 | }
|
---|
77 | }
|
---|
78 | }
|
---|
79 |
|
---|
80 | void MxQSlim::transform_quadrics(const Mat4& P)
|
---|
81 | {
|
---|
82 | for(uint j=0; j<quadrics.length(); j++)
|
---|
83 | quadrics(j).transform(P);
|
---|
84 | }
|
---|
85 |
|
---|
86 | void MxQSlim::discontinuity_constraint(MxVertexID i, MxVertexID j,
|
---|
87 | const MxFaceList& faces)
|
---|
88 | {
|
---|
89 | for(uint f=0; f<faces.length(); f++)
|
---|
90 | {
|
---|
91 | Vec3 org(m->vertex(i)), dest(m->vertex(j));
|
---|
92 | Vec3 e = dest - org;
|
---|
93 |
|
---|
94 | Vec3 n;
|
---|
95 | m->compute_face_normal(faces[f], n);
|
---|
96 |
|
---|
97 | Vec3 n2 = e ^ n;
|
---|
98 | unitize(n2);
|
---|
99 |
|
---|
100 | MxQuadric3 Q(n2, -(n2*org));
|
---|
101 | Q *= boundary_weight;
|
---|
102 |
|
---|
103 | if( weighting_policy == MX_WEIGHT_AREA ||
|
---|
104 | weighting_policy == MX_WEIGHT_AREA_AVG )
|
---|
105 | {
|
---|
106 | Q.set_area(norm2(e));
|
---|
107 | Q *= Q.area();
|
---|
108 | }
|
---|
109 |
|
---|
110 | quadrics(i) += Q;
|
---|
111 | quadrics(j) += Q;
|
---|
112 | }
|
---|
113 | }
|
---|
114 |
|
---|
115 | void MxQSlim::constrain_boundaries()
|
---|
116 | {
|
---|
117 | MxVertexList star;
|
---|
118 | MxFaceList faces;
|
---|
119 |
|
---|
120 | for(MxVertexID i=0; i<m->vert_count(); i++)
|
---|
121 | {
|
---|
122 | star.reset();
|
---|
123 | m->collect_vertex_star(i, star);
|
---|
124 |
|
---|
125 | for(uint j=0; j<star.length(); j++)
|
---|
126 | if( i < star(j) )
|
---|
127 | {
|
---|
128 | faces.reset();
|
---|
129 | m->collect_edge_neighbors(i, star(j), faces);
|
---|
130 | if( faces.length() == 1 )
|
---|
131 | discontinuity_constraint(i, star(j), faces);
|
---|
132 | }
|
---|
133 | }
|
---|
134 | }
|
---|
135 |
|
---|
136 |
|
---|
137 |
|
---|
138 |
|
---|
139 | MxEdgeQSlim::MxEdgeQSlim(MxStdModel& _m)
|
---|
140 | : MxQSlim(_m),
|
---|
141 | edge_links(_m.vert_count())
|
---|
142 | {
|
---|
143 | contraction_callback = NULL;
|
---|
144 | }
|
---|
145 |
|
---|
146 | MxEdgeQSlim::~MxEdgeQSlim()
|
---|
147 | {
|
---|
148 | // Delete everything remaining in the heap
|
---|
149 | for(uint i=0; i<heap.size(); i++)
|
---|
150 | delete ((MxQSlimEdge *)heap.item(i));
|
---|
151 | }
|
---|
152 |
|
---|
153 | ///////////////////////////////////////////////////////////////////////////
|
---|
154 | //
|
---|
155 | // IMPORTANT NOTE: These check_xxx functions assume that the local
|
---|
156 | // neighborhoods have been marked so that each face
|
---|
157 | // is marked with the number of involved vertices it has.
|
---|
158 | //
|
---|
159 |
|
---|
160 | double MxEdgeQSlim::check_local_compactness(uint v1, uint/*v2*/,
|
---|
161 | const float *vnew)
|
---|
162 | {
|
---|
163 | const MxFaceList& N1 = m->neighbors(v1);
|
---|
164 | double c_min = 1.0;
|
---|
165 |
|
---|
166 | for(uint i=0; i<N1.length(); i++)
|
---|
167 | if( m->face_mark(N1[i]) == 1 )
|
---|
168 | {
|
---|
169 | const MxFace& f = m->face(N1[i]);
|
---|
170 | Vec3 f_after[3];
|
---|
171 | for(uint j=0; j<3; j++)
|
---|
172 | f_after[j] = (f[j]==v1)?Vec3(vnew):Vec3(m->vertex(f[j]));
|
---|
173 |
|
---|
174 | double c=triangle_compactness(f_after[0], f_after[1], f_after[2]);
|
---|
175 |
|
---|
176 | if( c < c_min ) c_min = c;
|
---|
177 | }
|
---|
178 |
|
---|
179 | return c_min;
|
---|
180 | }
|
---|
181 |
|
---|
182 | double MxEdgeQSlim::check_local_inversion(uint v1,uint/*v2*/,const float *vnew)
|
---|
183 | {
|
---|
184 | double Nmin = 1.0;
|
---|
185 | const MxFaceList& N1 = m->neighbors(v1);
|
---|
186 |
|
---|
187 | for(uint i=0; i<N1.length(); i++)
|
---|
188 | if( m->face_mark(N1[i]) == 1 )
|
---|
189 | {
|
---|
190 | const MxFace& f = m->face(N1[i]);
|
---|
191 | Vec3 n_before;
|
---|
192 | m->compute_face_normal(N1[i], n_before);
|
---|
193 |
|
---|
194 | Vec3 f_after[3];
|
---|
195 | for(uint j=0; j<3; j++)
|
---|
196 | f_after[j] = (f[j]==v1)?Vec3(vnew):Vec3(m->vertex(f[j]));
|
---|
197 |
|
---|
198 | double delta = n_before *
|
---|
199 | triangle_normal(f_after[0], f_after[1], f_after[2]);
|
---|
200 |
|
---|
201 | if( delta < Nmin ) Nmin = delta;
|
---|
202 | }
|
---|
203 |
|
---|
204 | return Nmin;
|
---|
205 | }
|
---|
206 |
|
---|
207 | uint MxEdgeQSlim::check_local_validity(uint v1, uint /*v2*/, const float *vnew)
|
---|
208 |
|
---|
209 | {
|
---|
210 | const MxFaceList& N1 = m->neighbors(v1);
|
---|
211 | uint nfailed = 0;
|
---|
212 | uint i;
|
---|
213 |
|
---|
214 | for(i=0; i<N1.length(); i++)
|
---|
215 | if( m->face_mark(N1[i]) == 1 )
|
---|
216 | {
|
---|
217 | MxFace& f = m->face(N1[i]);
|
---|
218 | uint k = f.find_vertex(v1);
|
---|
219 | uint x = f[(k+1)%3];
|
---|
220 | uint y = f[(k+2)%3];
|
---|
221 |
|
---|
222 | float d_yx[3], d_vx[3], d_vnew[3], f_n[3], n[3];
|
---|
223 | mxv_sub(d_yx, m->vertex(y), m->vertex(x), 3); // d_yx = y-x
|
---|
224 | mxv_sub(d_vx, m->vertex(v1), m->vertex(x), 3); // d_vx = v-x
|
---|
225 | mxv_sub(d_vnew, vnew, m->vertex(x), 3); // d_vnew = vnew-x
|
---|
226 |
|
---|
227 | mxv_cross3(f_n, d_yx, d_vx);
|
---|
228 | mxv_cross3(n, f_n, d_yx); // n = ((y-x)^(v-x))^(y-x)
|
---|
229 | mxv_unitize(n, 3);
|
---|
230 |
|
---|
231 | // assert( mxv_dot(d_vx, n, 3) > -FEQ_EPS );
|
---|
232 | if(mxv_dot(d_vnew,n,3)<local_validity_threshold*mxv_dot(d_vx,n,3))
|
---|
233 | nfailed++;
|
---|
234 | }
|
---|
235 |
|
---|
236 | return nfailed;
|
---|
237 | }
|
---|
238 |
|
---|
239 | uint MxEdgeQSlim::check_local_degree(uint v1, uint v2, const float *)
|
---|
240 | {
|
---|
241 | const MxFaceList& N1 = m->neighbors(v1);
|
---|
242 | const MxFaceList& N2 = m->neighbors(v2);
|
---|
243 | uint i;
|
---|
244 | uint degree = 0;
|
---|
245 |
|
---|
246 | // Compute the degree of the vertex after contraction
|
---|
247 | //
|
---|
248 | for(i=0; i<N1.length(); i++)
|
---|
249 | if( m->face_mark(N1[i]) == 1 )
|
---|
250 | degree++;
|
---|
251 |
|
---|
252 | for(i=0; i<N2.length(); i++)
|
---|
253 | if( m->face_mark(N2[i]) == 1 )
|
---|
254 | degree++;
|
---|
255 |
|
---|
256 |
|
---|
257 | if( degree > vertex_degree_limit )
|
---|
258 | return degree - vertex_degree_limit;
|
---|
259 | else
|
---|
260 | return 0;
|
---|
261 | }
|
---|
262 |
|
---|
263 | void MxEdgeQSlim::apply_mesh_penalties(MxQSlimEdge *info)
|
---|
264 | {
|
---|
265 | uint i;
|
---|
266 |
|
---|
267 | const MxFaceList& N1 = m->neighbors(info->v1);
|
---|
268 | const MxFaceList& N2 = m->neighbors(info->v2);
|
---|
269 |
|
---|
270 | // Set up the face marks as the check_xxx() functions expect.
|
---|
271 | //
|
---|
272 | for(i=0; i<N2.length(); i++) m->face_mark(N2[i], 0);
|
---|
273 | for(i=0; i<N1.length(); i++) m->face_mark(N1[i], 1);
|
---|
274 | for(i=0; i<N2.length(); i++) m->face_mark(N2[i], m->face_mark(N2[i])+1);
|
---|
275 |
|
---|
276 | double base_error = info->heap_key();
|
---|
277 | double bias = 0.0;
|
---|
278 |
|
---|
279 | // Check for excess over degree bounds.
|
---|
280 | //
|
---|
281 | uint max_degree = MAX(N1.length(), N2.length());
|
---|
282 | if( max_degree > vertex_degree_limit )
|
---|
283 | bias += (max_degree-vertex_degree_limit) * meshing_penalty * 0.001;
|
---|
284 |
|
---|
285 | #if ALTERNATE_DEGREE_BIAS
|
---|
286 | // ??BUG: This code was supposed to be a slight improvement over
|
---|
287 | // the earlier version above. But it performs worse.
|
---|
288 | // Should check into why sometime.
|
---|
289 | //
|
---|
290 | uint degree_excess = check_local_degree(info->v1, info->v2, info->vnew);
|
---|
291 | if( degree_excess )
|
---|
292 | bias += degree_excess * meshing_penalty;
|
---|
293 | #endif
|
---|
294 |
|
---|
295 | // Local validity checks
|
---|
296 | //
|
---|
297 | uint nfailed = check_local_validity(info->v1, info->v2, info->vnew);
|
---|
298 | nfailed += check_local_validity(info->v2, info->v1, info->vnew);
|
---|
299 | if( nfailed )
|
---|
300 | bias += nfailed*meshing_penalty;
|
---|
301 |
|
---|
302 | if( compactness_ratio > 0.0 )
|
---|
303 | {
|
---|
304 | double c1_min=check_local_compactness(info->v1, info->v2, info->vnew);
|
---|
305 | double c2_min=check_local_compactness(info->v2, info->v1, info->vnew);
|
---|
306 | double c_min = MIN(c1_min, c2_min);
|
---|
307 |
|
---|
308 | // !!BUG: There's a small problem with this: it ignores the scale
|
---|
309 | // of the errors when adding the bias. For instance, enabling
|
---|
310 | // this on the cow produces bad results. I also tried
|
---|
311 | // += (base_error + FEQ_EPS) * (2-c_min), but that works
|
---|
312 | // poorly on the flat planar thing. A better solution is
|
---|
313 | // clearly needed.
|
---|
314 | //
|
---|
315 | // NOTE: The prior heuristic was
|
---|
316 | // if( ratio*cmin_before > cmin_after ) apply penalty;
|
---|
317 | //
|
---|
318 | if( c_min < compactness_ratio )
|
---|
319 | bias += (1-c_min);
|
---|
320 | }
|
---|
321 |
|
---|
322 | #if USE_OLD_INVERSION_CHECK
|
---|
323 | double Nmin1 = check_local_inversion(info->v1, info->v2, info->vnew);
|
---|
324 | double Nmin2 = check_local_inversion(info->v2, info->v1, info->vnew);
|
---|
325 | if( MIN(Nmin1, Nmin2) < 0.0 )
|
---|
326 | bias += meshing_penalty;
|
---|
327 | #endif
|
---|
328 |
|
---|
329 | info->heap_key(base_error - bias);
|
---|
330 | }
|
---|
331 |
|
---|
332 | void MxEdgeQSlim::compute_target_placement(MxQSlimEdge *info)
|
---|
333 | {
|
---|
334 | MxVertexID i=info->v1, j=info->v2;
|
---|
335 |
|
---|
336 | const Quadric &Qi=quadrics(i), &Qj=quadrics(j);
|
---|
337 |
|
---|
338 | Quadric Q = Qi; Q += Qj;
|
---|
339 | double e_min;
|
---|
340 |
|
---|
341 | if( placement_policy==MX_PLACE_OPTIMAL &&
|
---|
342 | Q.optimize(&info->vnew[X], &info->vnew[Y], &info->vnew[Z]) )
|
---|
343 | {
|
---|
344 | e_min = Q(info->vnew);
|
---|
345 | }
|
---|
346 | else
|
---|
347 | {
|
---|
348 | Vec3 vi(m->vertex(i)), vj(m->vertex(j));
|
---|
349 | Vec3 best;
|
---|
350 |
|
---|
351 | if( placement_policy>=MX_PLACE_LINE && Q.optimize(best, vi, vj) )
|
---|
352 | e_min = Q(best);
|
---|
353 | else
|
---|
354 | {
|
---|
355 | double ei=Q(vi), ej=Q(vj);
|
---|
356 |
|
---|
357 | if( ei < ej ) { e_min = ei; best = vi; }
|
---|
358 | else { e_min = ej; best = vj; }
|
---|
359 |
|
---|
360 | if( placement_policy>=MX_PLACE_ENDORMID )
|
---|
361 | {
|
---|
362 | Vec3 mid = (vi+vj)/2;
|
---|
363 | double e_mid = Q(mid);
|
---|
364 |
|
---|
365 | if( e_mid < e_min ) { e_min = e_mid; best = mid; }
|
---|
366 | }
|
---|
367 | }
|
---|
368 |
|
---|
369 | info->vnew[X] = best[X];
|
---|
370 | info->vnew[Y] = best[Y];
|
---|
371 | info->vnew[Z] = best[Z];
|
---|
372 | }
|
---|
373 |
|
---|
374 | if( weighting_policy == MX_WEIGHT_AREA_AVG )
|
---|
375 | e_min /= Q.area();
|
---|
376 |
|
---|
377 |
|
---|
378 | if ( weighting_policy == MX_WEIGHT_EDGE_ONLY )
|
---|
379 | {
|
---|
380 | Vec3 vi(m->vertex(i)), vj(m->vertex(j));
|
---|
381 | e_min = norm2(vi - vj);
|
---|
382 | }
|
---|
383 |
|
---|
384 |
|
---|
385 | info->heap_key(-e_min);
|
---|
386 | }
|
---|
387 |
|
---|
388 | /*
|
---|
389 | void MxEdgeQSlim::compute_target_placement(MxQSlimEdge *info)
|
---|
390 | {
|
---|
391 | MxVertexID i=info->v1, j=info->v2;
|
---|
392 |
|
---|
393 | const Quadric &Qi=quadrics(i), &Qj=quadrics(j);
|
---|
394 |
|
---|
395 | Quadric Q = Qi; Q += Qj;
|
---|
396 | double e_min;
|
---|
397 |
|
---|
398 | if( placement_policy==MX_PLACE_OPTIMAL &&
|
---|
399 | Q.optimize(&info->vnew[X], &info->vnew[Y], &info->vnew[Z]) )
|
---|
400 | {
|
---|
401 | e_min = Q(info->vnew);
|
---|
402 | }
|
---|
403 | else
|
---|
404 | {
|
---|
405 | Vec3 vi(m->vertex(i)), vj(m->vertex(j));
|
---|
406 | Vec3 best;
|
---|
407 |
|
---|
408 | if( placement_policy>=MX_PLACE_LINE && Q.optimize(best, vi, vj) )
|
---|
409 | e_min = Q(best);
|
---|
410 | else
|
---|
411 | {
|
---|
412 | double ei=Q(vi), ej=Q(vj);
|
---|
413 |
|
---|
414 | if( ei < ej ) { e_min = ei; best = vi; }
|
---|
415 | else { e_min = ej; best = vj; }
|
---|
416 |
|
---|
417 | if( placement_policy>=MX_PLACE_ENDORMID )
|
---|
418 | {
|
---|
419 | Vec3 mid = (vi+vj)/2;
|
---|
420 | double e_mid = Q(mid);
|
---|
421 |
|
---|
422 | if( e_mid < e_min ) { e_min = e_mid; best = mid; }
|
---|
423 | }
|
---|
424 | }
|
---|
425 |
|
---|
426 | info->vnew[X] = best[X];
|
---|
427 | info->vnew[Y] = best[Y];
|
---|
428 | info->vnew[Z] = best[Z];
|
---|
429 | }
|
---|
430 |
|
---|
431 | if( weighting_policy == MX_WEIGHT_AREA_AVG )
|
---|
432 | e_min /= Q.area();
|
---|
433 |
|
---|
434 | info->heap_key(-e_min);
|
---|
435 | }
|
---|
436 | */
|
---|
437 |
|
---|
438 | void MxEdgeQSlim::finalize_edge_update(MxQSlimEdge *info)
|
---|
439 | {
|
---|
440 | if( meshing_penalty > 1.0 )
|
---|
441 | apply_mesh_penalties(info);
|
---|
442 |
|
---|
443 | if( info->is_in_heap() )
|
---|
444 | heap.update(info);
|
---|
445 | else
|
---|
446 | heap.insert(info);
|
---|
447 | }
|
---|
448 |
|
---|
449 |
|
---|
450 | void MxEdgeQSlim::compute_edge_info(MxQSlimEdge *info)
|
---|
451 | {
|
---|
452 | compute_target_placement(info);
|
---|
453 |
|
---|
454 | finalize_edge_update(info);
|
---|
455 | }
|
---|
456 |
|
---|
457 | void MxEdgeQSlim::create_edge(MxVertexID i, MxVertexID j)
|
---|
458 | {
|
---|
459 | MxQSlimEdge *info = new MxQSlimEdge;
|
---|
460 |
|
---|
461 | edge_links(i).add(info);
|
---|
462 | edge_links(j).add(info);
|
---|
463 |
|
---|
464 | info->v1 = i;
|
---|
465 | info->v2 = j;
|
---|
466 |
|
---|
467 | compute_edge_info(info);
|
---|
468 | }
|
---|
469 |
|
---|
470 | void MxEdgeQSlim::collect_edges()
|
---|
471 | {
|
---|
472 | MxVertexList star;
|
---|
473 |
|
---|
474 | for(MxVertexID i=0; i<m->vert_count(); i++)
|
---|
475 | {
|
---|
476 | star.reset();
|
---|
477 | m->collect_vertex_star(i, star);
|
---|
478 |
|
---|
479 | for(uint j=0; j<star.length(); j++)
|
---|
480 | if( i < star(j) ) // Only add particular edge once
|
---|
481 | create_edge(i, star(j));
|
---|
482 | }
|
---|
483 | }
|
---|
484 |
|
---|
485 | void MxEdgeQSlim::initialize()
|
---|
486 | {
|
---|
487 | MxQSlim::initialize();
|
---|
488 | collect_edges();
|
---|
489 | }
|
---|
490 |
|
---|
491 | void MxEdgeQSlim::initialize(const MxEdge *edges, uint count)
|
---|
492 | {
|
---|
493 | MxQSlim::initialize();
|
---|
494 | for(uint i=0; i<count; i++)
|
---|
495 | create_edge(edges[i].v1, edges[i].v2);
|
---|
496 | }
|
---|
497 |
|
---|
498 | void MxEdgeQSlim::update_pre_contract(const MxPairContraction& conx)
|
---|
499 | {
|
---|
500 | MxVertexID v1=conx.v1, v2=conx.v2;
|
---|
501 | uint i, j;
|
---|
502 |
|
---|
503 | star.reset();
|
---|
504 | //
|
---|
505 | // Before, I was gathering the vertex "star" using:
|
---|
506 | // m->collect_vertex_star(v1, star);
|
---|
507 | // This doesn't work when we initially begin with a subset of
|
---|
508 | // the total edges. Instead, we need to collect the "star"
|
---|
509 | // from the edge links maintained at v1.
|
---|
510 | //
|
---|
511 | for(i=0; i<edge_links(v1).length(); i++)
|
---|
512 | star.add(edge_links(v1)[i]->opposite_vertex(v1));
|
---|
513 |
|
---|
514 | for(i=0; i<edge_links(v2).length(); i++)
|
---|
515 | {
|
---|
516 | MxQSlimEdge *e = edge_links(v2)(i);
|
---|
517 | MxVertexID u = (e->v1==v2)?e->v2:e->v1;
|
---|
518 | SanityCheck( e->v1==v2 || e->v2==v2 );
|
---|
519 | SanityCheck( u!=v2 );
|
---|
520 |
|
---|
521 | if( u==v1 || star.find(u) )
|
---|
522 | {
|
---|
523 | // This is a useless link --- kill it
|
---|
524 | bool found = edge_links(u).find(e, &j);
|
---|
525 | assert( found );
|
---|
526 | edge_links(u).remove(j);
|
---|
527 | heap.remove(e);
|
---|
528 | if( u!=v1 ) delete e; // (v1,v2) will be deleted later
|
---|
529 | }
|
---|
530 | else
|
---|
531 | {
|
---|
532 | // Relink this to v1
|
---|
533 | e->v1 = v1;
|
---|
534 | e->v2 = u;
|
---|
535 | edge_links(v1).add(e);
|
---|
536 | }
|
---|
537 | }
|
---|
538 |
|
---|
539 | edge_links(v2).reset();
|
---|
540 | }
|
---|
541 |
|
---|
542 | void MxEdgeQSlim::update_post_contract(const MxPairContraction& conx)
|
---|
543 | {
|
---|
544 | }
|
---|
545 |
|
---|
546 | void MxEdgeQSlim::apply_contraction(const MxPairContraction& conx)
|
---|
547 | {
|
---|
548 | //
|
---|
549 | // Pre-contraction update
|
---|
550 | valid_verts--;
|
---|
551 | valid_faces -= conx.dead_faces.length();
|
---|
552 | quadrics(conx.v1) += quadrics(conx.v2);
|
---|
553 |
|
---|
554 | update_pre_contract(conx);
|
---|
555 |
|
---|
556 | m->apply_contraction(conx);
|
---|
557 |
|
---|
558 | update_post_contract(conx);
|
---|
559 |
|
---|
560 | // Must update edge info here so that the meshing penalties
|
---|
561 | // will be computed with respect to the new mesh rather than the old
|
---|
562 | for(uint i=0; i<edge_links(conx.v1).length(); i++)
|
---|
563 | compute_edge_info(edge_links(conx.v1)[i]);
|
---|
564 | }
|
---|
565 |
|
---|
566 | void MxEdgeQSlim::update_pre_expand(const MxPairContraction&)
|
---|
567 | {
|
---|
568 | }
|
---|
569 |
|
---|
570 | void MxEdgeQSlim::update_post_expand(const MxPairContraction& conx)
|
---|
571 | {
|
---|
572 | MxVertexID v1=conx.v1, v2=conx.v2;
|
---|
573 | uint i;
|
---|
574 |
|
---|
575 | star.reset(); star2.reset();
|
---|
576 | PRECAUTION(edge_links(conx.v2).reset());
|
---|
577 | m->collect_vertex_star(conx.v1, star);
|
---|
578 | m->collect_vertex_star(conx.v2, star2);
|
---|
579 |
|
---|
580 | i = 0;
|
---|
581 | while( i<edge_links(v1).length() )
|
---|
582 | {
|
---|
583 | MxQSlimEdge *e = edge_links(v1)(i);
|
---|
584 | MxVertexID u = (e->v1==v1)?e->v2:e->v1;
|
---|
585 | SanityCheck( e->v1==v1 || e->v2==v1 );
|
---|
586 | SanityCheck( u!=v1 && u!=v2 );
|
---|
587 |
|
---|
588 | bool v1_linked = star.find(u);
|
---|
589 | bool v2_linked = star2.find(u);
|
---|
590 |
|
---|
591 | if( v1_linked )
|
---|
592 | {
|
---|
593 | if( v2_linked ) create_edge(v2, u);
|
---|
594 | i++;
|
---|
595 | }
|
---|
596 | else
|
---|
597 | {
|
---|
598 | // !! BUG: I expected this to be true, but it's not.
|
---|
599 | // Need to find out why, and whether it's my
|
---|
600 | // expectation or the code that's wrong.
|
---|
601 | // SanityCheck(v2_linked);
|
---|
602 | e->v1 = v2;
|
---|
603 | e->v2 = u;
|
---|
604 | edge_links(v2).add(e);
|
---|
605 | edge_links(v1).remove(i);
|
---|
606 | }
|
---|
607 |
|
---|
608 | compute_edge_info(e);
|
---|
609 | }
|
---|
610 |
|
---|
611 | if( star.find(v2) )
|
---|
612 | // ?? BUG: Is it legitimate for there not to be an edge here ??
|
---|
613 | create_edge(v1, v2);
|
---|
614 | }
|
---|
615 |
|
---|
616 |
|
---|
617 | void MxEdgeQSlim::apply_expansion(const MxPairContraction& conx)
|
---|
618 | {
|
---|
619 | update_pre_expand(conx);
|
---|
620 |
|
---|
621 | m->apply_expansion(conx);
|
---|
622 |
|
---|
623 | //
|
---|
624 | // Post-expansion update
|
---|
625 | valid_verts++;
|
---|
626 | valid_faces += conx.dead_faces.length();
|
---|
627 | quadrics(conx.v1) -= quadrics(conx.v2);
|
---|
628 |
|
---|
629 | update_post_expand(conx);
|
---|
630 | }
|
---|
631 |
|
---|
632 | bool MxEdgeQSlim::decimate(uint target)
|
---|
633 | {
|
---|
634 | MxPairContraction local_conx;
|
---|
635 |
|
---|
636 | while( valid_faces > target )
|
---|
637 | {
|
---|
638 | MxQSlimEdge *info = (MxQSlimEdge *)heap.extract();
|
---|
639 | if( !info ) { return false; }
|
---|
640 |
|
---|
641 | MxVertexID v1=info->v1, v2=info->v2;
|
---|
642 |
|
---|
643 | if( m->vertex_is_valid(v1) && m->vertex_is_valid(v2) )
|
---|
644 | {
|
---|
645 | MxPairContraction& conx = local_conx;
|
---|
646 |
|
---|
647 | m->compute_contraction(v1, v2, &conx, info->vnew);
|
---|
648 |
|
---|
649 | if( will_join_only && conx.dead_faces.length()>0 ) continue;
|
---|
650 |
|
---|
651 | if( contraction_callback )
|
---|
652 | (*contraction_callback)(conx, -info->heap_key());
|
---|
653 |
|
---|
654 | apply_contraction(conx);
|
---|
655 | }
|
---|
656 |
|
---|
657 | delete info;
|
---|
658 | }
|
---|
659 |
|
---|
660 | return true;
|
---|
661 | }
|
---|
662 |
|
---|
663 |
|
---|
664 |
|
---|
665 | void MxFaceQSlim::compute_face_info(MxFaceID f)
|
---|
666 | {
|
---|
667 | tri_info& info = f_info(f);
|
---|
668 | info.f = f;
|
---|
669 |
|
---|
670 | MxVertexID i = m->face(f)(0);
|
---|
671 | MxVertexID j = m->face(f)(1);
|
---|
672 | MxVertexID k = m->face(f)(2);
|
---|
673 |
|
---|
674 | const Quadric& Qi = quadrics(i);
|
---|
675 | const Quadric& Qj = quadrics(j);
|
---|
676 | const Quadric& Qk = quadrics(k);
|
---|
677 |
|
---|
678 | Quadric Q = Qi;
|
---|
679 | Q += Qj;
|
---|
680 | Q += Qk;
|
---|
681 |
|
---|
682 | if( placement_policy == MX_PLACE_OPTIMAL &&
|
---|
683 | Q.optimize(&info.vnew[X], &info.vnew[Y], &info.vnew[Z]) )
|
---|
684 | {
|
---|
685 | info.heap_key(-Q(info.vnew));
|
---|
686 | }
|
---|
687 | else
|
---|
688 | {
|
---|
689 | Vec3 v1(m->vertex(i)), v2(m->vertex(j)), v3(m->vertex(k));
|
---|
690 | double e1=Q(v1), e2=Q(v2), e3=Q(v3);
|
---|
691 |
|
---|
692 | Vec3 best;
|
---|
693 | double e_min;
|
---|
694 |
|
---|
695 | if( e1<=e2 && e1<=e3 ) { e_min=e1; best=v1; }
|
---|
696 | else if( e2<=e1 && e2<=e3 ) { e_min=e2; best=v2; }
|
---|
697 | else { e_min=e3; best=v3; }
|
---|
698 |
|
---|
699 | info.vnew[X] = best[X];
|
---|
700 | info.vnew[Y] = best[Y];
|
---|
701 | info.vnew[Z] = best[Z];
|
---|
702 | info.heap_key(-e_min);
|
---|
703 | }
|
---|
704 |
|
---|
705 | if( weighting_policy == MX_WEIGHT_AREA_AVG )
|
---|
706 | info.heap_key(info.heap_key() / Q.area());
|
---|
707 |
|
---|
708 | if( info.is_in_heap() )
|
---|
709 | heap.update(&info);
|
---|
710 | else
|
---|
711 | heap.insert(&info);
|
---|
712 | }
|
---|
713 |
|
---|
714 |
|
---|
715 | MxFaceQSlim::MxFaceQSlim(MxStdModel& _m)
|
---|
716 | : MxQSlim(_m), f_info(_m.face_count())
|
---|
717 | {
|
---|
718 | }
|
---|
719 |
|
---|
720 | void MxFaceQSlim::initialize()
|
---|
721 | {
|
---|
722 | MxQSlim::initialize();
|
---|
723 |
|
---|
724 | for(MxFaceID f=0; f<m->face_count(); f++)
|
---|
725 | compute_face_info(f);
|
---|
726 | }
|
---|
727 |
|
---|
728 | bool MxFaceQSlim::decimate(uint target)
|
---|
729 | {
|
---|
730 | unsigned int i;
|
---|
731 |
|
---|
732 | MxFaceList changed;
|
---|
733 |
|
---|
734 | while( valid_faces > target )
|
---|
735 | {
|
---|
736 | tri_info *info = (tri_info *)heap.extract();
|
---|
737 | if( !info ) { return false; }
|
---|
738 |
|
---|
739 | MxFaceID f = info->f;
|
---|
740 | MxVertexID v1 = m->face(f)(0),
|
---|
741 | v2 = m->face(f)(1),
|
---|
742 | v3 = m->face(f)(2);
|
---|
743 |
|
---|
744 | if( m->face_is_valid(f) )
|
---|
745 | {
|
---|
746 | //
|
---|
747 | // Perform the actual contractions
|
---|
748 | m->contract(v1, v2, v3, info->vnew, changed);
|
---|
749 |
|
---|
750 | quadrics(v1) += quadrics(v2); // update quadric of v1
|
---|
751 | quadrics(v1) += quadrics(v3);
|
---|
752 |
|
---|
753 | //
|
---|
754 | // Update valid counts
|
---|
755 | valid_verts -= 2;
|
---|
756 | for(i=0; i<changed.length(); i++)
|
---|
757 | if( !m->face_is_valid(changed(i)) ) valid_faces--;
|
---|
758 |
|
---|
759 | for(i=0; i<changed.length(); i++)
|
---|
760 | if( m->face_is_valid(changed(i)) )
|
---|
761 | compute_face_info(changed(i));
|
---|
762 | else
|
---|
763 | heap.remove(&f_info(changed(i)));
|
---|
764 | }
|
---|
765 | }
|
---|
766 |
|
---|
767 | return true;
|
---|
768 | }
|
---|