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 "caps_and_needles.h"

#include <CGLA/Vec3f.h>
#include <CGLA/Vec3d.h>
#include <Geometry/QEM.h>

#include "Manifold.h"

namespace HMesh
{
    using namespace std;
    using namespace CGLA;
    using namespace Geometry;

    namespace
    {
        /// get the minimum angle between 3 vertices
        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)));
        }

        /// need description
        float edge_error(const Manifold& m, HalfEdgeID h, const Vec3d& pa, const Vec3d& pb)
        {
            QEM q;
            Walker j = m.walker(h);

            FaceID f = j.face();
            if(f != InvalidFaceID)
                q += QEM(Vec3d(0), Vec3d(normal(m, f)));

            f = j.opp().face();
            if(f != InvalidFaceID)
                q += QEM(Vec3d(0), Vec3d(normal(m, f)));

            return q.error(pb - pa);
        }

        /// need description
        float vertex_error(const Manifold& m, VertexID v, const Vec3d& pb)
        {
            QEM q;
            Vec3d pa(m.pos(v));

            for(Walker vj = m.walker(v); !vj.full_circle(); vj = vj.circulate_vertex_cw()){
                FaceID f = vj.face();
                if(f != InvalidFaceID)
                    q += QEM(Vec3d(0), Vec3d(normal(m, f)));
            }
            return q.error(pb - pa);
        }
    }

    void remove_caps(Manifold& m, float thresh)
    {
        for(FaceIDIterator f = m.faces_begin(); f != m.faces_end(); ++f){
            Vec3d p[3];
            HalfEdgeID he[3];
            VertexID vh[3];

            // store ids of vertices and halfedges and vertex positions of face
            size_t n = 0;
            for(Walker fj = m.walker(*f); !fj.full_circle(); fj = fj.circulate_face_cw(), ++n){
                vh[n] = fj.vertex();
                p[n] = Vec3d(m.pos(vh[n]));

                // original face circulator implementation returned next halfedge, jumper doesn't. Can this be optimized? 
                he[n] = fj.halfedge();
            }
            assert(n == 3);

            // calculate the edge lengths of face
            bool is_collapsed = false;
            Vec3d edges[3];
            for(size_t 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;

            }
            // an edge length was close to 1e-20, thus collapsed
            if(is_collapsed)
                continue;

            for(size_t i = 0; i < 3; ++i){
                float ang = acos(max(-1.0, min(1.0, dot(-edges[(i+2)%3], edges[i]))));

                // flip long edge of current face if angle exceeds cap threshold and result is better than current cap
                if(ang > thresh){
                    size_t iplus1 = (i+1)%3;
                    Vec3d edge_dir = edges[iplus1];

                    Walker j = m.walker(he[iplus1]);
                    Vec3d v0(m.pos(j.vertex()));
                    Vec3d v1(m.pos(j.next().vertex()));
                    Vec3d v2(m.pos(j.opp().vertex()));
                    Vec3d v3(m.pos(j.opp().next().vertex()));

                    float m1 = min(min_angle(v0, v1, v2), min_angle(v0, v2, v3));
                    float m2 = min(min_angle(v0, v1, v3), min_angle(v1, v2, v3));

                    if(m1 < m2){
                                                // If the "cap vertex" projected onto the long edge is better in the 
                                                // sense that there is less geometric error after the flip, then we 
                                                // use the projected vertex. In other words, we see if it pays to snap
                                                // to the edge.
                                                Vec3d pprj = edge_dir * dot(edge_dir, p[i]-p[iplus1])+p[iplus1];
                        if(edge_error(m, he[iplus1], pprj, Vec3d(m.pos(vh[i]))) > vertex_error(m, vh[i], pprj))
                            m.pos(vh[i]) = pprj;
                        // flip if legal
                        if(precond_flip_edge(m, he[iplus1]))
                            m.flip_edge(he[iplus1]);
                        break;
                    }
                }
            }
        }

    }

    void remove_needles(Manifold& m, float thresh)
    {
        bool did_work = false;

        // remove needles until no more can be removed
        do{
            did_work = false;
            for(VertexIDIterator v = m.vertices_begin(); v != m.vertices_end(); ++v){
                // don't attempt to remove needles if vertex is boundary
                if(boundary(m, *v))
                    continue;

                for(Walker vj = m.walker(*v); !vj.full_circle(); vj = vj.circulate_vertex_cw()){
                    // don't attempt to remove needles if vertex of jumper halfedge is boundary
   //                 if(boundary(m, vj.vertex()))
//                        continue;

                    HalfEdgeID h = vj.opp().halfedge();
//                    VertexID n = vj.vertex();
                    float dist = length(m, h);

                    // collapse edge if allowed and needle is present
                    if(dist < thresh && precond_collapse_edge(m, h)){
//                        if(vertex_error(m, *v, Vec3d(m.pos(n))) < vertex_error(m, n, Vec3d(m.pos(*v))))
//                            m.pos(*v) = m.pos(n);
                        m.collapse_edge(h);
                        did_work = true;
                        break;
                    }
                }
            }
        }
        while(did_work);
    }
}