From e0a350b1484b67162abd260596a9ef325bfab44c Mon Sep 17 00:00:00 2001 From: DarioS Date: Tue, 11 Oct 2022 12:27:24 +0200 Subject: [PATCH] EgtGeomKernel 2.4i4 : - a superficie trimesh aggiunta CloneTriangle - modifiche a Triangulate per gestione MinAngle di triangolo. --- EgtGeomKernel.rc | Bin 11718 -> 11718 bytes SurfTriMesh.cpp | 35 ++++++++++++ SurfTriMesh.h | 5 +- Triangulate.cpp | 136 +++++++++++++++++++++++++++++++++-------------- Triangulate.h | 5 +- 5 files changed, 137 insertions(+), 44 deletions(-) diff --git a/EgtGeomKernel.rc b/EgtGeomKernel.rc index d36e9ac8f562b85936d4f12d245ebcc14159eb66..c432d0de78c31488b0fbc01f56ff1d1db9655970 100644 GIT binary patch delta 100 zcmX>WeJpyzFE&P#&G*I5Gfl4J(wY2?Q;pGNvLd(cW-qP@EMOU1W|)YS2jk{G;Q$~X Qr~yeBsu05nEu}f!0K9}B{r~^~ delta 100 zcmX>WeJpyzFE&Qw&G*I5Gfl4J(wY2?Q;pGhvLd(cW-qP@EMOU1W|)YS2jk{G;Q$~X Qr~yeBsu05nEu}f!0J`iS_5c6? diff --git a/SurfTriMesh.cpp b/SurfTriMesh.cpp index 6c3419d..b6525d9 100644 --- a/SurfTriMesh.cpp +++ b/SurfTriMesh.cpp @@ -827,6 +827,41 @@ SurfTriMesh::GetTriangleSmoothNormal( int nT, int nV, Vector3d& vtN) const return true ; } +//---------------------------------------------------------------------------- +SurfTriMesh* +SurfTriMesh::CloneTriangle( int nT) const +{ + // verifico validità superficie ed esistenza del triangolo + if ( ! IsValid() || nT < 0 || nT >= GetTriangleSize() || m_vTria[nT].nIdVert[0] == SVT_DEL) + return nullptr ; + + // Creo nuovo oggetto SurfTriMesh + PtrOwner pSurfTM( new( nothrow) SurfTriMesh) ; + if ( IsNull( pSurfTM)) + return nullptr ; + + // Copio il valore dei membri + pSurfTM->m_dLinTol = m_dLinTol ; + pSurfTM->m_dBoundaryAng = m_dBoundaryAng ; + pSurfTM->m_dCosBndAng = m_dCosBndAng ; + pSurfTM->m_dSmoothAng = m_dSmoothAng ; + pSurfTM->m_dCosSmAng = m_dCosSmAng ; + pSurfTM->m_nParts = 1 ; + + // Copio il triangolo + int nNewInd[3] = { pSurfTM->AddVertex( m_vVert[m_vTria[nT].nIdVert[0]].ptP), + pSurfTM->AddVertex( m_vVert[m_vTria[nT].nIdVert[1]].ptP), + pSurfTM->AddVertex( m_vVert[m_vTria[nT].nIdVert[2]].ptP)} ; + if ( pSurfTM->AddTriangle( nNewInd, m_vTria[nT].nTFlag) == SVT_NULL) + return nullptr ; + + // Aggiusto la superficie + pSurfTM->DoCompacting() ; + + // Restituisco la nuova superficie + return Release( pSurfTM) ; +} + //---------------------------------------------------------------------------- bool SurfTriMesh::GetLoops( POLYLINEVECTOR& vPL) const diff --git a/SurfTriMesh.h b/SurfTriMesh.h index 3f9b8f0..f5fa0d2 100644 --- a/SurfTriMesh.h +++ b/SurfTriMesh.h @@ -1,7 +1,7 @@ //---------------------------------------------------------------------------- -// EgalTech 2014-2021 +// EgalTech 2014-2022 //---------------------------------------------------------------------------- -// File : SurfTriMesh.h Data : 26.09.21 Versione : 2.3i2 +// File : SurfTriMesh.h Data : 10.10.22 Versione : 2.4i4 // Contenuto : Dichiarazione della classe Superficie TriMesh. // // @@ -258,6 +258,7 @@ class SurfTriMesh : public ISurfTriMesh, public IGeoObjRW bool GetVertexSmoothNormal( int nV, int nT, Vector3d& vtN) const override ; bool GetTriangleBoundaryEdges( int nId, TriFlags3d& TFlags) const override ; bool GetTriangleSmoothNormals( int nId, TriNormals3d& TNrms) const override ; + SurfTriMesh* CloneTriangle( int nT) const override ; bool GetLoops( POLYLINEVECTOR& vPL) const override ; bool GetSilhouette( const Vector3d& vtDir, double dTol, POLYLINEVECTOR& vPL) const override ; int GetFacetCount( void) const override ; diff --git a/Triangulate.cpp b/Triangulate.cpp index 279de26..27e3684 100644 --- a/Triangulate.cpp +++ b/Triangulate.cpp @@ -20,6 +20,7 @@ #include "/EgtDev/Include/EGkPolyLine.h" #include "/EgtDev/Include/EGkPlane3d.h" #include "/EgtDev/Include/EGkStringUtils3d.h" +#include "/EgtDev/Include/EgtNumUtils.h" #include using namespace std ; @@ -73,10 +74,17 @@ Triangulate::Make( const PolyLine& PL, PNTVECTOR& vPt, INTVECTOR& vTr) reverse( vPt.begin(), vPt.end()) ; // eseguo la triangolazione - if ( ! MakeByEC2( vPt, vTr) && - ! MakeByEC3( vPt, vTr)) { - LOG_INFO( GetEGkLogger(), "Info : problems in MakeByEC23(1)") - return MakeByEC_HPP( PL, bCCW, vPt, vTr) ; + 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 @@ -187,10 +195,17 @@ Triangulate::Make( const POLYLINEVECTOR& vPL, PNTVECTOR& vPt, INTVECTOR& vTr) return MakeByEC_HPP( vPL, bCCW, vPt, vTr) ; // eseguo la triangolazione - if ( ! MakeByEC2( vPt, vTr) && - ! MakeByEC3( vPt, vTr)) { - LOG_INFO( GetEGkLogger(), "Info : problems in MakeByEC23(N)") - return MakeByEC_HPP( vPL, bCCW, vPt, vTr) ; + 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 @@ -393,10 +408,11 @@ Triangulate::MakeByEC( const PNTVECTOR& vPt, INTVECTOR& vTr) // Ear Clipping algorithm enhanced, choose smaller diagonal //---------------------------------------------------------------------------- bool -Triangulate::MakeByEC2( const PNTVECTOR& vPt, INTVECTOR& vTr) +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()) ; @@ -453,11 +469,11 @@ Triangulate::MakeByEC2( const PNTVECTOR& vPt, INTVECTOR& vTr) int nLimN = min( 19, n / 2) ; for ( int h = 0 ; h < nLimN ; ++ h) { j = vNext[j] ; - if ( vEar[j] == EAS_NULL) - vEar[j] = ( TestTriangle( vPt, vPol, vPrev, vNext, j) ? EAS_OK : EAS_NO) ; - if ( vEar[j] == EAS_OK) { - double dMySqDist = SqDist( vPt[vPol[vPrev[j]]], vPt[vPol[vNext[j]]]) ; - if ( dMySqDist < 0.9 * dSqDist) { + 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 ; } @@ -467,11 +483,11 @@ Triangulate::MakeByEC2( const PNTVECTOR& vPt, INTVECTOR& vTr) int nLimP = min( 19, n / 2) ; for ( int h = 0 ; h < nLimP ; ++ h) { k = vPrev[k] ; - if ( vEar[k] == EAS_NULL) - vEar[k] = ( TestTriangle( vPt, vPol, vPrev, vNext, k) ? EAS_OK : EAS_NO) ; - if ( vEar[k] == EAS_OK) { - double dMySqDist = SqDist( vPt[vPol[vPrev[k]]], vPt[vPol[vNext[k]]]) ; - if ( dMySqDist < 0.9 * dSqDist) { + 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 ; } @@ -491,6 +507,8 @@ Triangulate::MakeByEC2( const PNTVECTOR& vPt, INTVECTOR& vTr) 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 ; @@ -516,14 +534,14 @@ Triangulate::MakeByEC2( const PNTVECTOR& vPt, INTVECTOR& vTr) //---------------------------------------------------------------------------- // Triangulate the CCW n-gon specified by the vertices vPt (Pt[n] != Pt[0]) -// Ear Clipping algorithm enhanced, choose smaller triangle aspect ratio -// AR = ( Lmax * Lmax) / ( 2 * Area) = SqLenMax / L1 * L2 +// Ear Clipping algorithm enhanced, choose max MinAng //---------------------------------------------------------------------------- bool -Triangulate::MakeByEC3( const PNTVECTOR& vPt, INTVECTOR& vTr) +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()) ; @@ -570,32 +588,50 @@ Triangulate::MakeByEC3( const PNTVECTOR& vPt, INTVECTOR& vTr) 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 & Aspect Ratio + // Square Length & MinAngle double dSqLen = SqDist(vPt[vPol[vPrev[i]]], vPt[vPol[vNext[i]]]) ; - double dAR = CalcTriangleAspectRatio( vPt[vPol[vPrev[i]]], vPt[vPol[i]], vPt[vPol[vNext[i]]]) ; - // Try with all other vertices - int ii = i ; - int j = vNext[i] ; - while ( j != ii) { - // New Square Length - double dNewSqLen = SqDist(vPt[vPol[vPrev[j]]], vPt[vPol[vNext[j]]]) ; - if ( dNewSqLen < 0.99 * dSqLen) { - // New Aspect Ratio - double dNewAR = CalcTriangleAspectRatio( vPt[vPol[vPrev[j]]], vPt[vPol[j]], vPt[vPol[vNext[j]]]) ; - // if better, test if triangle ok - if ( dNewAR < 30 || dNewAR < 1.2 * dAR) { + 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 = dNewSqLen ; - dAR = min( dAR, dNewAR) ; + dSqLen = dMySqLen ; + dMinAng = dMyMinAng ; i = j ; } } } - j = vNext[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 { @@ -610,6 +646,7 @@ Triangulate::MakeByEC3( const PNTVECTOR& vPt, INTVECTOR& vTr) 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 ; @@ -716,7 +753,7 @@ Triangulate::TestTriangle( const PNTVECTOR& vPt, const INTVECTOR& vPol, } //---------------------------------------------------------------------------- -// Ratio between max side and opposite height +// Ratio between Circumcircle Radius and Incircle Diameter (minimum 1 for equilateral) double Triangulate::CalcTriangleAspectRatio( const Point3d& ptPa, const Point3d& ptPb, const Point3d& ptPc) { @@ -730,6 +767,25 @@ Triangulate::CalcTriangleAspectRatio( const Point3d& ptPa, const Point3d& ptPb, 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) diff --git a/Triangulate.h b/Triangulate.h index ac19507..d5470b7 100644 --- a/Triangulate.h +++ b/Triangulate.h @@ -27,13 +27,14 @@ class Triangulate bool MakeByEC_HPP( const PolyLine& PL, bool bCCW, PNTVECTOR& vPt, INTVECTOR& vTr) ; bool MakeByEC_HPP( const POLYLINEVECTOR& vPL, bool bCCW, PNTVECTOR& vPt, INTVECTOR& vTr) ; bool MakeByEC( const PNTVECTOR& vPt, INTVECTOR& vTr) ; - bool MakeByEC2( const PNTVECTOR& vPt, INTVECTOR& vTr) ; - bool MakeByEC3( const PNTVECTOR& vPt, INTVECTOR& vTr) ; + bool MakeByEC2( const PNTVECTOR& vPt, INTVECTOR& vTr, double& dMinMinAng) ; + bool MakeByEC3( const PNTVECTOR& vPt, INTVECTOR& vTr, double& dMinMinAng) ; bool PrepareGrid( const PNTVECTOR& vPt, const INTVECTOR& vPol, const INTVECTOR& vPrev, const INTVECTOR& vNext) ; 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) ; + double CalcTriangleMinAngle( const Point3d& ptPa, const Point3d& ptPb, const Point3d& ptPc) ; double SquareDist( const Point3d& ptA, const Point3d& ptB) ; double TwoArea( const Point3d& ptA, const Point3d& ptB, const Point3d& ptC) ; bool AreSamePoint( const Point3d& ptA, const Point3d& ptB, double dToler = EPS_SMALL) ;