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

 //
 //  polarize.cpp
 //  GEL
 //
 //  Created by J. Andreas Bærentzen on 18/03/12.
 //  Copyright 2012 __MyCompanyName__. All rights reserved.
 //
 #include <queue>
+#include <iomanip>
 #include "polarize.h"
 #include <CGLA/Vec2d.h>
 #include <LinAlg/LapackFunc.h>
 #include <HMesh/triangulate.h>
 #include <HMesh/obj_save.h>
 #include <HMesh/curvature.h>
 #include <HMesh/quadric_simplify.h>
 #include <HMesh/mesh_optimization.h>
 #include <HMesh/smooth.h>
 #include <complex>
 using namespace CGLA;
 using namespace std;
 using namespace HMesh;
-inline bool same_level(float a, float b) {return abs(a-b) < 0.00001;}
+inline bool same_level(double a, double b) {return abs(a-b) < 1e-10;}
 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;
+            }
+        }
+}
+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,
                       vector<vector<int>>& boundaries)
+{
     segmentation = FaceAttributeVector<int>(m.no_faces(), -1);
     int SEG_NO = 0;
     for(auto f = m.faces_begin(); f != m.faces_end(); ++f)
+    {
         if (segmentation[*f] == -1)
+        {
             queue<FaceID> q;
             q.push(*f);
             vector<int> bound(0);
             while (!q.empty())
+            {
                 FaceID face = q.front();
                 q.pop();
                 segmentation[face] = SEG_NO;
                 for(Walker w = m.walker(face); !w.full_circle(); w = w.next())
+                {
                     if(ls_id[w.halfedge()] == -1)
+                    {
                         FaceID fopp = w.opp().face();
                         if(segmentation[fopp] == -1) {
                             q.push(fopp);
+                        }
+                    }
                     else {
                         bound.push_back(ls_id[w.halfedge()]);
+                    }
+                }
+            }
             sort(begin(bound), end(bound));
             auto new_end = unique(begin(bound), end(bound));
             boundaries.push_back(vector<int>(begin(bound),new_end));
             cout << "Boundaries of seg# " << SEG_NO << " : ";
             for(auto be:boundaries[SEG_NO])
                 cout << be << " ";
             cout << endl;
             SEG_NO += 1;
+        }
+    }
+}
 void shortest_edge_triangulate_face(Manifold& m, FaceID f0, VertexAttributeVector<int>& level_set_id_vertex)
