Subversion Repositories gelsvn

Rev

Rev 182 | Rev 229 | Go to most recent revision | Show entire file | Ignore whitespace | Details | Blame | Last modification | View Log | RSS feed

Rev 182 Rev 215
Line 1... Line 1...
1 
#include <algorithm>
 1 
#include <algorithm>
2 
#include "CGLA/Mat3x3f.h"
 2 
#include "CGLA/Mat3x3f.h"
3 
#include "smooth.h"
 3 
#include "smooth.h"
4 
#include "HMesh/VertexCirculator.h"
 4 
#include "HMesh/VertexCirculator.h"
- 
 
 5 
#include "HMesh/FaceCirculator.h"
5 
 
 6 
 
6 
 
 7 
 
7 
namespace HMesh
 8 
namespace HMesh
8 
{
 9 
{
9 
 
 10 
 
10 
	using namespace std;
 11 
		using namespace std;
11 
	using namespace CGLA;
 12 
		using namespace CGLA;
12 
	using namespace HMesh;
 13 
		using namespace HMesh;
13 
 
 14 
 
14 
	Vec3f laplacian(VertexIter vi)
 15 
		Vec3f laplacian(VertexIter vi)
15 
	{
 16 
		{
16 
		Vec3f avg_pos(0);
 17 
				Vec3f avg_pos(0);
17 
		int n=0;
 18 
				int n=0;
18
19
19 
		VertexCirculator vc(vi);
 20 
				VertexCirculator vc(vi);
20 
		for(;!vc.end();++vc)
 21 
				for(;!vc.end();++vc)
21 
			{
 22 
				{
22 
				avg_pos += vc.get_vertex()->pos;
 23 
						avg_pos += vc.get_vertex()->pos;
23 
				++n;
 24 
						++n;
24 
			}
 25 
				}
25 
		return avg_pos/n-vi->pos;
 26 
				return avg_pos/n-vi->pos;
26 
	}
 27 
		}
27 
 
 28 
 
28 
	void laplacian_smooth(Manifold& m, float t)
 29 
		void laplacian_smooth(Manifold& m, float t)
29 
	{
 - 
 
30 
		vector<Vec3f> pos(m.no_vertices());
 - 
 
31 
		int i=0;
 - 
 
32 
		for(VertexIter vi = m.vertices_begin();	vi != m.vertices_end(); ++vi,++i)
 - 
 
33 
			if(!is_boundary(vi))
 - 
 
34 
				pos[i] =  t*laplacian(vi) + vi->pos;
 - 
 
35 
		i=0;
 - 
 
36 
		for(VertexIter vi = m.vertices_begin();	vi != m.vertices_end(); ++vi,++i)
 - 
 
37 
			if(!is_boundary(vi))
 - 
 
38 
				vi->pos = pos[i];
 - 
 
39 
	}
 - 
 
40 
 
 - 
 
41 
	void taubin_smooth(Manifold& m, int max_iter)
 - 
 
42 
	{
 - 
 
43 
		for(int iter=0;iter<max_iter;++iter)
 - 
 
44 
			{
 30 
		{
45 
				vector<Vec3f> lap(m.no_vertices());
 31 
				vector<Vec3f> pos(m.no_vertices());
46 
				int i=0;
 32 
				int i=0;
47 
				for(VertexIter vi = m.vertices_begin();
33 
				for(VertexIter vi = m.vertices_begin();	vi != m.vertices_end(); ++vi,++i)
48 
						vi != m.vertices_end(); ++vi,++i)
 - 
 
49 
					if(!is_boundary(vi))
 34 
						if(!is_boundary(vi))
50 
							lap[i] =  laplacian(vi);
 35 
								pos[i] =  t*laplacian(vi) + vi->pos;
51 
				i=0;
 36 
				i=0;
- 
 
 37 
				for(VertexIter vi = m.vertices_begin();	vi != m.vertices_end(); ++vi,++i)
- 
 
 38 
						if(!is_boundary(vi))
- 
 
 39 
								vi->pos = pos[i];
- 
 
 40 
		}
- 
 
 41 
 
- 
 
 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;
52 
				for(VertexIter vi = m.vertices_begin();	vi != m.vertices_end();
53 
						for(VertexIter vi = m.vertices_begin();	vi != m.vertices_end();
53 
						++vi,++i)
 54 
								++vi,++i)
54 
					if(!is_boundary(vi))
 55 
								if(!is_boundary(vi))
55 
						vi->pos += (iter%2 == 0) ?
56 
										vi->pos += (iter%2 == 0) ?
56 
								 0.5  * lap[i] :
57 
												0.5  * lap[i] :
57 
								-0.52 * lap[i];
 58 
												-0.52 * lap[i];
- 
 
 59 
				}
- 
 
 60 
		}
- 
 
 61 
 
- 
 
 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 
		{
- 
 
 87 
				const float sigma = .1;
- 
 
 88 
 
- 
 
 89 
				vector<Vec3f> norms;
- 
 
 90 
				face_neighbourhood(f, norms);
- 
 
 91 
				float min_dist_sum=1e32;
- 
 
 92 
				int median=-1;
- 
 
 93
- 
 
 94 
				for(int i=0;i<norms.size();++i)
- 
 
 95 
				{
- 
 
 96 
						float dist_sum = 0;
- 
 
 97 
						for(int j=0;j<norms.size(); ++j)
- 
 
 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);
- 
 
 108 
				for(int i=0;i<norms.size();++i)
- 
 
 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];
58 
			}
 113 
				}
- 
 
 114 
				return normalize(avg_norm);
59 
	}
 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 
		}
60 
 
 168 
 
61 
 
 169 
 
62 
}
 170 
}