diff --git a/EgtGeomKernel.vcxproj b/EgtGeomKernel.vcxproj
index fc4f6b9..614a5ec 100644
--- a/EgtGeomKernel.vcxproj
+++ b/EgtGeomKernel.vcxproj
@@ -148,6 +148,7 @@ copy $(TargetPath) \EgtProg\DllD32
false
true
ProgramDatabase
+ \EgtDev\Extern\GeometricToolsEngine
Windows
@@ -239,6 +240,7 @@ copy $(TargetPath) \EgtProg\Dll32
true
false
None
+ \EgtDev\Extern\GeometricToolsEngine
Windows
diff --git a/Triangulate.cpp b/Triangulate.cpp
index d88fccd..17637cf 100644
--- a/Triangulate.cpp
+++ b/Triangulate.cpp
@@ -24,10 +24,12 @@
#include "/EgtDev/Include/EGkPolyLine.h"
#include "/EgtDev/Include/EGkPlane3d.h"
#include "/EgtDev/Include/EGkStringUtils3d.h"
+#include "Mathematics/TriangulateCDT.h"
#include
using namespace std ;
using namespace tpp ;
+using namespace gte ;
//----------------------------------------------------------------------------
enum EarStatus{ EAS_NULL = -1, EAS_NO = 0, EAS_OK = 1} ;
@@ -150,6 +152,13 @@ Triangulate::Make( const POLYLINEVECTOR& vPL, PNTVECTOR& vPt, INTVECTOR& vTr, Tr
}
return true ;
}
+ else if ( trgType == TRG_DEL_GTE) {
+ if ( ! MakeByDelaunayGTE( vPL, vPt, vTr)) {
+ LOG_ERROR( GetEGkLogger(), "Error in MakeByDelaunay ( Geometric Tools Engine)") ;
+ return false ;
+ }
+ return true ;
+ }
// ear clipping
// se non CCW inverto tutte le polilinee
@@ -233,6 +242,7 @@ Triangulate::Make( const POLYLINEVECTOR& vPL, PNTVECTOR& vPt, INTVECTOR& vTr, Tr
bool
Triangulate::MakeByDelaunay( const POLYLINEVECTOR& vPL, PNTVECTOR& vPt, INTVECTOR& vTr, bool bConforming, bool bQuality)
{
+ // sistema di riferimento locale sul piano della polyline
Plane3d plPlane ;
double dArea;
if ( ! vPL[0].IsClosedAndFlat( plPlane, dArea, 50 * EPS_SMALL))
@@ -243,41 +253,42 @@ Triangulate::MakeByDelaunay( const POLYLINEVECTOR& vPL, PNTVECTOR& vPt, INTVECTO
POLYLINEVECTOR vMyPL = vPL ;
for ( size_t t = 0 ; t < vPL.size() ; ++ t) {
vMyPL[t].ToLoc( frLoc) ;
+ if ( vMyPL[t].GetPointNbr() < 4) // per poter triangolare devo avere almeno tre punti distinti
+ return false ;
}
- // vertici e constraint per la triangolazione
+ // vertici e constraint per la triangolazione
vector vDelaunayVert ;
- INTVECTOR vDelaunayConstr ;
-
+ INTVECTOR vDelaunayConstr ;
for ( size_t i = 0 ; i < vMyPL.size() ; i++) {
Point3d ptFirst, pt ;
vMyPL[i].GetFirstPoint( ptFirst) ;
vDelaunayVert.push_back( Delaunay::Point( ptFirst.x, ptFirst.y)) ;
- int nFirst = vDelaunayVert.size() - 1 ; // indice del primo punto della polyline fra i vertici della triangolazione
+ int nFirst = ( int)vDelaunayVert.size() - 1 ; // indice del primo punto della polyline fra i vertici della triangolazione
vDelaunayConstr.push_back( nFirst) ;
while ( vMyPL[i].GetNextPoint( pt)) {
vDelaunayVert.push_back( Delaunay::Point( pt.x, pt.y)) ;
- // nei constraint gli indici vanno inseriti due volte perchè vengono letti a due a due per definire un segmento
- vDelaunayConstr.push_back( vDelaunayVert.size() - 1) ;
- vDelaunayConstr.push_back( vDelaunayVert.size() - 1) ;
+ // nei constraint gli indici vanno inseriti due volte perchè vengono letti a due a due per definire un segmento
+ vDelaunayConstr.push_back( ( int)vDelaunayVert.size() - 1) ;
+ vDelaunayConstr.push_back( ( int)vDelaunayVert.size() - 1) ;
}
- // impongo che l'ultimo vertice coincida con il primo ( curva chiusa)
+ // impongo che l'ultimo vertice coincida con il primo ( curva chiusa)
vDelaunayVert.pop_back() ;
vDelaunayConstr.erase(vDelaunayConstr.end() - 2, vDelaunayConstr.end()) ;
vDelaunayConstr.push_back( nFirst) ;
}
- // holes : sono definiti da un punto interno al buco
+ // holes : sono definiti da un punto interno al buco
std::vector vDelaunayHoles ;
for ( size_t i = 1 ; i < vMyPL.size() ; i++) {
Point3d pt ;
CurveComposite * pCrvHole = CreateBasicCurveComposite() ;
pCrvHole->FromPolyLine( vMyPL[i]) ;
pCrvHole->GetCentroid( pt) ;
- // se il centroide fosse esterno, cerco per tentativi un punto qualsiasi all'interno della curva
+ // se il centroide fosse esterno, cerco per tentativi un punto qualsiasi all'interno della curva
double dPar = 0.5 ;
while ( ! IsPointInsidePolyLine( pt, vMyPL[i])) {
double dParS, dParE ;
@@ -288,14 +299,14 @@ Triangulate::MakeByDelaunay( const POLYLINEVECTOR& vPL, PNTVECTOR& vPt, INTVECTO
pCrvHole->GetPointTang( dPar, ICurve::FROM_MINUS, pt, vtDir) ;
vtDir.Rotate( Z_AX, 0, -1) ;
pt += 2 * EPS_SMALL * vtDir ;
- // tento con un nuovo punto
+ // tento con un nuovo punto
dPar += 10 * EPS_SMALL ;
}
vDelaunayHoles.push_back( Delaunay::Point( pt.x, pt.y)) ;
}
- // parti concave
+ // parti concave
PNTVECTOR vPtConvexHull ;
vMyPL[0].GetConvexHullXY( vPtConvexHull) ;
CurveComposite* pCrv = CreateBasicCurveComposite() ;
@@ -306,13 +317,13 @@ Triangulate::MakeByDelaunay( const POLYLINEVECTOR& vPL, PNTVECTOR& vPt, INTVECTO
CurveLine * pLine = CreateBasicCurveLine() ;
pLine->Set( vPtConvexHull[i], vPtConvexHull[NextIdx]) ;
- // verifico se ho regioni da eliminare
+ // verifico se ho regioni da eliminare
IntersCurveCurve intCC( *pLine, *pCrv) ;
CRVCVECTOR ccClass ;
intCC.GetCurveClassification( 0, ccClass) ;
for ( size_t j = 0 ; j < ccClass.size() ; j ++) {
if ( ccClass[j].nClass == CRVC_OUT) {
- // cerco per tentativi un punto a caso all'interno della regione da eliminare
+ // cerco per tentativi un punto a caso all'interno della regione da eliminare
double dPar = ( ccClass[j].dParS + ccClass[j].dParE) / 2 ;
double dDist = 0 ;
Point3d pt ;
@@ -329,14 +340,14 @@ Triangulate::MakeByDelaunay( const POLYLINEVECTOR& vPL, PNTVECTOR& vPt, INTVECTO
pt += EPS_SMALL * ( vtDir) ;
vDelaunayHoles.push_back( Delaunay::Point( pt.x , pt.y)) ;
}
- // tento con nuovo punto
+ // tento con nuovo punto
dPar += 10 * EPS_SMALL ;
}
}
}
}
- // triangolazione
+ // triangolazione
Delaunay trGenerator( vDelaunayVert) ;
trGenerator.setSegmentConstraint( vDelaunayConstr) ;
trGenerator.setHolesConstraint( vDelaunayHoles) ;
@@ -345,11 +356,11 @@ Triangulate::MakeByDelaunay( const POLYLINEVECTOR& vPL, PNTVECTOR& vPt, INTVECTO
else
trGenerator.Triangulate( bQuality) ;
- // se non ho generato triangoli, errore
+ // se non ho generato triangoli, errore
if ( trGenerator.ntriangles() == 0)
return false ;
- // vertici
+ // vertici
std::vector< Delaunay::Point> vVertex = trGenerator.MyVertexTraverse() ;
for (size_t i = 0; i < vVertex.size(); i++) {
Point3d pt( vVertex[i][0], vVertex[i][1]) ;
@@ -357,13 +368,229 @@ Triangulate::MakeByDelaunay( const POLYLINEVECTOR& vPL, PNTVECTOR& vPt, INTVECTO
vPt.push_back( pt) ;
}
- // triangoli
+ // triangoli
vTr = trGenerator.MyTriangleTraverse() ;
return true ;
}
+//----------------------------------------------------------------------------
+// Delaunay Triangulation ( Geometric Tools Engine library)
+//----------------------------------------------------------------------------
+bool
+Triangulate::MakeByDelaunayGTE( const POLYLINEVECTOR& vPL, PNTVECTOR& vPt, INTVECTOR& vTr)
+{
+ // sistema di riferimento locale sul piano della polyline
+ Plane3d plPlane ;
+ double dArea;
+ if ( ! vPL[0].IsClosedAndFlat( plPlane, dArea, 50 * EPS_SMALL))
+ return false ;
+ Frame3d frLoc ;
+ frLoc.Set( plPlane.GetPoint(), plPlane.GetVersN()) ;
+
+ POLYLINEVECTOR vMyPL = vPL ;
+ for ( size_t t = 0 ; t < vPL.size() ; ++ t) {
+ vMyPL[t].ToLoc( frLoc) ;
+ if ( vMyPL[t].GetPointNbr() < 4)
+ return false ;
+ }
+
+ CurveComposite * pCrvPL = CreateBasicCurveComposite() ;
+ pCrvPL->FromPolyLine( vMyPL[0]) ;
+
+ // controllo che holes siano contenuti nel loop esterno
+ for ( size_t i = 1 ; i < vMyPL.size() ; i ++) {
+ CurveComposite * pCrvHole = CreateBasicCurveComposite() ;
+ pCrvHole->FromPolyLine( vMyPL[i]) ;
+
+ IntersCurveCurve intCC( *pCrvPL, *pCrvHole) ;
+ CRVCVECTOR ccClass ;
+ intCC.GetCurveClassification( 1, ccClass) ;
+
+ if ( ccClass.size() > 1 || ccClass[0].nClass != CRVC_IN)
+ return false ;
+ }
+
+ // Verifico non ci siano autointersezioni
+ SelfIntersCurve Inters( *pCrvPL) ;
+ if ( Inters.GetIntersCount() > 0 && Inters.GetOverlaps() == 0)
+ return false ;
+
+ // Gestione autointersezione
+ // if ( Inters.GetIntersCount() > 0 && Inters.GetOverlaps() == 0) {
+ //
+ // // vertici della triangolazione
+ // vector> vPoints ;
+ // Point3d pt ;
+ // vMyPL[0].GetFirstPoint( pt) ;
+ // vPoints.push_back( {pt.x, pt.y}) ;
+ // pt.ToGlob( frLoc) ;
+ // vPt.push_back( pt) ;
+ // while ( vMyPL[0].GetNextPoint( pt) ) {
+ // vPoints.push_back( {pt.x, pt.y}) ;
+ // pt.ToGlob( frLoc) ;
+ // vPt.push_back( pt) ;
+ // }
+ // vPt.pop_back() ;
+ // vPoints.pop_back() ;
+ //
+ // // i punti di autointersezione P sono salvati come pair: parametro di P, indice posizione di P in vPoints
+ // vector> dParams ;
+ // for ( int i = 0 ; i < Inters.GetIntersCount() ; i++) {
+ //
+ // IntCrvCrvInfo aInfo ;
+ // Inters.GetIntCrvCrvInfo( i, aInfo) ;
+ //
+ // Point3d pt ;
+ // pt = aInfo.IciA[0].ptI ;
+ // pt.ToGlob( frLoc) ;
+ //
+ // auto it = find_if( vPt.begin(), vPt.end(), [pt]( const Point3d& ptCrv){ return AreSamePointApprox( ptCrv, pt) ;} ) ;
+ // // se il punto di interezione è uno dei vertici
+ // if ( it != vPt.end()) {
+ // dParams.push_back( pair (aInfo.IciA[0].dU, it - vPt.begin())) ;
+ // dParams.push_back( pair (aInfo.IciB[0].dU, it - vPt.begin())) ;
+ // }
+ // // se non è uno dei vertici, lo aggiungo in vPt e in vPoints
+ // else {
+ // vPt.push_back( pt) ;
+ // pt.ToLoc( frLoc) ;
+ // vPoints.push_back( { pt.x, pt.y}) ;
+ // dParams.push_back( pair (aInfo.IciA[0].dU, vPoints.size() - 1)) ;
+ // dParams.push_back( pair (aInfo.IciB[0].dU, vPoints.size() - 1)) ;
+ // }
+ // }
+ // sort( dParams.begin(), dParams.end(), []( const pair& a, const pair& b){ return a.first < b.first ; }) ;
+ //
+ // // lati del poligono :
+ // INTVECTOR vPolygonIdx ;
+ // for ( int i = 0 ; i < vMyPL[0].GetPointNbr() ; i++)
+ // vPolygonIdx.push_back( i) ;
+ //
+ // // inserisco i nuovi lati determinati dai punti di intersezione
+ // for ( int i = 0 ; i < dParams.size() ; i++) {
+ // // se il parametro c'è già, non lo inserisco
+ // if ( find_if( vPolygonIdx.begin(), vPolygonIdx.end(),
+ // [dParams, i]( const int& val) { return abs( val - dParams[i].first) < 0 ; }) != vPolygonIdx.end())
+ // break ;
+ // // altrimenti lo inserisco in modo ordinato nel vettore
+ // auto it = lower_bound( vPolygonIdx.begin(), vPolygonIdx.end(), dParams[i].first,
+ // []( const int& a , const int& b) { return a < b ; }) ;
+ // vPolygonIdx.insert( it, dParams[i].second) ;
+ // }
+ //
+ // auto pPolygon = std::make_shared() ;
+ // pPolygon->polygon = vPolygonIdx ;
+ //
+ // // holes
+ // int k = vPoints.size() ;
+ // for ( size_t i = 1 ; i < vMyPL.size() ; i++) {
+ //
+ // auto pHole = std::make_shared() ;
+ // Point3d pt ;
+ // vMyPL[i].GetFirstPoint( pt) ;
+ // vPoints.push_back( { pt.x, pt.y}) ;
+ // pHole->polygon.push_back( k) ;
+ // k ++ ;
+ // while ( vMyPL[i].GetNextPoint( pt)) {
+ // vPoints.push_back( { pt.x, pt.y}) ;
+ // pHole->polygon.push_back( k) ;
+ // k ++ ;
+ // }
+ //
+ // vPoints.pop_back() ;
+ // pHole->polygon.pop_back() ;
+ // k -- ;
+ //
+ // pPolygon->child.push_back( pHole) ;
+ //
+ // }
+ //
+ // // Triangolazione
+ // PolygonTreeEx PTOutput ;
+ // TriangulateCDT triangulator ;
+ // triangulator( vPoints, pPolygon, PTOutput) ;
+ //
+ // for ( auto const& tri : PTOutput.interiorTriangles) {
+ // vTr.push_back( tri[0]) ;
+ // vTr.push_back( tri[1]) ;
+ // vTr.push_back( tri[2]) ;
+ // }
+ //
+ // return true ;
+ //
+ // }
+
+
+ // i vertici della triangolazione vPt sono i punti delle polyline nel sistema globale
+ for ( size_t i = 0 ; i < vPL.size() ; i++) {
+ Point3d pt ;
+ vPL[i].GetFirstPoint( pt) ;
+ vPt.push_back( Point3d( pt.x, pt.y, pt.z)) ;
+ while ( vPL[i].GetNextPoint( pt ) )
+ vPt.push_back( Point3d( pt.x, pt.y, pt.z)) ;
+ // ultimo punto della polyline coincide con il primo, quindi non serve
+ vPt.pop_back() ;
+ }
+
+ vector> vPoints ; // vertici per la triangolazione
+ auto pPolygon = std::make_shared() ; // poligono da triangolare
+
+ Point3d pt ;
+ int k = 0 ;
+ vMyPL[0].GetFirstPoint( pt) ;
+ vPoints.push_back( { pt.x, pt.y}) ;
+ pPolygon->polygon.push_back( k) ;
+ k ++ ;
+
+ while ( vMyPL[0].GetNextPoint( pt) ) {
+ vPoints.push_back( { pt.x, pt.y}) ;
+ pPolygon->polygon.push_back( k) ;
+ k ++ ;
+ }
+
+ // ultimo punto della polyline coincide con il primo, quindi non serve
+ vPoints.pop_back() ;
+ pPolygon->polygon.pop_back() ;
+ k -- ;
+
+ // holes
+ for ( size_t i = 1 ; i < vMyPL.size() ; i++) {
+
+ auto pHole = std::make_shared() ;
+ vMyPL[i].GetFirstPoint( pt) ;
+ vPoints.push_back( { pt.x, pt.y}) ;
+ pHole->polygon.push_back( k) ;
+ k ++ ;
+ while ( vMyPL[i].GetNextPoint( pt)) {
+ vPoints.push_back( { pt.x, pt.y}) ;
+ pHole->polygon.push_back( k) ;
+ k ++ ;
+ }
+
+ vPoints.pop_back() ;
+ pHole->polygon.pop_back() ;
+ k -- ;
+
+ pPolygon->child.push_back( pHole) ;
+
+ }
+
+ // Triangolazione
+ PolygonTreeEx PTOutput ;
+ TriangulateCDT triangulator ;
+ triangulator( vPoints, pPolygon, PTOutput) ;
+
+ for ( auto const & tri : PTOutput.interiorTriangles) {
+ vTr.push_back( tri[0]) ;
+ vTr.push_back( tri[1]) ;
+ vTr.push_back( tri[2]) ;
+ }
+
+ return true ;
+}
+
//----------------------------------------------------------------------------
bool
diff --git a/Triangulate.h b/Triangulate.h
index 2f630a6..cce6d13 100644
--- a/Triangulate.h
+++ b/Triangulate.h
@@ -18,9 +18,10 @@
//
enum TrgType { TRG_STANDARD, // ear clipping
- TRG_DEL_CONFORMING, // conforming Delaunay ( with quality constraint)
- TRG_DEL_QUALITY, // Delaunay with quality constraints ( no angle smaller than 20 degrees)
- TRG_DEL_NOQUALITY // Delaunay without quality constraints
+ TRG_DEL_CONFORMING, // conforming constrained Delaunay ( with quality constraint)
+ TRG_DEL_QUALITY, // constrained Delaunay with quality constraints ( no angle smaller than 20 degrees)
+ TRG_DEL_NOQUALITY, // constrained Delaunay without quality constraints
+ TRG_DEL_GTE // constrained Delaunay with Geomtric Tools Engine library
} ;
//----------------------------------------------------------------------------
@@ -37,6 +38,7 @@ class Triangulate
bool MakeByEC2( const PNTVECTOR& vPt, const INTVECTOR& vPol, INTVECTOR& vTr) ;
bool MakeByEC3( const PNTVECTOR& vPt, const INTVECTOR& vPol, INTVECTOR& vTr) ;
bool MakeByDelaunay( const POLYLINEVECTOR& vPL, PNTVECTOR& vPt, INTVECTOR& vTr, bool bConforming, bool bQuality) ;
+ bool MakeByDelaunayGTE( const POLYLINEVECTOR& vPL, PNTVECTOR& vPt, INTVECTOR& vTr) ;
bool TestTriangle( const PNTVECTOR& vPt, const INTVECTOR& vPol,
const INTVECTOR& vPrev, INTVECTOR& vNext, int i) ;
double CalcTriangleAspectRatio( const Point3d& ptPa, const Point3d& ptPb, const Point3d& ptPc) ;
diff --git a/triangle_impl.hpp b/triangle_impl.hpp
index fc5b2cf..e7f455f 100644
--- a/triangle_impl.hpp
+++ b/triangle_impl.hpp
@@ -10156,6 +10156,7 @@ struct behavior *b;
sortarray[j][0], sortarray[j][1]);
}
setvertextype(sortarray[j], UNDEADVERTEX);
+ if ( b->poly) setvertex2tri( sortarray[j], NULL) ;
m->undeads++;
} else {
i++;