Riccardo Elitropi
b3abcf73c3
- aggiunta funzione per il calcolo della curvatura in un vertice per una superficie TriMesh.
1169 lines
41 KiB
C++
1169 lines
41 KiB
C++
//----------------------------------------------------------------------------
|
||
// EgalTech 2015-2023
|
||
//----------------------------------------------------------------------------
|
||
// File : SurfTriMeshFaceting.cpp Data : 09.12.23 Versione : 2.5l2
|
||
// Contenuto : Implementazione della classe Superfici TriMesh.
|
||
//
|
||
//
|
||
//
|
||
// Modifiche : 26.03.14 DS Creazione modulo.
|
||
//
|
||
//
|
||
//----------------------------------------------------------------------------
|
||
|
||
//--------------------------- Include ----------------------------------------
|
||
#include "stdafx.h"
|
||
#include "SurfTriMesh.h"
|
||
#include "GeoConst.h"
|
||
#include "PolygonPlane.h"
|
||
#include "CurveLine.h"
|
||
#include "/EgtDev/Include/EgtPointerOwner.h"
|
||
#include <array>
|
||
#include <set>
|
||
#include <unordered_map>
|
||
#define EIGEN_NO_IO
|
||
#include "/EgtDev/Extern/Eigen/Dense"
|
||
|
||
using namespace std ;
|
||
|
||
//----------------------------------------------------------------------------
|
||
bool
|
||
SurfTriMesh::VerifyFaceting( void) const
|
||
{
|
||
if ( m_bFaceted)
|
||
return true ;
|
||
return (const_cast<SurfTriMesh*>(this))->UpdateFaceting() ;
|
||
}
|
||
|
||
//----------------------------------------------------------------------------
|
||
bool
|
||
SurfTriMesh::UpdateFaceting( void)
|
||
{
|
||
// salvo vecchio vettore facce
|
||
INTVECTOR vOldFacet = m_vFacet ;
|
||
|
||
// reset faceting
|
||
m_bFaceted = false ;
|
||
m_vFacet.clear() ;
|
||
for ( int i = 0 ; i < int( m_vTria.size()) ; ++ i)
|
||
m_vTria[i].nIdFacet = SVT_NULL ;
|
||
m_bFacEdged = false ;
|
||
|
||
// indice faccia corrente
|
||
int nFacet = -1 ;
|
||
|
||
// ricostruisco le sfaccettature come definite in precedenza (dove possibile)
|
||
bool bOk = true ;
|
||
for ( int j = 0 ; j < int( vOldFacet.size()) ; ++ j) {
|
||
int i = vOldFacet[j] ;
|
||
// salto triangoli inesistenti o già assegnati
|
||
if ( i >= int( m_vTria.size()) ||
|
||
m_vTria[i].nIdVert[0] == SVT_DEL ||
|
||
m_vTria[i].nIdFacet != SVT_NULL)
|
||
continue ;
|
||
// assegno indice di faccia al triangolo
|
||
m_vTria[i].nIdFacet = ++ nFacet ;
|
||
m_vFacet.push_back( i) ;
|
||
// aggiorno faccia
|
||
if ( ! UpdateOneFace( nFacet, i))
|
||
bOk = false ;
|
||
}
|
||
|
||
// ricostruisco le altre sfaccettature
|
||
for ( int i = 0 ; i < int( m_vTria.size()) ; ++ i) {
|
||
// salto triangoli cancellati o già assegnati
|
||
if ( m_vTria[i].nIdVert[0] == SVT_DEL ||
|
||
m_vTria[i].nIdFacet != SVT_NULL)
|
||
continue ;
|
||
// assegno indice di faccia al triangolo
|
||
m_vTria[i].nIdFacet = ++ nFacet ;
|
||
m_vFacet.push_back( i) ;
|
||
// aggiorno faccia
|
||
if ( ! UpdateOneFace( nFacet, i))
|
||
bOk = false ;
|
||
}
|
||
|
||
// se ci sono stati problemi, salvo nel log
|
||
if ( ! bOk)
|
||
LOG_ERROR( GetEGkLogger(), "SurfTM : UpdateFaceting error")
|
||
|
||
// calcolo facce piane effettuato
|
||
m_bFaceted = bOk ;
|
||
return bOk ;
|
||
}
|
||
|
||
//----------------------------------------------------------------------------
|
||
bool
|
||
SurfTriMesh::UpdateOneFace( int nFacet, int nT)
|
||
{
|
||
// piano del triangolo
|
||
Plane3d plPlane ;
|
||
if ( ! plPlane.Set( m_vVert[m_vTria[nT].nIdVert[0]].ptP, m_vTria[nT].vtN)) {
|
||
LOG_ERROR( GetEGkLogger(), "SurfTM : UpdateFaceting error in triangle data")
|
||
return false ;
|
||
}
|
||
// set di triangoli da aggiornare
|
||
set<int> stTria ;
|
||
stTria.insert( nT) ;
|
||
// finchè set non vuoto
|
||
bool bOk = true ;
|
||
while ( ! stTria.empty()) {
|
||
// tolgo un triangolo dal set
|
||
const auto iIt = stTria.begin() ;
|
||
int nT = *iIt ;
|
||
stTria.erase( iIt) ;
|
||
// aggiorno i triangoli adiacenti
|
||
for ( int j = 0 ; j < 3 ; ++ j) {
|
||
int nAdjT = m_vTria[nT].nIdAdjac[j] ;
|
||
if ( nAdjT != SVT_NULL && m_vTria[nAdjT].nIdFacet == SVT_NULL) {
|
||
if ( ! UpdateTriaFaceting( nT, nFacet, plPlane, nAdjT))
|
||
bOk = false ;
|
||
if ( m_vTria[nAdjT].nIdFacet == nFacet)
|
||
stTria.insert( nAdjT) ;
|
||
}
|
||
}
|
||
}
|
||
return bOk ;
|
||
}
|
||
|
||
//----------------------------------------------------------------------------
|
||
bool
|
||
SurfTriMesh::UpdateTriaFaceting( int nRefT, int nFacet, const Plane3d& plPlane, int nT)
|
||
{
|
||
// verifica scostamento della normale da quella del piano
|
||
if ( plPlane.GetVersN() * m_vTria[nT].vtN < m_dCosBndAng)
|
||
return true ;
|
||
// verifica scostamento del vertice opposto al lato in comune dal piano
|
||
int nV = SVT_NULL ;
|
||
for ( int i = 0 ; i < 3 ; ++i) {
|
||
if ( m_vTria[nT].nIdAdjac[i] == nRefT) {
|
||
nV = Prev( i) ;
|
||
break ;
|
||
}
|
||
}
|
||
if ( nV == SVT_NULL)
|
||
return true ;
|
||
double dTol = max( min( m_dLinTol, 20 * EPS_SMALL), 2 * EPS_SMALL) ;
|
||
if ( ! PointInPlaneEpsilon( m_vVert[m_vTria[nT].nIdVert[nV]].ptP, plPlane, dTol))
|
||
return true ;
|
||
// il triangolo fa parte della faccia
|
||
m_vTria[nT].nIdFacet = nFacet ;
|
||
return true ;
|
||
}
|
||
|
||
//----------------------------------------------------------------------------
|
||
int
|
||
SurfTriMesh::GetFacetCount( void) const
|
||
{
|
||
// la superficie deve essere validata
|
||
if ( m_nStatus != OK)
|
||
return 0 ;
|
||
// verifico stato sfaccettatura
|
||
if ( ! VerifyFaceting())
|
||
return 0 ;
|
||
// restituisco il numero
|
||
return int( m_vFacet.size()) ;
|
||
}
|
||
|
||
//----------------------------------------------------------------------------
|
||
bool
|
||
SurfTriMesh::SetFacet( int nInd, int nT)
|
||
{
|
||
// recupero la dimensione originale
|
||
int nPrevSize = int( m_vFacet.size()) ;
|
||
// determino la dimensione necessaria
|
||
int nNewSize = max( nInd + 1, nPrevSize) ;
|
||
// se necessaria dimensione maggiore
|
||
if ( nNewSize > nPrevSize) {
|
||
// espando vettore
|
||
try { m_vFacet.resize( nNewSize) ; }
|
||
catch (...) { return false ; }
|
||
// inizializzo a cancellate le eventuali facce intermedie
|
||
if ( ( nNewSize - nPrevSize) > 1) {
|
||
for ( int i = nPrevSize ; i < nNewSize ; ++ i)
|
||
m_vFacet[i] = SVT_DEL ;
|
||
}
|
||
}
|
||
// inserisco la faccia
|
||
m_vFacet[nInd] = nT ;
|
||
return true ;
|
||
}
|
||
|
||
//----------------------------------------------------------------------------
|
||
int
|
||
SurfTriMesh::GetFacetFromTria( int nT) const
|
||
{
|
||
// la superficie deve essere validata
|
||
if ( m_nStatus != OK)
|
||
return SVT_NULL ;
|
||
// verifico stato sfaccettatura
|
||
if ( ! VerifyFaceting())
|
||
return SVT_NULL ;
|
||
// l'indice del triangolo deve essere nei limiti
|
||
if ( nT < 0 || nT >= int( m_vTria.size()) || m_vTria[nT].nIdVert[0] == SVT_DEL)
|
||
return SVT_NULL ;
|
||
// restituisco l'indice di faccia
|
||
return m_vTria[nT].nIdFacet ;
|
||
}
|
||
|
||
//----------------------------------------------------------------------------
|
||
bool
|
||
SurfTriMesh::GetAllTriaInFacet( int nF, INTVECTOR& vT) const
|
||
{
|
||
// la superficie deve essere validata
|
||
if ( m_nStatus != OK)
|
||
return false ;
|
||
// verifico stato sfaccettatura
|
||
if ( ! VerifyFaceting())
|
||
return false ;
|
||
// l'indice della faccia deve essere nei limiti
|
||
if ( nF < 0 || nF >= int( m_vFacet.size()))
|
||
return false ;
|
||
// recupero l'indice del primo triangolo della faccia
|
||
int nT = m_vFacet[nF] ;
|
||
// incremento time stamp
|
||
++ m_nTimeStamp ;
|
||
// recupero i triangoli della faccia
|
||
vT.clear() ;
|
||
vT.reserve( 10) ;
|
||
vT.push_back( nT) ;
|
||
m_vTria[nT].nTemp = m_nTimeStamp ;
|
||
return VerifyAdjacTriaFacet( vT) ;
|
||
}
|
||
|
||
//----------------------------------------------------------------------------
|
||
bool
|
||
SurfTriMesh::VerifyAdjacTriaFacet( INTVECTOR& vT) const
|
||
{
|
||
INTVECTOR vStack{ vT.back()} ;
|
||
while ( ! vStack.empty()) {
|
||
array<int,3> vNew{ SVT_NULL, SVT_NULL, SVT_NULL} ;
|
||
for ( int j = 0 ; j < 3 ; ++ j) {
|
||
int nAdjT = m_vTria[vStack.back()].nIdAdjac[j] ;
|
||
if ( nAdjT != SVT_NULL &&
|
||
m_vTria[nAdjT].nTemp != m_nTimeStamp &&
|
||
m_vTria[nAdjT].nIdFacet == m_vTria[vStack.back()].nIdFacet) {
|
||
vT.push_back( nAdjT) ;
|
||
vNew[j] = nAdjT ;
|
||
m_vTria[nAdjT].nTemp = m_nTimeStamp ;
|
||
}
|
||
}
|
||
vStack.pop_back() ;
|
||
for ( int j = 0 ; j < 3 ; ++ j) {
|
||
if ( vNew[j] != SVT_NULL)
|
||
vStack.push_back( vNew[j]) ;
|
||
}
|
||
}
|
||
|
||
return true ;
|
||
}
|
||
|
||
//----------------------------------------------------------------------------
|
||
bool
|
||
SurfTriMesh::GetAllVertInFacet( int nF, INTVECTOR& vVert) const
|
||
{
|
||
// recupero tutti i triangoli della faccia
|
||
INTVECTOR vTria ;
|
||
if ( ! GetAllTriaInFacet( nF, vTria))
|
||
return false ;
|
||
// ne ricavo i vertici
|
||
vVert.clear() ;
|
||
INTUNORDSET usTempVert ;
|
||
for ( int nT : vTria) {
|
||
for ( int i = 0 ; i < 3 ; ++ i) {
|
||
int nV = m_vTria[nT].nIdVert[i] ;
|
||
if ( usTempVert.find( nV) == usTempVert.end()) {
|
||
usTempVert.insert( nV) ;
|
||
vVert.push_back( nV) ;
|
||
}
|
||
}
|
||
}
|
||
|
||
return true ;
|
||
}
|
||
|
||
//----------------------------------------------------------------------------
|
||
bool
|
||
SurfTriMesh::GetFacetNearestEndPoint( int nF, const Point3d& ptNear, Point3d& ptEnd, Vector3d& vtN) const
|
||
{
|
||
// recupero l'elenco dei triangoli della faccia
|
||
INTVECTOR vTria ;
|
||
if ( ! GetAllTriaInFacet( nF, vTria))
|
||
return false ;
|
||
// ciclo sui triangoli e sui loro lati di bordo
|
||
bool bFound = false ;
|
||
double dMinSqDist = SQ_INFINITO ;
|
||
for ( int i = 0 ; i < int( vTria.size()) ; ++i) {
|
||
int nT = vTria[i] ;
|
||
for ( int j = 0 ; j < 3 ; ++ j) {
|
||
int nAdjT = m_vTria[nT].nIdAdjac[j] ;
|
||
if ( nAdjT == SVT_NULL ||
|
||
m_vTria[nAdjT].nIdFacet != nF) {
|
||
Point3d ptTest = m_vVert[m_vTria[nT].nIdVert[j]].ptP ;
|
||
double dSqDist = SqDist( ptTest, ptNear) ;
|
||
if ( dSqDist < dMinSqDist) {
|
||
// salvo i dati del punto
|
||
dMinSqDist = dSqDist ;
|
||
ptEnd = ptTest ;
|
||
bFound = true ;
|
||
// recupero la normale della faccia (non quella mediata del vertice)
|
||
vtN = m_vTria[nT].vtN ;
|
||
}
|
||
}
|
||
}
|
||
}
|
||
return bFound ;
|
||
}
|
||
|
||
//----------------------------------------------------------------------------
|
||
bool
|
||
SurfTriMesh::GetFacetNearestMidPoint( int nF, const Point3d& ptNear, Point3d& ptMid, Vector3d& vtN) const
|
||
{
|
||
// recupero l'elenco dei triangoli della faccia
|
||
INTVECTOR vTria ;
|
||
if ( ! GetAllTriaInFacet( nF, vTria))
|
||
return false ;
|
||
// ciclo sui triangoli e sui loro lati di bordo
|
||
bool bFound = false ;
|
||
double dMinSqDist = SQ_INFINITO ;
|
||
for ( int i = 0 ; i < int( vTria.size()) ; ++i) {
|
||
int nT = vTria[i] ;
|
||
for ( int j = 0 ; j < 3 ; ++ j) {
|
||
int k = Next( j) ;
|
||
int nAdjT = m_vTria[nT].nIdAdjac[j] ;
|
||
if ( nAdjT == SVT_NULL ||
|
||
m_vTria[nAdjT].nIdFacet != nF) {
|
||
Point3d ptTest = Media( m_vVert[m_vTria[nT].nIdVert[j]].ptP,
|
||
m_vVert[m_vTria[nT].nIdVert[k]].ptP, 0.5) ;
|
||
double dSqDist = SqDist( ptTest, ptNear) ;
|
||
if ( dSqDist < dMinSqDist) {
|
||
// salvo i dati del punto
|
||
dMinSqDist = dSqDist ;
|
||
ptMid = ptTest ;
|
||
bFound = true ;
|
||
// recupero la normale della faccia (non quella mediata del vertice)
|
||
vtN = m_vTria[nT].vtN ;
|
||
}
|
||
}
|
||
}
|
||
}
|
||
return bFound ;
|
||
}
|
||
|
||
//----------------------------------------------------------------------------
|
||
bool
|
||
SurfTriMesh::GetFacetLoops( int nF, POLYLINEVECTOR& vPL) const
|
||
{
|
||
// recupero l'elenco dei triangoli della faccia
|
||
INTVECTOR vTria ;
|
||
if ( ! GetAllTriaInFacet( nF, vTria))
|
||
return false ;
|
||
// incremento time stamp
|
||
++ m_nTimeStamp ;
|
||
// ciclo sui triangoli
|
||
for ( int i = 0 ; i < int( vTria.size()) ; ++ i) {
|
||
int nT = vTria[i] ;
|
||
// se triangolo non ancora visitato
|
||
if ( m_vTria[nT].nTemp != m_nTimeStamp) {
|
||
// determino i triangoli adiacenti
|
||
int nAdjT[3] ;
|
||
for ( int j = 0 ; j < 3 ; ++ j)
|
||
nAdjT[j] = m_vTria[nT].nIdAdjac[j] ;
|
||
// determino i facet adiacenti
|
||
int nAdjF[3] ;
|
||
for ( int j = 0 ; j < 3 ; ++ j)
|
||
nAdjF[j] = ( nAdjT[j] != SVT_NULL ? m_vTria[nAdjT[j]].nIdFacet : SVT_NULL) ;
|
||
// se tutti e tre i lati sono di contorno
|
||
if ( nAdjF[0] != nF && nAdjF[1] != nF && nAdjF[2] != nF) {
|
||
// ho trovato un loop
|
||
vPL.emplace_back() ;
|
||
vPL.back().AddUPoint( nAdjF[0], m_vVert[m_vTria[nT].nIdVert[0]].ptP) ;
|
||
vPL.back().AddUPoint( nAdjF[1], m_vVert[m_vTria[nT].nIdVert[1]].ptP) ;
|
||
vPL.back().AddUPoint( nAdjF[2], m_vVert[m_vTria[nT].nIdVert[2]].ptP) ;
|
||
vPL.back().AddUPoint( SVT_NULL, m_vVert[m_vTria[nT].nIdVert[0]].ptP) ;
|
||
// marco il triangolo come verificato
|
||
m_vTria[nT].nTemp = m_nTimeStamp ;
|
||
}
|
||
// se i due lati 0 e 1 sono di contorno
|
||
else if ( nAdjF[0] != nF && nAdjF[1] != nF) {
|
||
// ho trovato l'inizio di un loop
|
||
vPL.emplace_back() ;
|
||
vPL.back().AddUPoint( nAdjF[0], m_vVert[m_vTria[nT].nIdVert[0]].ptP) ;
|
||
vPL.back().AddUPoint( nAdjF[1], m_vVert[m_vTria[nT].nIdVert[1]].ptP) ;
|
||
vPL.back().AddUPoint( SVT_NULL, m_vVert[m_vTria[nT].nIdVert[2]].ptP) ;
|
||
// marco il triangolo come verificato
|
||
m_vTria[nT].nTemp = m_nTimeStamp ;
|
||
// cammino lungo il loop fino a chiuderlo
|
||
if ( ! MarchAlongFacetLoop( nF, nT, 2, m_nTimeStamp, vPL.back()))
|
||
return false ;
|
||
}
|
||
// se i due lati 1 e 2 sono di contorno
|
||
else if ( nAdjF[1] != nF && nAdjF[2] != nF) {
|
||
// ho trovato l'inizio di un loop
|
||
vPL.emplace_back() ;
|
||
vPL.back().AddUPoint( nAdjF[1], m_vVert[m_vTria[nT].nIdVert[1]].ptP) ;
|
||
vPL.back().AddUPoint( nAdjF[2], m_vVert[m_vTria[nT].nIdVert[2]].ptP) ;
|
||
vPL.back().AddUPoint( SVT_NULL, m_vVert[m_vTria[nT].nIdVert[0]].ptP) ;
|
||
// marco il triangolo come verificato
|
||
m_vTria[nT].nTemp = m_nTimeStamp ;
|
||
// cammino lungo il loop fino a chiuderlo
|
||
if ( ! MarchAlongFacetLoop( nF, nT, 0, m_nTimeStamp, vPL.back()))
|
||
return false ;
|
||
}
|
||
// se i due lati 2 e 0 sono di contorno
|
||
else if ( nAdjF[2] != nF && nAdjF[0] != nF) {
|
||
// ho trovato l'inizio di un loop
|
||
vPL.emplace_back() ;
|
||
vPL.back().AddUPoint( nAdjF[2], m_vVert[m_vTria[nT].nIdVert[2]].ptP) ;
|
||
vPL.back().AddUPoint( nAdjF[0], m_vVert[m_vTria[nT].nIdVert[0]].ptP) ;
|
||
vPL.back().AddUPoint( SVT_NULL, m_vVert[m_vTria[nT].nIdVert[1]].ptP) ;
|
||
// marco il triangolo come verificato
|
||
m_vTria[nT].nTemp = m_nTimeStamp ;
|
||
// cammino lungo il loop fino a chiuderlo
|
||
if ( ! MarchAlongFacetLoop( nF, nT, 1, m_nTimeStamp, vPL.back()))
|
||
return false ;
|
||
}
|
||
// se il lato 0 è di contorno
|
||
else if ( nAdjF[0] != nF) {
|
||
// ho trovato l'inizio di un loop
|
||
vPL.emplace_back() ;
|
||
vPL.back().AddUPoint( nAdjF[0], m_vVert[m_vTria[nT].nIdVert[0]].ptP) ;
|
||
vPL.back().AddUPoint( SVT_NULL, m_vVert[m_vTria[nT].nIdVert[1]].ptP) ;
|
||
// marco il triangolo come verificato
|
||
m_vTria[nT].nTemp = m_nTimeStamp ;
|
||
// cammino lungo il loop fino a chiuderlo
|
||
if ( ! MarchAlongFacetLoop( nF, nT, 1, m_nTimeStamp, vPL.back()))
|
||
return false ;
|
||
}
|
||
// se il lato 1 è di contorno
|
||
else if ( nAdjF[1] != nF) {
|
||
// ho trovato l'inizio di un loop
|
||
vPL.emplace_back() ;
|
||
vPL.back().AddUPoint( nAdjF[1], m_vVert[m_vTria[nT].nIdVert[1]].ptP) ;
|
||
vPL.back().AddUPoint( SVT_NULL, m_vVert[m_vTria[nT].nIdVert[2]].ptP) ;
|
||
// marco il triangolo come verificato
|
||
m_vTria[nT].nTemp = m_nTimeStamp ;
|
||
// cammino lungo il loop fino a chiuderlo
|
||
if ( ! MarchAlongFacetLoop( nF, nT, 2, m_nTimeStamp, vPL.back()))
|
||
return false ;
|
||
}
|
||
// se il lato 2 è di contorno
|
||
else if ( nAdjF[2] != nF) {
|
||
// ho trovato l'inizio di un loop
|
||
vPL.emplace_back() ;
|
||
vPL.back().AddUPoint( nAdjF[2], m_vVert[m_vTria[nT].nIdVert[2]].ptP) ;
|
||
vPL.back().AddUPoint( SVT_NULL, m_vVert[m_vTria[nT].nIdVert[0]].ptP) ;
|
||
// marco il triangolo come verificato
|
||
m_vTria[nT].nTemp = m_nTimeStamp ;
|
||
// cammino lungo il loop fino a chiuderlo
|
||
if ( ! MarchAlongFacetLoop( nF, nT, 0, m_nTimeStamp, vPL.back()))
|
||
return false ;
|
||
}
|
||
// altrimenti non c'è contorno
|
||
else {
|
||
// marco il triangolo come verificato
|
||
m_vTria[nT].nTemp = m_nTimeStamp ;
|
||
}
|
||
}
|
||
}
|
||
|
||
// normale di riferimento della faccia
|
||
Vector3d vtN = m_vTria[vTria[0]].vtN ;
|
||
// in prima posizione ci deve essere il loop esterno
|
||
bool bOutFirst = false ;
|
||
for ( int i = 0 ; i < int( vPL.size()) ; ++ i) {
|
||
Plane3d plPlane ;
|
||
double dArea ;
|
||
if ( ! vPL[i].IsClosedAndFlat( plPlane, dArea, 100 * EPS_SMALL))
|
||
return false ;
|
||
// se loop esterno
|
||
if ( vtN * plPlane.GetVersN() > 0) {
|
||
// se non c'è ancora loop esterno in prima posizione
|
||
if ( ! bOutFirst) {
|
||
// lo sposto in prima posizione
|
||
if ( i != 0)
|
||
swap( vPL[0], vPL[i]) ;
|
||
bOutFirst = true ;
|
||
}
|
||
// altrimenti errore
|
||
else
|
||
return false ;
|
||
}
|
||
}
|
||
|
||
return bOutFirst ;
|
||
}
|
||
|
||
//----------------------------------------------------------------------------
|
||
bool
|
||
SurfTriMesh::MarchAlongFacetLoop( int nF, int nT, int nV, int nTimeStamp, PolyLine& PL) const
|
||
{
|
||
// mi muovo lungo il loop, un triangolo alla volta
|
||
bool bEnd = false ;
|
||
while ( ! bEnd) {
|
||
if ( ! MarchOneFacetTria( nF, nT, nV, nTimeStamp, PL, bEnd))
|
||
return false ;
|
||
}
|
||
return true ;
|
||
}
|
||
|
||
//----------------------------------------------------------------------------
|
||
bool
|
||
SurfTriMesh::MarchOneFacetTria( int nF, int& nT, int& nV, int nTimeStamp,
|
||
PolyLine& PL, bool& bEnd) const
|
||
{
|
||
// verifico esistenza triangolo adiacente, sul lato dopo il vertice
|
||
if ( m_vTria[nT].nIdAdjac[nV] == SVT_NULL)
|
||
return false ;
|
||
int nAdjT = m_vTria[nT].nIdAdjac[nV] ;
|
||
// verifico appartenga alla stessa faccia
|
||
if ( m_vTria[nAdjT].nIdFacet != nF)
|
||
return false ;
|
||
// recupero il suo lato di adiacenza (e verifico non abbia più adiacenze con il triangolo di partenza)
|
||
int nAdjS = SVT_NULL ;
|
||
for ( int i = 0 ; i < 3 ; ++ i) {
|
||
if ( m_vTria[nAdjT].nIdAdjac[i] == nT) {
|
||
if ( nAdjS == SVT_NULL)
|
||
nAdjS = i ;
|
||
else
|
||
return false ;
|
||
}
|
||
}
|
||
if ( nAdjS == SVT_NULL)
|
||
return false ;
|
||
// vertice di fine adiacenza e indice del successivo lato
|
||
int nAdjV = Next( nAdjS) ;
|
||
// verifico se il lato successivo è un bordo
|
||
int nNextT = m_vTria[nAdjT].nIdAdjac[nAdjV] ;
|
||
int nNextF = ( nNextT != SVT_NULL ? m_vTria[nNextT].nIdFacet : SVT_NULL) ;
|
||
if ( nNextF != nF) {
|
||
// se già recuperato
|
||
if ( m_vTria[nAdjT].nTemp == nTimeStamp) {
|
||
bEnd = true ;
|
||
return true ;
|
||
}
|
||
// dichiaro triangolo analizzato
|
||
m_vTria[nAdjT].nTemp = nTimeStamp ;
|
||
// modifico l'ultima adiacenza e aggiungo il lato al loop
|
||
PL.ModifyLastParam( nNextF) ;
|
||
nAdjV = Next( nAdjV) ;
|
||
PL.AddUPoint( SVT_NULL, m_vVert[m_vTria[nAdjT].nIdVert[nAdjV]].ptP) ;
|
||
// verifico anche il successivo
|
||
nNextT = m_vTria[nAdjT].nIdAdjac[nAdjV] ;
|
||
nNextF = ( nNextT != SVT_NULL ? m_vTria[nNextT].nIdFacet : SVT_NULL) ;
|
||
if ( nNextF != nF) {
|
||
// modifico l'ultima adiacenza e aggiungo il lato al loop
|
||
PL.ModifyLastParam( nNextF) ;
|
||
nAdjV = Next( nAdjV) ;
|
||
PL.AddUPoint( SVT_NULL, m_vVert[m_vTria[nAdjT].nIdVert[nAdjV]].ptP) ;
|
||
}
|
||
}
|
||
// devo passare al triangolo adiacente
|
||
nT = nAdjT ;
|
||
nV = nAdjV ;
|
||
return true ;
|
||
}
|
||
|
||
//----------------------------------------------------------------------------
|
||
bool
|
||
SurfTriMesh::GetFacetAdjacencies( int nF, INTMATRIX& vAdj) const
|
||
{
|
||
// recupero i loop di contorno
|
||
POLYLINEVECTOR vPL ;
|
||
if ( ! GetFacetLoops( nF, vPL))
|
||
return false ;
|
||
// ne ricavo le adiacenze
|
||
vAdj.resize( vPL.size()) ;
|
||
for ( int i = 0 ; i < int( vPL.size()) ; ++ i) {
|
||
vAdj[i].reserve( vPL[i].GetPointNbr()) ;
|
||
double dU ;
|
||
for ( bool bFound = vPL[i].GetFirstU( dU, true) ; bFound ; bFound = vPL[i].GetNextU( dU, true))
|
||
vAdj[i].emplace_back( int( round( dU))) ;
|
||
}
|
||
|
||
return true ;
|
||
}
|
||
|
||
//----------------------------------------------------------------------------
|
||
bool
|
||
SurfTriMesh::GetFacetsContact( int nF1, int nF2, bool& bAdjac, Point3d& ptP1, Point3d& ptP2, double& dAng) const
|
||
{
|
||
// recupero i loop di contorno della prima faccia
|
||
POLYLINEVECTOR vPL ;
|
||
if ( ! GetFacetLoops( nF1, vPL))
|
||
return false ;
|
||
// verifico esistenza seconda faccia
|
||
if ( nF2 < 0 || nF2 >= int( m_vFacet.size()))
|
||
return false ;
|
||
// verifico se c'è un contatto con la seconda faccia e recupero i punti estremi della eventuale linea di contatto
|
||
bAdjac = false ;
|
||
for ( int i = 0 ; i < int( vPL.size()) ; ++ i) {
|
||
double dUs, dUe ;
|
||
Point3d ptPs, ptPe ;
|
||
for ( bool bFound = vPL[i].GetFirstULine( &dUs, &ptPs, &dUe, &ptPe) ; bFound ; bFound = vPL[i].GetNextULine( &dUs, &ptPs, &dUe, &ptPe)) {
|
||
if ( int( round( dUs)) == nF2) {
|
||
if ( ! bAdjac) {
|
||
bAdjac = true ;
|
||
ptP1 = ptPs ;
|
||
ptP2 = ptPe ;
|
||
}
|
||
else {
|
||
// parametri del segmento già definito
|
||
Vector3d vtT ;
|
||
double dLen ;
|
||
DirDist( ptP1, ptP2, vtT, dLen) ;
|
||
// punto start rispetto al segmento
|
||
double dPs = ( ptPs - ptP1) * vtT ;
|
||
if ( dPs < - EPS_SMALL)
|
||
ptP1 = ptPs ;
|
||
else if ( dPs > dLen + EPS_SMALL)
|
||
ptP2 = ptPs ;
|
||
// punto end rispetto al segmento
|
||
double dPe = ( ptPe - ptP1) * vtT ;
|
||
if ( dPe < - EPS_SMALL)
|
||
ptP1 = ptPe ;
|
||
else if ( dPe > dLen + EPS_SMALL)
|
||
ptP2 = ptPe ;
|
||
}
|
||
}
|
||
}
|
||
}
|
||
// se facce in contatto, calcolo l'angolo che formano ( 0 = allineate, > 0 convesse, < 0 concave)
|
||
if ( bAdjac) {
|
||
Vector3d vtN1 = m_vTria[m_vFacet[nF1]].vtN ;
|
||
Vector3d vtN2 = m_vTria[m_vFacet[nF2]].vtN ;
|
||
Vector3d vtCm = ptP2 - ptP1 ;
|
||
bool bDet ;
|
||
if ( ! vtN1.GetRotation( vtN2, vtCm, dAng, bDet) || ! bDet)
|
||
return false ;
|
||
}
|
||
return true ;
|
||
}
|
||
|
||
//----------------------------------------------------------------------------
|
||
bool
|
||
SurfTriMesh::GetFacetCenter( int nF, Point3d& ptCen, Vector3d& vtN) const
|
||
{
|
||
// recupero i loop della faccia
|
||
POLYLINEVECTOR vPL ;
|
||
if ( ! GetFacetLoops( nF, vPL) || vPL.empty())
|
||
return false ;
|
||
// calcolo il centro del loop esterno (è il primo)
|
||
PolygonPlane PolyPlane ;
|
||
Point3d ptP ;
|
||
for ( bool bFound = vPL[0].GetFirstPoint( ptP) ; bFound ; bFound = vPL[0].GetNextPoint( ptP))
|
||
PolyPlane.AddPoint( ptP) ;
|
||
if ( ! PolyPlane.GetCentroid( ptCen))
|
||
return false ;
|
||
// recupero la normale di un triangolo della faccetta
|
||
vtN = m_vTria[m_vFacet[nF]].vtN ;
|
||
return true ;
|
||
}
|
||
|
||
//----------------------------------------------------------------------------
|
||
bool
|
||
SurfTriMesh::GetFacetNormal( int nF, Vector3d& vtN) const
|
||
{
|
||
// la superficie deve essere validata
|
||
if ( m_nStatus != OK)
|
||
return false ;
|
||
// verifico stato sfaccettatura
|
||
if ( ! VerifyFaceting())
|
||
return false ;
|
||
// l'indice della faccia deve essere nei limiti
|
||
if ( nF < 0 || nF >= int( m_vFacet.size()))
|
||
return false ;
|
||
// recupero la normale di un triangolo della faccetta
|
||
vtN = m_vTria[m_vFacet[nF]].vtN ;
|
||
return true ;
|
||
}
|
||
|
||
//----------------------------------------------------------------------------
|
||
bool
|
||
SurfTriMesh::GetFacetArea( int nF, double& dArea) const
|
||
{
|
||
// recupero i loop della faccia
|
||
POLYLINEVECTOR vPL ;
|
||
if ( ! GetFacetLoops( nF, vPL) || vPL.empty())
|
||
return false ;
|
||
// sommo le aree dei diversi loop (quelli interni hanno area negativa)
|
||
dArea = 0 ;
|
||
for ( size_t i = 0 ; i < vPL.size() ; ++ i) {
|
||
PolygonPlane PolyPlane ;
|
||
Point3d ptP ;
|
||
for ( bool bFound = vPL[i].GetFirstPoint( ptP) ; bFound ; bFound = vPL[i].GetNextPoint( ptP))
|
||
PolyPlane.AddPoint( ptP) ;
|
||
double dLoopArea ;
|
||
if ( ! PolyPlane.GetArea( dLoopArea))
|
||
return false ;
|
||
if ( i == 0)
|
||
dArea += dLoopArea ;
|
||
else
|
||
dArea -= dLoopArea ;
|
||
}
|
||
return true ;
|
||
}
|
||
|
||
//----------------------------------------------------------------------------
|
||
SurfTriMesh*
|
||
SurfTriMesh::CloneFacet( int nF) const
|
||
{
|
||
if ( ! IsValid())
|
||
return nullptr ;
|
||
// alloco la nuova superficie
|
||
PtrOwner<SurfTriMesh> pSTM( CreateBasicSurfTriMesh()) ;
|
||
if ( IsNull( pSTM))
|
||
return nullptr ;
|
||
// recupero l'elenco dei triangoli della faccia
|
||
INTVECTOR vTria ;
|
||
if ( ! GetAllTriaInFacet( nF, vTria))
|
||
return nullptr ;
|
||
// preparo lo spazio nella nuova superficie
|
||
int nTria = int( vTria.size()) ;
|
||
if ( ! pSTM->Init( nTria + 2, nTria, 1))
|
||
return nullptr ;
|
||
// tabella hash per indici vertici (vecchi e nuovi indici)
|
||
typedef unordered_map<int,int> VVMAP ;
|
||
VVMAP PntMap( nTria + 2) ;
|
||
// inserisco i triangoli nella nuova superficie
|
||
for ( int i = 0 ; i < nTria ; ++ i) {
|
||
// indice del triangolo
|
||
int nT = vTria[i] ;
|
||
// ciclo sui tre vertici
|
||
int nIdV[3] ;
|
||
for ( int j = 0 ; j < 3 ; ++ j) {
|
||
// verifico se vertice già presente
|
||
const auto it = PntMap.find( m_vTria[nT].nIdVert[j]) ;
|
||
if ( it == PntMap.end()) {
|
||
// aggiungo il vertice
|
||
if ( ( nIdV[j] = pSTM->AddVertex( m_vVert[m_vTria[nT].nIdVert[j]].ptP)) == SVT_NULL)
|
||
return nullptr ;
|
||
PntMap.insert( VVMAP::value_type( m_vTria[nT].nIdVert[j], nIdV[j])) ;
|
||
}
|
||
else
|
||
nIdV[j] = it->second ;
|
||
}
|
||
// inserisco il triangolo
|
||
if ( pSTM->AddTriangle( nIdV) == SVT_NULL)
|
||
return nullptr ;
|
||
}
|
||
// sistemazioni finali della superficie
|
||
if ( ! pSTM->AdjustTopology())
|
||
return nullptr ;
|
||
// restituisco la superficie
|
||
return Release( pSTM) ;
|
||
}
|
||
|
||
//----------------------------------------------------------------------------
|
||
bool
|
||
SurfTriMesh::RemoveFacet( int nF)
|
||
{
|
||
// verifico esistenza delle due facce
|
||
if ( ! IsValid() || nF < 0 || nF >= GetFacetCount())
|
||
return false ;
|
||
|
||
// recupero l'elenco dei triangoli della faccia
|
||
INTVECTOR vTria ;
|
||
if ( ! GetAllTriaInFacet( nF, vTria))
|
||
return false ;
|
||
|
||
// Rimuovo i triangoli della faccia
|
||
for ( int i : vTria)
|
||
RemoveTriangle( i) ;
|
||
|
||
// aggiorno le strutture di base dell'oggetto
|
||
if ( ! DoCompacting())
|
||
return false ;
|
||
|
||
// dichiaro necessità ricalcolo della grafica e di hashgrids3d
|
||
m_OGrMgr.Reset() ;
|
||
ResetHashGrids3d() ;
|
||
|
||
return true ;
|
||
}
|
||
|
||
//----------------------------------------------------------------------------
|
||
bool
|
||
SurfTriMesh::SwapFacets( int nF1, int nF2)
|
||
{
|
||
// verifico esistenza delle due facce
|
||
if ( ! IsValid() || nF1 < 0 || nF1 >= GetFacetCount() || nF2 < 0 || nF2 >= GetFacetCount())
|
||
return false ;
|
||
|
||
// recupero gli indici dei triangoli della prima e seconda faccia
|
||
INTVECTOR vTria1 ;
|
||
if ( ! GetAllTriaInFacet( nF1, vTria1))
|
||
return false ;
|
||
INTVECTOR vTria2 ;
|
||
if ( ! GetAllTriaInFacet( nF2, vTria2))
|
||
return false ;
|
||
|
||
// scambio gli indici di faccia dei triangoli
|
||
for ( int i = 0 ; i < int( vTria1.size()) ; ++ i)
|
||
m_vTria[vTria1[i]].nIdFacet = nF2 ;
|
||
for ( int i = 0 ; i < int( vTria2.size()) ; ++ i)
|
||
m_vTria[vTria2[i]].nIdFacet = nF1 ;
|
||
|
||
// scambio i riferimenti al primo triangolo delle facce
|
||
swap( m_vFacet[nF1], m_vFacet[nF2]) ;
|
||
|
||
return true ;
|
||
}
|
||
|
||
//----------------------------------------------------------------------------
|
||
bool
|
||
SurfTriMesh::GetFacetLocalBBox( int nF, BBox3d& b3Loc, int nFlag) const
|
||
{
|
||
// verifico esistenza della faccia
|
||
if ( nF < 0 || nF >= GetFacetCount())
|
||
return false ;
|
||
|
||
// recupero gli indici dei triangoli della faccia
|
||
INTVECTOR vTria ;
|
||
if ( ! GetAllTriaInFacet( nF, vTria))
|
||
return false ;
|
||
|
||
// assegno il box in locale
|
||
b3Loc.Reset() ;
|
||
for ( int i = 0 ; i < int( vTria.size()) ; ++ i) {
|
||
// indice del triangolo
|
||
int nT = vTria[i] ;
|
||
// ciclo sui tre vertici
|
||
for ( int j = 0 ; j < 3 ; ++ j)
|
||
b3Loc.Add( m_vVert[m_vTria[nT].nIdVert[j]].ptP) ;
|
||
}
|
||
return true ;
|
||
}
|
||
|
||
//----------------------------------------------------------------------------
|
||
bool
|
||
SurfTriMesh::GetFacetBBox( int nF, const Frame3d& frRef, BBox3d& b3Ref, int nFlag) const
|
||
{
|
||
// verifico esistenza della faccia
|
||
if ( nF < 0 || nF >= GetFacetCount())
|
||
return false ;
|
||
|
||
// recupero gli indici dei triangoli della faccia
|
||
INTVECTOR vTria ;
|
||
if ( ! GetAllTriaInFacet( nF, vTria))
|
||
return false ;
|
||
|
||
// assegno il box nel riferimento
|
||
b3Ref.Reset() ;
|
||
for ( int i = 0 ; i < int( vTria.size()) ; ++ i) {
|
||
// indice del triangolo
|
||
int nT = vTria[i] ;
|
||
// ciclo sui tre vertici
|
||
for ( int j = 0 ; j < 3 ; ++ j) {
|
||
Point3d ptTemp = m_vVert[m_vTria[nT].nIdVert[j]].ptP ;
|
||
ptTemp.ToGlob( frRef) ;
|
||
b3Ref.Add( ptTemp) ;
|
||
}
|
||
}
|
||
return true ;
|
||
}
|
||
|
||
//----------------------------------------------------------------------------
|
||
bool
|
||
SurfTriMesh::VerifyFacetEdging( void) const
|
||
{
|
||
if ( m_bFacEdged)
|
||
return true ;
|
||
return (const_cast<SurfTriMesh*>(this))->UpdateFacetEdging() ;
|
||
}
|
||
|
||
//----------------------------------------------------------------------------
|
||
bool
|
||
SurfTriMesh::UpdateFacetEdging( void)
|
||
{
|
||
// reset degli Edge
|
||
m_bFacEdged = false ;
|
||
m_vFacEdge.clear() ;
|
||
|
||
// verifico validità sfaccettatura
|
||
if ( ! VerifyFaceting())
|
||
return false ;
|
||
|
||
// reset flag temporaneo di edge dei triangoli
|
||
for ( int i = 0 ; i < int( m_vTria.size()) ; ++ i) {
|
||
m_vTria[i].nETempFlag[0] = 0 ;
|
||
m_vTria[i].nETempFlag[1] = 0 ;
|
||
m_vTria[i].nETempFlag[2] = 0 ;
|
||
}
|
||
|
||
// ricostruisco gli edge delle facce
|
||
for ( int nT = 0 ; nT < int( m_vTria.size()) ; ++ nT) {
|
||
StmTria& Tria = m_vTria[nT] ;
|
||
// salto triangoli cancellati
|
||
if ( Tria.nIdVert[0] == SVT_DEL)
|
||
continue ;
|
||
// ciclo sui lati del triangolo
|
||
for ( int nE = 0 ; nE < 3 ; ++ nE) {
|
||
int nE2 = ( nE + 1) % 3 ;
|
||
int nTAdj = Tria.nIdAdjac[nE] ;
|
||
if ( Tria.nETempFlag[nE] == 0) {
|
||
if ( nTAdj == SVT_NULL) {
|
||
m_vFacEdge.emplace_back( Tria.nIdVert[nE], Tria.nIdVert[nE2], Tria.nIdFacet, SVT_NULL, 180.) ;
|
||
}
|
||
else if ( Tria.nIdFacet != m_vTria[nTAdj].nIdFacet) {
|
||
// calcolo l'angolo tra le facce ( 0 = allineate, > 0 convesse, < 0 concave)
|
||
Vector3d vtN1 = Tria.vtN ;
|
||
Vector3d vtN2 = m_vTria[nTAdj].vtN ;
|
||
Vector3d vtCm = m_vVert[Tria.nIdVert[nE2]].ptP - m_vVert[Tria.nIdVert[nE]].ptP ;
|
||
double dAng ;
|
||
bool bDet ;
|
||
if ( ! vtN1.GetRotation( vtN2, vtCm, dAng, bDet) || ! bDet)
|
||
dAng = 0 ;
|
||
// inserisco in lista
|
||
m_vFacEdge.emplace_back( Tria.nIdVert[nE], Tria.nIdVert[nE2], Tria.nIdFacet, m_vTria[nTAdj].nIdFacet, dAng) ;
|
||
// marco i due semi edge
|
||
int nEAdj = 0 ;
|
||
if ( m_vTria[nTAdj].nIdAdjac[1] == nT)
|
||
nEAdj = 1 ;
|
||
else if ( m_vTria[nTAdj].nIdAdjac[2] == nT)
|
||
nEAdj = 2 ;
|
||
Tria.nETempFlag[nE] = 1 ;
|
||
m_vTria[nTAdj].nETempFlag[nEAdj] = 1 ;
|
||
}
|
||
}
|
||
}
|
||
}
|
||
|
||
// calcolo bordo delle facce piane effettuato
|
||
m_bFacEdged = true ;
|
||
|
||
return true ;
|
||
}
|
||
|
||
//----------------------------------------------------------------------------
|
||
int
|
||
SurfTriMesh::GetEdgeCount( void) const
|
||
{
|
||
// la superficie deve essere validata
|
||
if ( m_nStatus != OK)
|
||
return -1 ;
|
||
// verifico stato bordi sfaccettatura
|
||
if ( ! VerifyFacetEdging())
|
||
return -1 ;
|
||
// restituisco il numero
|
||
return int( m_vFacEdge.size()) ;
|
||
}
|
||
|
||
//----------------------------------------------------------------------------
|
||
bool
|
||
SurfTriMesh::GetEdge( int nInd, int& nV1, int& nV2, int& nFl, int& nFr, double& dAng) const
|
||
{
|
||
// la superficie deve essere validata
|
||
if ( m_nStatus != OK)
|
||
return false ;
|
||
// verifico stato bordi sfaccettatura
|
||
if ( ! VerifyFacetEdging())
|
||
return false ;
|
||
// verifico la validità dell'indice
|
||
if ( nInd < 0 || nInd > int( m_vFacEdge.size()))
|
||
return SVT_NULL ;
|
||
// recupero i dati
|
||
nV1 = m_vFacEdge[nInd].nIdVert[0] ;
|
||
nV2 = m_vFacEdge[nInd].nIdVert[1] ;
|
||
nFl = m_vFacEdge[nInd].nIdFacAdj[0] ;
|
||
nFr = m_vFacEdge[nInd].nIdFacAdj[1] ;
|
||
dAng = m_vFacEdge[nInd].dIntAng ;
|
||
// ritorno indice edge corrente
|
||
return true ;
|
||
}
|
||
|
||
//----------------------------------------------------------------------------
|
||
bool
|
||
SurfTriMesh::GetEdge( int nInd, Point3d& ptP1, Point3d& ptP2, double& dAng) const
|
||
{
|
||
// la superficie deve essere validata
|
||
if ( m_nStatus != OK)
|
||
return false ;
|
||
// verifico stato bordi sfaccettatura
|
||
if ( ! VerifyFacetEdging())
|
||
return false ;
|
||
// verifico la validità dell'indice
|
||
if ( nInd < 0 || nInd > int( m_vFacEdge.size()))
|
||
return SVT_NULL ;
|
||
// recupero i dati
|
||
ptP1 = m_vVert[m_vFacEdge[nInd].nIdVert[0]].ptP ;
|
||
ptP2 = m_vVert[m_vFacEdge[nInd].nIdVert[1]].ptP ;
|
||
dAng = m_vFacEdge[nInd].dIntAng ;
|
||
// ritorno indice edge corrente
|
||
return true ;
|
||
}
|
||
|
||
//----------------------------------------------------------------------------
|
||
bool
|
||
SurfTriMesh::GetEdges( ICURVEPOVECTOR& vpCurve) const
|
||
{
|
||
// la superficie deve essere validata
|
||
if ( m_nStatus != OK)
|
||
return false ;
|
||
// verifico stato bordi sfaccettatura
|
||
if ( ! VerifyFacetEdging())
|
||
return false ;
|
||
// ciclo sugli edge, creo le curve e le assegno al parametro da restituire
|
||
for ( int nE = 0 ; nE < int( m_vFacEdge.size()) ; ++ nE) {
|
||
Point3d ptS, ptE ;
|
||
if ( ! GetVertex( m_vFacEdge[nE].nIdVert[0], ptS) ||
|
||
! GetVertex( m_vFacEdge[nE].nIdVert[1], ptE))
|
||
return false ;
|
||
PtrOwner<CurveLine> pLine( CreateBasicCurveLine()) ;
|
||
if ( IsNull( pLine) || ! pLine->Set( ptS, ptE))
|
||
continue ;
|
||
pLine->SetTempParam( m_vFacEdge[nE].dIntAng) ;
|
||
vpCurve.emplace_back( Release( pLine)) ;
|
||
}
|
||
return true ;
|
||
}
|
||
|
||
//-----------------------------------------------------------------------------
|
||
// Funzione per il calcolo della curvatura massima e minima in un vertice della superficie TriMesh
|
||
// vtMinK : versore direzione curvatura minima ( dMinK è il suo modulo con segno)
|
||
// vtMaxK : versore direzione curvatura massima ( dMaxK è il suo modulo con segno)
|
||
// bPlanar : Flag per indicare le superficie localmente piana
|
||
// vtNorm : versore normale del piano tangente alla supericie nel vertice
|
||
// dGaussK : valore della curvatura gaussiana con segno ( da associarsi a vtNorm)
|
||
// dMeanK : valore della curvatura media con segno ( da associarsi a vtNorm)
|
||
// dMaxNormK : valore della curvatura normale massima con segno ( da associarsi a vtNorm)
|
||
// dMinNormK : valore della curvatura normale minima con segno ( da associarsi a vtNorm)
|
||
bool
|
||
SurfTriMesh::GetCurvature( int nV,
|
||
double& dMinK, Vector3d& vtMinK, double& dMaxK, Vector3d& vtMaxK, bool& bPlanar, Vector3d& vtNorm,
|
||
double* dGaussK, double* dMeanK, double* dMaxNormK, double* dMinNormK) const
|
||
{
|
||
// Controllo la validità della TriMesh e del Vertice
|
||
if ( ! IsValid())
|
||
return false ;
|
||
Point3d ptCurr ;
|
||
if ( ! GetVertex( nV, ptCurr))
|
||
return false ;
|
||
bPlanar = false ;
|
||
|
||
// Recupero tutti i vertici attorno al vertice corrente ( se presenza di lato libero, allora errore)
|
||
INTVECTOR vT ;
|
||
bool bCirc ;
|
||
if ( ! GetAllTriaAroundVertex( nV, vT, bCirc) || ! bCirc)
|
||
return false ;
|
||
|
||
// Calcolo la normale del vertice pesata mediante angolo sotteso
|
||
// [Meyer et al. (2003) – Discrete Differential-Geometry Operators for Triangulated 2-Manifolds
|
||
// Section 5.4.2. Tip-Angle Weights.]
|
||
// [https://www.cs.cmu.edu/~kmcrane/Projects/DDG/paper.pdf]
|
||
INTSET setIndTriaNeightbors ;
|
||
for ( const int& nT : vT) {
|
||
// Recupero le direzioni delle semirette per l'angolo al vertice corrente
|
||
Triangle3d Tria ;
|
||
int nIdVert[3] ;
|
||
if ( ! GetTriangle( nT, Tria) || ! GetTriangle( nT, nIdVert))
|
||
return false ;
|
||
Vector3d vtDir0 = V_NULL, vtDir1 = V_NULL ;
|
||
for ( int i = 0 ; i < 3 ; ++ i) {
|
||
if ( AreSamePointApprox( Tria.GetP( i), ptCurr)) {
|
||
vtDir0 = Tria.GetP( ( i + 1) % 3) - Tria.GetP( i) ;
|
||
vtDir0.Normalize() ;
|
||
vtDir1 = Tria.GetP( ( i + 2) % 3) - Tria.GetP( i) ;
|
||
vtDir1.Normalize() ;
|
||
break ;
|
||
}
|
||
}
|
||
for ( int i = 0 ; i < 3 ; ++ i) {
|
||
if ( nV == nIdVert[i]) {
|
||
setIndTriaNeightbors.insert( nIdVert[( i + 1) % 3]) ;
|
||
setIndTriaNeightbors.insert( nIdVert[( i + 2) % 3]) ;
|
||
break ;
|
||
}
|
||
}
|
||
// Calcolo l'angolo sotteso
|
||
double dCosTheta = max( -1., min( 1., ( vtDir0 * vtDir1))) ;
|
||
double dTheta = acos( dCosTheta) ;
|
||
// Aggiorno il contributo della normale al vertice
|
||
vtNorm += ( dTheta * Tria.GetN()) ;
|
||
}
|
||
vtNorm.Normalize() ;
|
||
|
||
// [ESTIMATING CURVATURE ON TRIANGULAR MESHES, cap. 2.1. Fitting Methods,
|
||
// par. 2.1.2. Quadric Fitting]
|
||
// [https://people.eecs.berkeley.edu/~jrs/meshpapers/GatzkeGrimm.pdf]
|
||
|
||
// Definisco il piano tangente al vertice corrente
|
||
Plane3d plTan ;
|
||
if ( ! plTan.Set( ptCurr, vtNorm))
|
||
return false ;
|
||
|
||
// Proietto i punti Neightbors(1) sul piano tangete ( senza ripeterli)
|
||
BIPNTVECTOR vPtPtProj ; vPtPtProj.reserve( setIndTriaNeightbors.size()) ;
|
||
for ( auto nIter = setIndTriaNeightbors.begin() ; nIter != setIndTriaNeightbors.end() ; ++ nIter) {
|
||
Point3d ptNeightbors ;
|
||
if ( ! GetVertex( *nIter, ptNeightbors))
|
||
return false ;
|
||
vPtPtProj.emplace_back( make_pair( ptNeightbors, ProjectPointOnPlane( ptNeightbors, plTan))) ;
|
||
}
|
||
|
||
// Recupero due versori perpendicolari nel piano definito
|
||
Vector3d vtTan1 = V_NULL ;
|
||
if ( abs( vtNorm.x) < 1./64. && abs( vtNorm.y) < 1./64.)
|
||
vtTan1 = Y_AX ^ vtNorm ;
|
||
else
|
||
vtTan1 = Z_AX ^ vtNorm ;
|
||
vtTan1.Normalize() ;
|
||
Vector3d vtTan2 = vtNorm ^ vtTan1 ;
|
||
vtTan2.Normalize() ;
|
||
|
||
// Sistema da risolvere mediante minimi quadrati : z(u,v) = Au^2 + Buv + Cv^2
|
||
// Questo sistema è molto più semplice e robusto rispetto a z(u,v) = Au^2 + Buv + Cv^2 + Du + Ev + F
|
||
// per il fatto che ora le coordinate sono in locale al piano tangente alla superficie ( mettendo
|
||
// quindi a 0 i coefficienti del primo ordine e il termine noto F)
|
||
// Definizione della matrice A(Nx3)
|
||
const int nPts = int( vPtPtProj.size()) ;
|
||
Eigen::MatrixXd mat_A( nPts, 3) ;
|
||
Eigen::VectorXd vec_b( nPts) ;
|
||
for ( int i = 0 ; i < nPts ; ++ i) {
|
||
Vector3d vtEdgeProj = vPtPtProj[i].second - ptCurr ;
|
||
Vector3d vtEdge = ( vPtPtProj[i].first - vPtPtProj[i].second) ;
|
||
double dU = vtEdgeProj * vtTan1 ;
|
||
double dV = vtEdgeProj * vtTan2 ;
|
||
double dW = vtEdge * vtNorm ;
|
||
mat_A( i, 0) = dU * dU ;
|
||
mat_A( i, 1) = dU * dV ;
|
||
mat_A( i, 2) = dV * dV ;
|
||
vec_b( i) = dW ;
|
||
}
|
||
// Risoluzione del sistema
|
||
Eigen::VectorXd vec_x = mat_A.colPivHouseholderQr().solve( vec_b) ;
|
||
|
||
// Costruzione della matrice hessiana
|
||
Eigen::Matrix2d mat_H ;
|
||
mat_H << 2. * vec_x( 0), vec_x( 1),
|
||
vec_x( 1), 2. * vec_x( 2) ;
|
||
|
||
// Calcolo gli autovalori ( quindi le curvature principali) e restituisco i risultati
|
||
Eigen::SelfAdjointEigenSolver<Eigen::Matrix2d> Solver( mat_H) ;
|
||
dMinK = Solver.eigenvalues()( 0) ;
|
||
dMaxK = Solver.eigenvalues()( 1) ;
|
||
Eigen::Vector2d dir_Min = Solver.eigenvectors().col( 0) ;
|
||
Eigen::Vector2d dir_Max = Solver.eigenvectors().col( 1) ;
|
||
vtMinK = dir_Min( 0) * vtTan1 + dir_Min( 1) * vtTan2 ;
|
||
vtMaxK = dir_Max( 0) * vtTan1 + dir_Max( 1) * vtTan2 ;
|
||
bPlanar = ( abs( dMinK) < EPS_SMALL && abs( dMaxK) < EPS_SMALL) ;
|
||
|
||
// Se rischiesto calcolo della curvatura Gaussiana
|
||
if ( dGaussK != nullptr)
|
||
*dGaussK = ( dMinK * dMaxK) ;
|
||
|
||
// Se richiesto calcolo della curvatura Media
|
||
if ( dMeanK != nullptr)
|
||
*dMeanK = ( dMinK + dMaxK) / 2. ;
|
||
|
||
// Se richiesto calcolo della curvatura normale lungo la direzione massima
|
||
if ( dMaxNormK != nullptr)
|
||
*dMaxNormK = dir_Max.transpose() * mat_H * dir_Max ;
|
||
|
||
// Se richiesto calcolo della curvatura normale lungo la direzione minima
|
||
if ( dMinNormK != nullptr)
|
||
*dMinNormK = dir_Min.transpose() * mat_H * dir_Min ;
|
||
|
||
return true ;
|
||
}
|