Subversion Repositories gelsvn

Rev

Rev 215 | Rev 600 | Go to most recent revision | Show entire file | Ignore whitespace | Details | Blame | Last modification | View Log | RSS feed

Rev 215 Rev 595
Line -... Line 1...
- 
 
 1 
/* ----------------------------------------------------------------------- *
1 
#include <queue>
 2 
 * This file is part of GEL, http://www.imm.dtu.dk/GEL
2 
#include "HMesh/FaceCirculator.h"
 3 
 * Copyright (C) the authors and DTU Informatics
3 
#include "HMesh/VertexCirculator.h"
 4 
 * For license and list of authors, see ../../doc/intro.pdf
- 
 
 5 
 * ----------------------------------------------------------------------- */
- 
 
 6 
 
4 
#include "caps_and_needles.h"
 7 
#include "caps_and_needles.h"
5 
#include "Geometry/QEM.h"
 - 
 
6 
 
 8 
 
7 
using namespace std;
 9 
#include <CGLA/Vec3f.h>
8 
using namespace CGLA;
 10 
#include <CGLA/Vec3d.h>
9 
using namespace GEO;
 - 
 
10 
using namespace HMesh;
 11 
#include <Geometry/QEM.h>
11 
 
 12 
 
- 
 
 13 
#include "Manifold.h"
12 
 
 14 
 
13 
namespace HMesh
 15 
