WebSVN – gelsvn – Diff – /trunk/src/HMesh/quadric_simplify.cpp


#include <queue>
#include <iostream>
#include "quadric_simplify.h"
#include "smooth.h"
#include "HMesh/VertexCirculator.h"
#include "Geometry/QEM.h"
 
 
namespace HMesh
{
	using namespace CGLA;
	using namespace std;
	using namespace GEO;
	using namespace HMesh;
	
	namespace
	{
		/* We create a record for each halfedge where we can keep its time
		stamp. If the time stamp on the halfedge record is bigger than
		the stamp on the simplification record, we cannot use the 
		simplification record (see below). */
		struct HalfEdgeRec
				{
					HalfEdgeIter he;
					int time_stamp;
					void halfedge_removed() {time_stamp = -1;}
					HalfEdgeRec(): time_stamp(0) {}
				};
		
		/* The simpliciation record contains information about a potential
		edge contraction */
		struct SimplifyRec
		{
			Vec3d opt_pos;  // optimal vertex position
			int he_index;   // Index (into HalfEdgeRec vector) of edge
							// we want to contract
			float err;      // Error associated with contraction
			int time_stamp; // Time stamp (see comment on HalfEdgeRec)
			int visits;     // Visits (number of times we considered this 
							// record).
		};
		
		bool operator<(const SimplifyRec& s1, const SimplifyRec& s2)
		{
			return s1.err > s2.err;
		}
		
		class QuadricSimplifier
		{
			typedef priority_queue<SimplifyRec> SimplifyQueue;
			typedef vector<QEM> QEMVec;
			typedef vector<HalfEdgeRec> HalfEdgeVec;
			
			Manifold& m;
			HalfEdgeVec halfedge_vec;
			QEMVec qem_vec;
			SimplifyQueue sim_queue;
			double singular_thresh;
			bool choose_optimal_positions;
			
			/* Compute the error associated with contraction of he and the
				optimal position of resulting vertex. */
			void push_simplify_rec(HalfEdgeIter he);
			
			/* Check whether the contraction is valid. See below for details*/
			bool check_consistency(HalfEdgeIter he, const Vec3d& opt_pos);
			
			/* Update the time stamp of a halfedge. A halfedge and its opp edge
				may have different stamps. We choose a stamp that is greater
				than either and assign to both.*/
			void update_time_stamp(HalfEdgeIter he)
			{
				int time_stamp = s_max(halfedge_vec[he->touched].time_stamp,
									   halfedge_vec[he->opp->touched].time_stamp);
				time_stamp++;
				halfedge_vec[he->touched].time_stamp = time_stamp;
				halfedge_vec[he->opp->touched].time_stamp = time_stamp;
			}
			
			/* Update time stamps for all halfedges in one ring of vi */
			void update_onering_timestamp(VertexIter vi);
			
			/* Perform a collapse - if conditions are met. Returns 1 or 0 
				accordingly. */
			int collapse(SimplifyRec& simplify_rec);
			
		public:
				
				/* Create a simplifier for a manifold */
				QuadricSimplifier(Manifold& _m, double _singular_thresh, bool _choose_optimal_positions): 
				m(_m), 
				halfedge_vec(m.no_halfedges()), 
				qem_vec(m.no_vertices()),
				singular_thresh(_singular_thresh),
				choose_optimal_positions(_choose_optimal_positions)
			{
					// Enumerate vertices
					m.enumerate_vertices();
					
					// Enumerate halfedges
					m.enumerate_halfedges();
			}
			
			/* Simplify doing at most max_work contractions */
			void reduce(int max_work);
		};
		
