Subversion Repositories gelsvn

/*

 *  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;

}

const Vec3d mean_curvature_normal(VertexIter v)

{

	if(!is_boundary(v))

	{

		Vec3d vertex(v->pos);

		Vec3d curv_normal(0);

		double a_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);

			a_sum += area(vc.get_face());

		}

		return curv_normal/(4*a_sum);

	}

	return Vec3d(0);

}

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;

}