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



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

        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.

    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()]));
    }
};

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

    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;

 
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;

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

    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(int 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;

    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()]))

    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;

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

        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());
            }

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

        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;

        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()]) &&

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

        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 &&

    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());
    }

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

    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()]))

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

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

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

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

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

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

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

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

Rev 586	Rev 587
Line 31...	Line 31...
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(HalfEdgeWalker w = m.halfedgewalker(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());`
Line 105...	Line 105...
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(HalfEdgeWalker w = m.halfedgewalker(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.`
Line 180...	Line 180...
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	`HalfEdgeWalker w = m.halfedgewalker(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	`};`
Line 196...	Line 196...
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	`HalfEdgeWalker w = m.halfedgewalker(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();`
Line 242...	Line 242...
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	`HalfEdgeWalker w = m.halfedgewalker(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	`Vec3f 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 Vec3f(0.0);`	258	`return Vec3d(0.0);`
259		259
260		260
261	`Vec3f n = normal(m, f);`	261	`Vec3d n = normal(m, f);`
262		262
263	`Vec3f gsum(0.0);`	263	`Vec3d gsum(0.0);`
264	`for(HalfEdgeWalker w = m.halfedgewalker(f); !w.full_circle(); w = w.next())`	264	`for(Walker w = m.walker(f); !w.full_circle(); w = w.next())`
265	`{`	265	`{`
266	`Vec3f 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;`
Line 285...	Line 285...
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	`HalfEdgeWalker w = m.halfedgewalker(h);`	290	`Walker w = m.walker(h);`
291	`FaceID f = w.face();`	291	`FaceID f = w.face();`
292	`Vec3f a = m.pos(w.opp().vertex());`	292	`Vec3d a = m.pos(w.opp().vertex());`
293	`Vec3f b = m.pos(w.vertex());`	293	`Vec3d b = m.pos(w.vertex());`
294	`Vec3f c = m.pos(w.next().vertex());`	294	`Vec3d c = m.pos(w.next().vertex());`
295	`Vec3f g = grad(m, fun, f);`	295	`Vec3d g = grad(m, fun, f);`
296	`Vec3f n = normal(m, f);`	296	`Vec3d n = normal(m, f);`
297	`Vec3f 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;`
Line 311...	Line 311...
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	`Vec3f laplacian(const Manifold& m, VertexID v)`	316	`Vec3d laplacian(const Manifold& m, VertexID v)`
317	`{`	317	`{`
318	`Vec3f avg_pos(0);`	318	`Vec3d avg_pos(0);`
319	`float asum = 0;`	319	`float asum = 0;`
320		320
321	`for(HalfEdgeWalker w = m.halfedgewalker(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<Vec3f> 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	`Vec3f n = normal(m, *v);`	336	`Vec3d n = normal(m, *v);`
337	`Vec3f 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	`}`
Line 359...	Line 359...
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	`HalfEdgeWalker w = m.halfedgewalker(*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(int i = 0; i<hidvec.size(); ++i)`
370	`{`	370	`{`
371	`HalfEdgeWalker w = m.halfedgewalker(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	`Vec3f 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	`Vec3f 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;`
Line 392...	Line 392...
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(HalfEdgeWalker w = m.halfedgewalker(*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	`HalfEdgeWalker 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()]))`
Line 430...	Line 430...
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	`HalfEdgeWalker w = m.halfedgewalker(*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;`
Line 468...	Line 468...
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	`HalfEdgeWalker w = m.halfedgewalker(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());`
Line 493...	Line 493...
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	`HalfEdgeWalker w = m.halfedgewalker(h);`	498	`Walker w = m.walker(h);`
499	`if(!touched[w.next().opp().halfedge()])`	499	`if(!touched[w.next().opp().halfedge()])`
500	`{`	500	`{`
501	`touched[w.next().opp().halfedge()] = 1;`	501	`touched[w.next().opp().halfedge()] = 1;`
502	`hq.push(w.next().opp().halfedge());`	502	`hq.push(w.next().opp().halfedge());`
503	`}`	503	`}`
Line 532...	Line 532...
532	`for(FaceIDIterator fid = m.faces_begin(); fid != m.faces_end(); ++fid)`	532	`for(FaceIDIterator fid = m.faces_begin(); fid != m.faces_end(); ++fid)`
533	`if(no_edges(m, *fid) > 3)`	533	`if(no_edges(m, *fid) > 3)`
534	`shortest_edge_triangulate_face(m, *fid, level_set_id_vertex);`	534	`shortest_edge_triangulate_face(m, *fid, level_set_id_vertex);`
535		535
536		536
537	`VertexAttributeVector<Vec3f> recalled_positions;`	537	`VertexAttributeVector<Vec3d> recalled_positions;`
538	`for(VertexIDIterator vid = m.vertices_begin(); vid != m.vertices_end(); ++vid)`	538	`for(VertexIDIterator vid = m.vertices_begin(); vid != m.vertices_end(); ++vid)`
539	`recalled_positions[vid] = m.pos(vid);`	539	`recalled_positions[vid] = m.pos(vid);`
540		540
541		541
542	`TriangleQuality tq_energy(level_set_id_vertex);`	542	`TriangleQuality tq_energy(level_set_id_vertex);`
Line 553...	Line 553...
553	`TriangleQualityValence tqv_energy(level_set_id_vertex);`	553	`TriangleQualityValence tqv_energy(level_set_id_vertex);`
554	`TriangleQuality tq_energy(level_set_id_vertex);`	554	`TriangleQuality tq_energy(level_set_id_vertex);`
555	`priority_queue_optimization(m, tqv_energy);`	555	`priority_queue_optimization(m, tqv_energy);`
556	`priority_queue_optimization(m, tq_energy);`	556	`priority_queue_optimization(m, tq_energy);`
557		557
558	`VertexAttributeVector<Vec3f> new_pos(m.allocated_vertices(), Vec3f(0));`	558	`VertexAttributeVector<Vec3d> new_pos(m.allocated_vertices(), Vec3d(0));`
559	`VertexAttributeVector<float> new_pos_wsum(m.allocated_vertices(), 0.0);`	559	`VertexAttributeVector<float> new_pos_wsum(m.allocated_vertices(), 0.0);`
560	`for(HalfEdgeIDIterator hid = m.halfedges_begin(); hid != m.halfedges_end(); ++hid)`	560	`for(HalfEdgeIDIterator hid = m.halfedges_begin(); hid != m.halfedges_end(); ++hid)`
561	`{`	561	`{`
562	`HalfEdgeWalker w = m.halfedgewalker(*hid);`	562	`Walker w = m.walker(*hid);`
563	`if(level_set_id[w.halfedge()] != 0)`	563	`if(level_set_id[w.halfedge()] != 0)`
564	`{`	564	`{`
565	`float weight = 1.0;//sqr(valency(m,w.opp().vertex())-2);`	565	`float weight = 1.0;//sqr(valency(m,w.opp().vertex())-2);`
566	`new_pos[w.vertex()] += weight*m.pos(w.opp().vertex());`	566	`new_pos[w.vertex()] += weight*m.pos(w.opp().vertex());`
567	`new_pos_wsum[w.vertex()] += weight;`	567	`new_pos_wsum[w.vertex()] += weight;`
Line 578...	Line 578...
578	`priority_queue_optimization(m, tqv_energy);`	578	`priority_queue_optimization(m, tqv_energy);`
579	`priority_queue_optimization(m, tq_energy);`	579	`priority_queue_optimization(m, tq_energy);`
580	`vector<HalfEdgeID> hidvec;`	580	`vector<HalfEdgeID> hidvec;`
581	`for(HalfEdgeIDIterator hid = m.halfedges_begin(); hid != m.halfedges_end(); ++hid)`	581	`for(HalfEdgeIDIterator hid = m.halfedges_begin(); hid != m.halfedges_end(); ++hid)`
582	`{`	582	`{`
583	`HalfEdgeWalker w = m.halfedgewalker(*hid);`	583	`Walker w = m.walker(*hid);`
584	`if(w.halfedge() < w.opp().halfedge() &&`	584	`if(w.halfedge() < w.opp().halfedge() &&`
585	`level_set_id[w.halfedge()] != 0 &&`	585	`level_set_id[w.halfedge()] != 0 &&`
586	`(`	586	`(`
587	`((level_set_id[w.opp().halfedge()] == level_set_id[w.opp().next().halfedge()] \|\|`	587	`((level_set_id[w.opp().halfedge()] == level_set_id[w.opp().next().halfedge()] \|\|`
588	`level_set_id[w.halfedge()] == level_set_id[w.next().halfedge()]) &&`	588	`level_set_id[w.halfedge()] == level_set_id[w.next().halfedge()]) &&`
Line 594...	Line 594...
594	`hidvec.push_back(w.halfedge());`	594	`hidvec.push_back(w.halfedge());`
595	`}`	595	`}`
596		596
597	`for(int i=0;i<hidvec.size(); ++i)`	597	`for(int i=0;i<hidvec.size(); ++i)`
598	`{`	598	`{`
599	`HalfEdgeWalker w = m.halfedgewalker(hidvec[i]);`	599	`Walker w = m.walker(hidvec[i]);`
600	`VertexID v = m.split_edge(hidvec[i]);`	600	`VertexID v = m.split_edge(hidvec[i]);`
601	`level_set_id_vertex[v] = level_set_id[w.halfedge()];`	601	`level_set_id_vertex[v] = level_set_id[w.halfedge()];`
602	`level_set_id[w.prev().halfedge()] = level_set_id[w.halfedge()];`	602	`level_set_id[w.prev().halfedge()] = level_set_id[w.halfedge()];`
603	`level_set_id[w.prev().opp().halfedge()] = level_set_id[w.halfedge()];`	603	`level_set_id[w.prev().opp().halfedge()] = level_set_id[w.halfedge()];`
604	`shortest_edge_triangulate_face(m, w.face(), level_set_id_vertex);`	604	`shortest_edge_triangulate_face(m, w.face(), level_set_id_vertex);`
Line 608...	Line 608...
608	`priority_queue_optimization(m, tqv_energy);`	608	`priority_queue_optimization(m, tqv_energy);`
609	`priority_queue_optimization(m, tq_energy);`	609	`priority_queue_optimization(m, tq_energy);`
610		610
611	`for(HalfEdgeIDIterator hid = m.halfedges_begin(); hid != m.halfedges_end(); ++hid)`	611	`for(HalfEdgeIDIterator hid = m.halfedges_begin(); hid != m.halfedges_end(); ++hid)`
612	`{`	612	`{`
613	`HalfEdgeWalker w = m.halfedgewalker(*hid);`	613	`Walker w = m.walker(*hid);`
614	`if(level_set_id[w.halfedge()] != 0 &&`	614	`if(level_set_id[w.halfedge()] != 0 &&`
615	`!(level_set_id[w.opp().halfedge()] == level_set_id[w.opp().next().halfedge()] \|\|`	615	`!(level_set_id[w.opp().halfedge()] == level_set_id[w.opp().next().halfedge()] \|\|`
616	`level_set_id[w.halfedge()] == level_set_id[w.next().halfedge()]) &&`	616	`level_set_id[w.halfedge()] == level_set_id[w.next().halfedge()]) &&`
617	`valency(m, w.vertex())<5 &&`	617	`valency(m, w.vertex())<5 &&`
618	`valency(m, w.opp().vertex())<5 &&`	618	`valency(m, w.opp().vertex())<5 &&`
Line 629...	Line 629...
629	`return;`	629	`return;`
630		630
631	`priority_queue<EdgeQElem> edge_queue;`	631	`priority_queue<EdgeQElem> edge_queue;`
632	`for(HalfEdgeIDIterator hid = m.halfedges_begin(); hid != m.halfedges_end(); ++hid)`	632	`for(HalfEdgeIDIterator hid = m.halfedges_begin(); hid != m.halfedges_end(); ++hid)`
633	`{`	633	`{`
634	`HalfEdgeWalker w = m.halfedgewalker(*hid);`	634	`Walker w = m.walker(*hid);`
635	`if(w.halfedge()<w.opp().halfedge() &&`	635	`if(w.halfedge()<w.opp().halfedge() &&`
636	`level_set_id[w.halfedge()] == 0)`	636	`level_set_id[w.halfedge()] == 0)`
637	`{`	637	`{`
638	`Vec3f v = normalize(m.pos(w.vertex()) - m.pos(w.opp().vertex()));`	638	`Vec3d v = normalize(m.pos(w.vertex()) - m.pos(w.opp().vertex()));`
639		639
640	`float weight = 0;`	640	`float weight = 0;`
641	`HalfEdgeWalker wo = m.halfedgewalker(w.opp().vertex());`	641	`Walker wo = m.walker(w.opp().vertex());`
642	`for(; !w.full_circle(); w = w.circulate_vertex_ccw())`	642	`for(; !w.full_circle(); w = w.circulate_vertex_ccw())`
643	`if(level_set_id[w.halfedge()] != 0)`	643	`if(level_set_id[w.halfedge()] != 0)`
644	`{`	644	`{`
645	`Vec3f e = normalize(m.pos(w.vertex()) - m.pos(w.opp().vertex()));`	645	`Vec3d e = normalize(m.pos(w.vertex()) - m.pos(w.opp().vertex()));`
646	`weight += abs(dot(v,e));`	646	`weight += abs(dot(v,e));`
647	`}`	647	`}`
648	`for(; !wo.full_circle(); wo = wo.circulate_vertex_ccw())`	648	`for(; !wo.full_circle(); wo = wo.circulate_vertex_ccw())`
649	`if(level_set_id[wo.halfedge()] != 0)`	649	`if(level_set_id[wo.halfedge()] != 0)`
650	`{`	650	`{`
651	`Vec3f e = normalize(m.pos(wo.vertex()) - m.pos(wo.opp().vertex()));`	651	`Vec3d e = normalize(m.pos(wo.vertex()) - m.pos(wo.opp().vertex()));`
652	`weight += abs(dot(v,e));`	652	`weight += abs(dot(v,e));`
653	`}`	653	`}`
654	`edge_queue.push(EdgeQElem(weight, *hid));`	654	`edge_queue.push(EdgeQElem(weight, *hid));`
655	`}`	655	`}`
656	`}`	656	`}`
Line 658...	Line 658...
658	`while(!edge_queue.empty())`	658	`while(!edge_queue.empty())`
659	`{`	659	`{`
660	`HalfEdgeID h = edge_queue.top().he;`	660	`HalfEdgeID h = edge_queue.top().he;`
661	`edge_queue.pop();`	661	`edge_queue.pop();`
662		662
663	`HalfEdgeWalker w = m.halfedgewalker(h);`	663	`Walker w = m.walker(h);`
664	`if(no_edges(m, w.face())==3 && no_edges(m,w.opp().face())==3 &&`	664	`if(no_edges(m, w.face())==3 && no_edges(m,w.opp().face())==3 &&`
665	`!(level_set_id[w.next().halfedge()] == level_set_id[w.opp().prev().halfedge()] \|\|`	665	`!(level_set_id[w.next().halfedge()] == level_set_id[w.opp().prev().halfedge()] \|\|`
666	`level_set_id[w.prev().halfedge()] == level_set_id[w.opp().next().halfedge()]))`	666	`level_set_id[w.prev().halfedge()] == level_set_id[w.opp().next().halfedge()]))`
667	`m.merge_faces(w.face(), w.halfedge());`	667	`m.merge_faces(w.face(), w.halfedge());`
668	`}`	668	`}`
Line 673...	Line 673...
673		673
674	`do{`	674	`do{`
675	`did_work = false;`	675	`did_work = false;`
676	`for(HalfEdgeIDIterator hid = m.halfedges_begin(); hid != m.halfedges_end(); ++hid)`	676	`for(HalfEdgeIDIterator hid = m.halfedges_begin(); hid != m.halfedges_end(); ++hid)`
677	`{`	677	`{`
678	`HalfEdgeWalker w = m.halfedgewalker(*hid);`	678	`Walker w = m.walker(*hid);`
679		679
680	`if(level_set_id[w.halfedge()] != 0 && no_edges(m, w.face())==3 &&`	680	`if(level_set_id[w.halfedge()] != 0 && no_edges(m, w.face())==3 &&`
681	`precond_collapse_edge(m, w.halfedge()))`	681	`precond_collapse_edge(m, w.halfedge()))`
682	`{`	682	`{`
683	`m.collapse_edge(w.halfedge(), true);`	683	`m.collapse_edge(w.halfedge(), true);`
Line 695...	Line 695...
695	`do {`	695	`do {`
696	`++k;`	696	`++k;`
697	`did_work = false;`	697	`did_work = false;`
698	`for(HalfEdgeIDIterator hid = m.halfedges_begin(); hid != m.halfedges_end(); ++hid)`	698	`for(HalfEdgeIDIterator hid = m.halfedges_begin(); hid != m.halfedges_end(); ++hid)`
699	`{`	699	`{`
700	`HalfEdgeWalker w0 = m.halfedgewalker(*hid);`	700	`Walker w0 = m.walker(*hid);`
701		701
702	`if(level_set_id[w0.halfedge()] != 0 &&`	702	`if(level_set_id[w0.halfedge()] != 0 &&`
703	`(level_set_id[w0.next().halfedge()] == 0 && level_set_id[w0.prev().halfedge()] == 0))`	703	`(level_set_id[w0.next().halfedge()] == 0 && level_set_id[w0.prev().halfedge()] == 0))`
704	`{`	704	`{`
705		705
706	`HalfEdgeWalker w = w0;`	706	`Walker w = w0;`
707	`bool do_split = false;`	707	`bool do_split = false;`
708	`for(;!w.full_circle(); w = w.next())`	708	`for(;!w.full_circle(); w = w.next())`
709	`{`	709	`{`
710	`if(level_set_id[w.halfedge()] != 0 &&`	710	`if(level_set_id[w.halfedge()] != 0 &&`
711	`(level_set_id[w.next().halfedge()] == level_set_id[w.halfedge()]))`	711	`(level_set_id[w.next().halfedge()] == level_set_id[w.halfedge()]))`
Line 752...	Line 752...
752	`// cout << "Remove short level set edges" << endl;`	752	`// cout << "Remove short level set edges" << endl;`
753	`// float avglen=0;`	753	`// float avglen=0;`
754	`// int n=0;`	754	`// int n=0;`
755	`// for(HalfEdgeIDIterator hid = m.halfedges_begin(); hid != m.halfedges_end(); ++hid)`	755	`// for(HalfEdgeIDIterator hid = m.halfedges_begin(); hid != m.halfedges_end(); ++hid)`
756	`// {`	756	`// {`
757	`// HalfEdgeWalker w = m.halfedgewalker(*hid);`	757	`// Walker w = m.walker(*hid);`
758	`// if(level_set_id[w.halfedge()] != 0)`	758	`// if(level_set_id[w.halfedge()] != 0)`
759	`// {`	759	`// {`
760	`// avglen += length(m, w.halfedge());`	760	`// avglen += length(m, w.halfedge());`
761	`// ++n;`	761	`// ++n;`
762	`// }`	762	`// }`
Line 769...	Line 769...
769		769
770		770
771	`// vector<HalfEdgeID> hidvec;`	771	`// vector<HalfEdgeID> hidvec;`
772	`// for(HalfEdgeIDIterator hid = m.halfedges_begin(); hid != m.halfedges_end(); ++hid)`	772	`// for(HalfEdgeIDIterator hid = m.halfedges_begin(); hid != m.halfedges_end(); ++hid)`
773	`// {`	773	`// {`
774	`// HalfEdgeWalker w = m.halfedgewalker(*hid);`	774	`// Walker w = m.walker(*hid);`
775	`// if(w.halfedge() < w.opp().halfedge() &&`	775	`// if(w.halfedge() < w.opp().halfedge() &&`
776	`// level_set_id[w.halfedge()] != 0 &&`	776	`// level_set_id[w.halfedge()] != 0 &&`
777	`// valency(m, w.vertex())+valency(m,w.opp().vertex())>12)`	777	`// valency(m, w.vertex())+valency(m,w.opp().vertex())>12)`
778	`// hidvec.push_back(w.halfedge());`	778	`// hidvec.push_back(w.halfedge());`
779	`// }`	779	`// }`
780	`//`	780	`//`
781	`// for(int i=0;i<hidvec.size(); ++i)`	781	`// for(int i=0;i<hidvec.size(); ++i)`
782	`// {`	782	`// {`
783	`// HalfEdgeWalker w = m.halfedgewalker(hidvec[i]);`	783	`// Walker w = m.walker(hidvec[i]);`
784	`// VertexID v = m.split_edge(hidvec[i]);`	784	`// VertexID v = m.split_edge(hidvec[i]);`
785	`// level_set_id_vertex[v] = level_set_id[w.halfedge()];`	785	`// level_set_id_vertex[v] = level_set_id[w.halfedge()];`
786	`// level_set_id[w.prev().halfedge()] = level_set_id[w.halfedge()];`	786	`// level_set_id[w.prev().halfedge()] = level_set_id[w.halfedge()];`
787	`// level_set_id[w.prev().opp().halfedge()] = level_set_id[w.halfedge()];`	787	`// level_set_id[w.prev().opp().halfedge()] = level_set_id[w.halfedge()];`
788	`// shortest_edge_triangulate_face(m, w.face(), level_set_id_vertex);`	788	`// shortest_edge_triangulate_face(m, w.face(), level_set_id_vertex);`
Line 797...	Line 797...
797	`//// ----------------------------`	797	`//// ----------------------------`
798	`//cout << "smooth level set curves" << endl;`	798	`//cout << "smooth level set curves" << endl;`
799	`//`	799	`//`
800	`//for(int iter=0;iter<2;++iter)`	800	`//for(int iter=0;iter<2;++iter)`
801	`//{`	801	`//{`
802	`// VertexAttributeVector<Vec3f> new_pos(m.allocated_vertices(), Vec3f(0));`	802	`// VertexAttributeVector<Vec3d> new_pos(m.allocated_vertices(), Vec3d(0));`
803	`// for(HalfEdgeIDIterator hid = m.halfedges_begin(); hid != m.halfedges_end(); ++hid)`	803	`// for(HalfEdgeIDIterator hid = m.halfedges_begin(); hid != m.halfedges_end(); ++hid)`
804	`// {`	804	`// {`
805	`// HalfEdgeWalker w = m.halfedgewalker(*hid);`	805	`// Walker w = m.walker(*hid);`
806	`// if(status[w.vertex()] == 1 && status[w.opp().vertex()] == 1 &&`	806	`// if(status[w.vertex()] == 1 && status[w.opp().vertex()] == 1 &&`
807	`// same_level(fun[w.vertex()], fun[w.opp().vertex()]))`	807	`// same_level(fun[w.vertex()], fun[w.opp().vertex()]))`
808	`// {`	808	`// {`
809	`// new_pos[w.vertex()] += m.pos(w.vertex()) + m.pos(w.opp().vertex());`	809	`// new_pos[w.vertex()] += m.pos(w.vertex()) + m.pos(w.opp().vertex());`
810	`// }`	810	`// }`
Line 853...	Line 853...
853	`// did_work = false;`	853	`// did_work = false;`
854	`// priority_queue<EdgeQElem> edge_queue;`	854	`// priority_queue<EdgeQElem> edge_queue;`
855	`//`	855	`//`
856	`// for(HalfEdgeIDIterator hid = m.halfedges_begin(); hid != m.halfedges_end(); ++hid)`	856	`// for(HalfEdgeIDIterator hid = m.halfedges_begin(); hid != m.halfedges_end(); ++hid)`
857	`// {`	857	`// {`
858	`// HalfEdgeWalker w = m.halfedgewalker(*hid);`	858	`// Walker w = m.walker(*hid);`
859	`// if(status[w.opp().vertex()] == 0 && status[w.vertex()] == 1)`	859	`// if(status[w.opp().vertex()] == 0 && status[w.vertex()] == 1)`
860	`// {`	860	`// {`
861	`// float weight = (abs(fun[w.vertex()] - fun[w.opp().vertex()]*length(m, w.halfedge())) + 1e-6);`	861	`// float weight = (abs(fun[w.vertex()] - fun[w.opp().vertex()]*length(m, w.halfedge())) + 1e-6);`
862	`// edge_queue.push(EdgeQElem(weight, w.halfedge()));`	862	`// edge_queue.push(EdgeQElem(weight, w.halfedge()));`
863	`// }`	863	`// }`
Line 886...	Line 886...
886	`//priority_queue<EdgeQElem> edge_queue;`	886	`//priority_queue<EdgeQElem> edge_queue;`
887	`//HalfEdgeAttributeVector<int> time_stamp(m.allocated_halfedges(), 0);`	887	`//HalfEdgeAttributeVector<int> time_stamp(m.allocated_halfedges(), 0);`
888	`//for(HalfEdgeIDIterator hid = m.halfedges_begin(); hid != m.halfedges_end(); ++hid)`	888	`//for(HalfEdgeIDIterator hid = m.halfedges_begin(); hid != m.halfedges_end(); ++hid)`
889	`//if(level_set_id[*hid]==0)`	889	`//if(level_set_id[*hid]==0)`
890	`//{`	890	`//{`
891	`// HalfEdgeWalker w = m.halfedgewalker(*hid);`	891	`// Walker w = m.walker(*hid);`
892	`// if(w.halfedge()<w.opp().halfedge() && !(status[w.vertex()]==1 && status[w.opp().vertex()]==1))`	892	`// if(w.halfedge()<w.opp().halfedge() && !(status[w.vertex()]==1 && status[w.opp().vertex()]==1))`
893	`// 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));`	893	`// 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));`
894	`// }`	894	`// }`
895	`//`	895	`//`
896	`//shortest_edge_triangulate(m);`	896	`//shortest_edge_triangulate(m);`
Line 898...	Line 898...
898	`//while(!edge_queue.empty())`	898	`//while(!edge_queue.empty())`
899	`//{`	899	`//{`
900	`// HalfEdgeID hid = edge_queue.top().he;`	900	`// HalfEdgeID hid = edge_queue.top().he;`
901	`// if(m.in_use(hid) && time_stamp[hid]== edge_queue.top().time_stamp)`	901	`// if(m.in_use(hid) && time_stamp[hid]== edge_queue.top().time_stamp)`
902	`// {`	902	`// {`
903	`// HalfEdgeWalker w = m.halfedgewalker(hid);`	903	`// Walker w = m.walker(hid);`
904	`// HalfEdgeWalker wa = w.next();`	904	`// Walker wa = w.next();`
905	`// HalfEdgeWalker wb = w.prev().opp();`	905	`// Walker wb = w.prev().opp();`
906	`//`	906	`//`
907	`// if(m.merge_faces(w.face(), hid))`	907	`// if(m.merge_faces(w.face(), hid))`
908	`// {`	908	`// {`
909	`// cout << ".";`	909	`// cout << ".";`
910	`// if(valency(m, wa.opp().vertex())==2 && precond_collapse_edge(m, wa.halfedge()))`	910	`// if(valency(m, wa.opp().vertex())==2 && precond_collapse_edge(m, wa.halfedge()))`
Line 939...	Line 939...
939	`// // bool did_work;`	939	`// // bool did_work;`
940	`// do{`	940	`// do{`
941	`// did_work = false;`	941	`// did_work = false;`
942	`// for(HalfEdgeIDIterator hid = m.halfedges_begin(); hid != m.halfedges_end(); ++hid)`	942	`// for(HalfEdgeIDIterator hid = m.halfedges_begin(); hid != m.halfedges_end(); ++hid)`
943	`// {`	943	`// {`
944	`// HalfEdgeWalker w = m.halfedgewalker(*hid);`	944	`// Walker w = m.walker(*hid);`
945	`// if(level_set_id[w.halfedge()] != 0 &&`	945	`// if(level_set_id[w.halfedge()] != 0 &&`
946	`// valency(m, w.vertex())==3 &&`	946	`// valency(m, w.vertex())==3 &&`
947	`// valency(m, w.opp().vertex())==3 &&`	947	`// valency(m, w.opp().vertex())==3 &&`
948	`// ((level_set_id[w.next().halfedge()] == 0 &&level_set_id[w.opp().next().halfedge()] == 0) \|\|`	948	`// ((level_set_id[w.next().halfedge()] == 0 &&level_set_id[w.opp().next().halfedge()] == 0) \|\|`
949	`// (level_set_id[w.prev().halfedge()] == 0 &&level_set_id[w.opp().prev().halfedge()] == 0)) &&`	949	`// (level_set_id[w.prev().halfedge()] == 0 &&level_set_id[w.opp().prev().halfedge()] == 0)) &&`
Line 959...	Line 959...
959		959
960	`// for(VertexIDIterator vid = m.vertices_begin(); vid != m.vertices_end(); ++vid)`	960	`// for(VertexIDIterator vid = m.vertices_begin(); vid != m.vertices_end(); ++vid)`
961	`// {`	961	`// {`
962	`// bool is_max = true;`	962	`// bool is_max = true;`
963	`// bool is_min = true;`	963	`// bool is_min = true;`
964	`// HalfEdgeWalker w = m.halfedgewalker(*vid);`	964	`// Walker w = m.walker(*vid);`
965	`// float f = fun[*vid];`	965	`// float f = fun[*vid];`
966	`// for(;!w.full_circle(); w = w.circulate_vertex_ccw())`	966	`// for(;!w.full_circle(); w = w.circulate_vertex_ccw())`
967	`// {`	967	`// {`
968	`// if(fun[w.vertex()] < f)`	968	`// if(fun[w.vertex()] < f)`
969	`// is_min = false;`	969	`// is_min = false;`

Subversion Repositories gelsvn

(root)/trunk/apps/MeshEdit/polarize.cpp @ 596 – Rev 586 → 587