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

 //
 //  Created by J. Andreas Bærentzen on 18/03/12.
 //  Copyright 2012 __MyCompanyName__. All rights reserved.
 //
 #include <queue>
+#include <iomanip>
 #include "polarize.h"
 #include <CGLA/Vec2d.h>
 #include <LinAlg/LapackFunc.h>
 #include <HMesh/triangulate.h>
 using namespace std;
 using namespace HMesh;
-inline bool same_level(float a, float b) {return abs(a-b) < 0.00001;}
+inline bool same_level(double a, double b) {return abs(a-b) < 1e-10;}
 struct LevelSetInfo
+{
     int id;
                 edge_weights[h] += tan(ang/2) / l;
                 edge_weights[wv.opp().halfedge()] += tan(ang_opp/2) / l;
+            }
+        }
+}
+template<typename T>
 void smooth_fun(Manifold& m,
-                    VertexAttributeVector<int>& nailed,
+                VertexAttributeVector<int>& nailed,
-                    VertexAttributeVector<double>& fun)
+                VertexAttributeVector<T>& fun, int iter=100)
+{
     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());
+    VertexAttributeVector<T> new_fun(m.no_vertices());
-    for(int i = 0; i < 1000; ++i)
+    for(int i = 0; i < iter; ++i)
+    {
         for(VertexIDIterator v = m.vertices_begin(); v != m.vertices_end(); ++v)
             if(!nailed[*v])
+            {
                 double w_sum = 0;
-                new_fun[*v] = 0;
+                new_fun[*v] = T(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;
+        fun = new_fun;
+    }
+}
 void segment_manifold(Manifold& m, HalfEdgeAttributeVector<int>& ls_id,
                       FaceAttributeVector<int>& segmentation,
             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)
         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]])
+                   (level_set_id_vertex[verts[i]] == 0 || level_set_id_vertex[verts[i]] != level_set_id_vertex[verts[j]]) &&
+                   valency(m, verts[i])<5 && valency(m,verts[j])<5
+                   )
                     vpairs.push_back(pair<int,int>(i, j));
+            }
+        }
         if(vpairs.empty()){
         gsum += gdir;
+    }
     return gsum;
