Ongoing Rearrangement of filepath

delete old trimesh content
This commit is contained in:
ganovelli 2011-04-01 17:18:15 +00:00
parent 6c0c32ecfe
commit 25b5f39dad
14 changed files with 0 additions and 6842 deletions

View File

@ -1,94 +0,0 @@
/****************************************************************************
* 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. *
* *
****************************************************************************/
// marco: removed types FaceType, FacePointer, FaceIterator to allow the use of this method from vertex meshes
/****************************************************************************
History
$Log: not supported by cvs2svn $
Revision 1.2 2004/09/15 11:16:27 ganovelli
changed P() to cP()
Revision 1.1 2004/04/05 11:56:13 cignoni
First working version!
Revision 1.2 2004/03/12 15:22:19 cignoni
Written some documentation and added to the trimes doxygen module
Revision 1.1 2004/03/05 10:59:24 cignoni
Changed name from plural to singular (normals->normal)
Revision 1.1 2004/03/04 00:05:50 cignoni
First working version!
Revision 1.1 2004/02/19 13:11:06 cignoni
Initial commit
****************************************************************************/
#ifndef __VCG_TRI_UPDATE_BOUNDING
#define __VCG_TRI_UPDATE_BOUNDING
namespace vcg {
namespace tri {
/// \ingroup trimesh
/// \headerfile bounding.h vcg/complex/trimesh/update/bounding.h
/// \brief 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 UpdateBounding
{
public:
typedef ComputeMeshType MeshType;
typedef typename MeshType::VertexType VertexType;
typedef typename MeshType::VertexPointer VertexPointer;
typedef typename MeshType::VertexIterator VertexIterator;
/// \brief Calculates the bounding box of the \code <ComputeMeshType> \endcode m
static void Box(ComputeMeshType &m)
{
m.bbox.SetNull();
VertexIterator vi;
for(vi = m.vert.begin(); vi != m.vert.end(); ++vi)
if( !(*vi).IsD() ) m.bbox.Add((*vi).cP());
}
}; // end class
} // End namespace
} // End namespace
#endif

View File

@ -1,688 +0,0 @@
/****************************************************************************
* 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.8 2008/05/14 10:03:29 ganovelli
Point3f->Coordtype
Revision 1.7 2008/04/23 16:37:15 onnis
VertexCurvature method added.
Revision 1.6 2008/04/04 10:26:12 cignoni
Cleaned up names, now Kg() gives back Gaussian Curvature (k1*k2), while Kh() gives back Mean Curvature 1/2(k1+k2)
Revision 1.5 2008/03/25 11:00:56 ganovelli
fixed bugs sign of principal direction and mean curvature value
Revision 1.4 2008/03/17 11:29:59 ganovelli
taubin and desbrun estimates added (-> see vcg/simplex/vertex/component.h [component_ocf.h|component_occ.h ]
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
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/space/index/grid_static_ptr.h>
#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>
#include <vcg/complex/trimesh/update/normal.h>
#include <vcg/complex/trimesh/point_sampling.h>
#include <vcg/complex/trimesh/append.h>
#include <vcg/complex/intersection.h>
#include <vcg/complex/trimesh/inertia.h>
#include <vcg/math/matrix33.h>
namespace vcg {
namespace tri {
/// \ingroup trimesh
/// \headerfile curvature.h vcg/complex/trimesh/update/curvature.h
/// \brief 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 MeshType>
class UpdateCurvature
{
public:
typedef typename MeshType::FaceType FaceType;
typedef typename MeshType::FacePointer FacePointer;
typedef typename MeshType::FaceIterator FaceIterator;
typedef typename MeshType::VertexIterator VertexIterator;
typedef typename MeshType::VertContainer VertContainer;
typedef typename MeshType::VertexType VertexType;
typedef typename MeshType::VertexPointer VertexPointer;
typedef vcg::face::VFIterator<FaceType> VFIteratorType;
typedef typename MeshType::CoordType CoordType;
typedef typename CoordType::ScalarType ScalarType;
private:
struct AdjVertex {
VertexType * vert;
float doubleArea;
bool isBorder;
};
public:
/// \brief Compute principal direction and magniuto of curvature.
/*
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);
// firstV is the first vertex of the 1ring neighboorhood
VertexType* firstV = pos.VFlip();
VertexType* tempV;
float totalDoubleAreaSize = 0.0f;
// compute the area of each triangle around the central vertex as well as their total area
do
{
// this bring the pos to the next triangle counterclock-wise
pos.FlipF();
pos.FlipE();
// tempV takes the next vertex in the 1ring neighborhood
tempV = pos.VFlip();
assert(tempV!=central_vertex);
AdjVertex v;
v.isBorder = pos.IsBorder();
v.vert = tempV;
v.doubleArea = vcg::DoubleArea(*pos.F());
totalDoubleAreaSize += v.doubleArea;
vertices.push_back(v);
}
while(tempV != firstV);
// compute the weights for the formula computing matrix M
for (size_t 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);
}
// compute I-NN^t to be used for computing the T_i's
Matrix33<ScalarType> Tp;
for (int i = 0; i < 3; ++i)
Tp[i][i] = 1.0f - powf(central_vertex->cN()[i],2);
Tp[0][1] = Tp[1][0] = -1.0f * (central_vertex->N()[0] * central_vertex->cN()[1]);
Tp[1][2] = Tp[2][1] = -1.0f * (central_vertex->cN()[1] * central_vertex->cN()[2]);
Tp[0][2] = Tp[2][0] = -1.0f * (central_vertex->cN()[0] * central_vertex->cN()[2]);
// for all neighbors vi compute the directional curvatures k_i and the T_i
// compute M by summing all w_i k_i T_i T_i^t
Matrix33<ScalarType> tempMatrix;
Matrix33<ScalarType> M;
M.SetZero();
for (size_t i = 0; i < vertices.size(); ++i) {
CoordType edge = (central_vertex->cP() - vertices[i].vert->cP());
float curvature = (2.0f * (central_vertex->cN().dot(edge)) ) / edge.SquaredNorm();
CoordType T = (Tp*edge).normalized();
tempMatrix.ExternalProduct(T,T);
M += tempMatrix * weights[i] * curvature ;
}
// compute vector W for the Householder matrix
CoordType W;
CoordType e1(1.0f,0.0f,0.0f);
if ((e1 - central_vertex->cN()).SquaredNorm() > (e1 + central_vertex->cN()).SquaredNorm())
W = e1 - central_vertex->cN();
else
W = e1 + central_vertex->cN();
W.Normalize();
// compute the Householder matrix I - 2WW^t
Matrix33<ScalarType> Q;
Q.SetIdentity();
tempMatrix.ExternalProduct(W,W);
Q -= tempMatrix * 2.0f;
// compute matrix Q^t M Q
Matrix33<ScalarType> QtMQ = (Q.transpose() * M * Q);
CoordType bl = Q.GetColumn(0);
CoordType T1 = Q.GetColumn(1);
CoordType T2 = Q.GetColumn(2);
// find sin and cos for the Givens rotation
float s,c;
// Gabriel Taubin hint and Valentino Fiorin impementation
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<ScalarType>::infinity();
for (int i=0; i<2; i++)
{
delta = sqrtf(powf(h[i], 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);
// diagonalize M
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 StMS(S.transpose() * minor2x2 * S);
// compute curvatures and curvature directions
float Principal_Curvature1 = (3.0f * StMS[0][0]) - StMS[1][1];
float Principal_Curvature2 = (3.0f * StMS[1][1]) - StMS[0][0];
CoordType Principal_Direction1 = T1 * c - T2 * s;
CoordType 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;
};
/** Curvature meseaure as described in the paper:
Robust principal curvatures on Multiple Scales, Yong-Liang Yang, Yu-Kun Lai, Shi-Min Hu Helmut Pottmann
SGP 2004
If pointVSfaceInt==true the covariance is computed by montecarlo sampling on the mesh (faster)
If pointVSfaceInt==false the covariance is computed by (analytic)integration over the surface (slower)
*/
typedef vcg::GridStaticPtr <FaceType, ScalarType > MeshGridType;
typedef vcg::GridStaticPtr <VertexType, ScalarType > PointsGridType;
static void PrincipalDirectionsPCA(MeshType &m, ScalarType r, bool pointVSfaceInt = true,vcg::CallBackPos * cb = NULL) {
std::vector<VertexType*> closests;
std::vector<ScalarType> distances;
std::vector<CoordType> points;
VertexIterator vi;
ScalarType area;
MeshType tmpM;
typename std::vector<CoordType>::iterator ii;
vcg::tri::TrivialSampler<MeshType> vs;
MeshGridType mGrid;
PointsGridType pGrid;
// Fill the grid used
if(pointVSfaceInt){
area = Stat<MeshType>::ComputeMeshArea(m);
vcg::tri::SurfaceSampling<MeshType,vcg::tri::TrivialSampler<MeshType> >::Montecarlo(m,vs,1000 * area / (2*M_PI*r*r ));
vi = vcg::tri::Allocator<MeshType>::AddVertices(tmpM,m.vert.size());
for(size_t y = 0; y < m.vert.size(); ++y,++vi) (*vi).P() = m.vert[y].P();
pGrid.Set(tmpM.vert.begin(),tmpM.vert.end());
} else{ mGrid.Set(m.face.begin(),m.face.end()); }
int jj = 0;
for(vi = m.vert.begin(); vi != m.vert.end(); ++vi){
vcg::Matrix33<ScalarType> A,eigenvectors;
vcg::Point3<ScalarType> bp,eigenvalues;
int nrot;
// sample the neighborhood
if(pointVSfaceInt)
{
vcg::tri::GetInSphereVertex<
MeshType,
PointsGridType,std::vector<VertexType*>,
std::vector<ScalarType>,
std::vector<CoordType> >(tmpM,pGrid, (*vi).cP(),r ,closests,distances,points);
A.Covariance(points,bp);
A*=area*area/1000;
}
else{
IntersectionBallMesh<MeshType,ScalarType>( m ,vcg::Sphere3<ScalarType>((*vi).cP(),r),tmpM );
vcg::Point3<ScalarType> _bary;
vcg::tri::Inertia<MeshType>::Covariance(tmpM,_bary,A);
}
Jacobi(A, eigenvalues , eigenvectors, nrot);
// get the estimate of curvatures from eigenvalues and eigenvectors
// find the 2 most tangent eigenvectors (by finding the one closest to the normal)
int best = 0; ScalarType bestv = fabs( (*vi).cN().dot(eigenvectors.GetColumn(0).normalized()) );
for(int i = 1 ; i < 3; ++i){
ScalarType prod = fabs((*vi).cN().dot(eigenvectors.GetColumn(i).normalized()));
if( prod > bestv){bestv = prod; best = i;}
}
(*vi).PD1() = eigenvectors.GetColumn( (best+1)%3).normalized();
(*vi).PD2() = eigenvectors.GetColumn( (best+2)%3).normalized();
// project them to the plane identified by the normal
vcg::Matrix33<ScalarType> rot;
ScalarType angle = acos((*vi).PD1().dot((*vi).N()));
rot.SetRotateRad( - (M_PI*0.5 - angle),(*vi).PD1()^(*vi).N());
(*vi).PD1() = rot*(*vi).PD1();
angle = acos((*vi).PD2().dot((*vi).N()));
rot.SetRotateRad( - (M_PI*0.5 - angle),(*vi).PD2()^(*vi).N());
(*vi).PD2() = rot*(*vi).PD2();
// copmutes the curvature values
const ScalarType r5 = r*r*r*r*r;
const ScalarType r6 = r*r5;
(*vi).K1() = (2.0/5.0) * (4.0*M_PI*r5 + 15*eigenvalues[(best+2)%3]-45.0*eigenvalues[(best+1)%3])/(M_PI*r6);
(*vi).K2() = (2.0/5.0) * (4.0*M_PI*r5 + 15*eigenvalues[(best+1)%3]-45.0*eigenvalues[(best+2)%3])/(M_PI*r6);
if((*vi).K1() < (*vi).K2()) { std::swap((*vi).K1(),(*vi).K2());
std::swap((*vi).PD1(),(*vi).PD2());
if (cb)
{
(*cb)(int(100.0f * (float)jj / (float)m.vn),"Vertices Analysis");
++jj;
} }
}
}
/// \brief Computes the discrete gaussian curvature.
/** For further details, please, refer to: \n
- <em> Discrete Differential-Geometry Operators for Triangulated 2-Manifolds Mark Meyer,
Mathieu Desbrun, Peter Schroder, Alan H. Barr VisMath '02, Berlin </em>
*/
static void MeanAndGaussian(MeshType & m)
{
assert(HasFFAdjacency(m));
float area0, area1, area2, angle0, angle1, angle2;
FaceIterator fi;
VertexIterator vi;
typename MeshType::CoordType e01v ,e12v ,e20v;
SimpleTempData<VertContainer, AreaData> TDAreaPtr(m.vert);
SimpleTempData<VertContainer, typename MeshType::CoordType> TDContr(m.vert);
vcg::tri::UpdateNormals<MeshType>::PerVertexNormalized(m);
//Compute AreaMix in H (vale anche per K)
for(vi=m.vert.begin(); vi!=m.vert.end(); ++vi) if(!(*vi).IsD())
{
(TDAreaPtr)[*vi].A = 0.0;
(TDContr)[*vi] =typename MeshType::CoordType(0.0,0.0,0.0);
(*vi).Kh() = 0.0;
(*vi).Kg() = (float)(2.0 * M_PI);
}
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)->cP() , (*fi).V(0)->cP() );
float e12 = SquaredDistance( (*fi).V(2)->cP() , (*fi).V(1)->cP() );
float e20 = SquaredDistance( (*fi).V(0)->cP() , (*fi).V(2)->cP() );
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;
}
else // obtuse
{
if(angle0 >= M_PI/2)
{
(TDAreaPtr)[(*fi).V(0)].A += vcg::DoubleArea<typename MeshType::FaceType>((*fi)) / 4.0;
(TDAreaPtr)[(*fi).V(1)].A += vcg::DoubleArea<typename MeshType::FaceType>((*fi)) / 8.0;
(TDAreaPtr)[(*fi).V(2)].A += vcg::DoubleArea<typename MeshType::FaceType>((*fi)) / 8.0;
}
else if(angle1 >= M_PI/2)
{
(TDAreaPtr)[(*fi).V(0)].A += vcg::DoubleArea<typename MeshType::FaceType>((*fi)) / 8.0;
(TDAreaPtr)[(*fi).V(1)].A += vcg::DoubleArea<typename MeshType::FaceType>((*fi)) / 4.0;
(TDAreaPtr)[(*fi).V(2)].A += vcg::DoubleArea<typename MeshType::FaceType>((*fi)) / 8.0;
}
else
{
(TDAreaPtr)[(*fi).V(0)].A += vcg::DoubleArea<typename MeshType::FaceType>((*fi)) / 8.0;
(TDAreaPtr)[(*fi).V(1)].A += vcg::DoubleArea<typename MeshType::FaceType>((*fi)) / 8.0;
(TDAreaPtr)[(*fi).V(2)].A += vcg::DoubleArea<typename MeshType::FaceType>((*fi)) / 4.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);
// Skip degenerate triangles.
if(angle0==0 || angle1==0 || angle1==0) continue;
e01v = ( (*fi).V(1)->cP() - (*fi).V(0)->cP() ) ;
e12v = ( (*fi).V(2)->cP() - (*fi).V(1)->cP() ) ;
e20v = ( (*fi).V(0)->cP() - (*fi).V(2)->cP() ) ;
TDContr[(*fi).V(0)] += ( e20v * (1.0/tan(angle1)) - e01v * (1.0/tan(angle2)) ) / 4.0;
TDContr[(*fi).V(1)] += ( e01v * (1.0/tan(angle2)) - e12v * (1.0/tan(angle0)) ) / 4.0;
TDContr[(*fi).V(2)] += ( e12v * (1.0/tan(angle0)) - e20v * (1.0/tan(angle1)) ) / 4.0;
(*fi).V(0)->Kg() -= angle0;
(*fi).V(1)->Kg() -= angle1;
(*fi).V(2)->Kg() -= 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->cP() - hp.v->cP();
hp1.FlipV();
hp1.NextB();
e2=hp1.v->cP() - hp.v->cP();
(*fi).V(i)->Kg() -= 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<ScalarType>::epsilon())
{
(*vi).Kh() = 0;
(*vi).Kg() = 0;
}
else
{
(*vi).Kh() = (((TDContr)[*vi].dot((*vi).cN())>0)?1.0:-1.0)*((TDContr)[*vi] / (TDAreaPtr) [*vi].A).Norm();
(*vi).Kg() /= (TDAreaPtr)[*vi].A;
}
}
}
/// \brief Update the mean and the gaussian curvature of a vertex.
/**
The function uses the VF adiacency to walk around the vertex.
\return It will return the voronoi area around the vertex. If (norm == true) the mean and the gaussian curvature are normalized.
Based on the paper <a href="http://www2.in.tu-clausthal.de/~hormann/papers/Dyn.2001.OTU.pdf"> <em> "Optimizing 3d triangulations using discrete curvature analysis" </em> </a>
*/
static float VertexCurvature(VertexPointer v, bool norm = true)
{
// VFAdjacency required!
assert(FaceType::HasVFAdjacency());
assert(VertexType::HasVFAdjacency());
VFIteratorType vfi(v);
float A = 0;
v->Kh() = 0;
v->Kg() = 2 * M_PI;
while (!vfi.End()) {
if (!vfi.F()->IsD()) {
FacePointer f = vfi.F();
int i = vfi.I();
VertexPointer v0 = f->V0(i), v1 = f->V1(i), v2 = f->V2(i);
float ang0 = math::Abs(Angle(v1->P() - v0->P(), v2->P() - v0->P() ));
float ang1 = math::Abs(Angle(v0->P() - v1->P(), v2->P() - v1->P() ));
float ang2 = M_PI - ang0 - ang1;
float s01 = SquaredDistance(v1->P(), v0->P());
float s02 = SquaredDistance(v2->P(), v0->P());
// voronoi cell of current vertex
if (ang0 >= M_PI/2)
A += (0.5f * DoubleArea(*f) - (s01 * tan(ang1) + s02 * tan(ang2)) / 8.0 );
else if (ang1 >= M_PI/2)
A += (s01 * tan(ang0)) / 8.0;
else if (ang2 >= M_PI/2)
A += (s02 * tan(ang0)) / 8.0;
else // non obctuse triangle
A += ((s02 / tan(ang1)) + (s01 / tan(ang2))) / 8.0;
// gaussian curvature update
v->Kg() -= ang0;
// mean curvature update
ang1 = math::Abs(Angle(f->N(), v1->N()));
ang2 = math::Abs(Angle(f->N(), v2->N()));
v->Kh() += ( (math::Sqrt(s01) / 2.0) * ang1 +
(math::Sqrt(s02) / 2.0) * ang2 );
}
++vfi;
}
v->Kh() /= 4.0f;
if(norm) {
if(A <= std::numeric_limits<float>::epsilon()) {
v->Kh() = 0;
v->Kg() = 0;
}
else {
v->Kh() /= A;
v->Kg() /= A;
}
}
return A;
}
static void VertexCurvature(MeshType & m){
for(VertexIterator vi = m.vert.begin(); vi != m.vert.end(); ++vi)
VertexCurvature(&*vi,false);
}
/*
Compute principal curvature directions and value with normal cycle:
@inproceedings{CohMor03,
author = {Cohen-Steiner, David and Morvan, Jean-Marie },
booktitle = {SCG '03: Proceedings of the nineteenth annual symposium on Computational geometry},
title - {Restricted delaunay triangulations and normal cycle}
year = {2003}
}
*/
static void PrincipalDirectionsNormalCycles(MeshType & m){
assert(VertexType::HasVFAdjacency());
assert(FaceType::HasFFAdjacency());
assert(FaceType::HasFaceNormal());
typename MeshType::VertexIterator vi;
for(vi = m.vert.begin(); vi != m.vert.end(); ++vi)
if(!((*vi).IsD())){
vcg::Matrix33<ScalarType> m33;m33.SetZero();
face::JumpingPos<typename MeshType::FaceType> p((*vi).VFp(),&(*vi));
p.FlipE();
typename MeshType::VertexType * firstv = p.VFlip();
assert(p.F()->V(p.VInd())==&(*vi));
do{
if( p.F() != p.FFlip()){
Point3<ScalarType> normalized_edge = p.F()->V(p.F()->Next(p.VInd()))->cP() - (*vi).P();
ScalarType edge_length = normalized_edge.Norm();
normalized_edge/=edge_length;
Point3<ScalarType> n1 = p.F()->cN();n1.Normalize();
Point3<ScalarType> n2 = p.FFlip()->cN();n2.Normalize();
ScalarType n1n2 = (n1 ^ n2).dot(normalized_edge);
n1n2 = std::max(std::min( ScalarType(1.0),n1n2),ScalarType(-1.0));
ScalarType beta = math::Asin(n1n2);
m33[0][0] += beta*edge_length*normalized_edge[0]*normalized_edge[0];
m33[0][1] += beta*edge_length*normalized_edge[1]*normalized_edge[0];
m33[1][1] += beta*edge_length*normalized_edge[1]*normalized_edge[1];
m33[0][2] += beta*edge_length*normalized_edge[2]*normalized_edge[0];
m33[1][2] += beta*edge_length*normalized_edge[2]*normalized_edge[1];
m33[2][2] += beta*edge_length*normalized_edge[2]*normalized_edge[2];
}
p.NextFE();
}while(firstv != p.VFlip());
if(m33.Determinant()==0.0){ // degenerate case
(*vi).K1() = (*vi).K2() = 0.0; continue;}
m33[1][0] = m33[0][1];
m33[2][0] = m33[0][2];
m33[2][1] = m33[1][2];
Point3<ScalarType> lambda;
Matrix33<ScalarType> vect;
int n_rot;
Jacobi(m33,lambda, vect,n_rot);
vect.transposeInPlace();
ScalarType normal = std::numeric_limits<ScalarType>::min();
int normI = 0;
for(int i = 0; i < 3; ++i)
if( fabs((*vi).N().Normalize().dot(vect.GetRow(i))) > normal )
{
normal= fabs((*vi).N().Normalize().dot(vect.GetRow(i)));
normI = i;
}
int maxI = (normI+2)%3;
int minI = (normI+1)%3;
if(fabs(lambda[maxI]) < fabs(lambda[minI])) std::swap(maxI,minI);
(*vi).PD1() = *(Point3<ScalarType>*)(& vect[maxI][0]);
(*vi).PD2() = *(Point3<ScalarType>*)(& vect[minI][0]);
(*vi).K1() = lambda[maxI];
(*vi).K2() = lambda[minI];
}
}
};
}
}
#endif

