Subversion Repositories gelsvn

Rev

Rev 182 | Rev 219 | Go to most recent revision | Show entire file | Regard whitespace | Details | Blame | Last modification | View Log | RSS feed

Rev 182 Rev 215
Line 1... Line 1...
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
Line 13... Line 14...
13 
	using namespace GEO;
 14 
		using namespace GEO;
14 
	using namespace HMesh;
 15 
		using namespace HMesh;
15 
 
 16 
 
16 
	namespace
 17 
		namespace
17 
	{
 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). */
18 
		struct HalfEdgeRec
 23 
				struct HalfEdgeRec
19 
		{
 24 
				{
20 
			HalfEdgeIter he;
 25 
						HalfEdgeIter he;
21 
			int time_stamp;
 26 
						int time_stamp;
22 
			void halfedge_removed() {time_stamp = -1;}
 27 
						void halfedge_removed() {time_stamp = -1;}
23 
			HalfEdgeRec(): time_stamp(0) {}
 28 
						HalfEdgeRec(): time_stamp(0) {}
24 
		};
 29 
				};
25 
 
 30 
 
- 
 
 31 
				/* The simpliciation record contains information about a potential
- 
 
 32 
					 edge contraction */
26 
		struct SimplifyRec
 33 
				struct SimplifyRec
27 
		{
 34 
				{
28 
			Vec3f opt_pos;
 35 
						Vec3d opt_pos;  // optimal vertex position
29 
			int he_index;
 36 
						int he_index;   // Index (into HalfEdgeRec vector) of edge
- 
 
 37 
						                // we want to contract
30 
			float err;
 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
31 
			int time_stamp;
 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;
Line 45... Line 55...
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);
51 
 
 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 
								}
- 
 
 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 
								{
Line 65... Line 96...
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 
 
- 
 
 105 
				bool QuadricSimplifier::check_consistency(HalfEdgeIter he,
- 
 
 106 
																									const Vec3d& opt_pos)
- 
 
 107 
				{
- 
 
 108 
						VertexIter v0 = he->vert;
- 
 
 109 
						VertexIter v1 = he->opp->vert;
- 
 
 110 
						Vec3d p0(v0->pos);
- 
 
 111 
 
- 
 
 112 
						/* This test is inspired by Garland's Ph.D. thesis. We try
- 
 
 113 
							 to detect whether flipped triangles will occur by sort of
- 
 
 114 
							 ensuring that the new vertex is in the hull of the one rings
- 
 
 115 
							 of the vertices at either end of the edge being contracted */
- 
 
 116
- 
 
 117 
						for(VertexCirculator vc(v0); !vc.end(); ++vc)
- 
 
 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);
- 
 
 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)
- 
 
 131 
												return false;
- 
 
 132 
								}
- 
 
 133 
						}
- 
 
 134 
 
- 
 
 135 
						/* Since the test above is not really enough, we also make sure
- 
 
 136 
							 that we do not introduce vertices of very high (>=10) or low
- 
 
 137 
							 (<=3) valency. Flipped faces seem to be always associated with
- 
 
 138 
							 very irregular valencies. It seems to get rid of most problems
- 
 
 139 
							 not fixed by the check above. */
- 
 
 140
- 
 
 141 
						if(valency(he->next->vert)<5)
- 
 
 142 
								return false;
- 
 
 143 
 
- 
 
 144 
						if((valency(he->vert) + valency(he->opp->vert) ) > 12)
- 
 
 145 
								return false;
- 
 
 146 
 
- 
 
 147 
						return true;
- 
 
 148 
				}
- 
 
 149 
 
- 
 
 150 
 
73 
		void QuadricSimplifier::push_simplify_rec(HalfEdgeIter he)
 151 
				void QuadricSimplifier::push_simplify_rec(HalfEdgeIter he)