+}
+Vec2d extend_fun2(HMesh::Manifold& m, HMesh::HalfEdgeID h,
+                  HMesh::VertexAttributeVector<double>& fun,
+                  HMesh::VertexAttributeVector<Vec2d>& 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);
+    Vec2d cval;
+    if(dot_nba > 1e-5)
+    {
+        float t = dot(N, c-a)/dot_nba;
+        if(t>0 && t < 1) {
+            Vec2d aval = fun2[w.opp().vertex()];
+            Vec2d bval = fun2[w.vertex()];
+            return (1-t)*aval + t*bval;
+        }
+    }
+    return Vec2d(0);
+}
 double solve_for_orthogonal_gradients(HMesh::Manifold& m, HMesh::VertexAttributeVector<double>& fun,
-                                     HalfEdgeID h, double beta, double gamma)
+                                      HalfEdgeID h, double beta, double gamma)
+{
     Vec3d n = normal(m, m.walker(h).face());
     Vec3d uvw[3];
     Vec3d abc;
     Walker w = m.walker(h);
     VertexID vid[] = {w.next().vertex(), w.opp().vertex(),w.vertex()};
     for(int i=0; !w.full_circle(); w = w.next(),++  i)
+    {
         uvw[i] = normalize(cross(n, m.pos(vid[(i+1)%3]) - m.pos(vid[(i+2)%3])));
         uvw[i] /= dot(m.pos(vid[i])-m.pos(vid[(i+1)%3]), uvw[i]);
         abc[i] = fun[vid[i]];
         c_new += cos(new_angle)*X + sin(new_angle)*Y;
+    }
     circle_pos[v] = normalize(c_new);
+}
+template<typename T>
+void orthogonal_trajectories(Manifold& m, VertexID v, VertexAttributeVector<double>& f, double min_length_fun,
+                             VertexAttributeVector<T>& h_in, VertexAttributeVector<T>& h_out)
+{
+    h_out[v] = T(0);
+    T num(0);
+    double denom = 0.0001;
+    for(Walker w = m.walker(v); !w.full_circle(); w = w.circulate_vertex_ccw())
+    {
+        double weight = area(m, w.face());
+        Vec3d a = m.pos(v);
+        Vec3d b = m.pos(w.vertex());
+        Vec3d c = m.pos(w.next().vertex());
+        Vec3d fvec(f[v],f[w.vertex()], f[w.next().vertex()]);
+        if(sqr_length(fvec)>0)
+            fvec.normalize();
+        double f_a = fvec[0];
+        double f_b = fvec[1];
+        double f_c = fvec[2];
+        T h_b = h_in[w.vertex()];
+        T h_c = h_in[w.next().vertex()];
+        Vec3d e_a = c-b;
+        Vec3d e_b = a-c;
+        Vec3d e_c = b-a;
+        Vec3d N = normalize(cross(e_b,e_c)+cross(e_c,e_a)+cross(e_a,e_b));
+        Vec3d g_a = cross(N, normalize(e_a));
+        g_a /= dot(g_a, e_b);
+        Vec3d g_b = cross(N, normalize(e_b));
+        g_b /= dot(g_b, e_c);
+        Vec3d g_c = cross(N, normalize(e_c));
+        g_c /= dot(g_c, e_a);
+        num -= weight * (h_b * (f_a*dot(g_a, g_b) + f_b*dot(g_b,g_b) + f_c*dot(g_b,g_c)) +
+                         h_c * (f_a*dot(g_a, g_c) + f_b*dot(g_b,g_c) + f_c*dot(g_c, g_c)));
+        denom += weight * (f_a * dot(g_a,g_a)+ f_b*dot(g_a,g_b) +f_c*dot(g_a, g_c));
+    }
+    h_out[v] = num / denom;
+}
 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;
+    double interval = (vmax-vmin)/(divisions+1);
     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)
+    for(int cut_no=1; cut_no<divisions+1;++cut_no)
+    {
-        Walker w = m.walker(hidvec[i]);
-        float b = (fun[w.vertex()]- vmin)/interval;
+        double cut_val = vmin + cut_no * 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());
+        vector<HalfEdgeID> hidvec;
-        for(int j=floor_b; j>floor_a; --j)
+        for(HalfEdgeID hid : m.halfedges())
+        {
-            float t = (j-a) / (b-a);
+            Walker w = m.walker(hid);
+            double bval = fun[w.vertex()];
+            double aval = fun[w.opp().vertex()];
+            if(aval<bval && aval <= cut_val && cut_val < bval)
+            {
+                double t = (cut_val-aval)/(bval-aval);
-            Vec3d p = t * pb + (1.0-t) * m.pos(w.opp().vertex());
+                Vec3d p = t*m.pos(w.vertex()) + (1.0-t)*m.pos(w.opp().vertex());
-            VertexID v_new = m.split_edge(w.halfedge());
+                VertexID vnew = m.split_edge(hid);
-            w = w.prev();
+                m.pos(vnew) = p;
-            status[v_new] = 1;
+                status[vnew] = 1;
-            fun[v_new] = j * interval + vmin;
+                fun[vnew] = cut_val;
-            m.pos(v_new) = p;
+            }
+        }
-    }
-    bool did_work;
-    do
-    {
-        did_work = false;
-        for(FaceIDIterator fid = m.faces_begin(); fid != m.faces_end(); ++fid)
+        for(FaceID fid : m.faces())
-            for(Walker w = m.walker(*fid);!w.full_circle(); w = w.next())
+                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()])))
+                    if(status[w.vertex()] == 1 && same_level(fun[w.vertex()],cut_val))
-                {
-                    Walker w0 = w;
-                    w = w.next().next();
-                    do
+                    {
-                        if(status[w.vertex()] == 1 && w.next().halfedge() != w0.halfedge() &&
+                        Walker w0 = w;
-                           same_level(fun[w0.vertex()],fun[w.vertex()]))
+                        w = w.next().next();
+                        do
+                        {
+                            if(status[w.vertex()] == 1 &&
+                               w.next().halfedge() != w0.halfedge() &&
+                               same_level(fun[w0.vertex()],cut_val))
+                            {
-                            m.split_face_by_edge(*fid, w0.vertex(), w.vertex());
+                                m.split_face_by_edge(fid, w0.vertex(), w.vertex());
-                            did_work = true;
+                                break;
-                            break;
+                            }
+                            w = w.next();
+                        }
+                        while(!w.full_circle());
-                        w = w.next();
+                        break;
+                    }
-                    while(!w.full_circle());
-                    break;
-                }
+    }
-    while(did_work);
     shortest_edge_triangulate(m);
     // ----------------------------
     cout << "Numbering the level sets" << endl;
