Subversion Repositories gelsvn

Rev

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

Rev 595 Rev 600
Line 2... Line 2...
2 
 * This file is part of GEL, http://www.imm.dtu.dk/GEL
 2 
 * This file is part of GEL, http://www.imm.dtu.dk/GEL
3 
 * Copyright (C) the authors and DTU Informatics
 3 
 * Copyright (C) the authors and DTU Informatics
4 
 * For license and list of authors, see ../../doc/intro.pdf
 4 
 * For license and list of authors, see ../../doc/intro.pdf
5 
 * ----------------------------------------------------------------------- */
 5 
 * ----------------------------------------------------------------------- */
6 
 
 6 
 
7 
#include "caps_and_needles.h"
 7 
#include "cleanup.h"
8 
 
 8 
 
9 
#include <CGLA/Vec3f.h>
 9 
#include <CGLA/Vec3f.h>
10 
#include <CGLA/Vec3d.h>
 10 
#include <CGLA/Vec3d.h>
11 
#include <Geometry/QEM.h>
 11 
#include <Geometry/QEM.h>
- 
 
 12 
#include <Geometry/KDTree.h>
12 
 
 13 
 
13 
#include "Manifold.h"
 14 
#include "Manifold.h"
14 
 
 15 
 
15 
namespace HMesh
 16 
