Subversion Repositories gelsvn

Rev

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

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