View File

@ -1,707 +0,0 @@
/****************************************************************************
* 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. *
* *
****************************************************************************/
/****************************************************************************
****************************************************************************/
#ifndef VCGLIB_UPDATE_CURVATURE_FITTING
#define VCGLIB_UPDATE_CURVATURE_FITTING
#include <vcg/space/index/grid_static_ptr.h>
#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>
#include <vcg/complex/trimesh/update/normal.h>
#include <vcg/complex/trimesh/point_sampling.h>
#include <vcg/complex/trimesh/append.h>
#include <vcg/complex/intersection.h>
#include <vcg/complex/trimesh/inertia.h>
#include <vcg/math/matrix33.h>
#include <vcg/Eigen/Core>
#include <vcg/Eigen/QR>
#include <vcg/Eigen/LU>
#include <vcg/Eigen/SVD>
// GG include
#include <vector>
#include <vcg/complex/trimesh/nring.h>
namespace vcg {
namespace tri {
/// \ingroup trimesh
/// \headerfile curvature_fitting.h vcg/complex/trimesh/update/curvature_fitting.h
/// \brief Computation of per-vertex directions and values of curvature.
/**
This class is used to compute the per-vertex directions and values of curvature using a quadric fitting method.
*/
template <class MeshType>
class UpdateCurvatureFitting
{
public:
typedef typename MeshType::FaceType FaceType;
typedef typename MeshType::FacePointer FacePointer;
typedef typename MeshType::FaceIterator FaceIterator;
typedef typename MeshType::VertexIterator VertexIterator;
typedef typename MeshType::VertContainer VertContainer;
typedef typename MeshType::VertexType VertexType;
typedef typename MeshType::VertexPointer VertexPointer;
typedef typename MeshType::VertexPointer VertexTypeP;
typedef vcg::face::VFIterator<FaceType> VFIteratorType;
typedef typename MeshType::CoordType CoordType;
typedef typename CoordType::ScalarType ScalarType;
class Quadric
{
public:
Quadric(double av, double bv, double cv, double dv, double ev)
{
a() = av;
b() = bv;
c() = cv;
d() = dv;
e() = ev;
}
double& a() { return data[0];}
double& b() { return data[1];}
double& c() { return data[2];}
double& d() { return data[3];}
double& e() { return data[4];}
double data[5];
double evaluate(double u, double v)
{
return a*u*u + b*u*v + c*v*v + d*u + e*v;
}
double du(double u, double v)
{
return 2.0*a()*u + b()*v + d();
}
double dv(double u, double v)
{
return 2.0*c()*v + b()*u + e();
}
double duv(double u, double v)
{
return b();
}
double duu(double u, double v)
{
return 2.0*a();
}
double dvv(double u, double v)
{
return 2.0*c();
}
static Quadric fit(std::vector<CoordType> VV)
{
assert(VV.size() >= 5);
Eigen::MatrixXf A(VV.size(),5);
Eigen::MatrixXf b(VV.size(),1);
Eigen::MatrixXf sol(VV.size(),1);
for(unsigned int c=0; c < VV.size(); ++c)
{
double u = VV[c].X();
double v = VV[c].Y();
double n = VV[c].Z();
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.transpose()*A).inverse()*A.transpose())*b;
return Quadric(sol[0],sol[1],sol[2],sol[3],sol[4]);
}
};
static CoordType project(VertexType* v, VertexType* vp)
{
return vp->P() - (v->N() * ((vp->P() - v->P()) * v->N()));
}
static std::vector<CoordType> computeReferenceFrames(VertexTypeP vi)
{
vcg::face::VFIterator<FaceType> vfi(vi);
int i = (vfi.I()+1)%3;
VertexTypeP vp = vfi.F()->V(i);
CoordType x = (project(&*vi,vp) - vi->P()).Normalize();
assert(fabs(x * vi->N()) < 0.1);
std::vector<CoordType> res(3);
res[0] = x;
res[1] = (vi->N() ^ res[0]).Normalize();
res[2] = (vi->N())/(vi->N()).Norm();
return res;
}
static std::set<CoordType> getSecondRing(VertexTypeP v)
{
std::set<VertexTypeP> ris;
std::set<CoordType> coords;
vcg::face::VFIterator<FaceType> vvi(v);
for(;!vvi.End();++vvi)
{
vcg::face::VFIterator<FaceType> vvi2(vvi.F()->V((vvi.I()+1)%3));
for(;!vvi2.End();++vvi2)
{
ris.insert(vvi2.F()->V((vvi2.I()+1)%3));
}
}
typename std::set<VertexTypeP>::iterator it;
for(it = ris.begin(); it != ris.end(); ++it)
coords.insert((*it)->P());
return coords;
}
static Quadric fitQuadric(VertexTypeP v, std::vector<CoordType>& ref)
{
std::set<CoordType> ring = getSecondRing(v);
if (ring.size() < 5)
return Quadric(1,1,1,1,1);
std::vector<CoordType> points;
typename std::set<CoordType>::iterator b = ring.begin();
typename std::set<CoordType>::iterator e = ring.end();
while(b != e)
{
// vtang non e` il v tangente!!!
CoordType vTang = *b - v->P();
double x = vTang * ref[0];
double y = vTang * ref[1];
double z = vTang * ref[2];
points.push_back(CoordType(x,y,z));
++b;
}
return Quadric::fit(points);
}
static void computeCurvature(MeshType & m)
{
assert(tri::HasPerVertexVFAdjacency(m) && tri::HasPerFaceVFAdjacency(m) );
vcg::tri::UpdateTopology<MeshType>::VertexFace(m);
vcg::tri::UpdateNormals<MeshType>::PerVertexAngleWeighted(m);
vcg::tri::UpdateNormals<MeshType>::NormalizeVertex(m);
VertexIterator vi;
for(vi = m.vert.begin(); vi!=m.vert.end(); ++vi )
{
std::vector<CoordType> ref = computeReferenceFrames(&*vi);
Quadric q = fitQuadric(&*vi,ref);
double a = q.a();
double b = q.b();
double c = q.c();
double d = q.d();
double e = q.e();
double E = 1.0 + d*d;
double F = d*e;
double G = 1.0 + e*e;
CoordType n = CoordType(-d,-e,1.0).Normalize();
vi->N() = ref[0] * n[0] + ref[1] * n[1] + ref[2] * n[2];
double L = 2.0 * a * n.Z();
double M = b * n.Z();
double N = 2 * c * n.Z();
// ----------------- 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();
c_val = -c_val;
CoordType v1, v2;
v1[0] = c_vec[0];
v1[1] = c_vec[1];
v1[2] = 0;
v2[0] = c_vec[2];
v2[1] = c_vec[3];
v2[2] = 0;
v1 = v1.Normalize();
v2 = v2.Normalize();
v1 = v1 * c_val[0];
v2 = v2 * c_val[1];
CoordType v1global = ref[0] * v1[0] + ref[1] * v1[1] + ref[2] * v1[2];
CoordType v2global = ref[0] * v2[0] + ref[1] * v2[1] + ref[2] * v2[2];
v1global.Normalize();
v2global.Normalize();
if (c_val[0] > c_val[1])
{
(*vi).PD1() = v1global;
(*vi).PD2() = v2global;
(*vi).K1() = c_val[0];
(*vi).K2() = c_val[1];
}
else
{
(*vi).PD1() = v2global;
(*vi).PD2() = v1global;
(*vi).K1() = c_val[1];
(*vi).K2() = c_val[0];
}
// ---- end Eigen stuff
}
}
// GG LOCAL CURVATURE
class QuadricLocal
{
public:
QuadricLocal ()
{
a() = b() = c() = d() = e() = 1.0;
}
QuadricLocal (double av, double bv, double cv, double dv, double ev)
{
a() = av;
b() = bv;
c() = cv;
d() = dv;
e() = ev;
}
double& a() { return data[0];}
double& b() { return data[1];}
double& c() { return data[2];}
double& d() { return data[3];}
double& e() { return data[4];}
double data[5];
double evaluate(double u, double v)
{
return a()*u*u + b()*u*v + c()*v*v + d()*u + e()*v;
}
double du(double u, double v)
{
return 2.0*a()*u + b()*v + d();
}
double dv(double u, double v)
{
return 2.0*c()*v + b()*u + e();
}
double duv(double u, double v)
{
return b();
}
double duu(double u, double v)
{
return 2.0*a();
}
double dvv(double u, double v)
{
return 2.0*c();
}
static QuadricLocal fit(std::vector<CoordType> &VV, bool svdRes, bool detCheck)
{
assert(VV.size() >= 5);
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].X();
double v = VV[c].Y();
double n = VV[c].Z();
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;
}
static int count = 0, index = 0;
double min = 0.000000000001; //1.0e-12
/*
if (!count)
printf("GNE %e\n", min);
*/
if (detCheck && ((A.transpose()*A).determinant() < min && (A.transpose()*A).determinant() > -min))
{
//A.svd().solve(b, &sol); A.svd().solve(b, &sol);
//cout << sol << endl;
printf("Quadric: unsolvable vertex %d %d\n", count, ++index);
//return Quadric (1, 1, 1, 1, 1);
A.svd().solve(b, &sol);
return QuadricLocal(sol[0],sol[1],sol[2],sol[3],sol[4]);
}
count++;
//for (int i = 0; i < 100; i++)
{
if (svdRes)
A.svd().solve(b, &sol);
else
sol = ((A.transpose()*A).inverse()*A.transpose())*b;
}
return QuadricLocal(sol[0],sol[1],sol[2],sol[3],sol[4]);
}
};
static void expandMaxLocal (MeshType & mesh, VertexType *v, int max, std::vector<VertexType*> *vv)
{
Nring<MeshType> rw = Nring<MeshType> (v, &mesh);
do rw.expand (); while (rw.allV.size() < max+1);
if (rw.allV[0] != v)
printf ("rw.allV[0] != *v\n");
vv->reserve ((size_t)max);
for (int i = 1; i < max+1; i++)
vv->push_back(rw.allV[i]);
rw.clear();
}
static void expandSphereLocal (MeshType & mesh, VertexType *v, float r, int min, std::vector<VertexType*> *vv)
{
Nring<MeshType> rw = Nring<MeshType> (v, &mesh);
bool isInside = true;
while (isInside)
{
rw.expand();
vv->reserve(rw.allV.size());
typename std::vector<VertexType*>::iterator b = rw.lastV.begin();
typename std::vector<VertexType*>::iterator e = rw.lastV.end();
isInside = false;
while(b != e)
{
if (((*b)->P() - v->P()).Norm() < r)
{
vv->push_back(*b);;
isInside = true;
}
++b;
}
}
//printf ("%d\n", vv->size());
rw.clear();
if (vv->size() < min)
{
vv->clear();
expandMaxLocal (mesh, v, min, vv);
}
}
static void getAverageNormal (VertexType *vp, std::vector<VertexType*> &vv, CoordType *ppn)
{
*ppn = CoordType (0,0,0);
for (typename std::vector<VertexType*>::iterator vpi = vv.begin(); vpi != vv.end(); ++vpi)
*ppn += (*vpi)->N();
*ppn += (*vp).N();
*ppn /= vv.size() + 1;
ppn->Normalize();
}
static void applyProjOnPlane (CoordType ppn, std::vector<VertexType*> &vin, std::vector<VertexType*> *vout)
{
for (typename std::vector<VertexType*>::iterator vpi = vin.begin(); vpi != vin.end(); ++vpi)
if ((*vpi)->N() * ppn > 0.0f)
vout->push_back (*vpi);
}
static CoordType projectLocal(VertexType* v, VertexType* vp, CoordType ppn)
{
return vp->P() - (ppn * ((vp->P() - v->P()) * ppn));
}
static void computeReferenceFramesLocal (VertexType *v, CoordType ppn, std::vector<CoordType> *ref)
{
vcg::face::VFIterator<FaceType> vfi (v);
int i = (vfi.I() + 1) % 3;
VertexTypeP vp = vfi.F()->V(i);
CoordType x = (projectLocal (v, vp, ppn) - v->P()).Normalize();
assert(fabs(x * ppn) < 0.1);
*ref = std::vector<CoordType>(3);
(*ref)[0] = x;
(*ref)[1] = (ppn ^ (*ref)[0]).Normalize();
(*ref)[2] = ppn.Normalize(); //ppn / ppn.Norm();
}
static void fitQuadricLocal (VertexType *v, std::vector<CoordType> ref, std::vector<VertexType*> &vv, QuadricLocal *q)
{
bool svdResolution = false;
bool zeroDeterminantCheck = false;
std::vector<CoordType> points;
points.reserve (vv.size());
typename std::vector<VertexType*>::iterator b = vv.begin();
typename std::vector<VertexType*>::iterator e = vv.end();
while(b != e)
{
CoordType cp = (*b)->P();
// vtang non e` il v tangente!!!
CoordType vTang = cp - v->P();
double x = vTang * ref[0];
double y = vTang * ref[1];
double z = vTang * ref[2];
points.push_back(CoordType(x,y,z));
++b;
}
*q = QuadricLocal::fit (points, svdResolution, zeroDeterminantCheck);
}
static void finalEigenStuff (VertexType *v, std::vector<CoordType> ref, QuadricLocal q)
{
double a = q.a();
double b = q.b();
double c = q.c();
double d = q.d();
double e = q.e();
double E = 1.0 + d*d;
double F = d*e;
double G = 1.0 + e*e;
CoordType n = CoordType(-d,-e,1.0).Normalize();
v->N() = ref[0] * n[0] + ref[1] * n[1] + ref[2] * n[2];
double L = 2.0 * a * n.Z();
double M = b * n.Z();
double N = 2 * c * n.Z();
// ----------------- 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();
c_val = -c_val;
CoordType v1, v2;
v1[0] = c_vec[0];
v1[1] = c_vec[1];
v1[2] = d * v1[0] + e * v1[1];
v2[0] = c_vec[2];
v2[1] = c_vec[3];
v2[2] = d * v2[0] + e * v2[1];
v1 = v1.Normalize();
v2 = v2.Normalize();
CoordType v1global = ref[0] * v1[0] + ref[1] * v1[1] + ref[2] * v1[2];
CoordType 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])
{
(*v).PD1() = v1global;
(*v).PD2() = v2global;
(*v).K1() = c_val[0];
(*v).K2() = c_val[1];
}
else
{
(*v).PD1() = v2global;
(*v).PD2() = v1global;
(*v).K1() = c_val[1];
(*v).K2() = c_val[0];
}
// ---- end Eigen stuff
}
static void updateCurvatureLocal (MeshType & mesh, float radiusSphere)
{
bool verbose = false;
bool projectionPlaneCheck = true;
int vertexesPerFit = 0;
int i = 0;
VertexIterator vi;
for(vi = mesh.vert.begin(); vi != mesh.vert.end(); ++vi, i++)
{
std::vector<VertexType*> vv;
std::vector<VertexType*> vvtmp;
int count;
if (verbose && !((count = (vi - mesh.vert.begin())) % 1000))
printf ("vertex %d of %d\n",count,mesh.vert.size());
// if (kRing != 0)
// expandRing (&*vi, kRing, 5, &vv);
// else
expandSphereLocal (mesh, &*vi, radiusSphere, 5, &vv);
assert (vv.size() >= 5);
CoordType ppn;
// if (averageNormalMode)
// //ppn = (*vi).N();
getAverageNormal (&*vi, vv, &ppn);
// else
// getProjPlaneNormal (&*vi, vv, &ppn);
if (projectionPlaneCheck)
{
vvtmp.reserve (vv.size ());
applyProjOnPlane (ppn, vv, &vvtmp);
if (vvtmp.size() >= 5)
vv = vvtmp;
}
vvtmp.clear();
// if (montecarloMaxVertexNum)
// {
// //printf ("P: %d\n", vv.size());
// vvtmp.reserve (vv.size ());
// //printf ("TP: %d\n", vvtmp.size());
// applyMontecarlo (montecarloMaxVertexNum, vv, &vvtmp);
// //printf ("TD: %d\n", vvtmp.size());
// vv = vvtmp;
// //printf ("D: %d\n", vv.size());
// //printf ("\n");
// }
assert (vv.size() >= 5);
std::vector<CoordType> ref;
computeReferenceFramesLocal (&*vi, ppn, &ref);
/*
printf ("%lf %lf %lf - %lf %lf %lf - %lf %lf %lf\n",
ref[0][0], ref[0][1], ref[0][2],
ref[1][0], ref[1][1], ref[1][2],
ref[2][0], ref[2][1], ref[2][2]);
*/
vertexesPerFit += vv.size();
//printf ("size: %d\n", vv.size());
QuadricLocal q;
fitQuadricLocal (&*vi, ref, vv, &q);
finalEigenStuff (&*vi, ref, q);
}
//if (verbose)
//printf ("average vertex num in each fit: %f, total %d, vn %d\n", ((float) vertexesPerFit) / mesh.vn, vertexesPerFit, mesh.vn);
if (verbose)
printf ("average vertex num in each fit: %f\n", ((float) vertexesPerFit) / mesh.vn);
}
};
}
}
#endif