-    VertexAttributeVector<int> nailed(m.no_vertices(),0);
     vector<LevelSetInfo> level_set_info;
     HalfEdgeAttributeVector<int> level_set_id(m.allocated_halfedges(), -1);
     VertexAttributeVector<int> level_set_id_vertex(m.allocated_vertices(), -1);
-    VertexAttributeVector<Vec2d> circle_pos;
-    for(auto vid = m.vertices_begin(); vid != m.vertices_end(); ++vid)
-        circle_pos[*vid] = Vec2d(0,0);
     int no_id=0;
-    for(HalfEdgeIDIterator hid = m.halfedges_begin(); hid != m.halfedges_end(); ++hid)
+    for(HalfEdgeID hid : m.halfedges())
+    {
-        Walker w = m.walker(*hid);
+        Walker w = m.walker(hid);
         if(status[w.vertex()] == 1 && status[w.opp().vertex()] == 1 &&
            same_level(fun[w.vertex()], fun[w.opp().vertex()]) &&
            level_set_id[w.halfedge()] == -1)
+        {
             LevelSetInfo lsi;
             lsi.id = no_id;
             lsi.fun_value = fun[w.vertex()];
             lsi.components = 1;
             lsi.no_vertices = 0;
-            float level_set_length = 0;
+            double 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()]))
+                while(!(status[w.vertex()] == 1 && same_level(fun[w.vertex()], fun[w.opp().vertex()])))
                     w = w.circulate_vertex_cw();
                 lsi.no_vertices += 1;
+            }
             lsi.h = w.halfedge();
             lsi.length = level_set_length;
-            double param = 0;
-            do
-            {
-                double angle = 2.0 * M_PI * param / lsi.length;
-                circle_pos[w.opp().vertex()] = Vec2d(cos(angle), sin(angle));
-                nailed[w.opp().vertex()] = 1;
-                param += length(m, w.halfedge());
-                w = w.next();
-                while(status[w.vertex()] != 1 || !same_level(fun[w.vertex()], fun[w.opp().vertex()]))
-                    w = w.circulate_vertex_cw();
-            }
-            while(w.halfedge() != lsi.h);
             level_set_info.push_back(lsi);
             level_set_info.back().print();
             ++no_id;
             assert(level_set_info.size()==no_id);
