Subversion Repositories gelsvn

Rev

Rev 149 | Rev 155 | Go to most recent revision | Details | Compare with Previous | Last modification | View Log | RSS feed

Rev Author Line No. Line
149 jab 1 
#include <queue>
2 
#include <iostream>
3 
#include "quadric_simplify.h"
4 
#include "smooth.h"
5 
#include "HMesh/VertexCirculator.h"
6 
#include "Geometry/QEM.h"
7 
 
8 
 
150 jab 9 
namespace HMesh
149 jab 10 
{
11 
  using namespace CGLA;
12 
  using namespace std;
13 
	using namespace GEO;
14 
	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 
 
26 
		struct SimplifyRec
27 
		{
28 
			Vec3f opt_pos;
29 
			int he_index;
30 
			float err;
31 
			int time_stamp;
32 
		};
33 
 
34 
		bool operator<(const SimplifyRec& s1, const SimplifyRec& s2)
35 
		{
36 
			return s1.err > s2.err;
37 
		}
38 
 
39 
		class QuadricSimplifier
40 
		{
41 
			typedef priority_queue<SimplifyRec> SimplifyQueue;
42 
			typedef vector<QEM> QEMVec;
43 
			typedef vector<HalfEdgeRec> HalfEdgeVec;
44 
 
45 
			Manifold& m;
46 
			HalfEdgeVec halfedge_vec;
47 
			QEMVec qem_vec;
48 
			SimplifyQueue sim_queue;
49 
 
50 
			void push_simplify_rec(HalfEdgeIter he);
51 
 
52 
			void put_onering_in_queue(VertexIter vi);
53 
 
54 
			int collapse(const SimplifyRec& simplify_rec);
55 
 
56 
		public:
57 
 
58 
			QuadricSimplifier(Manifold& _m):
59 
				m(_m),
60 
				halfedge_vec(m.no_halfedges()),
61 
				qem_vec(m.no_vertices())
62 
			{
63 
				// Enumerate vertices
64 
				m.enumerate_vertices();
65 
 
66 
				// Enumerate halfedges
67 
				m.enumerate_halfedges();
68 
			}
69 
 
70 
			void reduce(int max_work);
71 
		};
72 
 
73 
		void QuadricSimplifier::push_simplify_rec(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);
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
84 
			const QEM& Q1 = qem_vec[he->vert->touched];
85 
			const QEM& Q2 = qem_vec[he->opp->vert->touched];
86 
 
87 
			QEM q = Q1;
88 
			q += Q2;
89 
 
90 
			float err;
91 
			Vec3f opt_pos(0);
92 
 
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
102 
					// Random errors for comparison. Uncomment to use:
103 
					// 			float e0 = time_stamp + rand()/float(RAND_MAX);//q.error(v0);
104 
					// 			float e1 = time_stamp + rand()/float(RAND_MAX);//q.error(v1);
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 
						{
113 
							err = e0;
114 
							opt_pos = v0;
115 
						}
116 
					else
117 
						{
118 
							err = e1;
119 
							opt_pos = v1;
120 
						}
121 
				}
122 
 
123 
			// Create SimplifyRec
124 
			int he_index = he->touched;
125 
			SimplifyRec simplify_rec;
126 
			simplify_rec.opt_pos = opt_pos;
127 
			simplify_rec.err = err;
128 
			simplify_rec.he_index = he_index;
129 
			simplify_rec.time_stamp = time_stamp;
130 
 
131 
			// push it.
132 
			sim_queue.push(simplify_rec);
133 
		}
134 
 
135 
 
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)
141 
				{
142 
					HalfEdgeIter he = vc.get_halfedge();
143 
					push_simplify_rec(he);
144 
				}
145 
		}
146 
 
147 
		int QuadricSimplifier::collapse(const SimplifyRec& simplify_rec)
148 
		{
149 
			int he_index = simplify_rec.he_index;
150 
 
151 
			// Check the time stamp to verify that the simplification
152 
			// record is the newest. If the halfedge has been removed
153 
			// the time stamp is -1 and the comparison will also fail.
154 
			if(halfedge_vec[he_index].time_stamp == simplify_rec.time_stamp)
155 
				{
156 
					HalfEdgeIter he = halfedge_vec[he_index].he;
157 
					VertexIter v = he->opp->vert;
158 
					VertexIter n = he->vert;
159 
 
160 
					// If the edge is, in fact, collapsible
161 
					if(m.collapse_precond(v,he))
162 
						{
163 
							//cout << simplify_rec.err << " " << &(*he->vert) << endl;
164 
							// Get QEM for both end points
165 
							const QEM& Q1 = qem_vec[he->vert->touched];
166 
							const QEM& Q2 = qem_vec[he->opp->vert->touched];
167 
 
168 
							// Compute Q_new = Q_1 + Q_2
169 
							QEM q = Q1;
170 
							q += Q2;
171 
 
172 
							// Mark all halfedges that will be removed as dead
173 
							halfedge_vec[he->touched].halfedge_removed();
174 
							halfedge_vec[he->opp->touched].halfedge_removed();
175 
 
176 
							if(he->next->next->next == he)
177 
								{
178 
									halfedge_vec[he->next->touched].halfedge_removed();
179 
									halfedge_vec[he->next->next->touched].halfedge_removed();
180 
								}
181 
							if(he->opp->next->next->next == he->opp)
182 
								{
183 
									halfedge_vec[he->opp->next->touched].halfedge_removed();
184 
									halfedge_vec[he->opp->next->next->touched].halfedge_removed();
185 
								}
186 
 
187 
							// Do collapse
188 
							m.collapse_halfedge(v,he,false);
189 
							n->pos = simplify_rec.opt_pos;
190 
							qem_vec[n->touched] = q;
191 
 
192 
							put_onering_in_queue(n);
193 
							return 1;
194 
						}
195 
				}
196 
			return 0;
197 
		}
198 
 
199 
		void QuadricSimplifier::reduce(int max_work)
200 
		{
201 
 
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
209 
			for(VertexIter vi=m.vertices_begin(); vi != m.vertices_end(); ++vi)
210 
				{
211 
					Vec3d p(vi->pos);
212 
					QEM q;
213 
					for(VertexCirculator vc(vi); !vc.end(); ++vc)
214 
						{
215 
							FaceIter f = vc.get_face();
216 
							if(f != NULL_FACE_ITER)
217 
								{
218 
									Vec3d n(normal(f));
219 
									double a = area(f);
220 
									q += QEM(p, n, a / 3.0);
221 
								}
222 
						}
223 
					qem_vec[vi->touched] = q;
224 
				}
225 
 
226 
			for(HalfEdgeIter he = m.halfedges_begin();
227 
					he != m.halfedges_end(); ++he)
228 
				// For all halfedges,
229 
				{
230 
					if(halfedge_vec[he->touched].time_stamp == 0)
231 
						push_simplify_rec(he);
232 
				}
233 
 
234 
 
235 
			int work = 0;
236 
			while(!sim_queue.empty() && work < max_work)
237 
				{
238 
					int dwork = collapse(sim_queue.top());
239 
					if(dwork>0)
240 
						{
241 
							work += dwork;
242 
							if((work % 100) == 0)
243 
								cout << "work = " << work << endl;
244 
						}
245 
					sim_queue.pop();
246 
				}
247 
		}
248 
	}
249 
 
250 
	void quadric_simplify(Manifold& m, int max_work)
251 
	{
252 
		srand(1210);
253 
		QuadricSimplifier qsim(m);
254 
		qsim.reduce(max_work);
255 
	}
256 
 
257 
}