Subversion Repositories gelsvn

Rev

Rev 182 | Rev 219 | Go to most recent revision | Only display areas with differences | Ignore whitespace | Details | Blame | Last modification | View Log | RSS feed

Rev 182 Rev 215
1
#include <queue>
1
#include <queue>
2
#include <iostream>
2
#include <iostream>
3
#include "quadric_simplify.h"
3
#include "quadric_simplify.h"
4
#include "smooth.h"
4
#include "smooth.h"
-
 
5
#include "Util/BinaryHeap.h"
5
#include "HMesh/VertexCirculator.h"
6
#include "HMesh/VertexCirculator.h"
6
#include "Geometry/QEM.h"
7
#include "Geometry/QEM.h"
7
 
8
 
8
 
9
 
9
namespace HMesh
10
namespace HMesh
10
{
11
{
11
  using namespace CGLA;
12
		using namespace CGLA;
12
  using namespace std;
13
		using namespace std;
13
	using namespace GEO;
14
		using namespace GEO;
14
	using namespace HMesh;
15
		using namespace HMesh;
15
 
-
 
16
	namespace
-
 
17
	{
-
 
18
		struct HalfEdgeRec
-
 
19
		{
-
 
20
			HalfEdgeIter he;
-
 
21
			int time_stamp;
-
 
22
			void halfedge_removed() {time_stamp = -1;}
-
 
23
			HalfEdgeRec(): time_stamp(0) {}
-
 
24
		};
-
 
25
 
16
 
26
		struct SimplifyRec
17
		namespace
27
		{
18
		{
-
 
19
				/* We create a record for each halfedge where we can keep its time
-
 
20
					 stamp. If the time stamp on the halfedge record is bigger than
-
 
21
					 the stamp on the simplification record, we cannot use the 
-
 
22
					 simplification record (see below). */
-
 
23
				struct HalfEdgeRec
-
 
24
				{
28
			Vec3f opt_pos;
25
						HalfEdgeIter he;
29
			int he_index;
26
						int time_stamp;
-
 
27
						void halfedge_removed() {time_stamp = -1;}
-
 
28
						HalfEdgeRec(): time_stamp(0) {}
30
			float err;
29
				};
-
 
30
 
-
 
31
				/* The simpliciation record contains information about a potential
-
 
32
					 edge contraction */
31
			int time_stamp;
33
				struct SimplifyRec
-
 
34
				{
-
 
35
						Vec3d opt_pos;  // optimal vertex position
-
 
36
						int he_index;   // Index (into HalfEdgeRec vector) of edge
-
 
37
						                // we want to contract
-
 
38
						float err;      // Error associated with contraction
-
 
39
						int time_stamp; // Time stamp (see comment on HalfEdgeRec)
-
 
40
						int visits;     // Visits (number of times we considered this 
-
 
41
						                // record).
32
		};
42
				};
33
 
43
 
34
		bool operator<(const SimplifyRec& s1, const SimplifyRec& s2)
44
				bool operator<(const SimplifyRec& s1, const SimplifyRec& s2)
35
		{
45
				{
36
			return s1.err > s2.err;
46
						return s1.err > s2.err;
37
		}
47
				}
38
 
48
 
39
		class QuadricSimplifier
49
				class QuadricSimplifier
40
		{
50
				{
41
			typedef priority_queue<SimplifyRec> SimplifyQueue;
51
						typedef priority_queue<SimplifyRec> SimplifyQueue;
42
			typedef vector<QEM> QEMVec;
52
						typedef vector<QEM> QEMVec;
43
			typedef vector<HalfEdgeRec> HalfEdgeVec;
53
						typedef vector<HalfEdgeRec> HalfEdgeVec;
44
 
54
 
45
			Manifold& m;
55
						Manifold& m;
46
			HalfEdgeVec halfedge_vec;
56
						HalfEdgeVec halfedge_vec;
47
			QEMVec qem_vec;
57
						QEMVec qem_vec;
48
			SimplifyQueue sim_queue;
58
						SimplifyQueue sim_queue;
49
 
59
 
-
 
60
						/* Compute the error associated with contraction of he and the
-
 
61
							 optimal position of resulting vertex. */
50
			void push_simplify_rec(HalfEdgeIter he);
62
						void push_simplify_rec(HalfEdgeIter he);
-
 
63
 
-
 
64
						/* Check whether the contraction is valid. See below for details*/
-
 
65
						bool check_consistency(HalfEdgeIter he, const Vec3d& opt_pos);
-
 
66
 
-
 
67
						/* Update the time stamp of a halfedge. A halfedge and its opp edge
-
 
68
							 may have different stamps. We choose a stamp that is greater
-
 
69
							 than either and assign to both.*/
-
 
70
						void update_time_stamp(HalfEdgeIter he)
-
 
71
								{
-
 
72
										int time_stamp = s_max(halfedge_vec[he->touched].time_stamp,
-
 
73
																					 halfedge_vec[he->opp->touched].time_stamp);
-
 
74
										time_stamp++;
-
 
75
										halfedge_vec[he->touched].time_stamp = time_stamp;
-
 
76
										halfedge_vec[he->opp->touched].time_stamp = time_stamp;
-
 
77
								}
51
 
78
 
-
 
79
						/* Update time stamps for all halfedges in one ring of vi */
52
			void put_onering_in_queue(VertexIter vi);
80
						void update_onering_timestamp(VertexIter vi);
53
 
81
 
-
 
82
						/* Perform a collapse - if conditions are met. Returns 1 or 0 
-
 
83
							 accordingly. */
54
			int collapse(const SimplifyRec& simplify_rec);
84
						int collapse(SimplifyRec& simplify_rec);
55
 
85
 
56
		public:
86
				public:
57
 
87
 
-
 
88
						/* Create a simplifier for a manifold */
58
			QuadricSimplifier(Manifold& _m): 
89
						QuadricSimplifier(Manifold& _m): 
59
				m(_m), 
90
								m(_m), 
60
				halfedge_vec(m.no_halfedges()), 
91
								halfedge_vec(m.no_halfedges()), 
61
				qem_vec(m.no_vertices()) 
92
								qem_vec(m.no_vertices()) 
62
			{
93
								{
63
				// Enumerate vertices
94
										// Enumerate vertices
64
				m.enumerate_vertices();
95
										m.enumerate_vertices();
65
			
96
			
66
				// Enumerate halfedges
97
										// Enumerate halfedges
67
				m.enumerate_halfedges();
98
										m.enumerate_halfedges();
68
			}
99
								}
69
	
100
	
-
 
101
						/* Simplify doing at most max_work contractions */
70
			void reduce(int max_work);
102
						void reduce(int max_work);
71
		};
103
				};
72
 
104
 
73
		void QuadricSimplifier::push_simplify_rec(HalfEdgeIter he)
105
				bool QuadricSimplifier::check_consistency(HalfEdgeIter he, 
74
		{
-
 
75
			// Update time_stamp
-
 
76
			assert(halfedge_vec[he->touched].time_stamp >=0);
-
 
77
			assert(halfedge_vec[he->opp->touched].time_stamp>=0);
-
 
78
			int time_stamp = s_max(halfedge_vec[he->touched].time_stamp,
-
 
79
														 halfedge_vec[he->opp->touched].time_stamp);
106
																									const Vec3d& opt_pos)
80
			time_stamp++;
107
				{
81
			halfedge_vec[he->touched].time_stamp = time_stamp;
-
 
82
			halfedge_vec[he->opp->touched].time_stamp = time_stamp;
-
 
83
			// Get QEM for both end points
108
						VertexIter v0 = he->vert;
84
			const QEM& Q1 = qem_vec[he->vert->touched];
-
 
85
			const QEM& Q2 = qem_vec[he->opp->vert->touched];
109
						VertexIter v1 = he->opp->vert;
86
 
-
 
87
			QEM q = Q1;
-
 
88
			q += Q2;
-
 
89
 
-
 
90
			float err;
-
 
91
			Vec3f opt_pos(0);
110
						Vec3d p0(v0->pos);
92
 
111
 
93
			//cout << q.determinant() << endl;
-
 
94
			if (q.determinant()>= 1e-8)
-
 
95
				{
-
 
96
					opt_pos = q.opt_pos(0.001);
-
 
97
					err = q.error(opt_pos);
-
 
98
				}
-
 
99
			else
-
 
100
				{
-
 
101
					// To see the effectivenes of the scheme, we can create 
112
						/* This test is inspired by Garland's Ph.D. thesis. We try
102
					// Random errors for comparison. Uncomment to use:
113
							 to detect whether flipped triangles will occur by sort of
103
					// 			float e0 = time_stamp + rand()/float(RAND_MAX);//q.error(v0);
114
							 ensuring that the new vertex is in the hull of the one rings
104
					// 			float e1 = time_stamp + rand()/float(RAND_MAX);//q.error(v1);
115
							 of the vertices at either end of the edge being contracted */
105
 
-
 
106
					Vec3f v0 = he->vert->pos;
-
 
107
					Vec3f v1 = he->opp->vert->pos;
-
 
108
					float e0 = q.error(v0);
-
 
109
					float e1 = q.error(v1);
-
 
110
		
-
 
111
					if(e0 < e1)
-
 
112
						{
116
						
113
							err = e0;
-
 
114
							opt_pos = v0;
117
						for(VertexCirculator vc(v0); !vc.end(); ++vc)
115
						}
-
 
116
					else
-
 
117
						{
118
						{
-
 
119
								HalfEdgeIter h = vc.get_halfedge();
-
 
120
								if(h->vert != v1 && h->next->vert != v1)
-
 
121
								{
-
 
122
										Vec3d pa(h->vert->pos);
-
 
123
										Vec3d pb(h->next->vert->pos);
118
							err = e1;
124
										
-
 
125
										Vec3d dir = normalize(pb-pa);
-
 
126
 
-
 
127
										Vec3d n = p0 - pa;
-
 
128
										n = n - dir * dot(dir,n);
-
 
129
 
-
 
130
										if(dot(n,opt_pos - pa) <= 0)
119
							opt_pos = v1;
131
												return false;
-
 
132
								}
120
						}
133
						}
121
				}			
-
 
122
 
134
 
123
			// Create SimplifyRec
135
						/* Since the test above is not really enough, we also make sure 
124
			int he_index = he->touched;
136
							 that we do not introduce vertices of very high (>=10) or low 
125
			SimplifyRec simplify_rec;
137
							 (<=3) valency. Flipped faces seem to be always associated with
126
			simplify_rec.opt_pos = opt_pos;
138
							 very irregular valencies. It seems to get rid of most problems
127
			simplify_rec.err = err;
139
							 not fixed by the check above. */
-
 
140
							 
128
			simplify_rec.he_index = he_index;
141
						if(valency(he->next->vert)<5)
129
			simplify_rec.time_stamp = time_stamp;
142
								return false;
130
 
143
 
131
			// push it.
144
						if((valency(he->vert) + valency(he->opp->vert) ) > 12)
132
			sim_queue.push(simplify_rec);
145
								return false;
133
		}
-
 
134
		
-
 
135
 
146
 
136
		void QuadricSimplifier::put_onering_in_queue(VertexIter vi)
-
 
137
		{
-
 
138
			// For all emanating edges he
-
 
139
			for(VertexCirculator vc(vi);
-
 
140
					!vc.end(); ++vc)
147
						return true;
141
				{
-
 
142
					HalfEdgeIter he = vc.get_halfedge();
-
 
143
					push_simplify_rec(he);
-
 
144
				}
148
				}
145
		}
-
 
146
 
149
 
147
		int QuadricSimplifier::collapse(const SimplifyRec& simplify_rec)
-
 
148
		{
-
 
149
			int he_index = simplify_rec.he_index;
-
 
150
 
150
 
151
			// Check the time stamp to verify that the simplification 
151
				void QuadricSimplifier::push_simplify_rec(HalfEdgeIter he)
-
 
152
				{
152
			// record is the newest. If the halfedge has been removed
153
						// Get QEM for both end points
153
			// the time stamp is -1 and the comparison will also fail.
154
						const QEM& Q1 = qem_vec[he->vert->touched];
154
			if(halfedge_vec[he_index].time_stamp == simplify_rec.time_stamp)
155
						const QEM& Q2 = qem_vec[he->opp->vert->touched];
-
 
156
 
155
				{
157
						QEM q = Q1;
156
					HalfEdgeIter he = halfedge_vec[he_index].he;
158
						q += Q2;
-
 
159
 
157
					VertexIter v = he->opp->vert;
160
						float err;
158
					VertexIter n = he->vert;
161
						Vec3d opt_pos(0);
159
 
162
 
160
					// If the edge is, in fact, collapsible
163
						//cout << q.determinant() << endl;
161
					if(m.collapse_precond(he))
164
						if (q.determinant()>= 1e-5)
162
						{
165
						{
163
							//cout << simplify_rec.err << " " << &(*he->vert) << endl;
166
								opt_pos = q.opt_pos(1e-5);
164
							// Get QEM for both end points
167
								err = q.error(opt_pos);
-
 
168
						}
-
 
169
						else
-
 
170
						{
165
							const QEM& Q1 = qem_vec[he->vert->touched];
171
								// To see the effectivenes of the scheme, we can create 
166
							const QEM& Q2 = qem_vec[he->opp->vert->touched];
172
								// Random errors for comparison. Uncomment to use:
167
 
-
 
168
							// Compute Q_new = Q_1 + Q_2
173
								// 			float e0 = time_stamp + rand()/float(RAND_MAX);//q.error(v0);
169
							QEM q = Q1;
174
								// 			float e1 = time_stamp + rand()/float(RAND_MAX);//q.error(v1);
170
							q += Q2;
-
 
171
 
175
 
172
							// Mark all halfedges that will be removed as dead
176
								Vec3d v0(he->vert->pos);
173
							halfedge_vec[he->touched].halfedge_removed();
177
								Vec3d v1(he->opp->vert->pos);
-
 
178
								float e0 = q.error(v0);
174
							halfedge_vec[he->opp->touched].halfedge_removed();
179
								float e1 = q.error(v1);
175
						
180
		
176
							if(he->next->next->next == he)
181
								if(e0 < e1)
177
								{
182
								{
178
									halfedge_vec[he->next->touched].halfedge_removed();
183
										err = e0;
179
									halfedge_vec[he->next->next->touched].halfedge_removed();
184
										opt_pos = v0;
180
								}
185
								}
181
							if(he->opp->next->next->next == he->opp)
186
								else
182
								{
187
								{
183
									halfedge_vec[he->opp->next->touched].halfedge_removed();
188
										err = e1;
184
									halfedge_vec[he->opp->next->next->touched].halfedge_removed();
189
										opt_pos = v1;
185
								}
190
								}
-
 
191
						}			
186
 
192
 
187
							// Do collapse
193
						// Create SimplifyRec
188
							m.collapse_halfedge(he,false);
194
						int he_index = he->touched;
-
 
195
						SimplifyRec simplify_rec;
189
							n->pos = simplify_rec.opt_pos;
196
						simplify_rec.opt_pos = opt_pos;
-
 
197
						simplify_rec.err = err;
-
 
198
						simplify_rec.he_index = he_index;
-
 
199
						simplify_rec.time_stamp = halfedge_vec[he->touched].time_stamp;
190
							qem_vec[n->touched] = q;
200
						simplify_rec.visits = 0;
-
 
201
						// push it.
-
 
202
						sim_queue.push(simplify_rec);
191
						
203
						
192
							put_onering_in_queue(n);
-
 
193
							return 1;
-
 
194
						}
-
 
195
				}
204
				}
196
			return 0;
-
 
197
		}
205
		
198
 
206
 
199
		void QuadricSimplifier::reduce(int max_work)
207
				void QuadricSimplifier::update_onering_timestamp(VertexIter vi)
200
		{
208
				{
-
 
209
						// For all emanating edges he
-
 
210
						for(VertexCirculator vc(vi);
-
 
211
								!vc.end(); ++vc)
-
 
212
						{
-
 
213
								HalfEdgeIter he = vc.get_halfedge();
-
 
214
								update_time_stamp(he);
-
 
215
						}
-
 
216
				}
201
 
217
 
202
			// Set t = 0 for all halfedges
-
 
203
			int i=0;
-
 
204
			for(HalfEdgeIter he = m.halfedges_begin();
-
 
205
					he != m.halfedges_end(); ++he, ++i)
-
 
206
				halfedge_vec[i].he = he;
-
 
207
		
-
 
208
			// For all vertices, compute quadric and store in qem_vec
218
				int QuadricSimplifier::collapse(SimplifyRec& simplify_rec)
209
			for(VertexIter vi=m.vertices_begin(); vi != m.vertices_end(); ++vi)
-
 
210
				{
219
				{
211
					Vec3d p(vi->pos);
220
						int he_index = simplify_rec.he_index;
212
					QEM q;
221
						HalfEdgeIter he = halfedge_vec[he_index].he;
-
 
222
 
213
					for(VertexCirculator vc(vi); !vc.end(); ++vc)
223
						// Check the time stamp to verify that the simplification 
-
 
224
						// record is the newest. If the halfedge has been removed
-
 
225
						// the time stamp is -1 and the comparison will also fail.
-
 
226
						if(halfedge_vec[he_index].time_stamp == simplify_rec.time_stamp)
214
						{
227
						{
215
							FaceIter f = vc.get_face();
228
								VertexIter v = he->opp->vert;
-
 
229
								VertexIter n = he->vert;
-
 
230
 
-
 
231
								// If the edge is, in fact, collapsible
216
							if(f != NULL_FACE_ITER)
232
								if(m.collapse_precond(he))
217
								{
233
								{
-
 
234
										// If our consistency checks pass, we are relatively
-
 
235
										// sure that the contraction does not lead to a face flip.
-
 
236
										if(check_consistency(he, simplify_rec.opt_pos) && 
-
 
237
											 check_consistency(he->opp, simplify_rec.opt_pos))
-
 
238
										{
-
 
239
 
-
 
240
												//cout << simplify_rec.err << " " << &(*he->vert) << endl;
-
 
241
												// Get QEM for both end points
-
 
242
												const QEM& Q1 = qem_vec[he->vert->touched];
-
 
243
												const QEM& Q2 = qem_vec[he->opp->vert->touched];
-
 
244
												
-
 
245
												// Compute Q_new = Q_1 + Q_2
-
 
246
												QEM q = Q1;
-
 
247
												q += Q2;
-
 
248
												
-
 
249
												// Mark all halfedges that will be removed as dead
-
 
250
												halfedge_vec[he->touched].halfedge_removed();
-
 
251
												halfedge_vec[he->opp->touched].halfedge_removed();
-
 
252
						
-
 
253
												if(he->next->next->next == he)
-
 
254
												{
-
 
255
														halfedge_vec[he->next->touched].halfedge_removed();
-
 
256
														halfedge_vec[he->next->next->touched].halfedge_removed();
-
 
257
												}
-
 
258
												if(he->opp->next->next->next == he->opp)
-
 
259
												{
-
 
260
														halfedge_vec[he->opp->next->touched].halfedge_removed();
-
 
261
														halfedge_vec[he->opp->next->next->touched].halfedge_removed();
-
 
262
												}
-
 
263
												
218
									Vec3d n(normal(f));
264
												// Do collapse
-
 
265
												m.collapse_halfedge(he,false);
-
 
266
												n->pos = Vec3f(simplify_rec.opt_pos);
-
 
267
												qem_vec[n->touched] = q;
-
 
268
												
-
 
269
												update_onering_timestamp(n);
219
									double a = area(f);
270
												return 1;
-
 
271
										}
-
 
272
								}
-
 
273
								// If we are here, the collapse was not allowed. If we have
-
 
274
								// seen this simplify record less than 100 times, we try to
-
 
275
								// increase the error and store the record again. Maybe some
-
 
276
								// other contractions will make it more digestible later.
-
 
277
								if(simplify_rec.visits < 10)
-
 
278
								{
-
 
279
										simplify_rec.err += 5*simplify_rec.err + 1e-10;
220
									q += QEM(p, n, a / 3.0);
280
										++simplify_rec.visits;
-
 
281
										sim_queue.push(simplify_rec);
221
								}
282
								}
222
						}
283
						}
-
 
284
						// Ok the time stamp did not match. Create a new record.
-
 
285
						else if(halfedge_vec[he_index].time_stamp != -1)
223
					qem_vec[vi->touched] = q;
286
								push_simplify_rec(he);
-
 
287
 
-
 
288
						return 0;
224
				}
289
				}
225
 
290
 
226
			for(HalfEdgeIter he = m.halfedges_begin();
-
 
227
					he != m.halfedges_end(); ++he)
291
				void QuadricSimplifier::reduce(int max_work)
228
				// For all halfedges, 
-
 
229
				{
292
				{
-
 
293
						// Set t = 0 for all halfedges
-
 
294
						int i=0;
230
					if(halfedge_vec[he->touched].time_stamp == 0)
295
						for(HalfEdgeIter he = m.halfedges_begin();
231
						push_simplify_rec(he);
296
								he != m.halfedges_end(); ++he, ++i)
232
				}
297
								halfedge_vec[i].he = he;
233
 
298
 
-
 
299
						cout << "Computing quadrics" << endl;
-
 
300
		
-
 
301
						// For all vertices, compute quadric and store in qem_vec
-
 
302
						for(VertexIter vi=m.vertices_begin(); vi != m.vertices_end(); ++vi)
-
 
303
						{
-
 
304
								Vec3d p(vi->pos);
-
 
305
								QEM q;
-
 
306
								for(VertexCirculator vc(vi); !vc.end(); ++vc)
-
 
307
								{
-
 
308
										FaceIter f = vc.get_face();
-
 
309
										if(f != NULL_FACE_ITER)
-
 
310
										{
-
 
311
												Vec3d n(normal(f));
-
 
312
												double a = area(f);
-
 
313
												q += QEM(p, n, a / 3.0);
-
 
314
										}
-
 
315
								}
-
 
316
								qem_vec[vi->touched] = q;
-
 
317
						}
-
 
318
						cout << "Pushing initial halfedges" << endl;
234
 
319
 
-
 
320
						for(HalfEdgeIter he = m.halfedges_begin();
-
 
321
								he != m.halfedges_end(); ++he)
235
			int work = 0;
322
								// For all halfedges, 
-
 
323
						{
236
			while(!sim_queue.empty() && work < max_work)
324
								if(halfedge_vec[he->touched].time_stamp == 0)
-
 
325
								{
-
 
326
										update_time_stamp(he);
-
 
327
										push_simplify_rec(he);
-
 
328
								}
237
				{
329
						}
-
 
330
 
238
					int dwork = collapse(sim_queue.top());
331
						cout << "Simplify" << endl;
-
 
332
 
239
					if(dwork>0)
333
						int work = 0;
-
 
334
						while(!sim_queue.empty() && work < max_work)
240
						{
335
						{
-
 
336
								SimplifyRec simplify_record = sim_queue.top();
241
							work += dwork;
337
								sim_queue.pop();
-
 
338
 
-
 
339
								work += collapse(simplify_record);
242
							if((work % 100) == 0)
340
								if((work % 100) == 0)
-
 
341
								{
243
								cout << "work = " << work << endl;
342
										cout << "work = " << work << endl;
-
 
343
										cout << "sim Q size = " << sim_queue.size() << endl;
-
 
344
								}
244
						}
345
						}
245
					sim_queue.pop();
-
 
246
				}
346
				}
247
		}
347
		}
248
	}
-
 
249
 
348
 
250
	void quadric_simplify(Manifold& m, int max_work)
349
		void quadric_simplify(Manifold& m, int max_work)
251
	{
350
		{
252
		srand(1210);
351
				srand(1210);
253
		QuadricSimplifier qsim(m);
352
				QuadricSimplifier qsim(m);
254
		qsim.reduce(max_work);
353
				qsim.reduce(max_work);
255
	}
354
		}
256
 
355
 
257
}
356
}
258
 
357