principal_curvature.cpp

// Code from igl library
// Copyright (C) 2013 Alec Jacobson <[email protected]>
// 
// This Source Code Form is subject to the terms of the Mozilla Public License 
// v. 2.0. If a copy of the MPL was not distributed with this file, You can 
// obtain one at http://mozilla.org/MPL/2.0/.
#pragma once

#include <Eigen/Geometry>
#include <Eigen/Dense>

#include <vector>
#include <stdio.h>
#include <map>
#include <iostream>
#include <fstream>
#include <iomanip>
#include <iostream>
#include <queue>
#include <list>
#include <cmath>
#include <limits>

#include <Eigen/SparseCholesky>

// Compute the principal curvature directions and magnitude of the given triangle mesh
//   DerivedV derived from vertex positions matrix type: i.e. MatrixXd
//   DerivedF derived from face indices matrix type: i.e. MatrixXi
// Inputs:
//   V       eigen matrix #V by 3
//   F       #F by 3 list of mesh faces (must be triangles)
//   radius  controls the size of the neighbourhood used, 1 = average edge lenght
//
// Outputs:
//   PD1 #V by 3 maximal curvature direction for each vertex.
//   PD2 #V by 3 minimal curvature direction for each vertex.  
//  Note that to maximal/minimal curvature value is the rowise norm of PD1/PD2
//
// See also: moveVF, moveFV
//
// This function has been developed by: Nikolas De Giorgis, Luigi Rocca and Enrico Puppo.
// The algorithm is based on:
// Efficient Multi-scale Curvature and Crease Estimation
// Daniele Panozzo, Enrico Puppo, Luigi Rocca
// GraVisMa, 2010
namespace igl{

typedef enum
{
	SPHERE_SEARCH,
	K_RING_SEARCH
} searchType;

typedef enum
{
	AVERAGE,
	PROJ_PLANE
} normalType;

class CurvatureCalculator
{
public:
	/* Row number i represents the i-th vertex, whose columns are:
	curv[i][0] : K2
	curv[i][1] : K1
	curvDir[i][0] : PD1
	curvDir[i][1] : PD2
	*/
	std::vector< std::vector<double> > curv;
	std::vector< std::vector<Eigen::Vector3d> > curvDir;
	bool curvatureComputed;
	class Quadric
	{
	public:

		inline Quadric ()
		{
			a() = b() = c() = d() = e() = 1.0;
		}

		inline Quadric(double av, double bv, double cv, double dv, double ev)
		{
			a() = av;
			b() = bv;
			c() = cv;
			d() = dv;
			e() = ev;
		}

		inline double& a() { return data[0];}
		inline double& b() { return data[1];}
		inline double& c() { return data[2];}
		inline double& d() { return data[3];}
		inline double& e() { return data[4];}

		double data[5];

		inline double evaluate(double u, double v)
		{
			return a()*u*u + b()*u*v + c()*v*v + d()*u + e()*v;
		}

		inline double du(double u, double v)
		{
			return 2.0*a()*u + b()*v + d();
		}

		inline double dv(double u, double v)
		{
			return 2.0*c()*v + b()*u + e();
		}

		inline double duv(double , double )
		{
			return b();
		}

		inline double duu(double , double )
		{
			return 2.0*a();
		}

		inline double dvv(double , double )
		{
			return 2.0*c();
		}


		inline static Quadric fit(std::vector<Eigen::Vector3d> &VV, bool , bool)
		{
			using namespace std;
			assert(VV.size() >= 5);
			if (VV.size() < 5)
			{
				cerr << "ASSERT FAILED!" << endl;
				exit(0);
			}

			Eigen::MatrixXd A(VV.size(),5);
			Eigen::MatrixXd b(VV.size(),1);
			Eigen::MatrixXd sol(5,1);

			for(unsigned int c=0; c < VV.size(); ++c)
			{
				double u = VV[c][0];
				double v = VV[c][1];
				double n = VV[c][2];

				A(c,0) = u*u;
				A(c,1) = u*v;
				A(c,2) = v*v;
				A(c,3) = u;
				A(c,4) = v;

				b(c) = n;
			}

			sol=A.jacobiSvd(Eigen::ComputeThinU | Eigen::ComputeThinV).solve(b);

			return Quadric(sol(0),sol(1),sol(2),sol(3),sol(4));
		}
	};

public:

