WebSVN – gelsvn – Diff – /trunk/apps/MeshEdit/polarize.cpp


//
//  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>
 
#include <complex>
 
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(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())
    {
        Vec3d gdir = normalize(cross(n, m.pos(w.vertex()) - m.pos(w.opp().vertex())));
        float l = dot(m.pos(w.next().vertex())-m.pos(w.vertex()), gdir);
        gdir *= fun[w.next().vertex()]/l;
        gsum += gdir;
    }
    return gsum;
}
 
complex<double> complex_slerp(double t, const complex<double>& a, const complex<double>&b)
{
    double omega = (arg(a/b));
    cout << omega << endl;
    return (a*sin((1-t)*omega) + b*sin(t*omega))/sin(omega);
}
 
double mod_2pi(double x)
{
    return x-floor(x/(2.0 * M_PI))*2.0 * M_PI;
}
 
void extend_fun2(HMesh::Manifold& m, HMesh::HalfEdgeID h,
                 HMesh::VertexAttributeVector<double>& fun, 
                 HMesh::VertexAttributeVector<double>& 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);
    float v;
    if(dot_nba > 1e-5)
    {
        float t = dot(N, c-a)/dot_nba;
        double aval = fun2[w.opp().vertex()];
        double bval = fun2[w.vertex()];
        double val = (1-t)*aval + t*bval;
        v = mod_2pi(val);
        cout << t << " , " << v << endl; 
    }
    else
        v =  fun2[w.vertex()];
    fun2[w.next().vertex()] = v;
}
 
inline bool same_level(float a, float b) {return abs(a-b) < 0.00001;}
 