+{
     queue<FaceID> face_queue;
     face_queue.push(f0);
     while(!face_queue.empty())
+    {
         FaceID f = face_queue.front();
         face_queue.pop();
         // Create a vector of vertices.
         vector<VertexID> verts;
         for(Walker w = m.walker(f); !w.full_circle(); w = w.circulate_face_ccw())
+        {
             FaceID fa = w.face();
             FaceID fb = f;
             assert(fa==fb);
             verts.push_back(w.vertex());
+        }
         // If there are just three we are done.
         if(verts.size() == 3) continue;
         // Find vertex pairs that may be connected.
         vector<pair<int,int> > vpairs;
         const int N = verts.size();
         for(int i = 0; i < N - 2; ++i){
             for(int j = i + 2; j < N; ++j){
                 if(verts[i] != verts[j] &&
                    !connected(m, verts[i], verts[j]) &&
-                   (level_set_id_vertex[verts[i]] == 0 || level_set_id_vertex[verts[i]] != level_set_id_vertex[verts[j]])
+                   (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()){
             cout << "Warning: could not triangulate a face."
             << "Probably a vertex appears more than one time in other vertex's one-ring" << endl;
             continue;
+        }
         /* For all vertex pairs, find the edge lengths. Combine the
          vertices forming the shortest edge. */
         float min_len=FLT_MAX;
         int min_k = -1;
         for(size_t k = 0; k < vpairs.size(); ++k){
             int i = vpairs[k].first;
             int j = vpairs[k].second;
             float len = sqr_length(m.pos(verts[i]) - m.pos(verts[j]));
             if(len<min_len){
                 min_len = len;
                 min_k = k;
+            }
+        }
         assert(min_k != -1);
+        {
             // Split faces along edge whose midpoint is closest to isovalue
             int i = vpairs[min_k].first;
             int j = vpairs[min_k].second;
             FaceID f_new = m.split_face_by_edge(f, verts[i], verts[j]);
             if(no_edges(m, f)>3)
                 face_queue.push(f);
             if(no_edges(m, f_new)>3)
                 face_queue.push(f_new);
+        }
+    }
+}
 void shortest_edge_split_face(Manifold& m, FaceID f0, VertexAttributeVector<int>& level_set_id_vertex)
+{
     queue<FaceID> face_queue;
     face_queue.push(f0);
     while(!face_queue.empty())
+    {
         FaceID f = face_queue.front();
         face_queue.pop();
         // Create a vector of vertices.
         vector<VertexID> verts;
         for(Walker w = m.walker(f); !w.full_circle(); w = w.circulate_face_ccw())
+        {
             verts.push_back(w.vertex());
+        }
         // Find vertex pairs that may be connected.
         vector<pair<int,int> > vpairs;
         const int N = verts.size();
         for(int i = 0; i < N ; ++i){
             for(int j = 3; j < N-2; ++j){
                 int jj = (j+i)%N;
                 if(!connected(m, verts[i], verts[jj]) &&
                    (level_set_id_vertex[verts[i]] != level_set_id_vertex[verts[jj]]) &&
                    (level_set_id_vertex[verts[(i+1)%N]] == level_set_id_vertex[verts[i]]) &&
                    (level_set_id_vertex[verts[(i+N-1)%N]] == level_set_id_vertex[verts[i]]) &&
                    (level_set_id_vertex[verts[(jj+1)%N]] == level_set_id_vertex[verts[jj]]) &&
                    (level_set_id_vertex[verts[(jj+N-1)%N]] == level_set_id_vertex[verts[jj]]))
                     vpairs.push_back(pair<int,int>(i, jj));
+            }
+        }
         if(vpairs.empty()) continue;
         /* For all vertex pairs, find the edge lengths. Combine the
          vertices forming the shortest edge. */
         float min_len=FLT_MAX;
         int min_k = -1;
         for(size_t k = 0; k < vpairs.size(); ++k){
             int i = vpairs[k].first;
             int j = vpairs[k].second;
             float len = sqr_length(m.pos(verts[i]) - m.pos(verts[j]));
             if(len<min_len){
                 min_len = len;
                 min_k = k;
+            }
+        }
         assert(min_k != -1);
+        {
             // Split faces along edge whose midpoint is closest to isovalue
             int i = vpairs[min_k].first;
             int j = vpairs[min_k].second;
             FaceID f_new = m.split_face_by_edge(f, verts[i], verts[j]);
             if(no_edges(m, f)>5)
                 face_queue.push(f);
             if(no_edges(m, f_new)>5)
                 face_queue.push(f_new);
+        }
+    }
+}
 struct EdgeQElem {
     float priority;
     HalfEdgeID he;
     EdgeQElem(float _priority, HalfEdgeID _he): priority(_priority), he(_he) {}
 };
 bool operator<(const EdgeQElem& l, const EdgeQElem& r)
+{
     return l.priority < r.priority;
+}
 class FunctionalDifference: public EnergyFun
+{
     VertexAttributeVector<float>& fun;
     VertexAttributeVector<int>& status;
 public:
     FunctionalDifference(VertexAttributeVector<float>& _fun, VertexAttributeVector<int>& _status): fun(_fun), status(_status) {}
     virtual double delta_energy(const Manifold& m, HalfEdgeID h) const
+    {
         Walker w = m.walker(h);
         if(status[w.vertex()] == 1 && status[w.opp().vertex()]==1)
             return DBL_MAX;
         return sqr_length(m.pos(w.next().vertex())-m.pos(w.opp().next().vertex()))/(1e-6+abs(fun[w.next().vertex()]-fun[w.opp().next().vertex()])) - sqr_length(m.pos(w.vertex())-m.pos(w.opp().vertex()))/(1e-6+abs(fun[w.vertex()]-fun[w.opp().vertex()]));
+    }
 };
 class TriangleQualityValence: public EnergyFun
+{
     VertexAttributeVector<int>& idv;
     MinAngleEnergy mae;
     ValencyEnergy vae;
 public:
     TriangleQualityValence(VertexAttributeVector<int>& _idv): idv(_idv), mae(-1) {}
     virtual double delta_energy(const Manifold& m, HalfEdgeID h) const
+    {
         Walker w = m.walker(h);
         if(idv[w.next().vertex()] == idv[w.opp().next().vertex()] ||
            idv[w.vertex()] == idv[w.opp().vertex()])
             return DBL_MAX;
         VertexID v1 = w.opp().vertex();
         VertexID v2 = w.vertex();
         VertexID vo1 = w.next().vertex();
         VertexID vo2 = w.opp().next().vertex();
         int val1  = valency(m, v1);
         int val2  = valency(m, v2);
         int valo1 = valency(m, vo1);
         int valo2 = valency(m, vo2);
         // The optimal valency is four for a boundary vertex
         // and six elsewhere.
         int t1 = boundary(m, v1) ? 4 : 6;
         int t2 = boundary(m, v2) ? 4 : 6;
         int to1 = boundary(m, vo1) ? 4 : 6;
         int to2 = boundary(m, vo2) ? 4 : 6;
         int before =
         max(max(sqr(val1-t1),sqr(val2-t2)), max(sqr(valo1-to1), sqr(valo2-to2)));
         int after =
         max(max(sqr(valo1+1-to1),sqr(val1-1-t1)), max(sqr(val2-1-t2),sqr(valo2+1-to2)));
         return static_cast<double>(after-before);
         //        return vae.delta_energy(m,h);
         float la = length(m.pos(w.next().vertex())-m.pos(w.opp().next().vertex()));
         float lb = length(m.pos(w.vertex())-m.pos(w.opp().vertex()));
         return la-lb;
         return mae.delta_energy(m,h);
+    }
 };
 class TriangleQuality: public EnergyFun
+{
     VertexAttributeVector<int>& idv;
     MinAngleEnergy mae;
     ValencyEnergy vae;
 public:
     TriangleQuality(VertexAttributeVector<int>& _idv): idv(_idv), mae(-1) {}
     virtual double delta_energy(const Manifold& m, HalfEdgeID h) const
+    {
         Walker w = m.walker(h);
         if(idv[w.next().vertex()] == idv[w.opp().next().vertex()] ||
            idv[w.vertex()] == idv[w.opp().vertex()])
             return DBL_MAX;
         return mae.delta_energy(m,h);
+    }
 };
 Vec3d grad(HMesh::Manifold& m, HMesh::VertexAttributeVector<double>& fun, HMesh::FaceID f)
+{
     if(no_edges(m,f) != 3)
         return Vec3d(0.0);
     Vec3d n = normal(m, f);
     Vec3d gsum(0.0);
     for(Walker w = m.walker(f); !w.full_circle(); w = w.next())
+    {
         VertexID avid = w.opp().vertex();
         VertexID bvid = w.vertex();
         VertexID cvid = w.next().vertex();
         Vec3d gdir = normalize(cross(n, m.pos(bvid) - m.pos(avid)));
         double l = dot(m.pos(cvid)-m.pos(avid), gdir);
         gdir *= fun[cvid]/l;
         gsum += gdir;
+    }
     return gsum;
+}
+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]];
+    }
     Vec3d A(dot(uvw[0],uvw[0]),dot(uvw[0],uvw[1]),dot(uvw[0],uvw[2]));
     Vec3d B(dot(uvw[1],uvw[0]),dot(uvw[1],uvw[1]),dot(uvw[1],uvw[2]));
     Vec3d C(dot(uvw[2],uvw[0]),dot(uvw[2],uvw[1]),dot(uvw[2],uvw[2]));
     return - (beta * dot(B, abc) + gamma * dot(C, abc))/dot(A, abc);
