vcglib/vcg/complex/trimesh/update/curvature.h

/****************************************************************************
* VCGLib                                                            o o     *
* Visual and Computer Graphics Library                            o     o   *
*                                                                _   O  _   *
* Copyright(C) 2004                                                \/)\/    *
* Visual Computing Lab                                            /\/|      *
* ISTI - Italian National Research Council                           |      *
*                                                                    \      *
* All rights reserved.                                                      *
*                                                                           *
* This program is free software; you can redistribute it and/or modify      *   
* it under the terms of the GNU General Public License as published by      *
* the Free Software Foundation; either version 2 of the License, or         *
* (at your option) any later version.                                       *
*                                                                           *
* This program 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 General Public License (http://www.gnu.org/licenses/gpl.txt)          *
* for more details.                                                         *
*                                                                           *
****************************************************************************/
/****************************************************************************
  History
$Log: not supported by cvs2svn $
Revision 1.2  2005/10/25 09:17:41  spinelli
correct IsBorder

Revision 1.1  2005/02/22 16:40:29  ganovelli
created. This version writes the gaussian curvature on the Q() member of
the vertex

/****************************************************************************/

#ifndef VCGLIB_UPDATE_CURVATURE_
#define VCGLIB_UPDATE_CURVATURE_

#include <vcg/math/base.h>
#include <vcg/simplex/face/topology.h>
#include <vcg/simplex/face/pos.h>

namespace vcg {
namespace tri {

/** \addtogroup trimesh */
/*@{*/

/// Management, updating and computation of per-vertex and per-face normals.
/// This class is used to compute or update the normals that can be stored in the vertex or face component of a mesh.
template <class ComputeMeshType>
class UpdateCurvature
{

public:
typedef ComputeMeshType MeshType; 	
typedef typename MeshType::VertexType     VertexType;
typedef typename VertexType::NormalType     NormalType;
typedef typename VertexType::ScalarType ScalarType;
typedef typename MeshType::VertexPointer  VertexPointer;
typedef typename MeshType::VertexIterator VertexIterator;
typedef typename MeshType::FaceType       FaceType;
typedef typename MeshType::FacePointer    FacePointer;
typedef typename MeshType::FaceIterator   FaceIterator;


/** computes the discrete gaussian curvature as proposed in 
Discrete Differential-Geometry Operators for Triangulated 2-Manifolds Mark Meyer,
 Mathieu Desbrun, Peter Schroder, Alan H. Barr VisMath '02, Berlin
*/
static void Gaussian( MeshType &  m){
		assert(m.HasPerVertexQuality());

	typename MeshType::VertexIterator vi;   // iteratore vertice
	typename MeshType::FaceIterator fi;     // iteratore facce
	double *area;                  // areamix vector
	int i;												 // index
	double area0, area1, area2;
	double angle0, angle1, angle2; 
	
	//--- Initialization
	area = new double[m.vn];

	//reset the values to 0
	for(vi=m.vert.begin();vi!=m.vert.end();++vi) if(!(*vi).IsD())
		(*vi).Q() = 0.0;

	//--- compute Areamix
	for(fi=m.face.begin();fi!=m.face.end();++fi) if(!(*fi).IsD())
	{
		
		// angles
			 angle0 = math::Abs(Angle(	(*fi).V(1)->P()-(*fi).V(0)->P(),(*fi).V(2)->P()-(*fi).V(0)->P() ));
			 angle1 = math::Abs(Angle(	(*fi).V(0)->P()-(*fi).V(1)->P(),(*fi).V(2)->P()-(*fi).V(1)->P() ));
 			 angle2 = M_PI-(angle0+angle1);
		
		if((angle0 < M_PI/2) || (angle1 < M_PI/2) || (angle2 < M_PI/2))  // triangolo non ottuso
		{ 
			float e01 = SquaredDistance( (*fi).V(1)->P() , (*fi).V(0)->P() );
			float e12 = SquaredDistance( (*fi).V(2)->P() , (*fi).V(1)->P() );
			float e20 = SquaredDistance( (*fi).V(0)->P() , (*fi).V(2)->P() );
			
			// voronoi area v[0]
			area0 = ( e01*(1/tan(angle2)) + e20*(1/tan(angle1)) ) /8;
			// voronoi area v[1]
			area1 = ( e01*(1/tan(angle2)) + e12*(1/tan(angle0)) ) /8;
			// voronoi area v[2]
			area2 = ( e20*(1/tan(angle1)) + e20*(1/tan(angle0)) ) /8;
			
			(*fi).V(0)->Q()  += area0;
			(*fi).V(1)->Q()  += area1;
			(*fi).V(2)->Q()  += area2;
		}
		else // triangolo ottuso
		{  
			(*fi).V(0)->Q() += (*fi).DoubleArea() / 6;
			(*fi).V(1)->Q() += (*fi).DoubleArea() / 6;
			(*fi).V(2)->Q() += (*fi).DoubleArea() / 6;            
		}
	}

	i = 0;
	for(vi=m.vert.begin();vi!=m.vert.end();++vi,++i) if(!(*vi).IsD())
	{
		area[i] = (*vi).Q();
		(*vi).Q() = (float)(2.0 * M_PI);
	}
	
	if(false)
	for(fi=m.face.begin();fi!=m.face.end();++fi) if(!(*fi).IsD())
	{
		for(int i=0;i<3;i++)
		{
			if(vcg::face::IsBorder((*fi), i))
			{
				typename MeshType::CoordType e1,e2;
				vcg::face::Pos<FaceType> hp(&*fi,i,(*fi).V(i));
				//MeshType::hedgepos_type hp(&*fi,i,(*fi).V(i));
				vcg::face::Pos<FaceType> hp1=hp;
				//MeshType::hedgepos_type hp1=hp;

				hp1.FlipV();
	
				e1= hp1.v->P()-hp.v->P();
				hp1.FlipV();
				hp1.NextB();
				e2= hp1.v->P()-hp.v->P();
				(*fi).V(i)->Q() -=math::Abs(Angle(e1,e2));
			}
	  }
	}
	
	
	for(fi=m.face.begin();fi!=m.face.end();++fi)  if(!(*fi).IsD())
	{
		float angle0 = math::Abs(Angle(
			(*fi).V(1)->P()-(*fi).V(0)->P(),(*fi).V(2)->P()-(*fi).V(0)->P() ));
		float angle1 = math::Abs(Angle(
			(*fi).V(0)->P()-(*fi).V(1)->P(),(*fi).V(2)->P()-(*fi).V(1)->P() ));
		float angle2 = M_PI-(angle0+angle1);
		
		(*fi).V(0)->Q() -= angle0;
		(*fi).V(1)->Q() -= angle1;
		(*fi).V(2)->Q() -= angle2;
	}
	i=0;
	for(vi=m.vert.begin(); vi!=m.vert.end(); ++vi,++i) if(!(*vi).IsD())
	{
		(*vi).Q() /= area[i];
		(*vi).Q()=math::Clamp((*vi).Q(),-0.050f,0.050f);
		
 /*   if ( (*vi).Q() < 0 )
			(*vi).Q() = log( -(*vi).Q() );
		else if( (*vi).Q() > 0 )
			(*vi).Q() = log( (*vi).Q() );*/

	}
	
	//--- DeInit
	
	delete[] area;

}
};
}
}
#endif