Blame | Last modification | View Log | RSS feed
#include <queue>
#include <iostream>
#include "quadric_simplify.h"
#include "mesh_optimization.h"
#include "smooth.h"
#include "HMesh/VertexCirculator.h"
#include "Geometry/QEM.h"
namespace HMeshUtil
{
using namespace CGLA;
using namespace std;
using namespace GEO;
using namespace HMesh;
namespace
{
struct HalfEdgeRec
{
HalfEdgeIter he;
int time_stamp;
void halfedge_removed() {time_stamp = -1;}
HalfEdgeRec(): time_stamp(0) {}
};
struct SimplifyRec
{
Vec3f opt_pos;
int he_index;
float err;
int time_stamp;
};
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;
void push_simplify_rec(HalfEdgeIter he);
void put_onering_in_queue(VertexIter vi);
int collapse(const SimplifyRec& simplify_rec);
public:
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();
}
void reduce(int max_work);
};
void QuadricSimplifier::push_simplify_rec(HalfEdgeIter he)
{
// Update time_stamp
assert(halfedge_vec[he->touched].time_stamp >=0);
assert(halfedge_vec[he->opp->touched].time_stamp>=0);
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;
// 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;
Vec3f opt_pos(0);
//cout << q.determinant() << endl;
if (q.determinant()>= 1e-8)
{
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);
Vec3f v0 = he->vert->get_pos();
Vec3f v1 = he->opp->vert->get_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 = time_stamp;
// push it.
sim_queue.push(simplify_rec);
}
void QuadricSimplifier::put_onering_in_queue(VertexIter vi)
{
// For all emanating edges he
for(VertexCirculator vc(vi);
!vc.end(); ++vc)
{
HalfEdgeIter he = vc.get_halfedge();
push_simplify_rec(he);
}
}
int QuadricSimplifier::collapse(const SimplifyRec& simplify_rec)
{
int he_index = simplify_rec.he_index;
// 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)
{
HalfEdgeIter he = halfedge_vec[he_index].he;
VertexIter v = he->opp->vert;
VertexIter n = he->vert;
// If the edge is, in fact, collapsible
if(m.collapse_precond(v,he))
{
//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(v,he,false);
n->set_pos(simplify_rec.opt_pos);
qem_vec[n->touched] = q;
put_onering_in_queue(n);
return 1;
}
}
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;
// For all vertices, compute quadric and store in qem_vec
for(VertexIter vi=m.vertices_begin(); vi != m.vertices_end(); ++vi)
{
Vec3d p(vi->get_pos());
QEM q;
for(VertexCirculator vc(vi); !vc.end(); ++vc)
{
FaceIter f = vc.get_face();
if(f != NULL_FACE_ITER)
{
Vec3d n(get_normal(f));
double a = get_area(f);
q += QEM(p, n, a / 3.0);
}
}
qem_vec[vi->touched] = q;
}
for(HalfEdgeIter he = m.halfedges_begin();
he != m.halfedges_end(); ++he)
// For all halfedges,
{
if(halfedge_vec[he->touched].time_stamp == 0)
push_simplify_rec(he);
}
int work = 0;
while(!sim_queue.empty() && work < max_work)
{
int dwork = collapse(sim_queue.top());
if(dwork>0)
{
work += dwork;
if((work % 100) == 0)
cout << "work = " << work << endl;
}
sim_queue.pop();
}
}
}
void quadric_simplify(Manifold& m, int max_work)
{
srand(1210);
QuadricSimplifier qsim(m);
qsim.reduce(max_work);
}
}