		bool QuadricSimplifier::check_consistency(HalfEdgeIter he, 
												  const Vec3d& opt_pos)
		{
			VertexIter v0 = he->vert;
			VertexIter v1 = he->opp->vert;
			Vec3d p0(v0->pos);
			
			/* This test is inspired by Garland's Ph.D. thesis. We try
				to detect whether flipped triangles will occur by sort of
				ensuring that the new vertex is in the hull of the one rings
				of the vertices at either end of the edge being contracted 
				
				I also had an additional check intended to ensure that poor valencies
				would not be introduced, but it seemed to be unnecessary.
				
				*/
			
			for(VertexCirculator vc(v0); !vc.end(); ++vc)
			{
				HalfEdgeIter h = vc.get_halfedge();
				
				if(h->vert != v1 && h->next->vert != v1)
				{
					Vec3d pa(h->vert->pos);
					Vec3d pb(h->next->vert->pos);
					
					Vec3d dir = normalize(pb-pa);
					
					Vec3d n = p0 - pa;
					n = n - dir * dot(dir,n);
					
					if(dot(n,opt_pos - pa) <= 0)
						return false;
				}
			}
			
			return true;
		}
		
		
		void QuadricSimplifier::push_simplify_rec(HalfEdgeIter he)
		{
			// Get QEM for both end points
			const QEM& Q1 = qem_vec[he->vert->touched];
			const QEM& Q2 = qem_vec[he->opp->vert->touched];
			
			QEM q = Q1;
			q += Q2;
			
			float err;
			Vec3d opt_pos(0);
			Vec3d opt_origin = Vec3d(he->vert->pos+he->opp->vert->pos)/2.0;
			if(choose_optimal_positions)
				opt_pos = q.opt_pos(singular_thresh,opt_origin);
			else
			{
				Vec3d p1(he->vert->pos), p2(he->opp->vert->pos);
				float err1 = q.error(p1);
				float err2 = q.error(p2);
				if(err1<err2) 
				{
					err = err1;
					opt_pos = p1;
				}
				else
				{
					err = err2;
					opt_pos = p2;				
				}
			}
			
			err = q.error(opt_pos);
			
			// Create SimplifyRec
			int he_index = he->touched;
			SimplifyRec simplify_rec;
			simplify_rec.opt_pos = opt_pos;
			simplify_rec.err = err;
			simplify_rec.he_index = he_index;
			simplify_rec.time_stamp = halfedge_vec[he->touched].time_stamp;
			simplify_rec.visits = 0;
			// push it.
			sim_queue.push(simplify_rec);
			
		}
		
		
		void QuadricSimplifier::update_onering_timestamp(VertexIter vi)
		{
			// For all emanating edges he
			for(VertexCirculator vc(vi);
				!vc.end(); ++vc)
			{
				HalfEdgeIter he = vc.get_halfedge();
				update_time_stamp(he);
				push_simplify_rec(he);
			}
		}
		
		int QuadricSimplifier::collapse(SimplifyRec& simplify_rec)
		{
			int he_index = simplify_rec.he_index;
			HalfEdgeIter he = halfedge_vec[he_index].he;
			
			// Check the time stamp to verify that the simplification 
			// record is the newest. If the halfedge has been removed
			// the time stamp is -1 and the comparison will also fail.
			if(halfedge_vec[he_index].time_stamp == simplify_rec.time_stamp)
			{
				VertexIter v = he->opp->vert;
				VertexIter n = he->vert;
				
				// If the edge is, in fact, collapsible
				if(m.collapse_precond(he))
				{
					// If our consistency checks pass, we are relatively
					// sure that the contraction does not lead to a face flip.
					if(check_consistency(he, simplify_rec.opt_pos) && 
					   check_consistency(he->opp, simplify_rec.opt_pos))
					{
						
						//cout << simplify_rec.err << " " << &(*he->vert) << endl;
						// Get QEM for both end points
						const QEM& Q1 = qem_vec[he->vert->touched];
						const QEM& Q2 = qem_vec[he->opp->vert->touched];
						
						// Compute Q_new = Q_1 + Q_2
						QEM q = Q1;
						q += Q2;
						
						// Mark all halfedges that will be removed as dead
						halfedge_vec[he->touched].halfedge_removed();
						halfedge_vec[he->opp->touched].halfedge_removed();
						
						if(he->next->next->next == he)
						{
							halfedge_vec[he->next->touched].halfedge_removed();
							halfedge_vec[he->next->next->touched].halfedge_removed();
						}
						if(he->opp->next->next->next == he->opp)
						{
							halfedge_vec[he->opp->next->touched].halfedge_removed();
							halfedge_vec[he->opp->next->next->touched].halfedge_removed();
						}
						
						// Do collapse
						m.collapse_halfedge(he,false);
						n->pos = Vec3f(simplify_rec.opt_pos);
						qem_vec[n->touched] = q;
						
						update_onering_timestamp(n);
						return 1;
					}
				}
				// If we are here, the collapse was not allowed. If we have
				// seen this simplify record less than 100 times, we try to
				// increase the error and store the record again. Maybe some
				// other contractions will make it more digestible later.
				if(simplify_rec.visits < 100)
				{
					simplify_rec.err *= 1.01f;
					++simplify_rec.visits;
					sim_queue.push(simplify_rec);
				}
			}
			return 0;
		}
		