	Eigen::MatrixXd vertices;
	// Face list of current mesh    (#F x 3) or (#F x 4)
	// The i-th row contains the indices of the vertices that forms the i-th face in ccw order
	Eigen::MatrixXi faces;

	std::vector<std::vector<int> > vertex_to_vertices;
	std::vector<std::vector<int> > vertex_to_faces;
	std::vector<std::vector<int> > vertex_to_faces_index;
	Eigen::MatrixXd face_normals;
	Eigen::MatrixXd vertex_normals;

	/* Size of the neighborhood */
	double sphereRadius;
	int kRing;

	bool localMode; /* Use local mode */
	bool projectionPlaneCheck; /* Check collected vertices on tangent plane */
	bool montecarlo;
	bool svd; /* Use svd calculation instead of pseudoinverse */
	bool zeroDetCheck; /* Check if the determinant is close to zero */
	unsigned int montecarloN;

	searchType st; /* Use either a sphere search or a k-ring search */
	normalType nt;

	double lastRadius;
	double scaledRadius;
	std::string lastMeshName;

	/* Benchmark related variables */
	bool expStep; /* True if we want the radius to increase exponentially */
	int step;  /* If expStep==false, by how much rhe radius increases on every step */
	int maxSize; /* The maximum limit of the radius in the benchmark */

	inline CurvatureCalculator();

	inline void finalEigenStuff (size_t, std::vector<Eigen::Vector3d>, Quadric );
	inline void fitQuadric (Eigen::Vector3d, std::vector<Eigen::Vector3d> ref, const  std::vector<int>& , Quadric *);
	inline void applyProjOnPlane(Eigen::Vector3d, std::vector<int>, std::vector<int>&);
	inline void getSphere(const size_t, const double, std::vector<int>&, int min);
	inline void getKRing(const size_t, const double,std::vector<int>&);
	inline Eigen::Vector3d project(Eigen::Vector3d, Eigen::Vector3d, Eigen::Vector3d);
	inline void computeReferenceFrame(size_t, Eigen::Vector3d, std::vector<Eigen::Vector3d>&);
	inline void getAverageNormal(size_t, std::vector<int>, Eigen::Vector3d&);
	inline void getProjPlane(size_t, std::vector<int>, Eigen::Vector3d&);
	inline void applyMontecarlo(std::vector<int>&,std::vector<int>*);
	inline void computeCurvature();
	inline void printCurvature(std::string outpath);
	inline double getAverageEdge();

	inline static int rotateForward (double *v0, double *v1, double *v2)
	{
		double t;

		if (abs(*v2) >= abs(*v1) && abs(*v2) >= abs(*v0))
			return 0;

		t = *v0;
		*v0 = *v2;
		*v2 = *v1;
		*v1 = t;

		return 1 + rotateForward (v0, v1, v2);
	}

	inline static void rotateBackward (int nr, double *v0, double *v1, double *v2)
	{
		double t;

		if (nr == 0)
			return;

		t = *v2;
		*v2 = *v0;
		*v0 = *v1;
		*v1 = t;

		rotateBackward (nr - 1, v0, v1, v2);
	}

