Subversion Repositories gelsvn

#include <queue>
#include "HMesh/FaceCirculator.h"
#include "HMesh/VertexCirculator.h"
#include "caps_and_needles.h"
#include "Geometry/QEM.h"

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


namespace HMesh
{
        namespace
        {
                float min_angle(const Vec3d& v0,const Vec3d& v1,const Vec3d& v2)
                {
                        Vec3d a = normalize(v1-v0);
                        Vec3d b = normalize(v2-v1);
                        Vec3d c = normalize(v0-v2);
                        
                        return min(dot(a,-c), min(dot(b,-a), dot(c,-b)));
                }

                float edge_error(HalfEdgeIter h, 
                                                                                 const Vec3d& pa,
                                                                                 const Vec3d& pb)
                {
                        QEM q;
                        if(h->face != NULL_FACE_ITER)
                                q += QEM(Vec3d(0), Vec3d(normal(h->face)));
                        if(h->opp->face != NULL_FACE_ITER)
                                q += QEM(Vec3d(0), Vec3d(normal(h->opp->face)));
                        return q.error(pb-pa);
                }


                float vertex_error(VertexIter va, const Vec3d& pb)
                {
                        QEM q;
                        Vec3d pa(va->pos);
                        for(VertexCirculator vc(va); !vc.end(); ++vc)
                                {
                                        FaceIter f = vc.get_face();
                                        if(f != NULL_FACE_ITER)
                                                {
                                                        Vec3d n(normal(f));
                                                        q += QEM(Vec3d(0), Vec3d(n));
                                                }
                                }
                        return q.error(pb-pa);
                }

        }

        void remove_caps_from_trimesh(Manifold& mani, float ang_thresh)
        {
                for(FaceIter fi=mani.faces_begin(); fi != mani.faces_end(); ++fi)
                        {
                                Vec3d p[3];
                                HalfEdgeIter he[3];
                                VertexIter vi[3];

                                int n=0;
                                for(FaceCirculator fc(fi); !fc.end(); ++fc,++n)
                                        {
                                                vi[n] = fc.get_vertex();
                                                p[n]  = Vec3d(fc.get_vertex()->pos);
                                                he[n]= fc.get_halfedge();
                                        }
                                assert(n=3);
                                
                                bool is_collapsed = false;
                                Vec3d edges[3];
                                for(int i=0;i<3;++i)
                                        {
                                                edges[i] = p[(i+1)%3]-p[i];
                                                float l = length(edges[i]);
                                                if(l<1e-20)
                                                        is_collapsed = true;
                                                else
                                                        edges[i] /= l;
                                        }
                                if(is_collapsed) continue;
                                
                                for(int i=0;i<3;++i)
                                        {
                                                        float d = max(-1.0,min(1.0,dot(-edges[(i+2)%3],edges[i])));
                                                float ang = acos(d);
                                                if(ang > ang_thresh)
                                                        {
                                                                int iplus1 = (i+1)%3;
                                                                Vec3d edge_dir = edges[iplus1];
                                                                Vec3d pprj =  
                                                                        edge_dir * dot(edge_dir,p[i]-p[iplus1])+p[iplus1];

                                                                HalfEdgeIter h = he[iplus1];
                                                                Vec3d v0(h->vert->pos);
                                                                Vec3d v1(h->next->vert->pos);
                                                                Vec3d v2(h->opp->vert->pos);
                                                                Vec3d v3(h->opp->next->vert->pos);
                                                                
                                                                float m1 = min(min_angle(v0,v1,v2), min_angle(v0,v2,v3));
                                                                float m2 = min(min_angle(v0,pprj,v3), min_angle(pprj,v2,v3));
                                                                if(m1 < m2)
                                                                        {
                                                                                        if(edge_error(he[iplus1], pprj, Vec3d(vi[i]->pos)) >
                                                                                         vertex_error(vi[i], pprj))
                                                                                                vi[i]->pos = Vec3f(pprj);

                                                                                mani.flip(he[iplus1]);
                                                                                break;
                                                                        }
                                                        }
                                        }
                                
                        }
        }



        void remove_needles_from_trimesh(Manifold& mani, 
                                                                                                                                         float thresh)
        {
                bool did_work;
                do
                        {
                                did_work = false;
                                for(VertexIter vi=mani.vertices_begin(); 
                                                vi != mani.vertices_end(); ++vi)
                                        {
                                                        if(is_boundary(vi)) continue;
                                                for(VertexCirculator vc(vi);!vc.end();++vc)
                                                        {
                                                                if(is_boundary(vc.get_halfedge()->vert)) continue;
                                                                HalfEdgeIter he = vc.get_halfedge()->opp;

                                                                VertexIter n = vc.get_vertex();
                                                                float dist = length(he);
                                                                if((dist<thresh) && mani.collapse_precond(he))
                                                                        {
                                                                                        if(vertex_error(vi,Vec3d(n->pos)) < 
                                                                                                 vertex_error(n,Vec3d(vi->pos)))
                                                                                                        vi->pos = n->pos;

                                                                                mani.collapse_halfedge(he);
                                                                                did_work = true;
                                                                                break;
                                                                        }
                                                        }
                                        }
                        }
                while(did_work);
        }
}