+        }
+        {
             min_length_h = lsi.h;
             min_length_fun = lsi.fun_value;
             min_length = lsi.length;
             min_length_id = lsi.id;
+        }
+    }
-    smooth_fun(m, nailed, fun);
-    FaceAttributeVector<int> segmentation;
-    vector<vector<int>> boundaries;
-    segment_manifold(m, level_set_id, segmentation, boundaries);
-    VertexAttributeVector<Vec3d> cylinder_pos;
-    for(auto vid = m.vertices_begin();vid != m.vertices_end(); ++vid)
-        if (nailed[*vid])//level_set_id_vertex[*vid] == min_length_id)
-            cylinder_pos[*vid] = Vec3d(circle_pos[*vid][0], circle_pos[*vid][1],min_length_fun);
-        else
-            cylinder_pos[*vid] = Vec3d(0,0,fun[*vid]);
-    HalfEdgeAttributeVector<double> edge_weights(m.allocated_halfedges(), 0);
-    FaceAttributeVector<int> included(m.allocated_faces(), 1);
-    for(auto fid = m.faces_begin(); fid != m.faces_end(); ++fid)
-    {
-        if(boundaries[segmentation[*fid]].size() == 1)
-            included[*fid] = 0;
-    }
-    compute_edge_weights(m,edge_weights,included);
-    VertexAttributeVector<Vec3d> new_cylinder_pos(m.no_vertices());
-    for(int i = 0; i < 1500; ++i)
-    {
-        for(VertexIDIterator v = m.vertices_begin(); v != m.vertices_end(); ++v)
-            if(!nailed[*v])//level_set_id_vertex[*v] != min_length_id)
-            {
-//                recompute_circle_param(m, *v, fun, circle_pos, min_length_fun);
-                double w_sum = 0;
-                new_cylinder_pos[*v] = Vec3d(0);
-                for(Walker wv = m.walker(*v); !wv.full_circle(); wv = wv.circulate_vertex_ccw())
-                {
-                    double w = edge_weights[wv.halfedge()];
-                    new_cylinder_pos[*v] += w * cylinder_pos[wv.vertex()];
-                    w_sum += w;
-                }
-                if(w_sum> 0.0)
-                    new_cylinder_pos[*v] /= w_sum;
-                Vec2d v2d(new_cylinder_pos[*v][0], new_cylinder_pos[*v][1]);
-                double l =length(v2d);
-                if(l>0.1) {
-                    new_cylinder_pos[*v][0] /= l;
-                    new_cylinder_pos[*v][1] /= l;
-                }
-                if(level_set_id_vertex[*v] != -1)
-                    new_cylinder_pos[*v][2] = fun[*v];
-            }
-            else
-                new_cylinder_pos[*v] = cylinder_pos[*v];
-        for(auto v = m.vertices_begin();v != m.vertices_end(); ++v)
-            cylinder_pos[*v] = 0.9*cylinder_pos[*v]+ 0.1*new_cylinder_pos[*v];
-    }
-  for(auto vid = m.vertices_begin();vid != m.vertices_end(); ++vid)
-        parametrization[*vid] = Vec2d(cylinder_pos[*vid][0],cylinder_pos[*vid][1]);
-    return;
-//    for(auto vid = m.vertices_begin();vid != m.vertices_end(); ++vid)
-//        m.pos(*vid) = cylinder_pos[*vid];// Vec2d(cylinder_pos[*vid][0],cylinder_pos[*vid][1]);
-//    vector<pair<double,VertexID>> verts;
-//    for(auto vid = m.vertices_begin();vid != m.vertices_end(); ++vid)
-//        if(!same_level(min_length_fun, fun[*vid]))
-//            verts.push_back(pair<double,VertexID>(abs(fun[*vid]-min_length_fun),*vid));
-//    sort(verts.begin(), verts.end());
-//
-//    for(auto p:verts)
-//    {
-//        VertexID vid = p.second;
-//        recompute_circle_param(m, vid, fun, circle_pos, min_length_fun);
-//    }
-//    parametrization = circle_pos;
-//    return;
-//    // nail max length level set to circle.
-//    Walker w = m.walker(min_length_h);
-//    priority_queue<pair<double, HalfEdgeID>> hq;
-//    float param = 0;
-//    do
-//    {
-//        hq.push(pair<double,HalfEdgeID>(-abs(fun[w.next().vertex()]-min_length_fun),w.halfedge()));
-//        hq.push(pair<double,HalfEdgeID>(-abs(fun[w.opp().next().vertex()]-min_length_fun),w.opp().halfedge()));
-//
-//        nailed[w.opp().vertex()] = 1;
-//        w = w.next();
-//        while(level_set_id[w.halfedge()] != min_length_id)
-//            w = w.circulate_vertex_cw();
-//    }
-//    while(w.halfedge() != min_length_h);
-//    parametrize(m, nailed, circle_pos, parametrization, fun, min_length_fun);
-//    propagate_values(m, nailed, circle_pos, parametrization, level_set_id_vertex, level_set_info, fun, min_length_fun);
-////    for(int i=0;i<level_set_info.size(); ++i)
-////        level_set_info[i].print();
-//
-//    VertexAttributeVector<Vec3d> circle_center_attrib(m.no_vertices(),Vec3d(0,0,0));
-//    for(auto vid = m.vertices_begin(); vid != m.vertices_end(); ++vid)
-//    {
-//        if(level_set_id_vertex[*vid]>-1)
-//        {
-//            LevelSetInfo& lsi = level_set_info[level_set_id_vertex[*vid]];
-//            circle_center_attrib[*vid] = lsi.avg_pos;
-//            nailed[*vid] = 1;
-//        }
-//    }
-////    propagate(m, nailed, circle_center_attrib);
-//    Manifold m2 = m;
-//    for(auto vid = m2.vertices_begin(); vid != m2.vertices_end(); ++vid)
-//    {
-//        m2.pos(*vid) = circle_pos[*vid];
-//        m2.pos(*vid)[2] = 0.01 * fun[*vid];
-//    }
-//    obj_save("blob.obj", m2);
-//    for(auto vid = m.vertices_begin(); vid != m.vertices_end(); ++vid)
-//    {
-//        Vec3d p =circle_pos[*vid];//-circle_center_attrib[*vid];
-//        parametrization[*vid] = Vec2d(p[0],p[1]);
-//    }
     // ----------------------------
     cout << "Remove vertices not on level set curves" << endl;
