WebSVN – gelsvn – Diff – /branches/GEL_github/src/HMesh/smooth.cpp

 /* ----------------------------------------------------------------------- *
  * This file is part of GEL, http://www.imm.dtu.dk/GEL
  * Copyright (C) the authors and DTU Informatics
  * For license and list of authors, see ../../doc/intro.pdf
  * ----------------------------------------------------------------------- */
 #include <thread>
 #include "smooth.h"
 #include <future>
 #include <vector>
 #include <algorithm>
 #include "../CGLA/Mat3x3d.h"
 #include "../CGLA/Vec3d.h"
 #include "../CGLA/Quatd.h"
 #include "../Util/Timer.h"
 #include "Manifold.h"
 #include "AttributeVector.h"
 namespace HMesh
 {
     using namespace std;
     using namespace CGLA;
     class range {
     public:
         class iterator {
             friend class range;
         public:
             long int operator *() const { return i_; }
             const iterator &operator ++() { ++i_; return *this; }
             iterator operator ++(int) { iterator copy(*this); ++i_; return copy; }
             bool operator ==(const iterator &other) const { return i_ == other.i_; }
             bool operator !=(const iterator &other) const { return i_ != other.i_; }
         protected:
             iterator(long int start) : i_ (start) { }
         private:
             unsigned long i_;
         };
         iterator begin() const { return begin_; }
         iterator end() const { return end_; }
         range(long int  begin, long int end) : begin_(begin), end_(end) {}
     private:
         iterator begin_;
         iterator end_;
     };
-    typedef std::vector<std::vector<VertexID>> VertexIDBatches;
+    void for_each_vertex(Manifold& m, function<void(VertexID)> f) { for(auto v : m.vertices()) f(v); }
+    void for_each_face(Manifold& m, function<void(FaceID)> f) { for(auto p : m.faces()) f(p); }
+    void for_each_halfedge(Manifold& m, function<void(HalfEdgeID)> f) { for(auto h : m.halfedges()) f(h); }
-    void for_each_vertex(Manifold& m, function<void(VertexID)> f) { for(auto v : m.vertices()) f(v); }
-    void for_each_face(Manifold& m, function<void(FaceID)> f) { for(auto p : m.faces()) f(p); }
+    inline Vec3d laplacian(const Manifold& m, VertexID v)
-    void for_each_halfedge(Manifold& m, function<void(HalfEdgeID)> f) { for(auto h : m.halfedges()) f(h); }
+    {
+        Vec3d p(0);
+        int n = circulate_vertex_ccw(m, v, (std::function<void(VertexID)>)[&](VertexID v){ p += m.pos(v); });
-    inline Vec3d laplacian(const Manifold& m, VertexID v)
+        return p / n - m.pos(v);
-    {
+    }
-        Vec3d p(0);
-        int n = circulate_vertex_ccw(m, v, (std::function<void(VertexID)>)[&](VertexID v){ p += m.pos(v); });
-        return p / n - m.pos(v);
+    int CORES = 8;
-    }
+    void laplacian_smooth_example(Manifold& m)
+    {
-    int CORES = 8;
+        VertexAttributeVector<Vec3d> new_pos = m.positions_attribute_vector();
+        for(VertexID v : m.vertices())
-    void laplacian_smooth_example(Manifold& m)
+            if(!boundary(m, v))
-    {
+            {
-        VertexAttributeVector<Vec3d> new_pos = m.positions_attribute_vector();
+                Vec3d L(0);
-        for(VertexID v : m.vertices())
+                int n = circulate_vertex_ccw(m, v, (std::function<void(VertexID)>)[&](VertexID v){ L += m.pos(v); });
-            if(!boundary(m, v))
+                new_pos[v] = L/n;
-            {
+            }
-                Vec3d L(0);
+        m.positions_attribute_vector() = new_pos;
-                int n = circulate_vertex_ccw(m, v, (std::function<void(VertexID)>)[&](VertexID v){ L += m.pos(v); });
-                new_pos[v] = L/n;
+    }
-            }
-        m.positions_attribute_vector() = new_pos;
+    void laplacian_smooth0(Manifold& m,float weight, int max_iter)
+    {
-    }
+        for(int iter=0;iter<max_iter; ++iter) {
+            VertexAttributeVector<Vec3d> L_attr(m.no_vertices());
-    void laplacian_smooth0(Manifold& m,float weight, int max_iter)
-    {
+            for(VertexIDIterator v = m.vertices_begin(); v != m.vertices_end(); ++v)
-        for(int iter=0;iter<max_iter; ++iter) {
+                if(!boundary(m, *v))
-            VertexAttributeVector<Vec3d> L_attr(m.no_vertices());
+                    L_attr[*v] =laplacian(m, *v);
+            for(VertexIDIterator v = m.vertices_begin(); v != m.vertices_end(); ++v){
-            for(VertexIDIterator v = m.vertices_begin(); v != m.vertices_end(); ++v)
+                if(!boundary(m, *v))
-                if(!boundary(m, *v))
+                    m.pos(*v) += weight*L_attr[*v];
-                    L_attr[*v] =laplacian(m, *v);
+            }
-            for(VertexIDIterator v = m.vertices_begin(); v != m.vertices_end(); ++v){
+        }
-                if(!boundary(m, *v))
+    }
-                    m.pos(*v) += weight*L_attr[*v];
-            }
+    void laplacian_smooth1(Manifold& m,float weight, int max_iter)
-        }
+    {
-    }
+        for(int iter=0;iter<max_iter; ++iter) {
+            VertexAttributeVector<Vec3d> L_attr(m.no_vertices());
-    void laplacian_smooth1(Manifold& m,float weight, int max_iter)
-    {
+            for(VertexID v : m.vertices())
-        for(int iter=0;iter<max_iter; ++iter) {
+                if(!boundary(m, v))
-            VertexAttributeVector<Vec3d> L_attr(m.no_vertices());
+                    L_attr[v] =laplacian(m, v);
+            for(VertexID v : m.vertices()){
-            for(VertexID v : m.vertices())
+                if(!boundary(m, v))
-                if(!boundary(m, v))
+                    m.pos(v) += weight*L_attr[v];
-                    L_attr[v] =laplacian(m, v);
+            }
-            for(VertexID v : m.vertices()){
+        }
-                if(!boundary(m, v))
+    }
-                    m.pos(v) += weight*L_attr[v];
-            }
+    void laplacian_smooth2(Manifold& m,float weight, int max_iter)
-        }
+    {
-    }
+        auto new_pos = m.positions_attribute_vector();
+        for(int iter=0;iter<max_iter; ++iter) {
-    void laplacian_smooth2(Manifold& m,float weight, int max_iter)
+            for(auto v : m.vertices())
-    {
+                if(!boundary(m, v))
-        auto new_pos = m.positions_attribute_vector();
+                    new_pos[v] = weight*laplacian(m, v)+m.pos(v);
-        for(int iter=0;iter<max_iter; ++iter) {
+            m.positions_attribute_vector() = new_pos;
-            for(auto v : m.vertices())
+        }
-                if(!boundary(m, v))
+    }
-                    new_pos[v] = weight*laplacian(m, v)+m.pos(v);
-            m.positions_attribute_vector() = new_pos;
-        }
-    }
+    void laplacian_smooth3(Manifold& m, float weight, int max_iter)
+    {
+        for(int iter=0;iter<max_iter; ++iter) {
+            auto new_pos = m.positions_attribute_vector();
-    void laplacian_smooth3(Manifold& m, float weight, int max_iter)
+            for(auto v : m.vertices())
-    {
+                if(!boundary(m, v))
-        for(int iter=0;iter<max_iter; ++iter) {
+                    new_pos[v] = weight*laplacian(m, v)+m.pos(v);
-            auto new_pos = m.positions_attribute_vector();
+            m.positions_attribute_vector() = move(new_pos);
-            for(auto v : m.vertices())
+        }
-                if(!boundary(m, v))
+    }
-                    new_pos[v] = weight*laplacian(m, v)+m.pos(v);
-            m.positions_attribute_vector() = move(new_pos);
+    void laplacian_smooth4(Manifold& m,float weight, int max_iter)
-        }
+    {
-    }
+        auto new_pos = m.positions_attribute_vector();
+        for(int iter=0;iter<max_iter; ++iter) {
-    void laplacian_smooth4(Manifold& m,float weight, int max_iter)
+            for(VertexID v : m.vertices())
-    {
+                if(!boundary(m, v))
-        auto new_pos = m.positions_attribute_vector();
+                    new_pos[v] = weight*laplacian(m, v)+m.pos(v);
-        for(int iter=0;iter<max_iter; ++iter) {
+            swap(m.positions_attribute_vector(),new_pos);
-            for(VertexID v : m.vertices())
+        }
-                if(!boundary(m, v))
+    }
-                    new_pos[v] = weight*laplacian(m, v)+m.pos(v);
-            swap(m.positions_attribute_vector(),new_pos);
+    void laplacian_smooth4_5(Manifold& m,float weight, int max_iter)
-        }
+    {
-    }
+        auto new_pos = m.positions_attribute_vector();
+        for(int iter=0;iter<max_iter; ++iter) {
-    void laplacian_smooth4_5(Manifold& m,float weight, int max_iter)
+            for_each_vertex(m, [&](VertexID v) {new_pos[v] = weight*laplacian(m, v)+m.pos(v);});
-    {
+            swap(m.positions_attribute_vector(),new_pos);
-        auto new_pos = m.positions_attribute_vector();
+        }
-        for(int iter=0;iter<max_iter; ++iter) {
+    }
-            for_each_vertex(m, [&](VertexID v) {new_pos[v] = weight*laplacian(m, v)+m.pos(v);});
-            swap(m.positions_attribute_vector(),new_pos);
-        }
-    }
+    void laplacian_smooth5(Manifold& m, float weight, int max_iter)
+    {
+        for(int iter=0;iter<max_iter; ++iter) {
+            auto new_pos = m.positions_attribute_vector();
-    void laplacian_smooth5(Manifold& m, float weight, int max_iter)
+            vector<thread> t_vec;
-    {
+            for(auto v : m.vertices())
-        for(int iter=0;iter<max_iter; ++iter) {
+                if(!boundary(m, v))
-            auto new_pos = m.positions_attribute_vector();
+                    t_vec.push_back(thread([&](VertexID vid){
-            vector<thread> t_vec;
+                        if(!boundary(m, vid))
-            for(auto v : m.vertices())
+                            new_pos[vid] = weight*laplacian(m, vid)+ m.pos(vid);},v));
-                if(!boundary(m, v))
+            for(int i=0;i<t_vec.size();++i)
-                    t_vec.push_back(thread([&](VertexID vid){
+                t_vec[i].join();
-                        if(!boundary(m, vid))
+            m.positions_attribute_vector() = move(new_pos);
-                            new_pos[vid] = weight*laplacian(m, vid)+ m.pos(vid);},v));
+        }
-            for(int i=0;i<t_vec.size();++i)
+    }
-                t_vec[i].join();
-            m.positions_attribute_vector() = move(new_pos);
-        }
-    }
+    inline void laplacian_smooth_vertex(Manifold& m, const vector<VertexID>& vids,
+                                        VertexAttributeVector<Vec3d>& new_pos,
+                                        float weight){
+        for(VertexID v: vids)
-    inline void laplacian_smooth_vertex(Manifold& m, const vector<VertexID>& vids,
+            new_pos[v] = m.pos(v)+weight*laplacian(m, v);
-                                        VertexAttributeVector<Vec3d>& new_pos,
+    }
-                                        float weight){
-        for(VertexID v: vids)
-            new_pos[v] = m.pos(v)+weight*laplacian(m, v);
+    void laplacian_smooth6(Manifold& m, float weight, int max_iter)
-    }
+    {
+        vector<vector<VertexID>> vertex_ids(CORES);
+        auto batch_size = m.no_vertices()/CORES;
-    void laplacian_smooth6(Manifold& m, float weight, int max_iter)
+        int cnt = 0;
-    {
+        for_each_vertex(m, [&](VertexID v) {
-        vector<vector<VertexID>> vertex_ids(CORES);
+            if (!boundary(m, v))
-        auto batch_size = m.no_vertices()/CORES;
+                vertex_ids[(cnt++/batch_size)%CORES].push_back(v);
-        int cnt = 0;
+        });
-        for_each_vertex(m, [&](VertexID v) {
+        vector<thread> t_vec(CORES);
-            if (!boundary(m, v))
+        VertexAttributeVector<Vec3d> new_pos = m.positions_attribute_vector();
-                vertex_ids[(cnt++/batch_size)%CORES].push_back(v);
+        for(int iter=0;iter<max_iter; ++iter) {
-        });
+            for(int thread_no=0;thread_no<CORES;++thread_no)
-        vector<thread> t_vec(CORES);
+                t_vec[thread_no] = thread(laplacian_smooth_vertex,
-        VertexAttributeVector<Vec3d> new_pos = m.positions_attribute_vector();
+                                          ref(m), ref(vertex_ids[thread_no]),
-        for(int iter=0;iter<max_iter; ++iter) {
+                                          ref(new_pos), weight);
             for(int thread_no=0;thread_no<CORES;++thread_no)
-                t_vec[thread_no] = thread(laplacian_smooth_vertex,
+                t_vec[thread_no].join();
-                                          ref(m), ref(vertex_ids[thread_no]),
+            swap(m.positions_attribute_vector(), new_pos);
-                                          ref(new_pos), weight);
+        }
-            for(int thread_no=0;thread_no<CORES;++thread_no)
-                t_vec[thread_no].join();
+    }
-            swap(m.positions_attribute_vector(), new_pos);
-        }
+    typedef std::vector<std::vector<VertexID>> VertexIDBatches;
-    }
+//    template<typename  T>
+//    void for_each_vertex_parallel(int no_threads, const VertexIDBatches& batches, const T& f) {
-//    template<typename  T>
+//        vector<future<void>> f_vec(no_threads);
-//    void for_each_vertex_parallel(int no_threads, const VertexIDBatches& batches, const T& f) {
+//        for(auto t : range(0, no_threads))
-//        vector<future<void>> f_vec(no_threads);
+//            f_vec[t] = async(launch::async, f, ref(batches[t]));
-//        for(auto t : range(0, no_threads))
+//    }
-//            f_vec[t] = async(launch::async, f, ref(batches[t]));
-//    }
     template<typename  T>
     void for_each_vertex_parallel(int no_threads, const VertexIDBatches& batches, const T& f) {
         vector<thread> t_vec(no_threads);
         for(auto t : range(0, no_threads))
             t_vec[t] = thread(f, ref(batches[t]));
         for(auto t : range(0, no_threads))
             t_vec[t].join();
     }
     VertexIDBatches batch_vertices(Manifold& m) {
         VertexIDBatches vertex_ids(CORES);
         auto batch_size = m.no_vertices()/CORES;
         int cnt = 0;
         for_each_vertex(m, [&](VertexID v) {
             //if (!boundary(m, v))
                 vertex_ids[(cnt++/batch_size)%CORES].push_back(v);
         });
         return vertex_ids;
     }
     void laplacian_smooth7(Manifold& m, float weight, int max_iter)
     {
         auto vertex_ids = batch_vertices(m);
         auto new_pos = m.positions_attribute_vector();
         auto f = [&](const vector<VertexID>& vids) {
             for(VertexID v: vids)
                 new_pos[v] = m.pos(v)+weight*laplacian(m, v);
         };
         for(auto _ : range(0, max_iter)) {
             for_each_vertex_parallel(CORES, vertex_ids, f);
             swap(m.positions_attribute_vector(), new_pos);
         }
     }
     void laplacian_smooth(Manifold& m,float weight, int max_iter)
     {
 //        laplacian_smooth_example(m);
         laplacian_smooth7(m,weight,max_iter);
 //        Util::Timer tim;
 //        cout << "Method 0: ";
 //        for(int i=0;i<5;++i) {
 //            Manifold m2 = m;
 //            tim.start();
 //            laplacian_smooth0(m2, 1, 1000);
 //            cout << tim.get_secs() << ", ";
 //        }
 //        cout << endl;
 //
 //        cout << "Method 1: ";
 //        for(int i=0;i<5;++i) {
 //            Manifold m2 = m;
 //            tim.start();
 //            laplacian_smooth1(m2, 1, 1000);
 //            cout << tim.get_secs() << ", ";
 //        }
 //        cout << endl;
 //
 //        cout << "Method 2: ";
 //        for(int i=0;i<5;++i) {
 //            Manifold m2 = m;
 //            tim.start();
 //            laplacian_smooth2(m2, 1, 1000);
 //            cout << tim.get_secs() << ", ";
 //        }
 //        cout << endl;
 //
 //        cout << "Method 3: ";
 //        for(int i=0;i<5;++i) {
 //            Manifold m2 = m;
 //            tim.start();
 //            laplacian_smooth3(m2, 1, 1000);
 //            cout << tim.get_secs() << ", ";
 //        }
 //        cout << endl;
 //
 //        cout << "Method 4: ";
 //        for(int i=0;i<5;++i) {
 //            Manifold m2 = m;
 //            tim.start();
 //            laplacian_smooth4(m2, 1, 1000);
 //            cout << tim.get_secs() << ", ";
 //        }
 //        cout << endl;
 //
 //        cout << "Method 5: ";
 //        for(int i=0;i<3;++i) {
 //            Manifold m2 = m;
 //            tim.start();
 //            laplacian_smooth5(m2, 1, 1000);
 //            cout << " sec: "<< tim.get_secs() << ",";
 //        }
 //        cout << endl;
 //
 //        CORES = 2;
 //        cout << "Method 6 (2): ";
 //        for(int i=0;i<5;++i) {
 //            Manifold m2 = m;
 //            tim.start();
 //            laplacian_smooth6(m2, 1, 1000);
 //            cout << tim.get_secs() << ", ";
 //        }
 //        cout << endl;
 //
 //        CORES = 4;
 //        cout << "Method 7 (4): ";
 //        for(int i=0;i<5;++i) {
 //            Manifold m2 = m;
 //            tim.start();
 //            laplacian_smooth6(m2, 1, 1000);
 //            cout << tim.get_secs() << ", ";
 //        }
 //        cout << endl;
 //
 //        CORES = 8;
 //        cout << "Method 7 (8): ";
 //        for(int i=0;i<5;++i) {
 //            Manifold m2 = m;
 //            tim.start();
 //            laplacian_smooth6(m2, 1, 1000);
 //            cout << tim.get_secs() << ", ";
 //        }
 //        cout << endl;
     }
     void taubin_smooth(Manifold& m, int max_iter)
     {
         auto new_pos = m.positions_attribute_vector();
         for(int iter = 0; iter < 2*max_iter; ++iter) {
             for(VertexID v : m.vertices())
                 if(!boundary(m, v))
                     new_pos[v] = (iter%2 == 0 ? 0.5 : -0.52) * laplacian(m, v) + m.pos(v);
             swap(m.positions_attribute_vector(), new_pos);
         }
     }
     void face_neighbourhood(Manifold& m, FaceID f, vector<FaceID>& nbrs)
     {
         FaceAttributeVector<int> touched(m.allocated_faces(), -1);
         touched[f] = 1;
         nbrs.push_back(f);
         for(Walker wf = m.walker(f); !wf.full_circle(); wf = wf.circulate_face_cw()){
             for(Walker wv = m.walker(wf.vertex()); !wv.full_circle(); wv = wv.circulate_vertex_cw()){
                 FaceID fn = wv.face();
                 if(fn != InvalidFaceID && touched[fn] != touched[f]){
                     nbrs.push_back(fn);
                     touched[fn] = 1;
                 }
             }
         }
     }
     Vec3d fvm_filtered_normal(Manifold& m, FaceID f)
     {
         const float sigma = .1f;
         vector<FaceID> nbrs;
         face_neighbourhood(m, f, nbrs);
         float min_dist_sum=1e32f;
         long int median=-1;
         vector<Vec3d> normals(nbrs.size());
         for(size_t i=0;i<nbrs.size();++i)
         {
             normals[i] = normal(m, nbrs[i]);
             float dist_sum = 0;
             for(size_t j=0;j<nbrs.size(); ++j)
                 dist_sum += 1.0f - dot(normals[i], normals[j]);
             if(dist_sum < min_dist_sum)
             {
                 min_dist_sum = dist_sum;
                 median = i;
             }
         }
         assert(median != -1);
         Vec3d median_norm = normals[median];
         Vec3d avg_norm(0);
         for(size_t i=0;i<nbrs.size();++i)
         {
             float w = exp((dot(median_norm, normals[i])-1)/sigma);
             if(w<1e-2) w = 0;
             avg_norm += w*normals[i];
         }
         return normalize(avg_norm);
     }
     Vec3d bilateral_filtered_normal(Manifold& m, FaceID f, double avg_len)
     {
         vector<FaceID> nbrs;
         face_neighbourhood(m, f, nbrs);
         Vec3d p0 = centre(m, f);
         Vec3d n0 = normal(m, f);
         vector<Vec3d> normals(nbrs.size());
         vector<Vec3d> pos(nbrs.size());
         Vec3d fn(0);
         for(FaceID nbr : nbrs)
         {
             Vec3d n = normal(m, nbr);
             Vec3d p = centre(m, nbr);
             double w_a = exp(-acos(max(-1.0,min(1.0,dot(n,n0))))/(M_PI/32.0));
             double w_s = exp(-length(p-p0)/avg_len);
             fn += area(m, nbr)* w_a * w_s * n;
         }
         return normalize(fn);
     }
     void anisotropic_smooth(HMesh::Manifold& m, int max_iter, NormalSmoothMethod nsm)
     {
         double avg_len=0;
         for(HalfEdgeID hid : m.halfedges())
             avg_len += length(m, hid);
         avg_len /= 2.0;
         for(int iter = 0;iter<max_iter; ++iter)
         {
             FaceAttributeVector<Vec3d> filtered_norms(m.allocated_faces());
             for(FaceID f: m.faces()){
                 filtered_norms[f] = (nsm == BILATERAL_NORMAL_SMOOTH)?
                 bilateral_filtered_normal(m, f, avg_len):
                 fvm_filtered_normal(m, f);
             }
             VertexAttributeVector<Vec3d> vertex_positions(m.allocated_vertices(), Vec3d(0));
             VertexAttributeVector<int> count(m.allocated_vertices(), 0);
             for(int sub_iter=0;sub_iter<100;++sub_iter)
             {
                 for(HalfEdgeID hid : m.halfedges())
                 {
                     Walker w = m.walker(hid);
                     FaceID f = w.face();
                     if(f != InvalidFaceID){
                         VertexID v = w.vertex();
                         Vec3d dir = m.pos(w.opp().vertex()) - m.pos(v);
                         Vec3d n = filtered_norms[f];
                         vertex_positions[v] += m.pos(v) + 0.5 * n * dot(n, dir);
                         count[v] += 1;
                     }
                 }
                 double max_move= 0;
                 for(VertexID v : m.vertices())
                 {
                     Vec3d npos = vertex_positions[v] / double(count[v]);
                     double move = sqr_length(npos - m.pos(v));
                     if(move > max_move)
                         max_move = move;
                     m.pos(v) = npos;
                 }
                 if(max_move<sqr(1e-8*avg_len))
                 {
                     cout << "iters " << sub_iter << endl;
                     break;
                 }
             }
         }
     }
     void TAL_smoothing(Manifold& m, float w, int max_iter)
     {
         for(int iter=0;iter<max_iter;++iter) {
             VertexAttributeVector<float> vertex_areas;
             VertexAttributeVector<Vec3d> laplacians;
             for(VertexIDIterator vid = m.vertices_begin(); vid != m.vertices_end(); ++vid)
             {
                 vertex_areas[*vid] = 0;
                 for(Walker w = m.walker(*vid); !w.full_circle(); w = w.circulate_vertex_ccw())
                     if(w.face() != InvalidFaceID)
                         vertex_areas[*vid] += area(m, w.face());
             }
             for(VertexIDIterator vid = m.vertices_begin(); vid != m.vertices_end(); ++vid)
             {
                 laplacians[*vid] = Vec3d(0);
                 double weight_sum = 0.0;
                 if(boundary(m, *vid))
                 {
                     double angle_sum = 0;
                     for(Walker w = m.walker(*vid); !w.full_circle(); w = w.circulate_vertex_ccw())
                     {
                         if (w.face() != InvalidFaceID)
                         {
                             Vec3d vec_a = normalize(m.pos(w.vertex()) - m.pos(*vid));
                             Vec3d vec_b = normalize(m.pos(w.circulate_vertex_ccw().vertex()) -
                                                     m.pos(*vid));
                             angle_sum += acos(max(-1.0,min(1.0,dot(vec_a,vec_b))));
                         }
                         if(boundary(m,w.vertex()))
                         {
                             laplacians[*vid] += m.pos(w.vertex()) - m.pos(*vid);
                             weight_sum += 1.0;
                         }
                     }
                     laplacians[*vid] /= weight_sum;
                     laplacians[*vid] *= exp(-3.0*sqr(max(0.0, M_PI-angle_sum)));
                 }
                 else
                 {
                     for(Walker w = m.walker(*vid); !w.full_circle(); w = w.circulate_vertex_ccw())
                     {
                         float weight = vertex_areas[w.vertex()];
                         Vec3d l = m.pos(w.vertex()) - m.pos(*vid);
                         laplacians[*vid] +=  weight * l;
                         weight_sum += weight;
                     }
                     laplacians[*vid] /= weight_sum;
                     //                Vec3d n = normal(m, *vid);
                     //                if(sqr_length(n)>0.9)
                     //                    laplacians[*vid] -= n * dot(n, laplacians[*vid]);
                 }
             }
             for(VertexIDIterator vid = m.vertices_begin(); vid != m.vertices_end(); ++vid)
                 m.pos(*vid) += w*laplacians[*vid];
         }
     }
 }

Subversion Repositories gelsvn

(root)/branches/GEL_github/src/HMesh/smooth.cpp @ 684 – Rev 667 → 677