Subversion Repositories gelsvn

Rev

Rev 229 | Details | Compare with Previous | Last modification | View Log | RSS feed

Rev Author Line No. Line
149 jab 1 
#include <algorithm>
2 
#include "CGLA/Mat3x3f.h"
3 
#include "smooth.h"
4 
#include "HMesh/VertexCirculator.h"
215 jab 5 
#include "HMesh/FaceCirculator.h"
149 jab 6 
 
7 
 
150 jab 8 
namespace HMesh
149 jab 9 
{
10 
 
215 jab 11 
		using namespace std;
12 
		using namespace CGLA;
13 
		using namespace HMesh;
149 jab 14 
 
215 jab 15 
		Vec3f laplacian(VertexIter vi)
16 
		{
17 
				Vec3f avg_pos(0);
18 
				int n=0;
149 jab 19 
 
215 jab 20 
				VertexCirculator vc(vi);
21 
				for(;!vc.end();++vc)
22 
				{
23 
						avg_pos += vc.get_vertex()->pos;
24 
						++n;
25 
				}
26 
				return avg_pos/n-vi->pos;
27 
		}
149 jab 28 
 
215 jab 29 
		void laplacian_smooth(Manifold& m, float t)
30 
		{
31 
				vector<Vec3f> pos(m.no_vertices());
149 jab 32 
				int i=0;
215 jab 33 
				for(VertexIter vi = m.vertices_begin();	vi != m.vertices_end(); ++vi,++i)
34 
						if(!is_boundary(vi))
35 
								pos[i] =  t*laplacian(vi) + vi->pos;
149 jab 36 
				i=0;
215 jab 37 
				for(VertexIter vi = m.vertices_begin();	vi != m.vertices_end(); ++vi,++i)
38 
						if(!is_boundary(vi))
39 
								vi->pos = pos[i];
40 
		}
149 jab 41 
 
215 jab 42 
		void taubin_smooth(Manifold& m, int max_iter)
43 
		{
44 
				for(int iter=0;iter<max_iter;++iter)
45 
				{
46 
						vector<Vec3f> lap(m.no_vertices());
47 
						int i=0;
48 
						for(VertexIter vi = m.vertices_begin();
49 
								vi != m.vertices_end(); ++vi,++i)
50 
								if(!is_boundary(vi))
51 
										lap[i] =  laplacian(vi);
52 
						i=0;
53 
						for(VertexIter vi = m.vertices_begin();	vi != m.vertices_end();
54 
								++vi,++i)
55 
								if(!is_boundary(vi))
56 
										vi->pos += (iter%2 == 0) ?
57 
												0.5  * lap[i] :
58 
												-0.52 * lap[i];
59 
				}
60 
		}
149 jab 61 
 
215 jab 62 
		void face_neighbourhood(FaceIter f, vector<Vec3f>& nbrs)
63 
		{
64 
				nbrs.push_back(normal(f));
65 
				for(FaceCirculator fc(f); !fc.end(); ++fc)
66 
				{
67 
						for(VertexCirculator vc(fc.get_vertex()); !vc.end(); ++vc)
68 
						{
69 
								FaceIter fn = vc.get_face();
70 
								if(fn != NULL_FACE_ITER && fn->touched != f->touched)
71 
								{
72 
										Vec3f norm = normal(fn);
73 
										if(!isnan(sqr_length(norm)))
74 
												nbrs.push_back(norm);
75 
										else
76 
												cout << "bad normal detected" << endl;
77 
										fn->touched = f->touched;
78 
								}
79 
						}
80 
				}
81 
		}
82 
 
83 
 
84 
 
85 
		Vec3f filtered_normal(FaceIter f)
86 
		{
229 bj 87 
				const float sigma = .1f;
215 jab 88 
 
89 
				vector<Vec3f> norms;
90 
				face_neighbourhood(f, norms);
229 bj 91 
				float min_dist_sum=1e32f;
215 jab 92 
				int median=-1;
93 
 
229 bj 94 
				for(unsigned int i=0;i<norms.size();++i)
215 jab 95 
				{
96 
						float dist_sum = 0;
229 bj 97 
						for(unsigned int j=0;j<norms.size(); ++j)
215 jab 98 
								dist_sum += 1.0f - dot(norms[i], norms[j]);
99 
						if(dist_sum < min_dist_sum)
100 
						{
101 
								min_dist_sum = dist_sum;
102 
								median = i;
103 
						}
104 
				}
105 
				assert(median != -1);
106 
				Vec3f median_norm = norms[median];
107 
				Vec3f avg_norm(0);
229 bj 108 
				for(unsigned int i=0;i<norms.size();++i)
215 jab 109 
				{
110 
						float w = exp((dot(median_norm, norms[i])-1)/sigma);
111 
						if(w<1e-2) w = 0;
112 
						avg_norm += w*norms[i];
113 
				}
114 
				return normalize(avg_norm);
115 
		}
116 
 
117 
		void fvm_smooth(HMesh::Manifold& m, int max_iter)
118 
		{
119 
				for(int iter = 0;iter<max_iter; ++iter)
120 
				{
121 
						for(FaceIter f = m.faces_begin(); f != m.faces_end(); ++f)
122 
								f->touched = -1;
123 
 
124 
						vector<Vec3f> filtered_norms(m.no_faces());
125 
						int i=0;
126 
						for(FaceIter f = m.faces_begin(); f != m.faces_end(); ++f,++i)
127 
						{
128 
								f->touched = i;
129 
								filtered_norms[i] = filtered_normal(f);
130 
						}
131 
 
132 
						cout << "Updating vertices" << flush;
133 
						m.enumerate_faces();
134 
						m.enumerate_vertices();
135 
 
136 
						vector<Vec3f> vertex_positions(m.no_vertices());
137 
 
138 
						for(int iter=0;iter<10;++iter)
139 
						{
140 
								for(VertexIter vi=m.vertices_begin(); vi!=m.vertices_end(); ++vi)
141 
								{
142 
										Vec3f move(0);
143 
										for(VertexCirculator vc(vi); !vc.end(); ++vc)
144 
										{
145 
												HalfEdgeIter h = vc.get_halfedge();
146 
												FaceIter f1 = h->face;
147 
												FaceIter f2 = h->opp->face;
148 
												Vec3f dir = h->vert->pos - vi->pos;
149 
 
150 
												if(f1 != NULL_FACE_ITER)
151 
												{
152 
														Vec3f n1 = filtered_norms[f1->touched];
153 
														move += 0.01 * n1 * dot(n1, dir);
154 
												}
155 
												if(f2 != NULL_FACE_ITER)
156 
												{
157 
														Vec3f n2 = filtered_norms[f2->touched];
158 
														move += 0.01 * n2 * dot(n2, dir);
159 
												}
160 
										}
161 
										vertex_positions[vi->touched] = vi->pos + move;
162 
								}
163 
								for(VertexIter vi=m.vertices_begin(); vi!=m.vertices_end(); ++vi)
164 
										vi->pos = vertex_positions[vi->touched];
165 
						}
166 
				}
167 
		}
168 
 
169 
 
149 jab 170 
}