+}
 void recompute_circle_param(Manifold& m, VertexID v, VertexAttributeVector<double>& fun, VertexAttributeVector<Vec2d>& circle_pos, double f_origin)
+{
     Vec2d c_new(0);
     for(auto w = m.walker(v); !w.full_circle(); w = w.circulate_vertex_ccw())
+    {
         Walker w_base = w.next();
         Vec2d X = circle_pos[w_base.opp().vertex()];
         Vec2d Y(-X[1], X[0]);
         Vec2d vec = circle_pos[w_base.vertex()];
         double angle = atan2(max(-1.0,min(1.0,dot(vec,Y))),
                              max(-1.0,min(1.0,dot(vec,X))));
         double new_angle = solve_for_orthogonal_gradients(m, fun, w_base.halfedge(), 0, angle);
         new_angle = max(-M_PI/4, min(M_PI/4,new_angle));
         c_new += cos(new_angle)*X + sin(new_angle)*Y;
+    }
     circle_pos[v] = normalize(c_new);
+}
+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);
+        }
+    }
     cout << "ids" << no_id << " " << level_set_info.size() << endl;
     for(int i=0;i<level_set_info.size(); ++i)
         for(int j=0;j<level_set_info.size(); ++j)
             if(i != j && same_level(level_set_info[i].fun_value, level_set_info[j].fun_value))
                 level_set_info[i].components += 1;
     HalfEdgeID min_length_h;
     double min_length_fun;
     double min_length=FLT_MAX;
     int min_length_id = -1;
     for(int i=0;i<level_set_info.size(); ++i)
