Subversion Repositories gelsvn

#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;

						/* 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): 

								m(_m), 

								halfedge_vec(m.no_halfedges()), 

								qem_vec(m.no_vertices()) 

								{

										// 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 */

						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;

								}

						}

						/* Since the test above is not really enough, we also make sure 

							 that we do not introduce vertices of very high (>=10) or low 

							 (<=3) valency. Flipped faces seem to be always associated with

							 very irregular valencies. It seems to get rid of most problems

							 not fixed by the check above. */

						if(valency(he->next->vert)<5)

								return false;

						if((valency(he->vert) + valency(he->opp->vert) ) > 13)

								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);

						if (q.determinant()>= 1e-15)

						{

								opt_pos = q.opt_pos(0.001);

								err = q.error(opt_pos);

						}

						else

						{

								// To see the effectivenes of the scheme, we can create 

								// Random errors for comparison. Uncomment to use:

								// 			float e0 = time_stamp + rand()/float(RAND_MAX);//q.error(v0);

								// 			float e1 = time_stamp + rand()/float(RAND_MAX);//q.error(v1);

								Vec3d v0(he->vert->pos);

								Vec3d v1(he->opp->vert->pos);

								float e0 = q.error(v0);

								float e1 = q.error(v1);

								if(e0 < e1)

								{

										err = e0;

										opt_pos = v0;

								}

								else

								{

										err = e1;

										opt_pos = v1;

								}

						}			

						// 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);

						}

				}

				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 < 10)

								{

										simplify_rec.err += 5*simplify_rec.err + 1e-10;

										++simplify_rec.visits;

										sim_queue.push(simplify_rec);

								}

						}

						// Ok the time stamp did not match. Create a new record.

						else if(halfedge_vec[he_index].time_stamp != -1)

								push_simplify_rec(he);

						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);

								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);

										}

								}

								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 += collapse(simplify_record);

								if((work % 100) == 0)

								{

										cout << "work = " << work << endl;

										cout << "sim Q size = " << sim_queue.size() << endl;

								}

						}

				}

		}

		void quadric_simplify(Manifold& m, int max_work)

		{

				srand(1210);

				QuadricSimplifier qsim(m);

				qsim.reduce(max_work);

		}

}