Blame | Last modification | View Log | RSS feed
#include <cfloat>
#include "CGLA/Vec3d.h"
#include <queue>
#include "mesh_optimization.h"
#include "HMesh/FaceCirculator.h"
#include "HMesh/VertexCirculator.h"
using namespace std;
using namespace CGLA;
using namespace HMesh;
namespace HMeshUtil
{
// Small utility functions
namespace
{
class LineSeg
{
const CGLA::Vec3d p0;
const CGLA::Vec3d p1;
const float l;
const CGLA::Vec3d dir;
public:
LineSeg(const CGLA::Vec3d& _p0, const CGLA::Vec3d& _p1):
p0(_p0), p1(_p1), l(length(p1-p0)), dir((p1-p0)/l) {}
bool inseg(const CGLA::Vec3d& p) const
{
double t = dot(dir, p-p0);
if(t<0)
return false;
if(t>l)
return false;
return true;
}
};
Vec3d compute_normal(Vec3d* v)
{
Vec3d norm;
for(int i=0;i<4;++i)
{
norm[0] += (v[i][1]-v[(i+1)%4][1])*(v[i][2]+v[(i+1)%4][2]);
norm[1] += (v[i][2]-v[(i+1)%4][2])*(v[i][0]+v[(i+1)%4][0]);
norm[2] += (v[i][0]-v[(i+1)%4][0])*(v[i][1]+v[(i+1)%4][1]);
}
float l = norm.length();
if(l>0.0f)
norm /= l;
return norm;
}
bool would_flip(HalfEdgeIter h)
{
HalfEdgeIter ho = h->opp;
Vec3d v[4];
v[0]=Vec3d(h->vert->get_pos());
v[1]=Vec3d(ho->vert->get_pos());
v[2]=Vec3d(h->next->vert->get_pos());
v[3]=Vec3d(ho->next->vert->get_pos());
Vec3d dir(0,0,-1);//compute_normal(v);
Vec3d n1a = cross(v[3]-v[0], v[2]-v[0]);
Vec3d n2a = cross(v[2]-v[1], v[3]-v[1]);
if(dot(normalize(n1a),dir) < 0)
return true;
if(dot(normalize(n2a),dir) < 0)
return true;
return false;
}
}
// Energy functionals
namespace
{
class EnergyFun
{
public:
virtual double delta_energy(HalfEdgeIter) const = 0;
virtual double energy(HalfEdgeIter) const {return 0;}
};
class ValencyEnergy: public EnergyFun
{
public:
double delta_energy(HalfEdgeIter he) const
{
FaceIter f = he->face;
VertexIter v1 = he->opp->vert;
VertexIter v2 = he->vert;
VertexIter vo1 = he->next->vert;
VertexIter vo2 = he->opp->next->vert;
int val1 = valency(v1);
int val2 = valency(v2);
int valo1 = valency(vo1);
int valo2 = valency(vo2);
// The optimal valency is four for a boundary vertex
// and six elsewhere.
int t1 = v1->is_boundary() ? 4 : 6;
int t2 = v2->is_boundary() ? 4 : 6;
int to1 = vo1->is_boundary() ? 4 : 6;
int to2 = vo2->is_boundary() ? 4 : 6;
int before =
sqr(val1-t1)+sqr(val2-t2)+
sqr(valo1-to1)+sqr(valo2-to2);
int after =
sqr(valo1+1-to1)+sqr(val1-1-t1)+
sqr(val2-1-t2)+sqr(valo2+1-to2);
return static_cast<double>(after-before);
}
};
class RandomEnergy: public EnergyFun
{
public:
double delta_energy(HalfEdgeIter he) const
{
return static_cast<double>(rand()/static_cast<float>(RAND_MAX));
}
};
class MinAngleEnergy: public EnergyFun
{
double min_angle(const Vec3d& v0,const Vec3d& v1,const Vec3d& v2) const
{
Vec3d a = normalize(v1-v0);
Vec3d b = normalize(v2-v1);
Vec3d c = normalize(v0-v2);
return min(dot(a,-c), min(dot(b,-a), dot(c,-b)));
}
double thresh;
public:
MinAngleEnergy(double _thresh): thresh(_thresh) {}
double delta_energy(HalfEdgeIter he) const
{
FaceIter f = he->face;
FaceIter f2 = he->opp->face;
double t = dot(get_normal(f), get_normal(f2));
if(t> thresh || t < -thresh)
{
Vec3d v0(he->vert->get_pos());
Vec3d v1(he->next->vert->get_pos());
Vec3d v2(he->opp->vert->get_pos());
Vec3d v3(he->opp->next->vert->get_pos());
double before = min(min_angle(v0,v1,v2), min_angle(v0,v2,v3));
double after = min(min_angle(v0,v1,v3), min_angle(v1,v2,v3));
return -(after-before);
}
return DBL_MAX;
}
};
class DihedralEnergy: public EnergyFun
{
const double gamma;
const bool use_alpha;
double cos_ang(const Vec3d& n1, const Vec3d& n2) const
{
return s_max(-1.0, s_min(1.0, dot(n1, n2)));
}
double edge_alpha_energy(Vec3d v1, Vec3d v2, double ca) const
{
return pow(length(v1-v2)*(acos(ca)), 1.0f/gamma);
}
double edge_c_energy(Vec3d v1, Vec3d v2, double ca) const
{
return pow(length(v1-v2)*(1-ca), 1.0f/gamma);
}
mutable double ab_12;
mutable double ab_a1;
mutable double ab_b1;
mutable double ab_2c;
mutable double ab_2d;
mutable double aa_12;
mutable double aa_b1;
mutable double aa_c1;
mutable double aa_2a;
mutable double aa_2d;
public:
DihedralEnergy(double _gamma = 4.0, bool _use_alpha=false):
gamma(_gamma), use_alpha(_use_alpha) {}
void compute_angles(HalfEdgeIter he) const
{
Vec3d va(he->vert->get_pos());
Vec3d vb(he->opp->vert->get_pos());
Vec3d vc(he->next->vert->get_pos());
Vec3d vd(he->opp->next->vert->get_pos());
FaceIter fa = he->next->opp->face;
FaceIter fb = he->next->next->opp->face;
FaceIter fc = he->opp->next->opp->face;
FaceIter fd = he->opp->next->next->opp->face;
Vec3d n1 = normalize(cross(vc-va, vb-va));
Vec3d n2 = normalize(cross(vb-va, vd-va));
Vec3d na = fa==NULL_FACE_ITER? Vec3d(0) : Vec3d(get_normal(fa));
Vec3d nb = fb==NULL_FACE_ITER? Vec3d(0) : Vec3d(get_normal(fb));
Vec3d nc = fc==NULL_FACE_ITER? Vec3d(0) : Vec3d(get_normal(fc));
Vec3d nd = fd==NULL_FACE_ITER? Vec3d(0) : Vec3d(get_normal(fd));
Vec3d fn1 = normalize(cross(vb-vc, vd-vc));
Vec3d fn2 = normalize(cross(vd-vc, va-vc));
ab_12 = cos_ang(n1,n2);
ab_a1 = cos_ang(na,n1);
ab_b1 = cos_ang(nb,n1);
ab_2c = cos_ang(n2,nc);
ab_2d = cos_ang(n2,nd);
aa_12 = cos_ang(fn1,fn2);
aa_b1 = cos_ang(nb,fn1);
aa_c1 = cos_ang(nc, fn1);
aa_2a = cos_ang(fn2, na);
aa_2d = cos_ang(fn2,nd);
}
double energy(HalfEdgeIter he) const
{
double a = cos_ang(Vec3d(get_normal(he->face)),
Vec3d(get_normal(he->opp->face)));
Vec3d va(he->vert->get_pos());
Vec3d vb(he->opp->vert->get_pos());
if(use_alpha)
return edge_alpha_energy(va,vb,a);
return edge_c_energy(va,vb,a);
}
double delta_energy(HalfEdgeIter he) const
{
compute_angles(he);
Vec3d va(he->vert->get_pos());
Vec3d vb(he->opp->vert->get_pos());
Vec3d vc(he->next->vert->get_pos());
Vec3d vd(he->opp->next->vert->get_pos());
if(use_alpha)
{
double before =
edge_alpha_energy(va,vb,ab_12)
+edge_alpha_energy(va,vc,ab_a1)
+edge_alpha_energy(vc,vb,ab_b1)
+edge_alpha_energy(vd,vb,ab_2c)
+edge_alpha_energy(vd,va,ab_2d);
double after =
edge_alpha_energy(vd,vc,aa_12)
+edge_alpha_energy(vb,vc,aa_b1)
+edge_alpha_energy(vd,vb,aa_c1)
+edge_alpha_energy(va,vc,aa_2a)
+edge_alpha_energy(vd,va,aa_2d);
return after-before;
}
double before =
edge_c_energy(va,vb,ab_12)
+edge_c_energy(va,vc,ab_a1)
+edge_c_energy(vc,vb,ab_b1)
+edge_c_energy(vd,vb,ab_2c)
+edge_c_energy(vd,va,ab_2d);
double after =
edge_c_energy(vd,vc,aa_12)
+edge_c_energy(vb,vc,aa_b1)
+edge_c_energy(vd,vb,aa_c1)
+edge_c_energy(va,vc,aa_2a)
+edge_c_energy(vd,va,aa_2d);
return after-before;
}
double min_angle(HalfEdgeIter he) const
{
compute_angles(he);
return s_min(s_min(s_min(s_min(aa_12, aa_b1), aa_c1), aa_2a), aa_2d);
}
};
class CurvatureEnergy: public EnergyFun
{
// double abs_mean_curv(const Vec3d& v, const vector<Vec3d>& ring) const
// {
// const int N = ring.size();
// vector<Vec3d> norms(N,Vec3d(0));
// for(int i=0;i<N; ++i)
// {
// const Vec3d& vn = ring[i];
// const Vec3d& vnn = ring[(i+1)%N];
// Vec3d n = cross(vnn-v,vn-v);
// double ln = length(vn);
// if(ln> 1e-20)
// norms[i] = n/ln;
// }
// double H = 0.0f;
// for(int i=0;i<N; ++i)
// {
// const Vec3d& vn = ring[i];
// double e_len = length(v-vn);
// double beta =
// acos(max(-1.0, min(1.0, dot(norms[i], norms[(i+N-1)%N]))));
// H += e_len * beta;
// }
// return H/4.0f;
// }
double abs_mean_curv(const Vec3d& v, const vector<Vec3d>& ring) const
{
const int N = ring.size();
double H = 0;
for(int i=0;i<N;++i)
{
Vec3d vnim1 = ring[(i+N-1)%N] - v;
Vec3d vni = ring[i] - v;
Vec3d vnip1 = ring[(i+1)%N] - v;
Vec3d Nm = normalize(cross(vni, vnim1));
Vec3d Np = normalize(cross(vnip1, vni));
double beta = acos(max(-1.0, min(1.0, dot(Nm, Np))));
H += vni.length() * beta;
}
H /= 4;
return H;
}
public:
double delta_energy(HalfEdgeIter he) const
{
VertexIter va = he->vert;
VertexIter vb = he->opp->vert;
VertexIter vc = he->next->vert;
VertexIter vd = he->opp->next->vert;
Vec3d va_pos(va->get_pos());
Vec3d vb_pos(vb->get_pos());
Vec3d vc_pos(vc->get_pos());
Vec3d vd_pos(vd->get_pos());
vector<Vec3d> va_ring_bef;
vector<Vec3d> va_ring_aft;
vector<Vec3d> vb_ring_bef;
vector<Vec3d> vb_ring_aft;
vector<Vec3d> vc_ring_bef;
vector<Vec3d> vc_ring_aft;
vector<Vec3d> vd_ring_bef;
vector<Vec3d> vd_ring_aft;
for(VertexCirculator c(va); !c.end(); ++c)
{
va_ring_bef.push_back(Vec3d(c.get_vertex()->get_pos()));
if(c.get_vertex() != vb)
va_ring_aft.push_back(Vec3d(c.get_vertex()->get_pos()));
}
for(VertexCirculator c(vb); !c.end(); ++c)
{
vb_ring_bef.push_back(Vec3d(c.get_vertex()->get_pos()));
if(c.get_vertex() != va)
vb_ring_aft.push_back(Vec3d(c.get_vertex()->get_pos()));
}
for(VertexCirculator c(vc); !c.end(); ++c)
{
vc_ring_bef.push_back(Vec3d(c.get_vertex()->get_pos()));
vc_ring_aft.push_back(Vec3d(c.get_vertex()->get_pos()));
if(c.get_vertex() == va)
vc_ring_aft.push_back(Vec3d(vd->get_pos()));
}
for(VertexCirculator c(vd); !c.end(); ++c)
{
vd_ring_bef.push_back(Vec3d(c.get_vertex()->get_pos()));
vd_ring_aft.push_back(Vec3d(c.get_vertex()->get_pos()));
if(c.get_vertex() == vb)
vd_ring_aft.push_back(Vec3d(vc->get_pos()));
}
double before =
abs_mean_curv(va_pos,va_ring_bef) +
abs_mean_curv(vb_pos,vb_ring_bef) +
abs_mean_curv(vc_pos,vc_ring_bef) +
abs_mean_curv(vd_pos,vd_ring_bef);
double after =
abs_mean_curv(va_pos,va_ring_aft) +
abs_mean_curv(vb_pos,vb_ring_aft) +
abs_mean_curv(vc_pos,vc_ring_aft) +
abs_mean_curv(vd_pos,vd_ring_aft);
return after-before;
}
};
class GaussCurvatureEnergy: public EnergyFun
{
double gauss_curv(const Vec3d& v, const vector<Vec3d>& ring) const
{
const int N = ring.size();
double asum=0.0f;
double area_sum=0;
for(int i=0;i<N; ++i)
{
const Vec3d& v1 = ring[i];
const Vec3d& v2 = ring[(i+1)%N];
Vec3d a = v1-v;
Vec3d b = v2-v;
asum += acos(max(-1.0, min(1.0, dot(a,b)/(length(a)*length(b)))));
area_sum += 0.5 * length(cross(a,b));
}
return 3*abs(2 * M_PI - asum)/area_sum;
}
public:
double delta_energy(HalfEdgeIter he) const
{
VertexIter va = he->vert;
VertexIter vb = he->opp->vert;
VertexIter vc = he->next->vert;
VertexIter vd = he->opp->next->vert;
Vec3d va_pos(va->get_pos());
Vec3d vb_pos(vb->get_pos());
Vec3d vc_pos(vc->get_pos());
Vec3d vd_pos(vd->get_pos());
vector<Vec3d> va_ring_bef;
vector<Vec3d> va_ring_aft;
vector<Vec3d> vb_ring_bef;
vector<Vec3d> vb_ring_aft;
vector<Vec3d> vc_ring_bef;
vector<Vec3d> vc_ring_aft;
vector<Vec3d> vd_ring_bef;
vector<Vec3d> vd_ring_aft;
for(VertexCirculator c(va); !c.end(); ++c)
{
va_ring_bef.push_back(Vec3d(c.get_vertex()->get_pos()));
if(c.get_vertex() != vb)
va_ring_aft.push_back(Vec3d(c.get_vertex()->get_pos()));
}
for(VertexCirculator c(vb); !c.end(); ++c)
{
vb_ring_bef.push_back(Vec3d(c.get_vertex()->get_pos()));
if(c.get_vertex() != va)
vb_ring_aft.push_back(Vec3d(c.get_vertex()->get_pos()));
}
for(VertexCirculator c(vc); !c.end(); ++c)
{
vc_ring_bef.push_back(Vec3d(c.get_vertex()->get_pos()));
vc_ring_aft.push_back(Vec3d(c.get_vertex()->get_pos()));
if(c.get_vertex() == va)
vc_ring_aft.push_back(Vec3d(vd->get_pos()));
}
for(VertexCirculator c(vd); !c.end(); ++c)
{
vd_ring_bef.push_back(Vec3d(c.get_vertex()->get_pos()));
vd_ring_aft.push_back(Vec3d(c.get_vertex()->get_pos()));
if(c.get_vertex() == vb)
vd_ring_aft.push_back(Vec3d(vc->get_pos()));
}
double before =
gauss_curv(va_pos,va_ring_bef) +
gauss_curv(vb_pos,vb_ring_bef) +
gauss_curv(vc_pos,vc_ring_bef) +
gauss_curv(vd_pos,vd_ring_bef);
double after =
gauss_curv(va_pos,va_ring_aft) +
gauss_curv(vb_pos,vb_ring_aft) +
gauss_curv(vc_pos,vc_ring_aft) +
gauss_curv(vd_pos,vd_ring_aft);
return after-before;
}
};
struct PQElem
{
double pri;
HalfEdgeIter he;
int time;
PQElem() {}
PQElem(double _pri, HalfEdgeIter _he, int _time):
pri(_pri), he(_he), time(_time) {}
};
bool operator<(const PQElem& e0, const PQElem& e1)
{
return e0.pri > e1.pri;
}
void add_to_queue(priority_queue<PQElem>& Q,
HalfEdgeIter he,
const EnergyFun& efun)
{
if(!he->is_boundary())
{
double energy = efun.delta_energy(he);
int t = he->touched + 1;
he->touched = t;
he->opp->touched = t;
if((energy<0) && (t < 10000))
Q.push(PQElem(energy, he, t));
}
}
void add_one_ring_to_queue(priority_queue<PQElem>& Q,
VertexIter vi,
const EnergyFun& efun)
{
for(VertexCirculator vc(vi); !vc.end(); ++vc)
{
add_to_queue(Q, vc.get_halfedge(), efun);
}
}
void priority_queue_optimization(Manifold& m, const EnergyFun& efun)
{
priority_queue<PQElem> Q;
{
HalfEdgeIter he = m.halfedges_begin();
for(;he != m.halfedges_end(); ++he)
he->touched = 0;
he = m.halfedges_begin();
for(;he != m.halfedges_end(); ++he)
if(he->touched == 0)
add_to_queue(Q, he, efun);
}
double pri_old = -DBL_MAX;
while(!Q.empty())
{
PQElem elem = Q.top();
Q.pop();
HalfEdgeIter he = elem.he;
if(he->touched == elem.time)
{
if(m.flip(he))
{
add_one_ring_to_queue(Q, he->vert, efun);
add_one_ring_to_queue(Q, he->next->vert, efun);
add_one_ring_to_queue(Q, he->opp->vert, efun);
add_one_ring_to_queue(Q, he->opp->next->vert, efun);
}
}
}
}
void simulated_annealing_optimization(Manifold& m,
const EnergyFun& efun,
int max_iter=10000)
{
srand(0);
int swaps;
int iter = 0;
double T = 1;
double tmin=0;
double tsum=0;
{
HalfEdgeIter he = m.halfedges_begin();
for(;he != m.halfedges_end(); ++he)
{
if(!he->is_boundary())
{
double e = efun.delta_energy(he);
tmin = s_min(e, tmin);
tsum += efun.energy(he);
}
}
}
cout << "Initial energy:" << tsum << endl;
T = -2*tmin;
if(max_iter>0) do
{
vector<HalfEdgeIter> halfedges;
HalfEdgeIter he = m.halfedges_begin();
for(;he != m.halfedges_end(); ++he)
{
if(!he->is_boundary())
halfedges.push_back(he);
}
random_shuffle(halfedges.begin(), halfedges.end());
swaps = 0;
for(int i=0;i<halfedges.size();++i)
{
he = halfedges[i];
// if(!would_flip(he))
DihedralEnergy dih_en;
if(dih_en.min_angle(he) > -0.4)
{
double delta = efun.delta_energy(he);
if(delta < -1e-8)
{
if(m.flip(he))
++swaps;
}
else
{
delta = max(1e-8, delta);
double prob = min(0.9999, exp(-delta/T));
if(rand()/double(RAND_MAX) < prob)
{
if(m.flip(he))
++swaps;
}
}
}
}
cout << "Swaps = " << swaps << " T = " << T << endl;
if(iter % 5 == 0 && iter > 0)
T *= 0.9;
}
while(++iter < max_iter && swaps);
cout << "Iterations " << iter << endl;
}
}
void minimize_dihedral_angle(Manifold& m,
int iter,
bool anneal,
bool alpha,
double gamma)
{
DihedralEnergy energy_fun(gamma, alpha);
if(anneal)
simulated_annealing_optimization(m, energy_fun, iter);
else
priority_queue_optimization(m, energy_fun);
}
void randomize_mesh(Manifold& m, int max_iter)
{
RandomEnergy energy_fun;
simulated_annealing_optimization(m, energy_fun, max_iter);
}
void minimize_curvature(Manifold& m, bool anneal)
{
CurvatureEnergy energy_fun;
if(anneal)
simulated_annealing_optimization(m, energy_fun);
else
priority_queue_optimization(m, energy_fun);
}
void minimize_gauss_curvature(Manifold& m, bool anneal)
{
GaussCurvatureEnergy energy_fun;
if(anneal)
simulated_annealing_optimization(m, energy_fun);
else
priority_queue_optimization(m, energy_fun);
}
void maximize_min_angle(Manifold& m, float thresh, bool anneal)
{
MinAngleEnergy energy_fun(thresh);
if(anneal)
simulated_annealing_optimization(m, energy_fun);
else
priority_queue_optimization(m, energy_fun);
}
void optimize_valency(Manifold& m, bool anneal)
{
ValencyEnergy energy_fun;
if(anneal)
simulated_annealing_optimization(m, energy_fun);
else
priority_queue_optimization(m, energy_fun);
}
}