WebSVN – gelsvn – Diff – /branches/cpp11-devel/src/HMesh/subdivision.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 "dual.h"
#include "subdivision.h"
 
#include <vector>
#include "../CGLA/Vec3d.h"
 
#include "Manifold.h"
#include "AttributeVector.h"
 
namespace HMesh
{
    using namespace std;
    using namespace CGLA;
    
    void loop_split(Manifold& m_in, Manifold& m)
    {
        if(&m != &m_in)
            m = m_in;
        
        VertexAttributeVector<int> vtouched(m_in.allocated_vertices(), 0);
        
        vector<HalfEdgeID> hedges;
        for(HalfEdgeID h : m.halfedges())
            if(h<m.walker(h).opp().halfedge())
                hedges.push_back(h);
        
        for(HalfEdgeID h : hedges)
            vtouched[m.split_edge(h)] = 1;
        
        
        vector<FaceID> faces;
        for(FaceID fid : m.faces())
            faces.push_back(fid);
        for(FaceID fid : faces)
        {
            Walker w = m.walker(fid);
            
            if(vtouched[w.vertex()] == 0)
                w = w.next();
            
            assert(vtouched[w.vertex()] == 1);
            
            FaceID f = fid;
            VertexID v1, orig_vert = w.vertex();
            w = w.next();
            do
            {
                VertexID v0 = w.opp().vertex();
                w = w.next();
                v1  = w.vertex();
                w = w.next();
                assert(vtouched[v0] == 1);
                assert(vtouched[v1] == 1);
                f = m.split_face_by_edge(f, v0, v1);
            }
            while (v1 != orig_vert);
        }
    }
    
    void cc_split(Manifold& m_in, Manifold& m_out)
    {
        const int Invalid = -1;
 
        vector<Vec3d> new_points;
        new_points.reserve(m_in.no_vertices());
 
        VertexAttributeVector<int> vtouched(m_in.allocated_vertices(), Invalid);
        HalfEdgeAttributeVector<int> htouched(m_in.allocated_halfedges(), Invalid);
 
        int npsize = 0;
        for(VertexIDIterator v = m_in.vertices_begin(); v != m_in.vertices_end(); ++v){       
            vtouched[*v] = npsize;
            new_points.push_back(m_in.pos(*v));
            ++npsize;
        }
 
        for(HalfEdgeIDIterator h = m_in.halfedges_begin(); h != m_in.halfedges_end(); ++h){
            if(htouched[*h] != Invalid)
                continue;
 
            Walker w = m_in.walker(*h);
            htouched[*h] = htouched[w.opp().halfedge()] = npsize;
            new_points.push_back((m_in.pos(w.vertex()) + m_in.pos(w.opp().vertex())) * 0.5f);
            ++npsize;
 
        }
        vector<int> indices;
        vector<int> faces;
 
        for(FaceIDIterator f = m_in.faces_begin(); f != m_in.faces_end(); ++f){           
            for(Walker w = m_in.walker(*f); !w.full_circle(); w = w.circulate_face_cw()){
                indices.push_back(npsize);
                indices.push_back(htouched[w.halfedge()]);
                indices.push_back(vtouched[w.vertex()]);
                indices.push_back(htouched[w.next().halfedge()]);
                faces.push_back(4);
            }
            new_points.push_back(centre(m_in, *f));
            ++npsize;
        }
 
        m_out.clear();
        m_out.build(npsize, reinterpret_cast<double*>(&new_points[0]), faces.size(), &faces[0], &indices[0]);
    }
    
    void root3_subdivide(Manifold& m_in, Manifold& m)
    {
        if(&m != &m_in)
            m = m_in;
        
        VertexAttributeVector<int> vtouched(m.allocated_vertices(), 0);
        VertexAttributeVector<Vec3d> new_pos(m.allocated_vertices(), Vec3d(0));
        
        for (VertexIDIterator vid = m.vertices_begin(); vid != m.vertices_end(); ++vid) {
            int v = valency(m, *vid);
            double beta = (4.0 - 2.0 * cos(2.0 * M_PI / v))/(9.0*v);
            new_pos[*vid] = (1.0 - v * beta) * m.pos(*vid);
            for(Walker w = m.walker(*vid); !w.full_circle(); w = w.circulate_vertex_ccw())
            {
                new_pos[*vid] += beta * m.pos(w.vertex());
            }
        }
 
        vector<FaceID> faces;
        for(FaceIDIterator f = m.faces_begin(); f != m.faces_end(); ++f)
            faces.push_back(*f);
        for(int i=0;i<faces.size(); ++i)
            vtouched[m.split_face_by_vertex(faces[i])] = 1;
        
        for(HalfEdgeIDIterator h = m.halfedges_begin(); h != m.halfedges_end(); ++h)
        {
            Walker w = m.walker(*h);
            
            if(vtouched[w.vertex()] == 0 && vtouched[w.opp().vertex()] == 0 &&
                precond_flip_edge(m, *h))
                m.flip_edge(*h);
        }
        
        for (VertexIDIterator vid = m.vertices_begin(); vid != m.vertices_end(); ++vid)
            if(vtouched[*vid] == 0)
                m.pos(*vid) = new_pos[*vid];
    }
 
