Rev 640 | Go to most recent revision | Blame | Last modification | View Log | RSS feed
//
// polarize.cpp
// GEL
//
// Created by J. Andreas Bærentzen on 18/03/12.
// Copyright 2012 __MyCompanyName__. All rights reserved.
//
#include <queue>
#include <iomanip>
#include "harmonics.h"
#include "polarize.h"
#include <CGLA/Vec2d.h>
#include <LinAlg/LapackFunc.h>
#include <HMesh/triangulate.h>
#include <HMesh/obj_save.h>
#include <HMesh/curvature.h>
#include <HMesh/quadric_simplify.h>
#include <HMesh/mesh_optimization.h>
#include <HMesh/smooth.h>
#include <complex>
using namespace CGLA;
using namespace std;
using namespace HMesh;
inline bool same_level(double a, double b) {return abs(a-b) < 1e-10;}
enum VertexStatus {
RAW, POLE, SADDLE, COOKED
};
void kill_wedges(Manifold& m, const VertexAttributeVector<double>& fun, const VertexAttributeVector<VertexStatus>& status)
{
for(auto hid : m.halfedges())
if(m.in_use(hid)){
Walker w = m.walker(hid);
if(same_level(fun[w.vertex()], fun[w.opp().vertex()]) &&
status[w.vertex()] == COOKED &&
status[w.opp().vertex()] == COOKED &&
no_edges(m, w.face())==3 &&
status[w.next().vertex()] == COOKED &&
precond_collapse_edge(m, hid))
m.collapse_edge(hid);
}
}
void high_gradient_triangulate(Manifold& m, FaceID _f, const VertexAttributeVector<double>& fun)
{
queue<FaceID> fq;
fq.push(_f);
while(!fq.empty()) {
FaceID f = fq.front();
fq.pop();
if(no_edges(m, f)<=3)
continue;
// Create a vector of vertices.
vector<VertexID> verts;
circulate_face_ccw(m, f, [&](VertexID v){verts.push_back(v);});
// Find vertex pairs that may be connected.
vector<pair<int,int> > vpairs;
const int N = verts.size();
for(int i = 0; i < N - 2; ++i){
for(int j = i + 2; j < N; ++j){
if(verts[i] != verts[j] && !connected(m, verts[i], verts[j]))
vpairs.push_back(pair<int,int>(i, j));
}
}
if(vpairs.empty()){
cout << "Warning: could not triangulate a face."
<< "Probably a vertex appears more than one time in other vertex's one-ring" << endl;
continue;
}
/* For all vertex pairs, find the edge lengths. Combine the
vertices forming the shortest edge. */
float max_grad=-FLT_MAX;
size_t max_k = -1;
for(size_t k = 0; k < vpairs.size(); ++k){
int i = vpairs[k].first;
int j = vpairs[k].second;
double grad = abs(fun[verts[i]]-fun[verts[j]])/length(m.pos(verts[i]) - m.pos(verts[j]));
if(grad>max_grad){
max_grad = grad;
max_k = k;
}
}
assert(max_k != -1);
{
int i = vpairs[max_k].first;
int j = vpairs[max_k].second;
if(m.in_use(verts[i]) && m.in_use(verts[j])) {
FaceID f2 = m.split_face_by_edge(f, verts[i], verts[j]);
fq.push(f);
fq.push(f2);
}
}
}
}
FaceID merge_face_in_1_ring(Manifold& m, VertexID v)
{
Walker w=m.walker(v);
HalfEdgeID h = w.next().halfedge();
if(m.remove_vertex(v))
{
return m.close_hole(h);
}
return InvalidFaceID;
}
void compute_edge_weights(const Manifold& m, HalfEdgeAttributeVector<double>& edge_weights, FaceAttributeVector<int>& included)
{
edge_weights = HalfEdgeAttributeVector<double>(m.allocated_halfedges(), 0);
for(FaceIDIterator f = m.faces_begin(); f != m.faces_end(); ++f)
if(included[*f])
{
for(Walker wv = m.walker(*f); !wv.full_circle(); wv = wv.circulate_face_ccw())
{
HalfEdgeID h = wv.halfedge();
Vec3d p1(m.pos(wv.vertex()));
Vec3d p2(m.pos(wv.next().vertex()));
Vec3d p0(m.pos(wv.opp().vertex()));
double ang = acos(min(1.0, max(-1.0, dot(normalize(p1-p0), normalize(p2-p0)))));
double ang_opp = acos(min(1.0, max(-1.0, dot(normalize(p2-p1), normalize(p0-p1)))));
double l = (p1-p0).length();
edge_weights[h] += tan(ang/2) / l;
edge_weights[wv.opp().halfedge()] += tan(ang_opp/2) / l;
}
}
}
template<typename T>
void smooth_fun(const Manifold& m,
VertexAttributeVector<int>& nailed,
VertexAttributeVector<T>& fun, int iter=100)
{
HalfEdgeAttributeVector<double> edge_weights;
FaceAttributeVector<int> included(m.allocated_faces(),1);
compute_edge_weights(m,edge_weights, included);
VertexAttributeVector<T> new_fun(m.no_vertices());
for(int i = 0; i < iter; ++i)
{
for(VertexIDIterator v = m.vertices_begin(); v != m.vertices_end(); ++v)
if(!nailed[*v])
{
double w_sum = 0;
new_fun[*v] = T(0);
for(Walker wv = m.walker(*v); !wv.full_circle(); wv = wv.circulate_vertex_ccw())
{
double w = edge_weights[wv.halfedge()];
new_fun[*v] += w * fun[wv.vertex()];
w_sum += w;
}
new_fun[*v] /= w_sum;
}
else
new_fun[*v] = fun[*v];
fun = new_fun;
}
}
vector<VertexID> cut_mesh(Manifold& m, VertexAttributeVector<double>& fun, double cut_val, VertexAttributeVector<VertexStatus>& status)
{
cout << "Cutting at " << cut_val << endl;
vector<VertexID> new_verts;
vector<HalfEdgeID> hidvec;
for(HalfEdgeID hid : m.halfedges())
{
Walker w = m.walker(hid);
double bval = fun[w.vertex()];
double aval = fun[w.opp().vertex()];
if(aval<bval && aval < cut_val && cut_val <= bval)
{
double t = (cut_val-aval)/(bval-aval);
Vec3d p = t*m.pos(w.vertex()) + (1.0-t)*m.pos(w.opp().vertex());
VertexID vnew = m.split_edge(hid);
m.pos(vnew) = p;
status[vnew] = COOKED;
fun[vnew] = cut_val;
new_verts.push_back(vnew);
}
}
for(FaceID fid : m.faces())
for(Walker w = m.walker(fid);!w.full_circle(); w = w.next())
if(status[w.vertex()] == COOKED && same_level(fun[w.vertex()],cut_val))
{
Walker w0 = w;
w = w.next().next();
do
{
if(status[w.vertex()] == COOKED &&
w.next().halfedge() != w0.halfedge() &&
same_level(fun[w.vertex()],cut_val))
{
m.split_face_by_edge(fid, w0.vertex(), w.vertex());
break;
}
w = w.next();
}
while(!w.full_circle());
break;
}
return new_verts;
}
void smooth_cooked_vertices(Manifold& m, const vector<VertexID>& vertices, VertexAttributeVector<double>& fun,
VertexAttributeVector<VertexStatus>& status, int iter=1)
{
for(int i=0;i<iter;++i)
{
VertexAttributeVector<Vec3d> npos(m.allocated_vertices(), Vec3d(0));
VertexAttributeVector<int> cnt(m.allocated_vertices(), 0);
for(auto vid : vertices)
if(m.in_use(vid))
for(Walker w = m.walker(vid); !w.full_circle(); w = w.circulate_vertex_ccw())
if(status[vid]==COOKED) {
npos[vid] += m.pos(w.vertex());
++cnt[vid];
}
for(auto vid : vertices)
if(m.in_use(vid))
m.pos(vid) = (npos[vid])/cnt[vid];
}
}
int ssb(double a, double b)
{
const double THRESH = 1e-6;
double delta = a-b;
if(delta<-THRESH) return -1;
if(delta>THRESH) return 1;
return 0;
}
void classify_vertices(Manifold& m, VertexAttributeVector<double>& fun, VertexAttributeVector<VertexStatus>& status)
{
for(auto vid : m.vertices())
{
double val = fun[vid];
auto w = m.walker(vid);
int above_below = ssb(fun[w.vertex()],val);
int changes = 0;
w = w.circulate_vertex_ccw();
for(; !w.full_circle();w = w.circulate_vertex_ccw())
{
int above_below_new = ssb(fun[w.vertex()],val);
if(above_below != above_below_new)
++changes;
above_below=above_below_new;
}
if(changes==0)
status[vid]=POLE;
else if(changes>=4)
status[vid]=SADDLE;
else if(changes==2)
status[vid]=RAW;
else assert(0);
}
}
void remove_vertices(Manifold& m, VertexAttributeVector<double>& fun, VertexAttributeVector<VertexStatus>& status, double rmin, double rmax)
{
vector<VertexID> vertices_to_push;
for(auto vid: m.vertices())
if(status[vid]==RAW &&
fun[vid] >= rmin &&
fun[vid] <= rmax)
vertices_to_push.push_back(vid);
random_shuffle(begin(vertices_to_push), end(vertices_to_push));
for(auto vid : vertices_to_push)
{
FaceID fid = merge_face_in_1_ring(m,vid);
if(fid != InvalidFaceID)
high_gradient_triangulate(m, fid, fun);
}
}
void polarize_mesh(Manifold& m, VertexAttributeVector<double>& fun, double vmin, double vmax, const int divisions, VertexAttributeVector<Vec2d>& parametrization)
{
VertexAttributeVector<VertexStatus> status(m.allocated_vertices(), RAW);
classify_vertices(m, fun, status);
double interval = (vmax-vmin)/(divisions+1);
double base_cut_val = vmin + 0.02*(vmax-vmin);
vector<VertexID> new_verts = cut_mesh(m, fun, base_cut_val, status);
smooth_cooked_vertices(m, new_verts, fun, status, 1);
cout << "Base " << base_cut_val << endl;
for(int iter=0;iter<divisions+2;++iter)
{
int _iter=iter;//%2==0?iter/2:-1-iter/2;
double cut_val = base_cut_val + interval * _iter;
cout << "i,c " << _iter << "," << cut_val << endl;
remove_vertices(m, fun, status, cut_val, cut_val+interval);
vector<VertexID> new_verts = cut_mesh(m, fun,_iter<0?cut_val:cut_val+interval, status);
kill_wedges(m, fun, status);
smooth_cooked_vertices(m, new_verts, fun, status, 1);
}
int C=0, R=0, P=0,S=0;
for(auto vid : m.vertices())
{
switch(status[vid]) {
case RAW: ++R; break;
case COOKED: ++C; break;
case POLE: ++P; break;
case SADDLE: ++S; break;
}
}
cout << "R=" << R << " C=" << C << " P=" << P << " S=" << S << endl;
}
void make_height_fun(const HMesh::Manifold& m, HMesh::VertexAttributeVector<double>& fun,
double& vmin, double& vmax)
{
VertexIDIterator vid = m.vertices_begin();
VertexAttributeVector<int> nailed(m.allocated_vertices(), 0);
smooth_fun(m, nailed, fun,0);
vmin = vmax = m.pos(*vid)[1];
for(; vid != m.vertices_end(); ++vid)
{
double v = dot(m.pos(*vid),Vec3d(0.0,1,0.00));
fun[*vid] = v;
vmin = min(v, vmin);
vmax = max(v, vmax);
}
}
void make_adf_fun(HMesh::Manifold& m, double t, HMesh::VertexAttributeVector<double>& F,
double& vmin, double& vmax)
{
Harmonics harm(m);
harm.compute_adf(F, t);
VertexAttributeVector<int> nailed(m.allocated_vertices(), 0);
vmin = vmax = F[*m.vertices_begin()];
for(auto vid : m.vertices())
{
vmin = min(F[vid], vmin);
vmax = max(F[vid], vmax);
}
cout << vmin << vmax << endl;
}