Rev 631 | 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 "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;}
struct LevelSetInfo
{
int id;
int no_vertices;
Vec3d avg_pos;
double rad;
double length;
int components;
HalfEdgeID h;
double fun_value;
void print()
{
cout
<< " id : " << id
<< " no vertices : " << no_vertices
<< " avg pos : " << avg_pos
<< " rad : " << rad
<< " length : " << length
<< " # : " << components
<< " F : " << fun_value << endl;
}
};
void compute_edge_weights(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(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 segment_manifold(Manifold& m, HalfEdgeAttributeVector<int>& ls_id,
FaceAttributeVector<int>& segmentation,
vector<vector<int>>& boundaries)
{
segmentation = FaceAttributeVector<int>(m.no_faces(), -1);
int SEG_NO = 0;
for(auto f = m.faces_begin(); f != m.faces_end(); ++f)
{
if (segmentation[*f] == -1)
{
queue<FaceID> q;
q.push(*f);
vector<int> bound(0);
while (!q.empty())
{
FaceID face = q.front();
q.pop();
segmentation[face] = SEG_NO;
for(Walker w = m.walker(face); !w.full_circle(); w = w.next())
{
if(ls_id[w.halfedge()] == -1)
{
FaceID fopp = w.opp().face();
if(segmentation[fopp] == -1) {
q.push(fopp);
}
}
else {
bound.push_back(ls_id[w.halfedge()]);
}
}
}
sort(begin(bound), end(bound));
auto new_end = unique(begin(bound), end(bound));
boundaries.push_back(vector<int>(begin(bound),new_end));
cout << "Boundaries of seg# " << SEG_NO << " : ";
for(auto be:boundaries[SEG_NO])
cout << be << " ";
cout << endl;
SEG_NO += 1;
}
}
}
void shortest_edge_triangulate_face(Manifold& m, FaceID f0, 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]]) &&
valency(m, verts[i])<5 && valency(m,verts[j])<5
)
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 = sqr_length(m.pos(verts[i]) - m.pos(verts[j]));
if(len<min_len){
min_len = len;
min_k = k;
}
}
assert(min_k != -1);
{
// Split faces along edge whose midpoint is closest to isovalue
int i = vpairs[min_k].first;
int j = vpairs[min_k].second;
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);
}
}
}
void shortest_edge_split_face(Manifold& m, FaceID f0, 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())
{
verts.push_back(w.vertex());
}
// Find vertex pairs that may be connected.
vector<pair<int,int> > vpairs;
const int N = verts.size();
for(int i = 0; i < N ; ++i){
for(int j = 3; j < N-2; ++j){
int jj = (j+i)%N;
if(!connected(m, verts[i], verts[jj]) &&
(level_set_id_vertex[verts[i]] != level_set_id_vertex[verts[jj]]) &&
(level_set_id_vertex[verts[(i+1)%N]] == level_set_id_vertex[verts[i]]) &&
(level_set_id_vertex[verts[(i+N-1)%N]] == level_set_id_vertex[verts[i]]) &&
(level_set_id_vertex[verts[(jj+1)%N]] == level_set_id_vertex[verts[jj]]) &&
(level_set_id_vertex[verts[(jj+N-1)%N]] == level_set_id_vertex[verts[jj]]))
vpairs.push_back(pair<int,int>(i, jj));
}
}
if(vpairs.empty()) 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 = sqr_length(m.pos(verts[i]) - m.pos(verts[j]));
if(len<min_len){
min_len = len;
min_k = k;
}
}
assert(min_k != -1);
{
// Split faces along edge whose midpoint is closest to isovalue
int i = vpairs[min_k].first;
int j = vpairs[min_k].second;
FaceID f_new = m.split_face_by_edge(f, verts[i], verts[j]);
if(no_edges(m, f)>5)
face_queue.push(f);
if(no_edges(m, f_new)>5)
face_queue.push(f_new);
}
}
}
struct EdgeQElem {
float priority;
HalfEdgeID he;
EdgeQElem(float _priority, HalfEdgeID _he): priority(_priority), he(_he) {}
};
bool operator<(const EdgeQElem& l, const EdgeQElem& r)
{
return l.priority < r.priority;
}
class FunctionalDifference: public EnergyFun
{
VertexAttributeVector<float>& fun;
VertexAttributeVector<int>& status;
public:
FunctionalDifference(VertexAttributeVector<float>& _fun, VertexAttributeVector<int>& _status): fun(_fun), status(_status) {}
virtual double delta_energy(const Manifold& m, HalfEdgeID h) const
{
Walker w = m.walker(h);
if(status[w.vertex()] == 1 && status[w.opp().vertex()]==1)
return DBL_MAX;
return sqr_length(m.pos(w.next().vertex())-m.pos(w.opp().next().vertex()))/(1e-6+abs(fun[w.next().vertex()]-fun[w.opp().next().vertex()])) - sqr_length(m.pos(w.vertex())-m.pos(w.opp().vertex()))/(1e-6+abs(fun[w.vertex()]-fun[w.opp().vertex()]));
}
};
class TriangleQualityValence: public EnergyFun
{
VertexAttributeVector<int>& idv;
MinAngleEnergy mae;
ValencyEnergy vae;
public:
TriangleQualityValence(VertexAttributeVector<int>& _idv): idv(_idv), mae(-1) {}
virtual double delta_energy(const Manifold& m, HalfEdgeID h) const
{
Walker w = m.walker(h);
if(idv[w.next().vertex()] == idv[w.opp().next().vertex()] ||
idv[w.vertex()] == idv[w.opp().vertex()])
return DBL_MAX;
VertexID v1 = w.opp().vertex();
VertexID v2 = w.vertex();
VertexID vo1 = w.next().vertex();
VertexID vo2 = w.opp().next().vertex();
int val1 = valency(m, v1);
int val2 = valency(m, v2);
int valo1 = valency(m, vo1);
int valo2 = valency(m, vo2);
// The optimal valency is four for a boundary vertex
// and six elsewhere.
int t1 = boundary(m, v1) ? 4 : 6;
int t2 = boundary(m, v2) ? 4 : 6;
int to1 = boundary(m, vo1) ? 4 : 6;
int to2 = boundary(m, vo2) ? 4 : 6;
int before =
max(max(sqr(val1-t1),sqr(val2-t2)), max(sqr(valo1-to1), sqr(valo2-to2)));
int after =
max(max(sqr(valo1+1-to1),sqr(val1-1-t1)), max(sqr(val2-1-t2),sqr(valo2+1-to2)));
return static_cast<double>(after-before);
// return vae.delta_energy(m,h);
float la = length(m.pos(w.next().vertex())-m.pos(w.opp().next().vertex()));
float lb = length(m.pos(w.vertex())-m.pos(w.opp().vertex()));
return la-lb;
return mae.delta_energy(m,h);
}
};
class TriangleQuality: public EnergyFun
{
VertexAttributeVector<int>& idv;
MinAngleEnergy mae;
ValencyEnergy vae;
public:
TriangleQuality(VertexAttributeVector<int>& _idv): idv(_idv), mae(-1) {}
virtual double delta_energy(const Manifold& m, HalfEdgeID h) const
{
Walker w = m.walker(h);
if(idv[w.next().vertex()] == idv[w.opp().next().vertex()] ||
idv[w.vertex()] == idv[w.opp().vertex()])
return DBL_MAX;
return mae.delta_energy(m,h);
}
};
Vec3d grad(HMesh::Manifold& m, HMesh::VertexAttributeVector<double>& fun, HMesh::FaceID f)
{
if(no_edges(m,f) != 3)
return Vec3d(0.0);
Vec3d n = normal(m, f);
Vec3d gsum(0.0);
for(Walker w = m.walker(f); !w.full_circle(); w = w.next())
{
VertexID avid = w.opp().vertex();
VertexID bvid = w.vertex();
VertexID cvid = w.next().vertex();
Vec3d gdir = normalize(cross(n, m.pos(bvid) - m.pos(avid)));
double l = dot(m.pos(cvid)-m.pos(avid), gdir);
gdir *= fun[cvid]/l;
gsum += gdir;
}
return gsum;
}
Vec2d extend_fun2(HMesh::Manifold& m, HMesh::HalfEdgeID h,
HMesh::VertexAttributeVector<double>& fun,
HMesh::VertexAttributeVector<Vec2d>& fun2)
{
Walker w = m.walker(h);
FaceID f = w.face();
Vec3d a = m.pos(w.opp().vertex());
Vec3d b = m.pos(w.vertex());
Vec3d c = m.pos(w.next().vertex());
Vec3d g = grad(m, fun, f);
Vec3d n = normal(m, f);
Vec3d N = normalize(cross(g, n));
float dot_nba = dot(N,b-a);
Vec2d cval;
if(dot_nba > 1e-5)
{
float t = dot(N, c-a)/dot_nba;
if(t>0 && t < 1) {
Vec2d aval = fun2[w.opp().vertex()];
Vec2d bval = fun2[w.vertex()];
return (1-t)*aval + t*bval;
}
}
return Vec2d(0);
}
double solve_for_orthogonal_gradients(HMesh::Manifold& m, HMesh::VertexAttributeVector<double>& fun,
HalfEdgeID h, double beta, double gamma)
{
Vec3d n = normal(m, m.walker(h).face());
Vec3d uvw[3];
Vec3d abc;
Walker w = m.walker(h);
VertexID vid[] = {w.next().vertex(), w.opp().vertex(),w.vertex()};
for(int i=0; !w.full_circle(); w = w.next(),++ i)
{
uvw[i] = normalize(cross(n, m.pos(vid[(i+1)%3]) - m.pos(vid[(i+2)%3])));
uvw[i] /= dot(m.pos(vid[i])-m.pos(vid[(i+1)%3]), uvw[i]);
abc[i] = fun[vid[i]];
}
Vec3d A(dot(uvw[0],uvw[0]),dot(uvw[0],uvw[1]),dot(uvw[0],uvw[2]));
Vec3d B(dot(uvw[1],uvw[0]),dot(uvw[1],uvw[1]),dot(uvw[1],uvw[2]));
Vec3d C(dot(uvw[2],uvw[0]),dot(uvw[2],uvw[1]),dot(uvw[2],uvw[2]));
return - (beta * dot(B, abc) + gamma * dot(C, abc))/dot(A, abc);
}
void recompute_circle_param(Manifold& m, VertexID v, VertexAttributeVector<double>& fun, VertexAttributeVector<Vec2d>& circle_pos, double f_origin)
{
Vec2d c_new(0);
for(auto w = m.walker(v); !w.full_circle(); w = w.circulate_vertex_ccw())
{
Walker w_base = w.next();
Vec2d X = circle_pos[w_base.opp().vertex()];
Vec2d Y(-X[1], X[0]);
Vec2d vec = circle_pos[w_base.vertex()];
double angle = atan2(max(-1.0,min(1.0,dot(vec,Y))),
max(-1.0,min(1.0,dot(vec,X))));
double new_angle = solve_for_orthogonal_gradients(m, fun, w_base.halfedge(), 0, angle);
new_angle = max(-M_PI/4, min(M_PI/4,new_angle));
c_new += cos(new_angle)*X + sin(new_angle)*Y;
}
circle_pos[v] = normalize(c_new);
}
template<typename T>
void orthogonal_trajectories(Manifold& m, VertexID v, VertexAttributeVector<double>& f, double min_length_fun,
VertexAttributeVector<T>& h_in, VertexAttributeVector<T>& h_out)
{
h_out[v] = T(0);
T num(0);
double denom = 0.0001;
for(Walker w = m.walker(v); !w.full_circle(); w = w.circulate_vertex_ccw())
{
double weight = area(m, w.face());
Vec3d a = m.pos(v);
Vec3d b = m.pos(w.vertex());
Vec3d c = m.pos(w.next().vertex());
Vec3d fvec(f[v],f[w.vertex()], f[w.next().vertex()]);
if(sqr_length(fvec)>0)
fvec.normalize();
double f_a = fvec[0];
double f_b = fvec[1];
double f_c = fvec[2];
T h_b = h_in[w.vertex()];
T h_c = h_in[w.next().vertex()];
Vec3d e_a = c-b;
Vec3d e_b = a-c;
Vec3d e_c = b-a;
Vec3d N = normalize(cross(e_b,e_c)+cross(e_c,e_a)+cross(e_a,e_b));
Vec3d g_a = cross(N, normalize(e_a));
g_a /= dot(g_a, e_b);
Vec3d g_b = cross(N, normalize(e_b));
g_b /= dot(g_b, e_c);
Vec3d g_c = cross(N, normalize(e_c));
g_c /= dot(g_c, e_a);
num -= weight * (h_b * (f_a*dot(g_a, g_b) + f_b*dot(g_b,g_b) + f_c*dot(g_b,g_c)) +
h_c * (f_a*dot(g_a, g_c) + f_b*dot(g_b,g_c) + f_c*dot(g_c, g_c)));
denom += weight * (f_a * dot(g_a,g_a)+ f_b*dot(g_a,g_b) +f_c*dot(g_a, g_c));
}
h_out[v] = num / denom;
}
void polarize_mesh(Manifold& m, VertexAttributeVector<double>& fun, double vmin, double vmax, const int divisions, VertexAttributeVector<Vec2d>& parametrization)
{
double interval = (vmax-vmin)/(divisions+1);
VertexAttributeVector<int> status(m.allocated_vertices(), 0);
// ----------------------------------------
cout << "Tracing level set curves" << endl;
for(int cut_no=1; cut_no<divisions+1;++cut_no)
{
double cut_val = vmin + cut_no * interval;
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] = 1;
fun[vnew] = cut_val;
}
}
for(FaceID fid : m.faces())
for(Walker w = m.walker(fid);!w.full_circle(); w = w.next())
if(status[w.vertex()] == 1 && same_level(fun[w.vertex()],cut_val))
{
Walker w0 = w;
w = w.next().next();
do
{
if(status[w.vertex()] == 1 &&
w.next().halfedge() != w0.halfedge() &&
same_level(fun[w0.vertex()],cut_val))
{
m.split_face_by_edge(fid, w0.vertex(), w.vertex());
break;
}
w = w.next();
}
while(!w.full_circle());
break;
}
}
shortest_edge_triangulate(m);
// ----------------------------
cout << "Numbering the level sets" << endl;
vector<LevelSetInfo> level_set_info;
HalfEdgeAttributeVector<int> level_set_id(m.allocated_halfedges(), -1);
VertexAttributeVector<int> level_set_id_vertex(m.allocated_vertices(), -1);
int no_id=0;
for(HalfEdgeID hid : m.halfedges())
{
Walker w = m.walker(hid);
if(status[w.vertex()] == 1 && status[w.opp().vertex()] == 1 &&
same_level(fun[w.vertex()], fun[w.opp().vertex()]) &&
level_set_id[w.halfedge()] == -1)
{
LevelSetInfo lsi;
lsi.id = no_id;
lsi.fun_value = fun[w.vertex()];
lsi.components = 1;
lsi.no_vertices = 0;
double level_set_length = 0;
while(level_set_id[w.halfedge()] != no_id)
{
level_set_length += length(m,w.halfedge());
level_set_id[w.halfedge()] = no_id;
level_set_id[w.opp().halfedge()] = no_id;
level_set_id_vertex[w.vertex()] = no_id;
w = w.next();
while(!(status[w.vertex()] == 1 && same_level(fun[w.vertex()], fun[w.opp().vertex()])))
w = w.circulate_vertex_cw();
lsi.no_vertices += 1;
}
lsi.h = w.halfedge();
lsi.length = level_set_length;
level_set_info.push_back(lsi);
level_set_info.back().print();
++no_id;
assert(level_set_info.size()==no_id);
}
}
cout << "ids" << no_id << " " << level_set_info.size() << endl;
for(int i=0;i<level_set_info.size(); ++i)
for(int j=0;j<level_set_info.size(); ++j)
if(i != j && same_level(level_set_info[i].fun_value, level_set_info[j].fun_value))
level_set_info[i].components += 1;
HalfEdgeID min_length_h;
double min_length_fun;
double min_length=FLT_MAX;
int min_length_id = -1;
for(int i=0;i<level_set_info.size(); ++i)
{
LevelSetInfo& lsi = level_set_info[i];
if(lsi.components==1 && lsi.length<min_length)
{
min_length_h = lsi.h;
min_length_fun = lsi.fun_value;
min_length = lsi.length;
min_length_id = lsi.id;
}
}
// ----------------------------
cout << "Remove vertices not on level set curves" << endl;
}
void make_height_fun(const HMesh::Manifold& m, HMesh::VertexAttributeVector<double>& fun,
double& vmin, double& vmax)
{
VertexIDIterator vid = m.vertices_begin();
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);
}
}