Rev 572 | Go to most recent revision | Blame | Compare with Previous | 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 <HMesh/triangulate.h>
#include <HMesh/curvature.h>
#include <HMesh/quadric_simplify.h>
#include <HMesh/mesh_optimization.h>
#include <HMesh/smooth.h>
using namespace CGLA;
using namespace std;
using namespace HMesh;
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(HalfEdgeWalker w = m.halfedgewalker(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(HalfEdgeWalker w = m.halfedgewalker(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
{
HalfEdgeWalker w = m.halfedgewalker(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 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
{
HalfEdgeWalker w = m.halfedgewalker(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);
}
};
inline bool same_level(float a, float b) {return abs(a-b) < 0.00001;}
Vec3f laplacian(const Manifold& m, VertexID v)
{
Vec3f avg_pos(0);
float asum = 0;
for(HalfEdgeWalker w = m.halfedgewalker(v); !w.full_circle(); w = w.circulate_vertex_cw()){
float a = barycentric_area(m, w.vertex());
avg_pos += a * m.pos(w.vertex());
asum += a;
}
return avg_pos / asum - m.pos(v);
}
void aw_laplacian_smooth(Manifold& m, float t)
{
VertexAttributeVector<Vec3f> pos(m.total_vertices());
for(VertexIDIterator v = m.vertices_begin(); v != m.vertices_end(); ++v){
if(!boundary(m, *v))
{
Vec3f n = normal(m, *v);
Vec3f l = laplacian(m, *v);
if(sqr_length(n) > 0.8)
l -= n * dot(n,l);
pos[*v] = t * l + m.pos(*v);
}
}
for(VertexIDIterator v = m.vertices_begin(); v != m.vertices_end(); ++v){
if(!boundary(m, *v))
m.pos(*v) = pos[*v];
}
}
void polarize_mesh(Manifold& m, VertexAttributeVector<float>& fun, float vmin, float vmax, const int divisions)
{
float interval = (vmax-vmin)/divisions;
VertexAttributeVector<int> status(m.total_vertices(), 0);
// ----------------------------------------
cout << "Tracing level set curves" << endl;
vector<HalfEdgeID> hidvec;
for(HalfEdgeIDIterator hid = m.halfedges_begin(); hid != m.halfedges_end(); ++hid)
{
HalfEdgeWalker w = m.halfedgewalker(*hid);
if(fun[w.vertex()] > fun[w.opp().vertex()])
hidvec.push_back(*hid);
}
for(int i = 0; i<hidvec.size(); ++i)
{
HalfEdgeWalker w = m.halfedgewalker(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);
Vec3f pb = m.pos(w.vertex());
for(int j=floor_b; j>floor_a; --j)
{
float t = (j-a) / (b-a);
Vec3f 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(HalfEdgeWalker w = m.halfedgewalker(*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()])))
{
HalfEdgeWalker 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);
// ----------------------------
cout << "Numbering the level sets" << endl;
HalfEdgeAttributeVector<int> level_set_id(m.total_halfedges(), 0);
VertexAttributeVector<int> level_set_id_vertex(m.total_vertices(), 0);
int no_id=1;
for(HalfEdgeIDIterator hid = m.halfedges_begin(); hid != m.halfedges_end(); ++hid)
{
HalfEdgeWalker w = m.halfedgewalker(*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()] == 0)
{
while(level_set_id[w.halfedge()] != no_id)
{
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();
}
++no_id;
}
}
cout << "Number of level sets : " << (no_id-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 (int 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<Vec3f> 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);
for(int iter=0;iter< 100;++iter)
{
//// ----------------------------
cout << "smooth level set curves" << endl;
aw_laplacian_smooth(m,0.025);
VertexAttributeVector<Vec3f> new_pos(m.total_vertices(), Vec3f(0));
VertexAttributeVector<float> new_pos_wsum(m.total_vertices(), 0.0);
for(HalfEdgeIDIterator hid = m.halfedges_begin(); hid != m.halfedges_end(); ++hid)
{
HalfEdgeWalker w = m.halfedgewalker(*hid);
if(level_set_id[w.halfedge()] != 0)
{
float weight = 1;//(valency(m,w.opp().vertex()));
new_pos[w.vertex()] += weight*m.pos(w.opp().vertex());
new_pos_wsum[w.vertex()] += weight;
}
}
for(VertexIDIterator vid = m.vertices_begin(); vid != m.vertices_end(); ++vid)
{
float weight = 1;//(valency(m,*vid));
new_pos[*vid] += weight*m.pos(*vid);
new_pos_wsum[*vid] += weight;
m.pos(*vid) = new_pos[*vid] / new_pos_wsum[*vid];
}
}
TriangleQuality tq_energy(level_set_id_vertex);
priority_queue_optimization(m, tq_energy);
// vector<HalfEdgeID> hidvec;
// for(HalfEdgeIDIterator hid = m.halfedges_begin(); hid != m.halfedges_end(); ++hid)
// {
// HalfEdgeWalker w = m.halfedgewalker(*hid);
// if(w.halfedge() < w.opp().halfedge() &&
// level_set_id[w.halfedge()] != 0 &&
// valency(m, w.vertex())+valency(m,w.opp().vertex())>20)
// hidvec.push_back(w.halfedge());
// }
//
// for(int i=0;i<hidvec.size(); ++i)
// {
// HalfEdgeWalker w = m.halfedgewalker(hidvec[i]);
// VertexID v = m.split_edge(hidvec[i]);
// level_set_id_vertex[v] = level_set_id[w.halfedge()];
// level_set_id[w.prev().halfedge()] = level_set_id[w.halfedge()];
// level_set_id[w.prev().opp().halfedge()] = level_set_id[w.halfedge()];
// shortest_edge_triangulate_face(m, w.face(), level_set_id_vertex);
// shortest_edge_triangulate_face(m, w.opp().face(), level_set_id_vertex);
// }
//
// priority_queue_optimization(m, tq_energy);
// }
// for(VertexIDIterator vid = m.vertices_begin(); vid != m.vertices_end(); ++vid)
// m.pos(*vid) = recalled_positions[*vid];
return;
priority_queue<EdgeQElem> edge_queue;
for(HalfEdgeIDIterator hid = m.halfedges_begin(); hid != m.halfedges_end(); ++hid)
{
HalfEdgeWalker w = m.halfedgewalker(*hid);
if(w.halfedge()<w.opp().halfedge() &&
level_set_id[w.halfedge()] == 0)
{
Vec3f v = normalize(m.pos(w.vertex()) - m.pos(w.opp().vertex()));
float weight = 0;
HalfEdgeWalker wo = m.halfedgewalker(w.opp().vertex());
for(; !w.full_circle(); w = w.circulate_vertex_ccw())
if(level_set_id[w.halfedge()] != 0)
{
Vec3f e = normalize(m.pos(w.vertex()) - m.pos(w.opp().vertex()));
weight += abs(dot(v,e));
}
for(; !wo.full_circle(); wo = wo.circulate_vertex_ccw())
if(level_set_id[wo.halfedge()] != 0)
{
Vec3f e = normalize(m.pos(wo.vertex()) - m.pos(wo.opp().vertex()));
weight += abs(dot(v,e));
}
edge_queue.push(EdgeQElem(weight, *hid));
}
}
while(!edge_queue.empty())
{
HalfEdgeID h = edge_queue.top().he;
edge_queue.pop();
HalfEdgeWalker w = m.halfedgewalker(h);
if(level_set_id[w.next().halfedge()] == 0 ||
level_set_id[w.prev().halfedge()] == 0 ||
level_set_id[w.opp().next().halfedge()] == 0 ||
level_set_id[w.opp().prev().halfedge()] == 0)
m.merge_faces(w.face(), w.halfedge());
}
for(FaceIDIterator fid = m.faces_begin(); fid != m.faces_end(); ++fid)
if(no_edges(m,*fid) >= 6)
shortest_edge_split_face(m, *fid, level_set_id_vertex);
for(HalfEdgeIDIterator hid = m.halfedges_begin(); hid != m.halfedges_end(); ++hid)
{
HalfEdgeWalker w = m.halfedgewalker(*hid);
if(level_set_id[w.halfedge()] != 0 &&
valency(m, w.vertex())==3 &&
valency(m, w.opp().vertex())==3 &&
((level_set_id[w.next().halfedge()] == 0 &&level_set_id[w.opp().next().halfedge()] == 0) ||
(level_set_id[w.prev().halfedge()] == 0 &&level_set_id[w.opp().prev().halfedge()] == 0)) &&
precond_collapse_edge(m, w.halfedge()))
{
cout << "collapsing!!!" << endl;
m.collapse_edge(w.halfedge(), true);
did_work = true;
}
}
// bool did_work = false;
return;
int k=0;
do {
++k;
did_work = false;
for(HalfEdgeIDIterator hid = m.halfedges_begin(); hid != m.halfedges_end(); ++hid)
{
HalfEdgeWalker w0 = m.halfedgewalker(*hid);
if(level_set_id[w0.halfedge()] != 0 &&
(level_set_id[w0.next().halfedge()] == 0 && level_set_id[w0.prev().halfedge()] == 0))
{
HalfEdgeWalker w = w0;
bool do_split = false;
for(;!w.full_circle(); w = w.next())
{
if(level_set_id[w.halfedge()] != 0 &&
(level_set_id[w.next().halfedge()] == level_set_id[w.halfedge()]))
do_split = true;
}
if(do_split)
{
VertexID v = m.split_edge(w0.halfedge());
level_set_id_vertex[v] = level_set_id[w0.halfedge()];
level_set_id[w0.prev().halfedge()] = level_set_id[w0.halfedge()];
level_set_id[w0.prev().opp().halfedge()] = level_set_id[w0.halfedge()];
did_work = true;
}
}
for(FaceIDIterator fid = m.faces_begin(); fid != m.faces_end(); ++fid)
if(no_edges(m,*fid) >= 6)
{
shortest_edge_split_face(m, *fid, level_set_id_vertex);
did_work = true;
}
}
} while (did_work && k<1);
}
void make_height_fun(const HMesh::Manifold& m, HMesh::VertexAttributeVector<float>& fun,
float& vmin, float& vmax)
{
VertexIDIterator vid = m.vertices_begin();
vmin = vmax = m.pos(*vid)[2];
for(; vid != m.vertices_end(); ++vid)
{
float v = m.pos(*vid)[1];
fun[*vid] = v;
vmin = min(v, vmin);
vmax = max(v, vmax);
}
}
// //-------------------------
// cout << "Remove short level set edges" << endl;
// float avglen=0;
// int n=0;
// for(HalfEdgeIDIterator hid = m.halfedges_begin(); hid != m.halfedges_end(); ++hid)
// {
// HalfEdgeWalker w = m.halfedgewalker(*hid);
// if(level_set_id[w.halfedge()] != 0)
// {
// avglen += length(m, w.halfedge());
// ++n;
// }
// }
// avglen /= n;
// for(HalfEdgeIDIterator hid = m.halfedges_begin(); hid != m.halfedges_end(); ++hid)
// if (length(m,*hid)<0.25*avglen && precond_collapse_edge(m, *hid)) {
// m.collapse_edge(*hid);
// }
// vector<HalfEdgeID> hidvec;
// for(HalfEdgeIDIterator hid = m.halfedges_begin(); hid != m.halfedges_end(); ++hid)
// {
// HalfEdgeWalker w = m.halfedgewalker(*hid);
// if(w.halfedge() < w.opp().halfedge() &&
// level_set_id[w.halfedge()] != 0 &&
// valency(m, w.vertex())+valency(m,w.opp().vertex())>12)
// hidvec.push_back(w.halfedge());
// }
//
// for(int i=0;i<hidvec.size(); ++i)
// {
// HalfEdgeWalker w = m.halfedgewalker(hidvec[i]);
// VertexID v = m.split_edge(hidvec[i]);
// level_set_id_vertex[v] = level_set_id[w.halfedge()];
// level_set_id[w.prev().halfedge()] = level_set_id[w.halfedge()];
// level_set_id[w.prev().opp().halfedge()] = level_set_id[w.halfedge()];
// shortest_edge_triangulate_face(m, w.face(), level_set_id_vertex);
// shortest_edge_triangulate_face(m, w.opp().face(), level_set_id_vertex);
// }
//
// // priority_queue_optimization(m, tq_energy);
//// ----------------------------
//cout << "smooth level set curves" << endl;
//
//for(int iter=0;iter<2;++iter)
//{
// VertexAttributeVector<Vec3f> new_pos(m.total_vertices(), Vec3f(0));
// for(HalfEdgeIDIterator hid = m.halfedges_begin(); hid != m.halfedges_end(); ++hid)
// {
// HalfEdgeWalker w = m.halfedgewalker(*hid);
// if(status[w.vertex()] == 1 && status[w.opp().vertex()] == 1 &&
// same_level(fun[w.vertex()], fun[w.opp().vertex()]))
// {
// new_pos[w.vertex()] += m.pos(w.vertex()) + m.pos(w.opp().vertex());
// }
// }
// for(VertexIDIterator vid = m.vertices_begin(); vid != m.vertices_end(); ++vid)
// if(status[*vid] == 1)
// m.pos(*vid) = new_pos[*vid] / (4.0);
//}
//return;
//
//for (int iter=0; iter<1; ++iter) {
// cout << __FILE__ << __LINE__ << endl;
// cout << __FILE__ << __LINE__ << endl;
// TriangleQuality tq_energy(level_set_id, level_set_id_vertex);
// priority_queue_optimization(m, tq_energy);
// // FunctionalDifference energy(fun, status);
// // priority_queue_optimization(m, energy);
//
// // ----------------------------
// 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);
//
// for (int i=0; i<vid_vec.size(); ++i) {
// FaceID f = m.merge_one_ring(vid_vec[i]);
// }
// cout << "-" << endl;
// }
// --------------------------
// Triangulate polygons by connecting only vertices on different curves.
// shortest_edge_triangulate(m);
// TriangleQuality tq_energy(level_set_id);
// priority_queue_optimization(m, tq_energy);
// k=0;
// do {
// ++k;
// did_work = false;
// priority_queue<EdgeQElem> edge_queue;
//
// for(HalfEdgeIDIterator hid = m.halfedges_begin(); hid != m.halfedges_end(); ++hid)
// {
// HalfEdgeWalker w = m.halfedgewalker(*hid);
// if(status[w.opp().vertex()] == 0 && status[w.vertex()] == 1)
// {
// float weight = (abs(fun[w.vertex()] - fun[w.opp().vertex()]*length(m, w.halfedge())) + 1e-6);
// edge_queue.push(EdgeQElem(weight, w.halfedge()));
// }
//
//
// while(!edge_queue.empty())
// {
// HalfEdgeID he = edge_queue.top().he;
// edge_queue.pop();
// if(m.in_use(he))
// {
// if(precond_collapse_edge(m,he))
// {
// m.collapse_edge(he);
// did_work = true;
// }
// }
// }
// }
// }while (did_work && k < 100);
//
// cout << "k=" << k << endl;
// priority_queue_optimization(m, energy);
//priority_queue<EdgeQElem> edge_queue;
//HalfEdgeAttributeVector<int> time_stamp(m.total_halfedges(), 0);
//for(HalfEdgeIDIterator hid = m.halfedges_begin(); hid != m.halfedges_end(); ++hid)
//if(level_set_id[*hid]==0)
//{
// HalfEdgeWalker w = m.halfedgewalker(*hid);
// if(w.halfedge()<w.opp().halfedge() && !(status[w.vertex()]==1 && status[w.opp().vertex()]==1))
// edge_queue.push(EdgeQElem(-(sqr(fun[w.vertex()]-fun[w.opp().vertex()]))*sqr_length(m.pos(w.vertex())-m.pos(w.opp().vertex())),*hid,0));
// }
//
//shortest_edge_triangulate(m);
//int k=0;
//while(!edge_queue.empty())
//{
// HalfEdgeID hid = edge_queue.top().he;
// if(m.in_use(hid) && time_stamp[hid]== edge_queue.top().time_stamp)
// {
// HalfEdgeWalker w = m.halfedgewalker(hid);
// HalfEdgeWalker wa = w.next();
// HalfEdgeWalker wb = w.prev().opp();
//
// if(m.merge_faces(w.face(), hid))
// {
// cout << ".";
// if(valency(m, wa.opp().vertex())==2 && precond_collapse_edge(m, wa.halfedge()))
// {
// HalfEdgeID h = wa.halfedge();
// wa = wa.prev();
// m.collapse_edge(h, false);
// ++time_stamp[wa.halfedge()];
// ++time_stamp[wa.opp().halfedge()];
// if(level_set_id[wa.halfedge()]==0 && !(status[wa.vertex()]==1 && status[wa.opp().vertex()]==1))
// edge_queue.push(EdgeQElem(-(sqr(fun[wa.vertex()]-fun[wa.opp().vertex()]))*sqr_length(m.pos(wa.vertex())-m.pos(wa.opp().vertex())),wa.halfedge(),time_stamp[wa.halfedge()]));
// }
// if(valency(m, wb.opp().vertex())==2 && precond_collapse_edge(m, wb.halfedge()))
// {
// HalfEdgeID h = wb.halfedge();
// wb = wb.prev();
// m.collapse_edge(h, false);
// ++time_stamp[wb.halfedge()];
// ++time_stamp[wb.opp().halfedge()];
// if(level_set_id[wb.halfedge()]==0 && !(status[wb.vertex()]==1 && status[wb.opp().vertex()]==1))
// edge_queue.push(EdgeQElem(-(sqr(fun[wb.vertex()]-fun[wb.opp().vertex()]))*sqr_length(m.pos(wb.vertex())-m.pos(wb.opp().vertex())),wb.halfedge(),time_stamp[wb.halfedge()]));
// }
//
// }
//
// }
//
// edge_queue.pop();
// }
//
// // bool did_work;
// do{
// did_work = false;
// for(HalfEdgeIDIterator hid = m.halfedges_begin(); hid != m.halfedges_end(); ++hid)
// {
// HalfEdgeWalker w = m.halfedgewalker(*hid);
// if(level_set_id[w.halfedge()] != 0 &&
// valency(m, w.vertex())==3 &&
// valency(m, w.opp().vertex())==3 &&
// ((level_set_id[w.next().halfedge()] == 0 &&level_set_id[w.opp().next().halfedge()] == 0) ||
// (level_set_id[w.prev().halfedge()] == 0 &&level_set_id[w.opp().prev().halfedge()] == 0)) &&
// precond_collapse_edge(m, w.halfedge()))
// {
// cout << "collapsing!!!" << endl;
// m.collapse_edge(w.halfedge(), true);
// did_work = true;
// }
// }
// }while (did_work);
// for(VertexIDIterator vid = m.vertices_begin(); vid != m.vertices_end(); ++vid)
// {
// bool is_max = true;
// bool is_min = true;
// HalfEdgeWalker w = m.halfedgewalker(*vid);
// float f = fun[*vid];
// for(;!w.full_circle(); w = w.circulate_vertex_ccw())
// {
// if(fun[w.vertex()] < f)
// is_min = false;
// if(fun[w.vertex()] > f)
// is_max = false;
//
// }
// if(is_max || is_min)
// status[*vid] = 2;
// }