WebSVN – gelsvn – Diff – /branches/GEL_github/src/HMesh/cleanup.cpp

 /* ----------------------------------------------------------------------- *
  * 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 "cleanup.h"
 #include "../CGLA/Vec3f.h"
 #include "../CGLA/Vec3d.h"
 #include "../Geometry/QEM.h"
 #include "../Geometry/KDTree.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);
     }
-    int stitch_mesh(Manifold& m)
+    int remove_duplicates(Manifold& m, double rad)
     {
         KDTree<Vec3d, VertexID> vtree;
-        for(VertexIDIterator vid = m.vertices_begin(); vid != m.vertices_end(); ++vid)
+        for(VertexID v: m.vertices())
-            if(boundary(m, *vid))
+                vtree.insert(m.pos(v), v);
-                vtree.insert(m.pos(*vid), *vid);
+        vtree.build();
-        vtree.build();
+        VertexAttributeVector<int> vertex_dup(m.allocated_vertices(),0);
-        VertexAttributeVector<int> cluster_id(m.allocated_vertices(),-1);
+        vector<vector<HalfEdgeID> > clustered_halfedges;
-        vector<vector<HalfEdgeID> > clustered_halfedges;
+        int cnt = 0;
-        int cluster_ctr=0;
+        for(VertexID v: m.vertices())
-        for(VertexIDIterator vid = m.vertices_begin(); vid != m.vertices_end(); ++vid)
+            {
-            if(boundary(m, *vid) && cluster_id[*vid] == -1)
+                vector<Vec3d> keys;
-            {
+                vector<VertexID> vals;
-                vector<Vec3d> keys;
+                int n = vtree.in_sphere(m.pos(v), rad, keys, vals);
-                vector<VertexID> vals;
-                int n = vtree.in_sphere(m.pos(*vid), 1e-5, keys, vals);
+                if(n>1) {
+                    sort(vals.begin(), vals.end());
-                vector<HalfEdgeID> boundary_edges;
-                for(int i=0;i<n;++i)
+                    auto v2remove = vals.begin();
-                {
+                    for(v2remove++; v2remove != vals.end(); v2remove++) {
-                    cluster_id[vals[i]] = cluster_ctr;
+                        m.remove_vertex(*v2remove);
-                    Walker w = m.walker(vals[i]);
+                        ++cnt;
-                    boundary_edges.push_back(w.halfedge());
+                    }
                 }
-                clustered_halfedges.push_back(boundary_edges);
+            }
-                ++cluster_ctr;
+        return  cnt;
-            }
+    }
-        int unstitched=0;
-        for(HalfEdgeIDIterator hid = m.halfedges_begin(); hid != m.halfedges_end(); ++hid)
+    int stitch_mesh(Manifold& m, double rad)
-        {
+    {
-            HalfEdgeID h0 = *hid;
+        KDTree<Vec3d, VertexID> vtree;
-            Walker w = m.walker(h0);
-            if(w.face() == InvalidFaceID)
+        for(VertexIDIterator vid = m.vertices_begin(); vid != m.vertices_end(); ++vid)
-            {
+            if(boundary(m, *vid))
-                VertexID v0 = w.opp().vertex();
+                vtree.insert(m.pos(*vid), *vid);
-                VertexID v1 = w.vertex();
+        vtree.build();
-                int cid = cluster_id[v1];
+        VertexAttributeVector<int> cluster_id(m.allocated_vertices(),-1);
-                vector<HalfEdgeID>& stitch_candidates = clustered_halfedges[cid];
+        vector<vector<HalfEdgeID> > clustered_halfedges;
-                size_t i=0;
-                for(;i<stitch_candidates.size(); ++i)
+        int cluster_ctr=0;
-                {
+        for(VertexIDIterator vid = m.vertices_begin(); vid != m.vertices_end(); ++vid)
-                    HalfEdgeID h1 = stitch_candidates[i];
+            if(boundary(m, *vid) && cluster_id[*vid] == -1)
-                    if(m.in_use(h1))
+            {
-                    {
+                vector<Vec3d> keys;
-                        Walker w = m.walker(h1);
+                vector<VertexID> vals;
-                        if(cluster_id[w.vertex()] == cluster_id[v0])
+                int n = vtree.in_sphere(m.pos(*vid), rad, keys, vals);
-                            if(m.stitch_boundary_edges(h0,h1))
-                                break;
+                vector<HalfEdgeID> boundary_edges;
-                    }
+                for(int i=0;i<n;++i)
+                {
-                }
+                    cluster_id[vals[i]] = cluster_ctr;
-                if(i == stitch_candidates.size())
+                    Walker w = m.walker(vals[i]);
-                    ++unstitched;
+                    boundary_edges.push_back(w.halfedge());
-            }
+                }
-        }
+                clustered_halfedges.push_back(boundary_edges);
-        return unstitched;
+                ++cluster_ctr;
-    }
+            }
+        int unstitched=0;
-    void stitch_more(Manifold& mani)
+        for(HalfEdgeIDIterator hid = m.halfedges_begin(); hid != m.halfedges_end(); ++hid)
-    {
+        {
-        int unstitched, iter;
+            HalfEdgeID h0 = *hid;
-        for(iter=0;iter<2;++iter)
+            Walker w = m.walker(h0);
-        {
+            if(w.face() == InvalidFaceID)
-            unstitched = stitch_mesh(mani);
+            {
-            if(unstitched == 0)
+                VertexID v0 = w.opp().vertex();
-                break;
+                VertexID v1 = w.vertex();
-            vector<HalfEdgeID> hvec;
+                int cid = cluster_id[v1];
-            for(HalfEdgeIDIterator h = mani.halfedges_begin(); h != mani.halfedges_end();++h)
+                vector<HalfEdgeID>& stitch_candidates = clustered_halfedges[cid];
-                if(mani.walker(*h).face() == InvalidFaceID)
+                size_t i=0;
-                    hvec.push_back(*h);
+                for(;i<stitch_candidates.size(); ++i)
-            for(size_t i=0;i<hvec.size(); ++i)
+                {
-                mani.split_edge(hvec[i]);
+                    HalfEdgeID h1 = stitch_candidates[i];
+                    if(m.in_use(h1))
-        }
+                    {
-    }
+                        Walker w = m.walker(h1);
+                        if(cluster_id[w.vertex()] == cluster_id[v0])
+                            if(m.stitch_boundary_edges(h0,h1))
-    void close_holes(Manifold& m)
+                                break;
-    {
+                    }
-        for(HalfEdgeIDIterator h =  m.halfedges_begin(); h != m.halfedges_end(); ++h){
-            m.close_hole(*h);
+                }
-        }
+                if(i == stitch_candidates.size())
-    }
+                    ++unstitched;
+            }
+        }
+        return unstitched;
+    }
+    void remove_valence_two_vertices(Manifold & m)
+    {
+        for(VertexID v: m.vertices())
+            if(valency(m,v)==2)
+            {
+                Walker w = m.walker(v);
+                if(precond_collapse_edge(m, w.halfedge()))
+                    m.collapse_edge(w.halfedge());
+            }
+    }
+    void stitch_more(Manifold& mani, double rad)
+    {
+        int unstitched, iter;
+        for(iter=0;iter<2;++iter)
+        {
+            unstitched = stitch_mesh(mani, rad);
+            if(unstitched == 0)
+                break;
+            vector<HalfEdgeID> hvec;
+            for(HalfEdgeIDIterator h = mani.halfedges_begin(); h != mani.halfedges_end();++h)
+                if(mani.walker(*h).face() == InvalidFaceID)
+                    hvec.push_back(*h);
+            for(size_t i=0;i<hvec.size(); ++i)
+                mani.split_edge(hvec[i]);
+        }
+    }
+    void close_holes(Manifold& m)
+    {
+        for(HalfEdgeIDIterator h =  m.halfedges_begin(); h != m.halfedges_end(); ++h){
+            m.close_hole(*h);
+        }
+    }
+    void flip_orientation(Manifold& m)
+    {
+        vector<Vec3d> vertices;
+        VertexAttributeVector<int> idvec;
+        int k=0;
+        for(VertexID v: m.vertices()) {
+            idvec[v] = k++;
+            vertices.push_back(m.pos(v));
+        }
+        vector<int> indices;
+        vector<int> faces;
+        for(FaceID f: m.faces()) {
+            faces.push_back(no_edges(m, f));
+            circulate_face_cw(m, f, [&](VertexID v){
+                indices.push_back(idvec[v]);
+            });
+        }
+        m.clear();
+        m.build(vertices.size(), vertices[0].get(), faces.size(), &faces[0], &indices[0]);
+    }
+}

Subversion Repositories gelsvn

(root)/branches/GEL_github/src/HMesh/cleanup.cpp @ 677 – Rev 667 → 677