Subversion Repositories gelsvn

//

//  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);

    }

}