namespace HMesh
16 
{
 17 
{
17 
    using namespace std;
 18 
    using namespace std;
18 
    using namespace CGLA;
 19 
    using namespace CGLA;
19 
    using namespace Geometry;
 20 
    using namespace Geometry;
20 
 
 21
21 
    namespace
 22 
    namespace
22 
    {
 23 
    {
23 
        /// get the minimum angle between 3 vertices
 24 
        /// get the minimum angle between 3 vertices
24 
        float min_angle(const Vec3d& v0, const Vec3d& v1, const Vec3d& v2)
 25 
        float min_angle(const Vec3d& v0, const Vec3d& v1, const Vec3d& v2)
25 
        {
 26 
        {
26 
            Vec3d a = normalize(v1-v0);
 27 
            Vec3d a = normalize(v1-v0);
27 
            Vec3d b = normalize(v2-v1);
 28 
            Vec3d b = normalize(v2-v1);
28 
            Vec3d c = normalize(v0-v2);
 29 
            Vec3d c = normalize(v0-v2);
29 
 
 30
30 
            return min(dot(a, -c), min(dot(b, -a), dot(c, -b)));
 31 
            return min(dot(a, -c), min(dot(b, -a), dot(c, -b)));
31 
        }
 32 
        }
32 
 
 33
33 
        /// need description
 34 
        /// need description
34 
        float edge_error(const Manifold& m, HalfEdgeID h, const Vec3d& pa, const Vec3d& pb)
 35 
        float edge_error(const Manifold& m, HalfEdgeID h, const Vec3d& pa, const Vec3d& pb)
35 
        {
 36 
        {
36 
            QEM q;
 37 
            QEM q;
37 
            Walker j = m.walker(h);
 38 
            Walker j = m.walker(h);
38 
 
 39
39 
            FaceID f = j.face();
 40 
            FaceID f = j.face();
40 
            if(f != InvalidFaceID)
 41 
            if(f != InvalidFaceID)
41 
                q += QEM(Vec3d(0), Vec3d(normal(m, f)));
 42 
                q += QEM(Vec3d(0), Vec3d(normal(m, f)));
42 
 
 43
43 
            f = j.opp().face();
 44 
            f = j.opp().face();
44 
            if(f != InvalidFaceID)
 45 
            if(f != InvalidFaceID)
45 
                q += QEM(Vec3d(0), Vec3d(normal(m, f)));
 46 
                q += QEM(Vec3d(0), Vec3d(normal(m, f)));
46 
 
 47
47 
            return q.error(pb - pa);
 48 
            return q.error(pb - pa);
48 
        }
 49 
        }
49 
 
 50
50 
        /// need description
 51 
        /// need description
51 
        float vertex_error(const Manifold& m, VertexID v, const Vec3d& pb)
 52 
        float vertex_error(const Manifold& m, VertexID v, const Vec3d& pb)
52 
        {
 53 
        {
53 
            QEM q;
 54 
            QEM q;
54 
            Vec3d pa(m.pos(v));
 55 
            Vec3d pa(m.pos(v));
55 
 
 56
56 
            for(Walker vj = m.walker(v); !vj.full_circle(); vj = vj.circulate_vertex_cw()){
 57 
            for(Walker vj = m.walker(v); !vj.full_circle(); vj = vj.circulate_vertex_cw()){
57 
                FaceID f = vj.face();
 58 
                FaceID f = vj.face();
58 
                if(f != InvalidFaceID)
 59 
                if(f != InvalidFaceID)
59 
                    q += QEM(Vec3d(0), Vec3d(normal(m, f)));
 60 
                    q += QEM(Vec3d(0), Vec3d(normal(m, f)));
60 
            }
 61 
            }
61 
            return q.error(pb - pa);
 62 
            return q.error(pb - pa);
62 
        }
 63 
        }
63 
    }
 64 
    }
64 
 
 65
65 
    void remove_caps(Manifold& m, float thresh)
 66 
    void remove_caps(Manifold& m, float thresh)
66 
    {
 67 
    {
67 
        for(FaceIDIterator f = m.faces_begin(); f != m.faces_end(); ++f){
 68 
        for(FaceIDIterator f = m.faces_begin(); f != m.faces_end(); ++f){
68 
            Vec3d p[3];
 69 
            Vec3d p[3];
69 
            HalfEdgeID he[3];
 70 
            HalfEdgeID he[3];
70 
            VertexID vh[3];
 71 
            VertexID vh[3];
71 
 
 72
72 
            // store ids of vertices and halfedges and vertex positions of face
 73 
            // store ids of vertices and halfedges and vertex positions of face
73 
            size_t n = 0;
 74 
            size_t n = 0;
74 
            for(Walker fj = m.walker(*f); !fj.full_circle(); fj = fj.circulate_face_cw(), ++n){
 75 
            for(Walker fj = m.walker(*f); !fj.full_circle(); fj = fj.circulate_face_cw(), ++n){
75 
                vh[n] = fj.vertex();
 76 
                vh[n] = fj.vertex();
76 
                p[n] = Vec3d(m.pos(vh[n]));
 77 
                p[n] = Vec3d(m.pos(vh[n]));
77 
 
 - 
 
- 
 
 78
78 
                // original face circulator implementation returned next halfedge, jumper doesn't. Can this be optimized?
79 
                // original face circulator implementation returned next halfedge, jumper doesn't. Can this be optimized?
79 
                he[n] = fj.halfedge();
 80 
                he[n] = fj.halfedge();
80 
            }
 81 
            }
81 
            assert(n == 3);
 82 
            assert(n == 3);
82 
 
 83
83 
            // calculate the edge lengths of face
 84 
            // calculate the edge lengths of face
84 
            bool is_collapsed = false;
 85 
            bool is_collapsed = false;
85 
            Vec3d edges[3];
 86 
            Vec3d edges[3];
86 
            for(size_t i = 0; i < 3; ++i){
 87 
            for(size_t i = 0; i < 3; ++i){
87 
                edges[i] = p[(i+1)%3] - p[i];
 88 
                edges[i] = p[(i+1)%3] - p[i];
88 
                float l = length(edges[i]);
 89 
                float l = length(edges[i]);
89 
                if(l < 1e-20)
 90 
                if(l < 1e-20)
90 
                    is_collapsed = true;
91 
                    is_collapsed = true;
91 
                else
 92 
                else
92 
                    edges[i] /= l;
 93 
                    edges[i] /= l;
93 
 
 - 
 
- 
 
 94
94 
            }
 95 
            }
95 
            // an edge length was close to 1e-20, thus collapsed
 96 
            // an edge length was close to 1e-20, thus collapsed
96 
            if(is_collapsed)
 97 
            if(is_collapsed)
97 
                continue;
 98 
                continue;
98 
 
 99
99 
            for(size_t i = 0; i < 3; ++i){
 100 
            for(size_t i = 0; i < 3; ++i){
100 
                float ang = acos(max(-1.0, min(1.0, dot(-edges[(i+2)%3], edges[i]))));
 101 
                float ang = acos(max(-1.0, min(1.0, dot(-edges[(i+2)%3], edges[i]))));
101 
 
 - 
 
- 
 
 102
102 
                // flip long edge of current face if angle exceeds cap threshold and result is better than current cap
 103 
                // flip long edge of current face if angle exceeds cap threshold and result is better than current cap
103 
                if(ang > thresh){
 104 
                if(ang > thresh){
104 
                    size_t iplus1 = (i+1)%3;
 105 
                    size_t iplus1 = (i+1)%3;
105 
                    Vec3d edge_dir = edges[iplus1];
 106 
                    Vec3d edge_dir = edges[iplus1];
106 
 
 - 
 
- 
 
 107
107 
                    Walker j = m.walker(he[iplus1]);
 108 
                    Walker j = m.walker(he[iplus1]);
108 
                    Vec3d v0(m.pos(j.vertex()));
 109 
                    Vec3d v0(m.pos(j.vertex()));
109 
                    Vec3d v1(m.pos(j.next().vertex()));
 110 
                    Vec3d v1(m.pos(j.next().vertex()));
110 
                    Vec3d v2(m.pos(j.opp().vertex()));
 111 
                    Vec3d v2(m.pos(j.opp().vertex()));
111 
                    Vec3d v3(m.pos(j.opp().next().vertex()));
 112 
                    Vec3d v3(m.pos(j.opp().next().vertex()));
112 
 
 - 
 
- 
 
 113
113 
                    float m1 = min(min_angle(v0, v1, v2), min_angle(v0, v2, v3));
 114 
                    float m1 = min(min_angle(v0, v1, v2), min_angle(v0, v2, v3));
114 
                    float m2 = min(min_angle(v0, v1, v3), min_angle(v1, v2, v3));
 115 
                    float m2 = min(min_angle(v0, v1, v3), min_angle(v1, v2, v3));
115 
 
 - 
 
- 
 
 116
116 
                    if(m1 < m2){
 117 
                    if(m1 < m2){
117 
						// If the "cap vertex" projected onto the long edge is better in the
118 
						// If the "cap vertex" projected onto the long edge is better in the
118 
						// sense that there is less geometric error after the flip, then we
119 
						// sense that there is less geometric error after the flip, then we
119 
						// use the projected vertex. In other words, we see if it pays to snap
 120 
						// use the projected vertex. In other words, we see if it pays to snap
120 
						// to the edge.
 121 
						// to the edge.
121 
						Vec3d pprj = edge_dir * dot(edge_dir, p[i]-p[iplus1])+p[iplus1];
 122 
						Vec3d pprj = edge_dir * dot(edge_dir, p[i]-p[iplus1])+p[iplus1];
122 
                        if(edge_error(m, he[iplus1], pprj, Vec3d(m.pos(vh[i]))) > vertex_error(m, vh[i], pprj))
 123 
                        if(edge_error(m, he[iplus1], pprj, Vec3d(m.pos(vh[i]))) > vertex_error(m, vh[i], pprj))
123 
                            m.pos(vh[i]) = pprj;
 124 
                            m.pos(vh[i]) = pprj;
Line 127... Line 128...
127 
                        break;
 128 
                        break;
128 
                    }
 129 
                    }
129 
                }
 130 
                }
130 
            }
 131 
            }
131 
        }
 132 
        }
132 
 
 133
133 
    }
 134 
    }
134 
 
 135
135 
    void remove_needles(Manifold& m, float thresh)
 136 
    void remove_needles(Manifold& m, float thresh)
136 
    {
 137 
    {
137 
        bool did_work = false;
 138 
        bool did_work = false;
138 
 
 139
139 
        // remove needles until no more can be removed
 140 
        // remove needles until no more can be removed
140 
        do{
 141 
        do{
141 
            did_work = false;
 142 
            did_work = false;
142 
            for(VertexIDIterator v = m.vertices_begin(); v != m.vertices_end(); ++v){
 143 
            for(VertexIDIterator v = m.vertices_begin(); v != m.vertices_end(); ++v){
143 
                // don't attempt to remove needles if vertex is boundary
 144 
                // don't attempt to remove needles if vertex is boundary
144 
                if(boundary(m, *v))
 145 
                if(boundary(m, *v))
145 
                    continue;
 146 
                    continue;
146 
 
 - 
 
- 
 
 147
147 
                for(Walker vj = m.walker(*v); !vj.full_circle(); vj = vj.circulate_vertex_cw()){
 148 
                for(Walker vj = m.walker(*v); !vj.full_circle(); vj = vj.circulate_vertex_cw()){
148 
                    // don't attempt to remove needles if vertex of jumper halfedge is boundary
 149 
                    // don't attempt to remove needles if vertex of jumper halfedge is boundary
149 
   //                 if(boundary(m, vj.vertex()))
 150 
                    //                 if(boundary(m, vj.vertex()))
150 
//                        continue;
 151 
                    //                        continue;
151 
 
 - 
 
- 
 
 152
152 
                    HalfEdgeID h = vj.opp().halfedge();
 153 
                    HalfEdgeID h = vj.opp().halfedge();
153 
//                    VertexID n = vj.vertex();
 154 
                    //                    VertexID n = vj.vertex();
154 
                    float dist = length(m, h);
 155 
                    float dist = length(m, h);
155 
 
 - 
 
- 
 
 156
156 
                    // collapse edge if allowed and needle is present
 157 
                    // collapse edge if allowed and needle is present
157 
                    if(dist < thresh && precond_collapse_edge(m, h)){
 158 
                    if(dist < thresh && precond_collapse_edge(m, h)){
158 
//                        if(vertex_error(m, *v, Vec3d(m.pos(n))) < vertex_error(m, n, Vec3d(m.pos(*v))))
 159 
                        //                        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.pos(*v) = m.pos(n);
160 
                        m.collapse_edge(h);
 161 
                        m.collapse_edge(h);
161 
                        did_work = true;
 162 
                        did_work = true;
162 
                        break;
 163 
                        break;
163 
                    }
 164 
                    }
164 
                }
 165 
                }
165 
            }
 166 
            }
166 
        }
 167 
        }
167 
        while(did_work);
 168 
        while(did_work);
168 
    }
 169 
    }
- 
 
 170
- 
 
 171 
    int stitch_mesh(Manifold& m)
- 
 
 172 
    {
- 
 
 173 
        KDTree<Vec3d, VertexID> vtree;
- 
 
 174
- 
 
 175 
        for(VertexIDIterator vid = m.vertices_begin(); vid != m.vertices_end(); ++vid)
- 
 
 176 
            if(boundary(m, *vid))
- 
 
 177 
                vtree.insert(m.pos(*vid), *vid);
- 
 
 178 
        vtree.build();
- 
 
 179
- 
 
 180 
        VertexAttributeVector<int> cluster_id(m.allocated_vertices(),-1);
- 
 
 181 
        vector<vector<HalfEdgeID> > clustered_halfedges;
- 
 
 182
- 
 
 183 
        int cluster_ctr=0;
- 
 
 184 
        for(VertexIDIterator vid = m.vertices_begin(); vid != m.vertices_end(); ++vid)
- 
 
 185 
            if(boundary(m, *vid) && cluster_id[*vid] == -1)
- 
 
 186 
            {
- 
 
 187 
                vector<Vec3d> keys;
- 
 
 188 
                vector<VertexID> vals;
- 
 
 189 
                int n = vtree.in_sphere(m.pos(*vid), 1e-5, keys, vals);
- 
 
 190
- 
 
 191 
                vector<HalfEdgeID> boundary_edges;
- 
 
 192 
                for(int i=0;i<n;++i)
- 
 
 193 
                {
- 
 
 194 
                    cluster_id[vals[i]] = cluster_ctr;
- 
 
 195 
                    Walker w = m.walker(vals[i]);
- 
 
 196 
                    boundary_edges.push_back(w.halfedge());
- 
 
 197 
                }
- 
 
 198 
                clustered_halfedges.push_back(boundary_edges);
- 
 
 199 
                ++cluster_ctr;
- 
 
 200 
            }
- 
 
 201
- 
 
 202 
        int unstitched=0;
- 
 
 203 
        for(HalfEdgeIDIterator hid = m.halfedges_begin(); hid != m.halfedges_end(); ++hid)
- 
 
 204 
        {
- 
 
 205 
            HalfEdgeID h0 = *hid;
- 
 
 206 
            Walker w = m.walker(h0);
- 
 
 207 
            if(w.face() == InvalidFaceID)
- 
 
 208 
            {
- 
 
 209 
                VertexID v0 = w.opp().vertex();
- 
 
 210 
                VertexID v1 = w.vertex();
- 
 
 211
- 
 
 212 
                int cid = cluster_id[v1];
- 
 
 213 
                vector<HalfEdgeID>& stitch_candidates = clustered_halfedges[cid];
- 
 
 214 
                int i=0;
- 
 
 215 
                for(;i<stitch_candidates.size(); ++i)
- 
 
 216 
                {
- 
 
 217 
                    HalfEdgeID h1 = stitch_candidates[i];
- 
 
 218 
                    if(m.in_use(h1))
- 
 
 219 
                    {
- 
 
 220 
                        Walker w = m.walker(h1);
- 
 
 221 
                        if(cluster_id[w.vertex()] == cluster_id[v0])
- 
 
 222 
                            if(m.stitch_boundary_edges(h0,h1))
- 
 
 223 
                                break;
- 
 
 224 
                    }
- 
 
 225
- 
 
 226 
                }
- 
 
 227 
                if(i == stitch_candidates.size())
- 
 
 228 
                    ++unstitched;
- 
 
 229 
            }
- 
 
 230 
        }
- 
 
 231 
        return unstitched;
- 
 
 232 
    }
- 
 
 233
- 
 
 234
- 
 
 235 
    void stitch_more(Manifold& mani)
- 
 
 236 
    {
- 
 
 237 
        int unstitched, iter;
- 
 
 238 
        for(iter=0;iter<2;++iter)
- 
 
 239 
        {
- 
 
 240 
            unstitched = stitch_mesh(mani);
- 
 
 241 
            if(unstitched == 0)
- 
 
 242 
                break;
- 
 
 243
- 
 
 244 
            vector<HalfEdgeID> hvec;
- 
 
 245 
            for(HalfEdgeIDIterator h = mani.halfedges_begin(); h != mani.halfedges_end();++h)
- 
 
 246 
                if(mani.walker(*h).face() == InvalidFaceID)
- 
 
 247 
                    hvec.push_back(*h);
- 
 
 248 
            for(int i=0;i<hvec.size(); ++i)
- 
 
 249 
                mani.split_edge(hvec[i]);
- 
 
 250
- 
 
 251 
        }
- 
 
 252 
    }
- 
 
 253
- 
 
 254
- 
 
 255 
    void close_holes(Manifold& m)
- 
 
 256 
    {
- 
 
 257 
        for(HalfEdgeIDIterator h =  m.halfedges_begin(); h != m.halfedges_end(); ++h){
- 
 
 258 
            m.close_hole(*h);
- 
 
 259 
        }
- 
 
 260 
    }
- 
 
 261
169 
}
 262 
}
170 
 
 263