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;
    
 
 
    
    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); }
 
    
    Vec3d laplacian(const Manifold& m, VertexID v)
    {
        Vec3d avg_pos(0);
        int n = 0;
        
        circulate_vertex_ccw(m, v, [&](VertexID v){
            avg_pos += m.pos(v);
            ++n;
        });
        return avg_pos / 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,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);
        }
        
    }
    
    void laplacian_smooth(Manifold& m,float weight, int max_iter)
    {
        laplacian_smooth4(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 631	Rev 632
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		25
26		26
27		27
28	`void for_each_vertex(Manifold& m, function<void(VertexID)> f) { for(auto v : m.vertices()) f(v); }`	28	`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); }`	29	`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); }`	30	`void for_each_halfedge(Manifold& m, function<void(HalfEdgeID)> f) { for(auto h : m.halfedges()) f(h); }`
31		31
32		32
33	`Vec3d laplacian(const Manifold& m, VertexID v)`	33	`Vec3d laplacian(const Manifold& m, VertexID v)`
34	`{`	34	`{`
35	`Vec3d avg_pos(0);`	35	`Vec3d avg_pos(0);`
36	`int n = 0;`	36	`int n = 0;`
37		37
38	`circulate_vertex_ccw(m, v, [&](VertexID v){`	38	`circulate_vertex_ccw(m, v, [&](VertexID v){`
39	`avg_pos += m.pos(v);`	39	`avg_pos += m.pos(v);`
40	`++n;`	40	`++n;`
41	`});`	41	`});`
42	`return avg_pos / n - m.pos(v);`	42	`return avg_pos / n - m.pos(v);`
43	`}`	43	`}`
44		44
45		45
46	`int CORES = 8;`	46	`int CORES = 8;`
47		47
48	`void laplacian_smooth0(Manifold& m,float weight, int max_iter)`	48	`void laplacian_smooth0(Manifold& m,float weight, int max_iter)`
49	`{`	49	`{`
50	`for(int iter=0;iter<max_iter; ++iter) {`	50	`for(int iter=0;iter<max_iter; ++iter) {`
51	`VertexAttributeVector<Vec3d> L_attr(m.no_vertices());`	51	`VertexAttributeVector<Vec3d> L_attr(m.no_vertices());`
52		52
53	`for(VertexIDIterator v = m.vertices_begin(); v != m.vertices_end(); ++v)`	53	`for(VertexIDIterator v = m.vertices_begin(); v != m.vertices_end(); ++v)`
54	`if(!boundary(m, *v))`	54	`if(!boundary(m, *v))`
55	`L_attr[v] =laplacian(m, v);`	55	`L_attr[v] =laplacian(m, v);`
56	`for(VertexIDIterator v = m.vertices_begin(); v != m.vertices_end(); ++v){`	56	`for(VertexIDIterator v = m.vertices_begin(); v != m.vertices_end(); ++v){`
57	`if(!boundary(m, *v))`	57	`if(!boundary(m, *v))`
58	`m.pos(v) += weightL_attr[*v];`	58	`m.pos(v) += weightL_attr[*v];`
59	`}`	59	`}`
60	`}`	60	`}`
61	`}`	61	`}`
62		62
63	`void laplacian_smooth1(Manifold& m,float weight, int max_iter)`	63	`void laplacian_smooth1(Manifold& m,float weight, int max_iter)`
64	`{`	64	`{`
65	`for(int iter=0;iter<max_iter; ++iter) {`	65	`for(int iter=0;iter<max_iter; ++iter) {`
66	`VertexAttributeVector<Vec3d> L_attr(m.no_vertices());`	66	`VertexAttributeVector<Vec3d> L_attr(m.no_vertices());`
67		67
68	`for(VertexID v : m.vertices())`	68	`for(VertexID v : m.vertices())`
69	`if(!boundary(m, v))`	69	`if(!boundary(m, v))`
70	`L_attr[v] =laplacian(m, v);`	70	`L_attr[v] =laplacian(m, v);`
71	`for(VertexID v : m.vertices()){`	71	`for(VertexID v : m.vertices()){`
72	`if(!boundary(m, v))`	72	`if(!boundary(m, v))`
73	`m.pos(v) += weight*L_attr[v];`	73	`m.pos(v) += weight*L_attr[v];`
74	`}`	74	`}`
75	`}`	75	`}`
76	`}`	76	`}`
77		77
78	`void laplacian_smooth2(Manifold& m,float weight, int max_iter)`	78	`void laplacian_smooth2(Manifold& m,float weight, int max_iter)`
79	`{`	79	`{`
80	`auto new_pos = m.positions_attribute_vector();`	80	`auto new_pos = m.positions_attribute_vector();`
81	`for(int iter=0;iter<max_iter; ++iter) {`	81	`for(int iter=0;iter<max_iter; ++iter) {`
82	`for(auto v : m.vertices())`	82	`for(auto v : m.vertices())`
83	`if(!boundary(m, v))`	83	`if(!boundary(m, v))`
84	`new_pos[v] = weight*laplacian(m, v)+m.pos(v);`	84	`new_pos[v] = weight*laplacian(m, v)+m.pos(v);`
85	`m.positions_attribute_vector() = new_pos;`	85	`m.positions_attribute_vector() = new_pos;`
86	`}`	86	`}`
87	`}`	87	`}`
88		88
89		89
90		90
91	`void laplacian_smooth3(Manifold& m, float weight, int max_iter)`	91	`void laplacian_smooth3(Manifold& m, float weight, int max_iter)`
92	`{`	92	`{`
93	`for(int iter=0;iter<max_iter; ++iter) {`	93	`for(int iter=0;iter<max_iter; ++iter) {`
94	`auto new_pos = m.positions_attribute_vector();`	94	`auto new_pos = m.positions_attribute_vector();`
95	`for(auto v : m.vertices())`	95	`for(auto v : m.vertices())`
96	`if(!boundary(m, v))`	96	`if(!boundary(m, v))`
97	`new_pos[v] = weight*laplacian(m, v)+m.pos(v);`	97	`new_pos[v] = weight*laplacian(m, v)+m.pos(v);`
98	`m.positions_attribute_vector() = move(new_pos);`	98	`m.positions_attribute_vector() = move(new_pos);`
99	`}`	99	`}`
100	`}`	100	`}`
101		101
102	`void laplacian_smooth4(Manifold& m,float weight, int max_iter)`	102	`void laplacian_smooth4(Manifold& m,float weight, int max_iter)`
103	`{`	103	`{`
104	`auto new_pos = m.positions_attribute_vector();`	104	`auto new_pos = m.positions_attribute_vector();`
105	`for(int iter=0;iter<max_iter; ++iter) {`	105	`for(int iter=0;iter<max_iter; ++iter) {`
106	`for(auto v : m.vertices())`	106	`for(VertexID v : m.vertices())`
107	`if(!boundary(m, v))`	107	`if(!boundary(m, v))`
108	`new_pos[v] = weight*laplacian(m, v)+m.pos(v);`	108	`new_pos[v] = weight*laplacian(m, v)+m.pos(v);`
109	`swap(m.positions_attribute_vector(),new_pos);`	109	`swap(m.positions_attribute_vector(),new_pos);`
110	`}`	110	`}`
111	`}`	111	`}`
112		112
113	`void laplacian_smooth4_5(Manifold& m,float weight, int max_iter)`	113	`void laplacian_smooth4_5(Manifold& m,float weight, int max_iter)`
114	`{`	114	`{`
115	`auto new_pos = m.positions_attribute_vector();`	115	`auto new_pos = m.positions_attribute_vector();`
116	`for(int iter=0;iter<max_iter; ++iter) {`	116	`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);});`	117	`for_each_vertex(m, [&](VertexID v) {new_pos[v] = weight*laplacian(m, v)+m.pos(v);});`
118	`swap(m.positions_attribute_vector(),new_pos);`	118	`swap(m.positions_attribute_vector(),new_pos);`
119	`}`	119	`}`
120	`}`	120	`}`
121		121
122		122
123		123
124	`void laplacian_smooth5(Manifold& m, float weight, int max_iter)`	124	`void laplacian_smooth5(Manifold& m, float weight, int max_iter)`
125	`{`	125	`{`
126	`for(int iter=0;iter<max_iter; ++iter) {`	126	`for(int iter=0;iter<max_iter; ++iter) {`
127	`auto new_pos = m.positions_attribute_vector();`	127	`auto new_pos = m.positions_attribute_vector();`
128	`vector<thread> t_vec;`	128	`vector<thread> t_vec;`
129	`for(auto v : m.vertices())`	129	`for(auto v : m.vertices())`
130	`if(!boundary(m, v))`	130	`if(!boundary(m, v))`
131	`t_vec.push_back(thread([&](VertexID vid){`	131	`t_vec.push_back(thread([&](VertexID vid){`
132	`if(!boundary(m, vid))`	132	`if(!boundary(m, vid))`
133	`new_pos[vid] = weight*laplacian(m, vid)+ m.pos(vid);},v));`	133	`new_pos[vid] = weight*laplacian(m, vid)+ m.pos(vid);},v));`
134	`for(int i=0;i<t_vec.size();++i)`	134	`for(int i=0;i<t_vec.size();++i)`
135	`t_vec[i].join();`	135	`t_vec[i].join();`
136	`m.positions_attribute_vector() = move(new_pos);`	136	`m.positions_attribute_vector() = move(new_pos);`
137	`}`	137	`}`
138	`}`	138	`}`
139		139
140		140
141		141
142	`inline void laplacian_smooth_vertex(Manifold& m,vector<VertexID>& vids,`	142	`inline void laplacian_smooth_vertex(Manifold& m,vector<VertexID>& vids,`
143	`VertexAttributeVector<Vec3d>& new_pos,`	143	`VertexAttributeVector<Vec3d>& new_pos,`
144	`float weight){`	144	`float weight){`
145	`for(auto v: vids)`	145	`for(VertexID v: vids)`
146	`new_pos[v] = m.pos(v)+weight*laplacian(m, v);`	146	`new_pos[v] = m.pos(v)+weight*laplacian(m, v);`
147	`}`	147	`}`
148		148
149	`void laplacian_smooth6(Manifold& m, float weight, int max_iter)`	149	`void laplacian_smooth6(Manifold& m, float weight, int max_iter)`
150	`{`	150	`{`
151	`vector<vector<VertexID>> vertex_ids(CORES);`	151	`vector<vector<VertexID>> vertex_ids(CORES);`
152	`auto batch_size = m.no_vertices()/CORES;`	152	`auto batch_size = m.no_vertices()/CORES;`
153	`int cnt = 0;`	153	`int cnt = 0;`
154	`for_each_vertex(m, [&](VertexID v) {`	154	`for_each_vertex(m, [&](VertexID v) {`
155	`if (!boundary(m, v))`	155	`if (!boundary(m, v))`
156	`vertex_ids[(cnt++/batch_size)%CORES].push_back(v);`	156	`vertex_ids[(cnt++/batch_size)%CORES].push_back(v);`
157	`});`	157	`});`
158	`vector<thread> t_vec(CORES);`	158	`vector<thread> t_vec(CORES);`
159	`VertexAttributeVector<Vec3d> new_pos = m.positions_attribute_vector();`	159	`VertexAttributeVector<Vec3d> new_pos = m.positions_attribute_vector();`
160	`for(int iter=0;iter<max_iter; ++iter) {`	160	`for(int iter=0;iter<max_iter; ++iter) {`
161	`for(int thread_no=0;thread_no<CORES;++thread_no)`	161	`for(int thread_no=0;thread_no<CORES;++thread_no)`
162	`t_vec[thread_no] = thread(laplacian_smooth_vertex,`	162	`t_vec[thread_no] = thread(laplacian_smooth_vertex,`
163	`ref(m), ref(vertex_ids[thread_no]),`	163	`ref(m), ref(vertex_ids[thread_no]),`
164	`ref(new_pos), weight);`	164	`ref(new_pos), weight);`
165	`for(int thread_no=0;thread_no<CORES;++thread_no)`	165	`for(int thread_no=0;thread_no<CORES;++thread_no)`
166	`t_vec[thread_no].join();`	166	`t_vec[thread_no].join();`
167	`swap(m.positions_attribute_vector(), new_pos);`	167	`swap(m.positions_attribute_vector(), new_pos);`
168	`}`	168	`}`
169		169
170	`}`	170	`}`
171		171
172	`void laplacian_smooth(Manifold& m,float weight, int max_iter)`	172	`void laplacian_smooth(Manifold& m,float weight, int max_iter)`
173	`{`	173	`{`
174	`laplacian_smooth4(m,weight,max_iter);`	174	`laplacian_smooth4(m,weight,max_iter);`
175	`// Util::Timer tim;`	175	`// Util::Timer tim;`
176	`//`	176	`//`
177	`// cout << "Method 0: ";`	177	`// cout << "Method 0: ";`
178	`// for(int i=0;i<5;++i) {`	178	`// for(int i=0;i<5;++i) {`
179	`// Manifold m2 = m;`	179	`// Manifold m2 = m;`
180	`// tim.start();`	180	`// tim.start();`
181	`// laplacian_smooth0(m2, 1, 1000);`	181	`// laplacian_smooth0(m2, 1, 1000);`
182	`// cout << tim.get_secs() << ", ";`	182	`// cout << tim.get_secs() << ", ";`
183	`// }`	183	`// }`
184	`// cout << endl;`	184	`// cout << endl;`
185	`//`	185	`//`
186	`// cout << "Method 1: ";`	186	`// cout << "Method 1: ";`
187	`// for(int i=0;i<5;++i) {`	187	`// for(int i=0;i<5;++i) {`
188	`// Manifold m2 = m;`	188	`// Manifold m2 = m;`
189	`// tim.start();`	189	`// tim.start();`
190	`// laplacian_smooth1(m2, 1, 1000);`	190	`// laplacian_smooth1(m2, 1, 1000);`
191	`// cout << tim.get_secs() << ", ";`	191	`// cout << tim.get_secs() << ", ";`
192	`// }`	192	`// }`
193	`// cout << endl;`	193	`// cout << endl;`
194	`//`	194	`//`
195	`// cout << "Method 2: ";`	195	`// cout << "Method 2: ";`
196	`// for(int i=0;i<5;++i) {`	196	`// for(int i=0;i<5;++i) {`
197	`// Manifold m2 = m;`	197	`// Manifold m2 = m;`
198	`// tim.start();`	198	`// tim.start();`
199	`// laplacian_smooth2(m2, 1, 1000);`	199	`// laplacian_smooth2(m2, 1, 1000);`
200	`// cout << tim.get_secs() << ", ";`	200	`// cout << tim.get_secs() << ", ";`
201	`// }`	201	`// }`
202	`// cout << endl;`	202	`// cout << endl;`
203	`//`	203	`//`
204	`// cout << "Method 3: ";`	204	`// cout << "Method 3: ";`
205	`// for(int i=0;i<5;++i) {`	205	`// for(int i=0;i<5;++i) {`
206	`// Manifold m2 = m;`	206	`// Manifold m2 = m;`
207	`// tim.start();`	207	`// tim.start();`
208	`// laplacian_smooth3(m2, 1, 1000);`	208	`// laplacian_smooth3(m2, 1, 1000);`
209	`// cout << tim.get_secs() << ", ";`	209	`// cout << tim.get_secs() << ", ";`
210	`// }`	210	`// }`
211	`// cout << endl;`	211	`// cout << endl;`
212	`//`	212	`//`
213	`// cout << "Method 4: ";`	213	`// cout << "Method 4: ";`
214	`// for(int i=0;i<5;++i) {`	214	`// for(int i=0;i<5;++i) {`
215	`// Manifold m2 = m;`	215	`// Manifold m2 = m;`
216	`// tim.start();`	216	`// tim.start();`
217	`// laplacian_smooth4(m2, 1, 1000);`	217	`// laplacian_smooth4(m2, 1, 1000);`
218	`// cout << tim.get_secs() << ", ";`	218	`// cout << tim.get_secs() << ", ";`
219	`// }`	219	`// }`
220	`// cout << endl;`	220	`// cout << endl;`
221	`//`	221	`//`
222	`//// cout << "Method 5: ";`	222	`//// cout << "Method 5: ";`
223	`//// for(int i=0;i<3;++i) {`	223	`//// for(int i=0;i<3;++i) {`
224	`//// Manifold m2 = m;`	224	`//// Manifold m2 = m;`
225	`//// tim.start();`	225	`//// tim.start();`
226	`//// laplacian_smooth5(m2, 1, 1000);`	226	`//// laplacian_smooth5(m2, 1, 1000);`
227	`//// cout << " sec: "<< tim.get_secs() << ",";`	227	`//// cout << " sec: "<< tim.get_secs() << ",";`
228	`//// }`	228	`//// }`
229	`//// cout << endl;`	229	`//// cout << endl;`
230	`//`	230	`//`
231	`// CORES = 2;`	231	`// CORES = 2;`
232	`// cout << "Method 6 (2): ";`	232	`// cout << "Method 6 (2): ";`
233	`// for(int i=0;i<5;++i) {`	233	`// for(int i=0;i<5;++i) {`
234	`// Manifold m2 = m;`	234	`// Manifold m2 = m;`
235	`// tim.start();`	235	`// tim.start();`
236	`// laplacian_smooth6(m2, 1, 1000);`	236	`// laplacian_smooth6(m2, 1, 1000);`
237	`// cout << tim.get_secs() << ", ";`	237	`// cout << tim.get_secs() << ", ";`
238	`// }`	238	`// }`
239	`// cout << endl;`	239	`// cout << endl;`
240	`//`	240	`//`
241	`// CORES = 4;`	241	`// CORES = 4;`
242	`// cout << "Method 6 (4): ";`	242	`// cout << "Method 6 (4): ";`
243	`// for(int i=0;i<5;++i) {`	243	`// for(int i=0;i<5;++i) {`
244	`// Manifold m2 = m;`	244	`// Manifold m2 = m;`
245	`// tim.start();`	245	`// tim.start();`
246	`// laplacian_smooth6(m2, 1, 1000);`	246	`// laplacian_smooth6(m2, 1, 1000);`
247	`// cout << tim.get_secs() << ", ";`	247	`// cout << tim.get_secs() << ", ";`
248	`// }`	248	`// }`
249	`// cout << endl;`	249	`// cout << endl;`
250	`//`	250	`//`
251	`// CORES = 8;`	251	`// CORES = 8;`
252	`// cout << "Method 6 (8): ";`	252	`// cout << "Method 6 (8): ";`
253	`// for(int i=0;i<5;++i) {`	253	`// for(int i=0;i<5;++i) {`
254	`// Manifold m2 = m;`	254	`// Manifold m2 = m;`
255	`// tim.start();`	255	`// tim.start();`
256	`// laplacian_smooth6(m2, 1, 1000);`	256	`// laplacian_smooth6(m2, 1, 1000);`
257	`// cout << tim.get_secs() << ", ";`	257	`// cout << tim.get_secs() << ", ";`
258	`// }`	258	`// }`
259	`// cout << endl;`	259	`// cout << endl;`
260		260
261		261
262	`}`	262	`}`
263		263
264		264
265		265
266		266
267	`void taubin_smooth(Manifold& m, int max_iter)`	267	`void taubin_smooth(Manifold& m, int max_iter)`
268	`{`	268	`{`
269	`auto new_pos = m.positions_attribute_vector();`	269	`auto new_pos = m.positions_attribute_vector();`
270	`for(int iter = 0; iter < 2*max_iter; ++iter) {`	270	`for(int iter = 0; iter < 2*max_iter; ++iter) {`
271	`for(VertexID v : m.vertices())`	271	`for(VertexID v : m.vertices())`
272	`if(!boundary(m, v))`	272	`if(!boundary(m, v))`
273	`new_pos[v] = (iter%2 == 0 ? 0.5 : -0.52) * laplacian(m, v) + m.pos(v);`	273	`new_pos[v] = (iter%2 == 0 ? 0.5 : -0.52) * laplacian(m, v) + m.pos(v);`
274	`swap(m.positions_attribute_vector(), new_pos);`	274	`swap(m.positions_attribute_vector(), new_pos);`
275	`}`	275	`}`
276	`}`	276	`}`
277		277
278	`void face_neighbourhood(Manifold& m, FaceID f, vector<FaceID>& nbrs)`	278	`void face_neighbourhood(Manifold& m, FaceID f, vector<FaceID>& nbrs)`
279	`{`	279	`{`
280	`FaceAttributeVector<int> touched(m.allocated_faces(), -1);`	280	`FaceAttributeVector<int> touched(m.allocated_faces(), -1);`
281	`touched[f] = 1;`	281	`touched[f] = 1;`
282	`nbrs.push_back(f);`	282	`nbrs.push_back(f);`
283	`for(Walker wf = m.walker(f); !wf.full_circle(); wf = wf.circulate_face_cw()){`	283	`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()){`	284	`for(Walker wv = m.walker(wf.vertex()); !wv.full_circle(); wv = wv.circulate_vertex_cw()){`
285	`FaceID fn = wv.face();`	285	`FaceID fn = wv.face();`
286	`if(fn != InvalidFaceID && touched[fn] != touched[f]){`	286	`if(fn != InvalidFaceID && touched[fn] != touched[f]){`
287	`nbrs.push_back(fn);`	287	`nbrs.push_back(fn);`
288	`touched[fn] = 1;`	288	`touched[fn] = 1;`
289	`}`	289	`}`
290	`}`	290	`}`
291	`}`	291	`}`
292	`}`	292	`}`
293		293
294		294
295		295
296	`Vec3d fvm_filtered_normal(Manifold& m, FaceID f)`	296	`Vec3d fvm_filtered_normal(Manifold& m, FaceID f)`
297	`{`	297	`{`
298	`const float sigma = .1f;`	298	`const float sigma = .1f;`
299		299
300	`vector<FaceID> nbrs;`	300	`vector<FaceID> nbrs;`
301	`face_neighbourhood(m, f, nbrs);`	301	`face_neighbourhood(m, f, nbrs);`
302	`float min_dist_sum=1e32f;`	302	`float min_dist_sum=1e32f;`
303	`long int median=-1;`	303	`long int median=-1;`
304		304
305	`vector<Vec3d> normals(nbrs.size());`	305	`vector<Vec3d> normals(nbrs.size());`
306	`for(size_t i=0;i<nbrs.size();++i)`	306	`for(size_t i=0;i<nbrs.size();++i)`
307	`{`	307	`{`
308	`normals[i] = normal(m, nbrs[i]);`	308	`normals[i] = normal(m, nbrs[i]);`
309	`float dist_sum = 0;`	309	`float dist_sum = 0;`
310	`for(size_t j=0;j<nbrs.size(); ++j)`	310	`for(size_t j=0;j<nbrs.size(); ++j)`
311	`dist_sum += 1.0f - dot(normals[i], normals[j]);`	311	`dist_sum += 1.0f - dot(normals[i], normals[j]);`
312	`if(dist_sum < min_dist_sum)`	312	`if(dist_sum < min_dist_sum)`
313	`{`	313	`{`
314	`min_dist_sum = dist_sum;`	314	`min_dist_sum = dist_sum;`
315	`median = i;`	315	`median = i;`
316	`}`	316	`}`
317	`}`	317	`}`
318	`assert(median != -1);`	318	`assert(median != -1);`
319	`Vec3d median_norm = normals[median];`	319	`Vec3d median_norm = normals[median];`
320	`Vec3d avg_norm(0);`	320	`Vec3d avg_norm(0);`
321	`for(size_t i=0;i<nbrs.size();++i)`	321	`for(size_t i=0;i<nbrs.size();++i)`
322	`{`	322	`{`
323	`float w = exp((dot(median_norm, normals[i])-1)/sigma);`	323	`float w = exp((dot(median_norm, normals[i])-1)/sigma);`
324	`if(w<1e-2) w = 0;`	324	`if(w<1e-2) w = 0;`
325	`avg_norm += w*normals[i];`	325	`avg_norm += w*normals[i];`
326	`}`	326	`}`
327	`return normalize(avg_norm);`	327	`return normalize(avg_norm);`
328	`}`	328	`}`
329	`Vec3d bilateral_filtered_normal(Manifold& m, FaceID f, double avg_len)`	329	`Vec3d bilateral_filtered_normal(Manifold& m, FaceID f, double avg_len)`
330	`{`	330	`{`
331	`vector<FaceID> nbrs;`	331	`vector<FaceID> nbrs;`
332	`face_neighbourhood(m, f, nbrs);`	332	`face_neighbourhood(m, f, nbrs);`
333	`Vec3d p0 = centre(m, f);`	333	`Vec3d p0 = centre(m, f);`
334	`Vec3d n0 = normal(m, f);`	334	`Vec3d n0 = normal(m, f);`
335	`vector<Vec3d> normals(nbrs.size());`	335	`vector<Vec3d> normals(nbrs.size());`
336	`vector<Vec3d> pos(nbrs.size());`	336	`vector<Vec3d> pos(nbrs.size());`
337	`Vec3d fn(0);`	337	`Vec3d fn(0);`
338	`for(FaceID nbr : nbrs)`	338	`for(FaceID nbr : nbrs)`
339	`{`	339	`{`
340	`Vec3d n = normal(m, nbr);`	340	`Vec3d n = normal(m, nbr);`
341	`Vec3d p = centre(m, nbr);`	341	`Vec3d p = centre(m, nbr);`
342	`double w_a = exp(-acos(max(-1.0,min(1.0,dot(n,n0))))/(M_PI/32.0));`	342	`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);`	343	`double w_s = exp(-length(p-p0)/avg_len);`
344		344
345	`fn += area(m, nbr)* w_a * w_s * n;`	345	`fn += area(m, nbr)* w_a * w_s * n;`
346		346
347	`}`	347	`}`
348	`return normalize(fn);`	348	`return normalize(fn);`
349	`}`	349	`}`
350		350
351	`void anisotropic_smooth(HMesh::Manifold& m, int max_iter, NormalSmoothMethod nsm)`	351	`void anisotropic_smooth(HMesh::Manifold& m, int max_iter, NormalSmoothMethod nsm)`
352	`{`	352	`{`
353	`double avg_len=0;`	353	`double avg_len=0;`
354	`for(HalfEdgeID hid : m.halfedges())`	354	`for(HalfEdgeID hid : m.halfedges())`
355	`avg_len += length(m, hid);`	355	`avg_len += length(m, hid);`
356	`avg_len /= 2.0;`	356	`avg_len /= 2.0;`
357	`for(int iter = 0;iter<max_iter; ++iter)`	357	`for(int iter = 0;iter<max_iter; ++iter)`
358	`{`	358	`{`
359		359
360	`FaceAttributeVector<Vec3d> filtered_norms(m.allocated_faces());`	360	`FaceAttributeVector<Vec3d> filtered_norms(m.allocated_faces());`
361		361
362	`for(FaceID f: m.faces()){`	362	`for(FaceID f: m.faces()){`
363	`filtered_norms[f] = (nsm == BILATERAL_NORMAL_SMOOTH)?`	363	`filtered_norms[f] = (nsm == BILATERAL_NORMAL_SMOOTH)?`
364	`bilateral_filtered_normal(m, f, avg_len):`	364	`bilateral_filtered_normal(m, f, avg_len):`
365	`fvm_filtered_normal(m, f);`	365	`fvm_filtered_normal(m, f);`
366	`}`	366	`}`
367		367
368	`VertexAttributeVector<Vec3d> vertex_positions(m.allocated_vertices(), Vec3d(0));`	368	`VertexAttributeVector<Vec3d> vertex_positions(m.allocated_vertices(), Vec3d(0));`
369	`VertexAttributeVector<int> count(m.allocated_vertices(), 0);`	369	`VertexAttributeVector<int> count(m.allocated_vertices(), 0);`
370	`for(int sub_iter=0;sub_iter<100;++sub_iter)`	370	`for(int sub_iter=0;sub_iter<100;++sub_iter)`
371	`{`	371	`{`
372	`for(HalfEdgeID hid : m.halfedges())`	372	`for(HalfEdgeID hid : m.halfedges())`
373	`{`	373	`{`
374	`Walker w = m.walker(hid);`	374	`Walker w = m.walker(hid);`
375	`FaceID f = w.face();`	375	`FaceID f = w.face();`
376		376
377	`if(f != InvalidFaceID){`	377	`if(f != InvalidFaceID){`
378	`VertexID v = w.vertex();`	378	`VertexID v = w.vertex();`
379	`Vec3d dir = m.pos(w.opp().vertex()) - m.pos(v);`	379	`Vec3d dir = m.pos(w.opp().vertex()) - m.pos(v);`
380	`Vec3d n = filtered_norms[f];`	380	`Vec3d n = filtered_norms[f];`
381	`vertex_positions[v] += m.pos(v) + 0.5 * n * dot(n, dir);`	381	`vertex_positions[v] += m.pos(v) + 0.5 * n * dot(n, dir);`
382	`count[v] += 1;`	382	`count[v] += 1;`
383	`}`	383	`}`
384		384
385	`}`	385	`}`
386	`double max_move= 0;`	386	`double max_move= 0;`
387	`for(VertexID v : m.vertices())`	387	`for(VertexID v : m.vertices())`
388	`{`	388	`{`
389	`Vec3d npos = vertex_positions[v] / double(count[v]);`	389	`Vec3d npos = vertex_positions[v] / double(count[v]);`
390	`double move = sqr_length(npos - m.pos(v));`	390	`double move = sqr_length(npos - m.pos(v));`
391	`if(move > max_move)`	391	`if(move > max_move)`
392	`max_move = move;`	392	`max_move = move;`
393		393
394	`m.pos(v) = npos;`	394	`m.pos(v) = npos;`
395	`}`	395	`}`
396	`if(max_move<sqr(1e-8*avg_len))`	396	`if(max_move<sqr(1e-8*avg_len))`
397	`{`	397	`{`
398	`cout << "iters " << sub_iter << endl;`	398	`cout << "iters " << sub_iter << endl;`
399	`break;`	399	`break;`
400	`}`	400	`}`
401	`}`	401	`}`
402	`}`	402	`}`
403	`}`	403	`}`
404		404
405		405
406	`}`	406	`}`
407		407

Subversion Repositories gelsvn

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