Subversion Repositories gelsvn

#include "quadric_simplify.h"

#include "AttributeVector.h"

	using namespace GEO;

	using namespace CGLA;

	using namespace HMesh;

    using namespace Geometry;

	namespace

    namespace

	{

    {

		/* We create a record for each halfedge where we can keep its time

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

         stamp. If the time stamp on the halfedge record is bigger than

		the stamp on the simplification record, we cannot use the 

         the stamp on the simplification record, we cannot use the 

		simplification record (see below). */

         simplification record (see below). */

		struct HalfEdgeRec

        struct HalfEdgeRec

				{

        {

					HalfEdgeIter he;

            HalfEdgeID h;

					int time_stamp;

            int time_stamp;

					void halfedge_removed() {time_stamp = -1;}

            void halfedge_removed() {time_stamp = -1;}

					HalfEdgeRec(): time_stamp(0) {}

            HalfEdgeRec(): time_stamp(0) {}

				};

        };

		/* The simpliciation record contains information about a potential

        /* The simpliciation record contains information about a potential

		edge contraction */

         edge contraction */

		struct SimplifyRec

        struct SimplifyRec

		{

        {

			Vec3d opt_pos;  // optimal vertex position

            Vec3d opt_pos;  // optimal vertex position

			int he_index;   // Index (into HalfEdgeRec vector) of edge

            HalfEdgeID h;   // Index (into HalfEdgeRec vector) of edge

							// we want to contract

            // we want to contract

			float err;      // Error associated with contraction

            float err;      // Error associated with contraction

			int time_stamp; // Time stamp (see comment on HalfEdgeRec)

            int time_stamp; // Time stamp (see comment on HalfEdgeRec)

			int visits;     // Visits (number of times we considered this 

            int visits;     // Visits (number of times we considered this 

							// record).

            // record).

		};

        };

		bool operator<(const SimplifyRec& s1, const SimplifyRec& s2)

        bool operator<(const SimplifyRec& s1, const SimplifyRec& s2)

		{

        {

			return s1.err > s2.err;

            return s1.err > s2.err;

		}

        }

		class QuadricSimplifier

        class QuadricSimplifier

		{

        {

			typedef priority_queue<SimplifyRec> SimplifyQueue;

            typedef priority_queue<SimplifyRec> SimplifyQueue;

			typedef vector<QEM> QEMVec;

            typedef VertexAttributeVector<QEM> QEMVec;

			typedef vector<HalfEdgeRec> HalfEdgeVec;

            typedef HalfEdgeAttributeVector<HalfEdgeRec> HalfEdgeVec;

            typedef VertexAttributeVector<int> CollapseMask;

			Manifold& m;

			HalfEdgeVec halfedge_vec;

            Manifold& m;

			QEMVec qem_vec;

            HalfEdgeVec halfedge_vec;

			SimplifyQueue sim_queue;

            QEMVec qem_vec;

			double singular_thresh;

            CollapseMask collapse_mask;

			bool choose_optimal_positions;

            SimplifyQueue sim_queue;

            double singular_thresh;

			/* Compute the error associated with contraction of he and the

            bool choose_optimal_positions;

				optimal position of resulting vertex. */

			void push_simplify_rec(HalfEdgeIter he);

            /* Compute the error associated with contraction of he and the

             optimal position of resulting vertex. */

			/* Check whether the contraction is valid. See below for details*/

            void push_simplify_rec(HalfEdgeID h);

			bool check_consistency(HalfEdgeIter he, const Vec3d& opt_pos);

            /* Check whether the contraction is valid. See below for details*/

			/* Update the time stamp of a halfedge. A halfedge and its opp edge

            bool check_consistency(HalfEdgeID h, const Vec3d& opt_pos);

				may have different stamps. We choose a stamp that is greater

				than either and assign to both.*/

            /* Update the time stamp of a halfedge. A halfedge and its opp edge

			void update_time_stamp(HalfEdgeIter he)

             may have different stamps. We choose a stamp that is greater

			{

             than either and assign to both.*/

				int time_stamp = s_max(halfedge_vec[he->touched].time_stamp,

            void update_time_stamp(HalfEdgeID h)

									   halfedge_vec[he->opp->touched].time_stamp);

            {

				time_stamp++;

				Walker w = m.walker(h);

				halfedge_vec[he->touched].time_stamp = time_stamp;

				HalfEdgeID ho = w.opp().halfedge();

				halfedge_vec[he->opp->touched].time_stamp = time_stamp;

			}

                int time_stamp = s_max( halfedge_vec[h].time_stamp, halfedge_vec[ho].time_stamp);

                ++time_stamp;

			/* Update time stamps for all halfedges in one ring of vi */

                halfedge_vec[h].time_stamp = time_stamp;

			void update_onering_timestamp(VertexIter vi);

                halfedge_vec[ho].time_stamp = time_stamp;

            }

			/* Perform a collapse - if conditions are met. Returns 1 or 0 

				accordingly. */

            /* Update time stamps for all halfedges in one ring of vi */

			int collapse(SimplifyRec& simplify_rec);

            void update_onering_timestamp(VertexID v);

		public:

            /* Perform a collapse - if conditions are met. Returns 1 or 0 

             accordingly. */

				/* Create a simplifier for a manifold */

            int collapse(SimplifyRec& simplify_rec);

				QuadricSimplifier(Manifold& _m, double _singular_thresh, bool _choose_optimal_positions): 

				m(_m), 

        public:

				halfedge_vec(m.no_halfedges()), 

				qem_vec(m.no_vertices()),

            /* Create a simplifier for a manifold */

				singular_thresh(_singular_thresh),

            QuadricSimplifier(Manifold& _m, VertexAttributeVector<int> _collapse_mask,

				choose_optimal_positions(_choose_optimal_positions)

                              double _singular_thresh, bool _choose_optimal_positions):

			{

            m(_m), 

					// Enumerate vertices

            halfedge_vec(_m.allocated_halfedges()), 

					m.enumerate_vertices();

            qem_vec(_m.allocated_vertices()),

            collapse_mask(_collapse_mask),

					// Enumerate halfedges

            singular_thresh(_singular_thresh),

					m.enumerate_halfedges();

            choose_optimal_positions(_choose_optimal_positions)

			}

            {}

			/* Simplify doing at most max_work contractions */

            /* Simplify doing at most max_work contractions */

			void reduce(int max_work);

            void reduce(long int max_work);

		};

        };

		bool QuadricSimplifier::check_consistency(HalfEdgeIter he, 

        bool QuadricSimplifier::check_consistency(HalfEdgeID h, const Vec3d& opt_pos)

												  const Vec3d& opt_pos)

        {

		{

			Walker w = m.walker(h);

			VertexIter v0 = he->vert;

			VertexIter v1 = he->opp->vert;

			Vec3d p0(v0->pos);

            VertexID v0 = w.vertex();

            VertexID v1 = w.opp().vertex();

			/* This test is inspired by Garland's Ph.D. thesis. We try

            Vec3d p0(m.pos(v0));

				to detect whether flipped triangles will occur by sort of

				ensuring that the new vertex is in the hull of the one rings

            /* This test is inspired by Garland's Ph.D. thesis. We try

				of the vertices at either end of the edge being contracted 

             to detect whether flipped triangles will occur by sort of

             ensuring that the new vertex is in the hull of the one rings

				I also had an additional check intended to ensure that poor valencies

             of the vertices at either end of the edge being contracted 

				would not be introduced, but it seemed to be unnecessary.

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

			for(Walker w = m.walker(v0); !w.full_circle(); w = w.circulate_vertex_cw()){

                //ConstHalfEdgeHandle h = vc.halfedge();

				if(h->vert != v1 && h->next->vert != v1)

				{

                if(w.vertex()!= v1 && w.next().vertex() != v1){

					Vec3d pa(h->vert->pos);

                    Vec3d pa(m.pos(w.vertex()));

					Vec3d pb(h->next->vert->pos);

                    Vec3d pb(m.pos(w.next().vertex()));

					Vec3d dir = normalize(pb-pa);

                    Vec3d dir = normalize(pb - pa);

					Vec3d n = p0 - pa;

                    Vec3d n = p0 - pa;

					n = n - dir * dot(dir,n);

                    n = n - dir * dot(dir,n);

					if(dot(n,opt_pos - pa) <= 0)

                    if(dot(n,opt_pos - pa) <= 0)

						return false;

                        return false;

				}

                }

			}

            }

			return true;

            return true;

		}

        }

		void QuadricSimplifier::push_simplify_rec(HalfEdgeIter he)

        void QuadricSimplifier::push_simplify_rec(HalfEdgeID h)

		{

        {

			// Get QEM for both end points

			Walker w = m.walker(h);

			const QEM& Q1 = qem_vec[he->vert->touched];

            if(collapse_mask[w.opp().vertex()] == 0)

			const QEM& Q2 = qem_vec[he->opp->vert->touched];

            {

			QEM q = Q1;

                update_time_stamp(h);

			q += Q2;

                VertexID hv = w.vertex();

			float err;

                VertexID hov = w.opp().vertex();

			Vec3d opt_pos(0);

			Vec3d opt_origin = Vec3d(he->vert->pos+he->opp->vert->pos)/2.0;

                // Get QEM for both end points

			if(choose_optimal_positions)

                const QEM& Q1 = qem_vec[hv];

				opt_pos = q.opt_pos(singular_thresh,opt_origin);

                const QEM& Q2 = qem_vec[hov];

			else

			{

                QEM q = Q1;

				Vec3d p1(he->vert->pos), p2(he->opp->vert->pos);

                q += Q2;

				float err1 = q.error(p1);

				float err2 = q.error(p2);

                float err;

				if(err1<err2) 

                Vec3d opt_pos(0);

				{

                Vec3d opt_origin = Vec3d(m.pos(hv) + m.pos(hov)) * 0.5;

					err = err1;

                if(choose_optimal_positions)

					opt_pos = p1;

                    opt_pos = q.opt_pos(singular_thresh,opt_origin);

				}

                else

				else

                    opt_pos = Vec3d(m.pos(hv));

				{

					err = err2;

                err = q.error(opt_pos);

					opt_pos = p2;				

				}

                // Create SimplifyRec

			}

                SimplifyRec simplify_rec;

                simplify_rec.opt_pos = opt_pos;

			err = q.error(opt_pos);

                simplify_rec.err = err;

                simplify_rec.h = h;

			// Create SimplifyRec

                simplify_rec.time_stamp = halfedge_vec[h].time_stamp;

			int he_index = he->touched;

                simplify_rec.visits = 0;

			SimplifyRec simplify_rec;

                // push it.

			simplify_rec.opt_pos = opt_pos;

                sim_queue.push(simplify_rec);

			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.

        void QuadricSimplifier::update_onering_timestamp(VertexID v)

			sim_queue.push(simplify_rec);

        {

            // For all emanating edges h

		}

			for(Walker w = m.walker(v); !w.full_circle(); w = w.circulate_vertex_cw())

                push_simplify_rec(w.halfedge());

        }

		void QuadricSimplifier::update_onering_timestamp(VertexIter vi)

		{

        int QuadricSimplifier::collapse(SimplifyRec& simplify_rec)

			// For all emanating edges he

        { 

			for(VertexCirculator vc(vi);

            HalfEdgeID h = halfedge_vec[simplify_rec.h].h;

				!vc.end(); ++vc)

			{

			Walker w = m.walker(h);

				HalfEdgeIter he = vc.get_halfedge();

				update_time_stamp(he);

            // Check the time stamp to verify that the simplification 

				push_simplify_rec(he);

            // record is the newest. If the halfedge has been removed

			}

            // the time stamp is -1 and the comparison will also fail.

		}

            if(halfedge_vec[h].time_stamp == simplify_rec.time_stamp){

				Walker wo = w.opp();

		int QuadricSimplifier::collapse(SimplifyRec& simplify_rec)

                VertexID v = wo.vertex();

		{

                VertexID n = w.vertex();

			int he_index = simplify_rec.he_index;

			HalfEdgeIter he = halfedge_vec[he_index].he;

                // If the edge is, in fact, collapsible

                if(precond_collapse_edge(m, h)){      

			// Check the time stamp to verify that the simplification 

                    // If our consistency checks pass, we are relatively

			// record is the newest. If the halfedge has been removed

                    // sure that the contraction does not lead to a face flip.

			// the time stamp is -1 and the comparison will also fail.

                    if(check_consistency(h, simplify_rec.opt_pos) && check_consistency(wo.halfedge(), simplify_rec.opt_pos)){

			if(halfedge_vec[he_index].time_stamp == simplify_rec.time_stamp)

                        //cout << simplify_rec.err << " " << &(*he->vert) << endl;

			{

                        // Get QEM for both end points

				VertexIter v = he->opp->vert;

                        const QEM& Q1 = qem_vec[n];

				VertexIter n = he->vert;

                        const QEM& Q2 = qem_vec[v];

				// If the edge is, in fact, collapsible

                        // Compute Q_new = Q_1 + Q_2

				if(m.collapse_precond(he))

                        QEM q = Q1;

				{

                        q += Q2;

					// If our consistency checks pass, we are relatively

					// sure that the contraction does not lead to a face flip.

                        // Mark all halfedges that will be removed as dead

					if(check_consistency(he, simplify_rec.opt_pos) && 

                        halfedge_vec[w.halfedge()].halfedge_removed();

					   check_consistency(he->opp, simplify_rec.opt_pos))

                        halfedge_vec[wo.halfedge()].halfedge_removed();

					{

                        if(w.next().next().next().halfedge() == h){

                            halfedge_vec[w.next().halfedge()].halfedge_removed();

						// Get QEM for both end points

                            halfedge_vec[w.next().next().halfedge()].halfedge_removed();

						const QEM& Q1 = qem_vec[he->vert->touched];

                        }

						const QEM& Q2 = qem_vec[he->opp->vert->touched];

                        if(wo.next().next().next().halfedge() == wo.halfedge()){

                            halfedge_vec[wo.next().halfedge()].halfedge_removed();

						// Compute Q_new = Q_1 + Q_2

                            halfedge_vec[wo.next().next().halfedge()].halfedge_removed();

						QEM q = Q1;

                        }

						q += Q2;

                        // Do collapse

						// Mark all halfedges that will be removed as dead

                        m.collapse_edge(h);

						halfedge_vec[he->touched].halfedge_removed();

                        m.pos(n) = simplify_rec.opt_pos;

						halfedge_vec[he->opp->touched].halfedge_removed();

                        qem_vec[n] = q;

						if(he->next->next->next == he)

                        update_onering_timestamp(n);

						{

                        return 1;

							halfedge_vec[he->next->touched].halfedge_removed();

                    }

							halfedge_vec[he->next->next->touched].halfedge_removed();

                }

						}

						if(he->opp->next->next->next == he->opp)

						{

                // If we are here, the collapse was not allowed. If we have

							halfedge_vec[he->opp->next->touched].halfedge_removed();

                // seen this simplify record less than 100 times, we try to

							halfedge_vec[he->opp->next->next->touched].halfedge_removed();

                // increase the error and store the record again. Maybe some

						}

                // other contractions will make it more digestible later.

                if(simplify_rec.visits < 100){

						// Do collapse

                    simplify_rec.err *= 1.01f;

						m.collapse_halfedge(he,false);

                    ++simplify_rec.visits;

						n->pos = Vec3f(simplify_rec.opt_pos);

                    sim_queue.push(simplify_rec);

						qem_vec[n->touched] = q;

                }

            }

						update_onering_timestamp(n);

						return 1;

            return 0;

					}

        }

				}

				// If we are here, the collapse was not allowed. If we have

        void QuadricSimplifier::reduce(long int max_work)

				// seen this simplify record less than 100 times, we try to

        {

				// increase the error and store the record again. Maybe some

            // Set t = 0 for all halfedges

				// other contractions will make it more digestible later.

            for(HalfEdgeIDIterator h = m.halfedges_begin(); h != m.halfedges_end(); ++h){

				if(simplify_rec.visits < 100)

                halfedge_vec[*h].h = *h;

				{

            }

					simplify_rec.err *= 1.01f;

            cout << "Computing quadrics" << endl;

					++simplify_rec.visits;

					sim_queue.push(simplify_rec);

            // For all vertices, compute quadric and store in qem_vec

				}

            for(VertexIDIterator v = m.vertices_begin(); v != m.vertices_end(); ++v){

			}

                Vec3d p(m.pos(*v));

			return 0;

                Vec3d vn(normal(m, *v));

		}

                QEM q;

				for(Walker w = m.walker(*v); !w.full_circle(); w = w.circulate_vertex_cw()){

		void QuadricSimplifier::reduce(int max_work)

                    FaceID f = w.face();

		{

                    if(f != InvalidFaceID){

			// Set t = 0 for all halfedges

                        Vec3d n(normal(m, f));

			int i=0;

                        double a = area(m, f);

			for(HalfEdgeIter he = m.halfedges_begin();

                        q += QEM(p, n, a / 3.0);

				he != m.halfedges_end(); ++he, ++i)

                    }

				halfedge_vec[i].he = he;

                    if ((f == InvalidFaceID || w.opp().face() == InvalidFaceID ) && sqr_length(vn) > 0.0){

                        Vec3d edge = Vec3d(m.pos(w.vertex())) - p;

			cout << "Computing quadrics" << endl;

                        double edge_len = sqr_length(edge); 

                        if(edge_len > 0.0){

			// For all vertices, compute quadric and store in qem_vec

                            Vec3d n = cross(vn, edge);

			for(VertexIter vi=m.vertices_begin(); vi != m.vertices_end(); ++vi)

                            q += QEM(p, n, 2*edge_len);

			{

                        }

				Vec3d p(vi->pos);

                    }

				Vec3d vn(normal(vi));

                }

				QEM q;

                qem_vec[*v] = q;

				for(VertexCirculator vc(vi); !vc.end(); ++vc)

            }

				{

            cout << "Pushing initial halfedges" << endl;

					FaceIter f = vc.get_face();

					if(f != NULL_FACE_ITER)

            for(HalfEdgeIDIterator h = m.halfedges_begin(); h != m.halfedges_end(); ++h){

					{

                if(halfedge_vec[*h].time_stamp == 0)

						Vec3d n(normal(f));

                    push_simplify_rec(*h);

						double a = area(f);

            }

						q += QEM(p, n, a / 3.0);

            cout << "Simplify";

					}

					if ((f == NULL_FACE_ITER || 

            int work = 0;

						vc.get_opp_halfedge()->face ==  NULL_FACE_ITER) && sqr_length(vn) > 0.0)

            while(!sim_queue.empty() && work < max_work){

					{

                SimplifyRec simplify_record = sim_queue.top();

						Vec3d edge = Vec3d(vc.get_halfedge()->vert->pos)-p;

                sim_queue.pop();

						double edge_len = sqr_length(edge); 

						if(edge_len > 0.0)

                work += 2*collapse(simplify_record);

							{

                if((sim_queue.size() % 10000) == 0){

								Vec3d n = cross(vn, edge);

                    cout << ".";

								q += QEM(p, n, 2*edge_len);

                }

							}

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

					}

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

				}

            }

				qem_vec[vi->touched] = q;

            cout << endl;

			}

        }

			cout << "Pushing initial halfedges" << endl;

    }

			for(HalfEdgeIter he = m.halfedges_begin();

    void quadric_simplify(Manifold& m, double keep_fraction, double singular_thresh, bool choose_optimal_positions)

				he != m.halfedges_end(); ++he)

    {

				// For all halfedges, 

        gel_srand(1210);

			{

        long int F = m.no_faces();

				if(halfedge_vec[he->touched].time_stamp == 0)

        VertexAttributeVector<int> mask(m.no_faces(), 0);

				{

        long int max_work = max(static_cast<long int>(0), F- static_cast<long int>(keep_fraction * F));

					update_time_stamp(he);

        QuadricSimplifier qsim(m, mask, singular_thresh, choose_optimal_positions);

					push_simplify_rec(he);

        qsim.reduce(max_work);

				}

    }

			}

    void quadric_simplify(Manifold& m, VertexAttributeVector<int> mask, double keep_fraction, double singular_thresh, bool choose_optimal_positions)

			cout << "Simplify" << endl;

    {

        gel_srand(1210);

			int work = 0;

        long int F = m.no_faces();

			while(!sim_queue.empty() && work < max_work)

        long int max_work = keep_fraction == 0.0 ? INT_MAX : max(static_cast<long int>(0),

			{

                                                                 F- static_cast<long int>(keep_fraction * F));

				SimplifyRec simplify_record = sim_queue.top();

        QuadricSimplifier qsim(m, mask, singular_thresh, choose_optimal_positions);

				sim_queue.pop();

        qsim.reduce(max_work);

    }

				work += 2*collapse(simplify_record);

				if((work % 100) == 0)