+    {
         LevelSetInfo& lsi = level_set_info[i];
         if(lsi.components==1 && lsi.length<min_length)
+        {
             min_length_h = lsi.h;
             min_length_fun = lsi.fun_value;
             min_length = lsi.length;
             min_length_id = lsi.id;
+        }
+    }
-    smooth_fun(m, nailed, fun);
-    FaceAttributeVector<int> segmentation;
-    vector<vector<int>> boundaries;
-    segment_manifold(m, level_set_id, segmentation, boundaries);
-    VertexAttributeVector<Vec3d> cylinder_pos;
-    for(auto vid = m.vertices_begin();vid != m.vertices_end(); ++vid)
-        if (nailed[*vid])//level_set_id_vertex[*vid] == min_length_id)
-            cylinder_pos[*vid] = Vec3d(circle_pos[*vid][0], circle_pos[*vid][1],min_length_fun);
-        else
-            cylinder_pos[*vid] = Vec3d(0,0,fun[*vid]);
-    HalfEdgeAttributeVector<double> edge_weights(m.allocated_halfedges(), 0);
-    FaceAttributeVector<int> included(m.allocated_faces(), 1);
-    for(auto fid = m.faces_begin(); fid != m.faces_end(); ++fid)
-    {
-        if(boundaries[segmentation[*fid]].size() == 1)
-            included[*fid] = 0;
-    }
-    compute_edge_weights(m,edge_weights,included);
-    VertexAttributeVector<Vec3d> new_cylinder_pos(m.no_vertices());
-    for(int i = 0; i < 1500; ++i)
-    {
-        for(VertexIDIterator v = m.vertices_begin(); v != m.vertices_end(); ++v)
-            if(!nailed[*v])//level_set_id_vertex[*v] != min_length_id)
-            {
-//                recompute_circle_param(m, *v, fun, circle_pos, min_length_fun);
-                double w_sum = 0;
-                new_cylinder_pos[*v] = Vec3d(0);
-                for(Walker wv = m.walker(*v); !wv.full_circle(); wv = wv.circulate_vertex_ccw())
-                {
-                    double w = edge_weights[wv.halfedge()];
-                    new_cylinder_pos[*v] += w * cylinder_pos[wv.vertex()];
-                    w_sum += w;
-                }
-                if(w_sum> 0.0)
-                    new_cylinder_pos[*v] /= w_sum;
-                Vec2d v2d(new_cylinder_pos[*v][0], new_cylinder_pos[*v][1]);
-                double l =length(v2d);
-                if(l>0.1) {
-                    new_cylinder_pos[*v][0] /= l;
-                    new_cylinder_pos[*v][1] /= l;
-                }
-                if(level_set_id_vertex[*v] != -1)
-                    new_cylinder_pos[*v][2] = fun[*v];
-            }
-            else
-                new_cylinder_pos[*v] = cylinder_pos[*v];
-        for(auto v = m.vertices_begin();v != m.vertices_end(); ++v)
-            cylinder_pos[*v] = 0.9*cylinder_pos[*v]+ 0.1*new_cylinder_pos[*v];
-    }
-  for(auto vid = m.vertices_begin();vid != m.vertices_end(); ++vid)
-        parametrization[*vid] = Vec2d(cylinder_pos[*vid][0],cylinder_pos[*vid][1]);
-    return;
-//    for(auto vid = m.vertices_begin();vid != m.vertices_end(); ++vid)
-//        m.pos(*vid) = cylinder_pos[*vid];// Vec2d(cylinder_pos[*vid][0],cylinder_pos[*vid][1]);
-//    vector<pair<double,VertexID>> verts;
-//    for(auto vid = m.vertices_begin();vid != m.vertices_end(); ++vid)
-//        if(!same_level(min_length_fun, fun[*vid]))
-//            verts.push_back(pair<double,VertexID>(abs(fun[*vid]-min_length_fun),*vid));
-//    sort(verts.begin(), verts.end());
-//
-//    for(auto p:verts)
-//    {
-//        VertexID vid = p.second;
-//        recompute_circle_param(m, vid, fun, circle_pos, min_length_fun);
-//    }
-//    parametrization = circle_pos;
-//    return;
-//    // nail max length level set to circle.
-//    Walker w = m.walker(min_length_h);
-//    priority_queue<pair<double, HalfEdgeID>> hq;
-//    float param = 0;
-//    do
-//    {
-//        hq.push(pair<double,HalfEdgeID>(-abs(fun[w.next().vertex()]-min_length_fun),w.halfedge()));
-//        hq.push(pair<double,HalfEdgeID>(-abs(fun[w.opp().next().vertex()]-min_length_fun),w.opp().halfedge()));
-//
-//        nailed[w.opp().vertex()] = 1;
-//        w = w.next();
-//        while(level_set_id[w.halfedge()] != min_length_id)
-//            w = w.circulate_vertex_cw();
-//    }
-//    while(w.halfedge() != min_length_h);
-//    parametrize(m, nailed, circle_pos, parametrization, fun, min_length_fun);
-//    propagate_values(m, nailed, circle_pos, parametrization, level_set_id_vertex, level_set_info, fun, min_length_fun);
-////    for(int i=0;i<level_set_info.size(); ++i)
-////        level_set_info[i].print();
-//
-//    VertexAttributeVector<Vec3d> circle_center_attrib(m.no_vertices(),Vec3d(0,0,0));
-//    for(auto vid = m.vertices_begin(); vid != m.vertices_end(); ++vid)
-//    {
-//        if(level_set_id_vertex[*vid]>-1)
-//        {
-//            LevelSetInfo& lsi = level_set_info[level_set_id_vertex[*vid]];
-//            circle_center_attrib[*vid] = lsi.avg_pos;
-//            nailed[*vid] = 1;
-//        }
-//    }
-////    propagate(m, nailed, circle_center_attrib);
-//    Manifold m2 = m;
-//    for(auto vid = m2.vertices_begin(); vid != m2.vertices_end(); ++vid)
-//    {
-//        m2.pos(*vid) = circle_pos[*vid];
-//        m2.pos(*vid)[2] = 0.01 * fun[*vid];
-//    }
-//    obj_save("blob.obj", m2);
-//    for(auto vid = m.vertices_begin(); vid != m.vertices_end(); ++vid)
-//    {
-//        Vec3d p =circle_pos[*vid];//-circle_center_attrib[*vid];
-//        parametrization[*vid] = Vec2d(p[0],p[1]);
-//    }
     // ----------------------------
     cout << "Remove vertices not on level set curves" << endl;