View File

@ -1,109 +0,0 @@
/****************************************************************************
* 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.4 2006/05/16 21:36:54 cignoni
Removed unused box function and rewrote initial comment.
Revision 1.3 2006/05/15 13:12:36 pietroni
Updating of edge values id divided into 2 functions ( the first one update only a face...) added also resetting of edges flags.. (see first line of Set function)
Revision 1.2 2004/05/12 18:52:35 ganovelli
removed call to ComputeRT and put its body here
Revision 1.1 2004/05/12 10:39:45 ganovelli
created
****************************************************************************/
#ifndef __VCG_TRI_UPDATE_EDGES
#define __VCG_TRI_UPDATE_EDGES
#include <vcg/space/plane3.h>
namespace vcg {
namespace tri {
/// \ingroup trimesh
/// \headerfile edges.h vcg/complex/trimesh/update/edges.h
/// \brief This class is used to compute or update the precomputed data used to efficiently compute point-face distances.
template <class ComputeMeshType>
class UpdateEdges
{
public:
typedef ComputeMeshType MeshType;
typedef typename MeshType::VertexType VertexType;
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::FaceType::CoordType::ScalarType ScalarType;
static void Set(FaceType &f)
{
f.Flags() = f.Flags() & (~(FaceType::NORMX|FaceType::NORMY|FaceType::NORMZ));
// Primo calcolo degli edges
f.Edge(0) = f.V(1)->P(); f.Edge(0) -= f.V(0)->P();
f.Edge(1) = f.V(2)->P(); f.Edge(1) -= f.V(1)->P();
f.Edge(2) = f.V(0)->P(); f.Edge(2) -= f.V(2)->P();
// Calcolo di plane
f.Plane().SetDirection(f.Edge(0)^f.Edge(1));
f.Plane().SetOffset(f.Plane().Direction().dot(f.V(0)->P()));
f.Plane().Normalize();
// Calcolo migliore proiezione
ScalarType nx = math::Abs(f.Plane().Direction()[0]);
ScalarType ny = math::Abs(f.Plane().Direction()[1]);
ScalarType nz = math::Abs(f.Plane().Direction()[2]);
ScalarType d;
if(nx>ny && nx>nz) { f.Flags() |= FaceType::NORMX; d = 1/f.Plane().Direction()[0]; }
else if(ny>nz) { f.Flags() |= FaceType::NORMY; d = 1/f.Plane().Direction()[1]; }
else { f.Flags() |= FaceType::NORMZ; d = 1/f.Plane().Direction()[2]; }
// Scalatura spigoli
f.Edge(0)*=d;
f.Edge(1)*=d;
f.Edge(2)*=d;
}
static void Set(ComputeMeshType &m)
{
FaceIterator f;
for(f = m.face.begin(); f!=m.face.end(); ++f)
if(!(*f).IsD())
Set(*f);
}
}; // end class
} // End namespace
} // End namespace
#endif

View File

@ -1,488 +0,0 @@
/****************************************************************************
* 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. *
* *
****************************************************************************/
#ifndef FITMAPS_H
#define FITMAPS_H
#include <vcg/math/histogram.h>
#include <vcg/simplex/face/jumping_pos.h>
#include <vcg/complex/trimesh/allocate.h>
#include <vcg/complex/trimesh/update/flag.h>
#include <vector>
#include <set>
#include <vcg/complex/trimesh/update/normal.h>
#include <vcg/complex/trimesh/update/curvature.h>
#include <vcg/complex/trimesh/update/topology.h>
#include <vcg/complex/trimesh/update/bounding.h>
#include "vcg/complex/trimesh/update/curvature_fitting.h"
#include <vcg/Eigen/Core>
#include <vcg/Eigen/QR>
#include <vcg/Eigen/LU>
#include <vcg/Eigen/SVD>
#include <algorithm>
#include <vcg/complex/trimesh/nring.h>
#include <vcg/complex/trimesh/smooth.h>
using namespace Eigen;
namespace vcg { namespace tri {
template<class MeshType>
class Fitmaps
{
public:
typedef typename MeshType::FaceType FaceType;
typedef typename MeshType::VertexType VertexType;
typedef typename MeshType::ScalarType ScalarType;
typedef typename MeshType::FaceIterator FaceIterator;
typedef typename MeshType::VertexIterator VertexIterator;
typedef typename MeshType::CoordType CoordType;
typedef vcg::tri::Nring<MeshType> RingWalker;
class Bicubic
{
public:
Bicubic() {};
Bicubic(vector<double>& input)
{
data = input;
if (input.size() != 16)
{
assert(0);
}
}
// (u3 u2 u 1) (v3 v2 v 1)
//
// a u3 v3
// b u3 v2
// c u3 v1
// d u3 1
// e u2 v3
// f u2 v2
// g u2 v1
// h u2 1
// i u1 v3
// l u1 v2
// m u1 v1
// n u1 1
// o 1 v3
// p 1 v2
// q 1 v1
// r 1 1
double& a() { return data[0];}
double& b() { return data[1];}
double& c() { return data[2];}
double& d() { return data[3];}
double& e() { return data[4];}
double& f() { return data[5];}
double& g() { return data[6];}
double& h() { return data[7];}
double& i() { return data[8];}
double& l() { return data[9];}
double& m() { return data[10];}
double& n() { return data[11];}
double& o() { return data[12];}
double& p() { return data[13];}
double& q() { return data[14];}
double& r() { return data[15];}
vector<double> data;
double evaluate(double u, double v)
{
return
a() * u*u*u * v*v*v +
b() * u*u*u * v*v +
c() * u*u*u * v +
d() * u*u*u +
e() * u*u * v*v*v +
f() * u*u * v*v +
g() * u*u * v +
h() * u*u +
i() * u * v*v*v +
l() * u * v*v +
m() * u * v +
n() * u +
o() * v*v*v +
p() * v*v +
q() * v +
r();
}
double distanceRMS(std::vector<CoordType>& VV)
{
double error = 0;
for(typename std::vector<CoordType>::iterator it = VV.begin(); it != VV.end(); ++it)
{
double u = it->X();
double v = it->Y();
double n = it->Z();
double temp = evaluate(u,v);
error += (n - temp)*(n - temp);
}
error /= (double) VV.size();
return sqrt(error);
}
static Bicubic fit(std::vector<CoordType> VV)
{
assert(VV.size() >= 16);
Eigen::MatrixXf A(VV.size(),16);
Eigen::MatrixXf b(VV.size(),1);
Eigen::MatrixXf sol(16,1);
for(unsigned int c=0; c < VV.size(); ++c)
{
double u = VV[c].X();
double v = VV[c].Y();
double n = VV[c].Z();
A(c,0) = u*u*u * v*v*v;
A(c,1) = u*u*u * v*v;
A(c,2) = u*u*u * v;
A(c,3) = u*u*u;
A(c,4) = u*u * v*v*v;
A(c,5) = u*u * v*v;
A(c,6) = u*u * v;
A(c,7) = u*u;
A(c,8) = u * v*v*v;
A(c,9) = u * v*v;
A(c,10) = u * v;
A(c,11) = u;
A(c,12) = v*v*v;
A(c,13) = v*v;
A(c,14) = v;
A(c,15) = 1;
b[c] = n;
}
A.svd().solve(b, &sol);
vector<double> r(16);
for (int i=0; i < 16; ++i)
r.at(i) = sol[i];
return Bicubic(r);
}
};
Fitmaps()
{}
class radSorter
{
public:
radSorter(VertexType* v)
{
origin = v;
}
VertexType* origin;
bool operator() (VertexType* v1, VertexType* v2)
{
return (v1->P() - origin->P()).SquaredNorm() < (v2->P() - origin->P()).SquaredNorm();
}
};
float getMeanCurvature(VertexType* vp)
{
return (vp->K1() + vp->K2())/2.0;
}
static bool fitBicubicPoints(VertexType* v, std::vector<CoordType>& ref, Bicubic& ret, std::vector<CoordType>& points, std::vector<VertexType*>& ring)
{
points.clear();
if (ring.size() < 16)
{
return false;
}
typename std::vector<VertexType*>::iterator b = ring.begin();
typename std::vector<VertexType*>::iterator e = ring.end();
while(b != e)
{
CoordType vT = (*b)->P() - v->P();
double x = vT * ref[0];
double y = vT * ref[1];
double z = vT * ref[2];
points.push_back(CoordType(x,y,z));
++b;
}
ret = Bicubic::fit(points);
return true;
}
static double AverageEdgeLenght(MeshType& m)
{
double doubleA = 0;
for (FaceIterator fi = m.face.begin(); fi!=m.face.end(); fi++) if (!fi->IsD()) {
doubleA+=vcg::DoubleArea(*fi);
}
int nquads = m.fn / 2;
return sqrt( doubleA/(2*nquads) );
}
static void computeMFitmap(MeshType& m, float perc, int ringMax = 50)
{
vcg::tri::UpdateCurvatureFitting<MeshType>::computeCurvature(m);
vcg::tri::UpdateNormals<MeshType>::PerVertexAngleWeighted(m);
vcg::tri::UpdateTopology<MeshType>::FaceFace(m);
vcg::tri::UpdateTopology<MeshType>::VertexFace(m);
vcg::tri::UpdateBounding<MeshType>::Box(m);
int countTemp = 0;
RingWalker::clearFlags(&m);
for(VertexIterator it=m.vert.begin(); it!=m.vert.end();++it)
{
if ((countTemp++ % 100) == 0)
cerr << countTemp << "/" << m.vert.size() << endl;
RingWalker rw(&*it,&m);
CoordType nor = it->N();
float okfaces = 0;
float flipfaces = 0;
int count = 0;
do
{
count++;
rw.expand();
for(unsigned i=0; i<rw.lastF.size();++i)
{
CoordType vet1 = nor;
CoordType vet2 = rw.lastF[i]->N();
vet1.Normalize();
vet2.Normalize();
double scal = vet1 * vet2;
if ((scal) > 0)
okfaces += (vcg::DoubleArea(*rw.lastF[i]));
else
flipfaces += (vcg::DoubleArea(*rw.lastF[i]));
}
} while ((((flipfaces)/(okfaces + flipfaces))*100.0 < perc) && (count < ringMax));
std::sort(rw.lastV.begin(),rw.lastV.end(),radSorter(&*it));
it->Q() = ((*rw.lastV.begin())->P() - it->P()).Norm();
rw.clear();
}
vcg::tri::Smooth<MeshType>::VertexQualityLaplacian(m,2);
}
static vector<VertexType*> gatherNeighsSurface(VertexType* vt, float sigma, MeshType& m)
{
vector<VertexType*> current;
RingWalker rw(vt,&m);
bool exit = false;
do
{
rw.expand();
exit = true;
for(typename vector<VertexType*>::iterator it = rw.lastV.begin(); it != rw.lastV.end(); ++it)
{
if (((*it)->P() - vt->P()).Norm() < sigma)
{
current.push_back(*it);
exit = false;
}
}
} while (!exit);
rw.clear();
return current;
}
static void computeSFitmap(MeshType& m)//, float target = 1000)
{
vcg::tri::UpdateCurvatureFitting<MeshType>::computeCurvature(m);
vcg::tri::UpdateNormals<MeshType>::PerVertexAngleWeighted(m);
vcg::tri::UpdateTopology<MeshType>::FaceFace(m);
vcg::tri::UpdateTopology<MeshType>::VertexFace(m);
// update bounding box
vcg::tri::UpdateBounding<MeshType>::Box(m);
int countTemp = 0;
double e = AverageEdgeLenght(m);
int iteraz = 5; //2.0 * sqrt(m.vert.size()/target);
for(VertexIterator it=m.vert.begin(); it!=m.vert.end();++it)
{
if ((countTemp++ % 100) == 0)
cerr << countTemp << "/" << m.vert.size() << endl;
vector<float> oneX;
for (int iteration = 0; iteration<iteraz; ++iteration)
{
oneX.push_back((iteration+1)*(e));
}
std::vector<CoordType> ref(3);
ref[0] = it->PD1();
ref[1] = it->PD2();
ref[2] = it->PD1() ^ it->PD2();
ref[0].Normalize();
ref[1].Normalize();
ref[2].Normalize();
Bicubic b;
RingWalker::clearFlags(&m);
std::vector<VertexType*> pointsGlobal = gatherNeighsSurface(&*it,oneX.at(iteraz-1),m);
vector<float> onedimensional;
for (int iteration = 0; iteration<iteraz; ++iteration)
{
std::vector<VertexType*> points; // solo quelli nel raggio
std::vector<CoordType> projected; // come sopra ma in coord locali
for (typename std::vector<VertexType*>::iterator it2 = pointsGlobal.begin(); it2 != pointsGlobal.end(); ++it2)
{
if (((*it).P() - (*it2)->P()).Norm() < oneX.at(iteration))
points.push_back(*it2);
}
std::vector<VertexType*>& pointsFitting = points;
if (!fitBicubicPoints(&*it, ref, b, projected,pointsFitting))
{
onedimensional.push_back(0);
}
else
{
onedimensional.push_back(b.distanceRMS(projected));
}
}
// // vecchio fit ax^4
Eigen::MatrixXf Am(onedimensional.size(),1);
Eigen::MatrixXf bm(onedimensional.size(),1);
Eigen::MatrixXf sol(1,1);
for(unsigned int c=0; c < onedimensional.size(); ++c)
{
double x = oneX.at(c);
Am(c,0) = pow(x,4);
bm[c] = onedimensional[c];
}
Am.svd().solve(bm, &sol);
it->Q() = pow((double)sol[0],0.25);
// // nuovo fit ax^4 + b
// Eigen::MatrixXf Am(onedimensional.size()+1,2);
// Eigen::MatrixXf bm(onedimensional.size()+1,1);
// Eigen::MatrixXf sol(2,1);
//
// Am(0,0) = 0;
// Am(0,1) = 0;
// bm[0] = 0;
//
// for(unsigned int c=0; c < onedimensional.size(); ++c)
// {
// double x = oneX.at(c);
//
// Am(c,0) = pow(x,4);
// Am(c,1) = 1;
// bm[c] = onedimensional[c];
// }
//
// //sol = ((Am.transpose()*Am).inverse()*Am.transpose())*bm;
// Am.svd().solve(bm, &sol);
//
// cerr << "------" << sol[0] << " " << sol[1] << endl;
// if (sol[0] > 0)
// saliency[it] = pow((double)sol[0],0.25);
// else
// saliency[it] = 0;
}
vcg::tri::Smooth<MeshType>::VertexQualityLaplacian(m,1);
}
~Fitmaps(){};
};
}} // END NAMESPACES
#endif // FITMAPS_H

