Subversion Repositories gelsvn

Rev

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

Rev 630 Rev 631
Line 2... Line 2...
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 <thread>
7 
#include "smooth.h"
 8 
#include "smooth.h"
8 
 
 9 
 
9 
#include <vector>
 10 
#include <vector>
10 
#include <algorithm>
 11 
#include <algorithm>
11 
#include "../CGLA/Mat3x3f.h"
 12 
#include "../CGLA/Mat3x3d.h"
12 
#include "../CGLA/Vec3d.h"
 13 
#include "../CGLA/Vec3d.h"
- 
 
 14 
#include "../CGLA/Quatd.h"
- 
 
 15 
#include "../Util/Timer.h"
13 
 
 16 
 
14 
#include "Manifold.h"
 17 
#include "Manifold.h"
15 
#include "AttributeVector.h"
 18 
#include "AttributeVector.h"
16 
 
 19 
 
17 
namespace HMesh
 20 
namespace HMesh
18 
{
 21 
{
19 
    using namespace std;
 22 
    using namespace std;
20 
    using namespace CGLA;
 23 
    using namespace CGLA;
- 
 
 24
21 
 
 25 
 
- 
 
 26
- 
 
 27
- 
 
 28 
    void for_each_vertex(Manifold& m, function<void(VertexID)> f) { for(auto v : m.vertices()) f(v); }
- 
 
 29 
    void for_each_face(Manifold& m, function<void(FaceID)> f) { for(auto p : m.faces()) f(p); }
- 
 
 30 
    void for_each_halfedge(Manifold& m, function<void(HalfEdgeID)> f) { for(auto h : m.halfedges()) f(h); }
- 
 
 31 
 
- 
 
 32
22 
    Vec3d laplacian(const Manifold& m, VertexID v)
 33 
    Vec3d laplacian(const Manifold& m, VertexID v)
23 
    {
 34 
    {
24 
        Vec3d avg_pos(0);
 35 
        Vec3d avg_pos(0);
25 
        int n = 0;
 36 
        int n = 0;
26 
 
 37
27 
        for(Walker w = m.walker(v); !w.full_circle(); w = w.circulate_vertex_cw()){
 38 
        circulate_vertex_ccw(m, v, [&](VertexID v){
28 
            avg_pos += m.pos(w.vertex());
 39 
            avg_pos += m.pos(v);
29 
            ++n;
 40 
            ++n;
30 
        }
 41 
        });
31 
        return avg_pos / n - m.pos(v);
 42 
        return avg_pos / n - m.pos(v);
32 
    }
 43 
    }
- 
 
 44
- 
 
 45
- 
 
 46 
     int CORES = 8;
33 
 
 47
34 
    void laplacian_smooth(Manifold& m, float t)
 48 
    void laplacian_smooth0(Manifold& m,float weight, int max_iter)
35 
    {
 49 
    {
- 
 
 50 
        for(int iter=0;iter<max_iter; ++iter) {
36 
        VertexAttributeVector<Vec3d> pos(m.allocated_vertices());
 51 
            VertexAttributeVector<Vec3d> L_attr(m.no_vertices());
37 
 
 52
- 
 
 53 
            for(VertexIDIterator v = m.vertices_begin(); v != m.vertices_end(); ++v)
- 
 
 54 
                if(!boundary(m, *v))
- 
 
 55 
                    L_attr[*v] =laplacian(m, *v);
38 
        for(VertexIDIterator v = m.vertices_begin(); v != m.vertices_end(); ++v){
 56 
            for(VertexIDIterator v = m.vertices_begin(); v != m.vertices_end(); ++v){
39 
            if(!boundary(m, *v))
 57 
                if(!boundary(m, *v))
- 
 
 58 
                    m.pos(*v) += weight*L_attr[*v];
- 
 
 59 
            }
- 
 
 60 
        }
- 
 
 61 
    }
- 
 
 62
- 
 
 63 
    void laplacian_smooth1(Manifold& m,float weight, int max_iter)
- 
 
 64 
    {
- 
 
 65 
        for(int iter=0;iter<max_iter; ++iter) {
- 
 
 66 
            VertexAttributeVector<Vec3d> L_attr(m.no_vertices());
- 
 
 67
- 
 
 68 
            for(VertexID v : m.vertices())
- 
 
 69 
                if(!boundary(m, v))
- 
 
 70 
                    L_attr[v] =laplacian(m, v);
- 
 
 71 
            for(VertexID v : m.vertices()){
- 
 
 72 
                if(!boundary(m, v))
- 
 
 73 
                    m.pos(v) += weight*L_attr[v];
- 
 
 74 
            }
- 
 
 75 
        }
- 
 
 76 
    }
- 
 
 77
- 
 
 78 
    void laplacian_smooth2(Manifold& m,float weight, int max_iter)
- 
 
 79 
    {
- 
 
 80 
        auto new_pos = m.positions_attribute_vector();
- 
 
 81 
        for(int iter=0;iter<max_iter; ++iter) {
- 
 
 82 
            for(auto v : m.vertices())
- 
 
 83 
                if(!boundary(m, v))
40 
                pos[*v] =  t * laplacian(m, *v) + m.pos(*v);
 84 
                    new_pos[v] = weight*laplacian(m, v)+m.pos(v);
- 
 
 85 
            m.positions_attribute_vector() = new_pos;
41 
        }
 86 
        }
- 
 
 87 
    }
- 
 
 88 
 
42 
 
 89 
 
- 
 
 90
- 
 
 91 
    void laplacian_smooth3(Manifold& m, float weight, int max_iter)
- 
 
 92 
    {
- 
 
 93 
        for(int iter=0;iter<max_iter; ++iter) {
43 
        for(VertexIDIterator v = m.vertices_begin(); v != m.vertices_end(); ++v){
 94 
            auto new_pos = m.positions_attribute_vector();
- 
 
 95 
            for(auto v : m.vertices())
44 
            if(!boundary(m, *v))
 96 
                if(!boundary(m, v))
45 
                m.pos(*v) = pos[*v];
 97 
                    new_pos[v] = weight*laplacian(m, v)+m.pos(v);
- 
 
 98 
            m.positions_attribute_vector() = move(new_pos);
46 
        }
 99 
        }
47 
    }
 100 
    }
48 
 
 101 
 
49 
    void taubin_smooth(Manifold& m, int max_iter)
 102 
    void laplacian_smooth4(Manifold& m,float weight, int max_iter)
50 
    {
 103 
    {
- 
 
 104 
        auto new_pos = m.positions_attribute_vector();
51 
        for(int iter = 0; iter < max_iter; ++iter) {
 105 
        for(int iter=0;iter<max_iter; ++iter) {
- 
 
 106 
            for(auto v : m.vertices())
- 
 
 107 
                if(!boundary(m, v))
52 
            VertexAttributeVector<Vec3d> lap(m.allocated_vertices());
 108 
                    new_pos[v] = weight*laplacian(m, v)+m.pos(v);
- 
 
 109 
            swap(m.positions_attribute_vector(),new_pos);
- 
 
 110 
        }
- 
 
 111 
    }
- 
 
 112
- 
 
 113 
    void laplacian_smooth4_5(Manifold& m,float weight, int max_iter)
- 
 
 114 
    {
- 
 
 115 
        auto new_pos = m.positions_attribute_vector();
- 
 
 116 
        for(int iter=0;iter<max_iter; ++iter) {
- 
 
 117 
            for_each_vertex(m, [&](VertexID v) {new_pos[v] = weight*laplacian(m, v)+m.pos(v);});
- 
 
 118 
            swap(m.positions_attribute_vector(),new_pos);
- 
 
 119 
        }
- 
 
 120 
    }
53 
 
 121 
 
54 
            for(VertexIDIterator v = m.vertices_begin(); v != m.vertices_end(); ++v){
 - 
 
55 
                if(!boundary(m, *v))
 - 
 
56 
                    lap[*v] =  laplacian(m, *v);
 - 
 
57 
            }
 - 
 
58 
 
 122 
 
- 
 
 123
- 
 
 124 
    void laplacian_smooth5(Manifold& m, float weight, int max_iter)
- 
 
 125 
    {
- 
 
 126 
        for(int iter=0;iter<max_iter; ++iter) {
59 
            for(VertexIDIterator v = m.vertices_begin(); v != m.vertices_end(); ++v){
 127 
            auto new_pos = m.positions_attribute_vector();
- 
 
 128 
            vector<thread> t_vec;
- 
 
 129 
            for(auto v : m.vertices())
60 
                if(!boundary(m, *v))
 130 
                if(!boundary(m, v))
- 
 
 131 
                    t_vec.push_back(thread([&](VertexID vid){
- 
 
 132 
                        if(!boundary(m, vid))
61 
                    m.pos(*v) += (iter%2 == 0) ? 0.5  * lap[*v] : -0.52 * lap[*v];
 133 
                            new_pos[vid] = weight*laplacian(m, vid)+ m.pos(vid);},v));
- 
 
 134 
            for(int i=0;i<t_vec.size();++i)
- 
 
 135 
                t_vec[i].join();
- 
 
 136 
            m.positions_attribute_vector() = move(new_pos);
62 
            }
 137 
        }
- 
 
 138 
    }
- 
 
 139
- 
 
 140
- 
 
 141
- 
 
 142 
    inline void laplacian_smooth_vertex(Manifold& m,vector<VertexID>& vids,
- 
 
 143 
                                        VertexAttributeVector<Vec3d>& new_pos,
- 
 
 144 
                                        float weight){
- 
 
 145 
        for(auto v: vids)
- 
 
 146 
            new_pos[v] = m.pos(v)+weight*laplacian(m, v);
- 
 
 147 
    }
- 
 
 148
- 
 
 149 
    void laplacian_smooth6(Manifold& m, float weight, int max_iter)
- 
 
 150 
    {
- 
 
 151 
        vector<vector<VertexID>> vertex_ids(CORES);
- 
 
 152 
        auto batch_size = m.no_vertices()/CORES;
- 
 
 153 
        int cnt = 0;
- 
 
 154 
        for_each_vertex(m, [&](VertexID v) {
- 
 
 155 
            if (!boundary(m, v))
- 
 
 156 
                vertex_ids[(cnt++/batch_size)%CORES].push_back(v);
- 
 
 157 
        });
- 
 
 158 
        vector<thread> t_vec(CORES);
- 
 
 159 
        VertexAttributeVector<Vec3d> new_pos = m.positions_attribute_vector();
- 
 
 160 
        for(int iter=0;iter<max_iter; ++iter) {
- 
 
 161 
            for(int thread_no=0;thread_no<CORES;++thread_no)
- 
 
 162 
                t_vec[thread_no] = thread(laplacian_smooth_vertex,
- 
 
 163 
                                          ref(m), ref(vertex_ids[thread_no]),
- 
 
 164 
                                          ref(new_pos), weight);
- 
 
 165 
            for(int thread_no=0;thread_no<CORES;++thread_no)
- 
 
 166 
                t_vec[thread_no].join();
- 
 
 167 
            swap(m.positions_attribute_vector(), new_pos);
63 
        }
 168 
        }
- 
 
 169
- 
 
 170 
    }
- 
 
 171
- 
 
 172 
    void laplacian_smooth(Manifold& m,float weight, int max_iter)
- 
 
 173 
    {
- 
 
 174 
        laplacian_smooth4(m,weight,max_iter);
- 
 
 175 
//        Util::Timer tim;
- 
 
 176 
//
- 
 
 177 
//        cout << "Method 0: ";
- 
 
 178 
//        for(int i=0;i<5;++i) {
- 
 
 179 
//            Manifold m2 = m;
- 
 
 180 
//            tim.start();
- 
 
 181 
//            laplacian_smooth0(m2, 1, 1000);
- 
 
 182 
//            cout << tim.get_secs() << ", ";
- 
 
 183 
//        }
- 
 
 184 
//        cout << endl;
- 
 
 185 
//
- 
 
 186 
//        cout << "Method 1: ";
- 
 
 187 
//        for(int i=0;i<5;++i) {
- 
 
 188 
//            Manifold m2 = m;
- 
 
 189 
//            tim.start();
- 
 
 190 
//            laplacian_smooth1(m2, 1, 1000);
- 
 
 191 
//            cout << tim.get_secs() << ", ";
- 
 
 192 
//        }
- 
 
 193 
//        cout << endl;
- 
 
 194 
//
- 
 
 195 
//        cout << "Method 2: ";
- 
 
 196 
//        for(int i=0;i<5;++i) {
- 
 
 197 
//            Manifold m2 = m;
- 
 
 198 
//            tim.start();
- 
 
 199 
//            laplacian_smooth2(m2, 1, 1000);
- 
 
 200 
//            cout << tim.get_secs() << ", ";
- 
 
 201 
//        }
- 
 
 202 
//        cout << endl;
- 
 
 203 
//
- 
 
 204 
//        cout << "Method 3: ";
- 
 
 205 
//        for(int i=0;i<5;++i) {
- 
 
 206 
//            Manifold m2 = m;
- 
 
 207 
//            tim.start();
- 
 
 208 
//            laplacian_smooth3(m2, 1, 1000);
- 
 
 209 
//            cout << tim.get_secs() << ", ";
- 
 
 210 
//        }
- 
 
 211 
//        cout << endl;
- 
 
 212 
//
- 
 
 213 
//        cout << "Method 4: ";
- 
 
 214 
//        for(int i=0;i<5;++i) {
- 
 
 215 
//            Manifold m2 = m;
- 
 
 216 
//            tim.start();
- 
 
 217 
//            laplacian_smooth4(m2, 1, 1000);
- 
 
 218 
//            cout << tim.get_secs() << ", ";
- 
 
 219 
//        }
- 
 
 220 
//        cout << endl;
- 
 
 221 
//
- 
 
 222 
////        cout << "Method 5: ";
- 
 
 223 
////        for(int i=0;i<3;++i) {
- 
 
 224 
////            Manifold m2 = m;
- 
 
 225 
////            tim.start();
- 
 
 226 
////            laplacian_smooth5(m2, 1, 1000);
- 
 
 227 
////            cout << " sec: "<< tim.get_secs() << ",";
- 
 
 228 
////        }
- 
 
 229 
////        cout << endl;
- 
 
 230 
//
- 
 
 231 
//        CORES = 2;
- 
 
 232 
//        cout << "Method 6 (2): ";
- 
 
 233 
//        for(int i=0;i<5;++i) {
- 
 
 234 
//            Manifold m2 = m;
- 
 
 235 
//            tim.start();
- 
 
 236 
//            laplacian_smooth6(m2, 1, 1000);
- 
 
 237 
//            cout << tim.get_secs() << ", ";
- 
 
 238 
//        }
- 
 
 239 
//        cout << endl;
- 
 
 240 
//
- 
 
 241 
//        CORES = 4;
- 
 
 242 
//        cout << "Method 6 (4): ";
- 
 
 243 
//        for(int i=0;i<5;++i) {
- 
 
 244 
//            Manifold m2 = m;
- 
 
 245 
//            tim.start();
- 
 
 246 
//            laplacian_smooth6(m2, 1, 1000);
- 
 
 247 
//            cout << tim.get_secs() << ", ";
- 
 
 248 
//        }
- 
 
 249 
//        cout << endl;
- 
 
 250 
//
- 
 
 251 
//        CORES = 8;
- 
 
 252 
//        cout << "Method 6 (8): ";
- 
 
 253 
//        for(int i=0;i<5;++i) {
- 
 
 254 
//            Manifold m2 = m;
- 
 
 255 
//            tim.start();
- 
 
 256 
//            laplacian_smooth6(m2, 1, 1000);
- 
 
 257 
//            cout << tim.get_secs() << ", ";
- 
 
 258 
//        }
- 
 
 259 
//        cout << endl;
- 
 
 260 
 
- 
 
 261 
 
64 
    }
 262 
    }
65 
 
 263 
 
- 
 
 264
- 
 
 265
- 
 
 266
- 
 
 267 
    void taubin_smooth(Manifold& m, int max_iter)
- 
 
 268 
    {
- 
 
 269 
        auto new_pos = m.positions_attribute_vector();
- 
 
 270 
        for(int iter = 0; iter < 2*max_iter; ++iter) {
- 
 
 271 
            for(VertexID v : m.vertices())
- 
 
 272 
                if(!boundary(m, v))
- 
 
 273 
                    new_pos[v] = (iter%2 == 0 ? 0.5 : -0.52) * laplacian(m, v) + m.pos(v);
- 
 
 274 
            swap(m.positions_attribute_vector(), new_pos);
- 
 
 275 
        }
- 
 
 276 
    }
- 
 
 277
66 
    void face_neighbourhood(Manifold& m, FaceAttributeVector<int>& touched, FaceID f, vector<Vec3d>& nbrs)
 278 
    void face_neighbourhood(Manifold& m, FaceID f, vector<FaceID>& nbrs)
67 
    {
279 
    {
- 
 
 280 
        FaceAttributeVector<int> touched(m.allocated_faces(), -1);
- 
 
 281 
        touched[f] = 1;
68 
        nbrs.push_back(normal(m, f));
 282 
        nbrs.push_back(f);
69 
        for(Walker wf = m.walker(f); !wf.full_circle(); wf = wf.circulate_face_cw()){
 283 
        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()){
 284 
            for(Walker wv = m.walker(wf.vertex()); !wv.full_circle(); wv = wv.circulate_vertex_cw()){
71 
                FaceID fn = wv.face();
 285 
                FaceID fn = wv.face();
72 
                if(fn != InvalidFaceID && touched[fn] != touched[f]){
 286 
                if(fn != InvalidFaceID && touched[fn] != touched[f]){
73 
                    Vec3d norm = normal(m, fn);
 - 
 
74 
                    if(!isnan(sqr_length(norm)))
 - 
 
75 
                        nbrs.push_back(norm);
 287 
                    nbrs.push_back(fn);
76 
                    else
 - 
 
77 
                        cout << "bad normal detected" << endl;
 - 
 
78 
                    touched[fn] = touched[f];
 288 
                    touched[fn] = 1;
79 
                }
 289 
                }
80 
            }
 290 
            }
81 
        }
 291 
        }
82 
    }
 292 
    }