-    vector<VertexID> vid_vec;
-    for(VertexIDIterator vid = m.vertices_begin(); vid != m.vertices_end(); ++vid)
-        if(status[*vid]==0)
-            vid_vec.push_back(*vid);
-    random_shuffle(vid_vec.begin(), vid_vec.end());
-    for (size_t i=0; i<vid_vec.size(); ++i) {
-        FaceID f = m.merge_one_ring(vid_vec[i]);
-        if(f != InvalidFaceID)
-            shortest_edge_triangulate_face(m, f, level_set_id_vertex);
-        else
-            cout << "vertex not removed " << valency(m, vid_vec[i]) << endl;
-    }
-    for(FaceIDIterator fid = m.faces_begin(); fid != m.faces_end(); ++fid)
-        if(no_edges(m, *fid) > 3)
-            shortest_edge_triangulate_face(m, *fid, level_set_id_vertex);
-//    VertexAttributeVector<Vec3d> recalled_positions;
-//    for(VertexIDIterator vid = m.vertices_begin(); vid != m.vertices_end(); ++vid)
-//        recalled_positions[*vid] = m.pos(*vid);
-//
-//
-    TriangleQuality tq_energy(level_set_id_vertex);
-    priority_queue_optimization(m, tq_energy);
+}
 void make_height_fun(const HMesh::Manifold& m, HMesh::VertexAttributeVector<double>& fun,
                      double& vmin, double& vmax)
+{
     VertexIDIterator vid = m.vertices_begin();
-    vmin = vmax = m.pos(*vid)[2];
+    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 @ 640 – Rev 631 → 640