WebSVN – gelsvn – Diff – /trunk/apps/MeshEdit/curvature.cpp


/*
 *  curvature.cpp
 *  GEL
 *
 *  Created by J. Andreas Bærentzen on 23/09/08.
 *  Copyright 2008 __MyCompanyName__. All rights reserved.
 *
 */
 
#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 <HMesh/VertexCirculator.h>
#include <HMesh/FaceCirculator.h>
#include <HMesh/x3d_save.h>
#include <HMesh/x3d_load.h>
#include <HMesh/obj_load.h>
#include <HMesh/build_manifold.h>
#include <HMesh/mesh_optimization.h>
 
#include <LinAlg/Matrix.h>
#include <LinAlg/Vector.h>
#include <LinAlg/LapackFunc.h>
 
using namespace std;
using namespace HMesh;
using namespace LinAlg;
using namespace CGLA;
 
namespace {
	double scal = 0.001;
	double vector_scal = 0.001;
	
	template<class T> 
	void smooth_something_on_mesh(Manifold& m, vector<T>& vec, int smooth_steps)
	{
		for(int iter=0;iter<smooth_steps;++iter)
		{
			vector<T> new_vec(m.no_vertices());
			for(VertexIter vi=m.vertices_begin(); vi != m.vertices_end(); ++vi)
			{
				int i = vi->touched;
				new_vec[i] = vec[i];
				for(VertexCirculator vc(vi); !vc.end();++vc)
				{
					int j = vc.get_vertex()->touched;
					new_vec[i] += vec[j];
				}
				new_vec[i] /=(valency(vi)+ 1.0);
			}
			swap(vec,new_vec);
		}		
	}
}
 
