//==============================================================================
/*! \file
 * OpenMesh Toolkit for mesh analysis    \n
 * Copyright (c) 2010 by Rostislav Hulik     \n
 *
 * Author:  Rostislav Hulik, rosta.hulik@gmail.com  \n
 * Date:    2010/11/03                          \n
 *
 * This file is part of software developed for support of Rostislav Hulik's dissertation thesis at dcgm-robotics@FIT group.
 *
 * This file is free software: you can redistribute it and/or modify
 * it under the terms of the GNU Lesser General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 * 
 * This file is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU Lesser General Public License for more details.
 * 
 * You should have received a copy of the GNU Lesser General Public License
 * along with this file.  If not, see <http://www.gnu.org/licenses/>.
 * 
 * Description:
 * - Class computing curvatures on a mesh
 */

#ifndef _OM_CURVATURE_H_
#define _OM_CURVATURE_H_

#include <OpenMesh\Core\Utils\Property.hh>
#include <math.h>
#include <Eigen/QR>

namespace OMToolkit
{
	/**
	 * Class computes curvatures on a mesh
	 * @tparam Mesh Mesh type to work with
	 * @tparam Scalar type of computed curvature
	 */
	template <class Mesh, class Scalar>
	class OMCurvature
	{
		protected:
			/**
			 * Mesh vector type
			 */
			typedef typename Mesh::Normal Normal;

			/**
			 * Mesh vertex handle type
			 */
			typedef typename Mesh::VertexHandle VertexHandle;

			/**
			 * Mesh face handle type
			 */
			typedef typename Mesh::FaceHandle FaceHandle;

			/**
			 * Mesh halfedge handle type
			 */
			typedef typename Mesh::HalfedgeHandle HalfedgeHandle;

			/**
			 * Eigen matrix 3x3
			 */
			typedef Eigen::Matrix<Scalar, 3, 3> Matrix3x3;

			/**
			 * Eigen matrix 2x2
			 */
			typedef Eigen::Matrix<Scalar, 2, 2> Matrix2x2;

			/**
			 * Eigen row vector 3D
			 */
			typedef Eigen::Matrix<Scalar, 1, 3> RowVector3;

			/**
			 * Eigen row vector 2D
			 */
			typedef Eigen::Matrix<Scalar, 1, 2> RowVector2;

			/**
			 * Eigen column vector 3D
			 */
			typedef Eigen::Matrix<Scalar, 3, 1> ColumnVector3;

			/**
			 * Eigen column vector 2D
			 */
			typedef Eigen::Matrix<Scalar, 3, 1> ColumnVector2;

		public:
			
			/**
			 * Constructor - create an instance of this class 
			 * @param mesh Mesh to work with
			 * @param curvVectorHandle Handle to a curvature direction vertex property
			 * @param curvMagHandle Handle to a curvature magnitude vertex property
			 */
			OMCurvature(Mesh *mesh, OpenMesh::VPropHandleT<Normal> &curvVectorHandle, OpenMesh::VPropHandleT<Scalar> &curvMagHandle);

			/**
			 * Calculates mean curvature and its direction for all mesh vertices
			 * @return True, if calculated successfully
			 */
			bool calcMeanCurvature();

			/**
			 * Calculates Gauss curvature and for all mesh vertices
			 * @return True, if calculated successfully
			 */
			bool calcGaussCurvature();

			/**
			 * Calculates minimal curvature and for all mesh vertices
			 * @return True, if calculated successfully
			 * @todo Need to implement directions
			 */
			bool calcMinCurvature();

			/**
			 * Calculates maximal curvature and for all mesh vertices
			 * @return True, if calculated successfully
			 * @todo Need to implement directions
			 */
			bool calcMaxCurvature();
			
		protected:
			/**
			 * Function returns principal curvature directions for specified vertex
			 * Computation from paper of Gabriel Taubin - Estimating the tenspr of curvature of a surface from a polyhedral approximation
			 * @param vertex Vertex handle
			 * @param max maximum curvature direction vector
			 * @param min minimum curvature direction vector
			 */
			void getCurvatureDirections(VertexHandle vertex, Normal &max, Normal &min);

			/**
			 * Returns area of voronoi region aroun specified vertex
			 * @param vertex Vertex handle
			 * @return area of voronoi region
			 */
			Scalar getVoronoiArea(VertexHandle vertex);

			/**
			 * Computes mean curvature operator (vector, which direction equals mean curvature direction, its norm equals 2*mean curvature) for specified vertex
			 * @param vertex Vertex handle
			 * @return mean curvature operator
			 */
			Normal getMeanOperator(VertexHandle vertex);

			/**
			 * Computes Gauss curvature value for specified vertex
			 * @param vertex Vertex handle
			 * @return Gauss curvature magnitude
			 */
			Scalar getGaussOperator(VertexHandle vertex);
			
			/**
			 * Computes voronoi area of non-obtuse triangle from given vertex 
			 * @param face Given triangle
			 * @param face Vertex handle
			 * @return Area of voronoi region
			 */
			Scalar getVoronoiAreaTri(FaceHandle face, VertexHandle vertex);

			/**
			 * Computes voronoi area of obtuse triangle from given vertex 
			 * @param face Given triangle
			 * @param face Vertex handle
			 * @return Area of voronoi region
			 */
			Scalar getVoronoiAreaTriO(FaceHandle face, VertexHandle vertex);

			/**
			 * Predicate, is triangle obtuse? (Is one of its angles obtuse?)
			 * @param face Specified triangle
			 * @param half_edge Half edge handle pointing on found obtuse angle (if not obtuse, it is not specified)
			 * @return true, if triangle has an obtuse angle
			 */
			bool isObtuse(FaceHandle face, HalfedgeHandle &half_edge);

			/**
			 * Intern pointer to a mesh
			 */
			Mesh *m_mesh;

			/**
			 * Curvature direction property handle
			 */
			OpenMesh::VPropHandleT<Normal> m_curvVectorHandle;

			/**
			 * Curvature magnitude property handle
			 */
			OpenMesh::VPropHandleT<Scalar> m_curvMagHandle;

			/**
			 * Precomputed voronoi area vector for each face
			 */
			typedef OpenMesh::VectorT<Scalar, 3> AreasVector;

			/**
			 * Precomputed voronoi area property handle
			 */
			OpenMesh::FPropHandleT<AreasVector> m_areas;
	}; // class

#include <OMToolkit\OMCurvature.hxx>
} // namespace OMToolkit

#endif