Egalware/EgtGeomKernel
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity
Files
d92344f2bbc8da76b1e9be342d33cf2fd2d9f9e2
EgtGeomKernel/SurfTriMeshFaceting.cpp
T
Riccardo Elitropi d62d6946c3 EgtGeomKernel : 
- miglioramento della funzione per calcolo curvatura di una superficie TriMesh.
2025-11-06 13:04:40 +01:00

1166 lines
41 KiB
C++
Raw Blame History

//----------------------------------------------------------------------------
// 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
// dMinK : curvatura minima
// vtMinK : versore direzione curvatura minima
// dMaxK : curvatura massima
// vtMaxK : versore direzione curvatura massima
// bPlanar : Flag per indicare le superficie localmente piana
// vtNorm : versore normale del piano tangente alla supericie nel vertice
bool
SurfTriMesh::GetCurvature( int nV,
double& dMinK, Vector3d& vtMinK, double& dMaxK, Vector3d& vtMaxK, bool& bPlanar, Vector3d& vtNorm) 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 e distanza baricentrica
// ["Estimation Normal Vector of Triangular Mesh Vertex by Angle and Centroid Weights"]
// Controllo anche di non essere in presenza di uno spigolo vivo
INTSET setIndTriaNeightbors ;
bool bFirstTria = true ;
Vector3d vtNFirstTria = V_NULL ;
for ( const int& nT : vT) {
// Recupero il triangolo corrente
Triangle3d Tria ;
int nIdVert[3] ;
if ( ! GetTriangle( nT, Tria) || ! GetTriangle( nT, nIdVert))
return false ;
// Se il triangolo ha area troppo piccola, lo scarto
if ( Tria.GetArea() < EPS_SMALL)
continue ;
// Se primo triangolo salvo la sua normale di riferimento per successivi confronti
if ( bFirstTria) {
vtNFirstTria = Tria.GetN() ;
bFirstTria = false ;
}
else {
// Se spigolo vivo, la curvatura non può esistere
if ( Tria.GetN() * vtNFirstTria < m_dCosSmAng - EPS_SMALL)
return false ;
}
// Recupero le direzioni delle semirette per l'angolo al vertice corrente
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) ;
// Calcolo la distanza baricentrica
double dSqBarDist = max( EPS_SMALL, ( ptCurr - Tria.GetCentroid()).SqLen()) ;
// Aggiorno il contributo della normale al vertice
vtNorm += ( ( dTheta / dSqBarDist) * 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 {{ 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) ;
return true ;
}
Reference in New Issue View Git Blame Copy Permalink