View File

@ -1,453 +0,0 @@
/****************************************************************************
* 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.20 2007/05/31 15:24:50 ponchio
FIxed off-by-one error on FaceBorderFromNone.
Revision 1.19 2007/05/22 15:19:42 cignoni
Added VertexClear
Revision 1.18 2007/01/30 18:49:23 tarini
aggiunta la VertexBorderFromNone (flag bordo per vertici senza richiedere nulla)
Revision 1.17 2006/08/31 13:11:12 marfr960
corrected bounds of a vector scan
Revision 1.16 2006/08/30 12:59:49 marfr960
Added missing std:: to swap
Revision 1.15 2006/08/30 06:50:07 cignoni
Reverted to version 1.13. Version 1.14 was done on outdated version.
Revision 1.13 2006/06/18 20:49:30 cignoni
Added missing IsD tests
Revision 1.12 2006/05/03 21:23:25 cignoni
Corrected IsDeleted -> isD
Revision 1.11 2005/12/02 00:09:12 cignoni
Added assert(HasFlags) everywhere..
Revision 1.10 2005/07/06 08:16:34 ganovelli
set VertexBorderFromFace as static
Revision 1.9 2005/06/10 15:07:23 cignoni
Completed FaceBorderFromNone (and added a missing helper class)
Revision 1.8 2005/04/01 13:04:55 fiorin
Minor changes
Revision 1.7 2004/09/14 19:49:43 ganovelli
first compilation version
Revision 1.6 2004/07/15 00:13:39 cignoni
Better doxigen documentation
Revision 1.5 2004/07/06 06:27:02 cignoni
Added FaceBorderFromVF
Revision 1.4 2004/05/13 15:58:55 ganovelli
function Clear added
Revision 1.3 2004/03/12 15:22:19 cignoni
Written some documentation and added to the trimes doxygen module
Revision 1.2 2004/03/10 00:46:10 cignoni
changed to the face::IsBorder() style
Revision 1.1 2004/03/05 10:59:24 cignoni
Changed name from plural to singular (normals->normal)
Revision 1.1 2004/03/04 00:37:56 cignoni
First working version!
****************************************************************************/
#ifndef __VCG_TRI_UPDATE_FLAGS
#define __VCG_TRI_UPDATE_FLAGS
#include <vcg/simplex/face/pos.h>
namespace vcg {
namespace tri {
/// \ingroup trimesh
/// \headerfile flag.h vcg/complex/trimesh/update/flag.h
/// \brief Management, updating and computation of per-vertex and per-face flags (like border flags).
/**
This class is used to compute or update some of the flags that can be stored in the mesh components. For now just Border flags (e.g. the flag that tells if a given edge of a face belong to a border of the mesh or not).
*/
template <class UpdateMeshType>
class UpdateFlags
{
public:
typedef UpdateMeshType MeshType;
typedef vcg::face::Pos<typename UpdateMeshType::FaceType> PosType;
typedef typename MeshType::ScalarType ScalarType;
typedef typename MeshType::VertexType VertexType;
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;
/// \brief Reset all the mesh flags (both vertexes and faces) setting everithing to zero (the default value for flags)
static void Clear(MeshType &m)
{
assert(HasPerFaceFlags(m));
FaceIterator fi;
VertexIterator vi;
for(fi=m.face.begin(); fi!=m.face.end(); ++fi)
(*fi).Flags() = 0;
for(vi=m.vert.begin(); vi!=m.vert.end(); ++vi)
(*vi).Flags() = 0;
}
static void VertexClear(MeshType &m, unsigned int FlagMask = 0xffffffff)
{
VertexIterator vi;
int andMask = ~FlagMask;
for(vi=m.vert.begin(); vi!=m.vert.end(); ++vi)
if(!(*vi).IsD()) (*vi).Flags() &= andMask ;
}
static void FaceClear(MeshType &m, unsigned int FlagMask = 0xffffffff)
{
FaceIterator fi;
int andMask = ~FlagMask;
for(fi=m.face.begin(); fi!=m.face.end(); ++fi)
if(!(*fi).IsD()) (*fi).Flags() &= andMask ;
}
static void VertexSet(MeshType &m, unsigned int FlagMask)
{
VertexIterator vi;
for(vi=m.vert.begin(); vi!=m.vert.end(); ++vi)
if(!(*vi).IsD()) (*vi).Flags() |= FlagMask ;
}
static void FaceSet(MeshType &m, unsigned int FlagMask)
{
FaceIterator fi;
for(fi=m.face.begin(); fi!=m.face.end(); ++fi)
if(!(*fi).IsD()) (*fi).Flags() |= FlagMask ;
}
static void VertexClearV(MeshType &m) { VertexClear(m,VertexType::VISITED);}
static void VertexClearB(MeshType &m) { VertexClear(m,VertexType::BORDER);}
static void FaceClearV(MeshType &m) { FaceClear(m,FaceType::VISITED);}
static void FaceClearB(MeshType &m) { FaceClear(m,FaceType::BORDER012);}
static void FaceClearF(MeshType &m) { FaceClear(m,FaceType::FAUX012);}
static void VertexSetV(MeshType &m) { VertexSet(m,VertexType::VISITED);}
static void VertexSetB(MeshType &m) { VertexSet(m,VertexType::BORDER);}
static void FaceSetV(MeshType &m) { FaceSet(m,FaceType::VISITED);}
static void FaceSetB(MeshType &m) { FaceSet(m,FaceType::BORDER);}
static void FaceSetF(MeshType &m) { FaceSet(m,FaceType::FAUX012);}
/// \brief Compute the border flags for the faces using the Face-Face Topology.
/**
\warning Obviously it assumes that the topology has been correctly computed (see: UpdateTopology::FaceFace )
*/
static void FaceBorderFromFF(MeshType &m)
{
assert(HasPerFaceFlags(m));
// const int BORDERFLAG[3]={FaceType::BORDER0,FaceType::BORDER1,FaceType::BORDER2};
FaceIterator fi;
for(fi=m.face.begin();fi!=m.face.end();++fi)if(!(*fi).IsD())
for(int j=0;j<3;++j)
{
//if(!(*fi).IsManifold(j)) (*fi).SetCF(j);
//else
if(face::IsBorder(*fi,j)) (*fi).SetB(j);
else (*fi).ClearB(j);
}
}
static void FaceBorderFromVF(MeshType &m)
{
assert(HasPerFaceFlags(m));
VertexIterator vi;
assert(m.HasVFTopology());
int visitedBit=VertexType::NewBitFlag();
// Calcolo dei bordi
// per ogni vertice vi si cercano i vertici adiacenti che sono toccati da una faccia sola
// (o meglio da un numero dispari di facce)
const int BORDERFLAG[3]={FaceType::BORDER0, FaceType::BORDER1, FaceType::BORDER2};
for(vi=m.vert.begin();vi!=m.vert.end();++vi)
if(!(*vi).IsD())
{
for(face::VFIterator<FaceType> vfi(&*vi) ; !vfi.End(); ++vfi )
{
vfi.f->V1(vfi.z)->ClearUserBit(visitedBit);
vfi.f->V2(vfi.z)->ClearUserBit(visitedBit);
}
for(face::VFIterator<FaceType> vfi(&*vi) ; !vfi.End(); ++vfi )
{
if(vfi.f->V1(vfi.z)->IsUserBit(visitedBit)) vfi.f->V1(vfi.z)->ClearUserBit(visitedBit);
else vfi.f->V1(vfi.z)->SetUserBit(visitedBit);
if(vfi.f->V2(vfi.z)->IsUserBit(visitedBit)) vfi.f->V2(vfi.z)->ClearUserBit(visitedBit);
else vfi.f->V2(vfi.z)->SetUserBit(visitedBit);
}
for(face::VFIterator<FaceType> vfi(&*vi) ; !vfi.End(); ++vfi )
{
if(vfi.f->V(vfi.z)< vfi.f->V1(vfi.z) && vfi.f->V1(vfi.z)->IsUserBit(visitedBit))
vfi.f->Flags() |= BORDERFLAG[vfi.z];
if(vfi.f->V(vfi.z)< vfi.f->V2(vfi.z) && vfi.f->V2(vfi.z)->IsUserBit(visitedBit))
vfi.f->Flags() |= BORDERFLAG[(vfi.z+2)%3];
}
}
VertexType::DeleteBitFlag(VertexType::LastBitFlag());
}
class EdgeSorter
{
public:
VertexPointer v[2]; // Puntatore ai due vertici (Ordinati)
FacePointer f; // Puntatore alla faccia generatrice
int z; // Indice dell'edge nella faccia
EdgeSorter() {} // Nothing to do
void Set( const FacePointer pf, const int nz )
{
assert(pf!=0);
assert(nz>=0);
assert(nz<3);
v[0] = pf->V(nz);
v[1] = pf->V((nz+1)%3);
assert(v[0] != v[1]);
if( v[0] > v[1] ) std::swap(v[0],v[1]);
f = pf;
z = nz;
}
inline bool operator < ( const EdgeSorter & pe ) const {
if( v[0]<pe.v[0] ) return true;
else if( v[0]>pe.v[0] ) return false;
else return v[1] < pe.v[1];
}
inline bool operator == ( const EdgeSorter & pe ) const
{
return v[0]==pe.v[0] && v[1]==pe.v[1];
}
inline bool operator != ( const EdgeSorter & pe ) const
{
return v[0]!=pe.v[0] || v[1]!=pe.v[1];
}
};
// versione minimale che non calcola i complex flag.
static void VertexBorderFromNone(MeshType &m)
{
assert(HasPerVertexFlags(m));
std::vector<EdgeSorter> e;
typename UpdateMeshType::FaceIterator pf;
typename std::vector<EdgeSorter>::iterator p;
if( m.fn == 0 )
return;
e.resize(m.fn*3); // Alloco il vettore ausiliario
p = e.begin();
for(pf=m.face.begin();pf!=m.face.end();++pf) // Lo riempio con i dati delle facce
if( ! (*pf).IsD() )
for(int j=0;j<3;++j)
{
(*p).Set(&(*pf),j);
(*pf).ClearB(j);
++p;
}
assert(p==e.end());
sort(e.begin(), e.end()); // Lo ordino per vertici
typename std::vector<EdgeSorter>::iterator pe,ps;
for(ps = e.begin(), pe = e.begin(); pe < e.end(); ++pe) // Scansione vettore ausiliario
{
if( pe==e.end() || *pe != *ps ) // Trovo blocco di edge uguali
{
if(pe-ps==1) {
ps->v[0]->SetB();
ps->v[1]->SetB();
} else
if(pe-ps!=2) { // not twomanyfold!
for(;ps!=pe;++ps) {
ps->v[0]->SetB(); // Si settano border anche i complex.
ps->v[1]->SetB();
}
}
ps = pe;
}
}
}
/// This function fill the flags with the info on what is the best projection direction
/// for a given face. Used by the point-face distance function when do not exploiting pre-computed
/// per-face data (the so called edge component)
static void FaceProjection(MeshType &m)
{
FaceIterator fi;
for(fi=m.face.begin();fi!=m.face.end();++fi) // Lo riempio con i dati delle facce
if( ! (*fi).IsD() )
{
ScalarType nx = math::Abs((*fi).cN()[0]);
ScalarType ny = math::Abs((*fi).cN()[1]);
ScalarType nz = math::Abs((*fi).cN()[2]);
if(nx>ny && nx>nz) { (*fi).Flags() |= FaceType::NORMX; }
else if(ny>nz) { (*fi).Flags() |= FaceType::NORMY; }
else { (*fi).Flags() |= FaceType::NORMZ; }
}
}
/// Computes per-face border flags without requiring any kind of topology
/// It has a O(fn log fn) complexity.
static void FaceBorderFromNone(MeshType &m)
{
assert(HasPerFaceFlags(m));
std::vector<EdgeSorter> e;
typename UpdateMeshType::FaceIterator pf;
typename std::vector<EdgeSorter>::iterator p;
for(VertexIterator v=m.vert.begin();v!=m.vert.end();++v)
(*v).ClearB();
if( m.fn == 0 )
return;
FaceIterator fi;
int n_edges = 0;
for(fi = m.face.begin(); fi != m.face.end(); ++fi) if(! (*fi).IsD()) n_edges+=(*fi).VN();
e.resize(n_edges);
p = e.begin();
for(pf=m.face.begin();pf!=m.face.end();++pf) // Lo riempio con i dati delle facce
if( ! (*pf).IsD() )
for(int j=0;j<(*pf).VN();++j)
{
(*p).Set(&(*pf),j);
(*pf).ClearB(j);
++p;
}
assert(p==e.end());
sort(e.begin(), e.end()); // Lo ordino per vertici
typename std::vector<EdgeSorter>::iterator pe,ps;
ps = e.begin();pe=e.begin();
do
{
if( pe==e.end() || *pe != *ps ) // Trovo blocco di edge uguali
{
if(pe-ps==1) {
ps->f->SetB(ps->z);
} else
if(pe-ps!=2) { // Caso complex!!
for(;ps!=pe;++ps)
ps->f->SetB(ps->z); // Si settano border anche i complex.
}
ps = pe;
}
if(pe==e.end()) break;
++pe;
} while(true);
// TRACE("found %i border (%i complex) on %i edges\n",nborder,ncomplex,ne);
}
/// Compute the PerVertex Border flag deriving it from the faces
static void VertexBorderFromFace(MeshType &m)
{
assert(HasPerFaceFlags(m));
typename MeshType::VertexIterator v;
typename MeshType::FaceIterator f;
for(v=m.vert.begin();v!=m.vert.end();++v)
(*v).ClearB();
for(f=m.face.begin();f!=m.face.end();++f)
if(!(*f).IsD())
{
for(int z=0;z<(*f).VN();++z)
if( (*f).IsB(z) )
{
(*f).V(z)->SetB();
(*f).V((*f).Next(z))->SetB();
}
}
}
//
static void FaceFauxCrease(MeshType &m,float AngleRad)
{
assert(HasPerFaceFlags(m));
assert(HasFFAdjacency(m));
typename MeshType::FaceIterator f;
//initially everything is faux (e.g all internal)
FaceSetF(m);
for(f=m.face.begin();f!=m.face.end();++f)
{
if(!(*f).IsD())
{
for(int z=0;z<(*f).VN();++z)
{
if( face::IsBorder(*f,z) ) (*f).ClearF(z);
else
{
if(Angle((*f).N(), (*f).FFp(z)->N()) > AngleRad)
(*f).ClearF(z);
}
}
}
}
}
}; // end class
} // End namespace tri
} // End namespace vcg
#endif

View File

