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

Subversion Repositories gelsvn

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