Vec3d laplacian(const Manifold& m, VertexID v)
{
    Vec3d avg_pos(0);
    float asum = 0;
    
    for(Walker w = m.walker(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<Vec3d> pos(m.allocated_vertices());
    
    for(VertexIDIterator v = m.vertices_begin(); v != m.vertices_end(); ++v){
        if(!boundary(m, *v))
        {
            Vec3d n = normal(m, *v);
            Vec3d 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<double>& fun, double vmin, double vmax, const int divisions, VertexAttributeVector<double>& 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);
    
    
    
    
    // ----------------------------
    cout << "Numbering the level sets" << endl;
    
    float max_length = 0;
    double max_length_fun = 0;
    int max_length_id =-1;
    HalfEdgeID max_length_h;
    HalfEdgeAttributeVector<int> level_set_id(m.allocated_halfedges(), 0);
    VertexAttributeVector<int> level_set_id_vertex(m.allocated_vertices(), 0);
    int no_id=1;
    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()] == 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();
            }            
            if(level_set_length > max_length)
            {
                max_length = level_set_length;
                max_length_id = no_id;
                max_length_h = w.halfedge();
                max_length_fun = fun[w.vertex()];
            }
            ++no_id;
        }
    }
    cout << max_length << " " << max_length_id << endl;
    cout << "Number of level sets : " << (no_id-1);
    // ----------------------------
    
    shortest_edge_triangulate(m);
    
    
    // -------------
    cout << "Parametrize level sets " << endl;
    
    queue<HalfEdgeID> hq;
    HalfEdgeAttributeVector<int> touched(m.no_halfedges(), 0);
    Walker w = m.walker(max_length_h);
    float param = 0;
    do
    {
        parametrization[w.opp().vertex()] = 2.0 * M_PI * param / max_length;
        param += length(m, w.halfedge());
        w = w.next();
        while(level_set_id[w.halfedge()] != max_length_id)
            w = w.circulate_vertex_cw();
        hq.push(w.halfedge());
        hq.push(w.opp().halfedge());
        touched[w.halfedge()] = 1;
        touched[w.opp().halfedge()] = 1;
    }            
    while(w.halfedge() != max_length_h);
    
    while(!hq.empty())
    {
        HalfEdgeID h = hq.front();
        hq.pop();
        if(!touched[w.next().halfedge()])
        {
            extend_fun2(m, h, fun, parametrization);
            touched[w.next().halfedge()] = 1;
            touched[w.prev().halfedge()] = 1;
            Walker w = m.walker(h);
            if(!touched[w.next().opp().halfedge()])
            {
                touched[w.next().opp().halfedge()] = 1;
                hq.push(w.next().opp().halfedge());
            }
            if(!touched[w.prev().opp().halfedge()])
            {
                touched[w.prev().opp().halfedge()] = 1;
                hq.push(w.prev().opp().halfedge());
                
            }
        }
    }
    
    return; 
    
    // ----------------------------
    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);
    
    
    
    //// ----------------------------
    cout << "smooth level set curves" << endl;
    
    
    for(int iter=0;iter< 14;++iter)
    {
        TriangleQualityValence tqv_energy(level_set_id_vertex);
        TriangleQuality tq_energy(level_set_id_vertex);
        priority_queue_optimization(m, tqv_energy);
        priority_queue_optimization(m, tq_energy);
        
        VertexAttributeVector<Vec3d> new_pos(m.allocated_vertices(), Vec3d(0));
        VertexAttributeVector<float> new_pos_wsum(m.allocated_vertices(), 0.0);
        for(HalfEdgeIDIterator hid = m.halfedges_begin(); hid != m.halfedges_end(); ++hid)
        {
            Walker w = m.walker(*hid);
            if(level_set_id[w.halfedge()] != 0)
            {
                float weight = 1.0;//sqr(valency(m,w.opp().vertex())-2);
                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.0;//sqr(valency(m,*vid)-2);
            new_pos[*vid] += weight*m.pos(*vid);
            new_pos_wsum[*vid] += weight;
            m.pos(*vid) = 0.75*m.pos(*vid) + 0.25 *new_pos[*vid] / new_pos_wsum[*vid];
        }
        
        priority_queue_optimization(m, tqv_energy);
        priority_queue_optimization(m, tq_energy);
        vector<HalfEdgeID> hidvec;
        for(HalfEdgeIDIterator hid = m.halfedges_begin(); hid != m.halfedges_end(); ++hid)
        {
            Walker w = m.walker(*hid);
            if(w.halfedge() < w.opp().halfedge() && 
               level_set_id[w.halfedge()] != 0 &&
               (
                ((level_set_id[w.opp().halfedge()] == level_set_id[w.opp().next().halfedge()] ||
                  level_set_id[w.halfedge()] == level_set_id[w.next().halfedge()]) &&
                 valency(m, w.vertex())+valency(m,w.opp().vertex())>10) ||
                
                valency(m, w.vertex())+valency(m,w.opp().vertex())>14
                )
               )
                hidvec.push_back(w.halfedge());
        }
        
        for(size_t i=0;i<hidvec.size(); ++i)
        {
            Walker w = m.walker(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, tqv_energy);
        priority_queue_optimization(m, tq_energy);
        
        for(HalfEdgeIDIterator hid = m.halfedges_begin(); hid != m.halfedges_end(); ++hid)
        {
            Walker w = m.walker(*hid);
            if(level_set_id[w.halfedge()] != 0 &&
               !(level_set_id[w.opp().halfedge()] == level_set_id[w.opp().next().halfedge()] ||
                 level_set_id[w.halfedge()] == level_set_id[w.next().halfedge()])  &&
               valency(m, w.vertex())<5 &&
               valency(m, w.opp().vertex())<5 &&
               precond_collapse_edge(m, w.halfedge()))
            {
                m.collapse_edge(w.halfedge(), true);
                did_work = true;
            }
        }
        
        
    }
    
    return;
    
    priority_queue<EdgeQElem> edge_queue;
    for(HalfEdgeIDIterator hid = m.halfedges_begin(); hid != m.halfedges_end(); ++hid)
    {
        Walker w = m.walker(*hid);
        if(w.halfedge()<w.opp().halfedge() &&
           level_set_id[w.halfedge()] == 0)
        {
            Vec3d v = normalize(m.pos(w.vertex()) - m.pos(w.opp().vertex()));
            
            float weight = 0;
            Walker wo = m.walker(w.opp().vertex());
            for(; !w.full_circle(); w = w.circulate_vertex_ccw())
                if(level_set_id[w.halfedge()] != 0)
                {
                    Vec3d 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)
                {
                    Vec3d 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();
        
        Walker w = m.walker(h);
        if(no_edges(m, w.face())==3 && no_edges(m,w.opp().face())==3 &&
           !(level_set_id[w.next().halfedge()] == level_set_id[w.opp().prev().halfedge()] ||
             level_set_id[w.prev().halfedge()] == level_set_id[w.opp().next().halfedge()]))
            m.merge_faces(w.face(), w.halfedge());
    }
    
    
    return;
    
    
    do{
        did_work = false;
        for(HalfEdgeIDIterator hid = m.halfedges_begin(); hid != m.halfedges_end(); ++hid)
        {
            Walker w = m.walker(*hid);
            
            if(level_set_id[w.halfedge()] != 0 && no_edges(m, w.face())==3 &&
               precond_collapse_edge(m, w.halfedge()))
            {
                m.collapse_edge(w.halfedge(), true);
                did_work = true;
            }
        }
    }
    while(did_work);
    
    
    return;
    
    
    int k=0;
    do {
        ++k;
        did_work = false;
        for(HalfEdgeIDIterator hid = m.halfedges_begin(); hid != m.halfedges_end(); ++hid)
        {
            Walker w0 = m.walker(*hid);
            
            if(level_set_id[w0.halfedge()] != 0 &&
               (level_set_id[w0.next().halfedge()] == 0 && level_set_id[w0.prev().halfedge()] == 0))
            {
                
                Walker 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<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 = 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)
//        {
//            Walker w = m.walker(*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)
//        {
//            Walker w = m.walker(*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)
//        {
//            Walker w = m.walker(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<Vec3d> new_pos(m.allocated_vertices(), Vec3d(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()]))
//        {
//            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)
//            {
//                Walker w = m.walker(*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.allocated_halfedges(), 0);
//for(HalfEdgeIDIterator hid = m.halfedges_begin(); hid != m.halfedges_end(); ++hid)
//if(level_set_id[*hid]==0)
//{
//    Walker w = m.walker(*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)
//    {
//        Walker w = m.walker(hid);
//        Walker wa = w.next();
//        Walker 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)
//    {
//        Walker w = m.walker(*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;
//        Walker w = m.walker(*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;
//    }
 
 
 
 

Rev 596	Rev 607
1	`//`	1	`//`
2	`// polarize.cpp`	2	`// polarize.cpp`
3	`// GEL`	3	`// GEL`
4	`//`	4	`//`
5	`// Created by J. Andreas Bærentzen on 18/03/12.`	5	`// Created by J. Andreas Bærentzen on 18/03/12.`
6	`// Copyright 2012 __MyCompanyName__. All rights reserved.`	6	`// Copyright 2012 __MyCompanyName__. All rights reserved.`
7	`//`	7	`//`
8	`#include <queue>`	8	`#include <queue>`
9		9
10	`#include "polarize.h"`	10	`#include "polarize.h"`
11	`#include <HMesh/triangulate.h>`	11	`#include <HMesh/triangulate.h>`
12	`#include <HMesh/curvature.h>`	12	`#include <HMesh/curvature.h>`
13	`#include <HMesh/quadric_simplify.h>`	13	`#include <HMesh/quadric_simplify.h>`
14	`#include <HMesh/mesh_optimization.h>`	14	`#include <HMesh/mesh_optimization.h>`
15	`#include <HMesh/smooth.h>`	15	`#include <HMesh/smooth.h>`
16		16
17	`#include <complex>`	17	`#include <complex>`
18		18
19	`using namespace CGLA;`	19	`using namespace CGLA;`
20	`using namespace std;`	20	`using namespace std;`
21	`using namespace HMesh;`	21	`using namespace HMesh;`
22		22
23	`void shortest_edge_triangulate_face(Manifold& m, FaceID f0, VertexAttributeVector<int>& level_set_id_vertex)`	23	`void shortest_edge_triangulate_face(Manifold& m, FaceID f0, VertexAttributeVector<int>& level_set_id_vertex)`
24	`{`	24	`{`
25	`queue<FaceID> face_queue;`	25	`queue<FaceID> face_queue;`
26		26
27	`face_queue.push(f0);`	27	`face_queue.push(f0);`
28		28
29	`while(!face_queue.empty())`	29	`while(!face_queue.empty())`
30	`{`	30	`{`
31	`FaceID f = face_queue.front();`	31	`FaceID f = face_queue.front();`
32	`face_queue.pop();`	32	`face_queue.pop();`
33		33
34	`// Create a vector of vertices.`	34	`// Create a vector of vertices.`
35	`vector<VertexID> verts;`	35	`vector<VertexID> verts;`
36	`for(Walker w = m.walker(f); !w.full_circle(); w = w.circulate_face_ccw())`	36	`for(Walker w = m.walker(f); !w.full_circle(); w = w.circulate_face_ccw())`
37	`{`	37	`{`
38	`FaceID fa = w.face();`	38	`FaceID fa = w.face();`
39	`FaceID fb = f;`	39	`FaceID fb = f;`
40	`assert(fa==fb);`	40	`assert(fa==fb);`
41	`verts.push_back(w.vertex());`	41	`verts.push_back(w.vertex());`
42	`}`	42	`}`
43	`// If there are just three we are done.`	43	`// If there are just three we are done.`
44	`if(verts.size() == 3) continue;`	44	`if(verts.size() == 3) continue;`
45		45
46	`// Find vertex pairs that may be connected.`	46	`// Find vertex pairs that may be connected.`
47	`vector<pair<int,int> > vpairs;`	47	`vector<pair<int,int> > vpairs;`
48	`const int N = verts.size();`	48	`const int N = verts.size();`
49	`for(int i = 0; i < N - 2; ++i){`	49	`for(int i = 0; i < N - 2; ++i){`
50	`for(int j = i + 2; j < N; ++j){`	50	`for(int j = i + 2; j < N; ++j){`
51	`if(verts[i] != verts[j] &&`	51	`if(verts[i] != verts[j] &&`
52	`!connected(m, verts[i], verts[j]) &&`	52	`!connected(m, verts[i], verts[j]) &&`
53	`(level_set_id_vertex[verts[i]] == 0 \|\| level_set_id_vertex[verts[i]] != level_set_id_vertex[verts[j]])`	53	`(level_set_id_vertex[verts[i]] == 0 \|\| level_set_id_vertex[verts[i]] != level_set_id_vertex[verts[j]])`
54	`)`	54	`)`
55	`vpairs.push_back(pair<int,int>(i, j));`	55	`vpairs.push_back(pair<int,int>(i, j));`
56	`}`	56	`}`
57	`}`	57	`}`
58	`if(vpairs.empty()){`	58	`if(vpairs.empty()){`
59	`cout << "Warning: could not triangulate a face."`	59	`cout << "Warning: could not triangulate a face."`
60	`<< "Probably a vertex appears more than one time in other vertex's one-ring" << endl;`	60	`<< "Probably a vertex appears more than one time in other vertex's one-ring" << endl;`
61	`continue;`	61	`continue;`
62	`}`	62	`}`
63		63
64	`/* For all vertex pairs, find the edge lengths. Combine the`	64	`/* For all vertex pairs, find the edge lengths. Combine the`
65	`vertices forming the shortest edge. */`	65	`vertices forming the shortest edge. */`
66		66
67	`float min_len=FLT_MAX;`	67	`float min_len=FLT_MAX;`
68	`int min_k = -1;`	68	`int min_k = -1;`
69	`for(size_t k = 0; k < vpairs.size(); ++k){`	69	`for(size_t k = 0; k < vpairs.size(); ++k){`
70	`int i = vpairs[k].first;`	70	`int i = vpairs[k].first;`
71	`int j = vpairs[k].second;`	71	`int j = vpairs[k].second;`
72	`float len = sqr_length(m.pos(verts[i]) - m.pos(verts[j]));`	72	`float len = sqr_length(m.pos(verts[i]) - m.pos(verts[j]));`
73		73
74	`if(len<min_len){`	74	`if(len<min_len){`
75	`min_len = len;`	75	`min_len = len;`
76	`min_k = k;`	76	`min_k = k;`
77	`}`	77	`}`
78	`}`	78	`}`
79	`assert(min_k != -1);`	79	`assert(min_k != -1);`
80		80
81	`{`	81	`{`
82	`// Split faces along edge whose midpoint is closest to isovalue`	82	`// Split faces along edge whose midpoint is closest to isovalue`
83	`int i = vpairs[min_k].first;`	83	`int i = vpairs[min_k].first;`
84	`int j = vpairs[min_k].second;`	84	`int j = vpairs[min_k].second;`
85	`FaceID f_new = m.split_face_by_edge(f, verts[i], verts[j]);`	85	`FaceID f_new = m.split_face_by_edge(f, verts[i], verts[j]);`
86		86
87	`if(no_edges(m, f)>3)`	87	`if(no_edges(m, f)>3)`
88	`face_queue.push(f);`	88	`face_queue.push(f);`
89	`if(no_edges(m, f_new)>3)`	89	`if(no_edges(m, f_new)>3)`
90	`face_queue.push(f_new);`	90	`face_queue.push(f_new);`
91	`}`	91	`}`
92		92
93	`}`	93	`}`
94	`}`	94	`}`
95		95
96		96
97	`void shortest_edge_split_face(Manifold& m, FaceID f0, VertexAttributeVector<int>& level_set_id_vertex)`	97	`void shortest_edge_split_face(Manifold& m, FaceID f0, VertexAttributeVector<int>& level_set_id_vertex)`
98	`{`	98	`{`
99	`queue<FaceID> face_queue;`	99	`queue<FaceID> face_queue;`
100		100
101	`face_queue.push(f0);`	101	`face_queue.push(f0);`
102		102
103	`while(!face_queue.empty())`	103	`while(!face_queue.empty())`
104	`{`	104	`{`
105	`FaceID f = face_queue.front();`	105	`FaceID f = face_queue.front();`
106	`face_queue.pop();`	106	`face_queue.pop();`
107		107
108	`// Create a vector of vertices.`	108	`// Create a vector of vertices.`
109	`vector<VertexID> verts;`	109	`vector<VertexID> verts;`
110	`for(Walker w = m.walker(f); !w.full_circle(); w = w.circulate_face_ccw())`	110	`for(Walker w = m.walker(f); !w.full_circle(); w = w.circulate_face_ccw())`
111	`{`	111	`{`
112	`verts.push_back(w.vertex());`	112	`verts.push_back(w.vertex());`
113	`}`	113	`}`
114		114
115	`// Find vertex pairs that may be connected.`	115	`// Find vertex pairs that may be connected.`
116	`vector<pair<int,int> > vpairs;`	116	`vector<pair<int,int> > vpairs;`
117	`const int N = verts.size();`	117	`const int N = verts.size();`
118	`for(int i = 0; i < N ; ++i){`	118	`for(int i = 0; i < N ; ++i){`
119	`for(int j = 3; j < N-2; ++j){`	119	`for(int j = 3; j < N-2; ++j){`
120	`int jj = (j+i)%N;`	120	`int jj = (j+i)%N;`
121	`if(!connected(m, verts[i], verts[jj]) &&`	121	`if(!connected(m, verts[i], verts[jj]) &&`
122	`(level_set_id_vertex[verts[i]] != level_set_id_vertex[verts[jj]]) &&`	122	`(level_set_id_vertex[verts[i]] != level_set_id_vertex[verts[jj]]) &&`
123	`(level_set_id_vertex[verts[(i+1)%N]] == level_set_id_vertex[verts[i]]) &&`	123	`(level_set_id_vertex[verts[(i+1)%N]] == level_set_id_vertex[verts[i]]) &&`
124	`(level_set_id_vertex[verts[(i+N-1)%N]] == level_set_id_vertex[verts[i]]) &&`	124	`(level_set_id_vertex[verts[(i+N-1)%N]] == level_set_id_vertex[verts[i]]) &&`
125	`(level_set_id_vertex[verts[(jj+1)%N]] == level_set_id_vertex[verts[jj]]) &&`	125	`(level_set_id_vertex[verts[(jj+1)%N]] == level_set_id_vertex[verts[jj]]) &&`
126	`(level_set_id_vertex[verts[(jj+N-1)%N]] == level_set_id_vertex[verts[jj]]))`	126	`(level_set_id_vertex[verts[(jj+N-1)%N]] == level_set_id_vertex[verts[jj]]))`
127	`vpairs.push_back(pair<int,int>(i, jj));`	127	`vpairs.push_back(pair<int,int>(i, jj));`
128	`}`	128	`}`
129	`}`	129	`}`
130	`if(vpairs.empty()) continue;`	130	`if(vpairs.empty()) continue;`
131		131
132		132
133	`/* For all vertex pairs, find the edge lengths. Combine the`	133	`/* For all vertex pairs, find the edge lengths. Combine the`
134	`vertices forming the shortest edge. */`	134	`vertices forming the shortest edge. */`
135		135
136	`float min_len=FLT_MAX;`	136	`float min_len=FLT_MAX;`
137	`int min_k = -1;`	137	`int min_k = -1;`
138	`for(size_t k = 0; k < vpairs.size(); ++k){`	138	`for(size_t k = 0; k < vpairs.size(); ++k){`
139	`int i = vpairs[k].first;`	139	`int i = vpairs[k].first;`
140	`int j = vpairs[k].second;`	140	`int j = vpairs[k].second;`
141	`float len = sqr_length(m.pos(verts[i]) - m.pos(verts[j]));`	141	`float len = sqr_length(m.pos(verts[i]) - m.pos(verts[j]));`
142		142
143	`if(len<min_len){`	143	`if(len<min_len){`
144	`min_len = len;`	144	`min_len = len;`
145	`min_k = k;`	145	`min_k = k;`
146	`}`	146	`}`
147	`}`	147	`}`
148	`assert(min_k != -1);`	148	`assert(min_k != -1);`
149		149
150	`{`	150	`{`
151	`// Split faces along edge whose midpoint is closest to isovalue`	151	`// Split faces along edge whose midpoint is closest to isovalue`
152	`int i = vpairs[min_k].first;`	152	`int i = vpairs[min_k].first;`
153	`int j = vpairs[min_k].second;`	153	`int j = vpairs[min_k].second;`
154	`FaceID f_new = m.split_face_by_edge(f, verts[i], verts[j]);`	154	`FaceID f_new = m.split_face_by_edge(f, verts[i], verts[j]);`
155	`if(no_edges(m, f)>5)`	155	`if(no_edges(m, f)>5)`
156	`face_queue.push(f);`	156	`face_queue.push(f);`
157	`if(no_edges(m, f_new)>5)`	157	`if(no_edges(m, f_new)>5)`
158	`face_queue.push(f_new);`	158	`face_queue.push(f_new);`
159	`}`	159	`}`
160		160
161	`}`	161	`}`
162	`}`	162	`}`
163		163
164		164
165		165
166	`struct EdgeQElem {`	166	`struct EdgeQElem {`
167	`float priority;`	167	`float priority;`
168	`HalfEdgeID he;`	168	`HalfEdgeID he;`
169	`EdgeQElem(float _priority, HalfEdgeID _he): priority(_priority), he(_he) {}`	169	`EdgeQElem(float _priority, HalfEdgeID _he): priority(_priority), he(_he) {}`
170	`};`	170	`};`
171		171
172	`bool operator<(const EdgeQElem& l, const EdgeQElem& r)`	172	`bool operator<(const EdgeQElem& l, const EdgeQElem& r)`
173	`{`	173	`{`
174	`return l.priority < r.priority;`	174	`return l.priority < r.priority;`
175	`}`	175	`}`
176		176
177	`class FunctionalDifference: public EnergyFun`	177	`class FunctionalDifference: public EnergyFun`
178	`{`	178	`{`
179	`VertexAttributeVector<float>& fun;`	179	`VertexAttributeVector<float>& fun;`
180	`VertexAttributeVector<int>& status;`	180	`VertexAttributeVector<int>& status;`
181	`public:`	181	`public:`
182	`FunctionalDifference(VertexAttributeVector<float>& _fun, VertexAttributeVector<int>& _status): fun(_fun), status(_status) {}`	182	`FunctionalDifference(VertexAttributeVector<float>& _fun, VertexAttributeVector<int>& _status): fun(_fun), status(_status) {}`
183	`virtual double delta_energy(const Manifold& m, HalfEdgeID h) const`	183	`virtual double delta_energy(const Manifold& m, HalfEdgeID h) const`
184	`{`	184	`{`
185	`Walker w = m.walker(h);`	185	`Walker w = m.walker(h);`
186	`if(status[w.vertex()] == 1 && status[w.opp().vertex()]==1)`	186	`if(status[w.vertex()] == 1 && status[w.opp().vertex()]==1)`
187	`return DBL_MAX;`	187	`return DBL_MAX;`
188	`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()]));`	188	`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()]));`
189	`}`	189	`}`
190	`};`	190	`};`
191		191
192	`class TriangleQualityValence: public EnergyFun`	192	`class TriangleQualityValence: public EnergyFun`
193	`{`	193	`{`
194	`VertexAttributeVector<int>& idv;`	194	`VertexAttributeVector<int>& idv;`
195	`MinAngleEnergy mae;`	195	`MinAngleEnergy mae;`
196	`ValencyEnergy vae;`	196	`ValencyEnergy vae;`
197	`public:`	197	`public:`
198	`TriangleQualityValence(VertexAttributeVector<int>& _idv): idv(_idv), mae(-1) {}`	198	`TriangleQualityValence(VertexAttributeVector<int>& _idv): idv(_idv), mae(-1) {}`
199	`virtual double delta_energy(const Manifold& m, HalfEdgeID h) const`	199	`virtual double delta_energy(const Manifold& m, HalfEdgeID h) const`
200	`{`	200	`{`
201	`Walker w = m.walker(h);`	201	`Walker w = m.walker(h);`
202	`if(idv[w.next().vertex()] == idv[w.opp().next().vertex()] \|\|`	202	`if(idv[w.next().vertex()] == idv[w.opp().next().vertex()] \|\|`
203	`idv[w.vertex()] == idv[w.opp().vertex()])`	203	`idv[w.vertex()] == idv[w.opp().vertex()])`
204	`return DBL_MAX;`	204	`return DBL_MAX;`
205		205
206	`VertexID v1 = w.opp().vertex();`	206	`VertexID v1 = w.opp().vertex();`
207	`VertexID v2 = w.vertex();`	207	`VertexID v2 = w.vertex();`
208	`VertexID vo1 = w.next().vertex();`	208	`VertexID vo1 = w.next().vertex();`
209	`VertexID vo2 = w.opp().next().vertex();`	209	`VertexID vo2 = w.opp().next().vertex();`
210		210
211	`int val1 = valency(m, v1);`	211	`int val1 = valency(m, v1);`
212	`int val2 = valency(m, v2);`	212	`int val2 = valency(m, v2);`
213	`int valo1 = valency(m, vo1);`	213	`int valo1 = valency(m, vo1);`
214	`int valo2 = valency(m, vo2);`	214	`int valo2 = valency(m, vo2);`
215		215
216	`// The optimal valency is four for a boundary vertex`	216	`// The optimal valency is four for a boundary vertex`
217	`// and six elsewhere.`	217	`// and six elsewhere.`
218	`int t1 = boundary(m, v1) ? 4 : 6;`	218	`int t1 = boundary(m, v1) ? 4 : 6;`
219	`int t2 = boundary(m, v2) ? 4 : 6;`	219	`int t2 = boundary(m, v2) ? 4 : 6;`
220	`int to1 = boundary(m, vo1) ? 4 : 6;`	220	`int to1 = boundary(m, vo1) ? 4 : 6;`
221	`int to2 = boundary(m, vo2) ? 4 : 6;`	221	`int to2 = boundary(m, vo2) ? 4 : 6;`
222		222
223	`int before =`	223	`int before =`
224	`max(max(sqr(val1-t1),sqr(val2-t2)), max(sqr(valo1-to1), sqr(valo2-to2)));`	224	`max(max(sqr(val1-t1),sqr(val2-t2)), max(sqr(valo1-to1), sqr(valo2-to2)));`
225	`int after =`	225	`int after =`
226	`max(max(sqr(valo1+1-to1),sqr(val1-1-t1)), max(sqr(val2-1-t2),sqr(valo2+1-to2)));`	226	`max(max(sqr(valo1+1-to1),sqr(val1-1-t1)), max(sqr(val2-1-t2),sqr(valo2+1-to2)));`
227		227
228	`return static_cast<double>(after-before);`	228	`return static_cast<double>(after-before);`
229		229
230	`// return vae.delta_energy(m,h);`	230	`// return vae.delta_energy(m,h);`
231	`float la = length(m.pos(w.next().vertex())-m.pos(w.opp().next().vertex()));`	231	`float la = length(m.pos(w.next().vertex())-m.pos(w.opp().next().vertex()));`
232	`float lb = length(m.pos(w.vertex())-m.pos(w.opp().vertex()));`	232	`float lb = length(m.pos(w.vertex())-m.pos(w.opp().vertex()));`
233	`return la-lb;`	233	`return la-lb;`
234	`return mae.delta_energy(m,h);`	234	`return mae.delta_energy(m,h);`
235	`}`	235	`}`
236	`};`	236	`};`
237		237
238	`class TriangleQuality: public EnergyFun`	238	`class TriangleQuality: public EnergyFun`
239	`{`	239	`{`
240	`VertexAttributeVector<int>& idv;`	240	`VertexAttributeVector<int>& idv;`
241	`MinAngleEnergy mae;`	241	`MinAngleEnergy mae;`
242	`ValencyEnergy vae;`	242	`ValencyEnergy vae;`
243	`public:`	243	`public:`
244	`TriangleQuality(VertexAttributeVector<int>& _idv): idv(_idv), mae(-1) {}`	244	`TriangleQuality(VertexAttributeVector<int>& _idv): idv(_idv), mae(-1) {}`
245	`virtual double delta_energy(const Manifold& m, HalfEdgeID h) const`	245	`virtual double delta_energy(const Manifold& m, HalfEdgeID h) const`
246	`{`	246	`{`
247	`Walker w = m.walker(h);`	247	`Walker w = m.walker(h);`
248	`if(idv[w.next().vertex()] == idv[w.opp().next().vertex()] \|\|`	248	`if(idv[w.next().vertex()] == idv[w.opp().next().vertex()] \|\|`
249	`idv[w.vertex()] == idv[w.opp().vertex()])`	249	`idv[w.vertex()] == idv[w.opp().vertex()])`
250	`return DBL_MAX;`	250	`return DBL_MAX;`
251	`return mae.delta_energy(m,h);`	251	`return mae.delta_energy(m,h);`
252	`}`	252	`}`
253	`};`	253	`};`
254		254
255	`Vec3d grad(HMesh::Manifold& m, HMesh::VertexAttributeVector<double>& fun, HMesh::FaceID f)`	255	`Vec3d grad(HMesh::Manifold& m, HMesh::VertexAttributeVector<double>& fun, HMesh::FaceID f)`
256	`{`	256	`{`
257	`if(no_edges(m,f) != 3)`	257	`if(no_edges(m,f) != 3)`
258	`return Vec3d(0.0);`	258	`return Vec3d(0.0);`
259		259
260		260
261	`Vec3d n = normal(m, f);`	261	`Vec3d n = normal(m, f);`
262		262
263	`Vec3d gsum(0.0);`	263	`Vec3d gsum(0.0);`
264	`for(Walker w = m.walker(f); !w.full_circle(); w = w.next())`	264	`for(Walker w = m.walker(f); !w.full_circle(); w = w.next())`
265	`{`	265	`{`
266	`Vec3d gdir = normalize(cross(n, m.pos(w.vertex()) - m.pos(w.opp().vertex())));`	266	`Vec3d gdir = normalize(cross(n, m.pos(w.vertex()) - m.pos(w.opp().vertex())));`
267	`float l = dot(m.pos(w.next().vertex())-m.pos(w.vertex()), gdir);`	267	`float l = dot(m.pos(w.next().vertex())-m.pos(w.vertex()), gdir);`
268	`gdir *= fun[w.next().vertex()]/l;`	268	`gdir *= fun[w.next().vertex()]/l;`
269	`gsum += gdir;`	269	`gsum += gdir;`
270	`}`	270	`}`
271	`return gsum;`	271	`return gsum;`
272	`}`	272	`}`
273		273
274	`complex<double> complex_slerp(double t, const complex<double>& a, const complex<double>&b)`	274	`complex<double> complex_slerp(double t, const complex<double>& a, const complex<double>&b)`
275	`{`	275	`{`
276	`double omega = (arg(a/b));`	276	`double omega = (arg(a/b));`
277	`cout << omega << endl;`	277	`cout << omega << endl;`
278	`return (asin((1-t)omega) + bsin(tomega))/sin(omega);`	278	`return (asin((1-t)omega) + bsin(tomega))/sin(omega);`
279	`}`	279	`}`
280		280
281	`double mod_2pi(double x)`	281	`double mod_2pi(double x)`
282	`{`	282	`{`
283	`return x-floor(x/(2.0 * M_PI))2.0 M_PI;`	283	`return x-floor(x/(2.0 * M_PI))2.0 M_PI;`
284	`}`	284	`}`
285		285
286	`void extend_fun2(HMesh::Manifold& m, HMesh::HalfEdgeID h,`	286	`void extend_fun2(HMesh::Manifold& m, HMesh::HalfEdgeID h,`
287	`HMesh::VertexAttributeVector<double>& fun,`	287	`HMesh::VertexAttributeVector<double>& fun,`
288	`HMesh::VertexAttributeVector<double>& fun2)`	288	`HMesh::VertexAttributeVector<double>& fun2)`
289	`{`	289	`{`
290	`Walker w = m.walker(h);`	290	`Walker w = m.walker(h);`
291	`FaceID f = w.face();`	291	`FaceID f = w.face();`
292	`Vec3d a = m.pos(w.opp().vertex());`	292	`Vec3d a = m.pos(w.opp().vertex());`
293	`Vec3d b = m.pos(w.vertex());`	293	`Vec3d b = m.pos(w.vertex());`
294	`Vec3d c = m.pos(w.next().vertex());`	294	`Vec3d c = m.pos(w.next().vertex());`
295	`Vec3d g = grad(m, fun, f);`	295	`Vec3d g = grad(m, fun, f);`
296	`Vec3d n = normal(m, f);`	296	`Vec3d n = normal(m, f);`
297	`Vec3d N = -normalize(cross(g, n));`	297	`Vec3d N = -normalize(cross(g, n));`
298	`float dot_nba = dot(N,b-a);`	298	`float dot_nba = dot(N,b-a);`
299	`float v;`	299	`float v;`
300	`if(dot_nba > 1e-5)`	300	`if(dot_nba > 1e-5)`
301	`{`	301	`{`
302	`float t = dot(N, c-a)/dot_nba;`	302	`float t = dot(N, c-a)/dot_nba;`
303	`double aval = fun2[w.opp().vertex()];`	303	`double aval = fun2[w.opp().vertex()];`
304	`double bval = fun2[w.vertex()];`	304	`double bval = fun2[w.vertex()];`
305	`double val = (1-t)aval + tbval;`	305	`double val = (1-t)aval + tbval;`
306	`v = mod_2pi(val);`	306	`v = mod_2pi(val);`
307	`cout << t << " , " << v << endl;`	307	`cout << t << " , " << v << endl;`
308	`}`	308	`}`
309	`else`	309	`else`
310	`v = fun2[w.vertex()];`	310	`v = fun2[w.vertex()];`
311	`fun2[w.next().vertex()] = v;`	311	`fun2[w.next().vertex()] = v;`
312	`}`	312	`}`
313		313
314	`inline bool same_level(float a, float b) {return abs(a-b) < 0.00001;}`	314	`inline bool same_level(float a, float b) {return abs(a-b) < 0.00001;}`
315		315
316	`Vec3d laplacian(const Manifold& m, VertexID v)`	316	`Vec3d laplacian(const Manifold& m, VertexID v)`
317	`{`	317	`{`
318	`Vec3d avg_pos(0);`	318	`Vec3d avg_pos(0);`
319	`float asum = 0;`	319	`float asum = 0;`
320		320
321	`for(Walker w = m.walker(v); !w.full_circle(); w = w.circulate_vertex_cw()){`	321	`for(Walker w = m.walker(v); !w.full_circle(); w = w.circulate_vertex_cw()){`
322	`float a = barycentric_area(m, w.vertex());`	322	`float a = barycentric_area(m, w.vertex());`
323	`avg_pos += a * m.pos(w.vertex());`	323	`avg_pos += a * m.pos(w.vertex());`
324	`asum += a;`	324	`asum += a;`
325	`}`	325	`}`
326	`return avg_pos / asum - m.pos(v);`	326	`return avg_pos / asum - m.pos(v);`
327	`}`	327	`}`
328		328
329	`void aw_laplacian_smooth(Manifold& m, float t)`	329	`void aw_laplacian_smooth(Manifold& m, float t)`
330	`{`	330	`{`
331	`VertexAttributeVector<Vec3d> pos(m.allocated_vertices());`	331	`VertexAttributeVector<Vec3d> pos(m.allocated_vertices());`
332		332
333	`for(VertexIDIterator v = m.vertices_begin(); v != m.vertices_end(); ++v){`	333	`for(VertexIDIterator v = m.vertices_begin(); v != m.vertices_end(); ++v){`
334	`if(!boundary(m, *v))`	334	`if(!boundary(m, *v))`
335	`{`	335	`{`
336	`Vec3d n = normal(m, *v);`	336	`Vec3d n = normal(m, *v);`
337	`Vec3d l = laplacian(m, *v);`	337	`Vec3d l = laplacian(m, *v);`
338	`if(sqr_length(n) > 0.8)`	338	`if(sqr_length(n) > 0.8)`
339	`l -= n * dot(n,l);`	339	`l -= n * dot(n,l);`
340	`pos[v] = t l + m.pos(*v);`	340	`pos[v] = t l + m.pos(*v);`
341	`}`	341	`}`
342	`}`	342	`}`
343		343
344	`for(VertexIDIterator v = m.vertices_begin(); v != m.vertices_end(); ++v){`	344	`for(VertexIDIterator v = m.vertices_begin(); v != m.vertices_end(); ++v){`
345	`if(!boundary(m, *v))`	345	`if(!boundary(m, *v))`
346	`m.pos(v) = pos[v];`	346	`m.pos(v) = pos[v];`
347	`}`	347	`}`
348	`}`	348	`}`
349		349
350	`void polarize_mesh(Manifold& m, VertexAttributeVector<double>& fun, double vmin, double vmax, const int divisions, VertexAttributeVector<double>& parametrization)`	350	`void polarize_mesh(Manifold& m, VertexAttributeVector<double>& fun, double vmin, double vmax, const int divisions, VertexAttributeVector<double>& parametrization)`
351	`{`	351	`{`
352	`vmax -= 0.01 * (vmax-vmin);`	352	`vmax -= 0.01 * (vmax-vmin);`
353	`float interval = (vmax-vmin)/divisions;`	353	`float interval = (vmax-vmin)/divisions;`
354		354
355	`VertexAttributeVector<int> status(m.allocated_vertices(), 0);`	355	`VertexAttributeVector<int> status(m.allocated_vertices(), 0);`
356		356
357		357
358	`// ----------------------------------------`	358	`// ----------------------------------------`
359	`cout << "Tracing level set curves" << endl;`	359	`cout << "Tracing level set curves" << endl;`
360		360
361	`vector<HalfEdgeID> hidvec;`	361	`vector<HalfEdgeID> hidvec;`
362	`for(HalfEdgeIDIterator hid = m.halfedges_begin(); hid != m.halfedges_end(); ++hid)`	362	`for(HalfEdgeIDIterator hid = m.halfedges_begin(); hid != m.halfedges_end(); ++hid)`
363	`{`	363	`{`
364	`Walker w = m.walker(*hid);`	364	`Walker w = m.walker(*hid);`
365	`if(fun[w.vertex()] > fun[w.opp().vertex()])`	365	`if(fun[w.vertex()] > fun[w.opp().vertex()])`
366	`hidvec.push_back(*hid);`	366	`hidvec.push_back(*hid);`
367	`}`	367	`}`
368		368
369	`for(int i = 0; i<hidvec.size(); ++i)`	369	`for(size_t i = 0; i<hidvec.size(); ++i)`
370	`{`	370	`{`
371	`Walker w = m.walker(hidvec[i]);`	371	`Walker w = m.walker(hidvec[i]);`
372		372
373	`float b = (fun[w.vertex()]- vmin)/interval;`	373	`float b = (fun[w.vertex()]- vmin)/interval;`
374	`float a = (fun[w.opp().vertex()] - vmin)/interval;`	374	`float a = (fun[w.opp().vertex()] - vmin)/interval;`
375	`float floor_b = floor(b);`	375	`float floor_b = floor(b);`
376	`float floor_a = floor(a);`	376	`float floor_a = floor(a);`
377		377
378	`Vec3d pb = m.pos(w.vertex());`	378	`Vec3d pb = m.pos(w.vertex());`
379	`for(int j=floor_b; j>floor_a; --j)`	379	`for(int j=floor_b; j>floor_a; --j)`
380	`{`	380	`{`
381	`float t = (j-a) / (b-a);`	381	`float t = (j-a) / (b-a);`
382	`Vec3d p = t * pb + (1.0-t) * m.pos(w.opp().vertex());`	382	`Vec3d p = t * pb + (1.0-t) * m.pos(w.opp().vertex());`
383	`VertexID v_new = m.split_edge(w.halfedge());`	383	`VertexID v_new = m.split_edge(w.halfedge());`
384	`w = w.prev();`	384	`w = w.prev();`
385	`status[v_new] = 1;`	385	`status[v_new] = 1;`
386	`fun[v_new] = j * interval + vmin;`	386	`fun[v_new] = j * interval + vmin;`
387	`m.pos(v_new) = p;`	387	`m.pos(v_new) = p;`
388	`}`	388	`}`
389	`}`	389	`}`
390		390
391	`bool did_work;`	391	`bool did_work;`
392	`do`	392	`do`
393	`{`	393	`{`
394	`did_work = false;`	394	`did_work = false;`
395		395
396	`for(FaceIDIterator fid = m.faces_begin(); fid != m.faces_end(); ++fid)`	396	`for(FaceIDIterator fid = m.faces_begin(); fid != m.faces_end(); ++fid)`
397	`for(Walker w = m.walker(*fid);!w.full_circle(); w = w.next())`	397	`for(Walker w = m.walker(*fid);!w.full_circle(); w = w.next())`
398	`if(status[w.vertex()] == 1 && !(status[w.next().vertex()]==1 && same_level(fun[w.vertex()],fun[w.next().vertex()]))`	398	`if(status[w.vertex()] == 1 && !(status[w.next().vertex()]==1 && same_level(fun[w.vertex()],fun[w.next().vertex()]))`
399	`&& !(status[w.prev().vertex()]==1 && same_level(fun[w.vertex()],fun[w.prev().vertex()])))`	399	`&& !(status[w.prev().vertex()]==1 && same_level(fun[w.vertex()],fun[w.prev().vertex()])))`
400	`{`	400	`{`
401	`Walker w0 = w;`	401	`Walker w0 = w;`
402	`w = w.next().next();`	402	`w = w.next().next();`
403	`do`	403	`do`
404	`{`	404	`{`
405	`if(status[w.vertex()] == 1 && w.next().halfedge() != w0.halfedge() &&`	405	`if(status[w.vertex()] == 1 && w.next().halfedge() != w0.halfedge() &&`
406	`same_level(fun[w0.vertex()],fun[w.vertex()]))`	406	`same_level(fun[w0.vertex()],fun[w.vertex()]))`
407	`{`	407	`{`
408	`m.split_face_by_edge(*fid, w0.vertex(), w.vertex());`	408	`m.split_face_by_edge(*fid, w0.vertex(), w.vertex());`
409	`did_work = true;`	409	`did_work = true;`
410	`break;`	410	`break;`
411	`}`	411	`}`
412	`w = w.next();`	412	`w = w.next();`
413	`}`	413	`}`
414	`while(!w.full_circle());`	414	`while(!w.full_circle());`
415	`break;`	415	`break;`
416	`}`	416	`}`
417	`}`	417	`}`
418	`while(did_work);`	418	`while(did_work);`
419		419
420		420
421		421
422		422
423	`// ----------------------------`	423	`// ----------------------------`
424	`cout << "Numbering the level sets" << endl;`	424	`cout << "Numbering the level sets" << endl;`
425		425
426	`float max_length = 0;`	426	`float max_length = 0;`
427	`double max_length_fun = 0;`	427	`double max_length_fun = 0;`
428	`int max_length_id =-1;`	428	`int max_length_id =-1;`
429	`HalfEdgeID max_length_h;`	429	`HalfEdgeID max_length_h;`
430	`HalfEdgeAttributeVector<int> level_set_id(m.allocated_halfedges(), 0);`	430	`HalfEdgeAttributeVector<int> level_set_id(m.allocated_halfedges(), 0);`
431	`VertexAttributeVector<int> level_set_id_vertex(m.allocated_vertices(), 0);`	431	`VertexAttributeVector<int> level_set_id_vertex(m.allocated_vertices(), 0);`
432	`int no_id=1;`	432	`int no_id=1;`
433	`for(HalfEdgeIDIterator hid = m.halfedges_begin(); hid != m.halfedges_end(); ++hid)`	433	`for(HalfEdgeIDIterator hid = m.halfedges_begin(); hid != m.halfedges_end(); ++hid)`
434	`{`	434	`{`
435	`Walker w = m.walker(*hid);`	435	`Walker w = m.walker(*hid);`
436	`if(status[w.vertex()] == 1 && status[w.opp().vertex()] == 1 &&`	436	`if(status[w.vertex()] == 1 && status[w.opp().vertex()] == 1 &&`
437	`same_level(fun[w.vertex()], fun[w.opp().vertex()]) &&`	437	`same_level(fun[w.vertex()], fun[w.opp().vertex()]) &&`
438	`level_set_id[w.halfedge()] == 0)`	438	`level_set_id[w.halfedge()] == 0)`
439	`{`	439	`{`
440	`float level_set_length = 0;`	440	`float level_set_length = 0;`
441	`while(level_set_id[w.halfedge()] != no_id)`	441	`while(level_set_id[w.halfedge()] != no_id)`
442	`{`	442	`{`
443	`level_set_length += length(m,w.halfedge());`	443	`level_set_length += length(m,w.halfedge());`
444	`level_set_id[w.halfedge()] = no_id;`	444	`level_set_id[w.halfedge()] = no_id;`
445	`level_set_id[w.opp().halfedge()] = no_id;`	445	`level_set_id[w.opp().halfedge()] = no_id;`
446	`level_set_id_vertex[w.vertex()] = no_id;`	446	`level_set_id_vertex[w.vertex()] = no_id;`
447	`w = w.next();`	447	`w = w.next();`
448	`while(status[w.vertex()] != 1 \|\| !same_level(fun[w.vertex()], fun[w.opp().vertex()]))`	448	`while(status[w.vertex()] != 1 \|\| !same_level(fun[w.vertex()], fun[w.opp().vertex()]))`
449	`w = w.circulate_vertex_cw();`	449	`w = w.circulate_vertex_cw();`
450	`}`	450	`}`
451	`if(level_set_length > max_length)`	451	`if(level_set_length > max_length)`
452	`{`	452	`{`
453	`max_length = level_set_length;`	453	`max_length = level_set_length;`
454	`max_length_id = no_id;`	454	`max_length_id = no_id;`
455	`max_length_h = w.halfedge();`	455	`max_length_h = w.halfedge();`
456	`max_length_fun = fun[w.vertex()];`	456	`max_length_fun = fun[w.vertex()];`
457	`}`	457	`}`
458	`++no_id;`	458	`++no_id;`
459	`}`	459	`}`
460	`}`	460	`}`
461	`cout << max_length << " " << max_length_id << endl;`	461	`cout << max_length << " " << max_length_id << endl;`
462	`cout << "Number of level sets : " << (no_id-1);`	462	`cout << "Number of level sets : " << (no_id-1);`
463	`// ----------------------------`	463	`// ----------------------------`
464		464
465	`shortest_edge_triangulate(m);`	465	`shortest_edge_triangulate(m);`
466		466
467		467
468	`// -------------`	468	`// -------------`
469	`cout << "Parametrize level sets " << endl;`	469	`cout << "Parametrize level sets " << endl;`
470		470
471	`queue<HalfEdgeID> hq;`	471	`queue<HalfEdgeID> hq;`
472	`HalfEdgeAttributeVector<int> touched(m.no_halfedges(), 0);`	472	`HalfEdgeAttributeVector<int> touched(m.no_halfedges(), 0);`
473	`Walker w = m.walker(max_length_h);`	473	`Walker w = m.walker(max_length_h);`
474	`float param = 0;`	474	`float param = 0;`
475	`do`	475	`do`
476	`{`	476	`{`
477	`parametrization[w.opp().vertex()] = 2.0 * M_PI * param / max_length;`	477	`parametrization[w.opp().vertex()] = 2.0 * M_PI * param / max_length;`
478	`param += length(m, w.halfedge());`	478	`param += length(m, w.halfedge());`
479	`w = w.next();`	479	`w = w.next();`
480	`while(level_set_id[w.halfedge()] != max_length_id)`	480	`while(level_set_id[w.halfedge()] != max_length_id)`
481	`w = w.circulate_vertex_cw();`	481	`w = w.circulate_vertex_cw();`
482	`hq.push(w.halfedge());`	482	`hq.push(w.halfedge());`
483	`hq.push(w.opp().halfedge());`	483	`hq.push(w.opp().halfedge());`
484	`touched[w.halfedge()] = 1;`	484	`touched[w.halfedge()] = 1;`
485	`touched[w.opp().halfedge()] = 1;`	485	`touched[w.opp().halfedge()] = 1;`
486	`}`	486	`}`
487	`while(w.halfedge() != max_length_h);`	487	`while(w.halfedge() != max_length_h);`
488		488
489	`while(!hq.empty())`	489	`while(!hq.empty())`
490	`{`	490	`{`
491	`HalfEdgeID h = hq.front();`	491	`HalfEdgeID h = hq.front();`
492	`hq.pop();`	492	`hq.pop();`
493	`if(!touched[w.next().halfedge()])`	493	`if(!touched[w.next().halfedge()])`
494	`{`	494	`{`
495	`extend_fun2(m, h, fun, parametrization);`	495	`extend_fun2(m, h, fun, parametrization);`
496	`touched[w.next().halfedge()] = 1;`	496	`touched[w.next().halfedge()] = 1;`
497	`touched[w.prev().halfedge()] = 1;`	497	`touched[w.prev().halfedge()] = 1;`
498	`Walker w = m.walker(h);`	498	`Walker w = m.walker(h);`
499	`if(!touched[w.next().opp().halfedge()])`	499	`if(!touched[w.next().opp().halfedge()])`
500	`{`	500	`{`

Subversion Repositories gelsvn

(root)/trunk/apps/MeshEdit/polarize.cpp – Rev 596 → 607