@ -1,666 +0,0 @@
/****************************************************************************
* 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. *
* *
****************************************************************************/
#ifndef __VCGLIB_HALFEDGE_
#define __VCGLIB_HALFEDGE_
#include <vector>
#include <vcg/complex/trimesh/allocate.h>
#include <vcg/complex/trimesh/clean.h>
#include <vcg/complex/trimesh/update/topology.h>
#include <vcg/complex/trimesh/base.h>
#include <vcg/complex/trimesh/update/halfedge_topology.h>
namespace vcg
{
namespace tri{
/// \ingroup trimesh
/// \headerfile edge_support.h vcg/complex/trimesh/edge_support.h
/// \brief This class is used to build edge based data structure from indexed data structure and viceversa
/**
*/
template <class MeshType >
class UpdateHalfEdges{
public:
typedef typename MeshType::VertexType VertexType;
typedef typename MeshType::VertexPointer VertexPointer;
typedef typename MeshType::VertexIterator VertexIterator;
typedef typename MeshType::HEdgePointer HEdgePointer;
typedef typename MeshType::HEdgeType HEdgeType;
typedef typename MeshType::EdgePointer EdgePointer;
typedef typename MeshType::EdgeType EdgeType;
typedef typename MeshType::EdgeIterator EdgeIterator;
typedef typename MeshType::HEdgeIterator HEdgeIterator;
typedef typename MeshType::FaceIterator FaceIterator;
typedef typename MeshType::FaceType FaceType;
struct VertexPairEdgePtr{
VertexPairEdgePtr(VertexPointer _v0,VertexPointer _v1,HEdgePointer _ep):v0(_v0),v1(_v1),ep(_ep){if(v0>v1) std::swap(v0,v1);}
const bool operator <(const VertexPairEdgePtr &o) const {return (v0 == o.v0)? (v1<o.v1):(v0<o.v0);}
const bool operator ==(const VertexPairEdgePtr &o) const {return (v0 == o.v0)&& (v1==o.v1);}
VertexPointer v0,v1;
HEdgePointer ep;
};
struct FacePtrInt{
FacePtrInt ( FaceType * _f,int _i):f(_f),i(_i){}
FaceType * f;
int i;
};
typedef std::vector<bool> BitVector;
/**
build a half-edge data structure from an indexed data structure. Note that the half-edges are allocated here for the first time.
If you have a mesh where there are already edges, they will be removed and the data lost, so do not use this function
to just "update" the topology of half edges.
**/
static void FromIndexed(MeshType & m){
assert(HasFVAdjacency(m));
assert(HasHOppAdjacency(m));
assert(HasHNextAdjacency(m));
typename MeshType::template PerFaceAttributeHandle<BitVector> flagVisited =
vcg::tri::Allocator<MeshType>::template AddPerFaceAttribute<BitVector>(m,"");
std::vector<FacePtrInt > borderEdges;
// allocate all new half edges
FaceIterator fi;
unsigned int n_edges = 0;
// count how many half edge to allocate
for(fi = m.face.begin(); fi != m.face.end(); ++fi) if(! (*fi).IsD())
{n_edges+=(*fi).VN();
for(int i = 0; i < (*fi).VN(); ++i)
if(vcg::face::IsBorder<FaceType>((*fi),(i)))
++n_edges;
}
m.hedge.clear();
m.hn = 0;
// allocate the half edges
typename MeshType::HEdgeIterator ei = vcg::tri::Allocator<MeshType>::AddHEdges(m,n_edges);
std::vector<VertexPairEdgePtr> all;
int firstEdge = 0;
for(fi = m.face.begin(); fi != m.face.end(); ++fi)if(!(*fi).IsD()){
assert((*fi).VN()>2);
if(flagVisited[*fi].empty()) {flagVisited[*fi].resize((*fi).VN());}
for(int i = 0; i < (*fi).VN(); ++i,++ei)
{
(*ei).HVp() = (*fi).V(i);
(*ei).HNp() = &m.hedge[firstEdge + (i +1) % (*fi).VN()];
if(MeshType::HEdgeType::HasHFAdjacency())
(*ei).HFp() = &(*fi);
if( MeshType::FaceType::HasFHAdjacency())
(*fi).FHp() = &(*ei);
if(MeshType::HEdgeType::HasHPrevAdjacency())
(*ei).HPp() = &m.hedge[firstEdge + (i +(*fi).VN()-1) % (*fi).VN()];
if(HasVHAdjacency(m))
(*ei).HVp()->VHp() = &(*ei);
all.push_back(VertexPairEdgePtr((*fi).V(i), (*fi).V((*fi).Next(i)),&(*ei)));// it will be used to link the hedges
if( vcg::face::IsBorder<FaceType>((*fi),(i)))
borderEdges.push_back(FacePtrInt(&(*fi),i));
}
firstEdge += (*fi).VN();
}
// add all the border hedges
int borderLength;
typename std::vector<FacePtrInt >::iterator ebi;
for( ebi = borderEdges.begin(); ebi != borderEdges.end(); ++ebi)
if( !flagVisited[(*ebi).f][(*ebi).i])// not already inserted
{
borderLength = 0;
vcg::face::Pos<FaceType> bp((*ebi).f,(*ebi).i);
//FaceType * start = (*ebi).f;
VertexType * start = ((*ebi).f)->V((*ebi).i);
do{
all.push_back( VertexPairEdgePtr ( bp.f->V( bp.f->Next(bp.z) ),bp.f->V( bp.z ),&(*ei)));
(*ei).HVp() = bp.f->V(bp.f->Next(bp.z)) ;
flagVisited[bp.f][bp.z] = true;
++ei;
bp.NextB();
++borderLength;
}while (bp.v != start);
//}while (bp.f != start);
// run over the border edges to link the adjacencies
for(int be = 0; be < borderLength; ++be)
{
if(MeshType::HEdgeType::HasHFAdjacency())
m.hedge[firstEdge + be].HFp() = NULL;
if(MeshType::HEdgeType::HasHPrevAdjacency())
m.hedge[firstEdge + be].HPp() = &m.hedge[firstEdge + (be +borderLength-1) % borderLength];
m.hedge[firstEdge + be].HNp() = &m.hedge[firstEdge + (be +1) % borderLength];
}
firstEdge+=borderLength;
}
vcg::tri::Allocator<MeshType>:: template DeletePerFaceAttribute<BitVector>(m,flagVisited );
std::sort(all.begin(),all.end());
assert(all.size() == n_edges);
for(unsigned int i = 0 ; i < all.size(); )
if(all[i] == all[i+1])
{
all[i].ep->HOp() = all[i+1].ep;
all[i+1].ep->HOp() = all[i].ep;
i+=2;
}
else
{
all[i].ep->HOp() = all[i].ep;
i+=1;
}
if(HasEHAdjacency(m) && HasHEAdjacency(m))
{
assert(m.edge.size() == 0 || m.edge.size() == n_edges/2);
if ( m.edge.size() == 0 )
{
m.en = 0;
// allocate the edges
typename MeshType::EdgeIterator edge_i = vcg::tri::Allocator<MeshType>::AddEdges(m,n_edges/2);
for(ei = m.hedge.begin(); ei != m.hedge.end(); ++ei)
{
if((*ei).HEp() == NULL)
{
(*ei).HEp() = &(*edge_i);
(*ei).HOp()->HEp() = &(*edge_i);
(*edge_i).EHp() = &(*ei);
++edge_i;
}
}
}
else
{
if(HasEVAdjacency(m) && HasHEAdjacency(m) && HasEHAdjacency(m))
{
//update edge relations
for(typename MeshType::EdgeIterator ei1 = m.edge.begin(); ei1 != m.edge.end(); ++ei1 )
{
vector<HEdgePointer> hedges = HalfEdgeTopology<MeshType>::get_incident_hedges((*ei1).V(0));
for(typename vector<HEdgePointer>::iterator hi = hedges.begin(); hi != hedges.end(); ++hi)
{
if((*hi)->HOp()->HVp() == (*ei1).V(1))
{
assert((*hi)->HEp() == NULL);
assert((*hi)->HOp()->HEp() == NULL);
// EH
(*ei1).EHp() = *hi;
// HE
(*hi)->HEp() = &(*ei1);
(*hi)->HOp()->HEp() = &(*ei1);
break;
}
}
}
}
}
}
}
/**
Checks pointers FHEp() are valid
**/
static bool CheckConsistency_FHp(MeshType & m){
assert(MeshType::FaceType::HasFHAdjacency());
FaceIterator fi;
for(fi = m.face.begin(); fi != m.face.end(); ++fi)
if(!(*fi).IsD()){
if((*fi).FHp() < &(*m.hedge.begin())) return false;
if((*fi).FHp() > &(m.hedge.back())) return false;
}
return true;
}
/**
Checks that half edges and face relation are consistent
**/
static bool CheckConsistency(MeshType & m){
assert(MeshType::HEdgeType::HasHNextAdjacency());
assert(MeshType::HEdgeType::HasHOppAdjacency());
assert(MeshType::HEdgeType::HasHVAdjacency());
assert(MeshType::FaceType::HasFHAdjacency());
//bool hasHEF = ( MeshType::HEdgeType::HasHFAdjacency());
bool hasHP = ( MeshType::HEdgeType::HasHPrevAdjacency());
FaceIterator fi;
HEdgePointer ep,ep1;
int cnt = 0;
if( MeshType::HEdgeType::HasHFAdjacency() )
{
int iDb = 0;
for(fi = m.face.begin(); fi != m.face.end(); ++fi,++iDb)
if(!(*fi).IsD())
{
ep = ep1 = (*fi).FHp();
do{
if(ep->IsD())
return false; // the hedge should not be connected, it has been deleted
if( ! ep->HFp())
return false;
if(ep->HFp() != &(*fi))
return false;// hedge is not pointing to the rigth face
ep = ep->HNp();
if(cnt++ > m.hn)
return false; // hedges are ill connected (HENp())
}while(ep!=ep1);
}
}
HEdgePointer epPrev;
HEdgeIterator hi;
//bool extEdge ;
for( hi = m.hedge.begin(); hi != m.hedge.end(); ++hi)
if(!(*hi).IsD())
{
//cnt = 0;
epPrev = ep = ep1 = &(*hi);
//do{
//extEdge = (ep->HFp()==NULL);
if(hasHP)
{
if( !ep->HPp())
return false;
if( ep->HPp() == ep)
return false; // the previous of an edge cannot be the edge itself
if( ep->HNp()->HPp() != ep)
return false; // next and prev relation are not mutual
if( ep->HPp()->IsD())
return false; //
}
if( ! ep->HOp() )
return false;
if( ep->HOp() == ep)
return false; // opposite relation is not mutual
if( ep->HOp()->IsD())
return false;
if( ep->HOp()->HOp() != ep)
return false; // opposite relation is not mutual
if( HasHFAdjacency(m) )
{
if(ep->HFp())
{
if( ep->HFp()->IsD())
return false; // pointed face must not be deleted
}
}
if( HasHEAdjacency(m) && (m.en!=0))
{
if( ! ep->HEp())
return false; //halfedge must point to an edge
if( ep->HEp()->IsD())
return false; // pointed edge must not be deleted
if(ep->HEp() != ep->HOp()->HEp())
return false; // he and opposite he must point to the same edge
if(ep->HEp()->EHp() != ep && ep->HEp()->EHp() != ep->HOp() )
return false; // halfedge points to an edge not pointing it or its opposite
}
if( !ep->HNp() )
return false;
if( ep->HNp() == ep )
return false; // the next of an hedge cannot be the hedge itself
if( ep->HNp()->IsD())
return false; //
if(hasHP)
if( ep->HNp()->HPp() != ep)
return false; //
if( HasHVAdjacency(m) )
{
if( ! ep->HVp() )
return false; // halfedge must point to a vertex
if( ep->HVp()->IsD() )
return false; // pointed vertex must not be deleted
if( HasVHAdjacency(m) )
if( ! (ep->HVp()->VHp()) )
return false; // halfedge points to a vertex pointing NULL
}
ep = ep->HNp();
if( ep->HVp() != epPrev->HOp()->HVp())
return false;
epPrev = ep;
// if(cnt++ > m.hn)
// return false; // edges are ill connected (HENp())
//}while(ep!=ep1);
}
if(HasEHAdjacency(m) && HasHEAdjacency(m))
for(EdgeIterator ei = m.edge.begin(); ei != m.edge.end(); ++ei)
{
if(!(*ei).IsD())
{
if( !(*ei).EHp())
return false; //edge must have a valid pointer to his halfedge
if( (*ei).EHp()->HEp() != &(*ei) )
return false; // edge's halfedge must point to the edge itself
if( (*ei).EHp()->IsD())
return false;
}
}
if(HasVHAdjacency(m))
for(VertexIterator vi = m.vert.begin(); vi != m.vert.end(); ++vi)
{
if( !(*vi).IsD() )
if( (*vi).VHp() )
{
if( (*vi).VHp()->HVp() != &(*vi) )
return false;
if( (*vi).VHp()->IsD())
return false;
}
}
return true;
}
/** Set the relations HFp(), FHp() from a loop of edges to a face
*/
private:
static void SetRelationsLoopFace(HEdgeType * e0, FaceType * f){
assert(HEdgeType::HasHNextAdjacency());
assert(FaceType::HasFHAdjacency());
HEdgeType *e = e0;
assert(e!=NULL);
do{ e->HFp() = f; e = e->HNp(); } while(e != e0);
f->FHp() = e0;
}
/**
Merge the two faces. This will probably become a class template or a functor
*/
static void MergeFaces(FaceType *, FaceType *){}
/**
Find previous hedge in the loop
*/
static HEdgeType * PreviousEdge(HEdgeType * e0){
HEdgeType * ep = e0;
do{
if(ep->HNp() == e0) return ep;
ep = ep->HNp();
}while(ep!=e0);
assert(0); // degenerate loop
return 0;
}
public:
/** Adds an edge between the sources of e0 and e1 and set all the topology relations.
If the edges store the pointers to the faces then a new face is created.
<--- e1 ---- X <------e1_HEPp---
^
||
ei0 || ei1
||
v
----e0_HEPp-> X ----- e0 ------>
*/
static void AddHEdge(MeshType &m, HEdgeType * e0, HEdgeType * e1){
HEdgeType *iii =e0->HNp();
assert(e1!=e0->HNp());
assert(e0!=e1->HNp());
HEdgePointer tmp;
bool hasP = MeshType::HEdgeType::HasHPrevAdjacency();
assert(e0->HOp() != e1); // the hedge already exists
assert(e0!=e1->HNp());
std::vector<typename MeshType::HEdgePointer* > toUpdate;
toUpdate.push_back(&e0);
toUpdate.push_back(&e1);
HEdgeIterator ei0 = vcg::tri::Allocator<MeshType>::AddHEdges(m,2,toUpdate);
HEdgeIterator ei1 = ei0; ++ei1;
(*ei0).HNp() = e1;(*ei0).HVp() = e0->HVp();
(*ei1).HNp() = e0;(*ei1).HVp() = e1->HVp();
HEdgePointer e0_HEPp = 0,e1_HEPp = 0,ep =0;
if(hasP){
e0_HEPp = e0->HPp();
e1_HEPp = e1->HPp();
}else{// does not have pointer to previous, it must be computed
ep = e0;
do{
if(ep->HNp() == e0) e0_HEPp = ep;
if(ep->HNp() == e1) e1_HEPp = ep;
ep = ep->HNp();
}while(ep!=e0);
}
if(hasP){
(*ei0).HPp() = e0->HPp();
(*ei1).HPp() = e1->HPp();
e0->HPp() = &(*ei1);
e1->HPp() = &(*ei0);
}
e0_HEPp -> HNp() = &(*ei0);
e1_HEPp -> HNp() = &(*ei1);
(*ei0).HOp() = &(*ei1);
(*ei1).HOp() = &(*ei0);
if( HEdgeType::HasHFAdjacency() && FaceType::HasFHAdjacency()){
FaceIterator fi0 = vcg::tri::Allocator<MeshType>::AddFaces(m,1);
m.face.back().ImportData(*e0->HFp());
SetRelationsLoopFace(&(*ei0),e1->HFp()); // one loop to the old face
SetRelationsLoopFace(&(*ei1),&m.face.back()); // the other to the new face
}
}
/** Detach the topology relations of a given edge
<--- e->HENPp -X --- <---------eO_HEPp---
^
||
e || e->HEOp()
||
v
----e_HEPp--> X ----- e->HEOp->HENPp() ------>
*/
static void RemoveHEdge(MeshType &m, HEdgeType * e){
assert(MeshType::HEdgeType::HasHNextAdjacency());
assert(MeshType::HEdgeType::HasHOppAdjacency());
assert(MeshType::FaceType::HasFHAdjacency());
bool hasP = MeshType::HEdgeType::HasHPrevAdjacency();
HEdgePointer e_HEPp,eO_HEPp;
if(hasP){
e_HEPp = e->HPp();
eO_HEPp = e->HOp()->HPp();
}else{
e_HEPp = PreviousEdge(e);
eO_HEPp = PreviousEdge(e->HOp());
}
assert(e_HEPp->HNp() == e);
assert(eO_HEPp->HNp() == e->HOp());
e_HEPp->HNp() = e->HOp()->HNp();
eO_HEPp->HNp() = e-> HNp();
if(hasP) {
e->HOp()->HNp()->HPp() = e_HEPp;
e->HNp()->HPp() = eO_HEPp;
e->HPp() = NULL;
e-> HOp()->HPp() = NULL;
}
// take care of the faces
if(MeshType::HEdgeType::HasHFAdjacency()){
MergeFaces(e_HEPp->HFp(),eO_HEPp->HFp());
vcg::tri::Allocator<MeshType>::DeleteFace(m,*eO_HEPp->HFp());
SetRelationsLoopFace(e_HEPp,e_HEPp->HFp());
}
vcg::tri::Allocator<MeshType>::DeleteHEdge(m,*e->HOp());
vcg::tri::Allocator<MeshType>::DeleteHEdge(m,*e);
}
};// end class
template <class MeshType >
struct UpdateIndexed{
typedef typename MeshType::VertexType VertexType;
typedef typename MeshType::VertexPointer VertexPointer;
typedef typename MeshType::HEdgePointer HEdgePointer;
typedef typename MeshType::HEdgeType HEdgeType;
typedef typename MeshType::HEdgeIterator HEdgeIterator;
typedef typename MeshType::FaceIterator FaceIterator;
typedef typename MeshType::FaceType FaceType;
struct VertexPairEdgePtr{
VertexPairEdgePtr(VertexPointer _v0,VertexPointer _v1,HEdgePointer _ep):v0(_v0),v1(_v1),ep(_ep){if(v0>v1) std::swap(v0,v1);}
const bool operator <(const VertexPairEdgePtr &o) const {return (v0 == o.v0)? (v1<o.v1):(v0<o.v0);}
const bool operator ==(const VertexPairEdgePtr &o) const {return (v0 == o.v0)&& (v1==o.v1);}
VertexPointer v0,v1;
HEdgePointer ep;
};
/**
builds an indexed data structure from a half-edge data structure.
Note: if the half edge have the pointer to face
their relation FV (face-vertex) will be computed and the data possibly stored in the
face will be preserved.
**/
static void FromHalfEdges( MeshType & m ){
assert(HasFVAdjacency(m));
assert(MeshType::HEdgeType::HasHNextAdjacency());
assert(MeshType::HEdgeType::HasHVAdjacency());
assert(MeshType::HEdgeType::HasHOppAdjacency());
assert(MeshType::FaceType::HasFHAdjacency());
bool hasHEF;
//bool createFace,hasHEF,hasFHE;
// typename MeshType::template PerHEdgeAttributeHandle<bool> hV = Allocator<MeshType>::template AddPerHEdgeAttribute<bool>(m,"");
typename MeshType::HEdgeIterator ei;
typename MeshType::FacePointer fp;
typename MeshType::FaceIterator fi;
typename MeshType::HEdgePointer ep,epF;
//int vi = 0;
vcg::SimpleTempData<typename MeshType::HEdgeContainer,bool> hV(m.hedge);
hasHEF = (MeshType::HEdgeType::HasHFAdjacency());
assert( !hasHEF || (hasHEF && m.fn>0));
// if the edgetype has the pointer to face
// it is assumed the the edget2face pointer (HEFp) are correct
// and the faces are allocated
for ( ei = m.hedge.begin(); ei != m.hedge.end(); ++ei)
if(!(*ei).IsD()) // it has not been deleted
if(!hasHEF || ( hasHEF && (*ei).HFp()!=NULL)) // if it has a pointer to the face it is
// not null (i.e. it is not a border edge)
if(!hV[(*ei)] ) // it has not be visited yet
{
if(!hasHEF)// if it has
fp = &(* Allocator<MeshType>::AddFaces(m,1));
else
fp = (*ei).HFp();
ep = epF = &(*ei);
std::vector<VertexPointer> vpts;
do{vpts.push_back((*ep).HVp()); ep=ep->HNp();}while(ep!=epF);
//int idbg =fp->VN();
if(fp->VN() != vpts.size()){
fp->Dealloc();
fp ->Alloc(vpts.size());
}
//int idbg1 =fp->VN();
for(unsigned int i = 0; i < vpts.size();++i) fp ->V(i) = vpts[i];// set the pointer from face to vertex
hV[(*ei)] = true;
}
//Allocator<MeshType>::DeletePerHEdgeAttribute(m,hV);
}
};
} // end namespace vcg
}
#endif // __VCGLIB_EDGE_SUPPORT

