WebSVN – gelsvn – Diff – /branches/cpp11-devel/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 <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_;
    };
 
 
 
    
    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); }
 
    
    inline Vec3d laplacian(const Manifold& m, VertexID v)
    {
        Vec3d p(0);
 
 
        int n = circulate_vertex_ccw(m, v, [&](VertexID v){ p += m.pos(v); });
 
 
 
        return p / n - m.pos(v);
    }
    
    
    int CORES = 8;
    
    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());
            
            for(VertexIDIterator v = m.vertices_begin(); v != m.vertices_end(); ++v)
                if(!boundary(m, *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());
            
            for(VertexID v : m.vertices())
                if(!boundary(m, 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) {
            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();
            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) {
            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) {
            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();
            vector<thread> t_vec;
            for(auto v : m.vertices())
                if(!boundary(m, v))
                    t_vec.push_back(thread([&](VertexID vid){
                        if(!boundary(m, vid))
                            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)
            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;
        int cnt = 0;
        for_each_vertex(m, [&](VertexID v) {
            if (!boundary(m, v))
                vertex_ids[(cnt++/batch_size)%CORES].push_back(v);
        });
        vector<thread> t_vec(CORES);
        VertexAttributeVector<Vec3d> new_pos = m.positions_attribute_vector();
        for(int iter=0;iter<max_iter; ++iter) {
            for(int thread_no=0;thread_no<CORES;++thread_no)
                t_vec[thread_no] = thread(laplacian_smooth_vertex,
                                          ref(m), ref(vertex_ids[thread_no]),
                                          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 vector<vector<VertexID>> VertexIDBatches;
    
    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_smooth6(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 6 (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 6 (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;
                }
            }
        }
    }
    
    
}
 

Rev 632	Rev 633
1	`/* ----------------------------------------------------------------------- *`	1	`/* ----------------------------------------------------------------------- *`
2	`* This file is part of GEL, http://www.imm.dtu.dk/GEL`	2	`* This file is part of GEL, http://www.imm.dtu.dk/GEL`
3	`* Copyright (C) the authors and DTU Informatics`	3	`* Copyright (C) the authors and DTU Informatics`
4	`* For license and list of authors, see ../../doc/intro.pdf`	4	`* For license and list of authors, see ../../doc/intro.pdf`
5	`* ----------------------------------------------------------------------- */`	5	`* ----------------------------------------------------------------------- */`
6		6
7	`#include <thread>`	7	`#include <thread>`
8	`#include "smooth.h"`	8	`#include "smooth.h"`
9		9
10	`#include <vector>`	10	`#include <vector>`
11	`#include <algorithm>`	11	`#include <algorithm>`
12	`#include "../CGLA/Mat3x3d.h"`	12	`#include "../CGLA/Mat3x3d.h"`
13	`#include "../CGLA/Vec3d.h"`	13	`#include "../CGLA/Vec3d.h"`
14	`#include "../CGLA/Quatd.h"`	14	`#include "../CGLA/Quatd.h"`
15	`#include "../Util/Timer.h"`	15	`#include "../Util/Timer.h"`
16		16
17	`#include "Manifold.h"`	17	`#include "Manifold.h"`
18	`#include "AttributeVector.h"`	18	`#include "AttributeVector.h"`
19		19
20	`namespace HMesh`	20	`namespace HMesh`
21	`{`	21	`{`
22	`using namespace std;`	22	`using namespace std;`
23	`using namespace CGLA;`	23	`using namespace CGLA;`
24		24
-		25
-		26	`class range {`
-		27	`public:`
-		28	`class iterator {`
-		29	`friend class range;`
-		30	`public:`
-		31	`long int operator *() const { return i_; }`
-		32	`const iterator &operator ++() { ++i_; return *this; }`
-		33	`iterator operator ++(int) { iterator copy(*this); ++i_; return copy; }`
-		34
-		35	`bool operator ==(const iterator &other) const { return i_ == other.i_; }`
-		36	`bool operator !=(const iterator &other) const { return i_ != other.i_; }`
-		37
-		38	`protected:`
-		39	`iterator(long int start) : i_ (start) { }`
-		40
-		41	`private:`
-		42	`unsigned long i_;`
-		43	`};`
-		44
-		45	`iterator begin() const { return begin_; }`
-		46	`iterator end() const { return end_; }`
-		47	`range(long int begin, long int end) : begin_(begin), end_(end) {}`
-		48	`private:`
-		49	`iterator begin_;`
-		50	`iterator end_;`
-		51	`};`
-		52
25		53
26		54
27		55
28	`void for_each_vertex(Manifold& m, function<void(VertexID)> f) { for(auto v : m.vertices()) f(v); }`	56	`void for_each_vertex(Manifold& m, function<void(VertexID)> f) { for(auto v : m.vertices()) f(v); }`
29	`void for_each_face(Manifold& m, function<void(FaceID)> f) { for(auto p : m.faces()) f(p); }`	57	`void for_each_face(Manifold& m, function<void(FaceID)> f) { for(auto p : m.faces()) f(p); }`
30	`void for_each_halfedge(Manifold& m, function<void(HalfEdgeID)> f) { for(auto h : m.halfedges()) f(h); }`	58	`void for_each_halfedge(Manifold& m, function<void(HalfEdgeID)> f) { for(auto h : m.halfedges()) f(h); }`
31		59
32		60
33	`Vec3d laplacian(const Manifold& m, VertexID v)`	61	`inline Vec3d laplacian(const Manifold& m, VertexID v)`
34	`{`	62	`{`
35	`Vec3d avg_pos(0);`	63	`Vec3d p(0);`
36	`int n = 0;`	-
37		-
38	`circulate_vertex_ccw(m, v, [&](VertexID v){`	64	`int n = circulate_vertex_ccw(m, v, [&](VertexID v){ p += m.pos(v); });`
39	`avg_pos += m.pos(v);`	-
40	`++n;`	-
41	`});`	-
42	`return avg_pos / n - m.pos(v);`	65	`return p / n - m.pos(v);`
43	`}`	66	`}`
44		67
45		68
46	`int CORES = 8;`	69	`int CORES = 8;`
47		70
48	`void laplacian_smooth0(Manifold& m,float weight, int max_iter)`	71	`void laplacian_smooth0(Manifold& m,float weight, int max_iter)`
49	`{`	72	`{`
50	`for(int iter=0;iter<max_iter; ++iter) {`	73	`for(int iter=0;iter<max_iter; ++iter) {`
51	`VertexAttributeVector<Vec3d> L_attr(m.no_vertices());`	74	`VertexAttributeVector<Vec3d> L_attr(m.no_vertices());`
52		75
53	`for(VertexIDIterator v = m.vertices_begin(); v != m.vertices_end(); ++v)`	76	`for(VertexIDIterator v = m.vertices_begin(); v != m.vertices_end(); ++v)`
54	`if(!boundary(m, *v))`	77	`if(!boundary(m, *v))`
55	`L_attr[v] =laplacian(m, v);`	78	`L_attr[v] =laplacian(m, v);`
56	`for(VertexIDIterator v = m.vertices_begin(); v != m.vertices_end(); ++v){`	79	`for(VertexIDIterator v = m.vertices_begin(); v != m.vertices_end(); ++v){`
57	`if(!boundary(m, *v))`	80	`if(!boundary(m, *v))`
58	`m.pos(v) += weightL_attr[*v];`	81	`m.pos(v) += weightL_attr[*v];`
59	`}`	82	`}`
60	`}`	83	`}`
61	`}`	84	`}`
62		85
63	`void laplacian_smooth1(Manifold& m,float weight, int max_iter)`	86	`void laplacian_smooth1(Manifold& m,float weight, int max_iter)`
64	`{`	87	`{`
65	`for(int iter=0;iter<max_iter; ++iter) {`	88	`for(int iter=0;iter<max_iter; ++iter) {`
66	`VertexAttributeVector<Vec3d> L_attr(m.no_vertices());`	89	`VertexAttributeVector<Vec3d> L_attr(m.no_vertices());`
67		90
68	`for(VertexID v : m.vertices())`	91	`for(VertexID v : m.vertices())`
69	`if(!boundary(m, v))`	92	`if(!boundary(m, v))`
70	`L_attr[v] =laplacian(m, v);`	93	`L_attr[v] =laplacian(m, v);`
71	`for(VertexID v : m.vertices()){`	94	`for(VertexID v : m.vertices()){`
72	`if(!boundary(m, v))`	95	`if(!boundary(m, v))`
73	`m.pos(v) += weight*L_attr[v];`	96	`m.pos(v) += weight*L_attr[v];`
74	`}`	97	`}`
75	`}`	98	`}`
76	`}`	99	`}`
77		100
78	`void laplacian_smooth2(Manifold& m,float weight, int max_iter)`	101	`void laplacian_smooth2(Manifold& m,float weight, int max_iter)`
79	`{`	102	`{`
80	`auto new_pos = m.positions_attribute_vector();`	103	`auto new_pos = m.positions_attribute_vector();`
81	`for(int iter=0;iter<max_iter; ++iter) {`	104	`for(int iter=0;iter<max_iter; ++iter) {`
82	`for(auto v : m.vertices())`	105	`for(auto v : m.vertices())`
83	`if(!boundary(m, v))`	106	`if(!boundary(m, v))`
84	`new_pos[v] = weight*laplacian(m, v)+m.pos(v);`	107	`new_pos[v] = weight*laplacian(m, v)+m.pos(v);`
85	`m.positions_attribute_vector() = new_pos;`	108	`m.positions_attribute_vector() = new_pos;`
86	`}`	109	`}`
87	`}`	110	`}`
88		111
89		112
90		113
91	`void laplacian_smooth3(Manifold& m, float weight, int max_iter)`	114	`void laplacian_smooth3(Manifold& m, float weight, int max_iter)`
92	`{`	115	`{`
93	`for(int iter=0;iter<max_iter; ++iter) {`	116	`for(int iter=0;iter<max_iter; ++iter) {`
94	`auto new_pos = m.positions_attribute_vector();`	117	`auto new_pos = m.positions_attribute_vector();`
95	`for(auto v : m.vertices())`	118	`for(auto v : m.vertices())`
96	`if(!boundary(m, v))`	119	`if(!boundary(m, v))`
97	`new_pos[v] = weight*laplacian(m, v)+m.pos(v);`	120	`new_pos[v] = weight*laplacian(m, v)+m.pos(v);`
98	`m.positions_attribute_vector() = move(new_pos);`	121	`m.positions_attribute_vector() = move(new_pos);`
99	`}`	122	`}`
100	`}`	123	`}`
101		124
102	`void laplacian_smooth4(Manifold& m,float weight, int max_iter)`	125	`void laplacian_smooth4(Manifold& m,float weight, int max_iter)`
103	`{`	126	`{`
104	`auto new_pos = m.positions_attribute_vector();`	127	`auto new_pos = m.positions_attribute_vector();`
105	`for(int iter=0;iter<max_iter; ++iter) {`	128	`for(int iter=0;iter<max_iter; ++iter) {`
106	`for(VertexID v : m.vertices())`	129	`for(VertexID v : m.vertices())`
107	`if(!boundary(m, v))`	130	`if(!boundary(m, v))`
108	`new_pos[v] = weight*laplacian(m, v)+m.pos(v);`	131	`new_pos[v] = weight*laplacian(m, v)+m.pos(v);`
109	`swap(m.positions_attribute_vector(),new_pos);`	132	`swap(m.positions_attribute_vector(),new_pos);`
110	`}`	133	`}`
111	`}`	134	`}`
112		135
113	`void laplacian_smooth4_5(Manifold& m,float weight, int max_iter)`	136	`void laplacian_smooth4_5(Manifold& m,float weight, int max_iter)`
114	`{`	137	`{`
115	`auto new_pos = m.positions_attribute_vector();`	138	`auto new_pos = m.positions_attribute_vector();`
116	`for(int iter=0;iter<max_iter; ++iter) {`	139	`for(int iter=0;iter<max_iter; ++iter) {`
117	`for_each_vertex(m, [&](VertexID v) {new_pos[v] = weight*laplacian(m, v)+m.pos(v);});`	140	`for_each_vertex(m, [&](VertexID v) {new_pos[v] = weight*laplacian(m, v)+m.pos(v);});`
118	`swap(m.positions_attribute_vector(),new_pos);`	141	`swap(m.positions_attribute_vector(),new_pos);`
119	`}`	142	`}`
120	`}`	143	`}`
121		144
122		145
123		146
124	`void laplacian_smooth5(Manifold& m, float weight, int max_iter)`	147	`void laplacian_smooth5(Manifold& m, float weight, int max_iter)`
125	`{`	148	`{`
126	`for(int iter=0;iter<max_iter; ++iter) {`	149	`for(int iter=0;iter<max_iter; ++iter) {`
127	`auto new_pos = m.positions_attribute_vector();`	150	`auto new_pos = m.positions_attribute_vector();`
128	`vector<thread> t_vec;`	151	`vector<thread> t_vec;`
129	`for(auto v : m.vertices())`	152	`for(auto v : m.vertices())`
130	`if(!boundary(m, v))`	153	`if(!boundary(m, v))`
131	`t_vec.push_back(thread([&](VertexID vid){`	154	`t_vec.push_back(thread([&](VertexID vid){`
132	`if(!boundary(m, vid))`	155	`if(!boundary(m, vid))`
133	`new_pos[vid] = weight*laplacian(m, vid)+ m.pos(vid);},v));`	156	`new_pos[vid] = weight*laplacian(m, vid)+ m.pos(vid);},v));`
134	`for(int i=0;i<t_vec.size();++i)`	157	`for(int i=0;i<t_vec.size();++i)`
135	`t_vec[i].join();`	158	`t_vec[i].join();`
136	`m.positions_attribute_vector() = move(new_pos);`	159	`m.positions_attribute_vector() = move(new_pos);`
137	`}`	160	`}`
138	`}`	161	`}`
139		162
140		163
141		164
142	`inline void laplacian_smooth_vertex(Manifold& m,vector<VertexID>& vids,`	165	`inline void laplacian_smooth_vertex(Manifold& m, const vector<VertexID>& vids,`
143	`VertexAttributeVector<Vec3d>& new_pos,`	166	`VertexAttributeVector<Vec3d>& new_pos,`
144	`float weight){`	167	`float weight){`
145	`for(VertexID v: vids)`	168	`for(VertexID v: vids)`
146	`new_pos[v] = m.pos(v)+weight*laplacian(m, v);`	169	`new_pos[v] = m.pos(v)+weight*laplacian(m, v);`
147	`}`	170	`}`
148		171
-		172
149	`void laplacian_smooth6(Manifold& m, float weight, int max_iter)`	173	`void laplacian_smooth6(Manifold& m, float weight, int max_iter)`
150	`{`	174	`{`
151	`vector<vector<VertexID>> vertex_ids(CORES);`	175	`vector<vector<VertexID>> vertex_ids(CORES);`
152	`auto batch_size = m.no_vertices()/CORES;`	176	`auto batch_size = m.no_vertices()/CORES;`
153	`int cnt = 0;`	177	`int cnt = 0;`
154	`for_each_vertex(m, [&](VertexID v) {`	178	`for_each_vertex(m, [&](VertexID v) {`
155	`if (!boundary(m, v))`	179	`if (!boundary(m, v))`
156	`vertex_ids[(cnt++/batch_size)%CORES].push_back(v);`	180	`vertex_ids[(cnt++/batch_size)%CORES].push_back(v);`
157	`});`	181	`});`
158	`vector<thread> t_vec(CORES);`	182	`vector<thread> t_vec(CORES);`
159	`VertexAttributeVector<Vec3d> new_pos = m.positions_attribute_vector();`	183	`VertexAttributeVector<Vec3d> new_pos = m.positions_attribute_vector();`
160	`for(int iter=0;iter<max_iter; ++iter) {`	184	`for(int iter=0;iter<max_iter; ++iter) {`
161	`for(int thread_no=0;thread_no<CORES;++thread_no)`	185	`for(int thread_no=0;thread_no<CORES;++thread_no)`
162	`t_vec[thread_no] = thread(laplacian_smooth_vertex,`	186	`t_vec[thread_no] = thread(laplacian_smooth_vertex,`
163	`ref(m), ref(vertex_ids[thread_no]),`	187	`ref(m), ref(vertex_ids[thread_no]),`
164	`ref(new_pos), weight);`	188	`ref(new_pos), weight);`
165	`for(int thread_no=0;thread_no<CORES;++thread_no)`	189	`for(int thread_no=0;thread_no<CORES;++thread_no)`
166	`t_vec[thread_no].join();`	190	`t_vec[thread_no].join();`
167	`swap(m.positions_attribute_vector(), new_pos);`	191	`swap(m.positions_attribute_vector(), new_pos);`
168	`}`	192	`}`
169		193
170	`}`	194	`}`
171		195
-		196	`typedef vector<vector<VertexID>> VertexIDBatches;`
-		197
-		198	`template<typename T>`
-		199	`void for_each_vertex_parallel(int no_threads, const VertexIDBatches& batches, const T& f) {`
-		200	`vector<thread> t_vec(no_threads);`
-		201	`for(auto t : range(0, no_threads))`
-		202	`t_vec[t] = thread(f, ref(batches[t]));`
-		203	`for(auto t : range(0, no_threads))`
-		204	`t_vec[t].join();`
-		205	`}`
-		206
-		207	`VertexIDBatches batch_vertices(Manifold& m) {`
-		208	`VertexIDBatches vertex_ids(CORES);`
-		209	`auto batch_size = m.no_vertices()/CORES;`
-		210	`int cnt = 0;`
-		211	`for_each_vertex(m, [&](VertexID v) {`
-		212	`if (!boundary(m, v))`
-		213	`vertex_ids[(cnt++/batch_size)%CORES].push_back(v);`
-		214	`});`
-		215	`return vertex_ids;`
-		216	`}`
-		217
-		218	`void laplacian_smooth7(Manifold& m, float weight, int max_iter)`
-		219	`{`
-		220	`auto vertex_ids = batch_vertices(m);`
-		221	`auto new_pos = m.positions_attribute_vector();`
-		222	`auto f = [&](const vector<VertexID>& vids) {`
-		223	`for(VertexID v: vids)`
-		224	`new_pos[v] = m.pos(v)+weight*laplacian(m, v);`
-		225	`};`
-		226
-		227	`for(auto _ : range(0, max_iter)) {`
-		228	`for_each_vertex_parallel(CORES, vertex_ids, f);`
-		229	`swap(m.positions_attribute_vector(), new_pos);`
-		230	`}`
-		231	`}`
-		232
-		233
-		234
-		235
172	`void laplacian_smooth(Manifold& m,float weight, int max_iter)`	236	`void laplacian_smooth(Manifold& m,float weight, int max_iter)`
173	`{`	237	`{`
174	`laplacian_smooth4(m,weight,max_iter);`	238	`laplacian_smooth6(m,weight,max_iter);`
175	`// Util::Timer tim;`	239	`// Util::Timer tim;`
176	`//`	240	`//`
177	`// cout << "Method 0: ";`	241	`// cout << "Method 0: ";`
178	`// for(int i=0;i<5;++i) {`	242	`// for(int i=0;i<5;++i) {`
179	`// Manifold m2 = m;`	243	`// Manifold m2 = m;`
180	`// tim.start();`	244	`// tim.start();`
181	`// laplacian_smooth0(m2, 1, 1000);`	245	`// laplacian_smooth0(m2, 1, 1000);`
182	`// cout << tim.get_secs() << ", ";`	246	`// cout << tim.get_secs() << ", ";`
183	`// }`	247	`// }`
184	`// cout << endl;`	248	`// cout << endl;`
185	`//`	249	`//`
186	`// cout << "Method 1: ";`	250	`// cout << "Method 1: ";`
187	`// for(int i=0;i<5;++i) {`	251	`// for(int i=0;i<5;++i) {`
188	`// Manifold m2 = m;`	252	`// Manifold m2 = m;`
189	`// tim.start();`	253	`// tim.start();`
190	`// laplacian_smooth1(m2, 1, 1000);`	254	`// laplacian_smooth1(m2, 1, 1000);`
191	`// cout << tim.get_secs() << ", ";`	255	`// cout << tim.get_secs() << ", ";`
192	`// }`	256	`// }`
193	`// cout << endl;`	257	`// cout << endl;`
194	`//`	258	`//`
195	`// cout << "Method 2: ";`	259	`// cout << "Method 2: ";`
196	`// for(int i=0;i<5;++i) {`	260	`// for(int i=0;i<5;++i) {`
197	`// Manifold m2 = m;`	261	`// Manifold m2 = m;`
198	`// tim.start();`	262	`// tim.start();`
199	`// laplacian_smooth2(m2, 1, 1000);`	263	`// laplacian_smooth2(m2, 1, 1000);`
200	`// cout << tim.get_secs() << ", ";`	264	`// cout << tim.get_secs() << ", ";`
201	`// }`	265	`// }`
202	`// cout << endl;`	266	`// cout << endl;`
203	`//`	267	`//`
204	`// cout << "Method 3: ";`	268	`// cout << "Method 3: ";`
205	`// for(int i=0;i<5;++i) {`	269	`// for(int i=0;i<5;++i) {`
206	`// Manifold m2 = m;`	270	`// Manifold m2 = m;`
207	`// tim.start();`	271	`// tim.start();`
208	`// laplacian_smooth3(m2, 1, 1000);`	272	`// laplacian_smooth3(m2, 1, 1000);`
209	`// cout << tim.get_secs() << ", ";`	273	`// cout << tim.get_secs() << ", ";`
210	`// }`	274	`// }`
211	`// cout << endl;`	275	`// cout << endl;`
212	`//`	276	`//`
213	`// cout << "Method 4: ";`	277	`// cout << "Method 4: ";`
214	`// for(int i=0;i<5;++i) {`	278	`// for(int i=0;i<5;++i) {`
215	`// Manifold m2 = m;`	279	`// Manifold m2 = m;`
216	`// tim.start();`	280	`// tim.start();`
217	`// laplacian_smooth4(m2, 1, 1000);`	281	`// laplacian_smooth4(m2, 1, 1000);`
218	`// cout << tim.get_secs() << ", ";`	282	`// cout << tim.get_secs() << ", ";`
219	`// }`	283	`// }`
220	`// cout << endl;`	284	`// cout << endl;`
221	`//`	285	`//`
222	`//// cout << "Method 5: ";`	286	`//// cout << "Method 5: ";`
223	`//// for(int i=0;i<3;++i) {`	287	`//// for(int i=0;i<3;++i) {`
224	`//// Manifold m2 = m;`	288	`//// Manifold m2 = m;`
225	`//// tim.start();`	289	`//// tim.start();`
226	`//// laplacian_smooth5(m2, 1, 1000);`	290	`//// laplacian_smooth5(m2, 1, 1000);`
227	`//// cout << " sec: "<< tim.get_secs() << ",";`	291	`//// cout << " sec: "<< tim.get_secs() << ",";`
228	`//// }`	292	`//// }`
229	`//// cout << endl;`	293	`//// cout << endl;`
230	`//`	294	`//`
231	`// CORES = 2;`	295	`// CORES = 2;`
232	`// cout << "Method 6 (2): ";`	296	`// cout << "Method 6 (2): ";`
233	`// for(int i=0;i<5;++i) {`	297	`// for(int i=0;i<5;++i) {`
234	`// Manifold m2 = m;`	298	`// Manifold m2 = m;`
235	`// tim.start();`	299	`// tim.start();`
236	`// laplacian_smooth6(m2, 1, 1000);`	300	`// laplacian_smooth6(m2, 1, 1000);`
237	`// cout << tim.get_secs() << ", ";`	301	`// cout << tim.get_secs() << ", ";`
238	`// }`	302	`// }`
239	`// cout << endl;`	303	`// cout << endl;`
240	`//`	304	`//`
241	`// CORES = 4;`	305	`// CORES = 4;`
242	`// cout << "Method 6 (4): ";`	306	`// cout << "Method 6 (4): ";`
243	`// for(int i=0;i<5;++i) {`	307	`// for(int i=0;i<5;++i) {`
244	`// Manifold m2 = m;`	308	`// Manifold m2 = m;`
245	`// tim.start();`	309	`// tim.start();`
246	`// laplacian_smooth6(m2, 1, 1000);`	310	`// laplacian_smooth6(m2, 1, 1000);`
247	`// cout << tim.get_secs() << ", ";`	311	`// cout << tim.get_secs() << ", ";`
248	`// }`	312	`// }`
249	`// cout << endl;`	313	`// cout << endl;`
250	`//`	314	`//`
251	`// CORES = 8;`	315	`// CORES = 8;`
252	`// cout << "Method 6 (8): ";`	316	`// cout << "Method 6 (8): ";`
253	`// for(int i=0;i<5;++i) {`	317	`// for(int i=0;i<5;++i) {`
254	`// Manifold m2 = m;`	318	`// Manifold m2 = m;`
255	`// tim.start();`	319	`// tim.start();`
256	`// laplacian_smooth6(m2, 1, 1000);`	320	`// laplacian_smooth6(m2, 1, 1000);`
257	`// cout << tim.get_secs() << ", ";`	321	`// cout << tim.get_secs() << ", ";`
258	`// }`	322	`// }`
259	`// cout << endl;`	323	`// cout << endl;`
260		324
261		325
262	`}`	326	`}`
263		327
264		328
265		329
266		330
267	`void taubin_smooth(Manifold& m, int max_iter)`	331	`void taubin_smooth(Manifold& m, int max_iter)`
268	`{`	332	`{`
269	`auto new_pos = m.positions_attribute_vector();`	333	`auto new_pos = m.positions_attribute_vector();`
270	`for(int iter = 0; iter < 2*max_iter; ++iter) {`	334	`for(int iter = 0; iter < 2*max_iter; ++iter) {`
271	`for(VertexID v : m.vertices())`	335	`for(VertexID v : m.vertices())`
272	`if(!boundary(m, v))`	336	`if(!boundary(m, v))`
273	`new_pos[v] = (iter%2 == 0 ? 0.5 : -0.52) * laplacian(m, v) + m.pos(v);`	337	`new_pos[v] = (iter%2 == 0 ? 0.5 : -0.52) * laplacian(m, v) + m.pos(v);`
274	`swap(m.positions_attribute_vector(), new_pos);`	338	`swap(m.positions_attribute_vector(), new_pos);`
275	`}`	339	`}`
276	`}`	340	`}`
277		341
278	`void face_neighbourhood(Manifold& m, FaceID f, vector<FaceID>& nbrs)`	342	`void face_neighbourhood(Manifold& m, FaceID f, vector<FaceID>& nbrs)`
279	`{`	343	`{`
280	`FaceAttributeVector<int> touched(m.allocated_faces(), -1);`	344	`FaceAttributeVector<int> touched(m.allocated_faces(), -1);`
281	`touched[f] = 1;`	345	`touched[f] = 1;`
282	`nbrs.push_back(f);`	346	`nbrs.push_back(f);`
283	`for(Walker wf = m.walker(f); !wf.full_circle(); wf = wf.circulate_face_cw()){`	347	`for(Walker wf = m.walker(f); !wf.full_circle(); wf = wf.circulate_face_cw()){`
284	`for(Walker wv = m.walker(wf.vertex()); !wv.full_circle(); wv = wv.circulate_vertex_cw()){`	348	`for(Walker wv = m.walker(wf.vertex()); !wv.full_circle(); wv = wv.circulate_vertex_cw()){`
285	`FaceID fn = wv.face();`	349	`FaceID fn = wv.face();`
286	`if(fn != InvalidFaceID && touched[fn] != touched[f]){`	350	`if(fn != InvalidFaceID && touched[fn] != touched[f]){`
287	`nbrs.push_back(fn);`	351	`nbrs.push_back(fn);`
288	`touched[fn] = 1;`	352	`touched[fn] = 1;`
289	`}`	353	`}`
290	`}`	354	`}`
291	`}`	355	`}`
292	`}`	356	`}`
293		357
294		358
295		359
296	`Vec3d fvm_filtered_normal(Manifold& m, FaceID f)`	360	`Vec3d fvm_filtered_normal(Manifold& m, FaceID f)`
297	`{`	361	`{`
298	`const float sigma = .1f;`	362	`const float sigma = .1f;`
299		363
300	`vector<FaceID> nbrs;`	364	`vector<FaceID> nbrs;`
301	`face_neighbourhood(m, f, nbrs);`	365	`face_neighbourhood(m, f, nbrs);`
302	`float min_dist_sum=1e32f;`	366	`float min_dist_sum=1e32f;`
303	`long int median=-1;`	367	`long int median=-1;`
304		368
305	`vector<Vec3d> normals(nbrs.size());`	369	`vector<Vec3d> normals(nbrs.size());`
306	`for(size_t i=0;i<nbrs.size();++i)`	370	`for(size_t i=0;i<nbrs.size();++i)`
307	`{`	371	`{`
308	`normals[i] = normal(m, nbrs[i]);`	372	`normals[i] = normal(m, nbrs[i]);`
309	`float dist_sum = 0;`	373	`float dist_sum = 0;`
310	`for(size_t j=0;j<nbrs.size(); ++j)`	374	`for(size_t j=0;j<nbrs.size(); ++j)`
311	`dist_sum += 1.0f - dot(normals[i], normals[j]);`	375	`dist_sum += 1.0f - dot(normals[i], normals[j]);`
312	`if(dist_sum < min_dist_sum)`	376	`if(dist_sum < min_dist_sum)`
313	`{`	377	`{`
314	`min_dist_sum = dist_sum;`	378	`min_dist_sum = dist_sum;`
315	`median = i;`	379	`median = i;`
316	`}`	380	`}`
317	`}`	381	`}`
318	`assert(median != -1);`	382	`assert(median != -1);`
319	`Vec3d median_norm = normals[median];`	383	`Vec3d median_norm = normals[median];`
320	`Vec3d avg_norm(0);`	384	`Vec3d avg_norm(0);`
321	`for(size_t i=0;i<nbrs.size();++i)`	385	`for(size_t i=0;i<nbrs.size();++i)`
322	`{`	386	`{`
323	`float w = exp((dot(median_norm, normals[i])-1)/sigma);`	387	`float w = exp((dot(median_norm, normals[i])-1)/sigma);`
324	`if(w<1e-2) w = 0;`	388	`if(w<1e-2) w = 0;`
325	`avg_norm += w*normals[i];`	389	`avg_norm += w*normals[i];`
326	`}`	390	`}`
327	`return normalize(avg_norm);`	391	`return normalize(avg_norm);`
328	`}`	392	`}`
329	`Vec3d bilateral_filtered_normal(Manifold& m, FaceID f, double avg_len)`	393	`Vec3d bilateral_filtered_normal(Manifold& m, FaceID f, double avg_len)`
330	`{`	394	`{`
331	`vector<FaceID> nbrs;`	395	`vector<FaceID> nbrs;`
332	`face_neighbourhood(m, f, nbrs);`	396	`face_neighbourhood(m, f, nbrs);`
333	`Vec3d p0 = centre(m, f);`	397	`Vec3d p0 = centre(m, f);`
334	`Vec3d n0 = normal(m, f);`	398	`Vec3d n0 = normal(m, f);`
335	`vector<Vec3d> normals(nbrs.size());`	399	`vector<Vec3d> normals(nbrs.size());`
336	`vector<Vec3d> pos(nbrs.size());`	400	`vector<Vec3d> pos(nbrs.size());`
337	`Vec3d fn(0);`	401	`Vec3d fn(0);`
338	`for(FaceID nbr : nbrs)`	402	`for(FaceID nbr : nbrs)`
339	`{`	403	`{`
340	`Vec3d n = normal(m, nbr);`	404	`Vec3d n = normal(m, nbr);`
341	`Vec3d p = centre(m, nbr);`	405	`Vec3d p = centre(m, nbr);`
342	`double w_a = exp(-acos(max(-1.0,min(1.0,dot(n,n0))))/(M_PI/32.0));`	406	`double w_a = exp(-acos(max(-1.0,min(1.0,dot(n,n0))))/(M_PI/32.0));`
343	`double w_s = exp(-length(p-p0)/avg_len);`	407	`double w_s = exp(-length(p-p0)/avg_len);`
344		408
345	`fn += area(m, nbr)* w_a * w_s * n;`	409	`fn += area(m, nbr)* w_a * w_s * n;`
346		410
347	`}`	411	`}`
348	`return normalize(fn);`	412	`return normalize(fn);`
349	`}`	413	`}`
350		414
351	`void anisotropic_smooth(HMesh::Manifold& m, int max_iter, NormalSmoothMethod nsm)`	415	`void anisotropic_smooth(HMesh::Manifold& m, int max_iter, NormalSmoothMethod nsm)`
352	`{`	416	`{`
353	`double avg_len=0;`	417	`double avg_len=0;`
354	`for(HalfEdgeID hid : m.halfedges())`	418	`for(HalfEdgeID hid : m.halfedges())`
355	`avg_len += length(m, hid);`	419	`avg_len += length(m, hid);`
356	`avg_len /= 2.0;`	420	`avg_len /= 2.0;`
357	`for(int iter = 0;iter<max_iter; ++iter)`	421	`for(int iter = 0;iter<max_iter; ++iter)`
358	`{`	422	`{`
359		423
360	`FaceAttributeVector<Vec3d> filtered_norms(m.allocated_faces());`	424	`FaceAttributeVector<Vec3d> filtered_norms(m.allocated_faces());`
361		425
362	`for(FaceID f: m.faces()){`	426	`for(FaceID f: m.faces()){`
363	`filtered_norms[f] = (nsm == BILATERAL_NORMAL_SMOOTH)?`	427	`filtered_norms[f] = (nsm == BILATERAL_NORMAL_SMOOTH)?`
364	`bilateral_filtered_normal(m, f, avg_len):`	428	`bilateral_filtered_normal(m, f, avg_len):`
365	`fvm_filtered_normal(m, f);`	429	`fvm_filtered_normal(m, f);`
366	`}`	430	`}`
367		431
368	`VertexAttributeVector<Vec3d> vertex_positions(m.allocated_vertices(), Vec3d(0));`	432	`VertexAttributeVector<Vec3d> vertex_positions(m.allocated_vertices(), Vec3d(0));`
369	`VertexAttributeVector<int> count(m.allocated_vertices(), 0);`	433	`VertexAttributeVector<int> count(m.allocated_vertices(), 0);`
370	`for(int sub_iter=0;sub_iter<100;++sub_iter)`	434	`for(int sub_iter=0;sub_iter<100;++sub_iter)`
371	`{`	435	`{`
372	`for(HalfEdgeID hid : m.halfedges())`	436	`for(HalfEdgeID hid : m.halfedges())`
373	`{`	437	`{`
374	`Walker w = m.walker(hid);`	438	`Walker w = m.walker(hid);`
375	`FaceID f = w.face();`	439	`FaceID f = w.face();`
376		440
377	`if(f != InvalidFaceID){`	441	`if(f != InvalidFaceID){`
378	`VertexID v = w.vertex();`	442	`VertexID v = w.vertex();`
379	`Vec3d dir = m.pos(w.opp().vertex()) - m.pos(v);`	443	`Vec3d dir = m.pos(w.opp().vertex()) - m.pos(v);`
380	`Vec3d n = filtered_norms[f];`	444	`Vec3d n = filtered_norms[f];`
381	`vertex_positions[v] += m.pos(v) + 0.5 * n * dot(n, dir);`	445	`vertex_positions[v] += m.pos(v) + 0.5 * n * dot(n, dir);`
382	`count[v] += 1;`	446	`count[v] += 1;`
383	`}`	447	`}`
384		448
385	`}`	449	`}`
386	`double max_move= 0;`	450	`double max_move= 0;`
387	`for(VertexID v : m.vertices())`	451	`for(VertexID v : m.vertices())`
388	`{`	452	`{`
389	`Vec3d npos = vertex_positions[v] / double(count[v]);`	453	`Vec3d npos = vertex_positions[v] / double(count[v]);`
390	`double move = sqr_length(npos - m.pos(v));`	454	`double move = sqr_length(npos - m.pos(v));`
391	`if(move > max_move)`	455	`if(move > max_move)`
392	`max_move = move;`	456	`max_move = move;`
393		457
394	`m.pos(v) = npos;`	458	`m.pos(v) = npos;`
395	`}`	459	`}`
396	`if(max_move<sqr(1e-8*avg_len))`	460	`if(max_move<sqr(1e-8*avg_len))`
397	`{`	461	`{`
398	`cout << "iters " << sub_iter << endl;`	462	`cout << "iters " << sub_iter << endl;`
399	`break;`	463	`break;`
400	`}`	464	`}`
401	`}`	465	`}`
402	`}`	466	`}`
403	`}`	467	`}`
404		468
405		469
406	`}`	470	`}`
407		471

Subversion Repositories gelsvn

(root)/branches/cpp11-devel/src/HMesh/smooth.cpp – Rev 632 → 633