WebSVN – gelsvn – Diff – /branches/cpp11-devel/src/HMesh/curvature.cpp



        {
            for(int iter=0;iter<smooth_steps;++iter){
                VertexAttributeVector<T> new_vec(m.allocated_vertices());
                for(VertexIDIterator v = m.vertices_begin(); v != m.vertices_end(); ++v){
                    new_vec[*v] = vec[*v];
                    for(Walker w = m.walker(*v); !w.full_circle(); w = w.circulate_vertex_cw()){
                        new_vec[*v] += vec[w.vertex()];
                    }
                    new_vec[*v] /= (valency(m, *v) + 1.0);
                }
                swap(vec,new_vec);

 
    double voronoi_area(const Manifold& m, VertexID v)
    {
        double area_mixed = 0;
        //For each triangle T from the 1-ring neighborhood of x
        for(Walker w = m.walker(v); !w.full_circle(); w = w.circulate_vertex_cw()){
            double f_area = area(m, w.face());
 
            Vec3d v0(m.pos(v));
            Vec3d v1(m.pos(w.vertex()));
            Vec3d v2(m.pos(w.next().vertex()));

 
    double barycentric_area(const Manifold& m, VertexID v)
    {
        double barea = 0;
        //For each triangle T from the 1-ring neighborhood of x
        for(Walker w = m.walker(v); !w.full_circle(); w = w.circulate_vertex_cw()){
            barea += area(m, w.face())/3.0;
        }
        return barea;
    }
 

            return;
 
        Vec3d vertex(m.pos(v));
        curv_normal = Vec3d(0);
        w_sum = 0;
        for(Walker walker = m.walker(v); !walker.full_circle(); walker = walker.circulate_vertex_ccw()){
            Vec3d nbr(m.pos(walker.vertex()));
            Vec3d left(m.pos(walker.next().vertex()));
            Vec3d right(m.pos(walker.opp().next().vertex()));
 
            double d_left = dot(cond_normalize(nbr-left),cond_normalize(vertex-left));

    double gaussian_curvature_angle_defect(const Manifold& m, VertexID v)
    {
        if(boundary(m, v))
            return 0;
 
        Vec3d vertex(m.pos(v));
        vector<Vec3d> edges;
        for(Walker w = m.walker(v); !w.full_circle(); w = w.circulate_vertex_cw()){
            Vec3d e(normalize(m.pos(w.vertex()) - vertex));
            edges.push_back(e);
        }
        size_t N=edges.size();
        double angle_sum = 0;

    Mat3x3d curvature_tensor(const Manifold& m, HalfEdgeID h)
    {
        if(boundary(m, h))
            return Mat3x3d(0);
 
        Walker w = m.walker(h);
        Vec3d edge(m.pos(w.vertex()) - m.pos(w.opp().vertex()));
        double edge_len = length(edge);
        edge /= edge_len;
 
        Vec3d h_norm(normal(m, w.face()));

        Mat3x3d curv_tensor(0);
 
        if(boundary(m, v))
            return curv_tensor;
 
        for(Walker w = m.walker(v); !w.full_circle(); w = w.circulate_vertex_cw())
            curv_tensor += 0.5*curvature_tensor(m, w.halfedge());
 
        curv_tensor /= voronoi_area(m, v);
 
        return curv_tensor;

        orthogonal(Norm,X,Y);
        frame = Mat3x3d(X,Y,Norm);
        Vec3d centre(m.pos(v));
 
        vector<Vec3d> points;
        for(Walker w = m.walker(v); !w.full_circle(); w = w.circulate_vertex_cw())
            points.push_back(Vec3d(m.pos(w.vertex())));
 
        size_t N = points.size();
 
        CVector b(N);

            if(boundary(m, *v))
                continue;
            double a = voronoi_area(m, *v);
            tmp_curvature_tensors[*v] = curvature_tensors[*v] * a;
            tmp_area = a;
            for(Walker w = m.walker(*v); !w.full_circle(); w = w.circulate_vertex_cw()){
                if(!boundary(m, w.vertex())){
                    double a = voronoi_area(m, w.vertex());
                    tmp_curvature_tensors[*v] += curvature_tensors[w.vertex()]*a;
                    tmp_area += a;
                }

Rev 586	Rev 587
Line 45...	Line 45...
45	`{`	45	`{`
46	`for(int iter=0;iter<smooth_steps;++iter){`	46	`for(int iter=0;iter<smooth_steps;++iter){`
47	`VertexAttributeVector<T> new_vec(m.allocated_vertices());`	47	`VertexAttributeVector<T> new_vec(m.allocated_vertices());`
48	`for(VertexIDIterator v = m.vertices_begin(); v != m.vertices_end(); ++v){`	48	`for(VertexIDIterator v = m.vertices_begin(); v != m.vertices_end(); ++v){`
49	`new_vec[v] = vec[v];`	49	`new_vec[v] = vec[v];`
50	`for(HalfEdgeWalker w = m.halfedgewalker(*v); !w.full_circle(); w = w.circulate_vertex_cw()){`	50	`for(Walker w = m.walker(*v); !w.full_circle(); w = w.circulate_vertex_cw()){`
51	`new_vec[*v] += vec[w.vertex()];`	51	`new_vec[*v] += vec[w.vertex()];`
52	`}`	52	`}`
53	`new_vec[v] /= (valency(m, v) + 1.0);`	53	`new_vec[v] /= (valency(m, v) + 1.0);`
54	`}`	54	`}`
55	`swap(vec,new_vec);`	55	`swap(vec,new_vec);`
Line 59...	Line 59...
59		59
60	`double voronoi_area(const Manifold& m, VertexID v)`	60	`double voronoi_area(const Manifold& m, VertexID v)`
61	`{`	61	`{`
62	`double area_mixed = 0;`	62	`double area_mixed = 0;`
63	`//For each triangle T from the 1-ring neighborhood of x`	63	`//For each triangle T from the 1-ring neighborhood of x`
64	`for(HalfEdgeWalker w = m.halfedgewalker(v); !w.full_circle(); w = w.circulate_vertex_cw()){`	64	`for(Walker w = m.walker(v); !w.full_circle(); w = w.circulate_vertex_cw()){`
65	`double f_area = area(m, w.face());`	65	`double f_area = area(m, w.face());`
66		66
67	`Vec3d v0(m.pos(v));`	67	`Vec3d v0(m.pos(v));`
68	`Vec3d v1(m.pos(w.vertex()));`	68	`Vec3d v1(m.pos(w.vertex()));`
69	`Vec3d v2(m.pos(w.next().vertex()));`	69	`Vec3d v2(m.pos(w.next().vertex()));`
Line 93...	Line 93...
93		93
94	`double barycentric_area(const Manifold& m, VertexID v)`	94	`double barycentric_area(const Manifold& m, VertexID v)`
95	`{`	95	`{`
96	`double barea = 0;`	96	`double barea = 0;`
97	`//For each triangle T from the 1-ring neighborhood of x`	97	`//For each triangle T from the 1-ring neighborhood of x`
98	`for(HalfEdgeWalker w = m.halfedgewalker(v); !w.full_circle(); w = w.circulate_vertex_cw()){`	98	`for(Walker w = m.walker(v); !w.full_circle(); w = w.circulate_vertex_cw()){`
99	`barea += area(m, w.face())/3.0;`	99	`barea += area(m, w.face())/3.0;`
100	`}`	100	`}`
101	`return barea;`	101	`return barea;`
102	`}`	102	`}`
103		103
Line 107...	Line 107...
107	`return;`	107	`return;`
108		108
109	`Vec3d vertex(m.pos(v));`	109	`Vec3d vertex(m.pos(v));`
110	`curv_normal = Vec3d(0);`	110	`curv_normal = Vec3d(0);`
111	`w_sum = 0;`	111	`w_sum = 0;`
112	`for(HalfEdgeWalker walker = m.halfedgewalker(v); !walker.full_circle(); walker = walker.circulate_vertex_ccw()){`	112	`for(Walker walker = m.walker(v); !walker.full_circle(); walker = walker.circulate_vertex_ccw()){`
113	`Vec3d nbr(m.pos(walker.vertex()));`	113	`Vec3d nbr(m.pos(walker.vertex()));`
114	`Vec3d left(m.pos(walker.next().vertex()));`	114	`Vec3d left(m.pos(walker.next().vertex()));`
115	`Vec3d right(m.pos(walker.opp().next().vertex()));`	115	`Vec3d right(m.pos(walker.opp().next().vertex()));`
116		116
117	`double d_left = dot(cond_normalize(nbr-left),cond_normalize(vertex-left));`	117	`double d_left = dot(cond_normalize(nbr-left),cond_normalize(vertex-left));`
Line 152...	Line 152...
152	`double gaussian_curvature_angle_defect(const Manifold& m, VertexID v)`	152	`double gaussian_curvature_angle_defect(const Manifold& m, VertexID v)`
153	`{`	153	`{`
154	`if(boundary(m, v))`	154	`if(boundary(m, v))`
155	`return 0;`	155	`return 0;`
156		156
157	`Vec3f vertex(m.pos(v));`	157	`Vec3d vertex(m.pos(v));`
158	`vector<Vec3d> edges;`	158	`vector<Vec3d> edges;`
159	`for(HalfEdgeWalker w = m.halfedgewalker(v); !w.full_circle(); w = w.circulate_vertex_cw()){`	159	`for(Walker w = m.walker(v); !w.full_circle(); w = w.circulate_vertex_cw()){`
160	`Vec3d e(normalize(m.pos(w.vertex()) - vertex));`	160	`Vec3d e(normalize(m.pos(w.vertex()) - vertex));`
161	`edges.push_back(e);`	161	`edges.push_back(e);`
162	`}`	162	`}`
163	`size_t N=edges.size();`	163	`size_t N=edges.size();`
164	`double angle_sum = 0;`	164	`double angle_sum = 0;`
Line 175...	Line 175...
175	`Mat3x3d curvature_tensor(const Manifold& m, HalfEdgeID h)`	175	`Mat3x3d curvature_tensor(const Manifold& m, HalfEdgeID h)`
176	`{`	176	`{`
177	`if(boundary(m, h))`	177	`if(boundary(m, h))`
178	`return Mat3x3d(0);`	178	`return Mat3x3d(0);`
179		179
180	`HalfEdgeWalker w = m.halfedgewalker(h);`	180	`Walker w = m.walker(h);`
181	`Vec3d edge(m.pos(w.vertex()) - m.pos(w.opp().vertex()));`	181	`Vec3d edge(m.pos(w.vertex()) - m.pos(w.opp().vertex()));`
182	`double edge_len = length(edge);`	182	`double edge_len = length(edge);`
183	`edge /= edge_len;`	183	`edge /= edge_len;`
184		184
185	`Vec3d h_norm(normal(m, w.face()));`	185	`Vec3d h_norm(normal(m, w.face()));`
Line 200...	Line 200...
200	`Mat3x3d curv_tensor(0);`	200	`Mat3x3d curv_tensor(0);`
201		201
202	`if(boundary(m, v))`	202	`if(boundary(m, v))`
203	`return curv_tensor;`	203	`return curv_tensor;`
204		204
205	`for(HalfEdgeWalker w = m.halfedgewalker(v); !w.full_circle(); w = w.circulate_vertex_cw())`	205	`for(Walker w = m.walker(v); !w.full_circle(); w = w.circulate_vertex_cw())`
206	`curv_tensor += 0.5*curvature_tensor(m, w.halfedge());`	206	`curv_tensor += 0.5*curvature_tensor(m, w.halfedge());`
207		207
208	`curv_tensor /= voronoi_area(m, v);`	208	`curv_tensor /= voronoi_area(m, v);`
209		209
210	`return curv_tensor;`	210	`return curv_tensor;`
Line 222...	Line 222...
222	`orthogonal(Norm,X,Y);`	222	`orthogonal(Norm,X,Y);`
223	`frame = Mat3x3d(X,Y,Norm);`	223	`frame = Mat3x3d(X,Y,Norm);`
224	`Vec3d centre(m.pos(v));`	224	`Vec3d centre(m.pos(v));`
225		225
226	`vector<Vec3d> points;`	226	`vector<Vec3d> points;`
227	`for(HalfEdgeWalker w = m.halfedgewalker(v); !w.full_circle(); w = w.circulate_vertex_cw())`	227	`for(Walker w = m.walker(v); !w.full_circle(); w = w.circulate_vertex_cw())`
228	`points.push_back(Vec3d(m.pos(w.vertex())));`	228	`points.push_back(Vec3d(m.pos(w.vertex())));`
229		229
230	`size_t N = points.size();`	230	`size_t N = points.size();`
231		231
232	`CVector b(N);`	232	`CVector b(N);`
Line 268...	Line 268...
268	`if(boundary(m, *v))`	268	`if(boundary(m, *v))`
269	`continue;`	269	`continue;`
270	`double a = voronoi_area(m, *v);`	270	`double a = voronoi_area(m, *v);`
271	`tmp_curvature_tensors[v] = curvature_tensors[v] * a;`	271	`tmp_curvature_tensors[v] = curvature_tensors[v] * a;`
272	`tmp_area = a;`	272	`tmp_area = a;`
273	`for(HalfEdgeWalker w = m.halfedgewalker(*v); !w.full_circle(); w = w.circulate_vertex_cw()){`	273	`for(Walker w = m.walker(*v); !w.full_circle(); w = w.circulate_vertex_cw()){`
274	`if(!boundary(m, w.vertex())){`	274	`if(!boundary(m, w.vertex())){`
275	`double a = voronoi_area(m, w.vertex());`	275	`double a = voronoi_area(m, w.vertex());`
276	`tmp_curvature_tensors[v] += curvature_tensors[w.vertex()]a;`	276	`tmp_curvature_tensors[v] += curvature_tensors[w.vertex()]a;`
277	`tmp_area += a;`	277	`tmp_area += a;`
278	`}`	278	`}`

Subversion Repositories gelsvn

(root)/branches/cpp11-devel/src/HMesh/curvature.cpp @ 652 – Rev 586 → 587