Rev 642 | 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 <CGLA/Mat4x4d.h>
#include <CGLA/Vec3i.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 <HMesh/dual.h>
#include <Geometry/GridAlgorithm.h>
#include <Geometry/build_bbtree.h>
#include <Geometry/RGrid.h>
#include <Geometry/TrilinFilter.h>
#include <GLGraphics/ManifoldRenderer.h>
#include <complex>
using namespace CGLA;
using namespace Geometry;
using namespace std;
using namespace HMesh;
using namespace GLGraphics;
inline bool same_level(double a, double b) {return abs(a-b) < 1e-10;}
enum VertexStatus {
RAW, COOKED
};
void shortest_edge_triangulate_face(Manifold& m, FaceID f0, const VertexAttributeVector<int>& level_set_id_vertex)
{
queue<FaceID> face_queue;
face_queue.push(f0);
while(!face_queue.empty())
{
FaceID f = face_queue.front();
face_queue.pop();
// Create a vector of vertices.
vector<VertexID> verts;
for(Walker w = m.walker(f); !w.full_circle(); w = w.circulate_face_ccw())
{
FaceID fa = w.face();
FaceID fb = f;
assert(fa==fb);
verts.push_back(w.vertex());
}
// If there are just three we are done.
if(verts.size() == 3) continue;
// 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]) &&
(level_set_id_vertex[verts[i]] == 0 || level_set_id_vertex[verts[i]] != level_set_id_vertex[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 min_len=FLT_MAX;
int min_k = -1;
for(size_t k = 0; k < vpairs.size(); ++k){
int i = vpairs[k].first;
int j = vpairs[k].second;
float len = length(m.pos(verts[i]) - m.pos(verts[j]));
// cout << m.pos(verts[i]) << m.pos(verts[j]) << len << endl;
if(len<min_len){
min_len = len;
min_k = k;
}
}
assert(min_k != -1);
if(min_k>=0)
{
// Split faces along edge whose midpoint is closest to isovalue
int i = vpairs[min_k].first;
int j = vpairs[min_k].second;
// cout << i << " " << j << " " << min_k << " " << vpairs.size() << endl;
FaceID f_new = m.split_face_by_edge(f, verts[i], verts[j]);
if(no_edges(m, f)>3)
face_queue.push(f);
if(no_edges(m, f_new)>3)
face_queue.push(f_new);
}
else
cout << "failed to triangle. Suspect NaN vertex positions!!" << endl;
}
}
void flip_edges(Manifold& m, VertexAttributeVector<int> ls_id)
{
vector<HalfEdgeID> hvec;
for(HalfEdgeID hid : m.halfedges())
{
Walker w = m.walker(hid);
if(hid< w.opp().halfedge() &&
ls_id[w.vertex()] > 0 &&
ls_id[w.opp().vertex()] > 0 &&
ls_id[w.vertex()] != ls_id[w.opp().vertex()])
hvec.push_back(hid);
}
ValencyEnergy mae;
bool did_work = false;
int iter=0;
do {
did_work = false;
random_shuffle(begin(hvec), end(hvec));
for(auto hid : hvec)
{
Walker w =m.walker(hid);
if(mae.delta_energy(m, hid)<=0 &&
ls_id[w.next().vertex()] != ls_id[w.opp().next().vertex()]) {
m.flip_edge(hid);
did_work = true;
}
}
}
while (did_work && ++iter<1000);
cout << "flip iter" << iter << endl;
}
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;
}
}
void label_connected_components(Manifold&m, int& cur_id, const vector<VertexID>& vertices,
const VertexAttributeVector<double>& fun,
const VertexAttributeVector<VertexStatus>& status,
VertexAttributeVector<int>& level_id)
{
for(auto vid : vertices)
{
if(level_id[vid]>0)
continue;
queue<VertexID> vq;
vq.push(vid);
while(!vq.empty())
{
VertexID v = vq.front();
vq.pop();
level_id[v] = cur_id;
for(Walker w = m.walker(v); !w.full_circle(); w = w.circulate_vertex_ccw())
{
if(status[w.vertex()]==COOKED && same_level(fun[v], fun[w.vertex()]))
if(level_id[w.vertex()]==0)
vq.push(w.vertex());
else if(level_id[w.vertex()] != cur_id)
cout << "previously labeled vertex encountered" << endl;
}
}
++cur_id;
}
}
void connect_touched(Manifold& m, const VertexAttributeVector<int>& touched)
{
for(FaceID fid : m.faces())
for(Walker w = m.walker(fid);!w.full_circle(); w = w.next())
if(touched[w.vertex()])
{
Walker w0 = w;
w = w.next().next();
do
{
if(touched[w.vertex()] && w.next().halfedge() != w0.halfedge())
{
m.split_face_by_edge(fid, w0.vertex(), w.vertex());
break;
}
w = w.next();
}
while(!w.full_circle());
break;
}
}
void remove_shorts(Manifold& m, VertexAttributeVector<int> ls_id)
{
double avghlen=0;
for(auto hid : m.halfedges())
avghlen += length(m, hid);
avghlen /= m.no_halfedges();
for(auto hid : m.halfedges())
if(m.in_use(hid))
{
Walker w = m.walker(hid);
if(ls_id[w.vertex()] > 0 &&
ls_id[w.vertex()] == ls_id[w.opp().vertex()] &&
length(m, w.halfedge()) < 0.1 * avghlen &&
precond_collapse_edge(m, hid))
m.collapse_edge(hid);
}
}
vector<VertexID> separate_edges(Manifold& m, const vector<VertexID>& vertices, int id,
VertexAttributeVector<double>& fun,
VertexAttributeVector<VertexStatus>& status,
VertexAttributeVector<int>& ls_id)
{
VertexAttributeVector<int> touched(m.allocated_vertices(),0);
vector<VertexID> new_vertices;
// Cut all edges to vertices on other LS ID
for(auto vid: vertices)
if(ls_id[vid]==id)
{
vector<HalfEdgeID> edges_to_split;
for(Walker w = m.walker(vid); !w.full_circle(); w = w.circulate_vertex_ccw())
{
int id = ls_id[w.vertex()];
if(id > 0 & id != ls_id[vid])
edges_to_split.push_back(w.halfedge());
}
for (auto hid : edges_to_split)
{
VertexID vnew = m.split_edge(hid);
fun[vnew] = fun[vid];
m.pos(vnew) = m.pos(vid);
status[vnew] = COOKED;
ls_id[vnew] = ls_id[vid];
touched[vnew] = 1;
new_vertices.push_back(vnew);
}
}
// Label old vertices as plain (uncooked)
for(auto vid: vertices)
if(ls_id[vid]==id)
{
status[vid] = RAW;
ls_id[vid] = 0;
}
connect_touched(m, touched);
return new_vertices;
}
vector<VertexID> cut_mesh(Manifold& m, VertexAttributeVector<double>& fun, double cut_val,
VertexAttributeVector<VertexStatus>& status)
{
cout << "cutting @ " << cut_val << endl;
vector<VertexID> new_verts;
vector<HalfEdgeID> hidvec;
VertexAttributeVector<int> touched(m.allocated_vertices(), 0);
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;
touched[vnew] = 1;
new_verts.push_back(vnew);
}
}
connect_touched(m, touched);
return new_verts;
}
void smooth_cooked_vertices(Manifold& m, const vector<VertexID>& vertices,
VertexAttributeVector<VertexStatus>& status, int iter=7)
{
for(int i=0;i<iter;++i)
{
VertexAttributeVector<Vec3d> npos = m.positions_attribute_vector();
VertexAttributeVector<int> cnt(m.allocated_vertices(), 1);
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];
}
}
void remove_v5_vertices(Manifold& m, VertexAttributeVector<int>& ls_id, int id_ceiling)
{
for(auto vid: m.vertices())
{
int id = ls_id[vid];
if(id>0 && id<id_ceiling)
{
vector<HalfEdgeID> hit_list;
for(Walker w = m.walker(vid); !w.full_circle(); w=w.circulate_vertex_ccw())
if(ls_id[w.vertex()] >= id_ceiling)
hit_list.push_back(w.halfedge());
if(hit_list.size()>1)
{
sort(begin(hit_list),end(hit_list),
[&m](HalfEdgeID a, HalfEdgeID b){return length(m,a)>length(m, b);});
for(int i=1;i<hit_list.size();++i)
{
HalfEdgeID hid = hit_list[i];
Walker w = m.walker(hid);
HalfEdgeID surplus_h = w.opp().prev().halfedge();
if(precond_collapse_edge(m,surplus_h))
{
m.collapse_edge(surplus_h,true);
m.merge_faces(m.walker(hid).face(), hid);
}
}
}
}
}
}
void remove_wedges(Manifold& m, VertexAttributeVector<int>& ls_id)
{
for(auto hid : m.halfedges())
if(m.in_use(hid))
{
Walker w = m.walker(hid);
if(ls_id[w.vertex()]>0 &&
ls_id[w.vertex()] == ls_id[w.opp().vertex()] &&
ls_id[w.next().vertex()]>0 &&
no_edges(m, w.face())==3 &&
precond_collapse_edge(m, hid))
m.collapse_edge(hid, true);
}
}
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 remove_vertices(Manifold& m, VertexAttributeVector<double>& fun,
const VertexAttributeVector<VertexStatus>& status,
const VertexAttributeVector<int>& level_id, 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 = m.merge_one_ring(vid);
if(fid != InvalidFaceID)
shortest_edge_triangulate_face(m, fid, level_id);
}
}
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);
VertexAttributeVector<int> level_id(m.allocated_vertices(), 0);
int cur_id = 1;
double interval = (vmax-vmin)/(divisions+1);
double cut_val = vmin + 0.065*(vmax-vmin);
vector<VertexID> cut_verts = cut_mesh(m, fun, cut_val, status);
label_connected_components(m, cur_id, cut_verts, fun, status, level_id);
remove_vertices(m, fun, status, level_id, vmin, cut_val);
flip_edges(m, level_id);
for(int id=1;id<cur_id;++id) {
vector<VertexID> loop_verts = separate_edges(m, cut_verts, id, fun, status, level_id);
smooth_cooked_vertices(m, loop_verts, status, 50);
}
for(int iter=0;iter<divisions+2;++iter)
{
cut_val += interval;
vector<VertexID> cut_verts = cut_mesh(m, fun,cut_val, status);
int old_id = cur_id;
label_connected_components(m, cur_id, cut_verts, fun, status, level_id);
remove_vertices(m, fun, status, level_id, cut_val-1.01*interval, cut_val);
vector<VertexID> loop_verts;
for(int id=old_id;id<cur_id;++id)
{
vector<VertexID> tmp = separate_edges(m, cut_verts, id, fun, status, level_id);
loop_verts.insert(loop_verts.end(), tmp.begin(),tmp.end());
}
remove_wedges(m, level_id);
remove_v5_vertices(m, level_id, old_id);
smooth_cooked_vertices(m, loop_verts, status);
}
dual(m);
m.cleanup();
}
void simplify_polar_mesh(Manifold& m, double frac)
{
int N_orig = m.no_vertices();
HalfEdgeAttributeVector<int> touched(m.allocated_halfedges(),0);
vector<pair<double,vector<HalfEdgeID>>> work_vector;
for(auto hid : m.halfedges())
{
Walker w = m.walker(hid);
if(!touched[hid] &&
no_edges(m, w.face()) == 4 &&
no_edges(m, w.opp().face()) == 4)
{
double l_sum=0;
int n=0;
vector<HalfEdgeID> hvec;
bool abort=false;
for(int i=0;i<2;++i)
{
Walker w0 = (i==0?w:w.opp().next().next().opp());
while(!touched[w0.halfedge()])
{
touched[w0.halfedge()] = 1;
touched[w0.opp().halfedge()] = 1;
hvec.push_back(w0.halfedge());
++n;
l_sum += sqr(area(m,w0.face()));
if(valency(m,w0.vertex()) != 4 || valency(m, w0.opp().vertex()) != 4)
abort = true;
if(no_edges(m, w0.face())==4)
w0 = w0.next().next().opp();
else {
if (valency(m,w0.next().vertex())<5)
abort = true;
break;
}
}
}
if(abort)
continue;
else
work_vector.push_back(make_pair(l_sum/sqr(n),hvec));
}
}
sort(begin(work_vector),end(work_vector));
for(auto work : work_vector)
{
for(HalfEdgeID h :work.second)
if(m.in_use(h) && precond_collapse_edge(m, h))
m.collapse_edge(h,true);
else cout << "Collapse failed ..." << endl;
if(N_orig * frac > m.no_vertices())
break;
}
m.cleanup();
}
using namespace Geometry;
const Mat4x4d fit_bounding_volume(const Vec3d& p0,
const Vec3d& p7,
float buf_reg,
int vol_dim)
{
Vec3d sz = p7 - p0;
Vec3i dims(vol_dim);
Vec3d scal_vec = (Vec3d(dims)-Vec3d(2*buf_reg+2))/sz;
float scal = min(scal_vec[0], min(scal_vec[1],scal_vec[2]));
Mat4x4d m = translation_Mat4x4d(Vec3d(0)+Vec3d(buf_reg+1));
m *= scaling_Mat4x4d(Vec3d(scal));
m *= translation_Mat4x4d(-p0);
return m;
}
void smooth_and_refit(HMesh::Manifold& m_in, HMesh::Manifold& m_ref, int iter, int vol_dim)
{
AABBTree tree;
Manifold m_tmp = m_ref;
Vec3d p0,p7;
bbox(m_tmp, p0, p7);
Mat4x4d mat = fit_bounding_volume(p0, p7, 3, vol_dim);
Mat4x4d imat = invert(mat);
for(auto vid : m_tmp.vertices())
m_tmp.pos(vid) = mat.mul_3D_point(m_tmp.pos(vid));
build_AABBTree(m_tmp, tree);
cout << "Done building AABB tree" << endl;
RGrid<float> grid(Vec3i(vol_dim),FLT_MAX);
auto F = [&](const Vec3i& p,float& v)
{
v = tree.compute_signed_distance(Vec3f(p));
};
for_each_voxel(grid, F);
cout << "Made distance field" << endl;
TrilinFilter<RGrid<float>> df(&grid);
for (int i=0;i<iter; ++i)
{
auto vsafe = m_in.positions_attribute_vector();
laplacian_smooth(m_in,.25,10);
for(auto vid : m_in.vertices())
{
double w = exp(4-valency(m_in, vid));
m_in.pos(vid) *= w;
m_in.pos(vid) += (1-w) * vsafe[vid];
}
float avg = 0;
for(auto hid : m_in.halfedges())
avg += length(m_in, hid);
avg /= m_in.no_halfedges();
auto npos = m_in.positions_attribute_vector();
for(auto vid : m_in.vertices())
{
npos[vid] += 0.5 * avg * normal(m_in, vid);
}
m_in.positions_attribute_vector() = npos;
for(auto vid : m_in.vertices())
{
Vec3f pos = Vec3f(mat.mul_3D_point(m_in.pos(vid)));
if(df.in_domain(pos))
m_in.pos(vid) = imat.mul_3D_point(Vec3d(pos - df(pos) * df.grad(pos)));
else
cout << "outside domain" << 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;
}