Rev 489 | Blame | Last modification | View Log | RSS feed
/*
* 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/Manifold.h>
#include <HMesh/VertexHandle.h>
#include <HMesh/FaceHandle.h>
#include <HMesh/HalfEdgeHandle.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/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 {
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(VertexHandle v = m.vertices_begin(); v != m.vertices_end(); ++v){
int i = v.get().touched;
new_vec[i] = vec[i];
for(VertexCirculator vc(v); !vc.end(); ++vc){
int j = vc.vertex().get().touched;
new_vec[i] += vec[j];
}
new_vec[i] /= (valency(v) + 1.0);
}
swap(vec,new_vec);
}
}
}
double voronoi_area(const VertexHandle& v)
{
double area_mixed = 0;
//For each triangle T from the 1-ring neighborhood of x
for(VertexCirculator vc(v); !vc.end(); ++vc){
const FaceHandle& f = vc.face();
double f_area = area(f);
const HalfEdgeHandle he = vc.halfedge();
Vec3d v1(he.vert().get().pos);
Vec3d v2(he.next().vert().get().pos);
Vec3d v0(he.next().next().vert().get().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(const VertexHandle& v)
{
double barea = 0;
//For each triangle T from the 1-ring neighborhood of x
for(VertexCirculator vc(v); !vc.end(); ++vc){
const FaceHandle f = vc.face();
barea += area(f)/3.0;
}
return barea;
}
void unnormalized_mean_curvature_normal(const VertexHandle& v, Vec3d& curv_normal, double& w_sum)
{
if(!is_boundary(v)){
Vec3d vertex(v.get().pos);
curv_normal = Vec3d(0);
w_sum = 0;
for(VertexCirculator vc(v); !vc.end(); ++vc){
const HalfEdgeHandle h = vc.halfedge();
Vec3d nbr(h.vert().get().pos);
Vec3d left(h.next().vert().get().pos);
Vec3d right(h.opp().prev().opp().vert().get().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(const VertexHandle& 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(VertexHandle v = m.vertices_begin(); v != m.vertices_end(); ++v){
if(is_boundary(v))
continue;
sum += curvature[v.get().touched] * voronoi_area(v);
}
return sum;
}
const double gaussian_curvature_angle_defect(const VertexHandle& v)
{
if(!is_boundary(v)){
Vec3f vertex(v.get().pos);
vector<Vec3d> edges;
for(VertexCirculator vc(v); !vc.end(); ++vc){
Vec3d e(normalize(vc.vertex().get().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(const HalfEdgeHandle& h)
{
if(!is_boundary(h)){
Vec3d edge(h.vert().get().pos - h.opp().vert().get().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(const VertexHandle& v)
{
Mat3x3d curv_tensor(0);
if(!is_boundary(v)){
for(VertexCirculator vc(v); !vc.end(); ++vc){
curv_tensor += 0.5*curvature_tensor(vc.halfedge());
}
curv_tensor /= voronoi_area(v);
}
return curv_tensor;
}
void curvature_tensor_paraboloid(const VertexHandle& 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.get().pos);
vector<Vec3d> points;
for(VertexCirculator vc(v); !vc.end(); ++vc){
points.push_back(Vec3d(vc.vertex().get().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());
size_t i = 0;
for(VertexHandle v = m.vertices_begin(); v != m.vertices_end(); ++v, ++i){
v.get().touched = i;
curvature_tensors[i] = curvature_tensor_from_edges(v);
}
}
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;
size_t i = 0;
for(VertexHandle v = m.vertices_begin(); v != m.vertices_end(); ++v, ++i){
if(is_boundary(v))
continue;
double a = voronoi_area(v);
tmp_curvature_tensors[i] = curvature_tensors[i] * a;
tmp_area = a;
int n = 1;
for(VertexCirculator vc(v); !vc.end(); ++vc){
VertexHandle v1 = vc.vertex();
if(!is_boundary(v1)){
int j = v1.get().touched;
double a = voronoi_area(v1);
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(VertexHandle v = m.vertices_begin(); v != m.vertices_end(); ++v)
{
int i = v.get().touched;
curvature[i] = gaussian_curvature_angle_defect(v);
}
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(VertexHandle v = m.vertices_begin(); v != m.vertices_end(); ++v, ++i)
{
v.get().touched = i;
Vec3d N = mean_curvature_normal(v);
curvature[i] = length(N)*sign(dot(N,Vec3d(normal(v))));
}
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(VertexHandle v = m.vertices_begin(); v != m.vertices_end(); ++v, ++i)
{
v.get().touched = i;
Mat2x2d tensor;
Mat3x3d frame;
curvature_tensor_paraboloid(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[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(VertexHandle v = m.vertices_begin(); v != m.vertices_end(); ++v)
{
int i = v.get().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(v ));
orthogonal(n, dmin, dmax);
curvature[i] = 0;
cout << " rank 0 " << endl;
}
else if(s == 1)
{
Vec3d n(normal(v));
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;
}