-    vector<VertexID> vid_vec;
-    for(VertexIDIterator vid = m.vertices_begin(); vid != m.vertices_end(); ++vid)
-        if(status[*vid]==0)
-            vid_vec.push_back(*vid);
-    random_shuffle(vid_vec.begin(), vid_vec.end());
-    for (size_t i=0; i<vid_vec.size(); ++i) {
-        FaceID f = m.merge_one_ring(vid_vec[i]);
-        if(f != InvalidFaceID)
-            shortest_edge_triangulate_face(m, f, level_set_id_vertex);
-        else
-            cout << "vertex not removed " << valency(m, vid_vec[i]) << endl;
-    }
-    for(FaceIDIterator fid = m.faces_begin(); fid != m.faces_end(); ++fid)
-        if(no_edges(m, *fid) > 3)
-            shortest_edge_triangulate_face(m, *fid, level_set_id_vertex);
-//    VertexAttributeVector<Vec3d> recalled_positions;
-//    for(VertexIDIterator vid = m.vertices_begin(); vid != m.vertices_end(); ++vid)
-//        recalled_positions[*vid] = m.pos(*vid);
-//
-//
-    TriangleQuality tq_energy(level_set_id_vertex);
-    priority_queue_optimization(m, tq_energy);
+}
 void make_height_fun(const HMesh::Manifold& m, HMesh::VertexAttributeVector<double>& fun,
                      double& vmin, double& vmax)
+{
     VertexIDIterator vid = m.vertices_begin();
-    vmin = vmax = m.pos(*vid)[2];
+    vmin = vmax = m.pos(*vid)[1];
     for(; vid != m.vertices_end(); ++vid)
+    {
         double v = dot(m.pos(*vid),Vec3d(0.0,1,0.00));
         fun[*vid] = v;
         vmin = min(v, vmin);
         vmax = max(v, vmax);
+    }
+}

Subversion Repositories gelsvn

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