Egalware/EgtGeomKernel
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity
Files
master
EgtGeomKernel/Triangulate.cpp
T
Riccardo Elitropi 626d5b0e51 EgtGeomKernel 2.7f2 : 
- Aggiunte funzioni per calcolo di Offset per superfici chiuse TriMesh
- Piccola miglioria alla triangolazione (con SaraP)
- Migliorie per rimozioni TJunction, calcolo delle normali dei triangoli e creazione di una TriMesh a partire da uno ZMap (con SaraP)
- Aggiunte funzioni di SubtractMap e piccole modifiche per estensione dei Box di creazione per gli Zmap.
2025-06-16 11:34:23 +02:00

1647 lines
60 KiB
C++
Raw Permalink Blame History

//----------------------------------------------------------------------------
// EgalTech 2014-2021
//----------------------------------------------------------------------------
// File : Triangulate.cpp Data : 18.12.21 Versione : 2.3l
// Contenuto : Implementazione della classe Triangulate.
//
//
//
// Modifiche : 08.04.14 DS Creazione modulo.
// 23.06.14 DS Corretto errore con contorni CW.
//
//----------------------------------------------------------------------------
//--------------------------- Include ----------------------------------------
#include "stdafx.h"
#include "DllMain.h"
#include "Triangulate.h"
#include "ProjPlane.h"
#include "earcut.hpp"
#include "/EgtDev/Include/EGkSfrCreate.h"
#include "/EgtDev/Include/EGkPolyLine.h"
#include "/EgtDev/Include/EGkPlane3d.h"
#include "/EgtDev/Include/EGkStringUtils3d.h"
#include "/EgtDev/Include/EgtNumUtils.h"
#include "/EgtDev/Extern/fist/Include/api_fist.h"
#include <algorithm>
using namespace std ;
//----------------------------------------------------------------------------
enum EarStatus{ EAS_NULL = -1, EAS_NO = 0, EAS_OK = 1} ;
enum TrgType { TRG_STD = 0, TRG_CAP = 1, TRG_NEEDLE = 2} ;
//----------------------------------------------------------------------------
static bool ChangeStartPntVector( int nNewStart, PNTVECTOR& vPi) ;
static bool MakeByFist( const POLYLINEVECTOR& vPL, PNTVECTOR& vPt, INTVECTOR& vTr) ;
static bool RemoveFistInvalidTrg( PNTVECTOR& vPt, INTVECTOR& vTr) ;
static int CalcTriangleType( const PNTVECTOR& vPt, const INTVECTOR& vTr, int nTrg) ;
static bool FindAdjacentOnLongerEdge( PNTVECTOR& vPt, INTVECTOR& vTr, int nTrgA, int&nEdgeA, int& nTrgB, int& nEdgeB) ;
static bool FlipTrg( INTVECTOR& vTr, int nTrgA, int nTrgB, int nEdgeA, int nEdgeB) ;
static bool TestAdjacentOnEdge( INTVECTOR& vTr, int nTrgA, int nEdgeA, int nTrgTest1, int nTrgTest2, int& nTrgB, int& nEdgeB) ;
//----------------------------------------------------------------------------
static bool FORCE_EARCUT_HPP = false ;
static bool FORCE_FIST = false ;
//----------------------------------------------------------------------------
// In : PolyLine
// Out : PNTVECTOR (Point3d Vector) : points of the polyline
// INTVECTOR (int Vector) : 3*T indices of above points for T triangles
//----------------------------------------------------------------------------
bool
Triangulate::Make( const PolyLine& PL, PNTVECTOR& vPt, INTVECTOR& vTr)
{
// pulisco i vettori di ritorno
vPt.clear() ;
vTr.clear() ;
// verifico che la polilinea contenga almeno 4 punti (primo e ultimo coincidenti)
if ( PL.GetPointNbr() < 4 || ! PL.IsClosed())
return false ;
// se fist ( e geometria più complessa di un quadrilatero)
if ( FORCE_FIST && PL.GetPointNbr() > 5)
return MakeByFist( {PL}, vPt, vTr) ;
// verifico che la polilinea sia chiusa e piana e calcolo il piano medio del poligono
double dArea ;
Plane3d plPlane ;
if ( ! PL.IsClosedAndFlat( plPlane, dArea, 50 * EPS_SMALL))
return false ;
// determino il piano ottimale di proiezione e il relativo senso di rotazione
bool bCCW ;
if ( ! CalcProjPlane( plPlane.GetVersN(), m_nPlane, bCCW))
return false ;
// forzo esecuzione triangolazione con earcut.hpp
if ( FORCE_EARCUT_HPP) {
return MakeByEC_HPP( PL, bCCW, vPt, vTr) ;
}
// riempio il vettore con i vertici del poligono da triangolare
vPt.reserve( PL.GetPointNbr() - 1) ;
// salto il primo punto (coincide con l'ultimo)
Point3d ptP ;
if ( ! PL.GetFirstPoint( ptP))
return false ;
// inserisco i punti
while ( PL.GetNextPoint( ptP))
vPt.push_back( ptP) ;
// se non CCW inverto il vettore dei punti
if ( ! bCCW)
reverse( vPt.begin(), vPt.end()) ;
// eseguo la triangolazione
double dMinMinAng2, dMinMinAng3 ;
if ( ! MakeByEC2( vPt, vTr, dMinMinAng2)) {
if ( ! MakeByEC3( vPt, vTr, dMinMinAng3)) {
LOG_INFO( GetEGkLogger(), "Info : problems in MakeByEC23(1)")
return MakeByEC_HPP( PL, bCCW, vPt, vTr) ;
}
}
else if ( dMinMinAng2 < 5) {
INTVECTOR vMyTr ;
if ( MakeByEC3( vPt, vMyTr, dMinMinAng3) && dMinMinAng3 > dMinMinAng2)
vTr = vMyTr ;
}
// se era CW, devo invertire il senso dei triangoli
if ( ! bCCW)
reverse( vTr.begin(), vTr.end()) ;
return true ;
}
//----------------------------------------------------------------------------
// In : POLYLINEVECTOR : vector of polylines, the first outer, the others inner
// Out : PNTVECTOR (Point3d Vector) : points of the polyline
// INTVECTOR (int Vector) : 3*T indices of above points for T triangles
//----------------------------------------------------------------------------
bool
Triangulate::Make( const POLYLINEVECTOR& vPL, PNTVECTOR& vPt, INTVECTOR& vTr)
{
// pulisco i vettori di ritorno
vPt.clear() ;
vTr.clear() ;
// ci devono essere delle polilinee nel vettore
if ( &vPL == nullptr || vPL.empty())
return false ;
// se una sola polilinea mi riconduco al caso precedente
if ( vPL.size() == 1)
return Make( vPL[0], vPt, vTr) ;
// se fist
if ( FORCE_FIST)
return MakeByFist( vPL, vPt, vTr) ;
// verifico che la polilinea esterna sia chiusa e piana e calcolo il piano medio del poligono
double dArea ;
Plane3d plExtPlane ;
if ( ! vPL[0].IsClosedAndFlat( plExtPlane, dArea, 50 * EPS_SMALL))
return false ;
// determino il piano ottimale di proiezione e il relativo senso di rotazione
bool bCCW ;
if ( ! CalcProjPlane( plExtPlane.GetVersN(), m_nPlane, bCCW))
return false ;
// verifico le altre polilinee
for ( int i = 1 ; i < int( vPL.size()) ; ++i) {
// deve essere chiusa, giacere nello stesso piano ed essere orientata al contrario della esterna
double dArea2 ;
Plane3d plPlane ;
if ( ! vPL[i].IsClosedAndFlat( plPlane, dArea2, 50 * EPS_SMALL) ||
! AreOppositeVectorApprox( plExtPlane.GetVersN(), plPlane.GetVersN()))
return false ;
}
// forzo esecuzione triangolazione con earcut.hpp
if ( FORCE_EARCUT_HPP) {
return MakeByEC_HPP( vPL, bCCW, vPt, vTr) ;
}
// se non CCW inverto tutte le polilinee
if ( ! bCCW) {
for ( int i = 0 ; i < int( vPL.size()) ; ++i)
const_cast<PolyLine&>(vPL[i]).Invert( true) ;
}
// calcolo ordine decrescente su Xmax o Ymax o Zmax secondo m_nPlane per le curve interne
INTVECTOR vOrd ;
bool bOk = SortInternalLoops( vPL, vOrd) ;
// riempio il vettore con i punti dei poligoni da triangolare
// calcolo spazio totale e spazio massimo per un percorso interno
int nTot = vPL[0].GetPointNbr() - 1 ;
int nMax = 0 ;
for ( int i = 1 ; i < int( vPL.size()) ; ++i) {
nTot += vPL[i].GetPointNbr() + 1 ;
if ( vPL[i].GetPointNbr() > nMax)
nMax = vPL[i].GetPointNbr() ;
}
// riservo spazio pari al totale dei punti (anche quelli ripetuti di chiusura che servono per i link)
vPt.reserve( nTot) ;
// inserisco i punti del contorno esterno (tranne il primo che coincide con l'ultimo)
bOk = bOk && GetPntVectorFromPolyline( vPL[0], false, vPt) ;
// ciclo sui percorsi interni ordinati secondo Xmax decrescente
PNTVECTOR vPi ;
vPi.reserve( nMax) ;
for ( int i = 1 ; bOk && i < int( vPL.size()) ; ++ i) {
// riordino il percorso per avere punto iniziale con Xmax
if ( ! GetPntVectorFromPolyline( vPL[vOrd[i]], true, vPi)) {
bOk = false ;
break ;
}
// cerco un punto del percorso esterno visibile dal punto iniziale del precedente interno
int nI ;
if ( ! GetOuterPntToJoin( vPt, vPi[0], nI)) {
bOk = false ;
break ;
}
// riordino percorso esterno per avere questo punto all'inizio
if ( ! ChangeStartPntVector( nI, vPt)) {
bOk = false ;
break ;
}
// ripeto punto iniziale
vPt.push_back( vPt[0]) ;
// accodo percorso interno
for ( int j = 0 ; j < int( vPi.size()) ; ++j)
vPt.push_back( vPi[j]) ;
// ripeto punto iniziale di percorso interno
vPt.push_back( vPi[0]) ;
// cancello percorso interno
vPi.clear() ;
}
// ripristino l'orientamento delle polilinee originali per il caso non CCW
if ( ! bCCW) {
for ( int i = 0 ; i < int( vPL.size()) ; ++i)
const_cast<PolyLine&>(vPL[i]).Invert( true) ;
}
if ( ! bOk)
return MakeByEC_HPP( vPL, bCCW, vPt, vTr) ;
// eseguo la triangolazione
double dMinMinAng2, dMinMinAng3 ;
if ( ! MakeByEC2( vPt, vTr, dMinMinAng2)) {
if ( ! MakeByEC3( vPt, vTr, dMinMinAng3)) {
LOG_INFO( GetEGkLogger(), "Info : problems in MakeByEC23(N)")
return MakeByEC_HPP( vPL, bCCW, vPt, vTr) ;
}
}
else if ( dMinMinAng2 < 5) {
INTVECTOR vMyTr ;
if ( MakeByEC3( vPt, vMyTr, dMinMinAng3) && dMinMinAng3 > dMinMinAng2)
vTr = vMyTr ;
}
// se era CW, devo invertire il senso dei triangoli
if ( ! bCCW)
reverse( vTr.begin(), vTr.end()) ;
return true ;
}
//---------------------------------------------------------------------------
bool
Triangulate::MakeAdvanced( const POLYLINEVECTOR& vPLORIG, PNTVECTOR& vPt, INTVECTOR& vTr)
{
INTMATRIX vnPLIndMat ;
return MakeAdvanced( vPLORIG, vPt, vTr, vnPLIndMat) ;
}
//---------------------------------------------------------------------------
bool
Triangulate::MakeAdvanced( const POLYLINEVECTOR& vPLORIG, PNTVECTOR& vPt, INTVECTOR& vTr, const INTMATRIX& vnPLIndMatPre)
{
vPt.clear() ;
vTr.clear() ;
// se non ho PolyLine, allora non faccio nulla
if ( vPLORIG.empty())
return true ;
POLYLINEVECTOR vPL ;
Vector3d vtN ;
// se non sono stati passate le info per ordinare le polyline allora le ordino
INTMATRIX vnPLIndMat ;
BOOLVECTOR vbInvert ;
if ( vnPLIndMatPre.empty()) {
if ( ! CalcRegionPolyLines( vPLORIG, vtN, vnPLIndMat, vbInvert))
return false ;
vPL = vPLORIG ;
for ( int i = 0 ; i < int( vbInvert.size()) ; i++) {
if ( vbInvert[i])
vPL[i].Invert() ;
}
}
else {
// ho già calcolato e riordinato tutto, devo solo fare una copia delle polyline
// non serve fare le eventuali inversioni delle polyline, perché se è già stata calcolata la matrice dei chunck allora sono già state invertire
vPL = vPLORIG ;
vnPLIndMat = vnPLIndMatPre ;
}
// chiamo la Triangolazione per ogni riga della matrice ( quindi su ogni "Chunk")
for ( int i = 0 ; i < int( vnPLIndMat.size()) ; ++ i) {
PNTVECTOR vPt_tmp ; INTVECTOR vTr_tmp ;
POLYLINEVECTOR vPL_tmp ;
for ( int j = 0 ; j < int( vnPLIndMat[i].size()) ; ++ j)
vPL_tmp.push_back( vPL[vnPLIndMat[i][j]]) ;
if ( ! Make( vPL_tmp, vPt_tmp, vTr_tmp))
return false ;
int nSize = int( vPt.size()) ;
for ( int p = 0 ; p < int( vPt_tmp.size()) ; ++ p)
vPt.push_back( vPt_tmp[p]) ;
for ( int t = 0 ; t < int( vTr_tmp.size()) ; ++ t)
vTr.push_back( nSize + vTr_tmp[t]) ;
}
return true ;
}
//----------------------------------------------------------------------------
// Triangulate the CCW n-gon specified by the vertices vPt (Pt[n] != Pt[0])
// Ear Clipping algorithm from mapbox
//----------------------------------------------------------------------------
bool
Triangulate::MakeByEC_HPP( const PolyLine& PL, bool bCCW, PNTVECTOR& vPt, INTVECTOR& vTr)
{
// pulisco i vettori di ritorno
vPt.clear() ;
vTr.clear() ;
// riempio il vettore con i vertici del poligono da triangolare
vPt.reserve( PL.GetPointNbr() - 1) ;
// salto il primo punto (coincide con l'ultimo)
Point3d ptP ;
if ( ( bCCW && ! PL.GetFirstPoint( ptP)) ||
( ! bCCW && ! PL.GetLastPoint( ptP)))
return false ;
// inserisco i punti
while ( ( bCCW && PL.GetNextPoint( ptP)) ||
( ! bCCW && PL.GetPrevPoint( ptP)))
vPt.push_back( ptP) ;
// Create array
using Point = array<double, 2> ;
vector<vector<Point>> polygon ;
// Fill polygon structure with actual data. Any winding order works.
// The first polyline defines the main polygon.
polygon.resize( 1) ;
polygon[0].reserve( vPt.size()) ;
for ( int i = 0 ; i < int( vPt.size()) ; ++i) {
switch ( m_nPlane) {
default : // PL_XY
polygon[0].push_back( { vPt[i].x, vPt[i].y}) ;
break ;
case PL_YZ :
polygon[0].push_back( { vPt[i].y, vPt[i].z}) ;
break ;
case PL_ZX :
polygon[0].push_back( { vPt[i].z, vPt[i].x}) ;
break ;
}
}
// Run tessellation
// Returns array of indices that refer to the vertices of the input polygon.
// Three subsequent indices form a triangle. Output triangles are counterclockwise.
vTr = mapbox::earcut<int32_t>( polygon) ;
if ( vTr.empty()) {
LOG_INFO( GetEGkLogger(), "Info : problems in MakeByEC_HPP(1)")
return false ;
}
if ( ! bCCW)
reverse( vTr.begin(), vTr.end()) ;
return true ;
}
//----------------------------------------------------------------------------
bool
Triangulate::MakeByEC_HPP( const POLYLINEVECTOR& vPL, bool bCCW, PNTVECTOR& vPt, INTVECTOR& vTr)
{
// pulisco i vettori di ritorno
vPt.clear() ;
vTr.clear() ;
// Creo il vettore con i punti dei poligoni da triangolare
// calcolo spazio totale e spazio massimo per un percorso interno
int nTot = vPL[0].GetPointNbr() - 1 ;
for ( int i = 1 ; i < int( vPL.size()) ; ++i)
nTot += vPL[i].GetPointNbr() - 1 ;
// riservo spazio pari al totale dei punti (esclusi quelli ripetuti)
vPt.reserve( nTot) ;
// Create array
using Point = array<double, 2> ;
vector<vector<Point>> polygon ;
polygon.resize( vPL.size()) ;
// Eseguo riempimento di vettore e array
for ( int i = 0 ; i < int( vPL.size()) ; ++i) {
polygon[i].reserve( vPL[i].GetPointNbr() - 1) ;
// salto il primo punto (coincide con l'ultimo)
Point3d ptP ;
if ( ( bCCW && ! vPL[i].GetFirstPoint( ptP)) ||
( ! bCCW && ! vPL[i].GetLastPoint( ptP)))
return false ;
// inserisco i punti
while ( ( bCCW && vPL[i].GetNextPoint( ptP)) ||
( ! bCCW && vPL[i].GetPrevPoint( ptP))) {
vPt.push_back( ptP) ;
switch ( m_nPlane) {
default : // PL_XY
polygon[i].push_back( { ptP.x, ptP.y}) ;
break ;
case PL_YZ :
polygon[i].push_back( { ptP.y, ptP.z}) ;
break ;
case PL_ZX :
polygon[i].push_back( { ptP.z, ptP.x}) ;
break ;
}
}
}
// Run tessellation
// Returns array of indices that refer to the vertices of the input polygon.
// Three subsequent indices form a triangle. Output triangles are counterclockwise.
vTr = mapbox::earcut<int32_t>( polygon) ;
if ( vTr.empty()) {
LOG_INFO( GetEGkLogger(), "Info : problems in MakeByEC_HPP(N)")
return false ;
}
if ( ! bCCW)
reverse( vTr.begin(), vTr.end()) ;
return true ;
}
//----------------------------------------------------------------------------
// Triangulate the CCW n-gon specified by the vertices vPt (Pt[n] != Pt[0])
// Ear Clipping algorithm
//----------------------------------------------------------------------------
bool
Triangulate::MakeByEC( const PNTVECTOR& vPt, INTVECTOR& vTr)
{
// Clear triangle vector
vTr.clear() ;
// At least 3 points
int n = int( vPt.size()) ;
if ( n < 3)
return false ;
// Creo il vettore degli indici del Poligono
INTVECTOR vPol ;
vPol.reserve( n) ;
for ( int i = 0 ; i < n ; ++ i)
vPol.push_back( i) ;
// Preallocate triangle vector ( #triangles = n - 2)
vTr.reserve( 3 * ( n - 2)) ;
// Set up previous and next links to effectively form a double-linked vertex list
INTVECTOR vPrev( n) ;
INTVECTOR vNext( n) ;
for ( int j = 0 ; j < n ; ++ j) {
vPrev[j] = j - 1 ;
vNext[j] = j + 1 ;
}
vPrev[0] = n - 1 ;
vNext[n-1] = 0 ;
// Start at vertex 0
int i = 0 ;
int nCount = n ;
// Keep removing vertices until just a triangle left
while ( n >= 3) {
// To avoid infinite loop
if ( nCount <= 0)
return false ;
-- nCount ;
// Test if current vertex, v[i], is an ear
bool bIsEar = TestTriangle( vPt, vPol, vPrev, vNext, i) ;
bool bSave = bIsEar ;
// Verify if triangle is null (accept to discard)
if ( ! bIsEar) {
if ( Collinear( vPt[vPol[vPrev[i]]], vPt[vPol[i]], vPt[vPol[vNext[i]]]))
bIsEar = true ;
}
// If current vertex v[i] is an ear, delete it and visit the previous vertex
if ( bIsEar) {
// Triangle (v[prev[i]], v[i], v[next[i]]) is an ear
if ( bSave) {
vTr.push_back( vPol[vPrev[i]]) ;
vTr.push_back( vPol[i]) ;
vTr.push_back( vPol[vNext[i]]) ;
}
// Delete vertex v[i] by redirecting next and previous links
// of neighboring verts past it. Decrement vertex count
vNext[vPrev[i]] = vNext[i] ;
vPrev[vNext[i]] = vPrev[i] ;
n--;
// Visit the previous vertex next
i = vPrev[i] ;
// Reset Count
nCount = n ;
}
else {
// Current vertex is not an ear; visit the next vertex
i = vNext[i] ;
}
}
return true ;
}
//----------------------------------------------------------------------------
// Triangulate the CCW n-gon specified by the vertices vPt (Pt[n] != Pt[0])
// Ear Clipping algorithm enhanced, choose smaller diagonal
//----------------------------------------------------------------------------
bool
Triangulate::MakeByEC2( const PNTVECTOR& vPt, INTVECTOR& vTr, double& dMinMinAng)
{
// Clear triangle vector
vTr.clear() ;
dMinMinAng = ANG_STRAIGHT ;
// At least 3 points
int n = int( vPt.size()) ;
if ( n < 3)
return false ;
// Creo il vettore degli indici del Poligono
INTVECTOR vPol ;
vPol.reserve( n) ;
for ( int i = 0 ; i < n ; ++ i)
vPol.push_back( i) ;
// Preallocate triangle vector ( #triangles = n - 2)
vTr.reserve( 3 * ( n - 2)) ;
// Set up previous and next links to effectively form a double-linked vertex list
INTVECTOR vPrev( n) ;
INTVECTOR vNext( n) ;
for ( int j = 0 ; j < n ; ++ j) {
vPrev[j] = j - 1 ;
vNext[j] = j + 1 ;
}
vPrev[0] = n - 1 ;
vNext[n-1] = 0 ;
// Prepare Ear status vector
INTVECTOR vEar( n, EAS_NULL) ;
// Prepare PointGrid
if ( ! PrepareGrid( vPt, vPol, vPrev, vNext))
return false ;
// Start at vertex 0
int i = 0 ;
int nCount = n ;
// Keep removing vertices until just a triangle left
while ( n >= 3) {
// To avoid infinite loop
if ( nCount <= 0)
return false ;
-- nCount ;
// Test if current vertex, v[i], is an ear
if ( vEar[i] == EAS_NULL)
vEar[i] = ( TestTriangle( vPt, vPol, vPrev, vNext, i) ? EAS_OK : EAS_NO) ;
bool bIsEar = ( vEar[i] == EAS_OK) ;
bool bSave = bIsEar ;
if ( bIsEar) {
// Save square distance of diagonal
double dSqDist = SqDist(vPt[vPol[vPrev[i]]], vPt[vPol[vNext[i]]]) ;
// Start point for trials
int j = i ;
int k = i ;
// Try with 19 next
int nLimN = min( 19, n / 2) ;
for ( int h = 0 ; h < nLimN ; ++ h) {
j = vNext[j] ;
double dMySqDist = SqDist( vPt[vPol[vPrev[j]]], vPt[vPol[vNext[j]]]) ;
if ( dMySqDist < 0.9 * dSqDist) {
if ( vEar[j] == EAS_NULL)
vEar[j] = ( TestTriangle( vPt, vPol, vPrev, vNext, j) ? EAS_OK : EAS_NO) ;
if ( vEar[j] == EAS_OK) {
dSqDist = dMySqDist ;
i = j ;
}
}
}
// Try with 19 prev
int nLimP = min( 19, n / 2) ;
for ( int h = 0 ; h < nLimP ; ++ h) {
k = vPrev[k] ;
double dMySqDist = SqDist( vPt[vPol[vPrev[k]]], vPt[vPol[vNext[k]]]) ;
if ( dMySqDist < 0.9 * dSqDist) {
if ( vEar[k] == EAS_NULL)
vEar[k] = ( TestTriangle( vPt, vPol, vPrev, vNext, k) ? EAS_OK : EAS_NO) ;
if ( vEar[k] == EAS_OK) {
dSqDist = dMySqDist ;
i = k ;
}
}
}
}
// Verify if triangle is null (accept to discard)
else {
if ( Collinear( vPt[vPol[vPrev[i]]], vPt[vPol[i]], vPt[vPol[vNext[i]]]) &&
! SameDirection( vPt[vPol[vPrev[i]]], vPt[vPol[i]], vPt[vPol[vNext[i]]]))
bIsEar = true ;
}
// If current vertex v[i] is an ear, delete it and visit the previous vertex
if ( bIsEar) {
// Triangle (v[prev[i]], v[i], v[next[i]]) is an ear
if ( bSave) {
vTr.push_back( vPol[vPrev[i]]) ;
vTr.push_back( vPol[i]) ;
vTr.push_back( vPol[vNext[i]]) ;
double dMinAng = CalcTriangleMinAngle( vPt[vPol[vPrev[i]]], vPt[vPol[i]], vPt[vPol[vNext[i]]]) ;
dMinMinAng = min( dMinMinAng, dMinAng) ;
}
// Reset earity of diagonal endpoints
vEar[vPrev[i]] = EAS_NULL ;
vEar[vNext[i]] = EAS_NULL ;
// Delete vertex v[i] by redirecting next and previous links
// of neighboring verts past it. Decrement vertex count
vNext[vPrev[i]] = vNext[i] ;
vPrev[vNext[i]] = vPrev[i] ;
-- n ;
// Visit the previous vertex next
i = vPrev[i] ;
// Reset Count
nCount = n ;
}
else {
// Current vertex is not an ear; visit the next vertex
i = vNext[i] ;
}
}
return true ;
}
//----------------------------------------------------------------------------
// Triangulate the CCW n-gon specified by the vertices vPt (Pt[n] != Pt[0])
// Ear Clipping algorithm enhanced, choose max MinAng
//----------------------------------------------------------------------------
bool
Triangulate::MakeByEC3( const PNTVECTOR& vPt, INTVECTOR& vTr, double& dMinMinAng)
{
// Clear triangle vector
vTr.clear() ;
dMinMinAng = ANG_STRAIGHT ;
// At least 3 points
int n = int( vPt.size()) ;
if ( n < 3)
return false ;
// Creo il vettore degli indici del Poligono
INTVECTOR vPol ;
vPol.reserve( n) ;
for ( int i = 0 ; i < n ; ++ i)
vPol.push_back( i) ;
// Preallocate triangle vector ( #triangles = n - 2)
vTr.reserve( 3 * ( n - 2)) ;
// Set up previous and next links to effectively form a double-linked vertex list
INTVECTOR vPrev( n) ;
INTVECTOR vNext( n) ;
for ( int j = 0 ; j < n ; ++ j) {
vPrev[j] = j - 1 ;
vNext[j] = j + 1 ;
}
vPrev[0] = n - 1 ;
vNext[n-1] = 0 ;
// Prepare Ear status vector
INTVECTOR vEar( n, EAS_NULL) ;
// Prepare PointGrid
if ( ! PrepareGrid( vPt, vPol, vPrev, vNext))
return false ;
// Start at vertex 0
int i = 0 ;
int nCount = n ;
// Keep removing vertices until just a triangle left
while ( n >= 3) {
// To avoid infinite loop
if ( nCount <= 0)
return false ;
-- nCount ;
// Test if current vertex, v[i], is an ear
if ( vEar[i] == EAS_NULL)
vEar[i] = ( TestTriangle( vPt, vPol, vPrev, vNext, i) ? EAS_OK : EAS_NO) ;
bool bIsEar = ( vEar[i] == EAS_OK) ;
bool bSave = bIsEar ;
double dMinAng = -1 ;
if ( bIsEar) {
// Square Length & MinAngle
double dSqLen = SqDist(vPt[vPol[vPrev[i]]], vPt[vPol[vNext[i]]]) ;
dMinAng = CalcTriangleMinAngle( vPt[vPol[vPrev[i]]], vPt[vPol[i]], vPt[vPol[vNext[i]]]) ;
// Start point for trials
int j = i ;
int k = i ;
// Try with 29 next
int nLimN = min( 29, n / 2) ;
for ( int h = 0 ; h < nLimN ; ++ h) {
j = vNext[j] ;
double dMySqLen = SqDist( vPt[vPol[vPrev[j]]], vPt[vPol[vNext[j]]]) ;
if ( dMySqLen < 2 * dSqLen) {
double dMyMinAng = CalcTriangleMinAngle( vPt[vPol[vPrev[j]]], vPt[vPol[j]], vPt[vPol[vNext[j]]]) ;
if ( dMyMinAng > 1.01 * dMinAng) {
if ( vEar[j] == EAS_NULL)
vEar[j] = ( TestTriangle( vPt, vPol, vPrev, vNext, j) ? EAS_OK : EAS_NO) ;
if ( vEar[j] == EAS_OK) {
dSqLen = dMySqLen ;
dMinAng = dMyMinAng ;
i = j ;
}
}
}
}
// Try with 29 prev
int nLimP = min( 29, n / 2) ;
for ( int h = 0 ; h < nLimP ; ++ h) {
k = vPrev[k] ;
double dMySqLen = SqDist( vPt[vPol[vPrev[k]]], vPt[vPol[vNext[k]]]) ;
if ( dMySqLen < 2 * dSqLen) {
double dMyMinAng = CalcTriangleMinAngle( vPt[vPol[vPrev[k]]], vPt[vPol[k]], vPt[vPol[vNext[k]]]) ;
if ( dMyMinAng > 1.01 * dMinAng) {
if ( vEar[k] == EAS_NULL)
vEar[k] = ( TestTriangle( vPt, vPol, vPrev, vNext, k) ? EAS_OK : EAS_NO) ;
if ( vEar[k] == EAS_OK) {
dSqLen = dMySqLen ;
dMinAng = dMyMinAng ;
i = k ;
}
}
}
}
}
// Verify if triangle is null (accept to discard)
else {
if ( Collinear( vPt[vPol[vPrev[i]]], vPt[vPol[i]], vPt[vPol[vNext[i]]]) &&
! SameDirection( vPt[vPol[vPrev[i]]], vPt[vPol[i]], vPt[vPol[vNext[i]]]))
bIsEar = true ;
}
// If current vertex v[i] is an ear, delete it and visit the previous vertex
if ( bIsEar) {
// Triangle (v[prev[i]], v[i], v[next[i]]) is an ear to save
if ( bSave) {
vTr.push_back( vPol[vPrev[i]]) ;
vTr.push_back( vPol[i]) ;
vTr.push_back( vPol[vNext[i]]) ;
dMinMinAng = min( dMinMinAng, dMinAng) ;
}
// Reset earity of diagonal endpoints
vEar[vPrev[i]] = EAS_NULL ;
vEar[vNext[i]] = EAS_NULL ;
// Delete vertex v[i] by redirecting next and previous links
// of neighboring verts past it. Decrement vertex count
vNext[vPrev[i]] = vNext[i] ;
vPrev[vNext[i]] = vPrev[i] ;
-- n ;
// Visit the previous vertex next
i = vPrev[i] ;
// Reset Count
nCount = n ;
}
else {
// Current vertex is not an ear; visit the next vertex
i = vNext[i] ;
}
}
return true ;
}
//----------------------------------------------------------------------------
bool
Triangulate::PrepareGrid( const PNTVECTOR& vPt, const INTVECTOR& vPol,
const INTVECTOR& vPrev, const INTVECTOR& vNext)
{
// points number
int n = int( vPol.size()) ;
// overall box
BBox3d b3All ;
for ( int j = 0 ; j < n ; ++ j)
b3All.Add( vPt[vPol[j]]) ;
// grid cell dimension
double dCellDim ;
b3All.GetRadius( dCellDim) ;
dCellDim *= 2 / sqrt( n) ;
// grid
if ( ! m_VertGrid.Init( 2 * n, dCellDim))
return false ;
for ( int j = 0 ; j < n ; ++ j) {
// insert only reflex vertex
if ( ! TriangleIsCCW( vPt[vPol[vPrev[j]]], vPt[vPol[j]], vPt[vPol[vNext[j]]])) {
if ( ! m_VertGrid.InsertPoint( vPt[vPol[j]], j))
return false ;
}
}
m_vVert.reserve( n / 5) ;
return true ;
}
//----------------------------------------------------------------------------
bool
Triangulate::TestTriangle( const PNTVECTOR& vPt, const INTVECTOR& vPol,
const INTVECTOR& vPrev, INTVECTOR& vNext, int i)
{
// Test if current vertex, v[i], is an ear
bool bIsEar = true ;
// An ear must be convex (here counterclockwise)
if ( TriangleIsCCW( vPt[vPol[vPrev[i]]], vPt[vPol[i]], vPt[vPol[vNext[i]]]) &&
! Collinear( vPt[vPol[vPrev[i]]], vPt[vPol[i]], vPt[vPol[vNext[i]]])) {
// Loop over all vertices not part of the tentative ear
BBox3d b3Tria ;
b3Tria.Add( vPt[vPol[vPrev[i]]]) ;
b3Tria.Add( vPt[vPol[i]]) ;
b3Tria.Add( vPt[vPol[vNext[i]]]) ;
// espando per compensare piccole approssimazioni
b3Tria.Expand( 10 * EPS_SMALL) ;
m_VertGrid.Find( b3Tria, m_vVert) ;
for ( int j = 0 ; j < int( m_vVert.size()) ; ++ j) {
int k = m_vVert[j] ;
if ( k == i || k == vPrev[i] || k == vNext[i] )
continue ;
// If vertex k is on a triangle vertex
if ( AreSamePoint( vPt[vPol[k]], vPt[vPol[vPrev[i]]]) ||
AreSamePoint( vPt[vPol[k]], vPt[vPol[i]]) ||
AreSamePoint( vPt[vPol[k]], vPt[vPol[vNext[i]]])) {
// Test previous segment with triangle diagonal
if ( TestIntersection( vPt[vPol[vPrev[k]]], vPt[vPol[k]], vPt[vPol[vPrev[i]]], vPt[vPol[vNext[i]]])) {
bIsEar = false ;
break ;
}
// Test next segment with triangle diagonal
if ( TestIntersection( vPt[vPol[k]], vPt[vPol[vNext[k]]], vPt[vPol[vPrev[i]]], vPt[vPol[vNext[i]]])) {
bIsEar = false ;
break ;
}
}
// If vertex k is inside the ear triangle, then this is not an ear
else if ( TestPointInOrOnTriangle( vPt[vPol[k]], vPt[vPol[vPrev[i]]], vPt[vPol[i]], vPt[vPol[vNext[i]]])) {
bIsEar = false ;
break ;
}
}
}
else {
// The ear triangle is clockwise so v[i] is not an ear
bIsEar = false ;
}
return bIsEar ;
}
//----------------------------------------------------------------------------
// Ratio between Circumcircle Radius and Incircle Diameter (minimum 1 for equilateral)
double
Triangulate::CalcTriangleAspectRatio( const Point3d& ptPa, const Point3d& ptPb, const Point3d& ptPc)
{
double dSqDistA = SquareDist( ptPa, ptPb) ;
double dSqDistB = SquareDist( ptPb, ptPc) ;
double dSqDistC = SquareDist( ptPc, ptPa) ;
double dTwoArea = abs( TwoArea( ptPa, ptPb, ptPc)) ;
if ( dTwoArea < SQ_EPS_SMALL)
return INFINITO ;
else
return ( max( { dSqDistA, dSqDistB, dSqDistC}) / dTwoArea) ;
}
//----------------------------------------------------------------------------
double
Triangulate::CalcTriangleMinAngle( const Point3d& ptPa, const Point3d& ptPb, const Point3d& ptPc)
{
double dSqLenA = SquareDist( ptPa, ptPb) ;
double dSqLenB = SquareDist( ptPb, ptPc) ;
double dSqLenC = SquareDist( ptPc, ptPa) ;
if ( dSqLenA > dSqLenB)
swap( dSqLenA, dSqLenB) ;
if ( dSqLenA > dSqLenC)
swap( dSqLenA, dSqLenC) ;
double dLenB = sqrt( dSqLenB) ;
double dLenC = sqrt( dSqLenC) ;
double dCosA = ( dSqLenB + dSqLenC - dSqLenA) / ( 2 * dLenB * dLenC) ;
dCosA = Clamp( dCosA, -1., +1.) ;
double dAngA = acos( dCosA) * RADTODEG ;
return dAngA ;
}
//----------------------------------------------------------------------------
double
Triangulate::SquareDist( const Point3d& ptA, const Point3d& ptB)
{
switch ( m_nPlane) {
default : // PL_XY
return (( ptB.x - ptA.x) * ( ptB.x - ptA.x) + ( ptB.y - ptA.y) * ( ptB.y - ptA.y)) ;
case PL_YZ :
return (( ptB.y - ptA.y) * ( ptB.y - ptA.y) + ( ptB.z - ptA.z) * ( ptB.z - ptA.z)) ;
case PL_ZX :
return (( ptB.z - ptA.z) * ( ptB.z - ptA.z) + ( ptB.x - ptA.x) * ( ptB.x - ptA.x)) ;
}
}
//----------------------------------------------------------------------------
// Vector product is double of the area (positive if CCW)
double
Triangulate::TwoArea( const Point3d& ptA, const Point3d& ptB, const Point3d& ptC)
{
switch ( m_nPlane) {
default : // PL_XY
return ( ptB.x - ptA.x) * ( ptC.y - ptB.y) - ( ptB.y - ptA.y) * ( ptC.x - ptB.x) ;
case PL_YZ :
return ( ptB.y - ptA.y) * ( ptC.z - ptB.z) - ( ptB.z - ptA.z) * ( ptC.y - ptB.y) ;
case PL_ZX :
return ( ptB.z - ptA.z) * ( ptC.x - ptB.x) - ( ptB.x - ptA.x) * ( ptC.z - ptB.z) ;
}
}
//----------------------------------------------------------------------------
// Distance less than tolerance
bool
Triangulate::AreSamePoint( const Point3d& ptA, const Point3d& ptB, double dToler)
{
return ( SquareDist( ptA, ptB) < dToler * dToler) ;
}
//----------------------------------------------------------------------------
// Same Direction <--> Scalar Product Positive
bool
Triangulate::SameDirection( const Point3d& ptA, const Point3d& ptB, const Point3d& ptC)
{
switch ( m_nPlane) {
default : // PL_XY
return (( ptB.x - ptA.x) * ( ptC.x - ptB.x) + ( ptB.y - ptA.y) * ( ptC.y - ptB.y)) > 0 ;
case PL_YZ :
return (( ptB.y - ptA.y) * ( ptC.y - ptB.y) + ( ptB.z - ptA.z) * ( ptC.z - ptB.z)) > 0 ;
case PL_ZX :
return (( ptB.z - ptA.z) * ( ptC.z - ptB.z) + ( ptB.x - ptA.x) * ( ptC.x - ptB.x)) > 0 ;
}
}
//----------------------------------------------------------------------------
// Collinear <--> Null Area
bool
Triangulate::Collinear( const Point3d& ptA, const Point3d& ptB, const Point3d& ptC, double dToler)
{
double dSqDistA = SquareDist( ptA, ptB) ;
double dSqDistB = SquareDist( ptB, ptC) ;
double dSqDistC = SquareDist( ptC, ptA) ;
double dTwoArea = abs( TwoArea( ptA, ptB, ptC)) ;
return ( dTwoArea * dTwoArea < dToler * dToler * max( { dSqDistA, dSqDistB, dSqDistC})) ;
}
//----------------------------------------------------------------------------
// Positive Area means A -> B -> C counterclockwise
bool
Triangulate::TriangleIsCCW( const Point3d& ptA, const Point3d& ptB, const Point3d& ptC, double dToler)
{
return ( TwoArea( ptA, ptB, ptC) > dToler * dToler) ;
}
//----------------------------------------------------------------------------
// Proper intersection between line segments
bool
Triangulate::TestIntersection( const Point3d& ptA1, const Point3d& ptA2,
const Point3d& ptB1, const Point3d& ptB2)
{
// Collinearity is not considered intersection
if ( Collinear( ptA1, ptA2, ptB1) ||
Collinear( ptA1, ptA2, ptB2) ||
Collinear( ptB1, ptB2, ptA1) ||
Collinear( ptB1, ptB2, ptA2))
return false ;
// To intersect, the endpoints of a line segment must be on opposite side of the other line segment
return ( TriangleIsCCW( ptA1, ptA2, ptB1) != TriangleIsCCW( ptA1, ptA2, ptB2) &&
TriangleIsCCW( ptB1, ptB2, ptA1) != TriangleIsCCW( ptB1, ptB2, ptA2)) ;
}
//----------------------------------------------------------------------------
// Test if point p is contained in triangle (a, b, c)
bool
Triangulate::TestPointInTriangle( const Point3d& ptP, const Point3d& ptA, const Point3d& ptB, const Point3d& ptC)
{
// If P is on a vertex is considered outside
if ( AreSamePoint( ptP, ptA))
return false ;
if ( AreSamePoint( ptP, ptB))
return false ;
if ( AreSamePoint( ptP, ptC))
return false ;
// If P is on the right of at least one edge is outside
if ( TriangleIsCCW( ptA, ptP, ptB))
return false ;
if ( TriangleIsCCW( ptB, ptP, ptC))
return false ;
if ( TriangleIsCCW( ptC, ptP, ptA))
return false ;
// P is in triangle
return true ;
}
//----------------------------------------------------------------------------
// test if point p is inside or on the border of triangle (a, b, c)
bool
Triangulate::TestPointInOrOnTriangle( const Point3d& ptP, const Point3d& ptA, const Point3d& ptB, const Point3d& ptC)
{
// If P is on a vertex is considered inside
if ( AreSamePoint( ptP, ptA))
return true ;
if ( AreSamePoint( ptP, ptB))
return true ;
if ( AreSamePoint( ptP, ptC))
return true ;
// If P is on the right of at least one edge is outside
if ( TriangleIsCCW( ptA, ptP, ptB, EPS_SMALL))
return false ;
if ( TriangleIsCCW( ptB, ptP, ptC, EPS_SMALL))
return false ;
if ( TriangleIsCCW( ptC, ptP, ptA, EPS_SMALL))
return false ;
return true ;
}
//----------------------------------------------------------------------------
bool
Triangulate::SortInternalLoops( const POLYLINEVECTOR& vPL, INTVECTOR& vOrd)
{
// riempio vettore con indice e massimo specifico per ogni loop interno
typedef pair<int,double> INDMAX ; // coppia indice, massimo
typedef vector<INDMAX> INDMAXVECTOR ; // vettore di coppie indice, massimo
INDMAXVECTOR vMax ;
vMax.reserve( vPL.size()) ;
vMax.emplace_back( 0, 0.0) ;
for ( int i = 1 ; i < int( vPL.size()) ; ++ i) {
BBox3d b3Loc ;
vPL[i].GetLocalBBox( b3Loc) ;
switch ( m_nPlane) {
default : // PL_XY
vMax.emplace_back( i, b3Loc.GetMax().x) ;
break ;
case PL_YZ :
vMax.emplace_back( i, b3Loc.GetMax().y) ;
break ;
case PL_ZX :
vMax.emplace_back( i, b3Loc.GetMax().z) ;
break ;
}
}
// ordino vettore in senso decrescente rispetto al massimo
sort( vMax.begin() + 1, vMax.end(),
[]( const INDMAX& a, const INDMAX& b) { return ( a.second > b.second) ; }) ;
// copio indice nel vettore di ordine
vOrd.reserve( vPL.size()) ;
for ( int i = 0 ; i < int( vPL.size()) ; ++ i)
vOrd.push_back( vMax[i].first) ;
return true ;
}
//----------------------------------------------------------------------------
bool
Triangulate::GetPntVectorFromPolyline( const PolyLine& PL, bool bXmaxStart, PNTVECTOR& vPi)
{
// copio i punti nel vettore (tranne il primo che coincide con l'ultimo)
Point3d ptP ;
if ( ! PL.GetFirstPoint( ptP))
return false ;
while ( PL.GetNextPoint( ptP)) {
vPi.push_back( ptP) ;
}
// se non richiesto riordino per Xmax o Ymax o Zmax, ho finito
if ( ! bXmaxStart)
return true ;
// determino l'indice del punto con Xmax
int nI = - 1 ;
double dXmax = - INFINITO ;
for ( int i = 0 ; i < int( vPi.size()) ; ++ i) {
switch ( m_nPlane) {
default : // PL_XY
if ( vPi[i].x > dXmax) {
dXmax = vPi[i].x ;
nI = i ;
}
break ;
case PL_YZ :
if ( vPi[i].y > dXmax) {
dXmax = vPi[i].y ;
nI = i ;
}
break ;
case PL_ZX :
if ( vPi[i].z > dXmax) {
dXmax = vPi[i].z ;
nI = i ;
}
break ;
}
}
if ( nI == - 1)
return false ;
// riordino il vettore, per avere il punto con Xmax in prima posizione
return ChangeStartPntVector( nI, vPi) ;
}
//----------------------------------------------------------------------------
bool
Triangulate::GetOuterPntToJoin( const PNTVECTOR& vPt, const Point3d& ptP, int& nI)
{
// cerco prima intersezione del raggio dal punto con direzione X+ al contorno esterno
nI = - 1 ;
double dMinDist = INFINITO ;
Point3d ptInt ;
int nNumPt = int( vPt.size()) ;
for ( int i = 0 ; i < nNumPt ; ++ i) {
// indice punto precedente
int h = ( i - 1 >= 0 ) ? i - 1 : nNumPt - 1 ;
// indice punto successivo
int j = ( i + 1 < nNumPt) ? i + 1 : 0 ;
// mi metto nel piano principale
switch ( m_nPlane) {
default : // PL_XY
// se i punti coincidono
if ( abs( vPt[i].x - ptP.x) < EPS_SMALL && abs( vPt[i].y - ptP.y) < EPS_SMALL) {
dMinDist = 0 ;
nI = i ;
ptInt = vPt[i] ;
}
// se punto esattamente sul raggio e raggio interno al settore
else if ( vPt[i].x > ptP.x && abs( vPt[i].y - ptP.y) < EPS_SMALL &&
PointInSector( ptP, vPt[h], vPt[i], vPt[j])) {
// se distanza minore della minima, nuovo minimo
if ( ( vPt[i].x - ptP.x) < dMinDist) {
dMinDist = vPt[i].x - ptP.x ;
nI = i ;
ptInt = vPt[i] ;
}
}
// se segmento al punto successivo che attraversa il raggio e raggio interno ovvero a sinistra ( segmento crescente in Y)
else if ( vPt[i].y < ptP.y && vPt[j].y > ptP.y) {
// calcolo l'ascissa di intersezione
double dCoeff = ( ptP.y - vPt[i].y) / ( vPt[j].y - vPt[i].y) ;
double dX = vPt[i].x + ( vPt[j].x - vPt[i].x) * dCoeff ;
// se sta sul raggio e distanza minore della minima
if ( dX > ptP.x - EPS_SMALL && ( dX - ptP.x) < dMinDist) {
dMinDist = dX - ptP.x ;
nI = ( vPt[i].x >= vPt[j].x) ? i : j ;
double dZ = vPt[i].z + ( vPt[j].z - vPt[i].z) * dCoeff ;
ptInt.Set( dX, ptP.y, dZ) ;
}
}
break ;
case PL_YZ :
// se i punti coincidono
if ( abs( vPt[i].y - ptP.y) < EPS_SMALL && abs( vPt[i].z - ptP.z) < EPS_SMALL) {
dMinDist = 0 ;
nI = i ;
ptInt = vPt[i] ;
}
// se punto esattamente sul raggio e raggio interno al settore
else if ( vPt[i].y > ptP.y && abs( vPt[i].z - ptP.z) < EPS_SMALL &&
PointInSector( ptP, vPt[h], vPt[i], vPt[j])) {
// se distanza minore della minima, nuovo minimo
if ( ( vPt[i].y - ptP.y) < dMinDist) {
dMinDist = vPt[i].y - ptP.y ;
nI = i ;
ptInt = vPt[i] ;
}
}
// se segmento al punto successivo che attraversa il raggio e raggio interno ovvero a sinistra ( segmento crescente in Y)
else if ( vPt[i].z < ptP.z && vPt[j].z > ptP.z) {
// calcolo l'ascissa di intersezione
double dCoeff = ( ptP.z - vPt[i].z) / ( vPt[j].z - vPt[i].z) ;
double dY = vPt[i].y + ( vPt[j].y - vPt[i].y) * dCoeff ;
// se sta sul raggio e distanza minore della minima
if ( dY > ptP.y - EPS_SMALL && ( dY - ptP.y) < dMinDist) {
dMinDist = dY - ptP.y ;
nI = ( vPt[i].y >= vPt[j].y) ? i : j ;
double dX = vPt[i].x + ( vPt[j].x - vPt[i].x) * dCoeff ;
ptInt.Set( dX, dY, ptP.z) ;
}
}
break ;
case PL_ZX :
// se i punti coincidono
if ( abs( vPt[i].z - ptP.z) < EPS_SMALL && abs( vPt[i].x - ptP.x) < EPS_SMALL) {
dMinDist = 0 ;
nI = i ;
ptInt = vPt[i] ;
}
// se punto esattamente sul raggio e raggio interno al settore
else if ( vPt[i].z > ptP.z && abs( vPt[i].x - ptP.x) < EPS_SMALL &&
PointInSector( ptP, vPt[h], vPt[i], vPt[j])) {
// se distanza minore della minima, nuovo minimo
if ( ( vPt[i].z - ptP.z) < dMinDist) {
dMinDist = vPt[i].z - ptP.z ;
nI = i ;
ptInt = vPt[i] ;
}
}
// se segmento al punto successivo che attraversa il raggio e raggio interno ovvero a sinistra ( segmento crescente in Y)
else if ( vPt[i].x < ptP.x && vPt[j].x > ptP.x) {
// calcolo l'ascissa di intersezione
double dCoeff = ( ptP.x - vPt[i].x) / ( vPt[j].x - vPt[i].x) ;
double dZ = vPt[i].z + ( vPt[j].z - vPt[i].z) * dCoeff ;
// se sta sul raggio e distanza minore della minima
if ( dZ > ptP.z - EPS_SMALL && ( dZ - ptP.z) < dMinDist) {
dMinDist = dZ - ptP.z ;
nI = ( vPt[i].z >= vPt[j].z) ? i : j ;
double dY = vPt[i].y + ( vPt[j].y - vPt[i].y) * dCoeff ;
ptInt.Set( ptP.x, dY, dZ) ;
}
}
break ;
}
}
// non ho trovato alcunchè, errore
if ( nI == - 1)
return false ;
// se ho trovato un punto esatto del contorno, non devo fare altri controlli
if ( AreSamePointApprox( ptInt, vPt[nI]))
return true ;
// devo ora verificare che il segmento che unisce i punti non intersechi altri lati del contorno esterno
// altrimenti tengo il punto con raggio più vicino a X_AX o Y_AX o Z_AX secondo m_nPlane
int nJ = nI ;
Point3d ptPa = ptP ;
Point3d ptPb = vPt[nI] ;
Point3d ptPc = ptInt ;
bool bSwap = false ;
switch ( m_nPlane) {
default : /* PL_XY */ bSwap = ( ptPb.y > ptPc.y) ; break ;
case PL_YZ : bSwap = ( ptPb.z > ptPc.z) ; break ;
case PL_ZX : bSwap = ( ptPb.x > ptPc.x) ; break ;
}
if ( bSwap)
swap( ptPb, ptPc) ;
double dMinTan = INFINITO ;
double dMinSqDist = SQ_INFINITO ;
for ( int i = 0 ; i < nNumPt ; ++ i) {
// salto il punto già trovato
if ( i == nJ)
continue ;
// verifico se sta nel triangolo
if ( TestPointInTriangle( vPt[i], ptPa, ptPb, ptPc)) {
double dTan = INFINITO ;
switch ( m_nPlane) {
default : /* PL_XY */ dTan = abs(vPt[i].y - ptP.y) / (vPt[i].x - ptP.x) ; break ;
case PL_YZ : dTan = abs(vPt[i].z - ptP.z) / (vPt[i].y - ptP.y) ; break ;
case PL_ZX : dTan = abs(vPt[i].x - ptP.x) / (vPt[i].z - ptP.z) ; break ;
}
if ( dTan < dMinTan + EPS_ZERO) {
// indice punto precedente
int h = ( i - 1 >= 0) ? i - 1 : nNumPt - 1 ;
// indice punto successivo
int j = ( i + 1 < nNumPt) ? i + 1 : 0 ;
// verifico che il raggio che unisce i punti stia nel settore
if ( PointInSector( ptP, vPt[h], vPt[i], vPt[j])) {
double dSqDist = SqDist( vPt[i], ptP) ;
if ( dTan < dMinTan - EPS_ZERO || dSqDist < dMinSqDist) {
dMinTan = dTan ;
dMinSqDist = dSqDist ;
nI = i ;
}
}
}
}
}
return true ;
}
//----------------------------------------------------------------------------
bool
Triangulate::PointInSector( const Point3d& ptTest, const Point3d& ptPrev, const Point3d& ptCorn, const Point3d& ptNext)
{
// la parte valida del settore è a sinistra dei segmenti ptPrev --> ptCorn --> ptNext
// se corner convesso
if ( TriangleIsCCW( ptPrev, ptCorn, ptNext, 0))
return ( TriangleIsCCW( ptPrev, ptCorn, ptTest) &&
TriangleIsCCW( ptTest, ptCorn, ptNext)) ;
// altrimenti corner concavo ( reflex)
else
return ! ( TriangleIsCCW( ptTest, ptCorn, ptPrev, 0) &&
TriangleIsCCW( ptNext, ptCorn, ptTest, 0)) ;
}
//----------------------------------------------------------------------------
bool
ChangeStartPntVector( int nNewStart, PNTVECTOR& vPi)
{
// se il nuovo inizio coincide col vecchio, non devo fare alcunchè
if ( nNewStart == 0)
return true ;
// se il nuovo indice è oltre la dimensione del vettore, errore
if ( nNewStart >= int( vPi.size()))
return false ;
// ciclo di aggiustamento
rotate( vPi.begin(), vPi.begin() + nNewStart, vPi.end()) ;
return true ;
}
//----------------------------------------------------------------------------
// Funzioni per triangolare usando la libreria FIST
//----------------------------------------------------------------------------
//----------------------------------------------------------------------------
int
CalcTriangleType( const PNTVECTOR& vPt, const INTVECTOR& vTr, int nTrg)
{
Vector3d vtV1 = vPt[vTr[3*nTrg + 1]] - vPt[vTr[3*nTrg]] ;
Vector3d vtV2 = vPt[vTr[3*nTrg + 2]] - vPt[vTr[3*nTrg + 1]] ;
Vector3d vtN = vtV1 ^ vtV2 ;
double dSqN = vtN.SqLen() ;
double dSqLen1 = vtV1.SqLen() ;
double dSqLen2 = vtV2.SqLen() ;
double dSqLen3 = ( vtV1 + vtV2).SqLen() ;
// un triangolo è invalido ( needle o cap) se viene scartato dai controlli di AddTriangle
if ( dSqN < SQ_EPS_ZERO || dSqN < SQ_EPS_TRIA_H * max( {dSqLen1, dSqLen2, dSqLen3})) {
if ( dSqLen1 < SQ_EPS_ZERO || dSqLen2 < SQ_EPS_ZERO || dSqLen3 < SQ_EPS_ZERO)
return TRG_NEEDLE ;
else
return TRG_CAP ;
}
else
return TRG_STD ;
}
//----------------------------------------------------------------------------
bool
FlipTrg( INTVECTOR& vTr, int nTA, int nTB, int nEA, int nEB)
{
vTr[3*nTA + nEA] = vTr[3*nTB + ( nEB + 2) % 3] ;
vTr[3*nTB + nEB] = vTr[3*nTA + ( nEA + 2) % 3] ;
return true ;
}
//----------------------------------------------------------------------------
bool
TestAdjacentOnEdge( INTVECTOR& vTr, int nTA, int nEA, int nTTest1, int nTTest2, int& nTB, int& nEB)
{
// individuo quale triangolo tra nTTest1 e nTTest2 è quello adiacente a nTA lungo il lato nEA
nTB = -1 ;
nEB = -1 ;
// recupero i vertici di nEA
int nV1 = vTr[3*nTA + nEA] ;
int nV2 = vTr[3*nTA + ( nEA + 1) % 3] ;
// verifico se è TTest1 quello adiacente
if ( vTr[3*nTTest1] == nV2 && vTr[3*nTTest1+1] == nV1) {
nEB = 0 ;
nTB = nTTest1 ;
}
else if ( vTr[3*nTTest1+1] == nV2 && vTr[3*nTTest1+2] == nV1) {
nEB = 1 ;
nTB = nTTest1 ;
}
else if ( vTr[3*nTTest1] == nV1 && vTr[3*nTTest1+2] == nV2) {
nEB = 2 ;
nTB = nTTest1 ;
}
// verifico se è TTest2 quello adiacente
else if ( vTr[3*nTTest2] == nV2 && vTr[3*nTTest2+1] == nV1) {
nEB = 0 ;
nTB = nTTest2 ;
}
else if ( vTr[3*nTTest2+1] == nV2 && vTr[3*nTTest2+2] == nV1) {
nEB = 1 ;
nTB = nTTest2 ;
}
else if ( vTr[3*nTTest2] == nV1 && vTr[3*nTTest2+2] == nV2) {
nEB = 2 ;
nTB = nTTest2 ;
}
return true ;
}
//----------------------------------------------------------------------------
bool
FindAdjacentOnLongerEdge( PNTVECTOR& vPt, INTVECTOR& vTr, int nTA, int& nEA, int& nTB, int& nEB)
{
// recupero il lato più lungo del triangolo nTA
nEA = -1 ;
double dMaxLen = -1 ;
for ( int j = 0 ; j < 3 ; j++) {
double dCurrLen = SqDist( vPt[vTr[3*nTA + j]], vPt[vTr[3*nTA + ( j+1)%3]]) ;
if ( dCurrLen > dMaxLen) {
dMaxLen = dCurrLen ;
nEA = j ;
}
}
// cerco il triangolo nTB adiacente a nTA lungo nEA
nTB = -1 ;
nEB = -1 ;
int nV1 = vTr[3*nTA + nEA] ;
int nV2 = vTr[3*nTA + ( nEA + 1) % 3] ;
int nTria = vTr.size() / 3 ;
for ( int j = 0 ; j < nTria ; j++) {
if ( j == nTA)
continue ;
if ( vTr[3*j] == nV2 && vTr[3*j+1] == nV1) {
nEB = 0 ;
nTB = j ;
break ;
}
else if ( vTr[3*j+1] == nV2 && vTr[3*j+2] == nV1) {
nEB = 1 ;
nTB = j ;
break ;
}
else if ( vTr[3*j] == nV1 && vTr[3*j+2] == nV2) {
nEB = 2 ;
nTB = j ;
break ;
}
}
return true ;
}
//----------------------------------------------------------------------------
bool
RemoveFistInvalidTrg( PNTVECTOR& vPt, INTVECTOR& vTr)
{
// scorro tutti i triangoli alla ricerca di triangoli validi per FIST ma invalidi per le nostre tolleranze.
// I triangoli invalidi possono essere di due tipologie: cap o needle.
// I triangoli cap se eliminati danno origine a T-junctions, quindi devono essere gestiti opportunamente con dei flip.
// I triangoli needle se eliminati non sono problematici, ma i loro vertici coincidenti vanno gestiti opportunamente nel
// calcolo delle adiacenze dei triangoli cap.
// TO DO da capire e gestire casi in cui flip lascia triangoli invalidi
int nTria = int( vTr.size()) / 3 ;
INTVECTOR vCapTria ;
INTVECTOR vNeedleTria ;
BOOLVECTOR vbIsValidTria( nTria, true) ;
bool bRemovedTrg = false ;
for ( int i = 0 ; i < nTria ; i ++) {
int nTrgType = CalcTriangleType( vPt, vTr, i) ;
if ( nTrgType == TRG_CAP) {
vbIsValidTria[i] = false ;
vCapTria.emplace_back( i) ;
}
else if ( nTrgType == TRG_NEEDLE) {
vNeedleTria.emplace_back( i) ;
}
}
// se non ci sono triangoli cap non c'è bisogno di modifiche
if ( vCapTria.empty())
return true ;
// 1) elimino i triangoli di tipo needle gestendo i vertici coincidenti
if ( ! vNeedleTria.empty()) {
bRemovedTrg = true ;
for ( int i = 0 ; i < int( vNeedleTria.size()) ; i++) {
int nT = vNeedleTria[i] ;
if ( vTr[3*nT] != vTr[3*nT+1] && vTr[3*nT] != vTr[3*nT+2] && vTr[3*nT+1] != vTr[3*nT+2]) {
// individuo i vertici coincidenti
int nOldId = -1 ;
int nNewId = -1 ;
Vector3d vtV1 = vPt[vTr[3*nT+1]] - vPt[vTr[3*nT]] ;
if ( vtV1.SqLen() < SQ_EPS_ZERO) {
// vertici coincidenti 0 e 1
nOldId = vTr[3*nT] ;
nNewId = vTr[3*nT+1] ;
}
else {
Vector3d vtV2 = vPt[vTr[3*nT+2]] - vPt[vTr[3*nT+1]] ;
if ( vtV2.SqLen() < SQ_EPS_ZERO) {
// vertici coincidenti 1 e 2
nOldId = vTr[3*nT+1] ;
nNewId = vTr[3*nT+2] ;
}
else {
// vertici coincidenti 0 e 2
nOldId = vTr[3*nT] ;
nNewId = vTr[3*nT+2] ;
}
}
// aggiorno il vettore dei triangoli unificando i vertici coincidenti
for ( int k = 0 ; k < int( vTr.size()) ; k ++) {
if ( vTr[k] == nOldId)
vTr[k] = nNewId ;
}
}
// annullo il triangolo
vTr[3*nT] = -1 ;
vTr[3*nT+1] = -1 ;
vTr[3*nT+2] = -1 ;
}
}
// 2) sistemo i triangoli di tipo cap
for ( int i = 0 ; i < int( vCapTria.size()) ; i++) {
int nTA = vCapTria[i] ;
// verifico se il triangolo è già stato resto valido ( da una catena)
if ( vbIsValidTria[nTA])
continue ;
// trovo l'adiacenza sul suo lato più lungo
int nEA = -1 ;
int nEB = -1 ;
int nTB = -1 ;
FindAdjacentOnLongerEdge( vPt, vTr, nTA, nEA, nTB, nEB) ;
if ( nTB == -1) {
// se non ha un triangolo adiacente posso annullarlo
vTr[3*nTA] = -1 ;
vTr[3*nTA+1] = -1 ;
vTr[3*nTA+2] = -1 ;
vbIsValidTria[nTA] = true ;
bRemovedTrg = true ;
continue ;
}
else if ( vbIsValidTria[nTB]) {
// se è adiacente è valido, eseguo il flip
FlipTrg( vTr, nTA, nTB, nEA, nEB) ;
vbIsValidTria[nTA] = true ;
}
else {
// se adiacente è invalido, individuo una catena di triangoli invalidi da risolvere non appena si indentifica
// adiacenza con triangolo valido
INTVECTOR vChain, vChainEdges ;
vChain.emplace_back( nTA) ;
vChainEdges.emplace_back( nEA) ;
int nTCurr = nTB ;
int nEOther ;
while ( nTCurr != -1 && ! vbIsValidTria[nTCurr]) {
// calcolo il successivo
int nTOther, nECurr ;
FindAdjacentOnLongerEdge( vPt, vTr, nTCurr, nECurr, nTOther, nEOther) ;
if ( nTOther == vChain.back()) {
// se ho trovato un'adiacenza ambigua ( ovvero due triangoli invalidi adiacenti sui loro lati più lunghi)
// flip dei due triangoli per modificare il lato più lungo e togliere adiacenza ambigua
FlipTrg( vTr, nTCurr, nTOther, nECurr, nEOther) ;
if ( vChain.size() > 1) {
vChain.pop_back() ;
vChainEdges.pop_back() ;
// individuo il nuovo adiacente all'ultimo triangolo della catena tra i due appena flippati
TestAdjacentOnEdge( vTr, vChain.back(), vChainEdges.back(), nTCurr, nTOther, nTB, nEB) ;
nTOther = nTB ;
}
}
else {
vChain.emplace_back( nTCurr) ;
vChainEdges.emplace_back( nECurr) ;
}
// aggiorno per iterazione successiva
nTCurr = nTOther ;
}
// se la catena termina su triangolo nullo, annullo tutti i triangoli della catena
if ( nTCurr == -1) {
bRemovedTrg = true ;
for ( int k = 0 ; k < int( vChain.size()) ; k++) {
if ( vbIsValidTria[vChain[k]])
continue ;
vTr[3*vChain[k]] = -1 ;
vTr[3*vChain[k] + 1] = -1 ;
vTr[3*vChain[k] + 2] = -1 ;
vbIsValidTria[vChain[k]] = true ;
}
}
// se catena termina su un triangolo valido, applico il flip a cascata a partire dall'ultimo triangolo invalido
else {
FlipTrg( vTr, vChain.back(), nTCurr, vChainEdges.back(), nEOther) ;
vbIsValidTria[vChain.back()] = true ;
int nTrgTest1 = vChain.back() ;
int nTrgTest2 = nTCurr ;
for ( int j = int( vChain.size()-2) ; j >= 0 ; j--) {
// triangolo corrente
int nTA = vChain[j] ;
if ( vbIsValidTria[nTA])
continue ;
int nEA = vChainEdges[j] ;
// devo trovare il nuovo adiacente dopo il flip dei successivi nella catena
TestAdjacentOnEdge( vTr, nTA, nEA, nTrgTest1, nTrgTest2, nTB, nEB) ;
// flip per rendere valido il triangolo corrente
FlipTrg( vTr, nTA, nTB, nEA, nEB) ;
vbIsValidTria[nTA] = true ;
// aggiorno i trg di test per lo step successivo
nTrgTest1 = nTA ;
nTrgTest2 = nTB ;
}
}
}
}
// se sono stati annullati dei triangoli, li rimuovo dal vettore
if ( bRemovedTrg) {
INTVECTOR vTrTmp ;
vTrTmp.reserve( vTr.size()) ;
for ( int i = 0 ; i < int( vTr.size()) ; i++)
if ( vTr[i] != -1)
vTrTmp.emplace_back( vTr[i]) ;
swap( vTr, vTrTmp) ;
}
return true ;
}
//----------------------------------------------------------------------------
bool
MakeByFist( const POLYLINEVECTOR& vPL, PNTVECTOR& vPt, INTVECTOR& vTr)
{
// opzioni di triangolazione
rt_options rt_opt ;
InitDefaults( &rt_opt) ;
rt_opt.ears_top = false ;
rt_opt.ears_random = false ;
rt_opt.ears_sorted = true ;
rt_opt.ears_fancy = true ;
// creo oggetto di fist per triangolazione
global_struct fist ;
InitGlobalStruct( &fist, &rt_opt, true) ;
// allocazione ottimizzata degli array di fist
int nLoops = int( vPL.size()) ;
int nVertices = 0 ;
for ( int i = 0 ; i < nLoops ; i ++)
nVertices += ( vPL[i].GetPointNbr() - 1) ;
OptimizeMemoryAllocation( &fist, nLoops, nVertices) ;
// assegno i dati a fist
for ( int i = 0 ; i < nLoops ; i ++) {
// salvo i vertici saltando il primo punto ( coincide con l'ultimo)
Point3d pt ;
if ( ! vPL[i].GetFirstPoint( pt))
return false ;
while ( vPL[i].GetNextPoint( pt))
StoreVertex( &fist.c_vertex, pt.x, pt.y, pt.z) ;
// salvo il loop
int nPoints = vPL[i].GetPointNbr() ;
AddLoopInFace( &fist, nLoops - 1, i == 0, nPoints - 1) ;
}
StoreGroupNumber( &fist.c_list, &fist.c_vertex) ;
// eseguo triangolazione
Triangulate( &fist) ;
// recupero i vertici da fist
vPt.reserve( fist.c_vertex.num_vertices) ;
for ( int i = 0 ; i < fist.c_vertex.num_vertices ; i ++)
vPt.emplace_back( fist.c_vertex.vertices[i].x, fist.c_vertex.vertices[i].y, fist.c_vertex.vertices[i].z) ;
// recupero i triangoli da fist
vTr.reserve( 3 * fist.c_vertex.num_triangles) ;
for ( int i = 0 ; i < fist.c_vertex.num_triangles ; i ++) {
vTr.emplace_back( fist.c_vertex.triangles[i].v1) ;
vTr.emplace_back( fist.c_vertex.triangles[i].v2) ;
vTr.emplace_back( fist.c_vertex.triangles[i].v3) ;
}
// rimuovo eventuali triangoli validi per fist ma invalidi per le nostre tolleranze
RemoveFistInvalidTrg( vPt, vTr) ;
// chiudo fist e dealloco la sua memoria
FIST_TerminateProgram( &fist) ;
return true ;
}
Reference in New Issue View Git Blame Copy Permalink