    void rootCC_subdivide(Manifold& m_in, Manifold& m)
    {
        if(&m != &m_in)
            m = m_in;
        
        VertexAttributeVector<int> vtouched(m.allocated_vertices(), 0);
        VertexAttributeVector<Vec3d> new_pos(m.allocated_vertices(), Vec3d(0));
        
        for (VertexIDIterator vid = m.vertices_begin(); vid != m.vertices_end(); ++vid) {
            int v = valency(m, *vid);
            double beta = (4.0 - 2.0 * cos(2.0 * M_PI / v))/(9.0*v);
            new_pos[*vid] = (1.0 - v * beta) * m.pos(*vid);
            for(Walker w = m.walker(*vid); !w.full_circle(); w = w.circulate_vertex_ccw())
            {
                new_pos[*vid] += beta * m.pos(w.vertex());
            }
        }
        
        vector<FaceID> faces;
        for(FaceIDIterator f = m.faces_begin(); f != m.faces_end(); ++f)
            faces.push_back(*f);
        for(int i=0;i<faces.size(); ++i)
            vtouched[m.split_face_by_vertex(faces[i])] = 1;
        
        for(HalfEdgeIDIterator h = m.halfedges_begin(); h != m.halfedges_end(); ++h)
        {
            Walker w = m.walker(*h);
            
            if(vtouched[w.vertex()] == 0 && vtouched[w.opp().vertex()] == 0)
                m.merge_faces(w.face(), *h);
        }
        
        for (VertexIDIterator vid = m.vertices_begin(); vid != m.vertices_end(); ++vid)
            if(vtouched[*vid] == 0)
                m.pos(*vid) = new_pos[*vid];
    }
 
    
    void butterfly_subdivide(Manifold& m_in, Manifold& m)
    {
        const float S[4][6] = {{5.0/12.0, -1.0/12.0, -1.0/12.0, 0, 0, 0},
            {3.0/8.0, 0, -1.0/8.0, 0, 0, 0},
            {0.35,
                0.03090169943749475,
                -0.08090169943749474,
                -0.08090169943749474,
                0.03090169943749468,0 },
            {0, 1.0f/8, -1.0f/8, 0, -1.0f/8, 1.0f/8}};
 
        
        if(&m != &m_in)
            m = m_in;
 
        HalfEdgeAttributeVector<Vec3d> new_vertices_pos(m.allocated_halfedges(), Vec3d(0.0));
        HalfEdgeAttributeVector<int> htouched(m.allocated_halfedges(), 0);
        VertexAttributeVector<int> vtouched(m.allocated_vertices(), 0);
 
        for(HalfEdgeIDIterator hid = m.halfedges_begin(); hid != m.halfedges_end(); ++hid)
            {
                Walker w = m.walker(*hid);
                VertexID v0 = w.opp().vertex();
                
                int K = valency(m, v0);
                int Ko = valency(m, w.vertex());
                if((K==6 && Ko==6)|| K != 6)
                {
                    Vec3d pos((K==6 ? 1.0 : 3.0/4.0) * m.pos(v0));
                    for(int k=0;k<K; ++k, w = w.circulate_vertex_ccw())
                    {
                        double s = (K<=6) ? S[K-3][k]:(0.25+cos((2.0*M_PI*k)/K)+0.5*cos((4.0*M_PI*k)/K))/K;
                        pos += s * m.pos(w.vertex());                        
                    }
                    new_vertices_pos[*hid] = pos;
                    htouched[*hid] = 1;
                }
            }
        loop_split(m, m);
 
        for(HalfEdgeIDIterator hid = m.halfedges_begin(); hid != m.halfedges_end(); ++hid)
            if(htouched[*hid])
            {
                VertexID v = m.walker(*hid).opp().vertex();
                if(vtouched[v] == 0)
                    m.pos(v) =Vec3d(0);
                m.pos(v) += new_vertices_pos[*hid];
                vtouched[v]+=1;
            }
        for(VertexIDIterator vid = m.vertices_begin(); vid != m.vertices_end(); ++vid)
            if(vtouched[*vid])
                m.pos(*vid) /= vtouched[*vid];
    }
    
