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

EgtGeomKernel 2.4i4 :

Browse Source 
- a superficie trimesh aggiunta CloneTriangle
- modifiche a Triangulate per gestione MinAngle di triangolo.
This commit is contained in:
DarioS
2022-10-11 12:27:24 +02:00
parent 560a447349
commit e0a350b148
5 changed files with 137 additions and 44 deletions
Download Patch File Download Diff File Expand all files Collapse all files
Triangulate.cpp
+96 -40
View File
@@ -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 <algorithm>
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)