Subversion Repositories gelsvn

Rev

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

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