74 
		{
 152 
				{
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);
 - 
 
80 
			time_stamp++;
 - 
 
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
 153 
						// Get QEM for both end points
84 
			const QEM& Q1 = qem_vec[he->vert->touched];
 154 
						const QEM& Q1 = qem_vec[he->vert->touched];
85 
			const QEM& Q2 = qem_vec[he->opp->vert->touched];
 155 
						const QEM& Q2 = qem_vec[he->opp->vert->touched];
86 
 
 156 
 
87 
			QEM q = Q1;
 157 
						QEM q = Q1;
88 
			q += Q2;
 158 
						q += Q2;
89 
 
 159 
 
90 
			float err;
 160 
						float err;
91 
			Vec3f opt_pos(0);
 161 
						Vec3d opt_pos(0);
92 
 
 162 
 
93 
			//cout << q.determinant() << endl;
 163 
						//cout << q.determinant() << endl;
94 
			if (q.determinant()>= 1e-8)
 164 
						if (q.determinant()>= 1e-5)
95 
				{
 165 
						{
96 
					opt_pos = q.opt_pos(0.001);
 166 
								opt_pos = q.opt_pos(1e-5);
97 
					err = q.error(opt_pos);
 167 
								err = q.error(opt_pos);
98 
				}
 168 
						}
99 
			else
 169 
						else
100 
				{
 170 
						{
101 
					// To see the effectivenes of the scheme, we can create
171 
								// To see the effectivenes of the scheme, we can create
102 
					// Random errors for comparison. Uncomment to use:
 172 
								// Random errors for comparison. Uncomment to use:
103 
					// 			float e0 = time_stamp + rand()/float(RAND_MAX);//q.error(v0);
 173 
								// 			float e0 = time_stamp + rand()/float(RAND_MAX);//q.error(v0);
104 
					// 			float e1 = time_stamp + rand()/float(RAND_MAX);//q.error(v1);
 174 
								// 			float e1 = time_stamp + rand()/float(RAND_MAX);//q.error(v1);
105 
 
 175 
 
106 
					Vec3f v0 = he->vert->pos;
 176 
								Vec3d v0(he->vert->pos);
107 
					Vec3f v1 = he->opp->vert->pos;
 177 
								Vec3d v1(he->opp->vert->pos);
108 
					float e0 = q.error(v0);
 178 
								float e0 = q.error(v0);
109 
					float e1 = q.error(v1);
 179 
								float e1 = q.error(v1);
110
180
111 
					if(e0 < e1)
 181 
								if(e0 < e1)
112 
						{
 182 
								{
Line 124... Line 194...
124 
			int he_index = he->touched;
 194 
						int he_index = he->touched;
125 
			SimplifyRec simplify_rec;
 195 
						SimplifyRec simplify_rec;
126 
			simplify_rec.opt_pos = opt_pos;
 196 
						simplify_rec.opt_pos = opt_pos;
127 
			simplify_rec.err = err;
 197 
						simplify_rec.err = err;
128 
			simplify_rec.he_index = he_index;
 198 
						simplify_rec.he_index = he_index;
129 
			simplify_rec.time_stamp = time_stamp;
 199 
						simplify_rec.time_stamp = halfedge_vec[he->touched].time_stamp;
130 
 
 - 
 
- 
 
 200 
						simplify_rec.visits = 0;
131 
			// push it.
 201 
						// push it.
132 
			sim_queue.push(simplify_rec);
 202 
						sim_queue.push(simplify_rec);
- 
 
 203
133 
		}
 204 
				}
134
205
135 
 
 206 
 
136 
		void QuadricSimplifier::put_onering_in_queue(VertexIter vi)
 207 
				void QuadricSimplifier::update_onering_timestamp(VertexIter vi)
137 
		{
 208 
				{
138 
			// For all emanating edges he
 209 
						// For all emanating edges he
139 
			for(VertexCirculator vc(vi);
 210 
						for(VertexCirculator vc(vi);
140 
					!vc.end(); ++vc)
 211 
								!vc.end(); ++vc)
141 
				{
 212 
						{
142 
					HalfEdgeIter he = vc.get_halfedge();
 213 
								HalfEdgeIter he = vc.get_halfedge();
143 
					push_simplify_rec(he);
 214 
								update_time_stamp(he);
144 
				}
 215 
						}
145 
		}
 216 
				}
146 
 
 217 
 
147 
		int QuadricSimplifier::collapse(const SimplifyRec& simplify_rec)
 218 
				int QuadricSimplifier::collapse(SimplifyRec& simplify_rec)
148 
		{
 219 
				{
149 
			int he_index = simplify_rec.he_index;
 220 
						int he_index = simplify_rec.he_index;
- 
 
 221 
						HalfEdgeIter he = halfedge_vec[he_index].he;
150 
 
 222 
 
151 
			// Check the time stamp to verify that the simplification
223 
						// Check the time stamp to verify that the simplification
152 
			// record is the newest. If the halfedge has been removed
 224 
						// record is the newest. If the halfedge has been removed
153 
			// the time stamp is -1 and the comparison will also fail.
 225 
						// the time stamp is -1 and the comparison will also fail.
154 
			if(halfedge_vec[he_index].time_stamp == simplify_rec.time_stamp)
 226 
						if(halfedge_vec[he_index].time_stamp == simplify_rec.time_stamp)
155 
				{
 227 
						{
156 
					HalfEdgeIter he = halfedge_vec[he_index].he;
 - 
 
157 
					VertexIter v = he->opp->vert;
 228 
								VertexIter v = he->opp->vert;
158 
					VertexIter n = he->vert;
 229 
								VertexIter n = he->vert;
159 
 
 230 
 
160 
					// If the edge is, in fact, collapsible
 231 
								// If the edge is, in fact, collapsible
161 
					if(m.collapse_precond(he))
 232 
								if(m.collapse_precond(he))
162 
						{
 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 
 
163 
							//cout << simplify_rec.err << " " << &(*he->vert) << endl;
 240 
												//cout << simplify_rec.err << " " << &(*he->vert) << endl;
164 
							// Get QEM for both end points
 241 
												// Get QEM for both end points
165 
							const QEM& Q1 = qem_vec[he->vert->touched];
 242 
												const QEM& Q1 = qem_vec[he->vert->touched];
166 
							const QEM& Q2 = qem_vec[he->opp->vert->touched];
 243 
												const QEM& Q2 = qem_vec[he->opp->vert->touched];
167 
 
 244
Line 184... Line 261...
184 
									halfedge_vec[he->opp->next->next->touched].halfedge_removed();
 261 
														halfedge_vec[he->opp->next->next->touched].halfedge_removed();
185 
								}
 262 
												}
186 
 
 263
187 
							// Do collapse
 264 
												// Do collapse
188 
							m.collapse_halfedge(he,false);
 265 
												m.collapse_halfedge(he,false);
189 
							n->pos = simplify_rec.opt_pos;
 266 
												n->pos = Vec3f(simplify_rec.opt_pos);
190 
							qem_vec[n->touched] = q;
 267 
												qem_vec[n->touched] = q;
191
268
192 
							put_onering_in_queue(n);
 269 
												update_onering_timestamp(n);
193 
							return 1;
 270 
												return 1;
194 
						}
 271 
										}
195 
				}
 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;
- 
 
 280 
										++simplify_rec.visits;
- 
 
 281 
										sim_queue.push(simplify_rec);
- 
 
 282 
								}
- 
 
 283 
						}
- 
 
 284 
						// Ok the time stamp did not match. Create a new record.
- 
 
 285 
						else if(halfedge_vec[he_index].time_stamp != -1)
- 
 
 286 
								push_simplify_rec(he);
- 
 
 287 
 
196 
			return 0;
 288 
						return 0;
197 
		}
 289 
				}
