Egalware/EgtGeomKernel
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

EgtGeomKernel :

Browse Source 
- aggiunta triangolazione constrained Delaunay con libreria GeometricToolsEngine
- correzione errori triangolazione con TrianglePP.
This commit is contained in:
SaraP
2021-07-09 16:05:59 +02:00
parent 03815d7747
commit 8ac3fdab47
4 changed files with 254 additions and 22 deletions
Download Patch File Download Diff File Expand all files Collapse all files
Triangulate.cpp
+246 -19
View File
@@ -24,10 +24,12 @@
#include "/EgtDev/Include/EGkPolyLine.h"
#include "/EgtDev/Include/EGkPlane3d.h"
#include "/EgtDev/Include/EGkStringUtils3d.h"
#include "Mathematics/TriangulateCDT.h"
#include <algorithm>
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<Delaunay::Point> 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<Delaunay::Point> 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<Vector2<double>> 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<pair<double,int>> 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<double, int> (aInfo.IciA[0].dU, it - vPt.begin())) ;
// dParams.push_back( pair<double, int> (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<double, int> (aInfo.IciA[0].dU, vPoints.size() - 1)) ;
// dParams.push_back( pair<double, int> (aInfo.IciB[0].dU, vPoints.size() - 1)) ;
// }
// }
// sort( dParams.begin(), dParams.end(), []( const pair<double, int>& a, const pair<double, int>& 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<PolygonTree>() ;
// pPolygon->polygon = vPolygonIdx ;
//
// // holes
// int k = vPoints.size() ;
// for ( size_t i = 1 ; i < vMyPL.size() ; i++) {
//
// auto pHole = std::make_shared<PolygonTree>() ;
// 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<double> 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<Vector2<double>> vPoints ; // vertici per la triangolazione
auto pPolygon = std::make_shared<PolygonTree>() ; // 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<PolygonTree>() ;
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<double> 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