File diff suppressed because it is too large Load Diff

View File

@ -1,398 +0,0 @@
/****************************************************************************
* 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. *
* *
****************************************************************************/
#ifndef __VCG_TRI_UPDATE_NORMALS
#define __VCG_TRI_UPDATE_NORMALS
#include <vcg/space/triangle3.h>
#include <vcg/math/matrix33.h>
#include <vcg/complex/trimesh/update/flag.h>
namespace vcg {
namespace tri {
/// \ingroup trimesh
/// \headerfile normal.h vcg/complex/trimesh/update/normal.h
/// \brief 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 UpdateNormals
{
public:
typedef ComputeMeshType MeshType;
typedef typename MeshType::VertexType VertexType;
typedef typename MeshType::CoordType CoordType;
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;
/**
Set to zero all the normals. Usued by all the face averaging algorithms.
by default it does not clear the normals of unreferenced vertices because they could be still useful
*/
static void PerVertexClear(ComputeMeshType &m, bool ClearAllVertNormal=false)
{
assert(HasPerVertexNormal(m));
if(ClearAllVertNormal)
UpdateFlags<ComputeMeshType>::VertexClearV(m);
else
{
UpdateFlags<ComputeMeshType>::VertexSetV(m);
for(FaceIterator f=m.face.begin();f!=m.face.end();++f)
if( !(*f).IsD() )
for(int i=0;i<3;++i) (*f).V(i)->ClearV();
}
VertexIterator vi;
for(vi=m.vert.begin();vi!=m.vert.end();++vi)
if( !(*vi).IsD() && (*vi).IsRW() && (!(*vi).IsV()) )
(*vi).N() = NormalType((ScalarType)0,(ScalarType)0,(ScalarType)0);
}
/// \brief Calculates the face normal (if stored in the current face type)
static void PerFace(ComputeMeshType &m)
{
if( !m.HasPerFaceNormal()) return;
FaceIterator f;
for(f=m.face.begin();f!=m.face.end();++f)
if( !(*f).IsD() ) face::ComputeNormal(*f);
}
/// \brief Calculates the vertex normal. Exploiting or current face normals.
/**
The normal of a vertex v is the weigthed average of the normals of the faces incident on v.
*/
static void PerVertexFromCurrentFaceNormal(ComputeMeshType &m)
{
if( !m.HasPerVertexNormal()) return;
VertexIterator vi;
for(vi=m.vert.begin();vi!=m.vert.end();++vi)
if( !(*vi).IsD() && (*vi).IsRW() )
(*vi).N()=CoordType(0,0,0);
FaceIterator fi;
for(fi=m.face.begin();fi!=m.face.end();++fi)
if( !(*fi).IsD())
{
for(int j=0; j<3; ++j)
if( !(*fi).V(j)->IsD())
(*fi).V(j)->N() += (*fi).cN();
}
}
/// \brief Calculates the vertex normal. Exploiting or current face normals.
/**
The normal of a face f is the average of the normals of the vertices of f.
*/
static void PerFaceFromCurrentVertexNormal(ComputeMeshType &m)
{
for (FaceIterator fi=m.face.begin(); fi!=m.face.end(); ++fi)
if( !(*fi).IsD())
{
NormalType n;
n.SetZero();
for(int j=0; j<3; ++j)
n += fi->V(j)->cN();
n.Normalize();
fi->N() = n;
}
}
/// \brief Calculates the vertex normal. Without exploiting or touching face normals.
/**
The normal of a vertex v computed as a weighted sum f the incident face normals.
The weight is simlply the angle of the involved wedge. Described in:
G. Thurmer, C. A. Wuthrich
"Computing vertex normals from polygonal facets"
Journal of Graphics Tools, 1998
*/
static void PerVertexAngleWeighted(ComputeMeshType &m)
{
assert(HasPerVertexNormal(m));
PerVertexClear(m);
FaceIterator f;
for(f=m.face.begin();f!=m.face.end();++f)
if( !(*f).IsD() && (*f).IsR() )
{
typename FaceType::NormalType t = vcg::NormalizedNormal(*f);
NormalType e0 = ((*f).V1(0)->cP()-(*f).V0(0)->cP()).Normalize();
NormalType e1 = ((*f).V1(1)->cP()-(*f).V0(1)->cP()).Normalize();
NormalType e2 = ((*f).V1(2)->cP()-(*f).V0(2)->cP()).Normalize();
(*f).V(0)->N() += t*AngleN(e0,-e2);
(*f).V(1)->N() += t*AngleN(-e0,e1);
(*f).V(2)->N() += t*AngleN(-e1,e2);
}
}
/// \brief Calculates the vertex normal. Without exploiting or touching face normals.
/**
The normal of a vertex v is computed according to the formula described by Nelson Max in
Max, N., "Weights for Computing Vertex Normals from Facet Normals", Journal of Graphics Tools, 4(2) (1999)
The weight for each wedge is the cross product of the two edge over the product of the square of the two edge lengths.
According to the original paper it is perfect only for spherical surface, but it should perform well...
*/
static void PerVertexWeighted(ComputeMeshType &m)
{
assert(HasPerVertexNormal(m));
PerVertexClear(m);
FaceIterator f;
for(f=m.face.begin();f!=m.face.end();++f)
if( !(*f).IsD() && (*f).IsR() )
{
typename FaceType::NormalType t = vcg::Normal(*f);
ScalarType e0 = SquaredDistance((*f).V0(0)->cP(),(*f).V1(0)->cP());
ScalarType e1 = SquaredDistance((*f).V0(1)->cP(),(*f).V1(1)->cP());
ScalarType e2 = SquaredDistance((*f).V0(2)->cP(),(*f).V1(2)->cP());
(*f).V(0)->N() += t/(e0*e2);
(*f).V(1)->N() += t/(e0*e1);
(*f).V(2)->N() += t/(e1*e2);
}
}
/// \brief Calculates the vertex normal. Without exploiting or touching face normals.
/**
The normal of a vertex v is the classical area weigthed average of the normals of the faces incident on v.
*/
static void PerVertex(ComputeMeshType &m)
{
assert(HasPerVertexNormal(m));
PerVertexClear(m);
FaceIterator f;
for(f=m.face.begin();f!=m.face.end();++f)
if( !(*f).IsD() && (*f).IsR() )
{
//typename FaceType::NormalType t = (*f).Normal();
typename FaceType::NormalType t = vcg::Normal(*f);
for(int j=0; j<3; ++j)
if( !(*f).V(j)->IsD() && (*f).V(j)->IsRW() )
(*f).V(j)->N() += t;
}
}
/// \brief Calculates both vertex and face normals.
/**
The normal of a vertex v is the weigthed average of the normals of the faces incident on v.
*/
static void PerVertexPerFace(ComputeMeshType &m)
{
if( !m.HasPerVertexNormal() || !m.HasPerFaceNormal()) return;
PerFace(m);
PerVertexClear(m);
FaceIterator f;
for(f=m.face.begin();f!=m.face.end();++f)
if( !(*f).IsD() && (*f).IsR() )
{
for(int j=0; j<3; ++j)
if( !(*f).V(j)->IsD() && (*f).V(j)->IsRW() )
(*f).V(j)->N() += (*f).cN();
}
}
/// \brief Calculates both vertex and face normals.
/**
The normal of a vertex v is the weigthed average of the normals of the faces incident on v.
*/
static void PerVertexNormalizedPerFace(ComputeMeshType &m)
{
PerVertexPerFace(m);
NormalizeVertex(m);
}
/// \brief Normalize the lenght of the face normals.
static void NormalizeVertex(ComputeMeshType &m)
{
VertexIterator vi;
for(vi=m.vert.begin();vi!=m.vert.end();++vi)
if( !(*vi).IsD() && (*vi).IsRW() )
(*vi).N().Normalize();
}
/// \brief Normalize the lenght of the face normals.
static void NormalizeFace(ComputeMeshType &m)
{
FaceIterator fi;
for(fi=m.face.begin();fi!=m.face.end();++fi)
if( !(*fi).IsD() ) (*fi).N().Normalize();
}
static void AreaNormalizeFace(ComputeMeshType &m)
{
FaceIterator fi;
for(fi=m.face.begin();fi!=m.face.end();++fi)
if( !(*fi).IsD() )
{
(*fi).N().Normalize();
(*fi).N() = (*fi).N() * DoubleArea(*fi);
}
}
static void PerVertexNormalizedPerFaceNormalized(ComputeMeshType &m)
{
PerVertexNormalizedPerFace(m);
NormalizeFace(m);
}
static void PerFaceRW(ComputeMeshType &m, bool normalize=false)
{
if( !m.HasPerFaceNormal()) return;
FaceIterator f;
bool cn = true;
if(normalize)
{
for(f=m.m.face.begin();f!=m.m.face.end();++f)
if( !(*f).IsD() && (*f).IsRW() )
{
for(int j=0; j<3; ++j)
if( !(*f).V(j)->IsR()) cn = false;
if( cn ) face::ComputeNormalizedNormal(*f);
cn = true;
}
}
else
{
for(f=m.m.face.begin();f!=m.m.face.end();++f)
if( !(*f).IsD() && (*f).IsRW() )
{
for(int j=0; j<3; ++j)
if( !(*f).V(j)->IsR()) cn = false;
if( cn )
(*f).ComputeNormal();
cn = true;
}
}
}
static void PerFaceNormalized(ComputeMeshType &m)
{
if( !m.HasPerFaceNormal()) return;
FaceIterator f;
for(f=m.face.begin();f!=m.face.end();++f)
if( !(*f).IsD() ) face::ComputeNormalizedNormal(*f);
}
static void PerBitQuadFaceNormalized(ComputeMeshType &m)
{
if( !m.HasPerFaceNormal()) return;
PerFace(m);
FaceIterator f;
for(f=m.face.begin();f!=m.face.end();++f) {
if( !(*f).IsD() ) {
for (int k=0; k<3; k++) if (f->IsF(k))
if (&*f < f->FFp(k)) {
f->N() = f->FFp(k)->N() = (f->FFp(k)->N() + f->N()).Normalize();
}
}
}
}
/// \brief Calculates the vertex normal.
static void PerVertexNormalized(ComputeMeshType &m)
{
if( !m.HasPerVertexNormal()) return;
PerVertex(m);
for(VertexIterator vi=m.vert.begin();vi!=m.vert.end();++vi)
if( !(*vi).IsD() && (*vi).IsRW() )
(*vi).N().Normalize();
}
/// \brief Multiply the vertex normals by the matrix passed. By default, the scale component is removed.
static void PerVertexMatrix(ComputeMeshType &m, const Matrix44<ScalarType> &mat, bool remove_scaling= true){
float scale;
Matrix33<ScalarType> mat33(mat,3);
if( !m.HasPerVertexNormal()) return;
if(remove_scaling){
scale = pow(mat33.Determinant(),(ScalarType)(1.0/3.0));
mat33[0][0]/=scale;
mat33[1][1]/=scale;
mat33[2][2]/=scale;
}
for(VertexIterator vi=m.vert.begin();vi!=m.vert.end();++vi)
if( !(*vi).IsD() && (*vi).IsRW() )
(*vi).N() = mat33*(*vi).N();
}
/// \brief Multiply the face normals by the matrix passed. By default, the scale component is removed.
static void PerFaceMatrix(ComputeMeshType &m, const Matrix44<ScalarType> &mat, bool remove_scaling= true){
float scale;
Matrix33<ScalarType> mat33(mat,3);
if( !m.HasPerFaceNormal()) return;
if(remove_scaling){
scale = pow(mat33.Determinant(),ScalarType(1.0/3.0));
mat33[0][0]/=scale;
mat33[1][1]/=scale;
mat33[2][2]/=scale;
}
for(FaceIterator fi=m.face.begin();fi!=m.face.end();++fi)
if( !(*fi).IsD() && (*fi).IsRW() )
(*fi).N() = mat33* (*fi).N();
}
}; // end class
} // End namespace
} // End namespace
#endif

View File

@ -1,83 +0,0 @@
/****************************************************************************
* 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.1 2005/07/06 08:02:27 cignoni
Initial commit
****************************************************************************/
#ifndef __VCG_TRI_UPDATE_POSITION
#define __VCG_TRI_UPDATE_POSITION
#include "normal.h"
namespace vcg {
namespace tri {
/// \ingroup trimesh
/// \headerfile position.h vcg/complex/trimesh/update/position.h
/// \brief This class is used to update vertex position according to a transformation matrix.
template <class ComputeMeshType>
class UpdatePosition
{
public:
typedef ComputeMeshType MeshType;
typedef typename MeshType::ScalarType ScalarType;
typedef typename MeshType::VertexType VertexType;
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;
/// \brief Multiply
static void Matrix(ComputeMeshType &m, const Matrix44<ScalarType> &M, bool update_also_normals = true)
{
VertexIterator vi;
for(vi=m.vert.begin();vi!=m.vert.end();++vi)
if(!(*vi).IsD()) (*vi).P()=M*(*vi).cP();
if(update_also_normals){
if(m.HasPerVertexNormal()){
UpdateNormals<ComputeMeshType>::PerVertexMatrix(m,M);
}
if(m.HasPerFaceNormal()){
UpdateNormals<ComputeMeshType>::PerFaceMatrix(m,M);
}
}
}
}; // end class
} // End namespace
} // End namespace
#endif

View File