83 
 
 293
84 
 
 294
85 
 
 295
86 
    Vec3d filtered_normal(Manifold& m,  FaceAttributeVector<int>& touched, FaceID f)
 296 
    Vec3d fvm_filtered_normal(Manifold& m, FaceID f)
87 
    {
 297 
    {
88 
        const float sigma = .1f;
 298 
        const float sigma = .1f;
89 
 
 299
90 
        vector<Vec3d> norms;
 300 
        vector<FaceID> nbrs;
91 
        face_neighbourhood(m, touched, f, norms);
 301 
        face_neighbourhood(m, f, nbrs);
92 
        float min_dist_sum=1e32f;
 302 
        float min_dist_sum=1e32f;
93 
        long int median=-1;
 303 
        long int median=-1;
94 
 
 304
- 
 
 305 
        vector<Vec3d> normals(nbrs.size());
95 
        for(size_t i=0;i<norms.size();++i)
 306 
        for(size_t i=0;i<nbrs.size();++i)
96 
        {
 307 
        {
- 
 
 308 
            normals[i] = normal(m, nbrs[i]);
97 
            float dist_sum = 0;
 309 
            float dist_sum = 0;
98 
            for(size_t j=0;j<norms.size(); ++j)
 310 
            for(size_t j=0;j<nbrs.size(); ++j)
99 
                dist_sum += 1.0f - dot(norms[i], norms[j]);
 311 
                dist_sum += 1.0f - dot(normals[i], normals[j]);
100 
            if(dist_sum < min_dist_sum)
 312 
            if(dist_sum < min_dist_sum)
101 
            {
 313 
            {
102 
                min_dist_sum = dist_sum;
 314 
                min_dist_sum = dist_sum;
103 
                median = i;
 315 
                median = i;
104 
            }
 316 
            }
105 
        }
 317 
        }
106 
        assert(median != -1);
 318 
        assert(median != -1);
107 
        Vec3d median_norm = norms[median];
 319 
        Vec3d median_norm = normals[median];
108 
        Vec3d avg_norm(0);
 320 
        Vec3d avg_norm(0);
109 
        for(size_t i=0;i<norms.size();++i)
 321 
        for(size_t i=0;i<nbrs.size();++i)
110 
        {
 322 
        {
111 
            float w = exp((dot(median_norm, norms[i])-1)/sigma);
 323 
            float w = exp((dot(median_norm, normals[i])-1)/sigma);
112 
            if(w<1e-2) w = 0;
 324 
            if(w<1e-2) w = 0;
113 
            avg_norm += w*norms[i];
 325 
            avg_norm += w*normals[i];
114 
        }
 326 
        }
115 
        return normalize(avg_norm);
 327 
        return normalize(avg_norm);
116 
    }
 328 
    }
