WebSVN – gelsvn – Diff – /branches/cpp11-devel/src/HMesh/curvature.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 "curvature.h"
 
#include "../GLGraphics/gel_glut.h"
 
#include <iostream>
#include "../CGLA/eigensolution.h"
#include "../CGLA/Vec2d.h"
#include "../CGLA/Vec3d.h"
#include "../CGLA/Mat3x3d.h"
#include "../CGLA/Mat2x2d.h"
#include "../CGLA/Mat2x3d.h"
 
#include "Manifold.h"
#include "AttributeVector.h"
#include "x3d_save.h"
#include "x3d_load.h"
#include "obj_load.h"
#include "mesh_optimization.h"
 
#include "../LinAlg/Matrix.h"
#include "../LinAlg/Vector.h"
#include "../LinAlg/LapackFunc.h"
 
using namespace std;
 
using namespace LinAlg;
using namespace CGLA;
using namespace HMesh;
 
namespace HMesh
{
    namespace 
    {
        //double scal = 0.001;
        //double vector_scal = 0.001;
 
        template<class T> 
        void smooth_something_on_mesh(const Manifold& m, VertexAttributeVector<T>& vec, int smooth_steps)
        {
            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 a0 = acos(dot(v1-v0, v2-v0)/(length(v1-v0)*length(v2-v0)));
            double a1 = acos(dot(v2-v1, v0-v1)/(length(v2-v1)*length(v0-v1)));
            double a2 = acos(dot(v0-v2, v1-v2)/(length(v0-v2)*length(v1-v2)));
 
            if(a0>(M_PI/2.0) && a1>(M_PI/2.0) && a2>(M_PI/2.0)) // f is non-obtuse
            {
                // Add Voronoi formula (see Section 3.3)
                area_mixed += (1.0/8) * 
                    ((1.0/tan(a1)) * sqr_length(v2-v0) + 
                    (1.0/tan(a2)) * sqr_length(v1-v0));
            }
            else // Voronoi inappropriate
            {
                // Add either area(f)/4 or area(f)/2
                if(a0>M_PI/2.0)// the angle of f at x is obtuse
                    area_mixed += f_area/2;
                else
                    area_mixed += f_area/4;
            }
        }
        return area_mixed;
    }
 
    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;
    }
 
    void unnormalized_mean_curvature_normal(const Manifold& m, VertexID v, Vec3d& curv_normal, double& w_sum)
    {
        if(boundary(m, v))
            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 d_right = dot(cond_normalize(nbr-right),cond_normalize(vertex-right));
            double a_left  = acos(min(1.0, max(-1.0, d_left)));
            double a_right = acos(min(1.0, max(-1.0, d_right)));
 
            double w = 1.0/(1e-300+tan(a_left));
            w += 1.0/(1e-300+tan(a_right));
//            double w = sin(a_left + a_right) / (1e-300 + sin(a_left)*sin(a_right));
            curv_normal += w * (nbr-vertex);
            w_sum += w;
        }
 
    }
 
    Vec3d mean_curvature_normal(const Manifold& m, VertexID v)
    {
        Vec3d curv_normal;
        double w_sum;
        unnormalized_mean_curvature_normal(m, v, curv_normal, w_sum);
 
        return curv_normal / (4*voronoi_area(m, v));
    }
 
    double sum_curvatures(const Manifold& m, VertexAttributeVector<double>& curvature)
    {
        double sum = 0;
        for(VertexIDIterator v = m.vertices_begin(); v != m.vertices_end(); ++v){
            if(boundary(m, *v))
                continue;	
            sum += curvature[*v] * voronoi_area(m, *v);
        }
        return sum;
    }
 
 
    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;
        for(size_t i = 0; i < N; ++i)
        {
            double dot_prod = 
                std::max(-1.0, std::min(1.0, dot(edges[i],edges[(i+1)%N])));
            angle_sum += acos(dot_prod);
        }
        return (2*M_PI - angle_sum)/voronoi_area(m, v);
 
    }
 
    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()));
        Vec3d h_opp_norm(normal(m, w.opp().face()));
 
        Vec3d nc = cross(h_norm, h_opp_norm);
 
        double sign = (dot(nc, edge) >= 0) ? 1 : -1;
        double beta = asin(nc.length());
 
        Mat3x3d mat;
        outer_product(edge, edge, mat);
        return sign * edge_len * beta * mat;
    }
 
    Mat3x3d curvature_tensor_from_edges(const Manifold& m, VertexID v)
    {
        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;
    }
 
 
    void curvature_tensor_paraboloid(const Manifold& m, VertexID v, Mat2x2d& curv_tensor, Mat3x3d& frame)
    {
        if(boundary(m, v))
            return;
        // First estimate the normal and compute a transformation matrix
        // which takes us into tangent plane coordinates.
        Vec3d Norm = Vec3d(normal(m, v));
        Vec3d X,Y;
        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())));
 
        int N = int(points.size());
 
        CVector b(N);
        // Compute the matrix of parameter values
        CMatrix PMat(N, 3);
        for(int i = 0; i < N; ++i){
            Vec3d p = frame * (points[i]-centre);
            b[i] = p[2];
 
            PMat.set(i,0,0.5*sqr(p[0]));
            PMat.set(i,1,p[0]*p[1]);
            PMat.set(i,2,0.5*sqr(p[1]));
        }
 
        // Compute the coefficients of the polynomial surface
        CVector x(3);
        x = LinearLSSolve(PMat,b);
        if(isnan(x[0])) cout << __LINE__ << " " << PMat << b << endl ;
 
        // Finally compute the shape tensor from the coefficients
        // using the first and second fundamental forms.
        curv_tensor =Mat2x2d(x[0],x[1],x[1],x[2]);
 
    }
 
    void curvature_tensors_from_edges(const Manifold& m, VertexAttributeVector<Mat3x3d>& curvature_tensors)
    {
        for(VertexIDIterator v = m.vertices_begin(); v != m.vertices_end(); ++v)
            curvature_tensors[*v] = curvature_tensor_from_edges(m, *v);
    }
 
    void smooth_curvature_tensors(const Manifold& m, VertexAttributeVector<Mat3x3d>& curvature_tensors)
    {
        assert(curvature_tensors.size() == m.allocated_vertices());
        VertexAttributeVector<Mat3x3d> tmp_curvature_tensors(m.allocated_vertices());
        double tmp_area;
 
        for(VertexIDIterator v = m.vertices_begin(); v != m.vertices_end(); ++v){
            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;
                }
                tmp_curvature_tensors[*v] /= tmp_area;
            }
        }
        curvature_tensors = move(tmp_curvature_tensors);
    }
 
    void gaussian_curvature_angle_defects(const Manifold& m, VertexAttributeVector<double>& curvature, int smooth_steps)
    {
        for(VertexIDIterator v = m.vertices_begin(); v != m.vertices_end(); ++v)
            curvature[*v] = gaussian_curvature_angle_defect(m, *v);
 
        smooth_something_on_mesh(m, curvature, smooth_steps);
    }
 
    void mean_curvatures(const Manifold& m, VertexAttributeVector<double>& curvature, int smooth_steps)
    {
        for(VertexIDIterator v = m.vertices_begin(); v != m.vertices_end(); ++v)
			if(!boundary(m,*v))
			{
				Vec3d N = -mean_curvature_normal(m, *v);
				curvature[*v] = length(N) * sign(dot(N,Vec3d(normal(m, *v))));
			}	
        smooth_something_on_mesh(m, curvature, smooth_steps);	
    }
 
 
    void curvature_paraboloids( const Manifold& m, 
                                VertexAttributeVector<Vec3d>& min_curv_direction, 
                                VertexAttributeVector<Vec3d>& max_curv_direction,
                                VertexAttributeVector<double>& curvature)
    {
 
        for(VertexIDIterator v = m.vertices_begin(); v != m.vertices_end(); ++v){
            Mat2x2d tensor;
            Mat3x3d frame;
            curvature_tensor_paraboloid(m, *v, tensor, frame);
 
            Mat2x2d Q,L;
            int s = power_eigensolution(tensor, Q, L);
 
            if(s < 2)	
                cout << tensor << Q << L << endl;
 
            int max_idx = 0;
            int min_idx = 1;
 
            if(fabs(L[max_idx][max_idx])<fabs(L[min_idx][min_idx])) swap(max_idx, min_idx);
 
            Mat3x3d frame_t = transpose(frame);
 
            max_curv_direction[*v] = cond_normalize(frame_t * Vec3d(Q[max_idx][0], Q[max_idx][1], 0));
 
            min_curv_direction[*v] = cond_normalize(frame_t * Vec3d(Q[min_idx][0], Q[min_idx][1], 0));
 
            curvature[*v] = L[0][0]*L[1][1];
        }
    }
 
 
    void curvature_from_tensors(const Manifold& m,
                                const VertexAttributeVector<Mat3x3d>& curvature_tensors,
                                VertexAttributeVector<Vec3d>& min_curv_direction,
                                VertexAttributeVector<Vec3d>& max_curv_direction,
                                VertexAttributeVector<double>& curvature)
    {
        assert(curvature_tensors.size() == m.allocated_vertices());
 
        double max_val = -1e30;
        for(VertexIDIterator v = m.vertices_begin(); v != m.vertices_end(); ++v){
            Mat3x3d C,Q,L;
            C = curvature_tensors[*v];
            int s = power_eigensolution(C, Q, L);
            Vec3d dmin, dmax;
            if(s == 0)
            {
                Vec3d n(normal(m, *v));
                orthogonal(n, dmin, dmax);
                curvature[*v] = 0;
                cout << " rank 0 " << endl;
            }
            else if(s == 1)
            {
                Vec3d n(normal(m, *v));
                dmin = normalize(Q[0]);
                dmax = cross(n, dmin);
                curvature[*v] = 0;
                cout << " rank 1 " << endl;
            }
            else
            {
                /*				Vec3d l(fabs(L[0][0]), fabs(L[1][1]), fabs(L[2][2]));

                int z_idx=2;
                if(s==3)
                {
                if(l[0] < l[1])
                z_idx = l[0]<l[2] ? 0 : 2;
                else
                z_idx = l[1]<l[2] ? 1 : 2;
                }
                int max_idx = (z_idx + 1) % 3;
                int min_idx = (z_idx + 2) % 3;

                if(l[max_idx] < l[min_idx]) swap(max_idx, min_idx);
                */
                int max_idx = 0;
                int min_idx = 1;
                // Yes - the biggest eigenvalue corresponds to the min direction
                // and vice versa.
                dmin = normalize(Q[max_idx]);
                dmax = normalize(Q[min_idx]);
 
                curvature[*v] = L[max_idx][max_idx]*L[min_idx][min_idx];
 
            }
            min_curv_direction[*v] = dmin;
            max_curv_direction[*v] = dmax;
            max_val = max(fabs(curvature[*v]), max_val);
 
        }
        //scal = 1.0/max_val;
    }
}
 
 

