diff --git a/GdbExecutor.cpp b/GdbExecutor.cpp index d711108..ebc67bb 100644 --- a/GdbExecutor.cpp +++ b/GdbExecutor.cpp @@ -812,7 +812,7 @@ GdbExecutor::ExecuteSurfTriMesh( const string& sCmd2, const STRVECTOR& vsParams) pModCrv->ToGlob( frCrv) ; pModCrv->ToLoc( frDest) ; // ricavo l'approssimazione - if ( ! pCrv->ApproxWithLines( 0.01, 15, PL)) + if ( ! pModCrv->ApproxWithLines( 0.01, 15, PL)) return false ; } // creo la superficie @@ -861,7 +861,7 @@ GdbExecutor::ExecuteSurfTriMesh( const string& sCmd2, const STRVECTOR& vsParams) pModCrv->ToGlob( frCrv) ; pModCrv->ToLoc( frDest) ; // ricavo l'approssimazione - if ( ! pCrv->ApproxWithLines( 0.01, 15, PL)) + if ( ! pModCrv->ApproxWithLines( 0.01, 15, PL)) return false ; } // recupero il vettore di estrusione diff --git a/Triangulate.cpp b/Triangulate.cpp index a18f412..0dbb307 100644 --- a/Triangulate.cpp +++ b/Triangulate.cpp @@ -32,47 +32,47 @@ Triangulate::Make( const PolyLine& PL, PNTVECTOR& vPt, INTVECTOR& vTr) Plane3d plPlane ; if ( ! PL.IsPlanar( plPlane, 50 * EPS_SMALL)) return false ; - int nPlane ; bool bCCW ; if ( fabs( plPlane.vtN.z) >= fabs( plPlane.vtN.x) && fabs( plPlane.vtN.z) >= fabs( plPlane.vtN.y)) { - nPlane = PL_XY ; + m_nPlane = PL_XY ; bCCW = ( plPlane.vtN.z > 0) ; } else if ( fabs( plPlane.vtN.x) >= fabs( plPlane.vtN.y)) { - nPlane = PL_YZ ; + m_nPlane = PL_YZ ; bCCW = ( plPlane.vtN.x > 0) ; } else { - nPlane = PL_ZX ; + m_nPlane = PL_ZX ; bCCW = ( plPlane.vtN.y > 0) ; } // riempio il vettore con i vertici del poligono da triangolare vPt.clear() ; vPt.reserve( PL.GetPointNbr() - 1) ; + // salto il primo punto (coincide con l'ultimo) + Point3d ptP ; + if ( ! PL.GetFirstPoint( ptP)) + return false ; + // inserisco i punti + while ( PL.GetNextPoint( ptP)) + vPt.push_back( ptP) ; + // creo il vettore degli indici del Poligono + INTVECTOR vPol ; + int n = int( vPt.size()) ; + vPol.reserve( n) ; // se orientato correttamente (componente di N > 0) if ( bCCW) { - // salto il primo punto (coincide con l'ultimo) - Point3d ptP ; - if ( ! PL.GetFirstPoint( ptP)) - return false ; - // inserisco i punti - while ( PL.GetNextPoint( ptP)) - vPt.push_back( ptP) ; + for ( int i = 0 ; i < n ; ++ i) + vPol.push_back( i) ; } // altrimenti lo prendo al contrario else { - // salto l'ultimo punto (coincide con il primo) - Point3d ptP ; - if ( ! PL.GetLastPoint( ptP)) - return false ; - // inserisco i punti - while ( PL.GetPrevPoint( ptP)) - vPt.push_back( ptP) ; + for ( int i = n - 1 ; i >= 0 ; -- i) + vPol.push_back( i) ; } - // eseguo la triangolazione - return MakeByEC( vPt, nPlane, vTr) ; + // eseguo la triangolazione + return MakeByEC2( vPt, vPol, vTr) ; } //---------------------------------------------------------------------------- @@ -80,89 +80,192 @@ Triangulate::Make( const PolyLine& PL, PNTVECTOR& vPt, INTVECTOR& vTr) // Ear Clipping algorithm //---------------------------------------------------------------------------- bool -Triangulate::MakeByEC( const PNTVECTOR& vPt, int nPlane, INTVECTOR& vTr) +Triangulate::MakeByEC( const PNTVECTOR& vPt, const INTVECTOR& vPol, INTVECTOR& vTr) { + // Clear triangle vector + vTr.clear() ; + // At least 3 points - int n = int( vPt.size()) ; + int n = int( vPol.size()) ; if ( n < 3) return false ; - // Save the plane - if ( nPlane != PL_XY && nPlane != PL_YZ && nPlane != PL_ZX) - return false ; - m_nPlane = nPlane ; - - // Clear and preallocate triangle vector ( #triangles = n - 2) - vTr.clear() ; + // Preallocate triangle vector ( #triangles = n - 2) vTr.reserve( 3 * ( n - 2)) ; // Set up previous and next links to effectively form a double-linked vertex list - INTVECTOR prev( n) ; - INTVECTOR next( n) ; - for ( int i = 0 ; i < n ; ++ i) { - prev[i] = i - 1 ; - next[i] = i + 1 ; + INTVECTOR vPrev( n) ; + INTVECTOR vNext( n) ; + for ( int j = 0 ; j < n ; ++ j) { + vPrev[j] = j - 1 ; + vNext[j] = j + 1 ; } - prev[0] = n - 1 ; - next[n - 1] = 0 ; + vPrev[0] = n - 1 ; + vNext[n-1] = 0 ; // Start at vertex 0 int i = 0 ; int nCount = n ; // Keep removing vertices until just a triangle left - while ( n >= 3) { + while ( n > 3) { // To avoid infinite loop if ( nCount <= 0) return false ; -- nCount ; // Test if current vertex, v[i], is an ear - bool bIsEar = true ; - // An ear must be convex (here counterclockwise) - if ( TriangleIsCCW( vPt[prev[i]], vPt[i], vPt[next[i]])) { - // Loop over all vertices not part of the tentative ear - int k = next[next[i]] ; - do { - // If vertex k is inside the ear triangle, then this is not an ear - if ( TestPointInTriangle( vPt[k], vPt[prev[i]], vPt[i], vPt[next[i]])) { - bIsEar = false ; - break; - } - k = next[k] ; - } while (k != prev[i]) ; - } - else { - // The ‘ear’ triangle is clockwise so v[i] is not an ear - bIsEar = false ; - } - + bool bIsEar = TestTriangle( vPt, vPol, vPrev, vNext, i) ; // If current vertex v[i] is an ear, delete it and visit the previous vertex if ( bIsEar) { // Triangle (v[prev[i]], v[i], v[next[i]]) is an ear - vTr.push_back( prev[i]) ; - vTr.push_back( i) ; - vTr.push_back( next[i]) ; + vTr.push_back( vPol[vPrev[i]]) ; + vTr.push_back( vPol[i]) ; + vTr.push_back( vPol[vNext[i]]) ; // ‘Delete’ vertex v[i] by redirecting next and previous links // of neighboring verts past it. Decrement vertex count - next[prev[i]] = next[i] ; - prev[next[i]] = prev[i] ; + vNext[vPrev[i]] = vNext[i] ; + vPrev[vNext[i]] = vPrev[i] ; n--; // Visit the previous vertex next - i = prev[i] ; + i = vPrev[i] ; // Reset Count nCount = n ; } else { // Current vertex is not an ear; visit the next vertex - i = next[i] ; + i = vNext[i] ; } } + // Last triangle is an ear + vTr.push_back( vPol[vPrev[i]]) ; + vTr.push_back( vPol[i]) ; + vTr.push_back( vPol[vNext[i]]) ; + + return true ; +} + +//---------------------------------------------------------------------------- +// Triangulate the CCW n-gon specified by the vertices vPt (Pt[n] != Pt[0]) +// Ear Clipping algorithm enhanced +//---------------------------------------------------------------------------- +bool +Triangulate::MakeByEC2( const PNTVECTOR& vPt, const INTVECTOR& vPol, INTVECTOR& vTr) +{ + // Clear triangle vector + vTr.clear() ; + + // At least 3 points + int n = int( vPol.size()) ; + if ( n < 3) + return false ; + + // Preallocate triangle vector ( #triangles = n - 2) + vTr.reserve( 3 * ( n - 2)) ; + + // Set up previous and next links to effectively form a double-linked vertex list + INTVECTOR vPrev( n) ; + INTVECTOR vNext( n) ; + for ( int j = 0 ; j < n ; ++ j) { + vPrev[j] = j - 1 ; + vNext[j] = j + 1 ; + } + vPrev[0] = n - 1 ; + vNext[n-1] = 0 ; + + // Start at vertex 0 + int i = 0 ; + int nCount = n ; + // Keep removing vertices until just a triangle left + while ( n > 3) { + // To avoid infinite loop + if ( nCount <= 0) + return false ; + -- nCount ; + // Test if current vertex, v[i], is an ear + bool bIsEar = TestTriangle( vPt, vPol, vPrev, vNext, i) ; + if ( bIsEar) { + // Save square distance of diagonal + double dSqDist = SqDist(vPt[vPol[vPrev[i]]], vPt[vPol[vNext[i]]]) ; + // Try with next + int j = vNext[i] ; + bool bIsEar1 = TestTriangle( vPt, vPol, vPrev, vNext, j) ; + double dSqDist1 = INFINITO ; + if ( bIsEar1) + dSqDist1 = SqDist( vPt[vPol[vPrev[j]]], vPt[vPol[vNext[j]]]) ; + // Try with prev + int k = vPrev[i] ; + bool bIsEar2 = TestTriangle( vPt, vPol, vPrev, vNext, k) ; + double dSqDist2 = INFINITO ; + if ( bIsEar2) + dSqDist2 = SqDist( vPt[vPol[vPrev[k]]], vPt[vPol[vNext[k]]]) ; + // Choose the best + if ( dSqDist <= dSqDist1 && dSqDist <= dSqDist2) + ; // i is the better + else if ( dSqDist1 <= dSqDist2) + i = j ; + else + i = k ; + } + + // If current vertex v[i] is an ear, delete it and visit the previous vertex + if ( bIsEar) { + // Triangle (v[prev[i]], v[i], v[next[i]]) is an ear + vTr.push_back( vPol[vPrev[i]]) ; + vTr.push_back( vPol[i]) ; + vTr.push_back( vPol[vNext[i]]) ; + // ‘Delete’ vertex v[i] by redirecting next and previous links + // of neighboring verts past it. Decrement vertex count + vNext[vPrev[i]] = vNext[i] ; + vPrev[vNext[i]] = vPrev[i] ; + n--; + // Visit the previous vertex next + i = vPrev[i] ; + // Reset Count + nCount = n ; + } + else { + // Current vertex is not an ear; visit the next vertex + i = vNext[i] ; + } + } + // Last triangle is an ear + vTr.push_back( vPol[vPrev[i]]) ; + vTr.push_back( vPol[i]) ; + vTr.push_back( vPol[vNext[i]]) ; return true ; } //---------------------------------------------------------------------------- bool -Triangulate::TriangleIsCCW( const Point3d& ptA, const Point3d& ptB, const Point3d& ptC) +Triangulate::TestTriangle( const PNTVECTOR& vPt, const INTVECTOR& vPol, + const INTVECTOR& vPrev, INTVECTOR& vNext, int i) +{ + // Test if current vertex, v[i], is an ear + bool bIsEar = true ; + // An ear must be convex (here counterclockwise) + if ( TriangleIsCCW( vPt[vPol[vPrev[i]]], vPt[vPol[i]], vPt[vPol[vNext[i]]])) { + // Loop over all vertices not part of the tentative ear + int k = vNext[vNext[i]] ; + do { + // If vertex k is inside the ear triangle, then this is not an ear + if ( TestPointInTriangle( vPt[vPol[k]], vPt[vPol[vPrev[i]]], vPt[vPol[i]], vPt[vPol[vNext[i]]])) { + bIsEar = false ; + break ; + } + k = vNext[k] ; + } while (k != vPrev[i]) ; + } + else { + // The ‘ear’ triangle is clockwise so v[i] is not an ear + bIsEar = false ; + } + + return bIsEar ; +} + +//---------------------------------------------------------------------------- +bool +Triangulate::TriangleIsCCW( const Point3d& ptA, const Point3d& ptB, const Point3d& ptC, double dToler) { double dV11 ; double dV12 ; @@ -190,7 +293,7 @@ Triangulate::TriangleIsCCW( const Point3d& ptA, const Point3d& ptB, const Point3 break ; } - return ( ( dV11 * dV22 - dV12 * dV21) > 0) ; + return ( ( dV11 * dV22 - dV12 * dV21) > dToler * dToler) ; } //---------------------------------------------------------------------------- @@ -198,11 +301,19 @@ Triangulate::TriangleIsCCW( const Point3d& ptA, const Point3d& ptB, const Point3 bool Triangulate::TestPointInTriangle( const Point3d& ptP, const Point3d& ptA, const Point3d& ptB, const Point3d& ptC) { - if ( ! TriangleIsCCW( ptA, ptB, ptP)) + // if P is on a vertex is considered outside + if ( AreSamePointNear( ptP, ptA)) return false ; - if ( ! TriangleIsCCW( ptB, ptC, ptP)) + if ( AreSamePointNear( ptP, ptB)) return false ; - if ( ! TriangleIsCCW( ptC, ptA, ptP)) + if ( AreSamePointNear( ptP, ptC)) + return false ; + // if P is on the right of at least one edge is outside + if ( TriangleIsCCW( ptA, ptP, ptB)) + return false ; + if ( TriangleIsCCW( ptB, ptP, ptC)) + return false ; + if ( TriangleIsCCW( ptC, ptP, ptA)) return false ; return true ; } diff --git a/Triangulate.h b/Triangulate.h index 24ac50c..01cb650 100644 --- a/Triangulate.h +++ b/Triangulate.h @@ -24,10 +24,13 @@ class Triangulate { public : bool Make( const PolyLine& PL, PNTVECTOR& vPt, INTVECTOR& vTr) ; - bool MakeByEC( const PNTVECTOR& vPt, int nPlane, INTVECTOR& vTr) ; private : - bool TriangleIsCCW( const Point3d& ptA, const Point3d& ptB, const Point3d& ptC) ; + bool MakeByEC( const PNTVECTOR& vPt, const INTVECTOR& vPol, INTVECTOR& vTr) ; + bool MakeByEC2( const PNTVECTOR& vPt, const INTVECTOR& vPol, INTVECTOR& vTr) ; + bool TestTriangle( const PNTVECTOR& vPt, const INTVECTOR& vPol, + const INTVECTOR& vPrev, INTVECTOR& vNext, int i) ; + bool TriangleIsCCW( const Point3d& ptA, const Point3d& ptB, const Point3d& ptC, double dToler = EPS_SMALL) ; bool TestPointInTriangle( const Point3d& ptP, const Point3d& ptA, const Point3d& ptB, const Point3d& ptC) ; private :