WebSVN – gelsvn – Diff – /trunk/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 "smooth.h"
 
#include <vector>
#include <algorithm>
#include "../CGLA/Mat3x3f.h"
#include "../CGLA/Vec3d.h"
 
#include "Manifold.h"
#include "AttributeVector.h"
 
namespace HMesh
{
    using namespace std;
    using namespace CGLA;
 
    Vec3d laplacian(const Manifold& m, VertexID v)
    {
        Vec3d avg_pos(0);
        int n = 0;
 
        for(Walker w = m.walker(v); !w.full_circle(); w = w.circulate_vertex_cw()){
            avg_pos += m.pos(w.vertex());
            ++n;
        }
        return avg_pos / n - m.pos(v);
    }
 
    void laplacian_smooth(Manifold& m, float t)
    {
        VertexAttributeVector<Vec3d> pos(m.allocated_vertices());
 
        for(VertexIDIterator v = m.vertices_begin(); v != m.vertices_end(); ++v){
            if(!boundary(m, *v))
                pos[*v] =  t * laplacian(m, *v) + m.pos(*v);
        }
 
        for(VertexIDIterator v = m.vertices_begin(); v != m.vertices_end(); ++v){
            if(!boundary(m, *v))
                m.pos(*v) = pos[*v];
        }
    }
 
    void taubin_smooth(Manifold& m, int max_iter)
    {
        for(int iter = 0; iter < max_iter; ++iter) {
            VertexAttributeVector<Vec3d> lap(m.allocated_vertices());
 
            for(VertexIDIterator v = m.vertices_begin(); v != m.vertices_end(); ++v){
                if(!boundary(m, *v))
                    lap[*v] =  laplacian(m, *v);
            }
 
            for(VertexIDIterator v = m.vertices_begin(); v != m.vertices_end(); ++v){
                if(!boundary(m, *v))
                    m.pos(*v) += (iter%2 == 0) ? 0.5  * lap[*v] : -0.52 * lap[*v];
            }
        }
    }
 
    void face_neighbourhood(Manifold& m, FaceAttributeVector<int>& touched, FaceID f, vector<Vec3d>& nbrs)
    {	
        nbrs.push_back(normal(m, 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]){
                    Vec3d norm = normal(m, fn);
                    if(!isnan(sqr_length(norm)))
                        nbrs.push_back(norm);
                    else
                        cout << "bad normal detected" << endl;
                    touched[fn] = touched[f];
                }
            }
        }
    }
 
 
 
    Vec3d filtered_normal(Manifold& m,  FaceAttributeVector<int>& touched, FaceID f)
    {
        const float sigma = .1f;
 
        vector<Vec3d> norms;
        face_neighbourhood(m, touched, f, norms);
        float min_dist_sum=1e32f;
        long int median=-1;
 
        for(size_t i=0;i<norms.size();++i)
        {
            float dist_sum = 0;
            for(size_t j=0;j<norms.size(); ++j)
                dist_sum += 1.0f - dot(norms[i], norms[j]);
            if(dist_sum < min_dist_sum)
            {
                min_dist_sum = dist_sum;
                median = i;
            }
        }
        assert(median != -1);
        Vec3d median_norm = norms[median];
        Vec3d avg_norm(0);
        for(size_t i=0;i<norms.size();++i)
        {
            float w = exp((dot(median_norm, norms[i])-1)/sigma);
            if(w<1e-2) w = 0;
            avg_norm += w*norms[i];
        }
        return normalize(avg_norm);
    }
 
    void fvm_smooth(HMesh::Manifold& m, int max_iter)
    {
        for(int iter = 0;iter<max_iter; ++iter)
        {
            FaceAttributeVector<int> touched(m.allocated_faces(), -1);
            FaceAttributeVector<Vec3d> filtered_norms(m.allocated_faces());
 
            int i = 0;
            for(FaceIDIterator f = m.faces_begin(); f != m.faces_end(); ++f,++i){
                touched[*f] = i;
                filtered_norms[*f] = filtered_normal(m, touched, *f);
            }
 
            VertexAttributeVector<Vec3d> vertex_positions(m.allocated_vertices());
 
            for(int iter=0;iter<20;++iter)
            {
                for(VertexIDIterator v = m.vertices_begin(); v != m.vertices_end(); ++v){
                    Vec3d move(0);
                    for(Walker w = m.walker(*v); !w.full_circle(); w = w.circulate_vertex_cw()){
                        FaceID f1 = w.face();
                        FaceID f2 = w.opp().face();
                        Vec3d dir = m.pos(w.vertex()) - m.pos(*v);
 
                        if(f1 != InvalidFaceID){
                            Vec3d n1 = filtered_norms[f1];
                            move += 0.05 * n1 * dot(n1, dir);
                        }
                        if(f2 != InvalidFaceID)
                        {
                            Vec3d n2 = filtered_norms[f2];
                            move += 0.05 * n2 * dot(n2, dir);
                        }
                    }
                    vertex_positions[*v] = m.pos(*v) + move;
                }
                for(VertexIDIterator v = m.vertices_begin(); v != m.vertices_end(); ++v)
                    m.pos(*v) = vertex_positions[*v];
            }
        }
    }
 
 
}
 