namespace HMesh
14 
{
 16 
{
15 
	namespace
 17 
    using namespace std;
16 
	{
 18 
    using namespace CGLA;
17 
		float min_angle(const Vec3d& v0,const Vec3d& v1,const Vec3d& v2)
 19 
    using namespace Geometry;
18 
		{
 20 
 
19 
			Vec3d a = normalize(v1-v0);
 21 
    namespace
20 
			Vec3d b = normalize(v2-v1);
 22 
    {
21 
			Vec3d c = normalize(v0-v2);
 23 
        /// get the minimum angle between 3 vertices
22
24 
        float min_angle(const Vec3d& v0, const Vec3d& v1, const Vec3d& v2)
23 
			return min(dot(a,-c), min(dot(b,-a), dot(c,-b)));
 25 
        {
24 
		}
 26 
            Vec3d a = normalize(v1-v0);
25 
 
 27 
            Vec3d b = normalize(v2-v1);
26 
		float edge_error(HalfEdgeIter h,
28 
            Vec3d c = normalize(v0-v2);
27 
										 const Vec3d& pa,
 29 
 
28 
										 const Vec3d& pb)
 30 
            return min(dot(a, -c), min(dot(b, -a), dot(c, -b)));
29 
		{
 31 
        }
30 
			QEM q;
 32 
 
31 
			if(h->face != NULL_FACE_ITER)
 33 
        /// need description
32 
				q += QEM(Vec3d(0), Vec3d(normal(h->face)));
 34 
        float edge_error(const Manifold& m, HalfEdgeID h, const Vec3d& pa, const Vec3d& pb)
33 
			if(h->opp->face != NULL_FACE_ITER)
 35 
        {
34 
				q += QEM(Vec3d(0), Vec3d(normal(h->opp->face)));
 36 
            QEM q;
35 
			return q.error(pb-pa);
 37 
            Walker j = m.walker(h);
36 
		}
 38 
 
37 
 
 39 
            FaceID f = j.face();
38 
 
 40 
            if(f != InvalidFaceID)
39 
		float vertex_error(VertexIter va, const Vec3d& pb)
 41 
                q += QEM(Vec3d(0), Vec3d(normal(m, f)));
40 
		{
 42 
 
41 
			QEM q;
 43 
            f = j.opp().face();
42 
			Vec3d pa(va->pos);
 44 
            if(f != InvalidFaceID)
43 
			for(VertexCirculator vc(va); !vc.end(); ++vc)
 45 
                q += QEM(Vec3d(0), Vec3d(normal(m, f)));
44 
				{
 46 
 
45 
					FaceIter f = vc.get_face();
 47 
            return q.error(pb - pa);
46 
					if(f != NULL_FACE_ITER)
 48 
        }
47 
						{
 49 
 
48 
							Vec3d n(normal(f));
 50 
        /// need description
49 
							q += QEM(Vec3d(0), Vec3d(n));
 51 
        float vertex_error(const Manifold& m, VertexID v, const Vec3d& pb)
50 
						}
 52 
        {
51 
				}
 53 
            QEM q;
52 
			return q.error(pb-pa);
 54 
            Vec3d pa(m.pos(v));
53 
		}
 55 
 
54 
 
 56 
            for(Walker vj = m.walker(v); !vj.full_circle(); vj = vj.circulate_vertex_cw()){
55 
	}
 57 
                FaceID f = vj.face();
56 
 
 58 
                if(f != InvalidFaceID)
57 
	void remove_caps_from_trimesh(Manifold& mani, float ang_thresh)
 59 
                    q += QEM(Vec3d(0), Vec3d(normal(m, f)));
58 
	{
 60 
            }
59 
		for(FaceIter fi=mani.faces_begin(); fi != mani.faces_end(); ++fi)
 61 
            return q.error(pb - pa);
60 
			{
 62 
        }
61 
				Vec3d p[3];
 63 
    }
62 
				HalfEdgeIter he[3];
 64 
 
63 
				VertexIter vi[3];
 65 
    void remove_caps(Manifold& m, float thresh)
64 
 
 66 
    {
65 
				int n=0;
 67 
        for(FaceIDIterator f = m.faces_begin(); f != m.faces_end(); ++f){
66 
				for(FaceCirculator fc(fi); !fc.end(); ++fc,++n)
 68 
            Vec3d p[3];
67 
					{
 69 
            HalfEdgeID he[3];
68 
						vi[n] = fc.get_vertex();
 70 
            VertexID vh[3];
69 
						p[n]  = Vec3d(fc.get_vertex()->pos);
 71 
 
70 
						he[n]= fc.get_halfedge();
 72 
            // store ids of vertices and halfedges and vertex positions of face
71 
					}
 73 
            size_t n = 0;
72 
				assert(n=3);
 74 
            for(Walker fj = m.walker(*f); !fj.full_circle(); fj = fj.circulate_face_cw(), ++n){
73
75 
                vh[n] = fj.vertex();
74 
				bool is_collapsed = false;
 76 
                p[n] = Vec3d(m.pos(vh[n]));
75 
				Vec3d edges[3];
 77 
 
76 
				for(int i=0;i<3;++i)
 78 
                // original face circulator implementation returned next halfedge, jumper doesn't. Can this be optimized?
77 
					{
 79 
                he[n] = fj.halfedge();
78 
						edges[i] = p[(i+1)%3]-p[i];
 80 
            }
79 
						float l = length(edges[i]);
 81 
            assert(n == 3);
80 
						if(l<1e-20)
 82 
 
81 
							is_collapsed = true;
 83 
            // calculate the edge lengths of face
82 
						else
 84 
            bool is_collapsed = false;
83 
							edges[i] /= l;
 85 
            Vec3d edges[3];
84 
					}
 86 
            for(size_t i = 0; i < 3; ++i){
85 
				if(is_collapsed) continue;
 87 
                edges[i] = p[(i+1)%3] - p[i];
86
88 
                float l = length(edges[i]);
87 
				for(int i=0;i<3;++i)
 89 
                if(l < 1e-20)
88 
					{
 90 
                    is_collapsed = true;
89 
							float d = max(-1.0,min(1.0,dot(-edges[(i+2)%3],edges[i])));
 91 
                else
90 
						float ang = acos(d);
 92 
                    edges[i] /= l;
91 
						if(ang > ang_thresh)
 93 
 
92 
							{
 94 
            }
93 
								int iplus1 = (i+1)%3;
 95 
            // an edge length was close to 1e-20, thus collapsed
94 
								Vec3d edge_dir = edges[iplus1];
 96 
            if(is_collapsed)
95 
								Vec3d pprj =
97 
                continue;
96 
									edge_dir * dot(edge_dir,p[i]-p[iplus1])+p[iplus1];
 98 
 
97 
 
 99 
            for(size_t i = 0; i < 3; ++i){
98 
								HalfEdgeIter h = he[iplus1];
 100 
                float ang = acos(max(-1.0, min(1.0, dot(-edges[(i+2)%3], edges[i]))));
99 
								Vec3d v0(h->vert->pos);
 101 
 
100 
								Vec3d v1(h->next->vert->pos);
 102 
                // flip long edge of current face if angle exceeds cap threshold and result is better than current cap
101 
								Vec3d v2(h->opp->vert->pos);
 103 
                if(ang > thresh){
102 
								Vec3d v3(h->opp->next->vert->pos);
 104 
                    size_t iplus1 = (i+1)%3;
103
105 
                    Vec3d edge_dir = edges[iplus1];
104 
								float m1 = min(min_angle(v0,v1,v2), min_angle(v0,v2,v3));
 106 
 
105 
								float m2 = min(min_angle(v0,pprj,v3), min_angle(pprj,v2,v3));
 107 
                    Walker j = m.walker(he[iplus1]);
106 
								if(m1 < m2)
 108 
                    Vec3d v0(m.pos(j.vertex()));
107 
									{
 109 
                    Vec3d v1(m.pos(j.next().vertex()));
108 
											if(edge_error(he[iplus1], pprj, Vec3d(vi[i]->pos)) >
 110 
                    Vec3d v2(m.pos(j.opp().vertex()));
109 
											 vertex_error(vi[i], pprj))
 111 
                    Vec3d v3(m.pos(j.opp().next().vertex()));
110 
												vi[i]->pos = Vec3f(pprj);
 112 
 
111 
 
 113 
                    float m1 = min(min_angle(v0, v1, v2), min_angle(v0, v2, v3));
112 
										mani.flip(he[iplus1]);
 114 
                    float m2 = min(min_angle(v0, v1, v3), min_angle(v1, v2, v3));
113 
										break;
 115 
 
114 
									}
 116 
                    if(m1 < m2){
115 
							}
 117 
						// If the "cap vertex" projected onto the long edge is better in the
116 
					}
 118 
						// sense that there is less geometric error after the flip, then we
117
119 
						// use the projected vertex. In other words, we see if it pays to snap
118 
			}
 120 
						// to the edge.
119 
	}
 121 
						Vec3d pprj = edge_dir * dot(edge_dir, p[i]-p[iplus1])+p[iplus1];
120 
 
 122 
                        if(edge_error(m, he[iplus1], pprj, Vec3d(m.pos(vh[i]))) > vertex_error(m, vh[i], pprj))
121 
 
 123 
                            m.pos(vh[i]) = pprj;
122 
 
 124 
                        // flip if legal
123 
	void remove_needles_from_trimesh(Manifold& mani,
125 
                        if(precond_flip_edge(m, he[iplus1]))
124 
																	 float thresh)
 126 
                            m.flip_edge(he[iplus1]);
125 
	{
 127 
                        break;
126 
		bool did_work;
 128 
                    }
127 
		do
 129 
                }
128 
			{
 130 
            }
129 
				did_work = false;
 131 
        }
130 
				for(VertexIter vi=mani.vertices_begin();
132 
 
131 
						vi != mani.vertices_end(); ++vi)
 133 
    }
132 
					{
 134 
 
133 
							if(is_boundary(vi)) continue;
 135 
    void remove_needles(Manifold& m, float thresh)
134 
						for(VertexCirculator vc(vi);!vc.end();++vc)
 136 
    {
135 
							{
 137 
        bool did_work = false;
136 
								if(is_boundary(vc.get_halfedge()->vert)) continue;
 138 
 
137 
								HalfEdgeIter he = vc.get_halfedge()->opp;
 139 
        // remove needles until no more can be removed
138 
 
 140 
        do{
139 
								VertexIter n = vc.get_vertex();
 141 
            did_work = false;
140 
								float dist = length(he);
 142 
            for(VertexIDIterator v = m.vertices_begin(); v != m.vertices_end(); ++v){
141 
								if((dist<thresh) && mani.collapse_precond(he))
 143 
                // don't attempt to remove needles if vertex is boundary
142 
									{
 144 
                if(boundary(m, *v))
143 
											if(vertex_error(vi,Vec3d(n->pos)) <
145 
                    continue;
144 
												 vertex_error(n,Vec3d(vi->pos)))
 146 
 
145 
													vi->pos = n->pos;
 147 
                for(Walker vj = m.walker(*v); !vj.full_circle(); vj = vj.circulate_vertex_cw()){
146 
 
 148 
                    // don't attempt to remove needles if vertex of jumper halfedge is boundary
147 
										mani.collapse_halfedge(he);
 149 
   //                 if(boundary(m, vj.vertex()))
148 
										did_work = true;
 150 
//                        continue;
149 
										break;
 151 
 
150 
									}
 152 
                    HalfEdgeID h = vj.opp().halfedge();
151 
							}
 153 
//                    VertexID n = vj.vertex();
152 
					}
 154 
                    float dist = length(m, h);
153 
			}
 155 
 
154 
		while(did_work);
 156 
                    // collapse edge if allowed and needle is present
155 
	}
 157 
                    if(dist < thresh && precond_collapse_edge(m, h)){
156 
}
 158 
//                        if(vertex_error(m, *v, Vec3d(m.pos(n))) < vertex_error(m, n, Vec3d(m.pos(*v))))
- 
 
 159 
//                            m.pos(*v) = m.pos(n);
- 
 
 160 
                        m.collapse_edge(h);
- 
 
 161 
                        did_work = true;
- 
 
 162 
                        break;
- 
 
 163 
                    }
- 
 
 164 
                }
- 
 
 165 
            }
- 
 
 166 
        }
- 
 
 167 
        while(did_work);
- 
 
 168 
    }
- 
 
 169 
}
157 
 
 170