From 8ac3fdab478fffa2dc3dd10ae3bcdb7fab211c58 Mon Sep 17 00:00:00 2001 From: SaraP Date: Fri, 9 Jul 2021 16:05:59 +0200 Subject: [PATCH] EgtGeomKernel : - aggiunta triangolazione constrained Delaunay con libreria GeometricToolsEngine - correzione errori triangolazione con TrianglePP. --- EgtGeomKernel.vcxproj | 2 + Triangulate.cpp | 265 +++++++++++++++++++++++++++++++++++++++--- Triangulate.h | 8 +- triangle_impl.hpp | 1 + 4 files changed, 254 insertions(+), 22 deletions(-) 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++;