Rev 595	Rev 601
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 "smooth.h"`	7	`#include "smooth.h"`
8		8
9	`#include <vector>`	9	`#include <vector>`
10	`#include <algorithm>`	10	`#include <algorithm>`
11	`#include <CGLA/Mat3x3f.h>`	11	`#include "../CGLA/Mat3x3f.h"`
12	`#include <CGLA/Vec3d.h>`	12	`#include "../CGLA/Vec3d.h"`
13		13
14	`#include "Manifold.h"`	14	`#include "Manifold.h"`
15	`#include "AttributeVector.h"`	15	`#include "AttributeVector.h"`
16		16
17	`namespace HMesh`	17	`namespace HMesh`
18	`{`	18	`{`
19	`using namespace std;`	19	`using namespace std;`
20	`using namespace CGLA;`	20	`using namespace CGLA;`
21		21
22	`Vec3d laplacian(const Manifold& m, VertexID v)`	22	`Vec3d laplacian(const Manifold& m, VertexID v)`
23	`{`	23	`{`
24	`Vec3d avg_pos(0);`	24	`Vec3d avg_pos(0);`
25	`int n = 0;`	25	`int n = 0;`
26		26
27	`for(Walker w = m.walker(v); !w.full_circle(); w = w.circulate_vertex_cw()){`	27	`for(Walker w = m.walker(v); !w.full_circle(); w = w.circulate_vertex_cw()){`
28	`avg_pos += m.pos(w.vertex());`	28	`avg_pos += m.pos(w.vertex());`
29	`++n;`	29	`++n;`
30	`}`	30	`}`
31	`return avg_pos / n - m.pos(v);`	31	`return avg_pos / n - m.pos(v);`
32	`}`	32	`}`
33		33
34	`void laplacian_smooth(Manifold& m, float t)`	34	`void laplacian_smooth(Manifold& m, float t)`
35	`{`	35	`{`
36	`VertexAttributeVector<Vec3d> pos(m.allocated_vertices());`	36	`VertexAttributeVector<Vec3d> pos(m.allocated_vertices());`
37		37
38	`for(VertexIDIterator v = m.vertices_begin(); v != m.vertices_end(); ++v){`	38	`for(VertexIDIterator v = m.vertices_begin(); v != m.vertices_end(); ++v){`
39	`if(!boundary(m, *v))`	39	`if(!boundary(m, *v))`
40	`pos[v] = t laplacian(m, v) + m.pos(v);`	40	`pos[v] = t laplacian(m, v) + m.pos(v);`
41	`}`	41	`}`
42		42
43	`for(VertexIDIterator v = m.vertices_begin(); v != m.vertices_end(); ++v){`	43	`for(VertexIDIterator v = m.vertices_begin(); v != m.vertices_end(); ++v){`
44	`if(!boundary(m, *v))`	44	`if(!boundary(m, *v))`
45	`m.pos(v) = pos[v];`	45	`m.pos(v) = pos[v];`
46	`}`	46	`}`
47	`}`	47	`}`
48		48
49	`void taubin_smooth(Manifold& m, int max_iter)`	49	`void taubin_smooth(Manifold& m, int max_iter)`
50	`{`	50	`{`
51	`for(int iter = 0; iter < max_iter; ++iter) {`	51	`for(int iter = 0; iter < max_iter; ++iter) {`
52	`VertexAttributeVector<Vec3d> lap(m.allocated_vertices());`	52	`VertexAttributeVector<Vec3d> lap(m.allocated_vertices());`
53		53
54	`for(VertexIDIterator v = m.vertices_begin(); v != m.vertices_end(); ++v){`	54	`for(VertexIDIterator v = m.vertices_begin(); v != m.vertices_end(); ++v){`
55	`if(!boundary(m, *v))`	55	`if(!boundary(m, *v))`
56	`lap[v] = laplacian(m, v);`	56	`lap[v] = laplacian(m, v);`
57	`}`	57	`}`
58		58
59	`for(VertexIDIterator v = m.vertices_begin(); v != m.vertices_end(); ++v){`	59	`for(VertexIDIterator v = m.vertices_begin(); v != m.vertices_end(); ++v){`
60	`if(!boundary(m, *v))`	60	`if(!boundary(m, *v))`
61	`m.pos(v) += (iter%2 == 0) ? 0.5 lap[v] : -0.52 lap[*v];`	61	`m.pos(v) += (iter%2 == 0) ? 0.5 lap[v] : -0.52 lap[*v];`
62	`}`	62	`}`
63	`}`	63	`}`
64	`}`	64	`}`
65		65
66	`void face_neighbourhood(Manifold& m, FaceAttributeVector<int>& touched, FaceID f, vector<Vec3d>& nbrs)`	66	`void face_neighbourhood(Manifold& m, FaceAttributeVector<int>& touched, FaceID f, vector<Vec3d>& nbrs)`
67	`{`	67	`{`
68	`nbrs.push_back(normal(m, f));`	68	`nbrs.push_back(normal(m, f));`
69	`for(Walker wf = m.walker(f); !wf.full_circle(); wf = wf.circulate_face_cw()){`	69	`for(Walker wf = m.walker(f); !wf.full_circle(); wf = wf.circulate_face_cw()){`
70	`for(Walker wv = m.walker(wf.vertex()); !wv.full_circle(); wv = wv.circulate_vertex_cw()){`	70	`for(Walker wv = m.walker(wf.vertex()); !wv.full_circle(); wv = wv.circulate_vertex_cw()){`
71	`FaceID fn = wv.face();`	71	`FaceID fn = wv.face();`
72	`if(fn != InvalidFaceID && touched[fn] != touched[f]){`	72	`if(fn != InvalidFaceID && touched[fn] != touched[f]){`
73	`Vec3d norm = normal(m, fn);`	73	`Vec3d norm = normal(m, fn);`
74	`if(!isnan(sqr_length(norm)))`	74	`if(!isnan(sqr_length(norm)))`
75	`nbrs.push_back(norm);`	75	`nbrs.push_back(norm);`
76	`else`	76	`else`
77	`cout << "bad normal detected" << endl;`	77	`cout << "bad normal detected" << endl;`
78	`touched[fn] = touched[f];`	78	`touched[fn] = touched[f];`
79	`}`	79	`}`
80	`}`	80	`}`
81	`}`	81	`}`
82	`}`	82	`}`
83		83
84		84
85		85
86	`Vec3d filtered_normal(Manifold& m, FaceAttributeVector<int>& touched, FaceID f)`	86	`Vec3d filtered_normal(Manifold& m, FaceAttributeVector<int>& touched, FaceID f)`
87	`{`	87	`{`
88	`const float sigma = .1f;`	88	`const float sigma = .1f;`
89		89
90	`vector<Vec3d> norms;`	90	`vector<Vec3d> norms;`
91	`face_neighbourhood(m, touched, f, norms);`	91	`face_neighbourhood(m, touched, f, norms);`
92	`float min_dist_sum=1e32f;`	92	`float min_dist_sum=1e32f;`
93	`long int median=-1;`	93	`long int median=-1;`
94		94
95	`for(size_t i=0;i<norms.size();++i)`	95	`for(size_t i=0;i<norms.size();++i)`
96	`{`	96	`{`
97	`float dist_sum = 0;`	97	`float dist_sum = 0;`
98	`for(size_t j=0;j<norms.size(); ++j)`	98	`for(size_t j=0;j<norms.size(); ++j)`
99	`dist_sum += 1.0f - dot(norms[i], norms[j]);`	99	`dist_sum += 1.0f - dot(norms[i], norms[j]);`
100	`if(dist_sum < min_dist_sum)`	100	`if(dist_sum < min_dist_sum)`
101	`{`	101	`{`
102	`min_dist_sum = dist_sum;`	102	`min_dist_sum = dist_sum;`
103	`median = i;`	103	`median = i;`
104	`}`	104	`}`
105	`}`	105	`}`
106	`assert(median != -1);`	106	`assert(median != -1);`
107	`Vec3d median_norm = norms[median];`	107	`Vec3d median_norm = norms[median];`
108	`Vec3d avg_norm(0);`	108	`Vec3d avg_norm(0);`
109	`for(size_t i=0;i<norms.size();++i)`	109	`for(size_t i=0;i<norms.size();++i)`
110	`{`	110	`{`
111	`float w = exp((dot(median_norm, norms[i])-1)/sigma);`	111	`float w = exp((dot(median_norm, norms[i])-1)/sigma);`
112	`if(w<1e-2) w = 0;`	112	`if(w<1e-2) w = 0;`
113	`avg_norm += w*norms[i];`	113	`avg_norm += w*norms[i];`
114	`}`	114	`}`
115	`return normalize(avg_norm);`	115	`return normalize(avg_norm);`
116	`}`	116	`}`
117		117
118	`void fvm_smooth(HMesh::Manifold& m, int max_iter)`	118	`void fvm_smooth(HMesh::Manifold& m, int max_iter)`
119	`{`	119	`{`
120	`for(int iter = 0;iter<max_iter; ++iter)`	120	`for(int iter = 0;iter<max_iter; ++iter)`
121	`{`	121	`{`
122	`FaceAttributeVector<int> touched(m.allocated_faces(), -1);`	122	`FaceAttributeVector<int> touched(m.allocated_faces(), -1);`
123	`FaceAttributeVector<Vec3d> filtered_norms(m.allocated_faces());`	123	`FaceAttributeVector<Vec3d> filtered_norms(m.allocated_faces());`
124		124
125	`int i = 0;`	125	`int i = 0;`
126	`for(FaceIDIterator f = m.faces_begin(); f != m.faces_end(); ++f,++i){`	126	`for(FaceIDIterator f = m.faces_begin(); f != m.faces_end(); ++f,++i){`
127	`touched[*f] = i;`	127	`touched[*f] = i;`
128	`filtered_norms[f] = filtered_normal(m, touched, f);`	128	`filtered_norms[f] = filtered_normal(m, touched, f);`
129	`}`	129	`}`
130		130
131	`VertexAttributeVector<Vec3d> vertex_positions(m.allocated_vertices());`	131	`VertexAttributeVector<Vec3d> vertex_positions(m.allocated_vertices());`
132		132
133	`for(int iter=0;iter<20;++iter)`	133	`for(int iter=0;iter<20;++iter)`
134	`{`	134	`{`
135	`for(VertexIDIterator v = m.vertices_begin(); v != m.vertices_end(); ++v){`	135	`for(VertexIDIterator v = m.vertices_begin(); v != m.vertices_end(); ++v){`
136	`Vec3d move(0);`	136	`Vec3d move(0);`
137	`for(Walker w = m.walker(*v); !w.full_circle(); w = w.circulate_vertex_cw()){`	137	`for(Walker w = m.walker(*v); !w.full_circle(); w = w.circulate_vertex_cw()){`
138	`FaceID f1 = w.face();`	138	`FaceID f1 = w.face();`
139	`FaceID f2 = w.opp().face();`	139	`FaceID f2 = w.opp().face();`
140	`Vec3d dir = m.pos(w.vertex()) - m.pos(*v);`	140	`Vec3d dir = m.pos(w.vertex()) - m.pos(*v);`
141		141
142	`if(f1 != InvalidFaceID){`	142	`if(f1 != InvalidFaceID){`
143	`Vec3d n1 = filtered_norms[f1];`	143	`Vec3d n1 = filtered_norms[f1];`
144	`move += 0.05 * n1 * dot(n1, dir);`	144	`move += 0.05 * n1 * dot(n1, dir);`
145	`}`	145	`}`
146	`if(f2 != InvalidFaceID)`	146	`if(f2 != InvalidFaceID)`
147	`{`	147	`{`
148	`Vec3d n2 = filtered_norms[f2];`	148	`Vec3d n2 = filtered_norms[f2];`
149	`move += 0.05 * n2 * dot(n2, dir);`	149	`move += 0.05 * n2 * dot(n2, dir);`
150	`}`	150	`}`
151	`}`	151	`}`
152	`vertex_positions[v] = m.pos(v) + move;`	152	`vertex_positions[v] = m.pos(v) + move;`
153	`}`	153	`}`
154	`for(VertexIDIterator v = m.vertices_begin(); v != m.vertices_end(); ++v)`	154	`for(VertexIDIterator v = m.vertices_begin(); v != m.vertices_end(); ++v)`
155	`m.pos(v) = vertex_positions[v];`	155	`m.pos(v) = vertex_positions[v];`
156	`}`	156	`}`
157	`}`	157	`}`
158	`}`	158	`}`
159		159
160		160
161	`}`	161	`}`
162		162

Subversion Repositories gelsvn

(root)/trunk/src/HMesh/smooth.cpp – Rev 595 → 601