- 
 
 329 
    Vec3d bilateral_filtered_normal(Manifold& m, FaceID f, double avg_len)
- 
 
 330 
    {
- 
 
 331 
        vector<FaceID> nbrs;
- 
 
 332 
        face_neighbourhood(m, f, nbrs);
- 
 
 333 
        Vec3d p0 = centre(m, f);
- 
 
 334 
        Vec3d n0 = normal(m, f);
- 
 
 335 
        vector<Vec3d> normals(nbrs.size());
- 
 
 336 
        vector<Vec3d> pos(nbrs.size());
- 
 
 337 
        Vec3d fn(0);
- 
 
 338 
        for(FaceID nbr : nbrs)
- 
 
 339 
        {
- 
 
 340 
            Vec3d n = normal(m, nbr);
- 
 
 341 
            Vec3d p = centre(m, nbr);
- 
 
 342 
            double w_a = exp(-acos(max(-1.0,min(1.0,dot(n,n0))))/(M_PI/32.0));
- 
 
 343 
            double w_s = exp(-length(p-p0)/avg_len);
- 
 
 344
- 
 
 345 
            fn += area(m, nbr)* w_a * w_s * n;
- 
 
 346
- 
 
 347 
        }
- 
 
 348 
        return normalize(fn);
- 
 
 349 
    }