		void QuadricSimplifier::reduce(int max_work)
		{
			// Set t = 0 for all halfedges
			int i=0;
			for(HalfEdgeIter he = m.halfedges_begin();
				he != m.halfedges_end(); ++he, ++i)
				halfedge_vec[i].he = he;
			
			cout << "Computing quadrics" << endl;
			
			// For all vertices, compute quadric and store in qem_vec
			for(VertexIter vi=m.vertices_begin(); vi != m.vertices_end(); ++vi)
			{
				Vec3d p(vi->pos);
				Vec3d vn(normal(vi));
				QEM q;
				for(VertexCirculator vc(vi); !vc.end(); ++vc)
				{
					FaceIter f = vc.get_face();
					if(f != NULL_FACE_ITER)
					{
						Vec3d n(normal(f));
						double a = area(f);
						q += QEM(p, n, a / 3.0);
					}
					if ((f == NULL_FACE_ITER || 
						vc.get_opp_halfedge()->face ==  NULL_FACE_ITER) && sqr_length(vn) > 0.0)
					{
						Vec3d edge = Vec3d(vc.get_halfedge()->vert->pos)-p;
						double edge_len = sqr_length(edge); 
						if(edge_len > 0.0)
							{
								Vec3d n = cross(vn, edge);
								q += QEM(p, n, 2*edge_len);
							}
					}
				}
				qem_vec[vi->touched] = q;
			}
			cout << "Pushing initial halfedges" << endl;
			
			for(HalfEdgeIter he = m.halfedges_begin();
				he != m.halfedges_end(); ++he)
				// For all halfedges, 
			{
				if(halfedge_vec[he->touched].time_stamp == 0)
				{
					update_time_stamp(he);
					push_simplify_rec(he);
				}
			}
			
			cout << "Simplify" << endl;
			
			int work = 0;
			while(!sim_queue.empty() && work < max_work)
			{
				SimplifyRec simplify_record = sim_queue.top();
				sim_queue.pop();
				
				work += 2*collapse(simplify_record);
				if((work % 100) == 0)
				{
					cout << "work = " << work << endl;
					cout << "sim Q size = " << sim_queue.size() << endl;
				}
			}
		}
	}
	
	void quadric_simplify(Manifold& m, double keep_fraction, double singular_thresh, bool choose_optimal_positions)
	{
		gel_srand(1210);
		int F = m.no_faces();
		int max_work = max(0, F- static_cast<int>(keep_fraction * F));
		QuadricSimplifier qsim(m, singular_thresh, choose_optimal_positions);
		qsim.reduce(max_work);
	}
	
}
 

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

Subversion Repositories gelsvn

(root)/trunk/src/HMesh/quadric_simplify.cpp – Rev 373 → 417