Subversion Repositories gelsvn

Rev

Rev 229 | Go to most recent revision | Only display areas with differences | Ignore whitespace | Details | Blame | Last modification | View Log | RSS feed

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