    enum Subd {QUAD_SUBD, CC_SUBD, LOOP_SUBD, TRI_SUBD};
    
    void subd_smooth(Subd subd_method, Manifold& m)
    {
        
        VertexAttributeVector<Vec3d> new_vertices(m.no_vertices(), Vec3d(0,0,0));
        
        for(FaceIDIterator fid = m.faces_begin(); fid != m.faces_end(); ++fid)
        {
            FaceID f = *fid;
            for(Walker wlk = m.walker(f); !wlk.full_circle(); wlk = wlk.next())
            {
                float val = valency(m, wlk.vertex());
                float A,B;
                
                switch(subd_method)
                {
                    case QUAD_SUBD:
                        A = 1.0f / (4.0f * val);
                        B = 1.0f / (4.0f * val);
                        break;
                    case CC_SUBD:
                        A = (1.0f-3.0f/val)	* (1.0f/val);
                        B = sqr(1.0f/val);
                        break;
                    case TRI_SUBD:
                        A = 2.0f / (8.0f * val);
                        B = 3.0f / (8.0f * val);
                        break;
                    case LOOP_SUBD:
                        float w = 5.0f/8.0f - sqr(3.0f/8.0f + 0.25 * cos(2.0f*M_PI/val));
                        A = (1-2*w)/val;
                        B = w/val;
                        break;
                }
                for(Walker wlk2 = m.walker(f); !wlk2.full_circle(); wlk2 = wlk2.next())
                {
                    if(wlk2.vertex()==wlk.vertex())
                        new_vertices[wlk.vertex()] += A * m.pos(wlk2.vertex());
                    else
                        new_vertices[wlk.vertex()] += B * m.pos(wlk2.vertex());
                }
                
            }
        }
        for(VertexIDIterator vid = m.vertices_begin(); vid != m.vertices_end(); ++vid)
            m.pos(*vid) = new_vertices[*vid];
    }
 
    void cc_smooth(Manifold& m)
    {
        subd_smooth(CC_SUBD, m);
    }
    
    void loop_smooth(Manifold& m)
    {
        subd_smooth(LOOP_SUBD, m);
    }
 
 
}
 