@ -1,351 +0,0 @@
/****************************************************************************
* 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. *
* *
****************************************************************************/
#ifndef __VCG_TRI_UPDATE_QUALITY
#define __VCG_TRI_UPDATE_QUALITY
#include <vcg/simplex/face/pos.h>
#include <vcg/simplex/face/topology.h>
#include <vcg/complex/trimesh/update/flag.h>
#include <vcg/complex/trimesh/stat.h>
#include <algorithm>
#include <vector>
#include <stack>
#include <assert.h>
namespace vcg {
namespace tri {
/// \ingroup trimesh
/// \headerfile quality.h vcg/complex/trimesh/update/quality.h
/// \brief Generation of per-vertex and per-face qualities.
/**
It works according to various strategy, like geodesic distance from the border (UpdateQuality::VertexGeodesicFromBorder) or curvature ecc.
This class is templated over the mesh and (like all other Update* classes) has only static members; Typical usage:
\code
MyMeshType m;
UpdateQuality<MyMeshType>::VertexGeodesicFromBorder(m);
\endcode
*/
template <class UpdateMeshType>
class UpdateQuality
{
public:
typedef UpdateMeshType MeshType;
typedef typename MeshType::ScalarType ScalarType;
typedef typename MeshType::VertexType VertexType;
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;
class VQualityHeap
{
public:
float q;
VertexPointer p;
inline VQualityHeap( VertexPointer np )
{
q = np->Q();
p = np;
}
// Attenzione il minore e' maggiore
inline bool operator < ( const VQualityHeap & vq ) const { return q > vq.q; }
inline bool operator == ( const VQualityHeap & vq ) const { return q == vq.q; }
inline bool operator > ( const VQualityHeap & vq ) const { return q < vq.q; }
inline bool operator != ( const VQualityHeap & vq ) const { return q != vq.q; }
inline bool operator <= ( const VQualityHeap & vq ) const { return q >= vq.q; }
inline bool operator >= ( const VQualityHeap & vq ) const { return q <= vq.q; }
inline bool is_valid() const { return q==p->Q(); }
};
// *** IMPORTANT REQUIREMENTS
// VF topology
// Border FLags
// tri::UpdateTopology<SMesh>::VertexFace(sm);
// tri::UpdateFlags<SMesh>::FaceBorderFromVF(sm);
//
// Calcola la qualita' come distanza geodesica dal bordo della mesh.
// Robusta funziona anche per mesh non manifold.
// La qualita' memorizzata indica la distanza assoluta dal bordo della mesh.
// Nota prima del 13/11/03 in alcuni casi rari SPT andava in loop perche' poteva capitare
// che per approx numeriche ben strane pw->Q() > pv->Q()+d ma durante la memorizzazione
// della nuova distanza essa rimanesse uguale a prima. Patchato rimettendo i vertici nello
// heap solo se migliorano la distanza di un epsilon == 1/100000 della mesh diag.
/// \brief Compute, for each vertex of the mesh the geodesic distance from the border of the mesh itself.
/**
It uses the classical Dijkstra Shortest Path Tree algorithm.
The geodesic distance is approximated by allowing to walk only along edges of the mesh.
\warning VF topology, Per Vertex Quality and border flags already computed (see UpdateFlags::FaceBorderFromVF and UpdateTopology::VertexFace);
*/
static void VertexGeodesicFromBorder(MeshType &m) // R1
{
//Requirements
assert(m.HasVFTopology());
assert(m.HasPerVertexQuality());
std::vector< VQualityHeap > heap;
VertexIterator v;
FaceIterator f;
int j;
for(v=m.vert.begin();v!=m.vert.end();++v)
(*v).Q() = -1;
for(f=m.face.begin();f!=m.face.end();++f) // Inserisco nell'heap i v di bordo
if(!(*f).IsD())
for(j=0;j<3;++j)
if( (*f).IsB(j) )
{
for(int k=0;k<2;++k)
{
VertexPointer pv = (*f).V((j+k)%3);
if( pv->Q()==-1 )
{
pv->Q() = 0;
heap.push_back(VQualityHeap(pv));
}
}
}
const ScalarType loc_eps=m.bbox.Diag()/ScalarType(100000);
while( heap.size()!=0 ) // Shortest path tree
{
VertexPointer pv;
std::pop_heap(heap.begin(),heap.end());
if( ! heap.back().is_valid() )
{
heap.pop_back();
continue;
}
pv = heap.back().p;
heap.pop_back();
for(face::VFIterator<FaceType> vfi(pv) ; !vfi.End(); ++vfi )
{
for(int k=0;k<2;++k)
{
VertexPointer pw;
float d;
if(k==0) pw = vfi.f->V1(vfi.z);
else pw = vfi.f->V2(vfi.z);
d = Distance(pv->P(),pw->P());
if( pw->Q()==-1 || pw->Q() > pv->Q()+d + loc_eps)
{
pw->Q() = pv->Q()+d;
heap.push_back(VQualityHeap(pw));
std::push_heap(heap.begin(),heap.end());
}
}
}
}
for(v=m.vert.begin();v!=m.vert.end();++v)
if(v->Q()==-1)
v->Q() = 0;
}
/** Assign to each vertex of the mesh a constant quality value. Useful for initialization.
*/
static void VertexConstant(MeshType &m, float q)
{
VertexIterator vi;
for(vi=m.vert.begin();vi!=m.vert.end();++vi) if(!(*vi).IsD())
(*vi).Q()=q;
}
/** Clamp each vertex of the mesh with a range of values.
*/
static void VertexClamp(MeshType &m, float qmin, float qmax)
{
VertexIterator vi;
for(vi=m.vert.begin();vi!=m.vert.end();++vi) if(!(*vi).IsD())
(*vi).Q()=std::min(qmax, std::max(qmin,(*vi).Q()));
}
/** Normalize the vertex quality so that it fits in the specified range.
*/
static void VertexNormalize(MeshType &m, float qmin=0.0, float qmax=1.0)
{
ScalarType deltaRange = qmax-qmin;
std::pair<ScalarType,ScalarType> minmax = tri::Stat<MeshType>::ComputePerVertexQualityMinMax(m);
VertexIterator vi;
for(vi = m.vert.begin(); vi != m.vert.end(); ++vi)
(*vi).Q() = qmin+deltaRange*((*vi).Q() - minmax.first)/(minmax.second - minmax.first);
}
/** Normalize the face quality so that it fits in the specified range.
*/
static void FaceNormalize(MeshType &m, float qmin=0.0, float qmax=1.0)
{
ScalarType deltaRange = qmax-qmin;
std::pair<ScalarType,ScalarType> minmax = tri::Stat<MeshType>::ComputePerFaceQualityMinMax(m);
FaceIterator fi;
for(fi = m.face.begin(); fi != m.face.end(); ++fi)
(*fi).Q() = qmin+deltaRange*((*fi).Q() - minmax.first)/(minmax.second - minmax.first);
}
/** Assign to each face of the mesh a constant quality value. Useful for initialization.
*/
static void FaceConstant(MeshType &m, float q)
{
FaceIterator fi;
for(fi=m.face.begin();fi!=m.face.end();++fi)
(*fi).Q()=q;
}
static void VertexFromGaussianCurvature(MeshType &m)
{
VertexIterator vi;
for(vi=m.vert.begin();vi!=m.vert.end();++vi) if(!(*vi).IsD())
(*vi).Q() = (*vi).Kg();
}
static void VertexFromMeanCurvature(MeshType &m)
{
VertexIterator vi;
for(vi=m.vert.begin();vi!=m.vert.end();++vi) if(!(*vi).IsD())
(*vi).Q() = (*vi).Kh();
}
/*
* Absolute Curvature
*
* 2|H| if K >= 0
* |k1| + |k2| = <
* 2 * sqrt(|H|^2-K) otherwise
*
* defs and formulas taken from
*
* Improved curvature estimation for watershed segmentation of 3-dimensional meshes
* S Pulla, A Razdan, G Farin - Arizona State University, Tech. Rep, 2001
* and from
* Optimizing 3D triangulations using discrete curvature analysis
* N Dyn, K Hormann, SJ Kim, D Levin - Mathematical Methods for Curves and Surfaces: Oslo, 2000
*/
static void VertexFromAbsoluteCurvature(MeshType &m)
{
VertexIterator vi;
for(vi=m.vert.begin();vi!=m.vert.end();++vi) if(!(*vi).IsD())
{
if((*vi).Kg() >= 0)
(*vi).Q() = math::Abs( 2*(*vi).Kh() );
else
(*vi).Q() = 2*math::Sqrt(math::Abs( (*vi).Kh()*(*vi).Kh() - (*vi).Kg()));
}
}
/*
* RMS Curvature = sqrt(4H^2-2K)
* def and formula taken from
*
* Improved curvature estimation for watershed segmentation of 3-dimensional meshes
* S Pulla, A Razdan, G Farin - Arizona State University, Tech. Rep, 2001
*/
static void VertexFromRMSCurvature(MeshType &m)
{
VertexIterator vi;
for(vi=m.vert.begin();vi!=m.vert.end();++vi) if(!(*vi).IsD())
(*vi).Q() = math::Sqrt(math::Abs( 4*(*vi).Kh()*(*vi).Kh() - 2*(*vi).Kg()));
}
/*
Saturate the vertex quality so that for each vertex the gradient of the quality is lower than the given threshold value (in absolute value)
The saturation is done in a conservative way (quality is always decreased and never increased)
Note: requires VF adjacency.
*/
static void VertexSaturate(MeshType &m, ScalarType gradientThr=1.0)
{
UpdateFlags<MeshType>::VertexClearV(m);
std::stack<VertexPointer> st;
st.push(&*m.vert.begin());
while(!st.empty())
{
VertexPointer vc = st.top(); // the center
//printf("Stack size %i\n",st.size());
//printf("Pop elem %i %f\n",st.top() - &*m.vert.begin(), st.top()->Q());
st.pop();
vc->SetV();
std::vector<VertexPointer> star;
typename std::vector<VertexPointer>::iterator vvi;
face::VVStarVF<FaceType>(vc,star);
for(vvi=star.begin();vvi!=star.end();++vvi )
{
float &qi = (*vvi)->Q();
float distGeom = Distance((*vvi)->cP(),vc->cP()) / gradientThr;
// Main test if the quality varies more than the geometric displacement we have to lower something.
if( distGeom < fabs(qi - vc->Q()))
{
// center = 0 other=10 -> other =
// center = 10 other=0
if(vc->Q() > qi) // first case: the center of the star has to be lowered (and re-inserted in the queue).
{
//printf("Reinserting center %i \n",vc - &*m.vert.begin());
vc->Q() = qi+distGeom-0.00001f;
assert( distGeom > fabs(qi - vc->Q()));
st.push(vc);
break;
}
else
{
// second case: you have to lower qi, the vertex under examination.
assert( distGeom < fabs(qi - vc->Q()));
assert(vc->Q() < qi);
float newQi = vc->Q() + distGeom -0.00001f;
assert(newQi <= qi);
assert(vc->Q() < newQi);
assert( distGeom > fabs(newQi - vc->Q()) );
// printf("distGeom %f, qi %f, vc->Q() %f, fabs(qi - vc->Q()) %f\n",distGeom,qi,vc->Q(),fabs(qi - vc->Q()));
qi = newQi;
(*vvi)->ClearV();
}
}
if(!(*vvi)->IsV())
{
st.push( *vvi);
// printf("Reinserting side %i \n",*vvi - &*m.vert.begin());
(*vvi)->SetV();
}
}
}
}
}; //end class
} // end namespace
} // end namespace
#endif

View File

@ -1,446 +0,0 @@
/****************************************************************************
* 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. *
* *
****************************************************************************/
#ifndef __VCG_TRI_UPDATE_SELECTION
#define __VCG_TRI_UPDATE_SELECTION
#include <queue>
#include <vcg/complex/trimesh/update/flag.h>
namespace vcg {
namespace tri {
/// \ingroup trimesh
/// \brief A stack for saving and restoring selection.
/**
This class is used to save the current selection onto a stack for later use.
\todo it should be generalized to other attributes with a templated approach.
*/
template <class ComputeMeshType>
class SelectionStack
{
typedef typename ComputeMeshType::template PerVertexAttributeHandle< bool > vsHandle;
typedef typename ComputeMeshType::template PerFaceAttributeHandle< bool > fsHandle;
public:
SelectionStack(ComputeMeshType &m)
{
_m=&m;
}
bool push()
{
vsHandle vsH = Allocator<ComputeMeshType>::template AddPerVertexAttribute< bool >(*_m);
fsHandle fsH = Allocator<ComputeMeshType>::template AddPerFaceAttribute< bool > (*_m);
typename ComputeMeshType::VertexIterator vi;
for(vi = _m->vert.begin(); vi != _m->vert.end(); ++vi)
if( !(*vi).IsD() ) vsH[*vi] = (*vi).IsS() ;
typename ComputeMeshType::FaceIterator fi;
for(fi = _m->face.begin(); fi != _m->face.end(); ++fi)
if( !(*fi).IsD() ) fsH[*fi] = (*fi).IsS() ;
vsV.push_back(vsH);
fsV.push_back(fsH);
return true;
}
bool pop()
{
if(vsV.empty()) return false;
vsHandle vsH = vsV.back();
fsHandle fsH = fsV.back();
if(! (Allocator<ComputeMeshType>::template IsValidHandle(*_m, vsH))) return false;
typename ComputeMeshType::VertexIterator vi;
for(vi = _m->vert.begin(); vi != _m->vert.end(); ++vi)
if( !(*vi).IsD() )
if(vsH[*vi]) (*vi).SetS() ;
else (*vi).ClearS() ;
typename ComputeMeshType::FaceIterator fi;
for(fi = _m->face.begin(); fi != _m->face.end(); ++fi)
if( !(*fi).IsD() )
if(fsH[*fi]) (*fi).SetS() ;
else (*fi).ClearS() ;
Allocator<ComputeMeshType>::template DeletePerVertexAttribute<bool>(*_m,vsH);
Allocator<ComputeMeshType>::template DeletePerFaceAttribute<bool>(*_m,fsH);
vsV.pop_back();
fsV.pop_back();
return true;
}
private:
ComputeMeshType *_m;
std::vector<vsHandle> vsV;
std::vector<fsHandle> fsV;
};
/// \ingroup trimesh
/// \headerfile selection.h vcg/complex/trimesh/update/selection.h
/// \brief 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 UpdateSelection
{
public:
typedef ComputeMeshType MeshType;
typedef typename MeshType::ScalarType ScalarType;
typedef typename MeshType::VertexType VertexType;
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 vcg::Box3<ScalarType> Box3Type;
static size_t AllVertex(MeshType &m)
{
VertexIterator vi;
for(vi = m.vert.begin(); vi != m.vert.end(); ++vi)
if( !(*vi).IsD() ) (*vi).SetS();
return m.vn;
}
static size_t AllFace(MeshType &m)
{
FaceIterator fi;
for(fi = m.face.begin(); fi != m.face.end(); ++fi)
if( !(*fi).IsD() ) (*fi).SetS();
return m.fn;
}
static size_t ClearVertex(MeshType &m)
{
VertexIterator vi;
for(vi = m.vert.begin(); vi != m.vert.end(); ++vi)
if( !(*vi).IsD() ) (*vi).ClearS();
return 0;
}
static size_t ClearFace(MeshType &m)
{
FaceIterator fi;
for(fi = m.face.begin(); fi != m.face.end(); ++fi)
if( !(*fi).IsD() ) (*fi).ClearS();
return 0;
}
static void Clear(MeshType &m)
{
ClearVertex(m);
ClearFace(m);
}
static size_t CountFace(MeshType &m)
{
size_t selCnt=0;
FaceIterator fi;
for(fi=m.face.begin();fi!=m.face.end();++fi)
if(!(*fi).IsD() && (*fi).IsS()) ++selCnt;
return selCnt;
}
static size_t CountVertex(MeshType &m)
{
size_t selCnt=0;
VertexIterator vi;
for(vi=m.vert.begin();vi!=m.vert.end();++vi)
if(!(*vi).IsD() && (*vi).IsS()) ++selCnt;
return selCnt;
}
static size_t InvertFace(MeshType &m)
{
size_t selCnt=0;
FaceIterator fi;
for(fi=m.face.begin();fi!=m.face.end();++fi)
if(!(*fi).IsD())
{
if((*fi).IsS()) (*fi).ClearS();
else {
(*fi).SetS();
++selCnt;
}
}
return selCnt;
}
static size_t InvertVertex(MeshType &m)
{
size_t selCnt=0;
VertexIterator vi;
for(vi=m.vert.begin();vi!=m.vert.end();++vi)
if(!(*vi).IsD())
{
if((*vi).IsS()) (*vi).ClearS();
else {
(*vi).SetS();
++selCnt;
}
}
return selCnt;
}
/// \brief Select all the vertices that are touched by at least a single selected faces
static size_t VertexFromFaceLoose(MeshType &m)
{
size_t selCnt=0;
ClearVertex(m);
FaceIterator fi;
for(fi = m.face.begin(); fi != m.face.end(); ++fi)
if( !(*fi).IsD() && (*fi).IsS())
{
if( !(*fi).V(0)->IsS()) { (*fi).V(0)->SetS(); ++selCnt; }
if( !(*fi).V(1)->IsS()) { (*fi).V(1)->SetS(); ++selCnt; }
if( !(*fi).V(2)->IsS()) { (*fi).V(2)->SetS(); ++selCnt; }
}
return selCnt;
}
/// \brief Select ONLY the vertices that are touched ONLY by selected faces
/** In other words all the vertices having all the faces incident on them selected.
\warning Isolated vertices will not selected.
*/
static size_t VertexFromFaceStrict(MeshType &m)
{
VertexFromFaceLoose(m);
FaceIterator fi;
for(fi = m.face.begin(); fi != m.face.end(); ++fi)
if( !(*fi).IsD() && !(*fi).IsS())
{
(*fi).V(0)->ClearS();
(*fi).V(1)->ClearS();
(*fi).V(2)->ClearS();
}
return CountVertex(m);
}
/// \brief Select ONLY the faces with ALL the vertices selected
static size_t FaceFromVertexStrict(MeshType &m)
{
size_t selCnt=0;
ClearFace(m);
FaceIterator fi;
for(fi = m.face.begin(); fi != m.face.end(); ++fi)
if( !(*fi).IsD())
{
if((*fi).V(0)->IsS() && (*fi).V(1)->IsS() && (*fi).V(2)->IsS())
{
(*fi).SetS();
++selCnt;
}
}
return selCnt;
}
/// \brief Select all the faces with at least one selected vertex
static size_t FaceFromVertexLoose(MeshType &m)
{
size_t selCnt=0;
ClearFace(m);
FaceIterator fi;
for(fi = m.face.begin(); fi != m.face.end(); ++fi)
if( !(*fi).IsD() && !(*fi).IsS())
{
if((*fi).V(0)->IsS() || (*fi).V(1)->IsS() || (*fi).V(2)->IsS())
{
(*fi).SetS();
++selCnt;
}
}
return selCnt;
}
static size_t VertexFromBorderFlag(MeshType &m)
{
size_t selCnt=0;
ClearVertex(m);
VertexIterator vi;
for(vi = m.vert.begin(); vi != m.vert.end(); ++vi)
if( !(*vi).IsD() )
{
if((*vi).IsB() )
{
(*vi).SetS();
++selCnt;
}
}
return selCnt;
}
static size_t FaceFromBorderFlag(MeshType &m)
{
size_t selCnt=0;
ClearFace(m);
FaceIterator fi;
for(fi = m.face.begin(); fi != m.face.end(); ++fi)
if( !(*fi).IsD() )
{
if((*fi).IsB(0) || (*fi).IsB(1) || (*fi).IsB(2))
{
(*fi).SetS();
++selCnt;
}
}
return selCnt;
}
/// \brief This function select the faces that have an edge outside the given range.
static size_t FaceOutOfRangeEdge(MeshType &m, ScalarType MinEdgeThr=0, ScalarType MaxEdgeThr=(std::numeric_limits<ScalarType>::max)())
{
FaceIterator fi;
size_t count_fd = 0;
MinEdgeThr=MinEdgeThr*MinEdgeThr;
MaxEdgeThr=MaxEdgeThr*MaxEdgeThr;
for(fi=m.face.begin(); fi!=m.face.end();++fi)
if(!(*fi).IsD())
{
for(unsigned int i=0;i<3;++i)
{
const ScalarType squaredEdge=SquaredDistance((*fi).V0(i)->cP(),(*fi).V1(i)->cP());
if((squaredEdge<=MinEdgeThr) || (squaredEdge>=MaxEdgeThr) )
{
count_fd++;
(*fi).SetS();
break; // skip the rest of the edges of the tri
}
}
}
return count_fd;
}
/// \brief This function expand current selection to cover the whole connected component.
static size_t FaceConnectedFF(MeshType &m)
{
// it also assumes that the FF adjacency is well computed.
assert (HasFFAdjacency(m));
UpdateFlags<MeshType>::FaceClearV(m);
std::deque<FacePointer> visitStack;
size_t selCnt=0;
FaceIterator fi;
for(fi = m.face.begin(); fi != m.face.end(); ++fi)
if( !(*fi).IsD() && (*fi).IsS() && !(*fi).IsV() )
visitStack.push_back(&*fi);
while(!visitStack.empty())
{
FacePointer fp = visitStack.front();
visitStack.pop_front();
assert(!fp->IsV());
fp->SetV();
for(int i=0;i<3;++i) {
FacePointer ff = fp->FFp(i);
if(! ff->IsS())
{
ff->SetS();
++selCnt;
visitStack.push_back(ff);
assert(!ff->IsV());
}
}
}
return selCnt;
}
/// \brief Select ONLY the faces whose quality is in the specified closed interval.
static size_t FaceFromQualityRange(MeshType &m,float minq, float maxq)
{
size_t selCnt=0;
ClearFace(m);
FaceIterator fi;
assert(HasPerFaceQuality(m));
for(fi=m.face.begin();fi!=m.face.end();++fi)
if(!(*fi).IsD())
{
if( (*fi).Q()>=minq && (*fi).Q()<=maxq )
{
(*fi).SetS();
++selCnt;
}
}
return selCnt;
}
/// \brief Select ONLY the vertices whose quality is in the specified closed interval.
static size_t VertexFromQualityRange(MeshType &m,float minq, float maxq)
{
size_t selCnt=0;
ClearVertex(m);
VertexIterator vi;
assert(HasPerVertexQuality(m));
for(vi=m.vert.begin();vi!=m.vert.end();++vi)
if(!(*vi).IsD())
{
if( (*vi).Q()>=minq && (*vi).Q()<=maxq )
{
(*vi).SetS();
++selCnt;
}
}
return selCnt;
}
static int VertexInBox( MeshType & m, const Box3Type &bb)
{
int selCnt=0;
for (VertexIterator vi = m.vert.begin(); vi != m.vert.end(); ++vi) if(!(*vi).IsD())
{
if(bb.IsIn((*vi).cP()) ) {
(*vi).SetS();
++selCnt;
}
}
return selCnt;
}
void VertexNonManifoldEdges(MeshType &m)
{
assert(HasFFTopology(m));
VertexClear(m);
for (FaceIterator fi = m.face.begin(); fi != m.face.end(); ++fi) if (!fi->IsD())
{
for(int i=0;i<3;++i)
if(!IsManifold(*fi,i)){
(*fi).V0(i)->SetS();
(*fi).V1(i)->SetS();
}
}
}
}; // end class
} // End namespace
} // End namespace
#endif

