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);

    }

}