	inline static Eigen::Vector3d chooseMax (Eigen::Vector3d n, Eigen::Vector3d abc, double ab)
	{
		int i, max_i;
		double max_sp;
		Eigen::Vector3d nt[8];

		n.normalize ();
		abc.normalize ();

		max_sp = - std::numeric_limits<double>::max();

		for (i = 0; i < 4; i++)
		{
			nt[i] = n;
			if (ab > 0)
			{
				switch (i)
				{
				case 0:
					break;

				case 1:
					nt[i][2] = -n[2];
					break;

				case 2:
					nt[i][0] = -n[0];
					nt[i][1] = -n[1];
					break;

				case 3:
					nt[i][0] = -n[0];
					nt[i][1] = -n[1];
					nt[i][2] = -n[2];
					break;
				}
			}
			else
			{
				switch (i)
				{
				case 0:
					nt[i][0] = -n[0];
					break;

				case 1:
					nt[i][1] = -n[1];
					break;

				case 2:
					nt[i][0] = -n[0];
					nt[i][2] = -n[2];
					break;

				case 3:
					nt[i][1] = -n[1];
					nt[i][2] = -n[2];
					break;
				}
			}

			if (nt[i].dot(abc) > max_sp)
			{
				max_sp = nt[i].dot(abc);
				max_i = i;
			}
		}

		return nt[max_i];
	}
};


class comparer
{
public:
	inline bool operator() (const std::pair<int, double>& lhs, const std::pair<int, double>&rhs) const
	{
		return lhs.second>rhs.second;
	}
};

inline CurvatureCalculator::CurvatureCalculator()
{
	this->localMode=false;
	this->projectionPlaneCheck=false;
	this->sphereRadius=5;
	this->st=SPHERE_SEARCH;
	this->nt=PROJ_PLANE;
	this->montecarlo=false;
	this->montecarloN=0;
	this->kRing=3;
	this->svd=true;
	this->zeroDetCheck=true;
	this->curvatureComputed=false;
	this->expStep=true;
}

inline void CurvatureCalculator::fitQuadric (Eigen::Vector3d v, std::vector<Eigen::Vector3d> ref, const std::vector<int>& vv, Quadric *q)
{
	std::vector<Eigen::Vector3d> points;
	points.reserve (vv.size());

	for (unsigned int i = 0; i < vv.size(); ++i) {

		Eigen::Vector3d  cp = vertices.row(vv[i]);

		// vtang non e` il v tangente!!!
		Eigen::Vector3d  vTang = cp - v;

		double x = vTang.dot(ref[0]);
		double y = vTang.dot(ref[1]);
		double z = vTang.dot(ref[2]);
		points.push_back(Eigen::Vector3d (x,y,z));
	}
	*q = Quadric::fit (points, zeroDetCheck, svd);
}

inline void CurvatureCalculator::finalEigenStuff (size_t i, std::vector<Eigen::Vector3d> ref, Quadric q)
{

	double a = q.a();
	double b = q.b();
	double c = q.c();
	double d = q.d();
	double e = q.e();

	//  if (fabs(a) < 10e-8 || fabs(b) < 10e-8)
	//  {
	//    std::cout << "Degenerate quadric: " << i << std::endl;
	//  }

	double E = 1.0 + d*d;
	double F = d*e;
	double G = 1.0 + e*e;

	Eigen::Vector3d n = Eigen::Vector3d(-d,-e,1.0).normalized();

	double L = 2.0 * a * n[2];
	double M = b * n[2];
	double N = 2 * c * n[2];


	// ----------------- Eigen stuff
	Eigen::Matrix2d m;
	m << L*G - M*F, M*E-L*F, M*E-L*F, N*E-M*F;
	m = m / (E*G-F*F);
	Eigen::SelfAdjointEigenSolver<Eigen::Matrix2d> eig(m);

	Eigen::Vector2d c_val = eig.eigenvalues();
	Eigen::Matrix2d c_vec = eig.eigenvectors();

	// std::cerr << "c_val:" << c_val << std::endl;
	// std::cerr << "c_vec:" << c_vec << std::endl;

	// std::cerr << "c_vec:" << c_vec(0) << " "  << c_vec(1) << std::endl;

	c_val = -c_val;

	Eigen::Vector3d v1, v2;
	v1[0] = c_vec(0);
	v1[1] = c_vec(1);
	v1[2] = 0; //d * v1[0] + e * v1[1];

	v2[0] = c_vec(2);
	v2[1] = c_vec(3);
	v2[2] = 0; //d * v2[0] + e * v2[1];


	// v1 = v1.normalized();
	// v2 = v2.normalized();

	Eigen::Vector3d v1global = ref[0] * v1[0] + ref[1] * v1[1] + ref[2] * v1[2];
	Eigen::Vector3d v2global = ref[0] * v2[0] + ref[1] * v2[1] + ref[2] * v2[2];

	v1global.normalize();
	v2global.normalize();

	v1global *= c_val(0);
	v2global *= c_val(1);

	if (c_val[0] > c_val[1])
	{
		curv[i]=std::vector<double>(2);
		curv[i][0]=c_val(1);
		curv[i][1]=c_val(0);
		curvDir[i]=std::vector<Eigen::Vector3d>(2);
		curvDir[i][0]=v2global;
		curvDir[i][1]=v1global;
	}
	else
	{
		curv[i]=std::vector<double>(2);
		curv[i][0]=c_val(0);
		curv[i][1]=c_val(1);
		curvDir[i]=std::vector<Eigen::Vector3d>(2);
		curvDir[i][0]=v1global;
		curvDir[i][1]=v2global;
	}
	// ---- end Eigen stuff
}

inline void CurvatureCalculator::getKRing(const size_t start, const double r, std::vector<int>&vv)
{
	size_t bufsize = vertices.rows();
	vv.reserve(bufsize);
	std::list<std::pair<int,int> > queue;
	bool* visited = (bool*)calloc(bufsize,sizeof(bool));
	queue.push_back(std::pair<int,int>((int)start,0));
	visited[start]=true;
	while (!queue.empty())
	{
		int toVisit=queue.front().first;
		int distance=queue.front().second;
		queue.pop_front();
		vv.push_back(toVisit);
		if (distance<(int)r)
		{
			for (unsigned int i=0; i<vertex_to_vertices[toVisit].size(); i++)
			{
				int neighbor=vertex_to_vertices[toVisit][i];
				if (!visited[neighbor])
				{
					queue.push_back(std::pair<int,int> (neighbor,distance+1));
					visited[neighbor]=true;
				}
			}
		}
	}
	free(visited);
	return;
}


inline void CurvatureCalculator::getSphere(const size_t start, const double r, std::vector<int> &vv, int min)
{
	size_t bufsize=vertices.rows();
	vv.reserve(bufsize);
	std::list<size_t>* queue= new std::list<size_t>();
	bool* visited = (bool*)calloc(bufsize,sizeof(bool));
	queue->push_back(start);
	visited[start]=true;
	Eigen::Vector3d me=vertices.row(start);
	std::priority_queue<std::pair<int, double>, std::vector<std::pair<int, double> >, comparer >* extra_candidates= new  std::priority_queue<std::pair<int, double>, std::vector<std::pair<int, double> >, comparer >();
	while (!queue->empty())
	{
		size_t toVisit=queue->front();
		queue->pop_front();
		vv.push_back((int)toVisit);
		for (unsigned int i=0; i<vertex_to_vertices[toVisit].size(); i++)
		{
			int neighbor=vertex_to_vertices[toVisit][i];
			if (!visited[neighbor])
			{
				Eigen::Vector3d neigh=vertices.row(neighbor);
				double distance=(me-neigh).norm();
				if (distance<r)
					queue->push_back(neighbor);
				else if ((int)vv.size()<min)
					extra_candidates->push(std::pair<int,double>(neighbor,distance));
				visited[neighbor]=true;
			}
		}
	}
	while (!extra_candidates->empty() && (int)vv.size()<min)
	{
		std::pair<int, double> cand=extra_candidates->top();
		extra_candidates->pop();
		vv.push_back(cand.first);
		for (unsigned int i=0; i<vertex_to_vertices[cand.first].size(); i++)
		{
			int neighbor=vertex_to_vertices[cand.first][i];
			if (!visited[neighbor])
			{
				Eigen::Vector3d neigh=vertices.row(neighbor);
				double distance=(me-neigh).norm();
				extra_candidates->push(std::pair<int,double>(neighbor,distance));
				visited[neighbor]=true;
			}
		}
	}
	free(extra_candidates);
	free(queue);
	free(visited);
}

inline Eigen::Vector3d CurvatureCalculator::project(Eigen::Vector3d v, Eigen::Vector3d  vp, Eigen::Vector3d ppn)
{
	return (vp - (ppn * ((vp - v).dot(ppn))));
}

inline void CurvatureCalculator::computeReferenceFrame(size_t i, Eigen::Vector3d normal, std::vector<Eigen::Vector3d>& ref )
{

	Eigen::Vector3d longest_v=Eigen::Vector3d::Zero();
	longest_v=Eigen::Vector3d(vertices.row(vertex_to_vertices[i][0]));

	longest_v=(project(vertices.row(i),longest_v,normal)-Eigen::Vector3d(vertices.row(i))).normalized();

	/* L'ultimo asse si ottiene come prodotto vettoriale tra i due
	* calcolati */
	Eigen::Vector3d y_axis=(normal.cross(longest_v)).normalized();
	ref[0]=longest_v;
	ref[1]=y_axis;
	ref[2]=normal;
}

inline void CurvatureCalculator::getAverageNormal(size_t j, std::vector<int> vv, Eigen::Vector3d& normal)
{
	normal=(vertex_normals.row(j)).normalized();
	if (localMode)
		return;

	for (unsigned int i=0; i<vv.size(); i++)
	{
		normal+=vertex_normals.row(vv[i]).normalized();
	}
	normal.normalize();
}

inline void CurvatureCalculator::getProjPlane(size_t j, std::vector<int> vv, Eigen::Vector3d& ppn)
{
	int nr;
	double a, b, c;
	double nx, ny, nz;
	double abcq;

	a = b = c = 0;

	if (localMode)
	{
		for (unsigned int i=0; i<vertex_to_faces.at(j).size(); ++i)
		{
			Eigen::Vector3d faceNormal=face_normals.row(vertex_to_faces.at(j).at(i));
			a += faceNormal[0];
			b += faceNormal[1];
			c += faceNormal[2];
		}
	}
	else
	{
		for (unsigned int i=0; i<vv.size(); ++i)
		{
			a+= vertex_normals.row(vv[i])[0];
			b+= vertex_normals.row(vv[i])[1];
			c+= vertex_normals.row(vv[i])[2];
		}
	}
	nr = rotateForward (&a, &b, &c);
	abcq = a*a + b*b + c*c;
	nx = sqrt (a*a / abcq);
	ny = sqrt (b*b / abcq);
	nz = sqrt (1 - nx*nx - ny*ny);
	rotateBackward (nr, &a, &b, &c);
	rotateBackward (nr, &nx, &ny, &nz);

	ppn = chooseMax (Eigen::Vector3d(nx, ny, nz), Eigen::Vector3d (a, b, c), a * b);
	ppn.normalize();
}


inline double CurvatureCalculator::getAverageEdge()
{
	double sum = 0;
	int count = 0;

	for (int i = 0; i<faces.rows(); i++)
	{
		for (short unsigned j=0; j<3; j++)
		{
			Eigen::Vector3d p1=vertices.row(faces.row(i)[j]);
			Eigen::Vector3d p2=vertices.row(faces.row(i)[(j+1)%3]);

			double l = (p1-p2).norm();

			sum+=l;
			++count;
		}
	}

	return (sum/(double)count);
}


inline void CurvatureCalculator::applyProjOnPlane(Eigen::Vector3d ppn, std::vector<int> vin, std::vector<int> &vout)
{
	for (std::vector<int>::iterator vpi = vin.begin(); vpi != vin.end(); ++vpi)
		if (vertex_normals.row(*vpi) * ppn > 0.0f)
			vout.push_back (*vpi);
}

inline void CurvatureCalculator::applyMontecarlo(std::vector<int>& vin, std::vector<int> *vout)
{
	if (montecarloN >= vin.size ())
	{
		*vout = vin;
		return;
	}

	float p = ((float) montecarloN) / (float) vin.size();
	for (std::vector<int>::iterator vpi = vin.begin(); vpi != vin.end(); ++vpi)
	{
		float r;
		if ((r = ((float)rand () / RAND_MAX)) < p)
		{
			vout->push_back (*vpi);
		}
	}
}

inline void CurvatureCalculator::computeCurvature()
{
	using namespace std;

	//CHECK che esista la mesh
	size_t vertices_count=vertices.rows() ;

	if (vertices_count <=0)
		return;

	curvDir=std::vector< std::vector<Eigen::Vector3d> >(vertices_count);
	curv=std::vector<std::vector<double> >(vertices_count);


	scaledRadius=getAverageEdge()*sphereRadius;

	std::vector<int> vv;
	std::vector<int> vvtmp;
	Eigen::Vector3d normal;

	//double time_spent;
	//double searchtime=0, ref_time=0, fit_time=0, final_time=0;

	for (size_t i=0; i<vertices_count; ++i)
	{
		vv.clear();
		vvtmp.clear();
		Eigen::Vector3d me=vertices.row(i);
		switch (st)
		{
		case SPHERE_SEARCH:
			getSphere(i,scaledRadius,vv,6);
			break;
		case K_RING_SEARCH:
			getKRing(i,kRing,vv);
			break;
		default:
			fprintf(stderr,"Error: search type not recognized");
			return;
		}

		std::vector<Eigen::Vector3d> ref(3);
		if (vv.size()<6)
		{
			std::cerr << "Could not compute curvature of radius " << scaledRadius << endl;
			return;
		}
		switch (nt)
		{
		case AVERAGE:
			getAverageNormal(i,vv,normal);
			break;
		case PROJ_PLANE:
			getProjPlane(i,vv,normal);
			break;
		default:
			fprintf(stderr,"Error: normal type not recognized");
			return;
		}

		if (projectionPlaneCheck)
		{
			vvtmp.reserve (vv.size ());
			applyProjOnPlane (normal, vv, vvtmp);
			if (vvtmp.size() >= 6)
				vv = vvtmp;
		}

		if (vv.size()<6)
		{
			std::cerr << "Could not compute curvature of radius " << scaledRadius << endl;
			return;
		}
		if (montecarlo)
		{
			if(montecarloN<6)
				break;
			vvtmp.reserve(vv.size());
			applyMontecarlo(vv,&vvtmp);
			vv=vvtmp;
		}

		if (vv.size()<6)
			return;
		computeReferenceFrame(i,normal,ref);

		Quadric q;
		fitQuadric (me, ref, vv, &q);
		finalEigenStuff(i,ref,q);
	}

	lastRadius=sphereRadius;
	curvatureComputed=true;
}

inline void CurvatureCalculator::printCurvature(std::string outpath)
{
	using namespace std;
	if (!curvatureComputed)
		return;

	std::ofstream of;
	of.open(outpath.c_str());

	if (!of)
	{
		fprintf(stderr, "Error: could not open output file %s\n", outpath.c_str());
		return;
	}

	size_t vertices_count=vertices.rows();
	of << vertices_count << endl;
	for (int i=0; i<vertices_count; i++)
	{
		of << curv[i][0] << " " << curv[i][1] << " " << curvDir[i][0][0] << " " << curvDir[i][0][1] << " " << curvDir[i][0][2] << " " <<
			curvDir[i][1][0] << " " << curvDir[i][1][1] << " " << curvDir[i][1][2] << endl;
	}

	of.close();

}

template <typename FaceArray, typename Index>
inline void adjacency_list(const FaceArray & F,std::vector<std::vector<Index> >& A)
{
	A.clear(); 
	A.resize(F.maxCoeff()+1);

	// Loop over faces
	for(int i = 0;i<F.rows();i++)
	{
		// Loop over this face (triangle)
		for(int j = 0;j<3;j++)
		{
			// Get indices of edge: s --> d
			int s = F(i,j);
			int d = F(i,(j+1) % 3);
			A.at(s).push_back(d);
			A.at(d).push_back(s);
		}
	}

	// Remove duplicates
	for(int i=0; i<(int)A.size();++i)
	{
		std::sort(A[i].begin(), A[i].end());
		A[i].erase(std::unique(A[i].begin(), A[i].end()), A[i].end());
	}
}

template <typename DerivedV, typename DerivedF, typename IndexType>
inline void vf( const Eigen::PlainObjectBase<DerivedV>& V,const Eigen::PlainObjectBase<DerivedF>& F,
	std::vector<std::vector<IndexType> >& VF,std::vector<std::vector<IndexType> >& VFi)
{
	VF.clear(); VF.resize(V.rows());
	VFi.clear(); VFi.resize(V.rows());

	for(int fi=0; fi<F.rows(); ++fi){
		for(int i = 0; i < F.cols(); ++i){
			VF[F(fi,i)].push_back(fi);
			VFi[F(fi,i)].push_back(i);
		}
	}
}

template <typename DerivedV, typename DerivedF>
inline void principal_curvature(
	const Eigen::PlainObjectBase<DerivedV>& V,
	const Eigen::PlainObjectBase<DerivedF>& F,
	const Eigen::PlainObjectBase<DerivedV>& NV,
	const Eigen::PlainObjectBase<DerivedV>& NF,
	Eigen::PlainObjectBase<DerivedV>& PD1,
	Eigen::PlainObjectBase<DerivedV>& PD2,
	unsigned radius, double * maxValue = NULL)
{
	using namespace std;

	// Preallocate memory
	PD1.resize(V.rows(),3);
	PD2.resize(V.rows(),3);

	CurvatureCalculator cc;
	cc.vertices = V;
	cc.faces = F;
	cc.vertex_normals = NV;
	cc.face_normals = NF;
	cc.sphereRadius = radius;

	// Adjacency
	adjacency_list(F, cc.vertex_to_vertices);
	vf(V, F, cc.vertex_to_faces, cc.vertex_to_faces_index);

	// Compute
	cc.computeCurvature();

	if(maxValue) *maxValue = 0.0;

	// Copy it back
	for (unsigned i=0; i<V.rows(); i++)
	{
		Eigen::Vector3d d1;
		Eigen::Vector3d d2;
		d1 << cc.curvDir[i][0][0], cc.curvDir[i][0][1], cc.curvDir[i][0][2];
		d2 << cc.curvDir[i][1][0], cc.curvDir[i][1][1], cc.curvDir[i][1][2];
		d1.normalize();
		d2.normalize();
		d1 *= cc.curv[i][0];
		d2 *= cc.curv[i][1];
		PD1.row(i) = d1;
		PD2.row(i) = d2;

		// Keep track of maximum
		if(maxValue)
		{
			if(cc.curv[i][0] > *maxValue) *maxValue = cc.curv[i][0];
			if(cc.curv[i][1] > *maxValue) *maxValue = cc.curv[i][1];
		}

		if (PD1.row(i) * PD2.row(i).transpose() > 10e-6)
		{
			cerr << "Something is wrong, vertex: i" << endl;
			PD1.row(i) *= 0;
			PD2.row(i) *= 0;
		}
	}
}
}