WebSVN – gelsvn – Diff – /branches/cpp11-devel/apps/MeshEdit/polarize.cpp

 //  Copyright 2012 __MyCompanyName__. All rights reserved.
 //
 #include <queue>
 #include "polarize.h"
+#include <CGLA/Vec2d.h>
+#include <LinAlg/LapackFunc.h>
 #include <HMesh/triangulate.h>
+#include <HMesh/obj_save.h>
 #include <HMesh/curvature.h>
 #include <HMesh/quadric_simplify.h>
 #include <HMesh/mesh_optimization.h>
 #include <HMesh/smooth.h>
 using namespace CGLA;
 using namespace std;
 using namespace HMesh;
+inline bool same_level(float a, float b) {return abs(a-b) < 0.00001;}
+struct LevelSetInfo
+{
+    int id;
+    int no_vertices;
+    Vec3d avg_pos;
+    double rad;
+    double length;
+    int components;
+    HalfEdgeID h;
+    double fun_value;
+    void print()
+    {
+        cout
+        << " id : " << id
+        << " no vertices : " << no_vertices
+        << " avg pos : " << avg_pos
+        << " rad : " << rad
+        << " length : " << length
+        << " # : " << components
+        << " F : " << fun_value << endl;
+    }
+};
+void compute_edge_weights(Manifold& m, HalfEdgeAttributeVector<double>& edge_weights, FaceAttributeVector<int>& included)
+{
+    edge_weights = HalfEdgeAttributeVector<double>(m.allocated_halfedges(), 0);
+    for(FaceIDIterator f = m.faces_begin(); f != m.faces_end(); ++f)
+        if(included[*f])
+        {
+            for(Walker wv = m.walker(*f); !wv.full_circle(); wv = wv.circulate_face_ccw())
+            {
+                HalfEdgeID h = wv.halfedge();
+                Vec3d p1(m.pos(wv.vertex()));
+                Vec3d p2(m.pos(wv.next().vertex()));
+                Vec3d p0(m.pos(wv.opp().vertex()));
+                double ang = acos(min(1.0, max(-1.0, dot(normalize(p1-p0), normalize(p2-p0)))));
+                double ang_opp = acos(min(1.0, max(-1.0, dot(normalize(p2-p1), normalize(p0-p1)))));
+                double l = (p1-p0).length();
+                edge_weights[h] += tan(ang/2) / l;
+                edge_weights[wv.opp().halfedge()] += tan(ang_opp/2) / l;
+            }
+        }
+}
+void smooth_fun(Manifold& m,
+                    VertexAttributeVector<int>& nailed,
+                    VertexAttributeVector<double>& fun)
+{
+    HalfEdgeAttributeVector<double> edge_weights;
+    FaceAttributeVector<int> included(m.allocated_faces(),1);
+    compute_edge_weights(m,edge_weights, included);
+    VertexAttributeVector<double> new_fun(m.no_vertices());
+    for(int i = 0; i < 1000; ++i)
+    {
+        for(VertexIDIterator v = m.vertices_begin(); v != m.vertices_end(); ++v)
+            if(!nailed[*v])
+            {
+                double w_sum = 0;
+                new_fun[*v] = 0;
+                for(Walker wv = m.walker(*v); !wv.full_circle(); wv = wv.circulate_vertex_ccw())
+                {
+                    double w = edge_weights[wv.halfedge()];
+                    new_fun[*v] += w * fun[wv.vertex()];
+                    w_sum += w;
+                }
+                new_fun[*v] /= w_sum;
+            }
+            else
+                new_fun[*v] = fun[*v];
+            fun = new_fun;
+    }
+}
+void segment_manifold(Manifold& m, HalfEdgeAttributeVector<int>& ls_id,
+                      FaceAttributeVector<int>& segmentation,
+                      vector<vector<int>>& boundaries)
+{
+    segmentation = FaceAttributeVector<int>(m.no_faces(), -1);
+    int SEG_NO = 0;
+    for(auto f = m.faces_begin(); f != m.faces_end(); ++f)
+    {
+        if (segmentation[*f] == -1)
+        {
+            queue<FaceID> q;
+            q.push(*f);
+            vector<int> bound(0);
+            while (!q.empty())
+            {
+                FaceID face = q.front();
+                q.pop();
+                segmentation[face] = SEG_NO;
+                for(Walker w = m.walker(face); !w.full_circle(); w = w.next())
+                {
+                    if(ls_id[w.halfedge()] == -1)
+                    {
+                        FaceID fopp = w.opp().face();
+                        if(segmentation[fopp] == -1) {
+                            q.push(fopp);
+                        }
+                    }
+                    else {
+                        bound.push_back(ls_id[w.halfedge()]);
+                    }
+                }
+            }
+            sort(begin(bound), end(bound));
+            auto new_end = unique(begin(bound), end(bound));
+            boundaries.push_back(vector<int>(begin(bound),new_end));
+            cout << "Boundaries of seg# " << SEG_NO << " : ";
+            for(auto be:boundaries[SEG_NO])
+                cout << be << " ";
+            cout << endl;
+            SEG_NO += 1;
+        }
+    }
+}
 void shortest_edge_triangulate_face(Manifold& m, FaceID f0, VertexAttributeVector<int>& level_set_id_vertex)