Rev 631	Rev 632
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 "curvature.h"`	7	`#include "curvature.h"`
8		8
9	`#include "../GLGraphics/gel_glut.h"`	9	`#include "../GLGraphics/gel_glut.h"`
10		10
11	`#include <iostream>`	11	`#include <iostream>`
12	`#include "../CGLA/eigensolution.h"`	12	`#include "../CGLA/eigensolution.h"`
13	`#include "../CGLA/Vec2d.h"`	13	`#include "../CGLA/Vec2d.h"`
14	`#include "../CGLA/Vec3d.h"`	14	`#include "../CGLA/Vec3d.h"`
15	`#include "../CGLA/Mat3x3d.h"`	15	`#include "../CGLA/Mat3x3d.h"`
16	`#include "../CGLA/Mat2x2d.h"`	16	`#include "../CGLA/Mat2x2d.h"`
17	`#include "../CGLA/Mat2x3d.h"`	17	`#include "../CGLA/Mat2x3d.h"`
18		18
19	`#include "Manifold.h"`	19	`#include "Manifold.h"`
20	`#include "AttributeVector.h"`	20	`#include "AttributeVector.h"`
21	`#include "x3d_save.h"`	21	`#include "x3d_save.h"`
22	`#include "x3d_load.h"`	22	`#include "x3d_load.h"`
23	`#include "obj_load.h"`	23	`#include "obj_load.h"`
24	`#include "mesh_optimization.h"`	24	`#include "mesh_optimization.h"`
25		25
26	`#include "../LinAlg/Matrix.h"`	26	`#include "../LinAlg/Matrix.h"`
27	`#include "../LinAlg/Vector.h"`	27	`#include "../LinAlg/Vector.h"`
28	`#include "../LinAlg/LapackFunc.h"`	28	`#include "../LinAlg/LapackFunc.h"`
29		29
30	`using namespace std;`	30	`using namespace std;`
31		31
32	`using namespace LinAlg;`	32	`using namespace LinAlg;`
33	`using namespace CGLA;`	33	`using namespace CGLA;`
34	`using namespace HMesh;`	34	`using namespace HMesh;`
35		35
36	`namespace HMesh`	36	`namespace HMesh`
37	`{`	37	`{`
38	`namespace`	38	`namespace`
39	`{`	39	`{`
40	`//double scal = 0.001;`	40	`//double scal = 0.001;`
41	`//double vector_scal = 0.001;`	41	`//double vector_scal = 0.001;`
42		42
43	`template<class T>`	43	`template<class T>`
44	`void smooth_something_on_mesh(const Manifold& m, VertexAttributeVector<T>& vec, int smooth_steps)`	44	`void smooth_something_on_mesh(const Manifold& m, VertexAttributeVector<T>& vec, int smooth_steps)`
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(Walker w = m.walker(*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);`
56	`}`	56	`}`
57	`}`	57	`}`
58	`}`	58	`}`
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(Walker w = m.walker(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()));`
70		70
71	`double a0 = acos(dot(v1-v0, v2-v0)/(length(v1-v0)*length(v2-v0)));`	71	`double a0 = acos(dot(v1-v0, v2-v0)/(length(v1-v0)*length(v2-v0)));`
72	`double a1 = acos(dot(v2-v1, v0-v1)/(length(v2-v1)*length(v0-v1)));`	72	`double a1 = acos(dot(v2-v1, v0-v1)/(length(v2-v1)*length(v0-v1)));`
73	`double a2 = acos(dot(v0-v2, v1-v2)/(length(v0-v2)*length(v1-v2)));`	73	`double a2 = acos(dot(v0-v2, v1-v2)/(length(v0-v2)*length(v1-v2)));`
74		74
75	`if(a0>(M_PI/2.0) && a1>(M_PI/2.0) && a2>(M_PI/2.0)) // f is non-obtuse`	75	`if(a0>(M_PI/2.0) && a1>(M_PI/2.0) && a2>(M_PI/2.0)) // f is non-obtuse`
76	`{`	76	`{`
77	`// Add Voronoi formula (see Section 3.3)`	77	`// Add Voronoi formula (see Section 3.3)`
78	`area_mixed += (1.0/8) *`	78	`area_mixed += (1.0/8) *`
79	`((1.0/tan(a1)) * sqr_length(v2-v0) +`	79	`((1.0/tan(a1)) * sqr_length(v2-v0) +`
80	`(1.0/tan(a2)) * sqr_length(v1-v0));`	80	`(1.0/tan(a2)) * sqr_length(v1-v0));`
81	`}`	81	`}`
82	`else // Voronoi inappropriate`	82	`else // Voronoi inappropriate`
83	`{`	83	`{`
84	`// Add either area(f)/4 or area(f)/2`	84	`// Add either area(f)/4 or area(f)/2`
85	`if(a0>M_PI/2.0)// the angle of f at x is obtuse`	85	`if(a0>M_PI/2.0)// the angle of f at x is obtuse`
86	`area_mixed += f_area/2;`	86	`area_mixed += f_area/2;`
87	`else`	87	`else`
88	`area_mixed += f_area/4;`	88	`area_mixed += f_area/4;`
89	`}`	89	`}`
90	`}`	90	`}`
91	`return area_mixed;`	91	`return area_mixed;`
92	`}`	92	`}`
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(Walker w = m.walker(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
104	`void unnormalized_mean_curvature_normal(const Manifold& m, VertexID v, Vec3d& curv_normal, double& w_sum)`	104	`void unnormalized_mean_curvature_normal(const Manifold& m, VertexID v, Vec3d& curv_normal, double& w_sum)`
105	`{`	105	`{`
106	`if(boundary(m, v))`	106	`if(boundary(m, v))`
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(Walker walker = m.walker(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));`
118	`double d_right = dot(cond_normalize(nbr-right),cond_normalize(vertex-right));`	118	`double d_right = dot(cond_normalize(nbr-right),cond_normalize(vertex-right));`
119	`double a_left = acos(min(1.0, max(-1.0, d_left)));`	119	`double a_left = acos(min(1.0, max(-1.0, d_left)));`
120	`double a_right = acos(min(1.0, max(-1.0, d_right)));`	120	`double a_right = acos(min(1.0, max(-1.0, d_right)));`
121		121
122	`double w = 1.0/(1e-300+tan(a_left));`	122	`double w = 1.0/(1e-300+tan(a_left));`
123	`w += 1.0/(1e-300+tan(a_right));`	123	`w += 1.0/(1e-300+tan(a_right));`
124	`// double w = sin(a_left + a_right) / (1e-300 + sin(a_left)*sin(a_right));`	124	`// double w = sin(a_left + a_right) / (1e-300 + sin(a_left)*sin(a_right));`
125	`curv_normal += w * (nbr-vertex);`	125	`curv_normal += w * (nbr-vertex);`
126	`w_sum += w;`	126	`w_sum += w;`
127	`}`	127	`}`
128		128
129	`}`	129	`}`
130		130
131	`Vec3d mean_curvature_normal(const Manifold& m, VertexID v)`	131	`Vec3d mean_curvature_normal(const Manifold& m, VertexID v)`
132	`{`	132	`{`
133	`Vec3d curv_normal;`	133	`Vec3d curv_normal;`
134	`double w_sum;`	134	`double w_sum;`
135	`unnormalized_mean_curvature_normal(m, v, curv_normal, w_sum);`	135	`unnormalized_mean_curvature_normal(m, v, curv_normal, w_sum);`
136		136
137	`return curv_normal / (4*voronoi_area(m, v));`	137	`return curv_normal / (4*voronoi_area(m, v));`
138	`}`	138	`}`
139		139
140	`double sum_curvatures(const Manifold& m, VertexAttributeVector<double>& curvature)`	140	`double sum_curvatures(const Manifold& m, VertexAttributeVector<double>& curvature)`
141	`{`	141	`{`
142	`double sum = 0;`	142	`double sum = 0;`
143	`for(VertexIDIterator v = m.vertices_begin(); v != m.vertices_end(); ++v){`	143	`for(VertexIDIterator v = m.vertices_begin(); v != m.vertices_end(); ++v){`
144	`if(boundary(m, *v))`	144	`if(boundary(m, *v))`
145	`continue;`	145	`continue;`
146	`sum += curvature[v] voronoi_area(m, *v);`	146	`sum += curvature[v] voronoi_area(m, *v);`
147	`}`	147	`}`
148	`return sum;`	148	`return sum;`
149	`}`	149	`}`
150		150
151		151
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	`Vec3d vertex(m.pos(v));`	157	`Vec3d vertex(m.pos(v));`
158	`vector<Vec3d> edges;`	158	`vector<Vec3d> edges;`
159	`for(Walker w = m.walker(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;`
165	`for(size_t i = 0; i < N; ++i)`	165	`for(size_t i = 0; i < N; ++i)`
166	`{`	166	`{`
167	`double dot_prod =`	167	`double dot_prod =`
168	`s_max(-1.0, s_min(1.0, dot(edges[i],edges[(i+1)%N])));`	168	`std::max(-1.0, std::min(1.0, dot(edges[i],edges[(i+1)%N])));`
169	`angle_sum += acos(dot_prod);`	169	`angle_sum += acos(dot_prod);`
170	`}`	170	`}`
171	`return (2*M_PI - angle_sum)/voronoi_area(m, v);`	171	`return (2*M_PI - angle_sum)/voronoi_area(m, v);`
172		172
173	`}`	173	`}`
174		174
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	`Walker w = m.walker(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()));`
186	`Vec3d h_opp_norm(normal(m, w.opp().face()));`	186	`Vec3d h_opp_norm(normal(m, w.opp().face()));`
187		187
188	`Vec3d nc = cross(h_norm, h_opp_norm);`	188	`Vec3d nc = cross(h_norm, h_opp_norm);`
189		189
190	`double sign = (dot(nc, edge) >= 0) ? 1 : -1;`	190	`double sign = (dot(nc, edge) >= 0) ? 1 : -1;`
191	`double beta = asin(nc.length());`	191	`double beta = asin(nc.length());`
192		192
193	`Mat3x3d mat;`	193	`Mat3x3d mat;`
194	`outer_product(edge, edge, mat);`	194	`outer_product(edge, edge, mat);`
195	`return sign * edge_len * beta * mat;`	195	`return sign * edge_len * beta * mat;`
196	`}`	196	`}`
197		197
198	`Mat3x3d curvature_tensor_from_edges(const Manifold& m, VertexID v)`	198	`Mat3x3d curvature_tensor_from_edges(const Manifold& m, VertexID v)`
199	`{`	199	`{`
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(Walker w = m.walker(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;`
211	`}`	211	`}`
212		212
213		213
214	`void curvature_tensor_paraboloid(const Manifold& m, VertexID v, Mat2x2d& curv_tensor, Mat3x3d& frame)`	214	`void curvature_tensor_paraboloid(const Manifold& m, VertexID v, Mat2x2d& curv_tensor, Mat3x3d& frame)`
215	`{`	215	`{`
216	`if(boundary(m, v))`	216	`if(boundary(m, v))`
217	`return;`	217	`return;`
218	`// First estimate the normal and compute a transformation matrix`	218	`// First estimate the normal and compute a transformation matrix`
219	`// which takes us into tangent plane coordinates.`	219	`// which takes us into tangent plane coordinates.`
220	`Vec3d Norm = Vec3d(normal(m, v));`	220	`Vec3d Norm = Vec3d(normal(m, v));`
221	`Vec3d X,Y;`	221	`Vec3d X,Y;`
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(Walker w = m.walker(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	`int N = int(points.size());`	230	`int N = int(points.size());`
231		231
232	`CVector b(N);`	232	`CVector b(N);`
233	`// Compute the matrix of parameter values`	233	`// Compute the matrix of parameter values`
234	`CMatrix PMat(N, 3);`	234	`CMatrix PMat(N, 3);`
235	`for(int i = 0; i < N; ++i){`	235	`for(int i = 0; i < N; ++i){`
236	`Vec3d p = frame * (points[i]-centre);`	236	`Vec3d p = frame * (points[i]-centre);`
237	`b[i] = p[2];`	237	`b[i] = p[2];`
238		238
239	`PMat.set(i,0,0.5*sqr(p[0]));`	239	`PMat.set(i,0,0.5*sqr(p[0]));`
240	`PMat.set(i,1,p[0]*p[1]);`	240	`PMat.set(i,1,p[0]*p[1]);`
241	`PMat.set(i,2,0.5*sqr(p[1]));`	241	`PMat.set(i,2,0.5*sqr(p[1]));`
242	`}`	242	`}`
243		243
244	`// Compute the coefficients of the polynomial surface`	244	`// Compute the coefficients of the polynomial surface`
245	`CVector x(3);`	245	`CVector x(3);`
246	`x = LinearLSSolve(PMat,b);`	246	`x = LinearLSSolve(PMat,b);`
247	`if(isnan(x[0])) cout << __LINE__ << " " << PMat << b << endl ;`	247	`if(isnan(x[0])) cout << __LINE__ << " " << PMat << b << endl ;`
248		248
249	`// Finally compute the shape tensor from the coefficients`	249	`// Finally compute the shape tensor from the coefficients`
250	`// using the first and second fundamental forms.`	250	`// using the first and second fundamental forms.`
251	`curv_tensor =Mat2x2d(x[0],x[1],x[1],x[2]);`	251	`curv_tensor =Mat2x2d(x[0],x[1],x[1],x[2]);`
252		252
253	`}`	253	`}`
254		254
255	`void curvature_tensors_from_edges(const Manifold& m, VertexAttributeVector<Mat3x3d>& curvature_tensors)`	255	`void curvature_tensors_from_edges(const Manifold& m, VertexAttributeVector<Mat3x3d>& curvature_tensors)`
256	`{`	256	`{`
257	`for(VertexIDIterator v = m.vertices_begin(); v != m.vertices_end(); ++v)`	257	`for(VertexIDIterator v = m.vertices_begin(); v != m.vertices_end(); ++v)`
258	`curvature_tensors[v] = curvature_tensor_from_edges(m, v);`	258	`curvature_tensors[v] = curvature_tensor_from_edges(m, v);`
259	`}`	259	`}`
260		260
261	`void smooth_curvature_tensors(const Manifold& m, VertexAttributeVector<Mat3x3d>& curvature_tensors)`	261	`void smooth_curvature_tensors(const Manifold& m, VertexAttributeVector<Mat3x3d>& curvature_tensors)`
262	`{`	262	`{`
263	`assert(curvature_tensors.size() == m.allocated_vertices());`	263	`assert(curvature_tensors.size() == m.allocated_vertices());`
264	`VertexAttributeVector<Mat3x3d> tmp_curvature_tensors(m.allocated_vertices());`	264	`VertexAttributeVector<Mat3x3d> tmp_curvature_tensors(m.allocated_vertices());`
265	`double tmp_area;`	265	`double tmp_area;`
266		266
267	`for(VertexIDIterator v = m.vertices_begin(); v != m.vertices_end(); ++v){`	267	`for(VertexIDIterator v = m.vertices_begin(); v != m.vertices_end(); ++v){`
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(Walker w = m.walker(*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	`}`
279	`tmp_curvature_tensors[*v] /= tmp_area;`	279	`tmp_curvature_tensors[*v] /= tmp_area;`
280	`}`	280	`}`
281	`}`	281	`}`
282	`curvature_tensors = move(tmp_curvature_tensors);`	282	`curvature_tensors = move(tmp_curvature_tensors);`
283	`}`	283	`}`
284		284
285	`void gaussian_curvature_angle_defects(const Manifold& m, VertexAttributeVector<double>& curvature, int smooth_steps)`	285	`void gaussian_curvature_angle_defects(const Manifold& m, VertexAttributeVector<double>& curvature, int smooth_steps)`
286	`{`	286	`{`
287	`for(VertexIDIterator v = m.vertices_begin(); v != m.vertices_end(); ++v)`	287	`for(VertexIDIterator v = m.vertices_begin(); v != m.vertices_end(); ++v)`
288	`curvature[v] = gaussian_curvature_angle_defect(m, v);`	288	`curvature[v] = gaussian_curvature_angle_defect(m, v);`
289		289
290	`smooth_something_on_mesh(m, curvature, smooth_steps);`	290	`smooth_something_on_mesh(m, curvature, smooth_steps);`
291	`}`	291	`}`
292		292
293	`void mean_curvatures(const Manifold& m, VertexAttributeVector<double>& curvature, int smooth_steps)`	293	`void mean_curvatures(const Manifold& m, VertexAttributeVector<double>& curvature, int smooth_steps)`
294	`{`	294	`{`
295	`for(VertexIDIterator v = m.vertices_begin(); v != m.vertices_end(); ++v)`	295	`for(VertexIDIterator v = m.vertices_begin(); v != m.vertices_end(); ++v)`
296	`if(!boundary(m,*v))`	296	`if(!boundary(m,*v))`
297	`{`	297	`{`
298	`Vec3d N = -mean_curvature_normal(m, *v);`	298	`Vec3d N = -mean_curvature_normal(m, *v);`
299	`curvature[v] = length(N) sign(dot(N,Vec3d(normal(m, *v))));`	299	`curvature[v] = length(N) sign(dot(N,Vec3d(normal(m, *v))));`
300	`}`	300	`}`
301	`smooth_something_on_mesh(m, curvature, smooth_steps);`	301	`smooth_something_on_mesh(m, curvature, smooth_steps);`
302	`}`	302	`}`
303		303
304		304
305	`void curvature_paraboloids( const Manifold& m,`	305	`void curvature_paraboloids( const Manifold& m,`
306	`VertexAttributeVector<Vec3d>& min_curv_direction,`	306	`VertexAttributeVector<Vec3d>& min_curv_direction,`
307	`VertexAttributeVector<Vec3d>& max_curv_direction,`	307	`VertexAttributeVector<Vec3d>& max_curv_direction,`
308	`VertexAttributeVector<double>& curvature)`	308	`VertexAttributeVector<double>& curvature)`
309	`{`	309	`{`
310		310
311	`for(VertexIDIterator v = m.vertices_begin(); v != m.vertices_end(); ++v){`	311	`for(VertexIDIterator v = m.vertices_begin(); v != m.vertices_end(); ++v){`
312	`Mat2x2d tensor;`	312	`Mat2x2d tensor;`
313	`Mat3x3d frame;`	313	`Mat3x3d frame;`
314	`curvature_tensor_paraboloid(m, *v, tensor, frame);`	314	`curvature_tensor_paraboloid(m, *v, tensor, frame);`
315		315
316	`Mat2x2d Q,L;`	316	`Mat2x2d Q,L;`
317	`int s = power_eigensolution(tensor, Q, L);`	317	`int s = power_eigensolution(tensor, Q, L);`
318		318
319	`if(s < 2)`	319	`if(s < 2)`
320	`cout << tensor << Q << L << endl;`	320	`cout << tensor << Q << L << endl;`
321		321
322	`int max_idx = 0;`	322	`int max_idx = 0;`
323	`int min_idx = 1;`	323	`int min_idx = 1;`
324		324
325	`if(fabs(L[max_idx][max_idx])<fabs(L[min_idx][min_idx])) swap(max_idx, min_idx);`	325	`if(fabs(L[max_idx][max_idx])<fabs(L[min_idx][min_idx])) swap(max_idx, min_idx);`
326		326
327	`Mat3x3d frame_t = transpose(frame);`	327	`Mat3x3d frame_t = transpose(frame);`
328		328
329	`max_curv_direction[v] = cond_normalize(frame_t Vec3d(Q[max_idx][0], Q[max_idx][1], 0));`	329	`max_curv_direction[v] = cond_normalize(frame_t Vec3d(Q[max_idx][0], Q[max_idx][1], 0));`
330		330
331	`min_curv_direction[v] = cond_normalize(frame_t Vec3d(Q[min_idx][0], Q[min_idx][1], 0));`	331	`min_curv_direction[v] = cond_normalize(frame_t Vec3d(Q[min_idx][0], Q[min_idx][1], 0));`
332		332
333	`curvature[v] = L[0][0]L[1][1];`	333	`curvature[v] = L[0][0]L[1][1];`
334	`}`	334	`}`
335	`}`	335	`}`
336		336
337		337
338	`void curvature_from_tensors(const Manifold& m,`	338	`void curvature_from_tensors(const Manifold& m,`
339	`const VertexAttributeVector<Mat3x3d>& curvature_tensors,`	339	`const VertexAttributeVector<Mat3x3d>& curvature_tensors,`
340	`VertexAttributeVector<Vec3d>& min_curv_direction,`	340	`VertexAttributeVector<Vec3d>& min_curv_direction,`
341	`VertexAttributeVector<Vec3d>& max_curv_direction,`	341	`VertexAttributeVector<Vec3d>& max_curv_direction,`
342	`VertexAttributeVector<double>& curvature)`	342	`VertexAttributeVector<double>& curvature)`
343	`{`	343	`{`
344	`assert(curvature_tensors.size() == m.allocated_vertices());`	344	`assert(curvature_tensors.size() == m.allocated_vertices());`
345		345
346	`double max_val = -1e30;`	346	`double max_val = -1e30;`
347	`for(VertexIDIterator v = m.vertices_begin(); v != m.vertices_end(); ++v){`	347	`for(VertexIDIterator v = m.vertices_begin(); v != m.vertices_end(); ++v){`
348	`Mat3x3d C,Q,L;`	348	`Mat3x3d C,Q,L;`
349	`C = curvature_tensors[*v];`	349	`C = curvature_tensors[*v];`
350	`int s = power_eigensolution(C, Q, L);`	350	`int s = power_eigensolution(C, Q, L);`
351	`Vec3d dmin, dmax;`	351	`Vec3d dmin, dmax;`
352	`if(s == 0)`	352	`if(s == 0)`
353	`{`	353	`{`
354	`Vec3d n(normal(m, *v));`	354	`Vec3d n(normal(m, *v));`
355	`orthogonal(n, dmin, dmax);`	355	`orthogonal(n, dmin, dmax);`
356	`curvature[*v] = 0;`	356	`curvature[*v] = 0;`
357	`cout << " rank 0 " << endl;`	357	`cout << " rank 0 " << endl;`
358	`}`	358	`}`
359	`else if(s == 1)`	359	`else if(s == 1)`
360	`{`	360	`{`
361	`Vec3d n(normal(m, *v));`	361	`Vec3d n(normal(m, *v));`
362	`dmin = normalize(Q[0]);`	362	`dmin = normalize(Q[0]);`
363	`dmax = cross(n, dmin);`	363	`dmax = cross(n, dmin);`
364	`curvature[*v] = 0;`	364	`curvature[*v] = 0;`
365	`cout << " rank 1 " << endl;`	365	`cout << " rank 1 " << endl;`
366	`}`	366	`}`
367	`else`	367	`else`
368	`{`	368	`{`
369	`/* Vec3d l(fabs(L[0][0]), fabs(L[1][1]), fabs(L[2][2]));`	369	`/* Vec3d l(fabs(L[0][0]), fabs(L[1][1]), fabs(L[2][2]));`
370		370
371	`int z_idx=2;`	371	`int z_idx=2;`
372	`if(s==3)`	372	`if(s==3)`
373	`{`	373	`{`
374	`if(l[0] < l[1])`	374	`if(l[0] < l[1])`
375	`z_idx = l[0]<l[2] ? 0 : 2;`	375	`z_idx = l[0]<l[2] ? 0 : 2;`
376	`else`	376	`else`
377	`z_idx = l[1]<l[2] ? 1 : 2;`	377	`z_idx = l[1]<l[2] ? 1 : 2;`
378	`}`	378	`}`
379	`int max_idx = (z_idx + 1) % 3;`	379	`int max_idx = (z_idx + 1) % 3;`
380	`int min_idx = (z_idx + 2) % 3;`	380	`int min_idx = (z_idx + 2) % 3;`
381		381
382	`if(l[max_idx] < l[min_idx]) swap(max_idx, min_idx);`	382	`if(l[max_idx] < l[min_idx]) swap(max_idx, min_idx);`
383	`*/`	383	`*/`
384	`int max_idx = 0;`	384	`int max_idx = 0;`
385	`int min_idx = 1;`	385	`int min_idx = 1;`
386	`// Yes - the biggest eigenvalue corresponds to the min direction`	386	`// Yes - the biggest eigenvalue corresponds to the min direction`
387	`// and vice versa.`	387	`// and vice versa.`
388	`dmin = normalize(Q[max_idx]);`	388	`dmin = normalize(Q[max_idx]);`
389	`dmax = normalize(Q[min_idx]);`	389	`dmax = normalize(Q[min_idx]);`
390		390
391	`curvature[v] = L[max_idx][max_idx]L[min_idx][min_idx];`	391	`curvature[v] = L[max_idx][max_idx]L[min_idx][min_idx];`
392		392
393	`}`	393	`}`
394	`min_curv_direction[*v] = dmin;`	394	`min_curv_direction[*v] = dmin;`
395	`max_curv_direction[*v] = dmax;`	395	`max_curv_direction[*v] = dmax;`
396	`max_val = max(fabs(curvature[*v]), max_val);`	396	`max_val = max(fabs(curvature[*v]), max_val);`
397		397
398	`}`	398	`}`
399	`//scal = 1.0/max_val;`	399	`//scal = 1.0/max_val;`
400	`}`	400	`}`
401	`}`	401	`}`
402		402
403		403

Subversion Repositories gelsvn

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