Subversion Repositories gelsvn

/* ----------------------------------------------------------------------- *
 * This file is part of GEL, http://www.imm.dtu.dk/GEL
 * Copyright (C) the authors and DTU Informatics
 * For license and list of authors, see ../../doc/intro.pdf
 * ----------------------------------------------------------------------- */

#include "smooth.h"

#include <vector>
#include <algorithm>
#include "../CGLA/Mat3x3f.h"
#include "../CGLA/Vec3d.h"

#include "Manifold.h"
#include "AttributeVector.h"

namespace HMesh
{
    using namespace std;
    using namespace CGLA;

    Vec3d laplacian(const Manifold& m, VertexID v)
    {
        Vec3d avg_pos(0);
        int n = 0;

        for(Walker w = m.walker(v); !w.full_circle(); w = w.circulate_vertex_cw()){
            avg_pos += m.pos(w.vertex());
            ++n;
        }
        return avg_pos / n - m.pos(v);
    }

    void laplacian_smooth(Manifold& m, float t)
    {
        VertexAttributeVector<Vec3d> pos(m.allocated_vertices());

        for(VertexIDIterator v = m.vertices_begin(); v != m.vertices_end(); ++v){
            if(!boundary(m, *v))
                pos[*v] =  t * laplacian(m, *v) + m.pos(*v);
        }

        for(VertexIDIterator v = m.vertices_begin(); v != m.vertices_end(); ++v){
            if(!boundary(m, *v))
                m.pos(*v) = pos[*v];
        }
    }

    void taubin_smooth(Manifold& m, int max_iter)
    {
        for(int iter = 0; iter < max_iter; ++iter) {
            VertexAttributeVector<Vec3d> lap(m.allocated_vertices());

            for(VertexIDIterator v = m.vertices_begin(); v != m.vertices_end(); ++v){
                if(!boundary(m, *v))
                    lap[*v] =  laplacian(m, *v);
            }

            for(VertexIDIterator v = m.vertices_begin(); v != m.vertices_end(); ++v){
                if(!boundary(m, *v))
                    m.pos(*v) += (iter%2 == 0) ? 0.5  * lap[*v] : -0.52 * lap[*v];
            }
        }
    }

    void face_neighbourhood(Manifold& m, FaceAttributeVector<int>& touched, FaceID f, vector<Vec3d>& nbrs)
    {   
        nbrs.push_back(normal(m, f));
        for(Walker wf = m.walker(f); !wf.full_circle(); wf = wf.circulate_face_cw()){
            for(Walker wv = m.walker(wf.vertex()); !wv.full_circle(); wv = wv.circulate_vertex_cw()){
                FaceID fn = wv.face();
                if(fn != InvalidFaceID && touched[fn] != touched[f]){
                    Vec3d norm = normal(m, fn);
                    if(!isnan(sqr_length(norm)))
                        nbrs.push_back(norm);
                    else
                        cout << "bad normal detected" << endl;
                    touched[fn] = touched[f];
                }
            }
        }
    }



    Vec3d filtered_normal(Manifold& m,  FaceAttributeVector<int>& touched, FaceID f)
    {
        const float sigma = .1f;

        vector<Vec3d> norms;
        face_neighbourhood(m, touched, f, norms);
        float min_dist_sum=1e32f;
        long int median=-1;

        for(size_t i=0;i<norms.size();++i)
        {
            float dist_sum = 0;
            for(size_t j=0;j<norms.size(); ++j)
                dist_sum += 1.0f - dot(norms[i], norms[j]);
            if(dist_sum < min_dist_sum)
            {
                min_dist_sum = dist_sum;
                median = i;
            }
        }
        assert(median != -1);
        Vec3d median_norm = norms[median];
        Vec3d avg_norm(0);
        for(size_t i=0;i<norms.size();++i)
        {
            float w = exp((dot(median_norm, norms[i])-1)/sigma);
            if(w<1e-2) w = 0;
            avg_norm += w*norms[i];
        }
        return normalize(avg_norm);
    }

    void fvm_smooth(HMesh::Manifold& m, int max_iter)
    {
        for(int iter = 0;iter<max_iter; ++iter)
        {
            FaceAttributeVector<int> touched(m.allocated_faces(), -1);
            FaceAttributeVector<Vec3d> filtered_norms(m.allocated_faces());

            int i = 0;
            for(FaceIDIterator f = m.faces_begin(); f != m.faces_end(); ++f,++i){
                touched[*f] = i;
                filtered_norms[*f] = filtered_normal(m, touched, *f);
            }

            VertexAttributeVector<Vec3d> vertex_positions(m.allocated_vertices());

            for(int iter=0;iter<20;++iter)
            {
                for(VertexIDIterator v = m.vertices_begin(); v != m.vertices_end(); ++v){
                    Vec3d move(0);
                    for(Walker w = m.walker(*v); !w.full_circle(); w = w.circulate_vertex_cw()){
                        FaceID f1 = w.face();
                        FaceID f2 = w.opp().face();
                        Vec3d dir = m.pos(w.vertex()) - m.pos(*v);

                        if(f1 != InvalidFaceID){
                            Vec3d n1 = filtered_norms[f1];
                            move += 0.05 * n1 * dot(n1, dir);
                        }
                        if(f2 != InvalidFaceID)
                        {
                            Vec3d n2 = filtered_norms[f2];
                            move += 0.05 * n2 * dot(n2, dir);
                        }
                    }
                    vertex_positions[*v] = m.pos(*v) + move;
                }
                for(VertexIDIterator v = m.vertices_begin(); v != m.vertices_end(); ++v)
                    m.pos(*v) = vertex_positions[*v];
            }
        }
    }


}