View File

@ -1,105 +0,0 @@
/****************************************************************************
* 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: position.h,v $
****************************************************************************/
#ifndef __VCG_TRI_UPDATE_TEXTURE
#define __VCG_TRI_UPDATE_TEXTURE
//#include <vcg/space/plane.h>
namespace vcg {
namespace tri {
/// \ingroup trimesh
/// \headerfile texture.h vcg/complex/trimesh/update/texture.h
/// \brief This class is used to update vertex position according to a transformation matrix.
template <class ComputeMeshType>
class UpdateTexture
{
public:
typedef ComputeMeshType MeshType;
typedef typename MeshType::ScalarType ScalarType;
typedef typename MeshType::VertexType VertexType;
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;
static void WedgeTexFromPlanar(ComputeMeshType &m, Plane3<ScalarType> &pl)
{
FaceIterator fi;
for(fi=m.face.begin();fi!=m.face.end();++fi)
if(!(*fi).IsD())
{
}
}
static void WedgeTexFromCamera(ComputeMeshType &m, Plane3<ScalarType> &pl)
{
}
/// Currently texture coords are kept for ALL the triangles of a mesh. The texture id is stored with each face.
/// if a given face should not have tex coord it has the default -1 value for texture ID.
/// This function will add an new fake texture, add that to the list of textures and change all the -1 id to that value.
static void WedgeTexRemoveNull(ComputeMeshType &m, const std::string &texturename)
{
bool found=false;
FaceIterator fi;
// first loop lets check that there are -1 indexed textured face
for(fi=m.face.begin();fi!=m.face.end();++fi)
if(!(*fi).IsD()) if((*fi).WT(0).N()==-1) found = true;
if(!found) return;
m.textures.push_back(texturename);
int nullId=m.textures.size()-1;
for(fi=m.face.begin();fi!=m.face.end();++fi)
if(!(*fi).IsD()) if((*fi).WT(0).N()==-1)
{
(*fi).WT(0).N() = nullId;
(*fi).WT(1).N() = nullId;
(*fi).WT(2).N() = nullId;
}
}
}; // end class
} // End namespace
} // End namespace
#endif

View File

@ -1,489 +0,0 @@
/****************************************************************************
* 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.20 2008/04/04 10:27:34 cignoni
minor changes to the topology correctness checks
Revision 1.19 2007/05/29 00:07:06 ponchio
VFi++ -> ++VFi
Revision 1.18 2006/02/27 19:26:14 spinelli
minor bug in Face-Face topology loop fixed
Revision 1.17 2006/02/27 11:56:48 spinelli
minor bug in Face-Face topology loop fixed
Revision 1.16 2005/11/10 15:36:42 cignoni
Added clarifying comment in an assert
Revision 1.15 2004/10/20 07:33:10 cignoni
removed FaceBorderFlags (already present in update/flags.h)
Revision 1.14 2004/10/18 17:10:22 ganovelli
added ::FaceBorderFLags
Revision 1.13 2004/10/01 15:58:00 ponchio
Added include <vector>
Revision 1.12 2004/09/09 13:02:12 ponchio
Linux compatible path in #include
Revision 1.11 2004/08/07 16:18:20 pietroni
addet testFFTopology and testVFTopology functions used to test the rispective topology....
Revision 1.10 2004/07/15 11:35:08 ganovelli
Vfb to VFp
Revision 1.9 2004/07/15 00:13:39 cignoni
Better doxigen documentation
Revision 1.8 2004/06/02 16:42:44 ganovelli
typename for gcc compilation
Revision 1.7 2004/06/02 16:28:22 ganovelli
minor changes (swap =>> math::Swap)
Revision 1.6 2004/05/10 15:23:43 cignoni
Changed a FV -> VF in VertexFace topology computation
Revision 1.5 2004/05/06 15:24:38 pietroni
changed names to topology functions
Revision 1.4 2004/03/31 14:44:43 cignoni
Added Vertex-Face Topology
Revision 1.3 2004/03/12 15:22:19 cignoni
Written some documentation and added to the trimes doxygen module
Revision 1.2 2004/03/05 21:49:21 cignoni
First working version for face face
Revision 1.1 2004/03/04 00:53:24 cignoni
Initial commit
****************************************************************************/
#ifndef __VCG_TRI_UPDATE_TOPOLOGY
#define __VCG_TRI_UPDATE_TOPOLOGY
#include <algorithm>
#include <vector>
#include <vcg/simplex/face/pos.h>
#include <vcg/complex/trimesh/base.h>
namespace vcg {
namespace tri {
/// \ingroup trimesh
/// \headerfile topology.h vcg/complex/trimesh/update/topology.h
/// \brief Generation of per-vertex and per-face topological information.
template <class UpdateMeshType>
class UpdateTopology
{
public:
typedef UpdateMeshType MeshType;
typedef typename MeshType::VertexType VertexType;
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;
/// \headerfile topology.h vcg/complex/trimesh/update/topology.h
/// \brief Auxiliairy data structure for computing face face adjacency information.
/**
It identifies and edge storing two vertex pointer and a face pointer where it belong.
*/
class PEdge
{
public:
VertexPointer v[2]; // the two Vertex pointer are ordered!
FacePointer f; // the face where this edge belong
int z; // index in [0..2] of the edge of the face
PEdge() {}
void Set( FacePointer pf, const int nz )
{
assert(pf!=0);
assert(nz>=0);
assert(nz<pf->VN());
v[0] = pf->V(nz);
v[1] = pf->V(pf->Next(nz));
assert(v[0] != v[1]); // The face pointed by 'f' is Degenerate (two coincident vertexes)
if( v[0] > v[1] ) math::Swap(v[0],v[1]);
f = pf;
z = nz;
}
inline bool operator < ( const PEdge & pe ) const
{
if( v[0]<pe.v[0] ) return true;
else if( v[0]>pe.v[0] ) return false;
else return v[1] < pe.v[1];
}
inline bool operator == ( const PEdge & pe ) const
{
return v[0]==pe.v[0] && v[1]==pe.v[1];
}
};
// Fill a vector with all the edges of the mesh.
// each edge is stored in the vector the number of times that it appears in the mesh, with the referring face.
// optionally it can skip the faux edges (to retrieve only the real edges of a triangulated polygonal mesh)
static void FillEdgeVector(MeshType &m, std::vector<PEdge> &e, bool includeFauxEdge=true)
{
FaceIterator pf;
typename std::vector<PEdge>::iterator p;
// Alloco il vettore ausiliario
//e.resize(m.fn*3);
FaceIterator fi;
int n_edges = 0;
for(fi = m.face.begin(); fi != m.face.end(); ++fi) if(! (*fi).IsD()) n_edges+=(*fi).VN();
e.resize(n_edges);
p = e.begin();
for(pf=m.face.begin();pf!=m.face.end();++pf)
if( ! (*pf).IsD() )
for(int j=0;j<(*pf).VN();++j)
if(includeFauxEdge || !(*pf).IsF(j))
{
(*p).Set(&(*pf),j);
++p;
}
if(includeFauxEdge) assert(p==e.end());
else e.resize(p-e.begin());
}
static void FillUniqueEdgeVector(MeshType &m, std::vector<PEdge> &Edges, bool includeFauxEdge=true)
{
FillEdgeVector(m,Edges,includeFauxEdge);
sort(Edges.begin(), Edges.end()); // Lo ordino per vertici
typename std::vector< PEdge>::iterator newEnd = std::unique(Edges.begin(), Edges.end());
typename std::vector<PEdge>::iterator ei;
Edges.resize(newEnd-Edges.begin());
}
/// \brief Update the Face-Face topological relation by allowing to retrieve for each face what other faces shares their edges.
static void FaceFace(MeshType &m)
{
assert(HasFFAdjacency(m));
if( m.fn == 0 ) return;
std::vector<PEdge> e;
FillEdgeVector(m,e);
sort(e.begin(), e.end()); // Lo ordino per vertici
int ne = 0; // Numero di edge reali
typename std::vector<PEdge>::iterator pe,ps;
ps = e.begin();pe=e.begin();
//for(ps = e.begin(),pe=e.begin();pe<=e.end();++pe) // Scansione vettore ausiliario
do
{
if( pe==e.end() || !(*pe == *ps) ) // Trovo blocco di edge uguali
{
typename std::vector<PEdge>::iterator q,q_next;
for (q=ps;q<pe-1;++q) // Scansione facce associate
{
assert((*q).z>=0);
//assert((*q).z< 3);
q_next = q;
++q_next;
assert((*q_next).z>=0);
assert((*q_next).z< (*q_next).f->VN());
(*q).f->FFp(q->z) = (*q_next).f; // Collegamento in lista delle facce
(*q).f->FFi(q->z) = (*q_next).z;
}
assert((*q).z>=0);
assert((*q).z< (*q).f->VN());
(*q).f->FFp((*q).z) = ps->f;
(*q).f->FFi((*q).z) = ps->z;
ps = pe;
++ne; // Aggiorno il numero di edge
}
if(pe==e.end()) break;
++pe;
} while(true);
}
/// \brief Update the Vertex-Face topological relation.
/**
The function allows to retrieve for each vertex the list of faces sharing this vertex.
*/
static void VertexFace(MeshType &m)
{
if(!m.HasVFTopology()) return;
VertexIterator vi;
FaceIterator fi;
for(vi=m.vert.begin();vi!=m.vert.end();++vi)
{
(*vi).VFp() = 0;
(*vi).VFi() = 0;
}
for(fi=m.face.begin();fi!=m.face.end();++fi)
if( ! (*fi).IsD() )
{
for(int j=0;j<(*fi).VN();++j)
{
(*fi).VFp(j) = (*fi).V(j)->VFp();
(*fi).VFi(j) = (*fi).V(j)->VFi();
(*fi).V(j)->VFp() = &(*fi);
(*fi).V(j)->VFi() = j;
}
}
}
/// \headerfile topology.h vcg/complex/trimesh/update/topology.h
/// \brief Auxiliairy data structure for computing face face adjacency information.
/**
It identifies and edge storing two vertex pointer and a face pointer where it belong.
*/
class PEdgeTex
{
public:
typename FaceType::TexCoordType v[2]; // the two Vertex pointer are ordered!
FacePointer f; // the face where this edge belong
int z; // index in [0..2] of the edge of the face
PEdgeTex() {}
void Set( FacePointer pf, const int nz )
{
assert(pf!=0);
assert(nz>=0);
assert(nz<3);
v[0] = pf->WT(nz);
v[1] = pf->WT(pf->Next(nz));
assert(v[0] != v[1]); // The face pointed by 'f' is Degenerate (two coincident vertexes)
if( v[1] < v[0] ) std::swap(v[0],v[1]);
f = pf;
z = nz;
}
inline bool operator < ( const PEdgeTex & pe ) const
{
if( v[0]<pe.v[0] ) return true;
else if( pe.v[0]<v[0] ) return false;
else return v[1] < pe.v[1];
}
inline bool operator == ( const PEdgeTex & pe ) const
{
return (v[0]==pe.v[0]) && (v[1]==pe.v[1]);
}
inline bool operator != ( const PEdgeTex & pe ) const
{
return (v[0]!=pe.v[0]) || (v[1]!=pe.v[1]);
}
};
/// \brief Update the Face-Face topological relation
/**
The function allows to retrieve for each face what other faces shares their edges.
*/
static void FaceFaceFromTexCoord(MeshType &m)
{
// assert(HasFFTopology(m));
assert(HasPerWedgeTexCoord(m));
std::vector<PEdgeTex> e;
FaceIterator pf;
typename std::vector<PEdgeTex>::iterator p;
if( m.fn == 0 ) return;
// e.resize(m.fn*3); // Alloco il vettore ausiliario
FaceIterator fi;
int n_edges = 0;
for(fi = m.face.begin(); fi != m.face.end(); ++fi) if(! (*fi).IsD()) n_edges+=(*fi).VN();
e.resize(n_edges);
p = e.begin();
for(pf=m.face.begin();pf!=m.face.end();++pf) // Lo riempio con i dati delle facce
if( ! (*pf).IsD() )
for(int j=0;j<(*pf).VN();++j)
{
if( (*pf).WT(j) != (*pf).WT((*pf).Next(j)))
{
(*p).Set(&(*pf),j);
++p;
}
}
e.resize(p-e.begin()); // remove from the end of the edge vector the unitiailized ones
assert(p==e.end());
sort(e.begin(), e.end());
int ne = 0; // number of real edges
typename std::vector<PEdgeTex>::iterator pe,ps;
ps = e.begin();pe=e.begin();
//for(ps = e.begin(),pe=e.begin();pe<=e.end();++pe) // Scansione vettore ausiliario
do
{
if( pe==e.end() || (*pe) != (*ps) ) // Trovo blocco di edge uguali
{
typename std::vector<PEdgeTex>::iterator q,q_next;
for (q=ps;q<pe-1;++q) // Scansione facce associate
{
assert((*q).z>=0);
assert((*q).z< 3);
q_next = q;
++q_next;
assert((*q_next).z>=0);
assert((*q_next).z< (*q_next).f->VN());
(*q).f->FFp(q->z) = (*q_next).f; // Collegamento in lista delle facce
(*q).f->FFi(q->z) = (*q_next).z;
}
assert((*q).z>=0);
assert((*q).z< (*q).f->VN());
(*q).f->FFp((*q).z) = ps->f;
(*q).f->FFi((*q).z) = ps->z;
ps = pe;
++ne; // Aggiorno il numero di edge
}
if(pe==e.end()) break;
++pe;
} while(true);
}
/// \brief Test correctness of VFtopology
static void TestVertexFace(MeshType &m)
{
SimpleTempData<typename MeshType::VertContainer, int > numVertex(m.vert,0);
if(!m.HasVFTopology()) return;
FaceIterator fi;
for(fi=m.face.begin();fi!=m.face.end();++fi)
{
if (!(*fi).IsD())
{
numVertex[(*fi).V0(0)]++;
numVertex[(*fi).V1(0)]++;
numVertex[(*fi).V2(0)]++;
}
}
VertexIterator vi;
vcg::face::VFIterator<FaceType> VFi;
for(vi=m.vert.begin();vi!=m.vert.end();++vi)
{
if (!vi->IsD())
if(vi->VFp()!=0) // unreferenced vertices MUST have VF == 0;
{
int num=0;
assert(vi->VFp() >= &*m.face.begin());
assert(vi->VFp() <= &m.face.back());
VFi.f=vi->VFp();
VFi.z=vi->VFi();
while (!VFi.End())
{
num++;
assert(!VFi.F()->IsD());
assert((VFi.F()->V(VFi.I()))==&(*vi));
++VFi;
}
int num1=numVertex[&(*vi)];
assert(num==num1);
/*assert(num>1);*/
}
}
}
/// \brief Test correctness of FFtopology (only for 2Manifold Meshes!)
static void TestFaceFace(MeshType &m)
{
if(!m.HasFFTopology()) return;
for(FaceIterator fi=m.face.begin();fi!=m.face.end();++fi)
{
if (!fi->IsD())
{
for (int i=0;i<(*fi).VN();i++)
{
FaceType *ffpi=fi->FFp(i);
int e=fi->FFi(i);
//invariant property of FF topology for two manifold meshes
assert(ffpi->FFp(e) == &(*fi));
assert(ffpi->FFi(e) == i);
// Test that the two faces shares the same edge
// Vertices of the i-th edges of the first face
VertexPointer v0i= fi->V0(i);
VertexPointer v1i= fi->V1(i);
// Vertices of the corresponding edge on the other face
VertexPointer ffv0i= ffpi->V0(e);
VertexPointer ffv1i= ffpi->V1(e);
assert( (ffv0i==v0i) || (ffv0i==v1i) );
assert( (ffv1i==v0i) || (ffv1i==v1i) );
}
}
}
}
}; // end class
} // End namespace
} // End namespace
#endif