Subversion Repositories gelsvn

#include <algorithm>
#include "CGLA/Mat3x3f.h"
#include "smooth.h"
#include "HMesh/VertexCirculator.h"
#include "HMesh/FaceCirculator.h"


namespace HMesh
{

                using namespace std;
                using namespace CGLA;
                using namespace HMesh;

                Vec3f laplacian(VertexIter vi)
                {
                                Vec3f avg_pos(0);
                                int n=0;
                
                                VertexCirculator vc(vi);
                                for(;!vc.end();++vc)
                                {
                                                avg_pos += vc.get_vertex()->pos;
                                                ++n;
                                }
                                return avg_pos/n-vi->pos;
                }

                void laplacian_smooth(Manifold& m, float t)
                {
                                vector<Vec3f> pos(m.no_vertices());
                                int i=0;
                                for(VertexIter vi = m.vertices_begin(); vi != m.vertices_end(); ++vi,++i)
                                                if(!is_boundary(vi))
                                                                pos[i] =  t*laplacian(vi) + vi->pos;
                                i=0;
                                for(VertexIter vi = m.vertices_begin(); vi != m.vertices_end(); ++vi,++i)
                                                if(!is_boundary(vi))
                                                                vi->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(VertexIter vi = m.vertices_begin(); 
                                                                vi != m.vertices_end(); ++vi,++i)
                                                                if(!is_boundary(vi))
                                                                                lap[i] =  laplacian(vi);
                                                i=0;
                                                for(VertexIter vi = m.vertices_begin(); vi != m.vertices_end(); 
                                                                ++vi,++i)
                                                                if(!is_boundary(vi))
                                                                                vi->pos += (iter%2 == 0) ? 
                                                                                                0.5  * lap[i] : 
                                                                                                -0.52 * lap[i];
                                }
                }

                void face_neighbourhood(FaceIter f, vector<Vec3f>& nbrs)
                {       
                                nbrs.push_back(normal(f));
                                for(FaceCirculator fc(f); !fc.end(); ++fc)
                                {
                                                for(VertexCirculator vc(fc.get_vertex()); !vc.end(); ++vc)
                                                {
                                                                FaceIter fn = vc.get_face();
                                                                if(fn != NULL_FACE_ITER && fn->touched != f->touched)
                                                                {
                                                                                Vec3f norm = normal(fn);
                                                                                if(!isnan(sqr_length(norm)))
                                                                                                nbrs.push_back(norm);
                                                                                else
                                                                                                cout << "bad normal detected" << endl;
                                                                                fn->touched = f->touched;
                                                                }
                                                }
                                }
                }



                Vec3f filtered_normal(FaceIter 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(FaceIter f = m.faces_begin(); f != m.faces_end(); ++f)
                                                                f->touched = -1;

                                                vector<Vec3f> filtered_norms(m.no_faces());
                                                int i=0;
                                                for(FaceIter f = m.faces_begin(); f != m.faces_end(); ++f,++i)
                                                {
                                                                f->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(VertexIter vi=m.vertices_begin(); vi!=m.vertices_end(); ++vi)
                                                                {
                                                                                Vec3f move(0);
                                                                                for(VertexCirculator vc(vi); !vc.end(); ++vc)
                                                                                {
                                                                                                HalfEdgeIter h = vc.get_halfedge();
                                                                                                FaceIter f1 = h->face;
                                                                                                FaceIter f2 = h->opp->face;
                                                                                                Vec3f dir = h->vert->pos - vi->pos;

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


}