Subversion Repositories gelsvn

Rev

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

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