double voronoi_area(VertexIter v)
{
	double area_mixed = 0;
	//For each triangle T from the 1-ring neighborhood of x
	for(VertexCirculator vc(v); !vc.end(); ++vc)
	{
		FaceIter f = vc.get_face();
		double f_area = area(f);
		
		HalfEdgeIter he = vc.get_halfedge();
		Vec3d v1(he->vert->pos);
		Vec3d v2(he->next->vert->pos);
		Vec3d v0(he->next->next->vert->pos);
		
		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(VertexIter v)
{
	double barea = 0;
	//For each triangle T from the 1-ring neighborhood of x
	for(VertexCirculator vc(v); !vc.end(); ++vc)
	{
		FaceIter f = vc.get_face();
		barea += area(f)/3.0;
	}
	return barea;
}
 
void unnormalized_mean_curvature_normal(VertexIter v, Vec3d& curv_normal, double& w_sum)
{
	if(!is_boundary(v))
	{
		Vec3d vertex(v->pos);
		curv_normal = Vec3d(0);
		w_sum = 0;
		for(VertexCirculator vc(v); !vc.end(); ++vc)
		{
			HalfEdgeIter h = vc.get_halfedge();
			Vec3d nbr(h->vert->pos);
			Vec3d left(h->next->vert->pos);
			Vec3d right(h->opp->prev->opp->vert->pos);
			
			double d_left = dot(normalize(nbr-left),
								normalize(vertex-left));
			double d_right = dot(normalize(nbr-right),
								 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/tan(a_left) + 1.0/tan(a_right);
			
			curv_normal += w * (vertex-nbr);
 
			w_sum += w;
		}
 
 
	}
}
 
const Vec3d mean_curvature_normal(VertexIter v)
{
	Vec3d curv_normal;
	double w_sum;
	unnormalized_mean_curvature_normal(v, curv_normal, w_sum);
	
	return curv_normal / (4*voronoi_area(v));
}
 
double sum_curvatures(Manifold& m, vector<double>& curvature)
{
	double sum = 0;
	for(VertexIter v=m.vertices_begin(); v!=m.vertices_end(); ++v)
	{
		if(!is_boundary(v))
		{
			sum += curvature[v->touched] * voronoi_area(v);
		}
	}
	return sum;
}
 
 
const double gaussian_curvature_angle_defect(VertexIter v)
{
	if(!is_boundary(v))
	{
		Vec3f vertex(v->pos);
		vector<Vec3d> edges;
		for(VertexCirculator vc(v); !vc.end(); ++vc)
		{
			Vec3d e(normalize(vc.get_vertex()->pos-vertex));
			edges.push_back(e);
		}
		int N=edges.size();
		double angle_sum = 0;
		for(int i=0;i<N;++i)
		{
			double dot_prod = 
			s_max(-1.0, s_min(1.0, dot(edges[i],edges[(i+1)%N])));
			angle_sum += acos(dot_prod);
		}
		return (2*M_PI - angle_sum)/voronoi_area(v);
	}
	return 0;
}
 
const Mat3x3d curvature_tensor(HalfEdgeIter h)
{
	if(!is_boundary(h))
	{
		Vec3d edge(h->vert->pos - h->opp->vert->pos);
		double edge_len = length(edge);
		edge /= edge_len;
		
		Vec3d h_norm(normal(h->face));
		Vec3d h_opp_norm(normal(h->opp->face));
		
		Vec3d nc = cross(h_norm, h_opp_norm);
		
		double sign = (dot(nc, edge) >= 0) ? 1 : -1;
		double beta = asin(nc.length());
		
		Mat3x3d m;
		outer_product(edge, edge, m);
		return sign * edge_len * beta * m;
	}
	return Mat3x3d(0);
}
 
const Mat3x3d curvature_tensor_from_edges(VertexIter v)
{
	Mat3x3d curv_tensor(0);
	
	if(!is_boundary(v))
	{
		for(VertexCirculator vc(v); !vc.end(); ++vc)
		{
			curv_tensor += 0.5*curvature_tensor(vc.get_halfedge());
		}
		curv_tensor /= voronoi_area(v);
	}
	return curv_tensor;
}
 
 
void curvature_tensor_paraboloid(VertexIter v,
								 Mat2x2d& curv_tensor,
								 Mat3x3d& frame)
{
	if(!is_boundary(v))
	{
		// First estimate the normal and compute a transformation matrix
		// which takes us into tangent plane coordinates.
		Vec3d Norm = Vec3d(normal(v));
		Vec3d X,Y;
		orthogonal(Norm,X,Y);
		frame = Mat3x3d(X,Y,Norm);
		Vec3d centre(v->pos);
		
		vector<Vec3d> points;
		for(VertexCirculator vc(v); !vc.end(); ++vc)
			points.push_back(Vec3d(vc.get_vertex()->pos));
		
		int N = 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(Manifold& m, 
								  vector<Mat3x3d>& curvature_tensors)
{
	curvature_tensors.resize(m.no_vertices());
	int i=0;
	for(VertexIter vi=m.vertices_begin(); vi != m.vertices_end(); ++vi,++i)
	{
		vi->touched = i;
		curvature_tensors[i] = curvature_tensor_from_edges(vi);
	}
}
 
void smooth_curvature_tensors(Manifold& m,																	
							  vector<Mat3x3d>& curvature_tensors)
{
	assert(curvature_tensors.size() == m.no_vertices());
	vector<Mat3x3d> tmp_curvature_tensors(m.no_vertices());
	double tmp_area;
	int i=0;
	for(VertexIter vi=m.vertices_begin(); vi != m.vertices_end(); ++vi,++i)
		if(!is_boundary(vi))
		{
			double a = voronoi_area(vi);
			tmp_curvature_tensors[i] = curvature_tensors[i] * a;
			tmp_area = a;
			int n=1;
			for(VertexCirculator vc(vi); !vc.end(); ++vc)
				if(!is_boundary(vc.get_vertex()))
				{
					int j = vc.get_vertex()->touched;
					double a = voronoi_area(vc.get_vertex());
					tmp_curvature_tensors[i] += curvature_tensors[j]*a;
					tmp_area += a;
					++n;
				}
			tmp_curvature_tensors[i] /= tmp_area;
		}
	curvature_tensors = tmp_curvature_tensors;
}
 
void gaussian_curvature_angle_defects(Manifold& m,
									  vector<double>& curvature,
									  int smooth_steps)
{
	m.enumerate_vertices();
	curvature.resize(m.no_vertices());
	for(VertexIter vi=m.vertices_begin(); vi != m.vertices_end(); ++vi)
	{
		int i = vi->touched;
		curvature[i] = gaussian_curvature_angle_defect(vi);
	}
	smooth_something_on_mesh(m, curvature, smooth_steps);
}
 
void mean_curvatures(Manifold& m, vector<double>& curvature,int smooth_steps)
{
	curvature.resize(m.no_vertices());
	int i=0;
	for(VertexIter vi=m.vertices_begin(); vi != m.vertices_end(); ++vi,++i)
	{
		vi->touched = i;
		Vec3d N = mean_curvature_normal(vi);
		curvature[i] = length(N)*sign(dot(N,Vec3d(normal(vi))));
	}	
	smooth_something_on_mesh(m, curvature, smooth_steps);	
}
 
 
void curvature_paraboloids(Manifold& m,
						   vector<Vec3d>& min_curv_direction,
						   vector<Vec3d>& max_curv_direction,
						   vector<double>& curvature)
{
	int i=0;
	for(VertexIter vi=m.vertices_begin(); vi != m.vertices_end(); ++vi,++i)
	{
		vi->touched = i;
		Mat2x2d tensor;
		Mat3x3d frame;
		curvature_tensor_paraboloid(vi, 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[i] = 
		frame_t * Vec3d(Q[max_idx][0], Q[max_idx][1], 0);
		
		min_curv_direction[i] = 
		frame_t * Vec3d(Q[min_idx][0], Q[min_idx][1], 0);
		
		curvature[i] = L[0][0]*L[1][1];
	}
}
 
 
void curvature_from_tensors(Manifold& m,
							const vector<Mat3x3d>& curvature_tensors,
							vector<Vec3d>& min_curv_direction,
							vector<Vec3d>& max_curv_direction,
							vector<double>& curvature)
{
	assert(curvature_tensors.size() == m.no_vertices());
	min_curv_direction.resize(m.no_vertices());
	max_curv_direction.resize(m.no_vertices());
	curvature.resize(m.no_vertices());
	double max_val = -1e30;
	for(VertexIter vi=m.vertices_begin(); vi != m.vertices_end(); ++vi)
	{
		int i = vi->touched;
		Mat3x3d C,Q,L;
		C = curvature_tensors[i];
		int s = power_eigensolution(C, Q, L);
		Vec3d dmin, dmax;
		if(s==0)
		{
			Vec3d n(normal(vi));
			orthogonal(n,dmin, dmax);
			curvature[i] = 0;
			cout << " rank 0 " << endl;
		}
		else if(s==1)
		{
			Vec3d n(normal(vi));
			dmin = normalize(Q[0]);
			dmax = cross(n, dmin);
			curvature[i] = 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[i] = L[max_idx][max_idx]*L[min_idx][min_idx];
			
		}
		min_curv_direction[i] = dmin;
		max_curv_direction[i] = dmax;
		max_val = max(fabs(curvature[i]), max_val);
		
	}
	scal = 1.0/max_val;
}
 
 

Rev 403	Rev 412
1	`/*`	1	`/*`
2	`* curvature.cpp`	2	`* curvature.cpp`
3	`* GEL`	3	`* GEL`
4	`*`	4	`*`
5	`* Created by J. Andreas Bærentzen on 23/09/08.`	5	`* Created by J. Andreas Bærentzen on 23/09/08.`
6	`* Copyright 2008 __MyCompanyName__. All rights reserved.`	6	`* Copyright 2008 __MyCompanyName__. All rights reserved.`
7	`*`	7	`*`
8	`*/`	8	`*/`
9		9
10	`#include "curvature.h"`	10	`#include "curvature.h"`
11		11
12	`#include <GLGraphics/gel_glut.h>`	12	`#include <GLGraphics/gel_glut.h>`
13		13
14	`#include <iostream>`	14	`#include <iostream>`
15	`#include <CGLA/eigensolution.h>`	15	`#include <CGLA/eigensolution.h>`
16	`#include <CGLA/Vec2d.h>`	16	`#include <CGLA/Vec2d.h>`
17	`#include <CGLA/Vec3d.h>`	17	`#include <CGLA/Vec3d.h>`
18	`#include <CGLA/Mat3x3d.h>`	18	`#include <CGLA/Mat3x3d.h>`
19	`#include <CGLA/Mat2x2d.h>`	19	`#include <CGLA/Mat2x2d.h>`
20	`#include <CGLA/Mat2x3d.h>`	20	`#include <CGLA/Mat2x3d.h>`
21		21
22	`#include <HMesh/VertexCirculator.h>`	22	`#include <HMesh/VertexCirculator.h>`
23	`#include <HMesh/FaceCirculator.h>`	23	`#include <HMesh/FaceCirculator.h>`
24	`#include <HMesh/x3d_save.h>`	24	`#include <HMesh/x3d_save.h>`
25	`#include <HMesh/x3d_load.h>`	25	`#include <HMesh/x3d_load.h>`
26	`#include <HMesh/obj_load.h>`	26	`#include <HMesh/obj_load.h>`
27	`#include <HMesh/build_manifold.h>`	27	`#include <HMesh/build_manifold.h>`
28	`#include <HMesh/mesh_optimization.h>`	28	`#include <HMesh/mesh_optimization.h>`
29		29
30	`#include <LinAlg/Matrix.h>`	30	`#include <LinAlg/Matrix.h>`
31	`#include <LinAlg/Vector.h>`	31	`#include <LinAlg/Vector.h>`
32	`#include <LinAlg/LapackFunc.h>`	32	`#include <LinAlg/LapackFunc.h>`
33		33
34	`using namespace std;`	34	`using namespace std;`
35	`using namespace HMesh;`	35	`using namespace HMesh;`
36	`using namespace LinAlg;`	36	`using namespace LinAlg;`
37	`using namespace CGLA;`	37	`using namespace CGLA;`
38		38
39	`namespace {`	39	`namespace {`
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(Manifold& m, vector<T>& vec, int smooth_steps)`	44	`void smooth_something_on_mesh(Manifold& m, vector<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	`{`	47	`{`
48	`vector<T> new_vec(m.no_vertices());`	48	`vector<T> new_vec(m.no_vertices());`
49	`for(VertexIter vi=m.vertices_begin(); vi != m.vertices_end(); ++vi)`	49	`for(VertexIter vi=m.vertices_begin(); vi != m.vertices_end(); ++vi)`
50	`{`	50	`{`
51	`int i = vi->touched;`	51	`int i = vi->touched;`
52	`new_vec[i] = vec[i];`	52	`new_vec[i] = vec[i];`
53	`for(VertexCirculator vc(vi); !vc.end();++vc)`	53	`for(VertexCirculator vc(vi); !vc.end();++vc)`
54	`{`	54	`{`
55	`int j = vc.get_vertex()->touched;`	55	`int j = vc.get_vertex()->touched;`
56	`new_vec[i] += vec[j];`	56	`new_vec[i] += vec[j];`
57	`}`	57	`}`
58	`new_vec[i] /=(valency(vi)+ 1.0);`	58	`new_vec[i] /=(valency(vi)+ 1.0);`
59	`}`	59	`}`
60	`swap(vec,new_vec);`	60	`swap(vec,new_vec);`
61	`}`	61	`}`
62	`}`	62	`}`
63	`}`	63	`}`
64		64
65	`double voronoi_area(VertexIter v)`	65	`double voronoi_area(VertexIter v)`
66	`{`	66	`{`
67	`double area_mixed = 0;`	67	`double area_mixed = 0;`
68	`//For each triangle T from the 1-ring neighborhood of x`	68	`//For each triangle T from the 1-ring neighborhood of x`
69	`for(VertexCirculator vc(v); !vc.end(); ++vc)`	69	`for(VertexCirculator vc(v); !vc.end(); ++vc)`
70	`{`	70	`{`
71	`FaceIter f = vc.get_face();`	71	`FaceIter f = vc.get_face();`
72	`double f_area = area(f);`	72	`double f_area = area(f);`
73		73
74	`HalfEdgeIter he = vc.get_halfedge();`	74	`HalfEdgeIter he = vc.get_halfedge();`
75	`Vec3d v1(he->vert->pos);`	75	`Vec3d v1(he->vert->pos);`
76	`Vec3d v2(he->next->vert->pos);`	76	`Vec3d v2(he->next->vert->pos);`
77	`Vec3d v0(he->next->next->vert->pos);`	77	`Vec3d v0(he->next->next->vert->pos);`
78		78
79	`double a0 = acos(dot(v1-v0, v2-v0)/(length(v1-v0)*length(v2-v0)));`	79	`double a0 = acos(dot(v1-v0, v2-v0)/(length(v1-v0)*length(v2-v0)));`
80	`double a1 = acos(dot(v2-v1, v0-v1)/(length(v2-v1)*length(v0-v1)));`	80	`double a1 = acos(dot(v2-v1, v0-v1)/(length(v2-v1)*length(v0-v1)));`
81	`double a2 = acos(dot(v0-v2, v1-v2)/(length(v0-v2)*length(v1-v2)));`	81	`double a2 = acos(dot(v0-v2, v1-v2)/(length(v0-v2)*length(v1-v2)));`
82		82
83	`if(a0>(M_PI/2.0) && a1>(M_PI/2.0) && a2>(M_PI/2.0)) // f is non-obtuse`	83	`if(a0>(M_PI/2.0) && a1>(M_PI/2.0) && a2>(M_PI/2.0)) // f is non-obtuse`
84	`{`	84	`{`
85	`// Add Voronoi formula (see Section 3.3)`	85	`// Add Voronoi formula (see Section 3.3)`
86	`area_mixed += (1.0/8) *`	86	`area_mixed += (1.0/8) *`
87	`((1.0/tan(a1)) * sqr_length(v2-v0) +`	87	`((1.0/tan(a1)) * sqr_length(v2-v0) +`
88	`(1.0/tan(a2)) * sqr_length(v1-v0));`	88	`(1.0/tan(a2)) * sqr_length(v1-v0));`
89	`}`	89	`}`
90	`else // Voronoi inappropriate`	90	`else // Voronoi inappropriate`
91	`{`	91	`{`
92	`// Add either area(f)/4 or area(f)/2`	92	`// Add either area(f)/4 or area(f)/2`
93	`if(a0>M_PI/2.0)// the angle of f at x is obtuse`	93	`if(a0>M_PI/2.0)// the angle of f at x is obtuse`
94	`area_mixed += f_area/2;`	94	`area_mixed += f_area/2;`
95	`else`	95	`else`
96	`area_mixed += f_area/4;`	96	`area_mixed += f_area/4;`
97	`}`	97	`}`
98	`}`	98	`}`
99	`return area_mixed;`	99	`return area_mixed;`
100	`}`	100	`}`
101		101
102	`double barycentric_area(VertexIter v)`	102	`double barycentric_area(VertexIter v)`
103	`{`	103	`{`
104	`double barea = 0;`	104	`double barea = 0;`
105	`//For each triangle T from the 1-ring neighborhood of x`	105	`//For each triangle T from the 1-ring neighborhood of x`
106	`for(VertexCirculator vc(v); !vc.end(); ++vc)`	106	`for(VertexCirculator vc(v); !vc.end(); ++vc)`
107	`{`	107	`{`
108	`FaceIter f = vc.get_face();`	108	`FaceIter f = vc.get_face();`
109	`barea += area(f)/3.0;`	109	`barea += area(f)/3.0;`
110	`}`	110	`}`
111	`return barea;`	111	`return barea;`
112	`}`	112	`}`
113		113
114	`const Vec3d mean_curvature_normal(VertexIter v)`	114	`void unnormalized_mean_curvature_normal(VertexIter v, Vec3d& curv_normal, double& w_sum)`
115	`{`	115	`{`
116	`if(!is_boundary(v))`	116	`if(!is_boundary(v))`
117	`{`	117	`{`
118	`Vec3d vertex(v->pos);`	118	`Vec3d vertex(v->pos);`
119	`Vec3d curv_normal(0);`	119	`curv_normal = Vec3d(0);`
120	`double a_sum = 0;`	120	`w_sum = 0;`
121	`for(VertexCirculator vc(v); !vc.end(); ++vc)`	121	`for(VertexCirculator vc(v); !vc.end(); ++vc)`
122	`{`	122	`{`
123	`HalfEdgeIter h = vc.get_halfedge();`	123	`HalfEdgeIter h = vc.get_halfedge();`
124	`Vec3d nbr(h->vert->pos);`	124	`Vec3d nbr(h->vert->pos);`
125	`Vec3d left(h->next->vert->pos);`	125	`Vec3d left(h->next->vert->pos);`
126	`Vec3d right(h->opp->prev->opp->vert->pos);`	126	`Vec3d right(h->opp->prev->opp->vert->pos);`
127		127
128	`double d_left = dot(normalize(nbr-left),`	128	`double d_left = dot(normalize(nbr-left),`
129	`normalize(vertex-left));`	129	`normalize(vertex-left));`
130	`double d_right = dot(normalize(nbr-right),`	130	`double d_right = dot(normalize(nbr-right),`
131	`normalize(vertex-right));`	131	`normalize(vertex-right));`
132	`double a_left = acos(min(1.0, max(-1.0, d_left)));`	132	`double a_left = acos(min(1.0, max(-1.0, d_left)));`
133	`double a_right = acos(min(1.0, max(-1.0, d_right)));`	133	`double a_right = acos(min(1.0, max(-1.0, d_right)));`
134		134
135	`double w = 1.0/tan(a_left) + 1.0/tan(a_right);`	135	`double w = 1.0/tan(a_left) + 1.0/tan(a_right);`
136		136
137	`curv_normal += w * (vertex-nbr);`	137	`curv_normal += w * (vertex-nbr);`
138		-
139	`a_sum += area(vc.get_face());`	138	`w_sum += w;`
140	`}`	139	`}`
141		-
142	`return curv_normal/(4*a_sum);`	-
143	`}`	140	`}`
-		141	`}`
-		142
-		143	`const Vec3d mean_curvature_normal(VertexIter v)`
-		144	`{`
-		145	`Vec3d curv_normal;`
144	`return Vec3d(0);`	146	`double w_sum;`
-		147	`unnormalized_mean_curvature_normal(v, curv_normal, w_sum);`
-		148
-		149	`return curv_normal / (4*voronoi_area(v));`
145	`}`	150	`}`
146		151
147	`double sum_curvatures(Manifold& m, vector<double>& curvature)`	152	`double sum_curvatures(Manifold& m, vector<double>& curvature)`
148	`{`	153	`{`
149	`double sum = 0;`	154	`double sum = 0;`
150	`for(VertexIter v=m.vertices_begin(); v!=m.vertices_end(); ++v)`	155	`for(VertexIter v=m.vertices_begin(); v!=m.vertices_end(); ++v)`
151	`{`	156	`{`
152	`if(!is_boundary(v))`	157	`if(!is_boundary(v))`
153	`{`	158	`{`
154	`sum += curvature[v->touched] * voronoi_area(v);`	159	`sum += curvature[v->touched] * voronoi_area(v);`
155	`}`	160	`}`
156	`}`	161	`}`
157	`return sum;`	162	`return sum;`
158	`}`	163	`}`
159		164
160		165
161	`const double gaussian_curvature_angle_defect(VertexIter v)`	166	`const double gaussian_curvature_angle_defect(VertexIter v)`
162	`{`	167	`{`
163	`if(!is_boundary(v))`	168	`if(!is_boundary(v))`
164	`{`	169	`{`
165	`Vec3f vertex(v->pos);`	170	`Vec3f vertex(v->pos);`
166	`vector<Vec3d> edges;`	171	`vector<Vec3d> edges;`
167	`for(VertexCirculator vc(v); !vc.end(); ++vc)`	172	`for(VertexCirculator vc(v); !vc.end(); ++vc)`
168	`{`	173	`{`
169	`Vec3d e(normalize(vc.get_vertex()->pos-vertex));`	174	`Vec3d e(normalize(vc.get_vertex()->pos-vertex));`
170	`edges.push_back(e);`	175	`edges.push_back(e);`
171	`}`	176	`}`
172	`int N=edges.size();`	177	`int N=edges.size();`
173	`double angle_sum = 0;`	178	`double angle_sum = 0;`
174	`for(int i=0;i<N;++i)`	179	`for(int i=0;i<N;++i)`
175	`{`	180	`{`
176	`double dot_prod =`	181	`double dot_prod =`
177	`s_max(-1.0, s_min(1.0, dot(edges[i],edges[(i+1)%N])));`	182	`s_max(-1.0, s_min(1.0, dot(edges[i],edges[(i+1)%N])));`
178	`angle_sum += acos(dot_prod);`	183	`angle_sum += acos(dot_prod);`
179	`}`	184	`}`
180	`return (2*M_PI - angle_sum)/voronoi_area(v);`	185	`return (2*M_PI - angle_sum)/voronoi_area(v);`
181	`}`	186	`}`
182	`return 0;`	187	`return 0;`
183	`}`	188	`}`
184		189
185	`const Mat3x3d curvature_tensor(HalfEdgeIter h)`	190	`const Mat3x3d curvature_tensor(HalfEdgeIter h)`
186	`{`	191	`{`
187	`if(!is_boundary(h))`	192	`if(!is_boundary(h))`
188	`{`	193	`{`
189	`Vec3d edge(h->vert->pos - h->opp->vert->pos);`	194	`Vec3d edge(h->vert->pos - h->opp->vert->pos);`
190	`double edge_len = length(edge);`	195	`double edge_len = length(edge);`
191	`edge /= edge_len;`	196	`edge /= edge_len;`
192		197
193	`Vec3d h_norm(normal(h->face));`	198	`Vec3d h_norm(normal(h->face));`
194	`Vec3d h_opp_norm(normal(h->opp->face));`	199	`Vec3d h_opp_norm(normal(h->opp->face));`
195		200
196	`Vec3d nc = cross(h_norm, h_opp_norm);`	201	`Vec3d nc = cross(h_norm, h_opp_norm);`
197		202
198	`double sign = (dot(nc, edge) >= 0) ? 1 : -1;`	203	`double sign = (dot(nc, edge) >= 0) ? 1 : -1;`
199	`double beta = asin(nc.length());`	204	`double beta = asin(nc.length());`
200		205
201	`Mat3x3d m;`	206	`Mat3x3d m;`
202	`outer_product(edge, edge, m);`	207	`outer_product(edge, edge, m);`
203	`return sign * edge_len * beta * m;`	208	`return sign * edge_len * beta * m;`
204	`}`	209	`}`
205	`return Mat3x3d(0);`	210	`return Mat3x3d(0);`
206	`}`	211	`}`
207		212
208	`const Mat3x3d curvature_tensor_from_edges(VertexIter v)`	213	`const Mat3x3d curvature_tensor_from_edges(VertexIter v)`
209	`{`	214	`{`
210	`Mat3x3d curv_tensor(0);`	215	`Mat3x3d curv_tensor(0);`
211		216
212	`if(!is_boundary(v))`	217	`if(!is_boundary(v))`
213	`{`	218	`{`
214	`for(VertexCirculator vc(v); !vc.end(); ++vc)`	219	`for(VertexCirculator vc(v); !vc.end(); ++vc)`
215	`{`	220	`{`
216	`curv_tensor += 0.5*curvature_tensor(vc.get_halfedge());`	221	`curv_tensor += 0.5*curvature_tensor(vc.get_halfedge());`
217	`}`	222	`}`
218	`curv_tensor /= voronoi_area(v);`	223	`curv_tensor /= voronoi_area(v);`
219	`}`	224	`}`
220	`return curv_tensor;`	225	`return curv_tensor;`
221	`}`	226	`}`
222		227
223		228
224	`void curvature_tensor_paraboloid(VertexIter v,`	229	`void curvature_tensor_paraboloid(VertexIter v,`
225	`Mat2x2d& curv_tensor,`	230	`Mat2x2d& curv_tensor,`
226	`Mat3x3d& frame)`	231	`Mat3x3d& frame)`
227	`{`	232	`{`
228	`if(!is_boundary(v))`	233	`if(!is_boundary(v))`
229	`{`	234	`{`
230	`// First estimate the normal and compute a transformation matrix`	235	`// First estimate the normal and compute a transformation matrix`
231	`// which takes us into tangent plane coordinates.`	236	`// which takes us into tangent plane coordinates.`
232	`Vec3d Norm = Vec3d(normal(v));`	237	`Vec3d Norm = Vec3d(normal(v));`
233	`Vec3d X,Y;`	238	`Vec3d X,Y;`
234	`orthogonal(Norm,X,Y);`	239	`orthogonal(Norm,X,Y);`
235	`frame = Mat3x3d(X,Y,Norm);`	240	`frame = Mat3x3d(X,Y,Norm);`
236	`Vec3d centre(v->pos);`	241	`Vec3d centre(v->pos);`
237		242
238	`vector<Vec3d> points;`	243	`vector<Vec3d> points;`
239	`for(VertexCirculator vc(v); !vc.end(); ++vc)`	244	`for(VertexCirculator vc(v); !vc.end(); ++vc)`
240	`points.push_back(Vec3d(vc.get_vertex()->pos));`	245	`points.push_back(Vec3d(vc.get_vertex()->pos));`
241		246
242	`int N = points.size();`	247	`int N = points.size();`
243		248
244	`CVector b(N);`	249	`CVector b(N);`
245	`// Compute the matrix of parameter values`	250	`// Compute the matrix of parameter values`
246	`CMatrix PMat(N, 3);`	251	`CMatrix PMat(N, 3);`
247	`for(int i=0;i<N;++i)`	252	`for(int i=0;i<N;++i)`
248	`{`	253	`{`
249	`Vec3d p = frame * (points[i]-centre);`	254	`Vec3d p = frame * (points[i]-centre);`
250	`b[i] = p[2];`	255	`b[i] = p[2];`
251		256
252	`PMat.set(i,0,0.5*sqr(p[0]));`	257	`PMat.set(i,0,0.5*sqr(p[0]));`
253	`PMat.set(i,1,p[0]*p[1]);`	258	`PMat.set(i,1,p[0]*p[1]);`
254	`PMat.set(i,2,0.5*sqr(p[1]));`	259	`PMat.set(i,2,0.5*sqr(p[1]));`
255	`}`	260	`}`
256		261
257	`// Compute the coefficients of the polynomial surface`	262	`// Compute the coefficients of the polynomial surface`
258	`CVector x(3);`	263	`CVector x(3);`
259	`x = LinearLSSolve(PMat,b);`	264	`x = LinearLSSolve(PMat,b);`
260	`if(isnan(x[0])) cout << __LINE__ << " " << PMat << b << endl ;`	265	`if(isnan(x[0])) cout << __LINE__ << " " << PMat << b << endl ;`
261		266
262	`// Finally compute the shape tensor from the coefficients`	267	`// Finally compute the shape tensor from the coefficients`
263	`// using the first and second fundamental forms.`	268	`// using the first and second fundamental forms.`
264	`curv_tensor =Mat2x2d(x[0],x[1],x[1],x[2]);`	269	`curv_tensor =Mat2x2d(x[0],x[1],x[1],x[2]);`
265	`}`	270	`}`
266	`}`	271	`}`
267		272
268	`void curvature_tensors_from_edges(Manifold& m,`	273	`void curvature_tensors_from_edges(Manifold& m,`
269	`vector<Mat3x3d>& curvature_tensors)`	274	`vector<Mat3x3d>& curvature_tensors)`
270	`{`	275	`{`
271	`curvature_tensors.resize(m.no_vertices());`	276	`curvature_tensors.resize(m.no_vertices());`
272	`int i=0;`	277	`int i=0;`
273	`for(VertexIter vi=m.vertices_begin(); vi != m.vertices_end(); ++vi,++i)`	278	`for(VertexIter vi=m.vertices_begin(); vi != m.vertices_end(); ++vi,++i)`
274	`{`	279	`{`
275	`vi->touched = i;`	280	`vi->touched = i;`
276	`curvature_tensors[i] = curvature_tensor_from_edges(vi);`	281	`curvature_tensors[i] = curvature_tensor_from_edges(vi);`
277	`}`	282	`}`
278	`}`	283	`}`
279		284
280	`void smooth_curvature_tensors(Manifold& m,`	285	`void smooth_curvature_tensors(Manifold& m,`
281	`vector<Mat3x3d>& curvature_tensors)`	286	`vector<Mat3x3d>& curvature_tensors)`
282	`{`	287	`{`
283	`assert(curvature_tensors.size() == m.no_vertices());`	288	`assert(curvature_tensors.size() == m.no_vertices());`
284	`vector<Mat3x3d> tmp_curvature_tensors(m.no_vertices());`	289	`vector<Mat3x3d> tmp_curvature_tensors(m.no_vertices());`
285	`double tmp_area;`	290	`double tmp_area;`
286	`int i=0;`	291	`int i=0;`
287	`for(VertexIter vi=m.vertices_begin(); vi != m.vertices_end(); ++vi,++i)`	292	`for(VertexIter vi=m.vertices_begin(); vi != m.vertices_end(); ++vi,++i)`
288	`if(!is_boundary(vi))`	293	`if(!is_boundary(vi))`
289	`{`	294	`{`
290	`double a = voronoi_area(vi);`	295	`double a = voronoi_area(vi);`
291	`tmp_curvature_tensors[i] = curvature_tensors[i] * a;`	296	`tmp_curvature_tensors[i] = curvature_tensors[i] * a;`
292	`tmp_area = a;`	297	`tmp_area = a;`
293	`int n=1;`	298	`int n=1;`
294	`for(VertexCirculator vc(vi); !vc.end(); ++vc)`	299	`for(VertexCirculator vc(vi); !vc.end(); ++vc)`
295	`if(!is_boundary(vc.get_vertex()))`	300	`if(!is_boundary(vc.get_vertex()))`
296	`{`	301	`{`
297	`int j = vc.get_vertex()->touched;`	302	`int j = vc.get_vertex()->touched;`
298	`double a = voronoi_area(vc.get_vertex());`	303	`double a = voronoi_area(vc.get_vertex());`
299	`tmp_curvature_tensors[i] += curvature_tensors[j]*a;`	304	`tmp_curvature_tensors[i] += curvature_tensors[j]*a;`
300	`tmp_area += a;`	305	`tmp_area += a;`
301	`++n;`	306	`++n;`
302	`}`	307	`}`
303	`tmp_curvature_tensors[i] /= tmp_area;`	308	`tmp_curvature_tensors[i] /= tmp_area;`
304	`}`	309	`}`
305	`curvature_tensors = tmp_curvature_tensors;`	310	`curvature_tensors = tmp_curvature_tensors;`
306	`}`	311	`}`
307		312
308	`void gaussian_curvature_angle_defects(Manifold& m,`	313	`void gaussian_curvature_angle_defects(Manifold& m,`
309	`vector<double>& curvature,`	314	`vector<double>& curvature,`
310	`int smooth_steps)`	315	`int smooth_steps)`
311	`{`	316	`{`
312	`m.enumerate_vertices();`	317	`m.enumerate_vertices();`
313	`curvature.resize(m.no_vertices());`	318	`curvature.resize(m.no_vertices());`
314	`for(VertexIter vi=m.vertices_begin(); vi != m.vertices_end(); ++vi)`	319	`for(VertexIter vi=m.vertices_begin(); vi != m.vertices_end(); ++vi)`
315	`{`	320	`{`
316	`int i = vi->touched;`	321	`int i = vi->touched;`
317	`curvature[i] = gaussian_curvature_angle_defect(vi);`	322	`curvature[i] = gaussian_curvature_angle_defect(vi);`
318	`}`	323	`}`
319	`smooth_something_on_mesh(m, curvature, smooth_steps);`	324	`smooth_something_on_mesh(m, curvature, smooth_steps);`
320	`}`	325	`}`
321		326
322	`void mean_curvatures(Manifold& m, vector<double>& curvature,int smooth_steps)`	327	`void mean_curvatures(Manifold& m, vector<double>& curvature,int smooth_steps)`
323	`{`	328	`{`
324	`curvature.resize(m.no_vertices());`	329	`curvature.resize(m.no_vertices());`
325	`int i=0;`	330	`int i=0;`
326	`for(VertexIter vi=m.vertices_begin(); vi != m.vertices_end(); ++vi,++i)`	331	`for(VertexIter vi=m.vertices_begin(); vi != m.vertices_end(); ++vi,++i)`
327	`{`	332	`{`
328	`vi->touched = i;`	333	`vi->touched = i;`
329	`Vec3d N = mean_curvature_normal(vi);`	334	`Vec3d N = mean_curvature_normal(vi);`
330	`curvature[i] = length(N)*sign(dot(N,Vec3d(normal(vi))));`	335	`curvature[i] = length(N)*sign(dot(N,Vec3d(normal(vi))));`
331	`}`	336	`}`
332	`smooth_something_on_mesh(m, curvature, smooth_steps);`	337	`smooth_something_on_mesh(m, curvature, smooth_steps);`
333	`}`	338	`}`
334		339
335		340
336	`void curvature_paraboloids(Manifold& m,`	341	`void curvature_paraboloids(Manifold& m,`
337	`vector<Vec3d>& min_curv_direction,`	342	`vector<Vec3d>& min_curv_direction,`
338	`vector<Vec3d>& max_curv_direction,`	343	`vector<Vec3d>& max_curv_direction,`
339	`vector<double>& curvature)`	344	`vector<double>& curvature)`
340	`{`	345	`{`
341	`int i=0;`	346	`int i=0;`
342	`for(VertexIter vi=m.vertices_begin(); vi != m.vertices_end(); ++vi,++i)`	347	`for(VertexIter vi=m.vertices_begin(); vi != m.vertices_end(); ++vi,++i)`
343	`{`	348	`{`
344	`vi->touched = i;`	349	`vi->touched = i;`
345	`Mat2x2d tensor;`	350	`Mat2x2d tensor;`
346	`Mat3x3d frame;`	351	`Mat3x3d frame;`
347	`curvature_tensor_paraboloid(vi, tensor, frame);`	352	`curvature_tensor_paraboloid(vi, tensor, frame);`
348		353
349	`Mat2x2d Q,L;`	354	`Mat2x2d Q,L;`
350	`int s = power_eigensolution(tensor, Q, L);`	355	`int s = power_eigensolution(tensor, Q, L);`
351		356
352	`if(s<2)`	357	`if(s<2)`
353	`cout << tensor << Q << L << endl;`	358	`cout << tensor << Q << L << endl;`
354		359
355	`int max_idx=0;`	360	`int max_idx=0;`
356	`int min_idx=1;`	361	`int min_idx=1;`
357		362
358	`if(fabs(L[max_idx][max_idx])<fabs(L[min_idx][min_idx])) swap(max_idx, min_idx);`	363	`if(fabs(L[max_idx][max_idx])<fabs(L[min_idx][min_idx])) swap(max_idx, min_idx);`
359		364
360	`Mat3x3d frame_t = transpose(frame);`	365	`Mat3x3d frame_t = transpose(frame);`
361		366
362	`max_curv_direction[i] =`	367	`max_curv_direction[i] =`
363	`frame_t * Vec3d(Q[max_idx][0], Q[max_idx][1], 0);`	368	`frame_t * Vec3d(Q[max_idx][0], Q[max_idx][1], 0);`
364		369
365	`min_curv_direction[i] =`	370	`min_curv_direction[i] =`
366	`frame_t * Vec3d(Q[min_idx][0], Q[min_idx][1], 0);`	371	`frame_t * Vec3d(Q[min_idx][0], Q[min_idx][1], 0);`
367		372
368	`curvature[i] = L[0][0]*L[1][1];`	373	`curvature[i] = L[0][0]*L[1][1];`
369	`}`	374	`}`
370	`}`	375	`}`
371		376
372		377
373	`void curvature_from_tensors(Manifold& m,`	378	`void curvature_from_tensors(Manifold& m,`
374	`const vector<Mat3x3d>& curvature_tensors,`	379	`const vector<Mat3x3d>& curvature_tensors,`
375	`vector<Vec3d>& min_curv_direction,`	380	`vector<Vec3d>& min_curv_direction,`
376	`vector<Vec3d>& max_curv_direction,`	381	`vector<Vec3d>& max_curv_direction,`
377	`vector<double>& curvature)`	382	`vector<double>& curvature)`
378	`{`	383	`{`
379	`assert(curvature_tensors.size() == m.no_vertices());`	384	`assert(curvature_tensors.size() == m.no_vertices());`
380	`min_curv_direction.resize(m.no_vertices());`	385	`min_curv_direction.resize(m.no_vertices());`
381	`max_curv_direction.resize(m.no_vertices());`	386	`max_curv_direction.resize(m.no_vertices());`
382	`curvature.resize(m.no_vertices());`	387	`curvature.resize(m.no_vertices());`
383	`double max_val = -1e30;`	388	`double max_val = -1e30;`
384	`for(VertexIter vi=m.vertices_begin(); vi != m.vertices_end(); ++vi)`	389	`for(VertexIter vi=m.vertices_begin(); vi != m.vertices_end(); ++vi)`
385	`{`	390	`{`
386	`int i = vi->touched;`	391	`int i = vi->touched;`
387	`Mat3x3d C,Q,L;`	392	`Mat3x3d C,Q,L;`
388	`C = curvature_tensors[i];`	393	`C = curvature_tensors[i];`
389	`int s = power_eigensolution(C, Q, L);`	394	`int s = power_eigensolution(C, Q, L);`
390	`Vec3d dmin, dmax;`	395	`Vec3d dmin, dmax;`
391	`if(s==0)`	396	`if(s==0)`
392	`{`	397	`{`
393	`Vec3d n(normal(vi));`	398	`Vec3d n(normal(vi));`
394	`orthogonal(n,dmin, dmax);`	399	`orthogonal(n,dmin, dmax);`
395	`curvature[i] = 0;`	400	`curvature[i] = 0;`
396	`cout << " rank 0 " << endl;`	401	`cout << " rank 0 " << endl;`
397	`}`	402	`}`
398	`else if(s==1)`	403	`else if(s==1)`
399	`{`	404	`{`
400	`Vec3d n(normal(vi));`	405	`Vec3d n(normal(vi));`
401	`dmin = normalize(Q[0]);`	406	`dmin = normalize(Q[0]);`
402	`dmax = cross(n, dmin);`	407	`dmax = cross(n, dmin);`
403	`curvature[i] = 0;`	408	`curvature[i] = 0;`
404	`cout << " rank 1 " << endl;`	409	`cout << " rank 1 " << endl;`
405	`}`	410	`}`
406	`else`	411	`else`
407	`{`	412	`{`
408	`/* Vec3d l(fabs(L[0][0]), fabs(L[1][1]), fabs(L[2][2]));`	413	`/* Vec3d l(fabs(L[0][0]), fabs(L[1][1]), fabs(L[2][2]));`
409		414
410	`int z_idx=2;`	415	`int z_idx=2;`
411	`if(s==3)`	416	`if(s==3)`
412	`{`	417	`{`
413	`if(l[0] < l[1])`	418	`if(l[0] < l[1])`
414	`z_idx = l[0]<l[2] ? 0 : 2;`	419	`z_idx = l[0]<l[2] ? 0 : 2;`
415	`else`	420	`else`
416	`z_idx = l[1]<l[2] ? 1 : 2;`	421	`z_idx = l[1]<l[2] ? 1 : 2;`
417	`}`	422	`}`
418	`int max_idx = (z_idx + 1) % 3;`	423	`int max_idx = (z_idx + 1) % 3;`
419	`int min_idx = (z_idx + 2) % 3;`	424	`int min_idx = (z_idx + 2) % 3;`
420		425
421	`if(l[max_idx] < l[min_idx]) swap(max_idx, min_idx);`	426	`if(l[max_idx] < l[min_idx]) swap(max_idx, min_idx);`
422	`*/`	427	`*/`
423	`int max_idx = 0;`	428	`int max_idx = 0;`
424	`int min_idx = 1;`	429	`int min_idx = 1;`
425	`// Yes - the biggest eigenvalue corresponds to the min direction`	430	`// Yes - the biggest eigenvalue corresponds to the min direction`
426	`// and vice versa.`	431	`// and vice versa.`
427	`dmin = normalize(Q[max_idx]);`	432	`dmin = normalize(Q[max_idx]);`
428	`dmax = normalize(Q[min_idx]);`	433	`dmax = normalize(Q[min_idx]);`
429		434
430	`curvature[i] = L[max_idx][max_idx]*L[min_idx][min_idx];`	435	`curvature[i] = L[max_idx][max_idx]*L[min_idx][min_idx];`
431		436
432	`}`	437	`}`
433	`min_curv_direction[i] = dmin;`	438	`min_curv_direction[i] = dmin;`
434	`max_curv_direction[i] = dmax;`	439	`max_curv_direction[i] = dmax;`
435	`max_val = max(fabs(curvature[i]), max_val);`	440	`max_val = max(fabs(curvature[i]), max_val);`
436		441
437	`}`	442	`}`
438	`scal = 1.0/max_val;`	443	`scal = 1.0/max_val;`
439	`}`	444	`}`
440		445
441		446

Subversion Repositories gelsvn

(root)/trunk/apps/MeshEdit/curvature.cpp @ 488 – Rev 403 → 412