117 
 
 350
118 
    void fvm_smooth(HMesh::Manifold& m, int max_iter)
 351 
    void anisotropic_smooth(HMesh::Manifold& m, int max_iter, NormalSmoothMethod nsm)
119 
    {
 352 
    {
- 
 
 353 
        double avg_len=0;
- 
 
 354 
        for(HalfEdgeID hid : m.halfedges())
- 
 
 355 
            avg_len += length(m, hid);
- 
 
 356 
        avg_len /= 2.0;
120 
        for(int iter = 0;iter<max_iter; ++iter)
 357 
        for(int iter = 0;iter<max_iter; ++iter)
121 
        {
 358 
        {
122 
            FaceAttributeVector<int> touched(m.allocated_faces(), -1);
 359
123 
            FaceAttributeVector<Vec3d> filtered_norms(m.allocated_faces());
 360 
            FaceAttributeVector<Vec3d> filtered_norms(m.allocated_faces());
124 
 
 - 
 
125 
            int i = 0;
 361
126 
            for(FaceIDIterator f = m.faces_begin(); f != m.faces_end(); ++f,++i){
 362 
            for(FaceID f: m.faces()){
- 
 
 363 
                filtered_norms[f] = (nsm == BILATERAL_NORMAL_SMOOTH)?
127 
                touched[*f] = i;
 364 
                bilateral_filtered_normal(m, f, avg_len):
128 
                filtered_norms[*f] = filtered_normal(m, touched, *f);
 365 
                fvm_filtered_normal(m, f);
129 
            }
 366 
            }
130 
 
 367
131 
            VertexAttributeVector<Vec3d> vertex_positions(m.allocated_vertices());
 368 
            VertexAttributeVector<Vec3d> vertex_positions(m.allocated_vertices(), Vec3d(0));
132 
 
 - 
 
- 
 
 369 
            VertexAttributeVector<int> count(m.allocated_vertices(), 0);
133 
            for(int iter=0;iter<20;++iter)
 370 
            for(int sub_iter=0;sub_iter<100;++sub_iter)
134 
            {
 371 
            {
135 
                for(VertexIDIterator v = m.vertices_begin(); v != m.vertices_end(); ++v){
 372 
                for(HalfEdgeID hid : m.halfedges())
136 
                    Vec3d move(0);
 373 
                {
137 
                    for(Walker w = m.walker(*v); !w.full_circle(); w = w.circulate_vertex_cw()){
 - 
 
138 
                        FaceID f1 = w.face();
 374 
                    Walker w = m.walker(hid);
139 
                        FaceID f2 = w.opp().face();
 375 
                    FaceID f = w.face();
140 
                        Vec3d dir = m.pos(w.vertex()) - m.pos(*v);
 376
141 
 
 - 
 
142 
                        if(f1 != InvalidFaceID){
 377 
                    if(f != InvalidFaceID){
143 
                            Vec3d n1 = filtered_norms[f1];
 378 
                        VertexID v = w.vertex();
144 
                            move += 0.05 * n1 * dot(n1, dir);
 379 
                        Vec3d dir = m.pos(w.opp().vertex()) - m.pos(v);
145 
                        }
 - 
 
146 
                        if(f2 != InvalidFaceID)
 - 
 
147 
                        {
 - 
 
148 
                            Vec3d n2 = filtered_norms[f2];
 380 
                        Vec3d n = filtered_norms[f];
149 
                            move += 0.05 * n2 * dot(n2, dir);
 381 
                        vertex_positions[v] += m.pos(v) + 0.5 * n * dot(n, dir);
150 
                        }
 382 
                        count[v] += 1;
151 
                    }
 383 
                    }
- 
 
 384
- 
 
 385 
                }
- 
 
 386 
                double max_move= 0;
- 
 
 387 
                for(VertexID v : m.vertices())
- 
 
 388 
                {
152 
                    vertex_positions[*v] = m.pos(*v) + move;
 389 
                    Vec3d npos = vertex_positions[v] / double(count[v]);
- 
 
 390 
                    double move = sqr_length(npos - m.pos(v));
- 
 
 391 
                    if(move > max_move)
- 
 
 392 
                        max_move = move;
- 
 
 393
- 
 
 394 
                    m.pos(v) = npos;
- 
 
 395 
                }
- 
 
 396 
                if(max_move<sqr(1e-8*avg_len))
- 
 
 397 
                {
- 
 
 398 
                    cout << "iters " << sub_iter << endl;
- 
 
 399 
                    break;
153 
                }
 400 
                }
154 
                for(VertexIDIterator v = m.vertices_begin(); v != m.vertices_end(); ++v)
 - 
 
155 
                    m.pos(*v) = vertex_positions[*v];
 - 
 
156 
            }
 401 
            }
157 
        }
 402 
        }
158 
    }
 403 
    }
159 
 
 404
160 
 
 405
161 
}
 406 
}