198 
 
 290 
 
199 
		void QuadricSimplifier::reduce(int max_work)
 291 
				void QuadricSimplifier::reduce(int max_work)
200 
		{
 292 
				{
201 
 
 - 
 
202 
			// Set t = 0 for all halfedges
 293 
						// Set t = 0 for all halfedges
203 
			int i=0;
 294 
						int i=0;
204 
			for(HalfEdgeIter he = m.halfedges_begin();
 295 
						for(HalfEdgeIter he = m.halfedges_begin();
205 
					he != m.halfedges_end(); ++he, ++i)
 296 
								he != m.halfedges_end(); ++he, ++i)
206 
				halfedge_vec[i].he = he;
 297 
								halfedge_vec[i].he = he;
207
298 
 
- 
 
 299 
						cout << "Computing quadrics" << endl;
- 
 
 300
208 
			// For all vertices, compute quadric and store in qem_vec
 301 
						// For all vertices, compute quadric and store in qem_vec
209 
			for(VertexIter vi=m.vertices_begin(); vi != m.vertices_end(); ++vi)
 302 
						for(VertexIter vi=m.vertices_begin(); vi != m.vertices_end(); ++vi)
210 
				{
 303 
						{
211 
					Vec3d p(vi->pos);
 304 
								Vec3d p(vi->pos);
212 
					QEM q;
 305 
								QEM q;
Line 220... Line 313...
220 
									q += QEM(p, n, a / 3.0);
 313 
												q += QEM(p, n, a / 3.0);
221 
								}
 314 
										}
222 
						}
 315 
								}
223 
					qem_vec[vi->touched] = q;
 316 
								qem_vec[vi->touched] = q;
224 
				}
 317 
						}
- 
 
 318 
						cout << "Pushing initial halfedges" << endl;
225 
 
 319 
 
226 
			for(HalfEdgeIter he = m.halfedges_begin();
 320 
						for(HalfEdgeIter he = m.halfedges_begin();
227 
					he != m.halfedges_end(); ++he)
 321 
								he != m.halfedges_end(); ++he)
228 
				// For all halfedges,
322 
								// For all halfedges,
229 
				{
 323 
						{
230 
					if(halfedge_vec[he->touched].time_stamp == 0)
 324 
								if(halfedge_vec[he->touched].time_stamp == 0)
- 
 
 325 
								{
- 
 
 326 
										update_time_stamp(he);
231 
						push_simplify_rec(he);
 327 
										push_simplify_rec(he);
232 
				}
 328 
								}
- 
 
 329 
						}
233 
 
 330 
 
- 
 
 331 
						cout << "Simplify" << endl;
234 
 
 332 
 
235 
			int work = 0;
 333 
						int work = 0;
236 
			while(!sim_queue.empty() && work < max_work)
 334 
						while(!sim_queue.empty() && work < max_work)
237 
				{
 335 
						{
238 
					int dwork = collapse(sim_queue.top());
 336 
								SimplifyRec simplify_record = sim_queue.top();
239 
					if(dwork>0)
 337 
								sim_queue.pop();
240 
						{
 338 
 
241 
							work += dwork;
 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;
244 
						}
 344 
								}
245 
					sim_queue.pop();
 - 
 
246 
				}
 345 
						}
247 
		}
 346 
				}
248 
	}
 347 
		}
249 
 
 348 
 
250 
	void quadric_simplify(Manifold& m, int max_work)
 349 
		void quadric_simplify(Manifold& m, int max_work)