Subversion Repositories gelsvn

Rev

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

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