+{
     queue<FaceID> face_queue;
     face_queue.push(f0);
         // 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
+{
     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()]));
     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();
         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);
     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 n = normal(m, f);
     Vec3d gsum(0.0);
     for(Walker w = m.walker(f); !w.full_circle(); w = w.next())
+    {
+        VertexID avid = w.opp().vertex();
+        VertexID bvid = w.vertex();
+        VertexID cvid = w.next().vertex();
-        Vec3d gdir = normalize(cross(n, m.pos(w.vertex()) - m.pos(w.opp().vertex())));
+        Vec3d gdir = normalize(cross(n, m.pos(bvid) - m.pos(avid)));
-        float l = dot(m.pos(w.next().vertex())-m.pos(w.vertex()), gdir);
+        double l = dot(m.pos(cvid)-m.pos(avid), gdir);
-        gdir *= fun[w.next().vertex()]/l;
+        gdir *= fun[cvid]/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,
+double solve_for_orthogonal_gradients(HMesh::Manifold& m, HMesh::VertexAttributeVector<double>& fun,
-                 HMesh::VertexAttributeVector<double>& fun2)
+                                     HalfEdgeID h, double beta, double gamma)
+{
+    Vec3d n = normal(m, m.walker(h).face());
+    Vec3d uvw[3];
+    Vec3d abc;
     Walker w = m.walker(h);
-    FaceID f = w.face();
-    Vec3d a = m.pos(w.opp().vertex());
+    VertexID vid[] = {w.next().vertex(), w.opp().vertex(),w.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)
+    for(int i=0; !w.full_circle(); w = w.next(),++  i)
-{
-    Vec3d avg_pos(0);
-    float asum = 0;
+    {
-    for(Walker w = m.walker(v); !w.full_circle(); w = w.circulate_vertex_cw()){
+        uvw[i] = normalize(cross(n, m.pos(vid[(i+1)%3]) - m.pos(vid[(i+2)%3])));
-        float  a = barycentric_area(m, w.vertex());
-        avg_pos += a * m.pos(w.vertex());
+        uvw[i] /= dot(m.pos(vid[i])-m.pos(vid[(i+1)%3]), uvw[i]);
-        asum += a;
+        abc[i] = fun[vid[i]];
+    }
+    Vec3d A(dot(uvw[0],uvw[0]),dot(uvw[0],uvw[1]),dot(uvw[0],uvw[2]));
+    Vec3d B(dot(uvw[1],uvw[0]),dot(uvw[1],uvw[1]),dot(uvw[1],uvw[2]));
+    Vec3d C(dot(uvw[2],uvw[0]),dot(uvw[2],uvw[1]),dot(uvw[2],uvw[2]));
-    return avg_pos / asum - m.pos(v);
+    return - (beta * dot(B, abc) + gamma * dot(C, abc))/dot(A, abc);
+}
-void aw_laplacian_smooth(Manifold& m, float t)
+void recompute_circle_param(Manifold& m, VertexID v, VertexAttributeVector<double>& fun, VertexAttributeVector<Vec2d>& circle_pos, double f_origin)
+{
+    Vec2d c_new(0);
-    VertexAttributeVector<Vec3d> pos(m.allocated_vertices());
+    for(auto w = m.walker(v); !w.full_circle(); w = w.circulate_vertex_ccw())
+    {
-    for(VertexIDIterator v = m.vertices_begin(); v != m.vertices_end(); ++v){
-        if(!boundary(m, *v))
+        Walker w_base = w.next();
-        {
-            Vec3d n = normal(m, *v);
+        Vec2d X = circle_pos[w_base.opp().vertex()];
-            Vec3d l = laplacian(m, *v);
+        Vec2d Y(-X[1], X[0]);
-            if(sqr_length(n) > 0.8)
+        Vec2d vec = circle_pos[w_base.vertex()];
-                l -= n * dot(n,l);
+        double angle = atan2(max(-1.0,min(1.0,dot(vec,Y))),
-            pos[*v] =  t * l + m.pos(*v);
+                             max(-1.0,min(1.0,dot(vec,X))));
-        }
-    }
-    for(VertexIDIterator v = m.vertices_begin(); v != m.vertices_end(); ++v){
+        double new_angle = solve_for_orthogonal_gradients(m, fun, w_base.halfedge(), 0, angle);
-        if(!boundary(m, *v))
+        new_angle = max(-M_PI/4, min(M_PI/4,new_angle));
-            m.pos(*v) = pos[*v];
+        c_new += cos(new_angle)*X + sin(new_angle)*Y;
+    }
+    circle_pos[v] = normalize(c_new);
+}
-void polarize_mesh(Manifold& m, VertexAttributeVector<double>& fun, double vmin, double vmax, const int divisions, VertexAttributeVector<double>& parametrization)
+void polarize_mesh(Manifold& m, VertexAttributeVector<double>& fun, double vmin, double vmax, const int divisions, VertexAttributeVector<Vec2d>& parametrization)
+{
     vmax -= 0.01 * (vmax-vmin);
     float interval = (vmax-vmin)/divisions;
     VertexAttributeVector<int> status(m.allocated_vertices(), 0);
                     break;
+                }
+    }
     while(did_work);
+    shortest_edge_triangulate(m);
     // ----------------------------
     cout << "Numbering the level sets" << endl;
-    float max_length = 0;
-    double max_length_fun = 0;
+    VertexAttributeVector<int> nailed(m.no_vertices(),0);
-    int max_length_id =-1;
-    HalfEdgeID max_length_h;
+    vector<LevelSetInfo> level_set_info;
-    HalfEdgeAttributeVector<int> level_set_id(m.allocated_halfedges(), 0);
+    HalfEdgeAttributeVector<int> level_set_id(m.allocated_halfedges(), -1);
-    VertexAttributeVector<int> level_set_id_vertex(m.allocated_vertices(), 0);
+    VertexAttributeVector<int> level_set_id_vertex(m.allocated_vertices(), -1);
+    VertexAttributeVector<Vec2d> circle_pos;
+    for(auto vid = m.vertices_begin(); vid != m.vertices_end(); ++vid)
+        circle_pos[*vid] = Vec2d(0,0);
-    int no_id=1;
+    int no_id=0;
     for(HalfEdgeIDIterator hid = m.halfedges_begin(); hid != m.halfedges_end(); ++hid)
+    {
         Walker w = m.walker(*hid);
         if(status[w.vertex()] == 1 && status[w.opp().vertex()] == 1 &&
            same_level(fun[w.vertex()], fun[w.opp().vertex()]) &&
-           level_set_id[w.halfedge()] == 0)
+           level_set_id[w.halfedge()] == -1)
+        {
+            LevelSetInfo lsi;
+            lsi.id = no_id;
+            lsi.fun_value = fun[w.vertex()];
+            lsi.components = 1;
+            lsi.no_vertices = 0;
             float level_set_length = 0;
             while(level_set_id[w.halfedge()] != no_id)
+            {
                 level_set_length += length(m,w.halfedge());
                 level_set_id[w.halfedge()] = no_id;
                 level_set_id[w.opp().halfedge()] = no_id;
                 level_set_id_vertex[w.vertex()] = no_id;
                 w = w.next();
                 while(status[w.vertex()] != 1 || !same_level(fun[w.vertex()], fun[w.opp().vertex()]))
                     w = w.circulate_vertex_cw();
+                lsi.no_vertices += 1;
-            }
+            }
+            lsi.h = w.halfedge();
-            if(level_set_length > max_length)
+            lsi.length = level_set_length;
+            double param = 0;
+            do
+            {
-                max_length = level_set_length;
+                double angle = 2.0 * M_PI * param / lsi.length;
+                circle_pos[w.opp().vertex()] = Vec2d(cos(angle), sin(angle));
-                max_length_id = no_id;
+                nailed[w.opp().vertex()] = 1;
-                max_length_h = w.halfedge();
+                param += length(m, w.halfedge());
+                w = w.next();
+                while(status[w.vertex()] != 1 || !same_level(fun[w.vertex()], fun[w.opp().vertex()]))
-                max_length_fun = fun[w.vertex()];
+                    w = w.circulate_vertex_cw();
+            }
+            while(w.halfedge() != lsi.h);
+            level_set_info.push_back(lsi);
+            level_set_info.back().print();
             ++no_id;
+            assert(level_set_info.size()==no_id);
+        }
+    }
-    cout << max_length << " " << max_length_id << endl;
-    cout << "Number of level sets : " << (no_id-1);
-    // ----------------------------
+    cout << "ids" << no_id << " " << level_set_info.size() << endl;
+    for(int i=0;i<level_set_info.size(); ++i)
+        for(int j=0;j<level_set_info.size(); ++j)
+            if(i != j && same_level(level_set_info[i].fun_value, level_set_info[j].fun_value))
+                level_set_info[i].components += 1;
-    shortest_edge_triangulate(m);
+    HalfEdgeID min_length_h;
+    double min_length_fun;
+    double min_length=FLT_MAX;
+    int min_length_id = -1;
+    for(int i=0;i<level_set_info.size(); ++i)
+    {
+        LevelSetInfo& lsi = level_set_info[i];
+        if(lsi.components==1 && lsi.length<min_length)
+        {
+            min_length_h = lsi.h;
+            min_length_fun = lsi.fun_value;
+            min_length = lsi.length;
+            min_length_id = lsi.id;
+        }
+    }
+    smooth_fun(m, nailed, fun);
+    FaceAttributeVector<int> segmentation;
-    // -------------
+    vector<vector<int>> boundaries;
-    cout << "Parametrize level sets " << endl;
+    segment_manifold(m, level_set_id, segmentation, boundaries);
-    queue<HalfEdgeID> hq;
+    VertexAttributeVector<Vec3d> cylinder_pos;
-    HalfEdgeAttributeVector<int> touched(m.no_halfedges(), 0);
+    for(auto vid = m.vertices_begin();vid != m.vertices_end(); ++vid)
-    Walker w = m.walker(max_length_h);
+        if (nailed[*vid])//level_set_id_vertex[*vid] == min_length_id)
+            cylinder_pos[*vid] = Vec3d(circle_pos[*vid][0], circle_pos[*vid][1],min_length_fun);
-    float param = 0;
+        else
+            cylinder_pos[*vid] = Vec3d(0,0,fun[*vid]);
-    do
+    HalfEdgeAttributeVector<double> edge_weights(m.allocated_halfedges(), 0);
+    FaceAttributeVector<int> included(m.allocated_faces(), 1);
+    for(auto fid = m.faces_begin(); fid != m.faces_end(); ++fid)
+    {
-        parametrization[w.opp().vertex()] = 2.0 * M_PI * param / max_length;
-        param += length(m, w.halfedge());
+        if(boundaries[segmentation[*fid]].size() == 1)
-        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;
+            included[*fid] = 0;
-        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);
+    compute_edge_weights(m,edge_weights,included);
-            touched[w.next().halfedge()] = 1;
+    VertexAttributeVector<Vec3d> new_cylinder_pos(m.no_vertices());
-            touched[w.prev().halfedge()] = 1;
-            Walker w = m.walker(h);
+    for(int i = 0; i < 1500; ++i)
-            if(!touched[w.next().opp().halfedge()])
-            {
+    {
-                touched[w.next().opp().halfedge()] = 1;
+        for(VertexIDIterator v = m.vertices_begin(); v != m.vertices_end(); ++v)
-                hq.push(w.next().opp().halfedge());
-            }
-            if(!touched[w.prev().opp().halfedge()])
+            if(!nailed[*v])//level_set_id_vertex[*v] != min_length_id)
+            {
+//                recompute_circle_param(m, *v, fun, circle_pos, min_length_fun);
+                double w_sum = 0;
+                new_cylinder_pos[*v] = Vec3d(0);
+                for(Walker wv = m.walker(*v); !wv.full_circle(); wv = wv.circulate_vertex_ccw())
+                {
-                touched[w.prev().opp().halfedge()] = 1;
+                    double w = edge_weights[wv.halfedge()];
+                    new_cylinder_pos[*v] += w * cylinder_pos[wv.vertex()];
+                    w_sum += w;
+                }
+                if(w_sum> 0.0)
+                    new_cylinder_pos[*v] /= w_sum;
+                Vec2d v2d(new_cylinder_pos[*v][0], new_cylinder_pos[*v][1]);
+                double l =length(v2d);
+                if(l>0.1) {
-                hq.push(w.prev().opp().halfedge());
+                    new_cylinder_pos[*v][0] /= l;
+                    new_cylinder_pos[*v][1] /= l;
+                }
+                if(level_set_id_vertex[*v] != -1)
+                    new_cylinder_pos[*v][2] = fun[*v];
+            }
-        }
+            else
+                new_cylinder_pos[*v] = cylinder_pos[*v];
+        for(auto v = m.vertices_begin();v != m.vertices_end(); ++v)
+            cylinder_pos[*v] = 0.9*cylinder_pos[*v]+ 0.1*new_cylinder_pos[*v];
+    }
+  for(auto vid = m.vertices_begin();vid != m.vertices_end(); ++vid)
+        parametrization[*vid] = Vec2d(cylinder_pos[*vid][0],cylinder_pos[*vid][1]);
+    return;
+//    for(auto vid = m.vertices_begin();vid != m.vertices_end(); ++vid)
+//        m.pos(*vid) = cylinder_pos[*vid];// Vec2d(cylinder_pos[*vid][0],cylinder_pos[*vid][1]);
+//    vector<pair<double,VertexID>> verts;
+//    for(auto vid = m.vertices_begin();vid != m.vertices_end(); ++vid)
+//        if(!same_level(min_length_fun, fun[*vid]))
+//            verts.push_back(pair<double,VertexID>(abs(fun[*vid]-min_length_fun),*vid));
+//    sort(verts.begin(), verts.end());
+//
+//    for(auto p:verts)
+//    {
+//        VertexID vid = p.second;
+//        recompute_circle_param(m, vid, fun, circle_pos, min_length_fun);
+//    }
+//    parametrization = circle_pos;
+//    return;
+//    // nail max length level set to circle.
+//    Walker w = m.walker(min_length_h);
+//    priority_queue<pair<double, HalfEdgeID>> hq;
+//    float param = 0;
+//    do
+//    {
+//        hq.push(pair<double,HalfEdgeID>(-abs(fun[w.next().vertex()]-min_length_fun),w.halfedge()));
+//        hq.push(pair<double,HalfEdgeID>(-abs(fun[w.opp().next().vertex()]-min_length_fun),w.opp().halfedge()));
+//
+//        nailed[w.opp().vertex()] = 1;
+//        w = w.next();
+//        while(level_set_id[w.halfedge()] != min_length_id)
+//            w = w.circulate_vertex_cw();
+//    }
+//    while(w.halfedge() != min_length_h);
+//    parametrize(m, nailed, circle_pos, parametrization, fun, min_length_fun);
+//    propagate_values(m, nailed, circle_pos, parametrization, level_set_id_vertex, level_set_info, fun, min_length_fun);
+////    for(int i=0;i<level_set_info.size(); ++i)
+////        level_set_info[i].print();
+//
+//    VertexAttributeVector<Vec3d> circle_center_attrib(m.no_vertices(),Vec3d(0,0,0));
+//    for(auto vid = m.vertices_begin(); vid != m.vertices_end(); ++vid)
+//    {
+//        if(level_set_id_vertex[*vid]>-1)
+//        {
+//            LevelSetInfo& lsi = level_set_info[level_set_id_vertex[*vid]];
+//            circle_center_attrib[*vid] = lsi.avg_pos;
+//            nailed[*vid] = 1;
+//        }
+//    }
+////    propagate(m, nailed, circle_center_attrib);
+//    Manifold m2 = m;
+//    for(auto vid = m2.vertices_begin(); vid != m2.vertices_end(); ++vid)
+//    {
+//        m2.pos(*vid) = circle_pos[*vid];
+//        m2.pos(*vid)[2] = 0.01 * fun[*vid];
+//    }
+//    obj_save("blob.obj", m2);
+//    for(auto vid = m.vertices_begin(); vid != m.vertices_end(); ++vid)
+//    {
+//        Vec3d p =circle_pos[*vid];//-circle_center_attrib[*vid];
+//        parametrization[*vid] = Vec2d(p[0],p[1]);
+//    }
-    return;
     // ----------------------------
     cout << "Remove vertices not on level set curves" << endl;
     vector<VertexID> vid_vec;
     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;
+//    VertexAttributeVector<Vec3d> recalled_positions;
-    for(VertexIDIterator vid = m.vertices_begin(); vid != m.vertices_end(); ++vid)
+//    for(VertexIDIterator vid = m.vertices_begin(); vid != m.vertices_end(); ++vid)
-        recalled_positions[*vid] = m.pos(*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];
+        double v = dot(m.pos(*vid),Vec3d(0.0,1,0.00));
         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;
-//    }

Subversion Repositories gelsvn

(root)/branches/cpp11-devel/apps/MeshEdit/polarize.cpp – Rev 630 → 631