Subversion Repositories gelsvn

Rev

Rev 315 | Rev 362 | Go to most recent revision | Show entire file | Ignore whitespace | Details | Blame | Last modification | View Log | RSS feed

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