Rev 631	Rev 666
1	`/* ----------------------------------------------------------------------- *`	1	`/* ----------------------------------------------------------------------- *`
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 "dual.h"`	7	`#include "dual.h"`
8	`#include "subdivision.h"`	8	`#include "subdivision.h"`
9		9
10	`#include <vector>`	10	`#include <vector>`
11	`#include "../CGLA/Vec3d.h"`	11	`#include "../CGLA/Vec3d.h"`
12		12
13	`#include "Manifold.h"`	13	`#include "Manifold.h"`
14	`#include "AttributeVector.h"`	14	`#include "AttributeVector.h"`
15		15
16	`namespace HMesh`	16	`namespace HMesh`
17	`{`	17	`{`
18	`using namespace std;`	18	`using namespace std;`
19	`using namespace CGLA;`	19	`using namespace CGLA;`
20		20
21	`void loop_split(Manifold& m_in, Manifold& m)`	21	`void loop_split(Manifold& m_in, Manifold& m)`
22	`{`	22	`{`
23	`if(&m != &m_in)`	23	`if(&m != &m_in)`
24	`m = m_in;`	24	`m = m_in;`
25		25
26	`VertexAttributeVector<int> vtouched(m_in.allocated_vertices(), 0);`	26	`VertexAttributeVector<int> vtouched(m_in.allocated_vertices(), 0);`
27		27
28	`vector<HalfEdgeID> hedges;`	28	`vector<HalfEdgeID> hedges;`
29	`for(HalfEdgeID h : m.halfedges())`	29	`for(HalfEdgeID h : m.halfedges())`
30	`if(h<m.walker(h).opp().halfedge())`	30	`if(h<m.walker(h).opp().halfedge())`
31	`hedges.push_back(h);`	31	`hedges.push_back(h);`
32		32
33	`for(HalfEdgeID h : hedges)`	33	`for(HalfEdgeID h : hedges)`
34	`vtouched[m.split_edge(h)] = 1;`	34	`vtouched[m.split_edge(h)] = 1;`
35		35
36		36
37	`vector<FaceID> faces;`	37	`vector<FaceID> faces;`
38	`for(FaceID fid : m.faces())`	38	`for(FaceID fid : m.faces())`
39	`faces.push_back(fid);`	39	`faces.push_back(fid);`
40	`for(FaceID fid : faces)`	40	`for(FaceID fid : faces)`
41	`{`	41	`{`
42	`Walker w = m.walker(fid);`	42	`Walker w = m.walker(fid);`
43		43
44	`if(vtouched[w.vertex()] == 0)`	44	`if(vtouched[w.vertex()] == 0)`
45	`w = w.next();`	45	`w = w.next();`
46		46
47	`assert(vtouched[w.vertex()] == 1);`	47	`assert(vtouched[w.vertex()] == 1);`
48		48
49	`FaceID f = fid;`	49	`FaceID f = fid;`
50	`VertexID v1, orig_vert = w.vertex();`	50	`VertexID v1, orig_vert = w.vertex();`
51	`w = w.next();`	51	`w = w.next();`
52	`do`	52	`do`
53	`{`	53	`{`
54	`VertexID v0 = w.opp().vertex();`	54	`VertexID v0 = w.opp().vertex();`
55	`w = w.next();`	55	`w = w.next();`
56	`v1 = w.vertex();`	56	`v1 = w.vertex();`
57	`w = w.next();`	57	`w = w.next();`
58	`assert(vtouched[v0] == 1);`	58	`assert(vtouched[v0] == 1);`
59	`assert(vtouched[v1] == 1);`	59	`assert(vtouched[v1] == 1);`
60	`f = m.split_face_by_edge(f, v0, v1);`	60	`f = m.split_face_by_edge(f, v0, v1);`
61	`}`	61	`}`
62	`while (v1 != orig_vert);`	62	`while (v1 != orig_vert);`
63	`}`	63	`}`
64	`}`	64	`}`
65		65
66	`void cc_split(Manifold& m_in, Manifold& m_out)`	66	`void cc_split(Manifold& m_in, Manifold& m_out)`
67	`{`	67	`{`
68	`const int Invalid = -1;`	68	`const int Invalid = -1;`
69		69
70	`vector<Vec3d> new_points;`	70	`vector<Vec3d> new_points;`
71	`new_points.reserve(m_in.no_vertices());`	71	`new_points.reserve(m_in.no_vertices());`
72		72
73	`VertexAttributeVector<int> vtouched(m_in.allocated_vertices(), Invalid);`	73	`VertexAttributeVector<int> vtouched(m_in.allocated_vertices(), Invalid);`
74	`HalfEdgeAttributeVector<int> htouched(m_in.allocated_halfedges(), Invalid);`	74	`HalfEdgeAttributeVector<int> htouched(m_in.allocated_halfedges(), Invalid);`
75		75
76	`int npsize = 0;`	76	`int npsize = 0;`
77	`for(VertexIDIterator v = m_in.vertices_begin(); v != m_in.vertices_end(); ++v){`	77	`for(VertexIDIterator v = m_in.vertices_begin(); v != m_in.vertices_end(); ++v){`
78	`vtouched[*v] = npsize;`	78	`vtouched[*v] = npsize;`
79	`new_points.push_back(m_in.pos(*v));`	79	`new_points.push_back(m_in.pos(*v));`
80	`++npsize;`	80	`++npsize;`
81	`}`	81	`}`
82		82
83	`for(HalfEdgeIDIterator h = m_in.halfedges_begin(); h != m_in.halfedges_end(); ++h){`	83	`for(HalfEdgeIDIterator h = m_in.halfedges_begin(); h != m_in.halfedges_end(); ++h){`
84	`if(htouched[*h] != Invalid)`	84	`if(htouched[*h] != Invalid)`
85	`continue;`	85	`continue;`
86		86
87	`Walker w = m_in.walker(*h);`	87	`Walker w = m_in.walker(*h);`
88	`htouched[*h] = htouched[w.opp().halfedge()] = npsize;`	88	`htouched[*h] = htouched[w.opp().halfedge()] = npsize;`
89	`new_points.push_back((m_in.pos(w.vertex()) + m_in.pos(w.opp().vertex())) * 0.5f);`	89	`new_points.push_back((m_in.pos(w.vertex()) + m_in.pos(w.opp().vertex())) * 0.5f);`
90	`++npsize;`	90	`++npsize;`
91		91
92	`}`	92	`}`
93	`vector<int> indices;`	93	`vector<int> indices;`
94	`vector<int> faces;`	94	`vector<int> faces;`
95		95
96	`for(FaceIDIterator f = m_in.faces_begin(); f != m_in.faces_end(); ++f){`	96	`for(FaceIDIterator f = m_in.faces_begin(); f != m_in.faces_end(); ++f){`
97	`for(Walker w = m_in.walker(*f); !w.full_circle(); w = w.circulate_face_cw()){`	97	`for(Walker w = m_in.walker(*f); !w.full_circle(); w = w.circulate_face_cw()){`
98	`indices.push_back(npsize);`	98	`indices.push_back(npsize);`
99	`indices.push_back(htouched[w.halfedge()]);`	99	`indices.push_back(htouched[w.halfedge()]);`
100	`indices.push_back(vtouched[w.vertex()]);`	100	`indices.push_back(vtouched[w.vertex()]);`
101	`indices.push_back(htouched[w.next().halfedge()]);`	101	`indices.push_back(htouched[w.next().halfedge()]);`
102	`faces.push_back(4);`	102	`faces.push_back(4);`
103	`}`	103	`}`
104	`new_points.push_back(centre(m_in, *f));`	104	`new_points.push_back(centre(m_in, *f));`
105	`++npsize;`	105	`++npsize;`
106	`}`	106	`}`
107		107
108	`m_out.clear();`	108	`m_out.clear();`
109	`m_out.build(npsize, reinterpret_cast<double*>(&new_points[0]), faces.size(), &faces[0], &indices[0]);`	109	`m_out.build(npsize, reinterpret_cast<double*>(&new_points[0]), faces.size(), &faces[0], &indices[0]);`
110	`}`	110	`}`
111		111
112	`void root3_subdivide(Manifold& m_in, Manifold& m)`	112	`void root3_subdivide(Manifold& m_in, Manifold& m)`
113	`{`	113	`{`
114	`if(&m != &m_in)`	114	`if(&m != &m_in)`
115	`m = m_in;`	115	`m = m_in;`
116		116
117	`VertexAttributeVector<int> vtouched(m.allocated_vertices(), 0);`	117	`VertexAttributeVector<int> vtouched(m.allocated_vertices(), 0);`
118	`VertexAttributeVector<Vec3d> new_pos(m.allocated_vertices(), Vec3d(0));`	118	`VertexAttributeVector<Vec3d> new_pos(m.allocated_vertices(), Vec3d(0));`
119		119
120	`for (VertexIDIterator vid = m.vertices_begin(); vid != m.vertices_end(); ++vid) {`	120	`for (VertexIDIterator vid = m.vertices_begin(); vid != m.vertices_end(); ++vid) {`
121	`int v = valency(m, *vid);`	121	`int v = valency(m, *vid);`
122	`double beta = (4.0 - 2.0 * cos(2.0 * M_PI / v))/(9.0*v);`	122	`double beta = (4.0 - 2.0 * cos(2.0 * M_PI / v))/(9.0*v);`
123	`new_pos[vid] = (1.0 - v beta) * m.pos(*vid);`	123	`new_pos[vid] = (1.0 - v beta) * m.pos(*vid);`
124	`for(Walker w = m.walker(*vid); !w.full_circle(); w = w.circulate_vertex_ccw())`	124	`for(Walker w = m.walker(*vid); !w.full_circle(); w = w.circulate_vertex_ccw())`
125	`{`	125	`{`
126	`new_pos[vid] += beta m.pos(w.vertex());`	126	`new_pos[vid] += beta m.pos(w.vertex());`
127	`}`	127	`}`
128	`}`	128	`}`
129		129
130	`vector<FaceID> faces;`	130	`vector<FaceID> faces;`
131	`for(FaceIDIterator f = m.faces_begin(); f != m.faces_end(); ++f)`	131	`for(FaceIDIterator f = m.faces_begin(); f != m.faces_end(); ++f)`
132	`faces.push_back(*f);`	132	`faces.push_back(*f);`
133	`for(int i=0;i<faces.size(); ++i)`	133	`for(int i=0;i<faces.size(); ++i)`
134	`vtouched[m.split_face_by_vertex(faces[i])] = 1;`	134	`vtouched[m.split_face_by_vertex(faces[i])] = 1;`
135		135
136	`for(HalfEdgeIDIterator h = m.halfedges_begin(); h != m.halfedges_end(); ++h)`	136	`for(HalfEdgeIDIterator h = m.halfedges_begin(); h != m.halfedges_end(); ++h)`
137	`{`	137	`{`
138	`Walker w = m.walker(*h);`	138	`Walker w = m.walker(*h);`
139		139
140	`if(vtouched[w.vertex()] == 0 && vtouched[w.opp().vertex()] == 0 &&`	140	`if(vtouched[w.vertex()] == 0 && vtouched[w.opp().vertex()] == 0 &&`
141	`precond_flip_edge(m, *h))`	141	`precond_flip_edge(m, *h))`
142	`m.flip_edge(*h);`	142	`m.flip_edge(*h);`
143	`}`	143	`}`
144		144
145	`for (VertexIDIterator vid = m.vertices_begin(); vid != m.vertices_end(); ++vid)`	145	`for (VertexIDIterator vid = m.vertices_begin(); vid != m.vertices_end(); ++vid)`
146	`if(vtouched[*vid] == 0)`	146	`if(vtouched[*vid] == 0)`
147	`m.pos(vid) = new_pos[vid];`	147	`m.pos(vid) = new_pos[vid];`
148	`}`	148	`}`
-		149
-		150	`void rootCC_subdivide(Manifold& m_in, Manifold& m)`
149		151	`{`
-		152	`if(&m != &m_in)`
-		153	`m = m_in;`
-		154
-		155	`VertexAttributeVector<int> vtouched(m.allocated_vertices(), 0);`
-		156	`VertexAttributeVector<Vec3d> new_pos(m.allocated_vertices(), Vec3d(0));`
-		157
-		158	`for (VertexIDIterator vid = m.vertices_begin(); vid != m.vertices_end(); ++vid) {`
-		159	`int v = valency(m, *vid);`
-		160	`double beta = (4.0 - 2.0 * cos(2.0 * M_PI / v))/(9.0*v);`
-		161	`new_pos[vid] = (1.0 - v beta) * m.pos(*vid);`
-		162	`for(Walker w = m.walker(*vid); !w.full_circle(); w = w.circulate_vertex_ccw())`
-		163	`{`
-		164	`new_pos[vid] += beta m.pos(w.vertex());`
-		165	`}`
-		166	`}`
-		167
-		168	`vector<FaceID> faces;`
-		169	`for(FaceIDIterator f = m.faces_begin(); f != m.faces_end(); ++f)`
-		170	`faces.push_back(*f);`
-		171	`for(int i=0;i<faces.size(); ++i)`
-		172	`vtouched[m.split_face_by_vertex(faces[i])] = 1;`
-		173
-		174	`for(HalfEdgeIDIterator h = m.halfedges_begin(); h != m.halfedges_end(); ++h)`
-		175	`{`
-		176	`Walker w = m.walker(*h);`
-		177
-		178	`if(vtouched[w.vertex()] == 0 && vtouched[w.opp().vertex()] == 0)`
-		179	`m.merge_faces(w.face(), *h);`
-		180	`}`
-		181
-		182	`for (VertexIDIterator vid = m.vertices_begin(); vid != m.vertices_end(); ++vid)`
-		183	`if(vtouched[*vid] == 0)`
-		184	`m.pos(vid) = new_pos[vid];`
-		185	`}`
-		186
150		187
151	`void butterfly_subdivide(Manifold& m_in, Manifold& m)`	188	`void butterfly_subdivide(Manifold& m_in, Manifold& m)`
152	`{`	189	`{`
153	`const float S[4][6] = {{5.0/12.0, -1.0/12.0, -1.0/12.0, 0, 0, 0},`	190	`const float S[4][6] = {{5.0/12.0, -1.0/12.0, -1.0/12.0, 0, 0, 0},`
154	`{3.0/8.0, 0, -1.0/8.0, 0, 0, 0},`	191	`{3.0/8.0, 0, -1.0/8.0, 0, 0, 0},`
155	`{0.35,`	192	`{0.35,`
156	`0.03090169943749475,`	193	`0.03090169943749475,`
157	`-0.08090169943749474,`	194	`-0.08090169943749474,`
158	`-0.08090169943749474,`	195	`-0.08090169943749474,`
159	`0.03090169943749468,0 },`	196	`0.03090169943749468,0 },`
160	`{0, 1.0f/8, -1.0f/8, 0, -1.0f/8, 1.0f/8}};`	197	`{0, 1.0f/8, -1.0f/8, 0, -1.0f/8, 1.0f/8}};`
161		198
162		199
163	`if(&m != &m_in)`	200	`if(&m != &m_in)`
164	`m = m_in;`	201	`m = m_in;`
165		202
166	`HalfEdgeAttributeVector<Vec3d> new_vertices_pos(m.allocated_halfedges(), Vec3d(0.0));`	203	`HalfEdgeAttributeVector<Vec3d> new_vertices_pos(m.allocated_halfedges(), Vec3d(0.0));`
167	`HalfEdgeAttributeVector<int> htouched(m.allocated_halfedges(), 0);`	204	`HalfEdgeAttributeVector<int> htouched(m.allocated_halfedges(), 0);`
168	`VertexAttributeVector<int> vtouched(m.allocated_vertices(), 0);`	205	`VertexAttributeVector<int> vtouched(m.allocated_vertices(), 0);`
169		206
170	`for(HalfEdgeIDIterator hid = m.halfedges_begin(); hid != m.halfedges_end(); ++hid)`	207	`for(HalfEdgeIDIterator hid = m.halfedges_begin(); hid != m.halfedges_end(); ++hid)`
171	`{`	208	`{`
172	`Walker w = m.walker(*hid);`	209	`Walker w = m.walker(*hid);`
173	`VertexID v0 = w.opp().vertex();`	210	`VertexID v0 = w.opp().vertex();`
174		211
175	`int K = valency(m, v0);`	212	`int K = valency(m, v0);`
176	`int Ko = valency(m, w.vertex());`	213	`int Ko = valency(m, w.vertex());`
177	`if((K==6 && Ko==6)\|\| K != 6)`	214	`if((K==6 && Ko==6)\|\| K != 6)`
178	`{`	215	`{`
179	`Vec3d pos((K==6 ? 1.0 : 3.0/4.0) * m.pos(v0));`	216	`Vec3d pos((K==6 ? 1.0 : 3.0/4.0) * m.pos(v0));`
180	`for(int k=0;k<K; ++k, w = w.circulate_vertex_ccw())`	217	`for(int k=0;k<K; ++k, w = w.circulate_vertex_ccw())`
181	`{`	218	`{`
182	`double s = (K<=6) ? S[K-3][k]:(0.25+cos((2.0M_PIk)/K)+0.5cos((4.0M_PI*k)/K))/K;`	219	`double s = (K<=6) ? S[K-3][k]:(0.25+cos((2.0M_PIk)/K)+0.5cos((4.0M_PI*k)/K))/K;`
183	`pos += s * m.pos(w.vertex());`	220	`pos += s * m.pos(w.vertex());`
184	`}`	221	`}`
185	`new_vertices_pos[*hid] = pos;`	222	`new_vertices_pos[*hid] = pos;`
186	`htouched[*hid] = 1;`	223	`htouched[*hid] = 1;`
187	`}`	224	`}`
188	`}`	225	`}`
189	`loop_split(m, m);`	226	`loop_split(m, m);`
190		227
191	`for(HalfEdgeIDIterator hid = m.halfedges_begin(); hid != m.halfedges_end(); ++hid)`	228	`for(HalfEdgeIDIterator hid = m.halfedges_begin(); hid != m.halfedges_end(); ++hid)`
192	`if(htouched[*hid])`	229	`if(htouched[*hid])`
193	`{`	230	`{`
194	`VertexID v = m.walker(*hid).opp().vertex();`	231	`VertexID v = m.walker(*hid).opp().vertex();`
195	`if(vtouched[v] == 0)`	232	`if(vtouched[v] == 0)`
196	`m.pos(v) =Vec3d(0);`	233	`m.pos(v) =Vec3d(0);`
197	`m.pos(v) += new_vertices_pos[*hid];`	234	`m.pos(v) += new_vertices_pos[*hid];`
198	`vtouched[v]+=1;`	235	`vtouched[v]+=1;`
199	`}`	236	`}`
200	`for(VertexIDIterator vid = m.vertices_begin(); vid != m.vertices_end(); ++vid)`	237	`for(VertexIDIterator vid = m.vertices_begin(); vid != m.vertices_end(); ++vid)`
201	`if(vtouched[*vid])`	238	`if(vtouched[*vid])`
202	`m.pos(vid) /= vtouched[vid];`	239	`m.pos(vid) /= vtouched[vid];`
203	`}`	240	`}`
204		241
205	`enum Subd {QUAD_SUBD, CC_SUBD, LOOP_SUBD, TRI_SUBD};`	242	`enum Subd {QUAD_SUBD, CC_SUBD, LOOP_SUBD, TRI_SUBD};`
206		243
207	`void subd_smooth(Subd subd_method, Manifold& m)`	244	`void subd_smooth(Subd subd_method, Manifold& m)`
208	`{`	245	`{`
209		246
210	`VertexAttributeVector<Vec3d> new_vertices(m.no_vertices(), Vec3d(0,0,0));`	247	`VertexAttributeVector<Vec3d> new_vertices(m.no_vertices(), Vec3d(0,0,0));`
211		248
212	`for(FaceIDIterator fid = m.faces_begin(); fid != m.faces_end(); ++fid)`	249	`for(FaceIDIterator fid = m.faces_begin(); fid != m.faces_end(); ++fid)`
213	`{`	250	`{`
214	`FaceID f = *fid;`	251	`FaceID f = *fid;`
215	`for(Walker wlk = m.walker(f); !wlk.full_circle(); wlk = wlk.next())`	252	`for(Walker wlk = m.walker(f); !wlk.full_circle(); wlk = wlk.next())`
216	`{`	253	`{`
217	`float val = valency(m, wlk.vertex());`	254	`float val = valency(m, wlk.vertex());`
218	`float A,B;`	255	`float A,B;`
219		256
220	`switch(subd_method)`	257	`switch(subd_method)`
221	`{`	258	`{`
222	`case QUAD_SUBD:`	259	`case QUAD_SUBD:`
223	`A = 1.0f / (4.0f * val);`	260	`A = 1.0f / (4.0f * val);`
224	`B = 1.0f / (4.0f * val);`	261	`B = 1.0f / (4.0f * val);`
225	`break;`	262	`break;`
226	`case CC_SUBD:`	263	`case CC_SUBD:`
227	`A = (1.0f-3.0f/val) * (1.0f/val);`	264	`A = (1.0f-3.0f/val) * (1.0f/val);`
228	`B = sqr(1.0f/val);`	265	`B = sqr(1.0f/val);`
229	`break;`	266	`break;`
230	`case TRI_SUBD:`	267	`case TRI_SUBD:`
231	`A = 2.0f / (8.0f * val);`	268	`A = 2.0f / (8.0f * val);`
232	`B = 3.0f / (8.0f * val);`	269	`B = 3.0f / (8.0f * val);`
233	`break;`	270	`break;`
234	`case LOOP_SUBD:`	271	`case LOOP_SUBD:`
235	`float w = 5.0f/8.0f - sqr(3.0f/8.0f + 0.25 * cos(2.0f*M_PI/val));`	272	`float w = 5.0f/8.0f - sqr(3.0f/8.0f + 0.25 * cos(2.0f*M_PI/val));`
236	`A = (1-2*w)/val;`	273	`A = (1-2*w)/val;`
237	`B = w/val;`	274	`B = w/val;`
238	`break;`	275	`break;`
239	`}`	276	`}`
240	`for(Walker wlk2 = m.walker(f); !wlk2.full_circle(); wlk2 = wlk2.next())`	277	`for(Walker wlk2 = m.walker(f); !wlk2.full_circle(); wlk2 = wlk2.next())`
241	`{`	278	`{`
242	`if(wlk2.vertex()==wlk.vertex())`	279	`if(wlk2.vertex()==wlk.vertex())`
243	`new_vertices[wlk.vertex()] += A * m.pos(wlk2.vertex());`	280	`new_vertices[wlk.vertex()] += A * m.pos(wlk2.vertex());`
244	`else`	281	`else`
245	`new_vertices[wlk.vertex()] += B * m.pos(wlk2.vertex());`	282	`new_vertices[wlk.vertex()] += B * m.pos(wlk2.vertex());`
246	`}`	283	`}`
247		284
248	`}`	285	`}`
249	`}`	286	`}`
250	`for(VertexIDIterator vid = m.vertices_begin(); vid != m.vertices_end(); ++vid)`	287	`for(VertexIDIterator vid = m.vertices_begin(); vid != m.vertices_end(); ++vid)`
251	`m.pos(vid) = new_vertices[vid];`	288	`m.pos(vid) = new_vertices[vid];`
252	`}`	289	`}`
253		290
254	`void cc_smooth(Manifold& m)`	291	`void cc_smooth(Manifold& m)`
255	`{`	292	`{`
256	`subd_smooth(CC_SUBD, m);`	293	`subd_smooth(CC_SUBD, m);`
257	`}`	294	`}`
258		295
259	`void loop_smooth(Manifold& m)`	296	`void loop_smooth(Manifold& m)`
260	`{`	297	`{`
261	`subd_smooth(LOOP_SUBD, m);`	298	`subd_smooth(LOOP_SUBD, m);`
262	`}`	299	`}`
263		300
264		301
265	`}`	302	`}`
266		303

Subversion Repositories gelsvn

(root)/branches/cpp11-devel/src/HMesh/subdivision.cpp – Rev 631 → 666