Rev 630 | 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 "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(float a, float b) {return abs(a-b) < 0.00001;}
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;
}
}
}
void smooth_fun(Manifold& m,
VertexAttributeVector<int>& nailed,
VertexAttributeVector<double>& fun)
{
HalfEdgeAttributeVector<double> edge_weights;
FaceAttributeVector<int> included(m.allocated_faces(),1);
compute_edge_weights(m,edge_weights, included);
VertexAttributeVector<double> new_fun(m.no_vertices());
for(int i = 0; i < 1000; ++i)
{
for(VertexIDIterator v = m.vertices_begin(); v != m.vertices_end(); ++v)
if(!nailed[*v])
{
double w_sum = 0;
new_fun[*v] = 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]])
)
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;
}
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);
}
void polarize_mesh(Manifold& m, VertexAttributeVector<double>& fun, double vmin, double vmax, const int divisions, VertexAttributeVector<Vec2d>& parametrization)
{
vmax -= 0.01 * (vmax-vmin);
float interval = (vmax-vmin)/divisions;
VertexAttributeVector<int> status(m.allocated_vertices(), 0);
// ----------------------------------------
cout << "Tracing level set curves" << endl;
vector<HalfEdgeID> hidvec;
for(HalfEdgeIDIterator hid = m.halfedges_begin(); hid != m.halfedges_end(); ++hid)
{
Walker w = m.walker(*hid);
if(fun[w.vertex()] > fun[w.opp().vertex()])
hidvec.push_back(*hid);
}
for(size_t i = 0; i<hidvec.size(); ++i)
{
Walker w = m.walker(hidvec[i]);
float b = (fun[w.vertex()]- vmin)/interval;
float a = (fun[w.opp().vertex()] - vmin)/interval;
float floor_b = floor(b);
float floor_a = floor(a);
Vec3d pb = m.pos(w.vertex());
for(int j=floor_b; j>floor_a; --j)
{
float t = (j-a) / (b-a);
Vec3d p = t * pb + (1.0-t) * m.pos(w.opp().vertex());
VertexID v_new = m.split_edge(w.halfedge());
w = w.prev();
status[v_new] = 1;
fun[v_new] = j * interval + vmin;
m.pos(v_new) = p;
}
}
bool did_work;
do
{
did_work = false;
for(FaceIDIterator fid = m.faces_begin(); fid != m.faces_end(); ++fid)
for(Walker w = m.walker(*fid);!w.full_circle(); w = w.next())
if(status[w.vertex()] == 1 && !(status[w.next().vertex()]==1 && same_level(fun[w.vertex()],fun[w.next().vertex()]))
&& !(status[w.prev().vertex()]==1 && same_level(fun[w.vertex()],fun[w.prev().vertex()])))
{
Walker w0 = w;
w = w.next().next();
do
{
if(status[w.vertex()] == 1 && w.next().halfedge() != w0.halfedge() &&
same_level(fun[w0.vertex()],fun[w.vertex()]))
{
m.split_face_by_edge(*fid, w0.vertex(), w.vertex());
did_work = true;
break;
}
w = w.next();
}
while(!w.full_circle());
break;
}
}
while(did_work);
shortest_edge_triangulate(m);
// ----------------------------
cout << "Numbering the level sets" << endl;
VertexAttributeVector<int> nailed(m.no_vertices(),0);
vector<LevelSetInfo> level_set_info;
HalfEdgeAttributeVector<int> level_set_id(m.allocated_halfedges(), -1);
VertexAttributeVector<int> level_set_id_vertex(m.allocated_vertices(), -1);
VertexAttributeVector<Vec2d> circle_pos;
for(auto vid = m.vertices_begin(); vid != m.vertices_end(); ++vid)
circle_pos[*vid] = Vec2d(0,0);
int no_id=0;
for(HalfEdgeIDIterator hid = m.halfedges_begin(); hid != m.halfedges_end(); ++hid)
{
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;
float 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;
double param = 0;
do
{
double angle = 2.0 * M_PI * param / lsi.length;
circle_pos[w.opp().vertex()] = Vec2d(cos(angle), sin(angle));
nailed[w.opp().vertex()] = 1;
param += length(m, w.halfedge());
w = w.next();
while(status[w.vertex()] != 1 || !same_level(fun[w.vertex()], fun[w.opp().vertex()]))
w = w.circulate_vertex_cw();
}
while(w.halfedge() != lsi.h);
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;
}
}
smooth_fun(m, nailed, fun);
FaceAttributeVector<int> segmentation;
vector<vector<int>> boundaries;
segment_manifold(m, level_set_id, segmentation, boundaries);
VertexAttributeVector<Vec3d> cylinder_pos;
for(auto vid = m.vertices_begin();vid != m.vertices_end(); ++vid)
if (nailed[*vid])//level_set_id_vertex[*vid] == min_length_id)
cylinder_pos[*vid] = Vec3d(circle_pos[*vid][0], circle_pos[*vid][1],min_length_fun);
else
cylinder_pos[*vid] = Vec3d(0,0,fun[*vid]);
HalfEdgeAttributeVector<double> edge_weights(m.allocated_halfedges(), 0);
FaceAttributeVector<int> included(m.allocated_faces(), 1);
for(auto fid = m.faces_begin(); fid != m.faces_end(); ++fid)
{
if(boundaries[segmentation[*fid]].size() == 1)
included[*fid] = 0;
}
compute_edge_weights(m,edge_weights,included);
VertexAttributeVector<Vec3d> new_cylinder_pos(m.no_vertices());
for(int i = 0; i < 1500; ++i)
{
for(VertexIDIterator v = m.vertices_begin(); v != m.vertices_end(); ++v)
if(!nailed[*v])//level_set_id_vertex[*v] != min_length_id)
{
// recompute_circle_param(m, *v, fun, circle_pos, min_length_fun);
double w_sum = 0;
new_cylinder_pos[*v] = Vec3d(0);
for(Walker wv = m.walker(*v); !wv.full_circle(); wv = wv.circulate_vertex_ccw())
{
double w = edge_weights[wv.halfedge()];
new_cylinder_pos[*v] += w * cylinder_pos[wv.vertex()];
w_sum += w;
}
if(w_sum> 0.0)
new_cylinder_pos[*v] /= w_sum;
Vec2d v2d(new_cylinder_pos[*v][0], new_cylinder_pos[*v][1]);
double l =length(v2d);
if(l>0.1) {
new_cylinder_pos[*v][0] /= l;
new_cylinder_pos[*v][1] /= l;
}
if(level_set_id_vertex[*v] != -1)
new_cylinder_pos[*v][2] = fun[*v];
}
else
new_cylinder_pos[*v] = cylinder_pos[*v];
for(auto v = m.vertices_begin();v != m.vertices_end(); ++v)
cylinder_pos[*v] = 0.9*cylinder_pos[*v]+ 0.1*new_cylinder_pos[*v];
}
for(auto vid = m.vertices_begin();vid != m.vertices_end(); ++vid)
parametrization[*vid] = Vec2d(cylinder_pos[*vid][0],cylinder_pos[*vid][1]);
return;
// for(auto vid = m.vertices_begin();vid != m.vertices_end(); ++vid)
// m.pos(*vid) = cylinder_pos[*vid];// Vec2d(cylinder_pos[*vid][0],cylinder_pos[*vid][1]);
// vector<pair<double,VertexID>> verts;
// for(auto vid = m.vertices_begin();vid != m.vertices_end(); ++vid)
// if(!same_level(min_length_fun, fun[*vid]))
// verts.push_back(pair<double,VertexID>(abs(fun[*vid]-min_length_fun),*vid));
// sort(verts.begin(), verts.end());
//
// for(auto p:verts)
// {
// VertexID vid = p.second;
// recompute_circle_param(m, vid, fun, circle_pos, min_length_fun);
// }
// parametrization = circle_pos;
// return;
// // nail max length level set to circle.
// Walker w = m.walker(min_length_h);
// priority_queue<pair<double, HalfEdgeID>> hq;
// float param = 0;
// do
// {
// hq.push(pair<double,HalfEdgeID>(-abs(fun[w.next().vertex()]-min_length_fun),w.halfedge()));
// hq.push(pair<double,HalfEdgeID>(-abs(fun[w.opp().next().vertex()]-min_length_fun),w.opp().halfedge()));
//
// nailed[w.opp().vertex()] = 1;
// w = w.next();
// while(level_set_id[w.halfedge()] != min_length_id)
// w = w.circulate_vertex_cw();
// }
// while(w.halfedge() != min_length_h);
// parametrize(m, nailed, circle_pos, parametrization, fun, min_length_fun);
// propagate_values(m, nailed, circle_pos, parametrization, level_set_id_vertex, level_set_info, fun, min_length_fun);
//// for(int i=0;i<level_set_info.size(); ++i)
//// level_set_info[i].print();
//
// VertexAttributeVector<Vec3d> circle_center_attrib(m.no_vertices(),Vec3d(0,0,0));
// for(auto vid = m.vertices_begin(); vid != m.vertices_end(); ++vid)
// {
// if(level_set_id_vertex[*vid]>-1)
// {
// LevelSetInfo& lsi = level_set_info[level_set_id_vertex[*vid]];
// circle_center_attrib[*vid] = lsi.avg_pos;
// nailed[*vid] = 1;
// }
// }
//// propagate(m, nailed, circle_center_attrib);
// Manifold m2 = m;
// for(auto vid = m2.vertices_begin(); vid != m2.vertices_end(); ++vid)
// {
// m2.pos(*vid) = circle_pos[*vid];
// m2.pos(*vid)[2] = 0.01 * fun[*vid];
// }
// obj_save("blob.obj", m2);
// for(auto vid = m.vertices_begin(); vid != m.vertices_end(); ++vid)
// {
// Vec3d p =circle_pos[*vid];//-circle_center_attrib[*vid];
// parametrization[*vid] = Vec2d(p[0],p[1]);
// }
// ----------------------------
cout << "Remove vertices not on level set curves" << endl;
vector<VertexID> vid_vec;
for(VertexIDIterator vid = m.vertices_begin(); vid != m.vertices_end(); ++vid)
if(status[*vid]==0)
vid_vec.push_back(*vid);
random_shuffle(vid_vec.begin(), vid_vec.end());
for (size_t i=0; i<vid_vec.size(); ++i) {
FaceID f = m.merge_one_ring(vid_vec[i]);
if(f != InvalidFaceID)
shortest_edge_triangulate_face(m, f, level_set_id_vertex);
else
cout << "vertex not removed " << valency(m, vid_vec[i]) << endl;
}
for(FaceIDIterator fid = m.faces_begin(); fid != m.faces_end(); ++fid)
if(no_edges(m, *fid) > 3)
shortest_edge_triangulate_face(m, *fid, level_set_id_vertex);
// VertexAttributeVector<Vec3d> recalled_positions;
// for(VertexIDIterator vid = m.vertices_begin(); vid != m.vertices_end(); ++vid)
// recalled_positions[*vid] = m.pos(*vid);
//
//
TriangleQuality tq_energy(level_set_id_vertex);
priority_queue_optimization(m, tq_energy);
}
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)[2];
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);
}
}