Subversion Repositories gelsvn

#include "smooth.h"

#include <vector>
#include <algorithm>

#include <CGLA/Mat3x3f.h>
#include <CGLA/Vec3f.h>

#include "Manifold.h"
#include "VertexHandle.h"
#include "FaceHandle.h"
#include "HalfEdgeHandle.h"
#include "VertexCirculator.h"
#include "FaceCirculator.h"


namespace HMesh
{
    using std::vector;
    using std::cout;
    using std::endl;
    using std::flush;
    using CGLA::Vec3f;
    using CGLA::isnan;

    Vec3f laplacian(const VertexHandle& vi)
    {
        Vec3f avg_pos(0);
        int n = 0;

        for(VertexCirculator vc(vi); !vc.end(); ++vc){
            avg_pos += vc.vertex().get().pos;
            ++n;
        }
        return avg_pos / n - vi.get().pos;
    }

    void laplacian_smooth(Manifold& m, float t)
    {
        vector<Vec3f> pos(m.no_vertices());
        int i=0;
        for(VertexHandle vi = m.vertices_begin(); vi != m.vertices_end(); ++vi, ++i){
            if(!is_boundary(vi))
                pos[i] =  t*laplacian(vi) + vi.get().pos;
        }
        i=0;
        for(VertexHandle vi = m.vertices_begin();       vi != m.vertices_end(); ++vi, ++i){
            if(!is_boundary(vi))
                vi.get().pos = pos[i];
        }
    }

    void taubin_smooth(Manifold& m, int max_iter)
    {
        for(int iter = 0; iter < max_iter; ++iter) {
            vector<Vec3f> lap(m.no_vertices());
            int i=0;
            for(VertexHandle vi = m.vertices_begin(); vi != m.vertices_end(); ++vi, ++i){
                if(!is_boundary(vi))
                    lap[i] =  laplacian(vi);
            }
            i=0;
            for(VertexHandle vi = m.vertices_begin(); vi != m.vertices_end(); ++vi, ++i){
                if(!is_boundary(vi))
                    vi.get().pos += (iter%2 == 0) ? 0.5  * lap[i] : -0.52 * lap[i];
            }
        }
    }

    void face_neighbourhood(const FaceHandle& f, vector<Vec3f>& nbrs)
    {   
        nbrs.push_back(normal(f));
        for(FaceCirculator fc(f); !fc.end(); ++fc){
            for(VertexCirculator vc(fc.vertex()); !vc.end(); ++vc){
                FaceHandle fn = vc.face();
                if(fn != NULL_FACE_HANDLE && fn.get().touched != f.get().touched){
                    Vec3f norm = normal(fn);
                    if(!isnan(sqr_length(norm)))
                        nbrs.push_back(norm);
                    else
                        cout << "bad normal detected" << endl;
                    fn.get().touched = f.get().touched;
                }
            }
        }
    }



    Vec3f filtered_normal(const FaceHandle& f)
    {
        const float sigma = .1f;

        vector<Vec3f> norms;
        face_neighbourhood(f, norms);
        float min_dist_sum=1e32f;
        int median=-1;

        for(unsigned int i=0;i<norms.size();++i)
        {
            float dist_sum = 0;
            for(unsigned int 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);
        Vec3f median_norm = norms[median];
        Vec3f avg_norm(0);
        for(unsigned int 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)
        {
            for(FaceHandle f = m.faces_begin(); f != m.faces_end(); ++f){
                f.get().touched = -1;
            }

            vector<Vec3f> filtered_norms(m.no_faces());
            int i=0;
            for(FaceHandle f = m.faces_begin(); f != m.faces_end(); ++f,++i){
                f.get().touched = i;
                filtered_norms[i] = filtered_normal(f);
            }

            cout << "Updating vertices" << flush;
            m.enumerate_faces();
            m.enumerate_vertices();

            vector<Vec3f> vertex_positions(m.no_vertices());

            for(int iter=0;iter<10;++iter)
            {
                for(VertexHandle vi = m.vertices_begin(); vi != m.vertices_end(); ++vi){
                    Vec3f move(0);
                    for(VertexCirculator vc(vi); !vc.end(); ++vc){
                        HalfEdgeHandle h = vc.halfedge();
                        FaceHandle f1 = h.face();
                        FaceHandle f2 = h.opp().face();
                        Vec3f dir = h.vert().get().pos - vi.get().pos;

                        if(f1 == NULL_FACE_HANDLE){
                            Vec3f n1 = filtered_norms[f1.get().touched];
                            move += 0.01 * n1 * dot(n1, dir);
                        }
                        if(f2 == NULL_FACE_HANDLE)
                        {
                            Vec3f n2 = filtered_norms[f2.get().touched];
                            move += 0.01 * n2 * dot(n2, dir);
                        }
                    }
                    vertex_positions[vi.get().touched] = vi.get().pos + move;
                }
                for(VertexHandle vi = m.vertices_begin(); vi!=m.vertices_end(); ++vi)
                    vi.get().pos = vertex_positions[vi.get().touched];
            }
        }
    }


}