taubin and desbrun estimates added (-> see vcg/simplex/vertexplus/component.h [component_ocf.h|component_occ.h ]
This commit is contained in:
ganovelli
parent
4e7d6a2765
commit
8a265c9d32
|
|
@ -23,6 +23,11 @@
|
|||
/****************************************************************************
|
||||
History
|
||||
$Log: not supported by cvs2svn $
|
||||
Revision 1.3 2006/02/27 18:02:57 ponchio
|
||||
Area -> doublearea/2
|
||||
|
||||
added some typename
|
||||
|
||||
Revision 1.2 2005/10/25 09:17:41 spinelli
|
||||
correct IsBorder
|
||||
|
||||
|
|
@ -36,8 +41,11 @@ the vertex
|
|||
#define VCGLIB_UPDATE_CURVATURE_
|
||||
|
||||
#include <vcg/math/base.h>
|
||||
#include <vcg/math/matrix.h>
|
||||
#include <vcg/simplex/face/topology.h>
|
||||
#include <vcg/simplex/face/pos.h>
|
||||
#include <vcg/simplex/face/jumping_pos.h>
|
||||
#include <vcg/container/simple_temporary_data.h>
|
||||
|
||||
namespace vcg {
|
||||
namespace tri {
|
||||
|
|
@ -47,140 +55,330 @@ namespace tri {
|
|||
|
||||
/// 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>
|
||||
template <class MeshType>
|
||||
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;
|
||||
typedef typename MeshType::FaceIterator FaceIterator;
|
||||
typedef typename MeshType::VertexIterator VertexIterator;
|
||||
typedef typename MeshType::VertContainer VertContainer;
|
||||
typedef typename MeshType::VertexType VertexType;
|
||||
typedef typename MeshType::FaceType FaceType;
|
||||
typedef typename MeshType::CoordType CoordType;
|
||||
typedef typename CoordType::ScalarType ScalarType;
|
||||
|
||||
|
||||
private:
|
||||
typedef struct AdjVertex {
|
||||
VertexType * vert;
|
||||
float doubleArea;
|
||||
bool isBorder;
|
||||
};
|
||||
public:
|
||||
/*
|
||||
Compute principal direction and magniuto of curvature as describe in the paper:
|
||||
@InProceedings{bb33922,
|
||||
author = "G. Taubin",
|
||||
title = "Estimating the Tensor of Curvature of a Surface from a
|
||||
Polyhedral Approximation",
|
||||
booktitle = "International Conference on Computer Vision",
|
||||
year = "1995",
|
||||
pages = "902--907",
|
||||
URL = "http://dx.doi.org/10.1109/ICCV.1995.466840",
|
||||
bibsource = "http://www.visionbib.com/bibliography/describe440.html#TT32253",
|
||||
}
|
||||
*/
|
||||
static void PrincipalDirections(MeshType &m) {
|
||||
|
||||
assert(m.HasVFTopology());
|
||||
|
||||
vcg::tri::UpdateNormals<MeshType>::PerVertexNormalized(m);
|
||||
|
||||
VertexIterator vi;
|
||||
for (vi =m.vert.begin(); vi !=m.vert.end(); ++vi) {
|
||||
if ( ! (*vi).IsD() && (*vi).VFp() != NULL) {
|
||||
|
||||
VertexType * central_vertex = &(*vi);
|
||||
|
||||
std::vector<float> weights;
|
||||
std::vector<AdjVertex> vertices;
|
||||
|
||||
vcg::face::JumpingPos<FaceType> pos((*vi).VFp(), central_vertex);
|
||||
|
||||
VertexType* firstV = pos.VFlip();
|
||||
VertexType* tempV;
|
||||
float totalDoubleAreaSize = 0.0f;
|
||||
|
||||
if (((firstV->P()-central_vertex->P())^(pos.VFlip()->P()-central_vertex->P()))*central_vertex->N()<=0.0f)
|
||||
{
|
||||
pos.Set(central_vertex->VFp(), central_vertex);
|
||||
pos.FlipE();
|
||||
firstV = pos.VFlip();
|
||||
}
|
||||
else pos.Set(central_vertex->VFp(), central_vertex);
|
||||
|
||||
do
|
||||
{
|
||||
pos.NextE();
|
||||
tempV = pos.VFlip();
|
||||
|
||||
AdjVertex v;
|
||||
|
||||
v.isBorder = pos.IsBorder();
|
||||
v.vert = tempV;
|
||||
v.doubleArea = ((pos.F()->V(1)->P() - pos.F()->V(0)->P()) ^ (pos.F()->V(2)->P()- pos.F()->V(0)->P())).Norm();;
|
||||
totalDoubleAreaSize += v.doubleArea;
|
||||
|
||||
vertices.push_back(v);
|
||||
}
|
||||
while(tempV != firstV);
|
||||
|
||||
for (int i = 0; i < vertices.size(); ++i) {
|
||||
if (vertices[i].isBorder) {
|
||||
weights.push_back(vertices[i].doubleArea / totalDoubleAreaSize);
|
||||
} else {
|
||||
weights.push_back(0.5f * (vertices[i].doubleArea + vertices[(i-1)%vertices.size()].doubleArea) / totalDoubleAreaSize);
|
||||
}
|
||||
assert(weights.back() < 1.0f);
|
||||
}
|
||||
|
||||
Matrix33f Tp;
|
||||
for (int i = 0; i < 3; ++i)
|
||||
Tp[i][i] = 1.0f - powf(central_vertex->N()[i],2);
|
||||
Tp[0][1] = Tp[1][0] = -1.0f * (central_vertex->N()[0] * central_vertex->N()[1]);
|
||||
Tp[1][2] = Tp[2][1] = -1.0f * (central_vertex->N()[1] * central_vertex->N()[2]);
|
||||
Tp[0][2] = Tp[2][0] = -1.0f * (central_vertex->N()[0] * central_vertex->N()[2]);
|
||||
|
||||
Matrix33f tempMatrix;
|
||||
Matrix33f M;
|
||||
M.SetZero();
|
||||
for (int i = 0; i < vertices.size(); ++i) {
|
||||
Point3f edge = (central_vertex->P() - vertices[i].vert->P());
|
||||
float curvature = (2.0f * (central_vertex->N() * edge) ) / edge.SquaredNorm();
|
||||
Point3f T = (Tp*edge).Normalize();
|
||||
tempMatrix.ExternalProduct(T,T);
|
||||
M += tempMatrix * weights[i] * curvature ;
|
||||
}
|
||||
|
||||
Point3f W;
|
||||
Point3f e1(1.0f,0.0f,0.0f);
|
||||
if ((e1 - central_vertex->N()).SquaredNorm() > (e1 + central_vertex->N()).SquaredNorm())
|
||||
W = e1 - central_vertex->N();
|
||||
else
|
||||
W = e1 + central_vertex->N();
|
||||
W.Normalize();
|
||||
|
||||
Matrix33f Q;
|
||||
Q.SetIdentity();
|
||||
tempMatrix.ExternalProduct(W,W);
|
||||
Q -= tempMatrix * 2.0f;
|
||||
|
||||
Matrix33f Qt(Q);
|
||||
Qt.Transpose();
|
||||
|
||||
Matrix33f QtMQ = (Qt * M * Q);
|
||||
|
||||
Point3f bl = Q.GetColumn(0);
|
||||
Point3f T1 = Q.GetColumn(1);
|
||||
Point3f T2 = Q.GetColumn(2);
|
||||
|
||||
float s,c;
|
||||
// Gabriel Taubin hint and Valentino Fiorin impementation
|
||||
float qt21 = QtMQ[2][1];
|
||||
float qt12 = QtMQ[1][2];
|
||||
float alpha = QtMQ[1][1]-QtMQ[2][2];
|
||||
float beta = QtMQ[2][1];
|
||||
|
||||
float h[2];
|
||||
float delta = sqrtf(4.0f*powf(alpha, 2) +16.0f*powf(beta, 2));
|
||||
h[0] = (2.0f*alpha + delta) / (2.0f*beta);
|
||||
h[1] = (2.0f*alpha - delta) / (2.0f*beta);
|
||||
|
||||
float t[2];
|
||||
float best_c, best_s;
|
||||
float min_error = std::numeric_limits<float>::infinity();
|
||||
for (int i=0; i<2; i++)
|
||||
{
|
||||
delta = sqrtf(powf(h[1], 2) + 4.0f);
|
||||
t[0] = (h[i]+delta) / 2.0f;
|
||||
t[1] = (h[i]-delta) / 2.0f;
|
||||
|
||||
for (int j=0; j<2; j++)
|
||||
{
|
||||
float squared_t = powf(t[j], 2);
|
||||
float denominator = 1.0f + squared_t;
|
||||
s = (2.0f*t[j]) / denominator;
|
||||
c = (1-squared_t) / denominator;
|
||||
|
||||
float approximation = c*s*alpha + (powf(c, 2) - powf(s, 2))*beta;
|
||||
float angle_similarity = fabs(acosf(c)/asinf(s));
|
||||
float error = fabs(1.0f-angle_similarity)+fabs(approximation);
|
||||
if (error<min_error)
|
||||
{
|
||||
min_error = error;
|
||||
best_c = c;
|
||||
best_s = s;
|
||||
}
|
||||
}
|
||||
}
|
||||
c = best_c;
|
||||
s = best_s;
|
||||
|
||||
vcg::ndim::MatrixMNf minor2x2 (2,2);
|
||||
vcg::ndim::MatrixMNf S (2,2);
|
||||
|
||||
|
||||
/** computes the discrete gaussian curvature as proposed in
|
||||
minor2x2[0][0] = QtMQ[1][1];
|
||||
minor2x2[0][1] = QtMQ[1][2];
|
||||
minor2x2[1][0] = QtMQ[2][1];
|
||||
minor2x2[1][1] = QtMQ[2][2];
|
||||
|
||||
S[0][0] = S[1][1] = c;
|
||||
S[0][1] = s;
|
||||
S[1][0] = -1.0f * s;
|
||||
|
||||
vcg::ndim::MatrixMNf St (S);
|
||||
St.Transpose();
|
||||
|
||||
vcg::ndim::MatrixMNf StMS(St * minor2x2 * S);
|
||||
|
||||
float Principal_Curvature1 = (3.0f * StMS[0][0]) - StMS[1][1];
|
||||
float Principal_Curvature2 = (3.0f * StMS[1][1]) - StMS[0][0];
|
||||
|
||||
Point3f Principal_Direction1 = T1 * c - T2 * s;
|
||||
Point3f Principal_Direction2 = T1 * s + T2 * c;
|
||||
|
||||
(*vi).PD1() = Principal_Direction1 ;
|
||||
(*vi).PD2() = Principal_Direction2 ;
|
||||
(*vi).K1() = -Principal_Curvature1;
|
||||
(*vi).K2() = -Principal_Curvature2;
|
||||
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
|
||||
class AreaData
|
||||
{
|
||||
public:
|
||||
float A;
|
||||
};
|
||||
|
||||
|
||||
/** 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());
|
||||
*/
|
||||
static void MeanAndGaussian(MeshType & m)
|
||||
{
|
||||
float area0, area1, area2, angle0, angle1, angle2, e01, e12, e20;
|
||||
FaceIterator fi;
|
||||
VertexIterator vi;
|
||||
|
||||
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];
|
||||
SimpleTempData<VertContainer, AreaData> TDAreaPtr(m.vert); TDAreaPtr.Start();
|
||||
|
||||
//reset the values to 0
|
||||
for(vi=m.vert.begin();vi!=m.vert.end();++vi) if(!(*vi).IsD())
|
||||
(*vi).Q() = 0.0;
|
||||
//Calcola AreaMix in H (vale anche per K)
|
||||
for(vi=m.vert.begin(); vi!=m.vert.end(); ++vi) if(!(*vi).IsD())
|
||||
{
|
||||
(TDAreaPtr)[*vi].A = 0;
|
||||
(*vi).H() = 0;
|
||||
(*vi).K() = (float)(2.0 * M_PI);
|
||||
}
|
||||
|
||||
//--- 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;
|
||||
}
|
||||
}
|
||||
for(fi=m.face.begin();fi!=m.face.end();++fi) if( !(*fi).IsD())
|
||||
{
|
||||
// angles
|
||||
angle0 = math::Abs(Angle( (*fi).P(1)-(*fi).P(0),(*fi).P(2)-(*fi).P(0) ));
|
||||
angle1 = math::Abs(Angle( (*fi).P(0)-(*fi).P(1),(*fi).P(2)-(*fi).P(1) ));
|
||||
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() );
|
||||
|
||||
area0 = ( e20*(1.0/tan(angle1)) + e01*(1.0/tan(angle2)) ) / 8.0;
|
||||
area1 = ( e01*(1.0/tan(angle2)) + e12*(1.0/tan(angle0)) ) / 8.0;
|
||||
area2 = ( e12*(1.0/tan(angle0)) + e20*(1.0/tan(angle1)) ) / 8.0;
|
||||
|
||||
(TDAreaPtr)[(*fi).V(0)].A += area0;
|
||||
(TDAreaPtr)[(*fi).V(1)].A += area1;
|
||||
(TDAreaPtr)[(*fi).V(2)].A += area2;
|
||||
|
||||
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;
|
||||
}
|
||||
else // triangolo ottuso
|
||||
{
|
||||
(TDAreaPtr)[(*fi).V(0)].A += vcg::DoubleArea<typename MeshType::FaceType>((*fi)) / 6.0;
|
||||
(TDAreaPtr)[(*fi).V(1)].A += vcg::DoubleArea<typename MeshType::FaceType>((*fi)) / 6.0;
|
||||
(TDAreaPtr)[(*fi).V(2)].A += vcg::DoubleArea<typename MeshType::FaceType>((*fi)) / 6.0;
|
||||
}
|
||||
}
|
||||
|
||||
for(fi=m.face.begin();fi!=m.face.end();++fi) if( !(*fi).IsD() )
|
||||
{
|
||||
angle0 = math::Abs(Angle( (*fi).P(1)-(*fi).P(0),(*fi).P(2)-(*fi).P(0) ));
|
||||
angle1 = math::Abs(Angle( (*fi).P(0)-(*fi).P(1),(*fi).P(2)-(*fi).P(1) ));
|
||||
angle2 = M_PI-(angle0+angle1);
|
||||
|
||||
e01 = ( (*fi).V(1)->P() - (*fi).V(0)->P() ) * (*fi).V(0)->N();
|
||||
e12 = ( (*fi).V(2)->P() - (*fi).V(1)->P() ) * (*fi).V(1)->N();
|
||||
e20 = ( (*fi).V(0)->P() - (*fi).V(2)->P() ) * (*fi).V(2)->N();
|
||||
|
||||
area0 = ( e20 * (1.0/tan(angle1)) + e01 * (1.0/tan(angle2)) ) / 4.0;
|
||||
area1 = ( e01 * (1.0/tan(angle2)) + e12 * (1.0/tan(angle0)) ) / 4.0;
|
||||
area2 = ( e12 * (1.0/tan(angle0)) + e20 * (1.0/tan(angle1)) ) / 4.0;
|
||||
|
||||
(*fi).V(0)->H() += area0;
|
||||
(*fi).V(1)->H() += area1;
|
||||
(*fi).V(2)->H() += area2;
|
||||
|
||||
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));
|
||||
(*fi).V(0)->K() -= angle0;
|
||||
(*fi).V(1)->K() -= angle1;
|
||||
(*fi).V(2)->K() -= angle2;
|
||||
|
||||
|
||||
for(int i=0;i<3;i++)
|
||||
{
|
||||
if(vcg::face::IsBorder((*fi), i))
|
||||
{
|
||||
CoordType e1,e2;
|
||||
vcg::face::Pos<FaceType> hp(&*fi, i, (*fi).V(i));
|
||||
vcg::face::Pos<FaceType> 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)->K() -= math::Abs(Angle(e1,e2));
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
for(vi=m.vert.begin(); vi!=m.vert.end(); ++vi) if(!(*vi).IsD() /*&& !(*vi).IsB()*/)
|
||||
{
|
||||
if((TDAreaPtr)[*vi].A<=std::numeric_limits<float>::epsilon())
|
||||
{
|
||||
(*vi).H() = 0;
|
||||
(*vi).K() = 0;
|
||||
}
|
||||
else
|
||||
{
|
||||
(*vi).H() /= (TDAreaPtr)[*vi].A;
|
||||
(*vi).K() /= (TDAreaPtr)[*vi].A;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
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;
|
||||
TDAreaPtr.Stop();
|
||||
|
||||
}
|
||||
}
|
||||
|
||||
};
|
||||
|
||||
|
||||
}
|
||||
}
|
||||
#endif
|
||||
|
|
|
|||
Loading…
Reference in New Issue