diff --git a/ArcXxTgArc.cpp b/ArcXxTgArc.cpp index 0c9b527..f098c3f 100644 --- a/ArcXxTgArc.cpp +++ b/ArcXxTgArc.cpp @@ -266,7 +266,7 @@ CalcCircleCenTgCircle( const Point3d& ptC, const Point3d& ptCen, double dRad, BI if ( dRad < EPS_SMALL) { // se la distanza tra i punti è significativa, c'è una soluzione if ( dDist > EPS_SMALL) { - vCenPtg.push_back( make_pair( ptC, ptCen)) ; + vCenPtg.emplace_back( ptC, ptCen) ; return 1 ; } // altrimenti, nessuna soluzione @@ -279,13 +279,13 @@ CalcCircleCenTgCircle( const Point3d& ptC, const Point3d& ptCen, double dRad, BI return 0 ; // se è minore del raggio, c'è una sola soluzione else if ( dDist < dRad) { - vCenPtg.push_back( make_pair( ptC, ptCen - vtDir * dRad)) ; + vCenPtg.emplace_back( ptC, ptCen - vtDir * dRad) ; return 1 ; } // altrimenti ci sono due soluzioni else { - vCenPtg.push_back( make_pair( ptC, ptCen - vtDir * dRad)) ; - vCenPtg.push_back( make_pair( ptC, ptCen + vtDir * dRad)) ; + vCenPtg.emplace_back( ptC, ptCen - vtDir * dRad) ; + vCenPtg.emplace_back( ptC, ptCen + vtDir * dRad) ; return 2 ; } } diff --git a/ChainCurves.cpp b/ChainCurves.cpp index 761eb47..8e902da 100644 --- a/ChainCurves.cpp +++ b/ChainCurves.cpp @@ -44,7 +44,7 @@ ChainCurves::AddCurve( int nId, const Point3d& ptStart, const Vector3d& vtStart, if ( ! m_sCrvId.insert( nId).second) return true ; // inserisco i dati della curva nel vettore - m_vCrvData.push_back( CrvData( nId, ptStart, vtStart, ptEnd, vtEnd)) ; + m_vCrvData.emplace_back( nId, ptStart, vtStart, ptEnd, vtEnd) ; int nV = int( m_vCrvData.size()) ; // li inserisco nel PointGrid3d m_PointGrid.InsertPoint( ptStart, nV) ; diff --git a/CreateCurveAux.cpp b/CreateCurveAux.cpp index 42ddc36..e5251a4 100644 --- a/CreateCurveAux.cpp +++ b/CreateCurveAux.cpp @@ -39,7 +39,7 @@ CalcLinePointTgCircle( const Point3d& ptP, const Point3d& ptCen, double dRad, BI if ( dRad < EPS_SMALL) { // se la distanza tra i punti è significativa, c'è una soluzione : la retta per i due punti if ( dDist > EPS_SMALL) { - vBiPnt.push_back( make_pair( ptP, ptCen)) ; + vBiPnt.emplace_back( ptP, ptCen) ; return 1 ; } // altrimenti, nessuna soluzione @@ -64,12 +64,12 @@ CalcLinePointTgCircle( const Point3d& ptP, const Point3d& ptCen, double dRad, BI // tangente a destra Point3d ptT = ptK + vtOrtho ; ptT.z = ptCen.z ; - vBiPnt.push_back( make_pair( ptP, ptT)) ; + vBiPnt.emplace_back( ptP, ptT) ; // tangente a sinistra ptT = ptK - vtOrtho ; ptT.z = ptCen.z ; - vBiPnt.push_back( make_pair( ptP, ptT)) ; + vBiPnt.emplace_back( ptP, ptT) ; return 2 ; } diff --git a/CurveBezier.cpp b/CurveBezier.cpp index 9ae69bb..6fdc178 100644 --- a/CurveBezier.cpp +++ b/CurveBezier.cpp @@ -1237,7 +1237,7 @@ CurveBezier::GetParamAtLength( double dLen, double& dU) const double dUFin = 1 ; while ( ( bOk = GetSegmentParam( dLen, dCurrLen, dSegLen, dUIni, dUFin)) && dSegLen > MIN_SEG_LEN) { // allungo di un poco l'intervallo di ricerca per compensare l'approssimazione della lunghezza - dUFin = min( dUFin + ( dUFin - dUIni) * 0.05, 1) ; + dUFin = min( dUFin + ( dUFin - dUIni) * 0.05, 1.) ; } // aggiustamento finale parametro if ( bOk) diff --git a/CurveByInterp.cpp b/CurveByInterp.cpp index 3ebb065..a21bee3 100644 --- a/CurveByInterp.cpp +++ b/CurveByInterp.cpp @@ -136,11 +136,11 @@ CurveByInterp::CalcTangents( int nMethod) // se ci sono solo 2 punti, le tangenti devono essere dirette lungo la linea che li unisce if ( nSize == 2) { // non esiste derivata prima del primo punto - m_vPrevDer.push_back( Vector3d( 0, 0, 0)) ; + m_vPrevDer.emplace_back( 0, 0, 0) ; m_vNextDer.push_back( ( m_vPnt[1] - m_vPnt[0]) / ( m_vPar[1] - m_vPar[0])) ; m_vPrevDer.push_back( m_vNextDer[0]) ; // non esiste derivata dopo il secondo e ultimo punto - m_vNextDer.push_back( Vector3d( 0, 0, 0)) ; + m_vNextDer.emplace_back( 0, 0, 0) ; return true ; } // verifico se curva chiusa (primo e ultimo punto coincidono) diff --git a/DistPointArc.cpp b/DistPointArc.cpp index 9b98d60..1f6de07 100644 --- a/DistPointArc.cpp +++ b/DistPointArc.cpp @@ -58,7 +58,7 @@ DistPointArc::DistPointCircle( const Point3d& ptP, const ICurveArc& arArc) // calcolo del valore di minima distanza m_dDist = Dist( ptP, ptMinDist) ; // salvo i dati - m_Info.push_back( MinDistPCInfo( MDPCI_NORMAL, dParam, ptMinDist)) ; + m_Info.emplace_back( MDPCI_NORMAL, dParam, ptMinDist) ; } // altrimenti tutti i punti della circonferenza sono a minima distanza else { @@ -67,9 +67,9 @@ DistPointArc::DistPointCircle( const Point3d& ptP, const ICurveArc& arArc) // salvo iniziale e finale, come estremi del range Point3d ptMinDist ; arArc.GetStartPoint( ptMinDist) ; - m_Info.push_back( MinDistPCInfo( MDPCI_START_CONT, 0, ptMinDist)) ; + m_Info.emplace_back( MDPCI_START_CONT, 0, ptMinDist) ; arArc.GetEndPoint( ptMinDist) ; - m_Info.push_back( MinDistPCInfo( MDPCI_END_CONT, 1, ptMinDist)) ; + m_Info.emplace_back( MDPCI_END_CONT, 1, ptMinDist) ; } } @@ -92,7 +92,7 @@ DistPointArc::DistPointFlatArc( const Point3d& ptP, const ICurveArc& arArc) // altro punto con la stessa minima distanza if ( i == 1 && fabs( dDist - m_dDist) < EPS_SMALL) { // lo aggiungo - m_Info.push_back( MinDistPCInfo( MDPCI_NORMAL, dU, ptTest)) ; + m_Info.emplace_back( MDPCI_NORMAL, dU, ptTest) ; } // primo punto o punto con minima distanza più bassa else if ( i == 0 || dDist < m_dDist) { @@ -100,7 +100,7 @@ DistPointArc::DistPointFlatArc( const Point3d& ptP, const ICurveArc& arArc) m_dDist = dDist ; // il nuovo vettore deve contenere solo quest'ultimo minimo m_Info.clear() ; - m_Info.push_back( MinDistPCInfo( MDPCI_NORMAL, dU, ptTest)) ; + m_Info.emplace_back( MDPCI_NORMAL, dU, ptTest) ; } } } diff --git a/DistPointCrvAux.cpp b/DistPointCrvAux.cpp index 45e09c6..f749693 100644 --- a/DistPointCrvAux.cpp +++ b/DistPointCrvAux.cpp @@ -169,7 +169,7 @@ FilterMinDistPointCurve( const Point3d& ptP, const ICurve& cCurve, // se abbastanza lontano e non su bordi intervallino lo aggiungo if ( SqDist( (*Iter).ptQ, Info.back().ptQ) > MIN_LEN_CONT_MDPC && ! bLastOnEnd && fabs( (*Iter).dPar - (*Iter).dParMin) > EPS_ZERO) { - Info.push_back( MinDistPCInfo( MDPCI_NORMAL, (*Iter).dPar, (*Iter).ptQ)) ; + Info.emplace_back( MDPCI_NORMAL, (*Iter).dPar, (*Iter).ptQ) ; bLastOnEnd = fabs( (*Iter).dPar - (*Iter).dParMax) < EPS_ZERO ; } // altrimenti lo sostituisco se distanza minore @@ -185,7 +185,7 @@ FilterMinDistPointCurve( const Point3d& ptP, const ICurve& cCurve, dMinDist = (*Iter).dDist ; // il nuovo vettore deve contenere solo quest'ultimo minimo Info.clear() ; - Info.push_back( MinDistPCInfo( MDPCI_NORMAL, (*Iter).dPar, (*Iter).ptQ)) ; + Info.emplace_back( MDPCI_NORMAL, (*Iter).dPar, (*Iter).ptQ) ; bLastOnEnd = fabs( (*Iter).dPar - (*Iter).dParMax) < EPS_ZERO ; } } diff --git a/DistPointCurve.cpp b/DistPointCurve.cpp index 1a78825..458c7c6 100644 --- a/DistPointCurve.cpp +++ b/DistPointCurve.cpp @@ -62,7 +62,7 @@ DistPointCurve::LineCalculate( const Point3d& ptP, const ICurve& Curve, bool bIs if ( dstPtLn.m_dSqDist >= 0) { m_dDist = sqrt( dstPtLn.m_dSqDist) ; - m_Info.push_back( MinDistPCInfo( MDPCI_NORMAL, dstPtLn.m_dParam, dstPtLn.m_ptMinDist)) ; + m_Info.emplace_back( MDPCI_NORMAL, dstPtLn.m_dParam, dstPtLn.m_ptMinDist) ; } } diff --git a/EgtGeomKernel.rc b/EgtGeomKernel.rc index 15e1d95..23d25d9 100644 Binary files a/EgtGeomKernel.rc and b/EgtGeomKernel.rc differ diff --git a/EgtGeomKernel.vcxproj b/EgtGeomKernel.vcxproj index c484bd8..6e4ec49 100644 --- a/EgtGeomKernel.vcxproj +++ b/EgtGeomKernel.vcxproj @@ -169,7 +169,7 @@ copy $(TargetPath) \EgtProg\DllD64 Level3 Use - Full + MaxSpeed true true WIN32;_WINDOWS;I_AM_EGK;NDEBUG;%(PreprocessorDefinitions) @@ -179,7 +179,7 @@ copy $(TargetPath) \EgtProg\DllD64 AnySuitable true StreamingSIMDExtensions2 - false + true Precise @@ -210,7 +210,7 @@ copy $(TargetPath) \EgtProg\Dll32 Level3 Use - Full + MaxSpeed true true WIN32;_WINDOWS;I_AM_EGK;NDEBUG;%(PreprocessorDefinitions) @@ -220,7 +220,7 @@ copy $(TargetPath) \EgtProg\Dll32 AnySuitable true NotSet - false + true Precise diff --git a/FontNfe.cpp b/FontNfe.cpp index bb132c1..be39b6e 100644 --- a/FontNfe.cpp +++ b/FontNfe.cpp @@ -360,7 +360,7 @@ NfeFont::ApproxWithLines( const string& sText, int nInsPos, double dLinTol, doub pCrvNew->Scale( GLOB_FRM, dScaX, dScaY, dScaZ) ; pCrvNew->Translate( vtMove) ; // approssimo con linee - lstPL.push_back( PolyLine()) ; + lstPL.emplace_back() ; pCrvNew->ApproxWithLines( dLinTol, dAngTolDeg, lstPL.back()) ; } bIter = pIter->GoToNext() ; @@ -443,7 +443,7 @@ NfeFont::ApproxWithArcs( const string& sText, int nInsPos, double dLinTol, doubl pCrvNew->Scale( GLOB_FRM, dScaX, dScaY, dScaZ) ; pCrvNew->Translate( vtMove) ; // approssimo con archi - lstPA.push_back( PolyArc()) ; + lstPA.emplace_back() ; pCrv->ApproxWithArcs( dLinTol, dAngTolDeg, lstPA.back()) ; } bIter = pIter->GoToNext() ; @@ -507,7 +507,7 @@ NfeFont::GetTextLines( const string& sText, int nInsPos, PNTVECTOR& vPt, STRVECT SetCodePoints( vTmpCode, sLine) ; vLine.push_back( sLine) ; vTmpCode.clear() ; - vPt.push_back( Point3d( 0, vtMove.y, 0)) ; + vPt.emplace_back( 0, vtMove.y, 0) ; } // sistemo la posizione vtMove.Set( 0, vtMove.y - dH, 0) ; @@ -537,7 +537,7 @@ NfeFont::GetTextLines( const string& sText, int nInsPos, PNTVECTOR& vPt, STRVECT SetCodePoints( vTmpCode, sLine) ; vLine.push_back( sLine) ; vTmpCode.clear() ; - vPt.push_back( Point3d( 0, vtMove.y, 0)) ; + vPt.emplace_back( 0, vtMove.y, 0) ; } // sistemazioni per la posizione del punto di inserimento diff --git a/FontOs.cpp b/FontOs.cpp index 1a531ad..785a083 100644 --- a/FontOs.cpp +++ b/FontOs.cpp @@ -374,7 +374,7 @@ OsFont::ApproxWithLines( const string& sText, int nInsPos, double dLinTol, doubl (*iIter)->Scale( GLOB_FRM, dScaX, dScaY, dScaZ) ; (*iIter)->Translate( vtMove) ; // approssimo con rette - lstPL.push_back( PolyLine()) ; + lstPL.emplace_back() ; if ( ! (*iIter)->ApproxWithLines( dLinTol, dAngTolDeg, lstPL.back())) return false ; } @@ -452,7 +452,7 @@ OsFont::ApproxWithArcs( const string& sText, int nInsPos, double dLinTol, double (*iIter)->Scale( GLOB_FRM, dScaX, dScaY, dScaZ) ; (*iIter)->Translate( vtMove) ; // approssimo con archi - lstPA.push_back( PolyArc()) ; + lstPA.emplace_back() ; if ( ! (*iIter)->ApproxWithArcs( dLinTol, dAngTolDeg, lstPA.back())) return false ; } @@ -522,7 +522,7 @@ OsFont::GetTextLines( const string& sText, int nInsPos, PNTVECTOR& vPt, STRVECTO SetCodePoints( vTmpCode, sLine) ; vLine.push_back( sLine) ; vTmpCode.clear() ; - vPt.push_back( Point3d( 0, vtMove.y, 0)) ; + vPt.emplace_back( 0, vtMove.y, 0) ; } // sistemo la posizione vtMove.Set( 0, vtMove.y - dH, 0) ; @@ -549,7 +549,7 @@ OsFont::GetTextLines( const string& sText, int nInsPos, PNTVECTOR& vPt, STRVECTO SetCodePoints( vTmpCode, sLine) ; vLine.push_back( sLine) ; vTmpCode.clear() ; - vPt.push_back( Point3d( 0, vtMove.y, 0)) ; + vPt.emplace_back( 0, vtMove.y, 0) ; } // sistemazioni per la posizione del punto di inserimento diff --git a/FontOs.h b/FontOs.h index 51398ab..4210c21 100644 --- a/FontOs.h +++ b/FontOs.h @@ -15,6 +15,7 @@ #include "/EgtDev/Include/EGkCurve.h" #include "/EgtDev/Include/EGnStringBase.h" +#define NOMINMAX #include #include diff --git a/GdbExecutor.cpp b/GdbExecutor.cpp index 984ca5f..ca21f97 100644 --- a/GdbExecutor.cpp +++ b/GdbExecutor.cpp @@ -2124,7 +2124,7 @@ GdbExecutor::ExecuteSurfTriMesh( const string& sCmd2, const STRVECTOR& vsParams) } // se creazione per triangolazione di un contorno chiuso e piano else if ( sCmd2 == "CONT" || sCmd2 == "BYCONTOUR") { - return SurfTriMeshByContour( vsParams) ; + return SurfTriMeshByFlatContour( vsParams) ; } // se creazione per estrusione else if ( sCmd2 == "EXTR" || sCmd2 == "BYEXTRUSION") { @@ -2321,7 +2321,7 @@ GdbExecutor::SurfTriMeshByTriangleSoup( const STRVECTOR& vsParams) //---------------------------------------------------------------------------- bool -GdbExecutor::SurfTriMeshByContour( const STRVECTOR& vsParams) +GdbExecutor::SurfTriMeshByFlatContour( const STRVECTOR& vsParams) { // 3 o 4 parametri : Id, ParentId, IdCurve[, dLinTol] if ( vsParams.size() != 3 && vsParams.size() != 4) @@ -2340,7 +2340,7 @@ GdbExecutor::SurfTriMeshByContour( const STRVECTOR& vsParams) if ( vsParams.size() == 4) FromString( vsParams[3], dLinTol) ; // calcolo la superficie - ISurfTriMesh* pSTM = GetSurfTriMeshByContour( *CrvLoc.Get(), dLinTol) ; + ISurfTriMesh* pSTM = GetSurfTriMeshByFlatContour( CrvLoc.Get(), dLinTol) ; if ( pSTM == nullptr) return false ; // inserisco la superficie trimesh nel DB @@ -2372,7 +2372,7 @@ GdbExecutor::SurfTriMeshByExtrusion( const STRVECTOR& vsParams) if ( vsParams.size() == 5) FromString( vsParams[4], dLinTol) ; // calcolo la superficie - ISurfTriMesh* pSTM = GetSurfTriMeshByExtrusion( *CrvLoc.Get(), vtExtr, false, dLinTol) ; + ISurfTriMesh* pSTM = GetSurfTriMeshByExtrusion( CrvLoc.Get(), vtExtr, false, dLinTol) ; if ( pSTM == nullptr) return false ; // inserisco la superficie trimesh nel DB @@ -2405,7 +2405,7 @@ GdbExecutor::SurfTriMeshByTwoPaths( const STRVECTOR& vsParams) if ( vsParams.size() == 5) FromString( vsParams[4], dLinTol) ; // calcolo la superficie - ISurfTriMesh* pSTM = GetSurfTriMeshRuled( *CrvLoc1.Get(), *CrvLoc2.Get(), dLinTol) ; + ISurfTriMesh* pSTM = GetSurfTriMeshRuled( CrvLoc1.Get(), CrvLoc2.Get(), dLinTol) ; if ( pSTM == nullptr) return false ; // inserisco la superficie trimesh nel DB @@ -2463,7 +2463,7 @@ GdbExecutor::SurfTriMeshByScrewing( bool bMove, const STRVECTOR& vsParams) if ( ! FromString( vsParams[5], dAngRotDeg)) return false ; // calcolo la superficie - ISurfTriMesh* pSTM = GetSurfTriMeshByScrewing( *CrvLoc.Get(), ptAx, vtAx, dAngRotDeg, dMove, dLinTol) ; + ISurfTriMesh* pSTM = GetSurfTriMeshByScrewing( CrvLoc.Get(), ptAx, vtAx, dAngRotDeg, dMove, dLinTol) ; if ( pSTM == nullptr) return false ; // inserisco la superficie trimesh nel DB diff --git a/GdbExecutor.h b/GdbExecutor.h index 046cb46..4c816b3 100644 --- a/GdbExecutor.h +++ b/GdbExecutor.h @@ -106,7 +106,7 @@ class GdbExecutor : public IGdbExecutor bool SurfTriMeshAddTriangle( const STRVECTOR& vsParams) ; bool SurfTriMeshEnd( const STRVECTOR& vsParams) ; bool SurfTriMeshByTriangleSoup( const STRVECTOR& vsParams) ; - bool SurfTriMeshByContour( const STRVECTOR& vsParams) ; + bool SurfTriMeshByFlatContour( const STRVECTOR& vsParams) ; bool SurfTriMeshByExtrusion( const STRVECTOR& vsParams) ; bool SurfTriMeshByTwoPaths( const STRVECTOR& vsParams) ; bool SurfTriMeshByScrewing( bool bMove, const STRVECTOR& vsParams) ; diff --git a/LineTgTwoCurves.cpp b/LineTgTwoCurves.cpp index 788c6de..155e344 100644 --- a/LineTgTwoCurves.cpp +++ b/LineTgTwoCurves.cpp @@ -252,12 +252,12 @@ CalcLineExtTg2Circles( const Point3d& ptC1, double dRad1, const Point3d& ptC2, d // offset a sinistra ptT1 = ptC1 + vtDir ; ptT2 = ptC2 + vtDir ; - vBiPnt.push_back( make_pair( ptT1, ptT2)) ; + vBiPnt.emplace_back( ptT1, ptT2) ; ++ nSol ; // offset a destra ptT1 = ptC1 - vtDir ; ptT2 = ptC2 - vtDir ; - vBiPnt.push_back( make_pair( ptT1, ptT2)) ; + vBiPnt.emplace_back( ptT1, ptT2) ; ++ nSol ; } else if ( nSol1 == 2) { @@ -269,7 +269,7 @@ CalcLineExtTg2Circles( const Point3d& ptC1, double dRad1, const Point3d& ptC2, d // punti della prima retta ptT1 = ptC1 + vtDir ; ptT2 = vAuxBiPnt[0].second + vtDir ; - vBiPnt.push_back( make_pair( ptT1, ptT2)) ; + vBiPnt.emplace_back( ptT1, ptT2) ; ++ nSol ; // direzione di offset della seconda retta vtDir = vAuxBiPnt[1].second - ptC2 ; @@ -277,7 +277,7 @@ CalcLineExtTg2Circles( const Point3d& ptC1, double dRad1, const Point3d& ptC2, d // punti della seconda retta ptT1 = ptC1 + vtDir ; ptT2 = vAuxBiPnt[1].second + vtDir ; - vBiPnt.push_back( make_pair( ptT1, ptT2)) ; + vBiPnt.emplace_back( ptT1, ptT2) ; ++ nSol ; } @@ -309,7 +309,7 @@ CalcLineIntTg2Circles( const Point3d& ptC1, double dRad1, const Point3d& ptC2, d // punti della prima retta ptT1 = ptC1 + vtDir ; ptT2 = vAuxBiPnt[0].second + vtDir ; - vBiPnt.push_back( make_pair( ptT1, ptT2)) ; + vBiPnt.emplace_back( ptT1, ptT2) ; ++ nSol ; // direzione di offset della seconda retta vtDir = ptC2 - vAuxBiPnt[1].second ; @@ -317,7 +317,7 @@ CalcLineIntTg2Circles( const Point3d& ptC1, double dRad1, const Point3d& ptC2, d // punti della seconda retta ptT1 = ptC1 + vtDir ; ptT2 = vAuxBiPnt[1].second + vtDir ; - vBiPnt.push_back( make_pair( ptT1, ptT2)) ; + vBiPnt.emplace_back( ptT1, ptT2) ; ++ nSol ; } diff --git a/PointGrid3d.cpp b/PointGrid3d.cpp index 6b2fd65..3912d80 100644 --- a/PointGrid3d.cpp +++ b/PointGrid3d.cpp @@ -146,7 +146,7 @@ PointGrid3d::InsertPoint( const Point3d& ptP, int nId) { // inserimento nel map int nKey = PointHash( Get1dCellNbr( ptP.x), Get1dCellNbr( ptP.y), Get1dCellNbr( ptP.z)) ; - m_MMap.insert( make_pair( nKey, make_pair( ptP, nId))) ; + m_MMap.emplace( nKey, make_pair( ptP, nId)) ; // aggiornamento del bbox complessivo m_BBox.Add( ptP) ; return true ; @@ -198,6 +198,32 @@ PointGrid3d::Find( const Point3d& ptTest, double dTol, INTVECTOR& vnIds) return ( vnIds.size() > 0) ; } +//---------------------------------------------------------------------------- +bool +PointGrid3d::Find( const BBox3d& b3Test, INTVECTOR& vnIds) +{ + // pulisco il risultato + vnIds.clear() ; + // determino il range di celle sui tre assi + IBox iBox ; + if ( ! Get3dRangeNbr( b3Test, iBox)) + return false ; + // ciclo su tutte le celle del range + for ( int i = iBox.nXmin ; i <= iBox.nXmax ; ++ i) { + for ( int j = iBox.nYmin ; j <= iBox.nYmax ; ++ j) { + for ( int k = iBox.nZmin ; k <= iBox.nZmax ; ++ k) { + IPNTI_UMMAP_CRANGE MMrange = m_MMap.equal_range( PointHash( i, j, k)) ; + for ( ; MMrange.first != MMrange.second ; ++ MMrange.first) { + if ( b3Test.Encloses( (*MMrange.first).second.first)) + vnIds.push_back( (*MMrange.first).second.second) ; + } + } + } + } + + return ( vnIds.size() > 0) ; +} + //---------------------------------------------------------------------------- bool PointGrid3d::Find( const Point3d& ptTest, double dTol, int& nId) @@ -315,6 +341,24 @@ PointGrid3d::Get1dCellNbr( double dCoord) return static_cast( floor( dCoord * m_dInvCellDim)) ; } +//---------------------------------------------------------------------------- +bool +PointGrid3d::Get3dRangeNbr( const BBox3d& b3Test, IBox& iBox) +{ + // calcolo intersezione tra box + BBox3d b3Int ; + if ( ! m_BBox.FindIntersection( b3Test, b3Int)) + return false ; + // ricavo gli indici sui tre assi degli estremi del box + iBox.nXmin = Get1dCellNbr( b3Int.GetMin().x) ; + iBox.nYmin = Get1dCellNbr( b3Int.GetMin().y) ; + iBox.nZmin = Get1dCellNbr( b3Int.GetMin().z) ; + iBox.nXmax = Get1dCellNbr( b3Int.GetMax().x) ; + iBox.nYmax = Get1dCellNbr( b3Int.GetMax().y) ; + iBox.nZmax = Get1dCellNbr( b3Int.GetMax().z) ; + return true ; +} + //---------------------------------------------------------------------------- bool PointGrid3d::Get3dRangeNbr( const Point3d& ptTest, double dTol, IBox& iBox) @@ -337,5 +381,5 @@ PointGrid3d::Get3dRangeNbr( const Point3d& ptTest, double dTol, IBox& iBox) int PointGrid3d::PointHash( int nX, int nY, int nZ) { - return ( nX * 73856093 ^ nY * 19349663 ^ nZ * 83492791) ; + return ( nX * 73856093 ^ nY * 19349653 ^ nZ * 83492791) ; } diff --git a/PolyArc.cpp b/PolyArc.cpp index b9f6428..b364b92 100644 --- a/PolyArc.cpp +++ b/PolyArc.cpp @@ -64,7 +64,7 @@ PolyArc::AddUPoint( double dPar, const Point3d& ptP, double dBulge) return true ; } try { - m_lUPointBs.push_back( UPointB( dPar, ptP, dBulge)) ; + m_lUPointBs.emplace_back( dPar, ptP, dBulge) ; } catch (...) { return false ; diff --git a/PolyLine.cpp b/PolyLine.cpp index 727bdea..c2c8035 100644 --- a/PolyLine.cpp +++ b/PolyLine.cpp @@ -55,7 +55,7 @@ PolyLine::AddUPoint( double dPar, const Point3d& ptP) } // eseguo inserimento try { - m_lUPoints.push_back( POINTU( ptP, dPar)) ; + m_lUPoints.emplace_back( ptP, dPar) ; } catch (...) { return false ; @@ -258,6 +258,20 @@ PolyLine::Split( double dU, PolyLine& PL) return true ; } +//---------------------------------------------------------------------------- +bool +PolyLine::GetLocalBBox( BBox3d& b3Loc) const +{ + // assegno il box in locale, scorrendo tutti i punti + b3Loc.Reset() ; + for ( PNTULIST::const_iterator iter = m_lUPoints.begin() ; + iter != m_lUPoints.end() ; + ++ iter) + b3Loc.Add( iter->first) ; + + return true ; +} + //---------------------------------------------------------------------------- bool PolyLine::IsClosed( void) const diff --git a/PolynomialPoint3d.cpp b/PolynomialPoint3d.cpp index 36b18b9..6520dd5 100644 --- a/PolynomialPoint3d.cpp +++ b/PolynomialPoint3d.cpp @@ -29,7 +29,7 @@ PolynomialPoint3d::SetDegree( int nDegree) m_Coeff.reserve( nDegree + 1) ; m_nDegree = nDegree ; for ( int i = 0 ; i <= m_nDegree ; ++ i) - m_Coeff.push_back( Point3d( 0, 0, 0)) ; + m_Coeff.emplace_back( 0, 0, 0) ; return true ; } catch (...) { @@ -48,7 +48,7 @@ PolynomialPoint3d::EnsureDegree( int nDegree) try { m_Coeff.reserve( nDegree + 1) ; for ( int i = m_nDegree + 1 ; i <= nDegree ; ++ i) - m_Coeff.push_back( Point3d()) ; + m_Coeff.emplace_back() ; m_nDegree = nDegree ; return true ; } diff --git a/StmFromCurves.cpp b/StmFromCurves.cpp index 0504c5e..6e50707 100644 --- a/StmFromCurves.cpp +++ b/StmFromCurves.cpp @@ -20,16 +20,20 @@ using namespace std ; +//------------------------------------------------------------------------------- +static bool CalcRegionPolyLines( const ICURVEPVECTOR& vpCurve, double dLinTol, + POLYLINEVECTOR& vPL, Vector3d& vtN) ; + //------------------------------------------------------------------------------- ISurfTriMesh* -GetSurfTriMeshByContour( const ICurve& Curve, double dLinTol) +GetSurfTriMeshByFlatContour( const ICurve* pCurve, double dLinTol) { // verifica parametri - if ( &Curve == nullptr) + if ( pCurve == nullptr) return nullptr ; // calcolo la polilinea che approssima la curva PolyLine PL ; - if ( ! Curve.ApproxWithLines( dLinTol, ANG_TOL_STD_DEG, PL)) + if ( ! pCurve->ApproxWithLines( dLinTol, ANG_TOL_STD_DEG, PL)) return nullptr ; // creo e setto la superficie trimesh PtrOwner pSTM( CreateSurfTriMesh()) ; @@ -41,15 +45,35 @@ GetSurfTriMeshByContour( const ICurve& Curve, double dLinTol) //------------------------------------------------------------------------------- ISurfTriMesh* -GetSurfTriMeshByExtrusion( const ICurve& Curve, const Vector3d& vtExtr, +GetSurfTriMeshByRegion( const ICURVEPVECTOR& vpCurve, double dLinTol) +{ + // verifica parametri + if ( &vpCurve == nullptr || vpCurve.empty()) + return nullptr ; + // calcolo le polilinee che approssimano le curve della regione + POLYLINEVECTOR vPL ; + Vector3d vtN ; + if ( ! CalcRegionPolyLines( vpCurve, dLinTol, vPL, vtN)) + return nullptr ; + // creo e setto la superficie trimesh + PtrOwner pSTM( CreateSurfTriMesh()) ; + if ( IsNull( pSTM) || ! pSTM->CreateByRegion( vPL)) + return nullptr ; + // restituisco la superficie + return Release( pSTM) ; +} + +//------------------------------------------------------------------------------- +ISurfTriMesh* +GetSurfTriMeshByExtrusion( const ICurve* pCurve, const Vector3d& vtExtr, bool bCapEnds, double dLinTol) { // verifica parametri - if ( &Curve == nullptr || &vtExtr == nullptr) + if ( pCurve == nullptr || &vtExtr == nullptr) return nullptr ; // calcolo la polilinea che approssima la curva PolyLine PL ; - if ( ! Curve.ApproxWithLines( dLinTol, ANG_TOL_STD_DEG, PL)) + if ( ! pCurve->ApproxWithLines( dLinTol, ANG_TOL_STD_DEG, PL)) return nullptr ; // se richiesta chiusura agli estremi bool bDoCapEnds = false ; @@ -62,6 +86,7 @@ GetSurfTriMeshByExtrusion( const ICurve& Curve, const Vector3d& vtExtr, double dOrthoExtr = plPlane.vtN * vtExtr ; if ( ( fabs( dOrthoExtr) > EPS_SMALL)) { bDoCapEnds = true ; + // se negativa, inverto il senso del contorno if ( dOrthoExtr < 0) PL.Invert() ; } @@ -95,17 +120,75 @@ GetSurfTriMeshByExtrusion( const ICurve& Curve, const Vector3d& vtExtr, //------------------------------------------------------------------------------- ISurfTriMesh* -GetSurfTriMeshByRevolve( const ICurve& Curve, const Point3d& ptAx, const Vector3d& vtAx, +GetSurfTriMeshByRegionExtrusion( const ICURVEPVECTOR& vpCurve, const Vector3d& vtExtr, double dLinTol) +{ + // verifica parametri + if ( &vpCurve == nullptr || vpCurve.empty() || &vtExtr == nullptr) + return nullptr ; + // se una sola curva, uso la funzione precedente + if ( vpCurve.size() == 1 ) + return GetSurfTriMeshByExtrusion( vpCurve[0], vtExtr, true, dLinTol) ; + // calcolo le polilinee che approssimano le curve della regione + POLYLINEVECTOR vPL ; + Vector3d vtN ; + if ( ! CalcRegionPolyLines( vpCurve, dLinTol, vPL, vtN)) + return nullptr ; + // verifico la direzione di estrusione + double dOrthoExtr = vtN * vtExtr ; + if ( ( fabs( dOrthoExtr) < EPS_SMALL)) + return nullptr ; + // se componente estrusione negativa, inverto tutti i percorsi + if ( dOrthoExtr < 0) { + for ( int i = 0 ; i < int( vPL.size()) ; ++ i) + vPL[i].Invert() ; + } + // creo la prima superficie di estremità + PtrOwner pSTM( CreateSurfTriMesh()) ; + if ( IsNull( pSTM) || ! pSTM->CreateByRegion( vPL)) + return nullptr ; + // creo la seconda superficie e la unisco alla prima + { // copio la prima superficie + PtrOwner pSTM2( GetSurfTriMesh( pSTM->Clone())) ; + if ( IsNull( pSTM2)) + return nullptr ; + // inverto la prima superficie + pSTM->Invert() ; + // traslo la seconda + pSTM2->Translate( vtExtr) ; + // la unisco alla prima + if ( ! pSTM->DoSewing( *pSTM2)) + return nullptr ; + } + // creo e unisco le diverse superfici di estrusione + for ( int i = 0 ; i < int( vPL.size()) ; ++ i) { + // estrusione + PtrOwner pSTM2( CreateSurfTriMesh()) ; + if ( IsNull( pSTM2) || ! pSTM2->CreateByExtrusion( vPL[i], vtExtr)) + return nullptr ; + // la unisco alla superficie principale + if ( ! pSTM->DoSewing( *pSTM2)) + return nullptr ; + } + // compatto la superficie + if ( ! pSTM->DoCompacting()) + return nullptr ; + // restituisco la superficie + return Release( pSTM) ; +} + +//------------------------------------------------------------------------------- +ISurfTriMesh* +GetSurfTriMeshByRevolve( const ICurve* pCurve, const Point3d& ptAx, const Vector3d& vtAx, bool bCapEnds, double dLinTol) { // verifica parametri - if ( &Curve == nullptr || &ptAx == nullptr || &vtAx == nullptr) + if ( pCurve == nullptr || &ptAx == nullptr || &vtAx == nullptr) return nullptr ; // limite minimo su tolleranza dLinTol = max( dLinTol, EPS_SMALL) ; // calcolo la polilinea che approssima la curva PolyLine PL ; - if ( ! Curve.ApproxWithLines( dLinTol, ANG_TOL_STD_DEG, PL)) + if ( ! pCurve->ApproxWithLines( dLinTol, ANG_TOL_STD_DEG, PL)) return nullptr ; // calcolo lo step di rotazione double dMaxRad = 0 ; @@ -156,17 +239,17 @@ GetSurfTriMeshByRevolve( const ICurve& Curve, const Point3d& ptAx, const Vector3 //------------------------------------------------------------------------------- ISurfTriMesh* -GetSurfTriMeshByScrewing( const ICurve& Curve, const Point3d& ptAx, const Vector3d& vtAx, +GetSurfTriMeshByScrewing( const ICurve* pCurve, const Point3d& ptAx, const Vector3d& vtAx, double dAngRotDeg, double dMove, double dLinTol) { // verifica parametri - if ( &Curve == nullptr || &ptAx == nullptr || &vtAx == nullptr) + if ( pCurve == nullptr || &ptAx == nullptr || &vtAx == nullptr) return nullptr ; // limite minimo su tolleranza dLinTol = max( dLinTol, EPS_SMALL) ; // calcolo la polilinea che approssima la curva PolyLine PL ; - if ( ! Curve.ApproxWithLines( dLinTol, ANG_TOL_STD_DEG, PL)) + if ( ! pCurve->ApproxWithLines( dLinTol, ANG_TOL_STD_DEG, PL)) return nullptr ; // calcolo lo step di rotazione double dMaxRad = 0 ; @@ -189,18 +272,18 @@ GetSurfTriMeshByScrewing( const ICurve& Curve, const Point3d& ptAx, const Vector //------------------------------------------------------------------------------- ISurfTriMesh* -GetSurfTriMeshRuled( const ICurve& Curve1, const ICurve& Curve2, double dLinTol) +GetSurfTriMeshRuled( const ICurve* pCurve1, const ICurve* pCurve2, double dLinTol) { // verifica parametri - if ( &Curve1 == nullptr || &Curve2 == nullptr) + if ( pCurve1 == nullptr || pCurve2 == nullptr) return nullptr ; // calcolo la polilinea che approssima la prima curva PolyLine PL1 ; - if ( ! Curve1.ApproxWithLines( dLinTol, ANG_TOL_STD_DEG, PL1)) + if ( ! pCurve1->ApproxWithLines( dLinTol, ANG_TOL_STD_DEG, PL1)) return nullptr ; // calcolo la polilinea che approssima la seconda curva PolyLine PL2 ; - if ( ! Curve2.ApproxWithLines( dLinTol, ANG_TOL_STD_DEG, PL2)) + if ( ! pCurve2->ApproxWithLines( dLinTol, ANG_TOL_STD_DEG, PL2)) return nullptr ; // creo e setto la superficie trimesh PtrOwner pSTM( CreateSurfTriMesh()) ; @@ -209,3 +292,62 @@ GetSurfTriMeshRuled( const ICurve& Curve1, const ICurve& Curve2, double dLinTol) // restituisco la superficie return Release( pSTM) ; } + +//------------------------------------------------------------------------------- +bool +CalcRegionPolyLines( const ICURVEPVECTOR& vpCurve, double dLinTol, + POLYLINEVECTOR& vPL, Vector3d& vtN) +{ + // calcolo le polilinee che approssimano le curve + POLYLINEVECTOR vPLtmp ; + vPLtmp.resize( vpCurve.size()) ; + for ( int i = 0 ; i < int( vpCurve.size()) ; ++ i) { + if ( ! vpCurve[i]->ApproxWithLines( dLinTol, ANG_TOL_STD_DEG, vPLtmp[i])) + return nullptr ; + } + // ne calcolo l'area e genero un ordine in senso decrescente + typedef std::pair INDAREA ; // coppia indice, area + typedef std::vector INDAREAVECTOR ; // vettore di coppie indice, area + INDAREAVECTOR vArea ; + vArea.reserve( vPLtmp.size()) ; + Vector3d vtN0 ; + for ( int i = 0 ; i < int( vPLtmp.size()) ; ++ i) { + // verifico chiusura, calcolo piano medio e area + Plane3d plPlane ; + double dArea ; + if ( ! vPLtmp[i].IsClosedAndFlat( plPlane, dArea, 50 * EPS_SMALL)) + return false ; + // imposto la normale del primo contorno come riferimento + if ( i == 0) + vtN0 = plPlane.vtN ; + // verifico che le normali siano molto vicine + if ( ! AreSameOrOppositeVectorApprox( plPlane.vtN, vtN0)) + return false ; + // assegno il segno all'area secondo il verso della normale + if ( ( plPlane.vtN * vtN0) > 0) + vArea.emplace_back( i, dArea) ; + else + vArea.emplace_back( i, - dArea) ; + } + sort( vArea.begin(), vArea.end(), + []( const INDAREA& a, const INDAREA&b) { return fabs( a.second) > fabs( b.second) ; }) ; + // sposto le polilinee nel vettore da restituire secondo l'ordine + vPL.clear() ; + vPL.resize( vPLtmp.size()) ; + bool bCCW = true ; + for ( int i = 0 ; i < int( vPLtmp.size()) ; ++ i) { + // scambio + swap( vPL[i], vPLtmp[vArea[i].first]) ; + // verifico senso di rotazione del contorno esterno + if ( i == 0) + bCCW = ( vArea[i].second > 0) ; + // aggiusto gli altri contorni + else { + if ( ( bCCW && vArea[i].second > 0) || ( ! bCCW && vArea[i].second < 0)) + vPL[i].Invert() ; + } + } + // restituisco la normale positiva alla regione + vtN = ( bCCW ? vtN0 : - vtN0) ; + return true ; +} \ No newline at end of file diff --git a/SurfTriMesh.cpp b/SurfTriMesh.cpp index 42f5885..8a6015e 100644 --- a/SurfTriMesh.cpp +++ b/SurfTriMesh.cpp @@ -75,7 +75,7 @@ SurfTriMesh::AddVertex( const Point3d& ptVert) m_nStatus = TO_VERIFY ; m_OGrMgr.Reset() ; // inserisco il vertice - try { m_vVert.push_back( StmVert( ptVert)) ;} + try { m_vVert.emplace_back( ptVert) ;} catch(...) { return SVT_NULL ;} // ne determino l'indice return int( m_vVert.size() - 1) ; @@ -138,7 +138,7 @@ SurfTriMesh::AddTriangle( const int nIdVert[3]) m_nStatus = TO_VERIFY ; m_OGrMgr.Reset() ; // inserisco il triangolo - try { m_vTria.push_back( StmTria( nIdVert)) ;} + try { m_vTria.emplace_back( nIdVert) ;} catch(...) { return SVT_NULL ;} // ne determino l'indice return int( m_vTria.size() - 1) ; @@ -871,6 +871,43 @@ SurfTriMesh::CreateByFlatContour( const PolyLine& PL) return AdjustTopology() ; } +//---------------------------------------------------------------------------- +bool +SurfTriMesh::CreateByRegion( const POLYLINEVECTOR& vPL) +{ + // eseguo la triangolazione, dopo aver verificato che l'insieme di contorni costituisca una regione + PNTVECTOR vPnt ; + INTVECTOR vTria ; + Triangulate Tri ; + if ( ! Tri.Make( vPL, vPnt, vTria)) + return false ; + + // inizializzo la superficie + int nVert = int( vPnt.size()) ; + int nTria = int( vTria.size()) / 3 ; + if ( ! Init( nVert, nTria)) + return false ; + + // inserisco i vertici nella superficie + for ( int i = 0 ; i < int( vPnt.size()) ; ++ i) { + if ( AddVertex( vPnt[i]) == SVT_NULL) + return false ; + } + + // recupero i triangoli e li inserisco nella superficie + int vV[3] ; + for ( int i = 0 ; i < nTria ; ++i) { + vV[0] = vTria[3*i] ; + vV[1] = vTria[3*i+1] ; + vV[2] = vTria[3*i+2] ; + if ( AddTriangle( vV) == SVT_NULL) + return false ; + } + + // compatto e sistemo la topologia + return AdjustTopology() ; +} + //---------------------------------------------------------------------------- bool SurfTriMesh::CreateByExtrusion( const PolyLine& PL, const Vector3d& vtExtr) diff --git a/SurfTriMesh.h b/SurfTriMesh.h index c8c6129..6a423f0 100644 --- a/SurfTriMesh.h +++ b/SurfTriMesh.h @@ -17,7 +17,6 @@ #include "DllMain.h" #include "GeoObjRW.h" #include "/EgtDev/Include/EGkSurfTriMesh.h" -#include "/EgtDev/Include/EgtNumCollection.h" #include //---------------------------------------------------------------------------- @@ -110,6 +109,7 @@ class SurfTriMesh : public ISurfTriMesh, public IGeoObjRW virtual int AddTriangle( const int nIdVert[3]) ; virtual bool AdjustTopology( void) ; virtual bool CreateByFlatContour( const PolyLine& PL) ; + virtual bool CreateByRegion( const POLYLINEVECTOR& vPL) ; virtual bool CreateByExtrusion( const PolyLine& PL, const Vector3d& vtExtr) ; virtual bool CreateByTwoCurves( const PolyLine& PL1, const PolyLine& PL2) ; virtual bool CreateByRevolution( const PolyLine& PL, const Point3d& ptAx, const Vector3d& vtAx, diff --git a/Triangulate.cpp b/Triangulate.cpp index db0e986..e60c70f 100644 --- a/Triangulate.cpp +++ b/Triangulate.cpp @@ -13,12 +13,24 @@ //--------------------------- Include ---------------------------------------- #include "stdafx.h" +#include #include "Triangulate.h" #include "\EgtDev\Include\EGkPolyLine.h" #include "\EgtDev\Include\EGkPlane3d.h" +#include "DllMain.h" +#include "/EgtDev/Include/EGkStringUtils3d.h" +#include "/EgtDev/Include/EgtLogger.h" + using namespace std ; +//---------------------------------------------------------------------------- +static enum PlaneType { PL_XY = 1, PL_YZ = 2, PL_ZX = 3} ; +static enum EarStatus{ EAS_NULL = -1, EAS_NO = 0, EAS_OK = 1} ; + +//---------------------------------------------------------------------------- +static bool ChangeStartPntVector( int nNewStart, PNTVECTOR& vPi) ; + //---------------------------------------------------------------------------- // In : PolyLine // Out : PNTVECTOR (Point3d Vector) : points of the polyline @@ -56,20 +68,15 @@ Triangulate::Make( const PolyLine& PL, PNTVECTOR& vPt, INTVECTOR& vTr) // inserisco i punti while ( PL.GetNextPoint( ptP)) vPt.push_back( ptP) ; + // se non CCW inverto il vettore dei punti + if ( ! bCCW) + reverse( vPt.begin(), vPt.end()) ; // 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) { - for ( int i = 0 ; i < n ; ++ i) - vPol.push_back( i) ; - } - // altrimenti lo prendo al contrario - else { - for ( int i = n - 1 ; i >= 0 ; -- i) - vPol.push_back( i) ; - } + for ( int i = 0 ; i < n ; ++ i) + vPol.push_back( i) ; // eseguo la triangolazione if ( ! MakeByEC2( vPt, vPol, vTr)) @@ -82,6 +89,153 @@ Triangulate::Make( const PolyLine& PL, PNTVECTOR& vPt, INTVECTOR& vTr) return true ; } +//---------------------------------------------------------------------------- +// In : POLYLINEVECTOR : vector of polylines, the first outer, the others inner +// Out : PNTVECTOR (Point3d Vector) : points of the polyline +// INTVECTOR (int Vector) : 3*T indices of above points for T triangles +//---------------------------------------------------------------------------- +bool +Triangulate::Make( const POLYLINEVECTOR& vPL, PNTVECTOR& vPt, INTVECTOR& vTr) +{ + // pulisco i vettori di ritorno + vPt.clear() ; + vTr.clear() ; + // ci devono essere delle polilinee nel vettore + if ( &vPL == nullptr || vPL.empty()) + return false ; + // se una sola polilinea mi riconduco al caso precedente + if ( vPL.size() == 1) + return Make( vPL[0], vPt, vTr) ; + // verifico che la polilinea esterna sia chiusa e piana e calcolo il piano medio del poligono + double dArea ; + Plane3d plExtPlane ; + if ( ! vPL[0].IsClosedAndFlat( plExtPlane, dArea, 50 * EPS_SMALL)) + return false ; + bool bCCW ; + if ( fabs( plExtPlane.vtN.z) >= fabs( plExtPlane.vtN.x) && + fabs( plExtPlane.vtN.z) >= fabs( plExtPlane.vtN.y)) { + m_nPlane = PL_XY ; + bCCW = ( plExtPlane.vtN.z > 0) ; + } + else if ( fabs( plExtPlane.vtN.x) >= fabs( plExtPlane.vtN.y)) { + m_nPlane = PL_YZ ; + bCCW = ( plExtPlane.vtN.x > 0) ; + } + else { + m_nPlane = PL_ZX ; + bCCW = ( plExtPlane.vtN.y > 0) ; + } + // verifico le altre polilinee + for ( int i = 1 ; i < int( vPL.size()) ; ++i) { + // deve essere chiusa, giacere nello stesso piano ed essere orientata al contrario della esterna + double dArea ; + Plane3d plPlane ; + if ( ! vPL[i].IsClosedAndFlat( plPlane, dArea, 50 * EPS_SMALL) || + ! AreOppositeVectorApprox( plExtPlane.vtN, plPlane.vtN)) + return false ; + } + // se non CCW inverto tutte le polilinee + if ( ! bCCW) { + for ( int i = 0 ; i < int( vPL.size()) ; ++i) + const_cast(vPL[i]).Invert( true) ; + } + // calcolo ordine decrescente su Xmax o Ymax o Zmax secondo m_nPlane per le curve interne + INTVECTOR vOrd ; + if ( ! SortInternalLoops( vPL, vOrd)) + return false ; + // riempio il vettore con i punti dei poligoni da triangolare + // calcolo spazio totale e spazio massimo per un percorso interno + int nTot = vPL[0].GetPointNbr() - 1 ; + int nMax = 0 ; + for ( int i = 1 ; i < int( vPL.size()) ; ++i) { + nTot += vPL[i].GetPointNbr() + 1 ; + if ( vPL[i].GetPointNbr() > nMax) + nMax = vPL[i].GetPointNbr() ; + } + // riservo spazio pari al totale dei punti (anche quelli ripetuti di chiusura che servono per i link) + vPt.reserve( nTot) ; + // inserisco i punti del contorno esterno (tranne il primo che coincide con l'ultimo) + if ( ! GetPntVectorFromPolyline( vPL[0], false, vPt)) + return false ; + // ciclo sui percorsi interni ordinati secondo Xmax decrescente + PNTVECTOR vPi ; + vPi.reserve( nMax) ; + for ( int i = 1 ; i < int( vPL.size()) ; ++ i) { + // riordino il percorso per avere punto iniziale con Xmax + if ( ! GetPntVectorFromPolyline( vPL[vOrd[i]], true, vPi)) + return false ; + // cerco un punto del percorso esterno visibile dal punto iniziale del precedente interno + int nI ; + if ( ! GetOuterPntToJoin( vPt, vPi[0], nI)) + return false ; + // riordino percorso esterno per avere questo punto all'inizio + if ( ! ChangeStartPntVector( nI, vPt)) + return false ; + // ripeto punto iniziale + vPt.push_back( vPt[0]) ; + // accodo percorso interno + for ( int j = 0 ; j < int( vPi.size()) ; ++j) + vPt.push_back( vPi[j]) ; + // ripeto punto iniziale di percorso interno + vPt.push_back( vPi[0]) ; + // cancello percorso interno + vPi.clear() ; + } + // ripristino l'orientamento delle polilinee originali per il caso non CCW + if ( ! bCCW) { + for ( int i = 0 ; i < int( vPL.size()) ; ++i) + const_cast(vPL[i]).Invert( true) ; + } + + // creo il vettore degli indici del Poligono + INTVECTOR vPol ; + int n = int( vPt.size()) ; + vPol.reserve( n) ; + // non devo gestire separatamente CCW perchè ho già invertito i punti + for ( int i = 0 ; i < n ; ++ i) + vPol.push_back( i) ; + + // eseguo la triangolazione + if ( ! MakeByEC2( vPt, vPol, vTr)) + return false ; + + // se era CW, devo invertire il senso dei triangoli + if ( ! bCCW) + reverse( vTr.begin(), vTr.end()) ; + + return true ; +} + +//---------------------------------------------------------------------------- +bool +Triangulate::PrepareGrid( const PNTVECTOR& vPt, const INTVECTOR& vPol, + const INTVECTOR& vPrev, const INTVECTOR& vNext) +{ + // points number + int n = int( vPol.size()) ; + // overall box + BBox3d b3All ; + for ( int j = 0 ; j < n ; ++ j) + b3All.Add( vPt[vPol[j]]) ; + // grid cell dimension + double dCellDim ; + b3All.GetRadius( dCellDim) ; + dCellDim *= 2 / sqrt( n) ; + // grid + if ( ! m_VertGrid.Init( 2 * n, dCellDim)) + return false ; + for ( int j = 0 ; j < n ; ++ j) { + // insert only reflex vertex + if ( ! TriangleIsCCW( vPt[vPol[vPrev[j]]], vPt[vPol[j]], vPt[vPol[vNext[j]]])) { + if ( ! m_VertGrid.InsertPoint( vPt[vPol[j]], j)) + return false ; + } + } + m_vVert.reserve( n / 5) ; + + return true ; +} + //---------------------------------------------------------------------------- // Triangulate the CCW n-gon specified by the vertices vPt (Pt[n] != Pt[0]) // Ear Clipping algorithm @@ -114,19 +268,27 @@ Triangulate::MakeByEC( const PNTVECTOR& vPt, const INTVECTOR& vPol, INTVECTOR& v 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 = TestTriangle( vPt, vPol, vPrev, vNext, i) ; + bool bSave = bIsEar ; + // Verify if triangle is null (accept to discard) + if ( ! bIsEar) { + if ( Collinear( vPt[vPol[vPrev[i]]], vPt[vPol[i]], vPt[vPol[vNext[i]]])) + bIsEar = true ; + } // 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]]) ; + if ( bSave) { + 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] ; @@ -142,17 +304,13 @@ Triangulate::MakeByEC( const PNTVECTOR& vPt, const INTVECTOR& vPol, INTVECTOR& v 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 +// Ear Clipping algorithm enhanced, choose smaller diagonal //---------------------------------------------------------------------------- bool Triangulate::MakeByEC2( const PNTVECTOR& vPt, const INTVECTOR& vPol, INTVECTOR& vTr) @@ -178,39 +336,53 @@ Triangulate::MakeByEC2( const PNTVECTOR& vPt, const INTVECTOR& vPol, INTVECTOR& vPrev[0] = n - 1 ; vNext[n-1] = 0 ; + // Prepare Ear status vector + INTVECTOR vEar( n, EAS_NULL) ; + + // Prepare PointGrid + if ( ! PrepareGrid( vPt, vPol, vPrev, vNext)) + return false ; + // 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 = TestTriangle( vPt, vPol, vPrev, vNext, i) ; + if ( vEar[i] == EAS_NULL) + vEar[i] = ( TestTriangle( vPt, vPol, vPrev, vNext, i) ? EAS_OK : EAS_NO) ; + bool bIsEar = ( vEar[i] == EAS_OK) ; + bool bSave = bIsEar ; if ( bIsEar) { // Save square distance of diagonal double dSqDist = SqDist(vPt[vPol[vPrev[i]]], vPt[vPol[vNext[i]]]) ; - // Try with next or next^2 - int j = vNext[i] ; + // Try with 3 next + int j = i ; double dSqDist1 = INFINITO ; - if ( TestTriangle( vPt, vPol, vPrev, vNext, j)) - dSqDist1 = SqDist( vPt[vPol[vPrev[j]]], vPt[vPol[vNext[j]]]) ; - else { + for ( int h = 0 ; h < 3 ; ++ h) { j = vNext[j] ; - if ( TestTriangle( vPt, vPol, vPrev, vNext, j)) + if ( vEar[j] == EAS_NULL) + vEar[j] = ( TestTriangle( vPt, vPol, vPrev, vNext, j) ? EAS_OK : EAS_NO) ; + if ( vEar[j] == EAS_OK) { dSqDist1 = SqDist( vPt[vPol[vPrev[j]]], vPt[vPol[vNext[j]]]) ; + break ; + } } - // Try with prev or prev^2 - int k = vPrev[i] ; + // Try with 3 prev + int k = i ; double dSqDist2 = INFINITO ; - if ( TestTriangle( vPt, vPol, vPrev, vNext, k)) - dSqDist2 = SqDist( vPt[vPol[vPrev[k]]], vPt[vPol[vNext[k]]]) ; - else { + for ( int h = 0 ; h < 3 ; ++ h) { k = vPrev[k] ; - if ( TestTriangle( vPt, vPol, vPrev, vNext, k)) + if ( vEar[k] == EAS_NULL) + vEar[k] = ( TestTriangle( vPt, vPol, vPrev, vNext, k) ? EAS_OK : EAS_NO) ; + if ( vEar[k] == EAS_OK) { dSqDist2 = SqDist( vPt[vPol[vPrev[k]]], vPt[vPol[vNext[k]]]) ; + break ; + } } // Choose the better (EPS_ZERO to have a difference) if ( dSqDist < dSqDist1 + EPS_ZERO && dSqDist < dSqDist2 + EPS_ZERO) @@ -220,32 +392,151 @@ Triangulate::MakeByEC2( const PNTVECTOR& vPt, const INTVECTOR& vPol, INTVECTOR& else i = k ; } - + // Verify if triangle is null (accept to discard) + else { + if ( Collinear( vPt[vPol[vPrev[i]]], vPt[vPol[i]], vPt[vPol[vNext[i]]]) && + ! Aligned( vPt[vPol[vPrev[i]]], vPt[vPol[i]], vPt[vPol[vNext[i]]])) + bIsEar = true ; + } // 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 ; + // Triangle (v[prev[i]], v[i], v[next[i]]) is an ear + if ( bSave) { + vTr.push_back( vPol[vPrev[i]]) ; + vTr.push_back( vPol[i]) ; + vTr.push_back( vPol[vNext[i]]) ; + } + // Reset earity of diagonal endpoints + vEar[vPrev[i]] = EAS_NULL ; + vEar[vNext[i]] = EAS_NULL ; + // ‘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 ; +} + +//---------------------------------------------------------------------------- +// 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 +//---------------------------------------------------------------------------- +bool +Triangulate::MakeByEC3( 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 ; + + // Prepare Ear status vector + INTVECTOR vEar( n, EAS_NULL) ; + + // Prepare PointGrid + if ( ! PrepareGrid( vPt, vPol, vPrev, vNext)) + return false ; + + // 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 + if ( vEar[i] == EAS_NULL) + vEar[i] = ( TestTriangle( vPt, vPol, vPrev, vNext, i) ? EAS_OK : EAS_NO) ; + bool bIsEar = ( vEar[i] == EAS_OK) ; + bool bSave = bIsEar ; + if ( bIsEar) { + // Square Length & Aspect Ratio + 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 verticis + 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) { + 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) ; + i = j ; + } + } + } + j = vNext[j] ; + } + } + // Verify if triangle is null (accept to discard) + else { + if ( Collinear( vPt[vPol[vPrev[i]]], vPt[vPol[i]], vPt[vPol[vNext[i]]]) && + ! Aligned( vPt[vPol[vPrev[i]]], vPt[vPol[i]], vPt[vPol[vNext[i]]])) + bIsEar = true ; + } + // 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 to save + if ( bSave) { + vTr.push_back( vPol[vPrev[i]]) ; + vTr.push_back( vPol[i]) ; + vTr.push_back( vPol[vNext[i]]) ; + } + // Reset earity of diagonal endpoints + vEar[vPrev[i]] = EAS_NULL ; + vEar[vNext[i]] = EAS_NULL ; + // ‘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] ; + } + } return true ; } @@ -260,15 +551,36 @@ Triangulate::TestTriangle( const PNTVECTOR& vPt, const INTVECTOR& vPol, // 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]) ; + BBox3d b3Tria ; + b3Tria.Add( vPt[vPol[vPrev[i]]]) ; + b3Tria.Add( vPt[vPol[i]]) ; + b3Tria.Add( vPt[vPol[vNext[i]]]) ; + m_VertGrid.Find( b3Tria, m_vVert) ; + for ( int j = 0 ; j < int( m_vVert.size()) ; ++ j) { + int k = m_vVert[j] ; + if ( k == i || k == vPrev[i] || k == vNext[i] ) + continue ; + // If vertex k is on a triangle vertex + if ( AreSamePoint( vPt[vPol[k]], vPt[vPol[vPrev[i]]]) || + AreSamePoint( vPt[vPol[k]], vPt[vPol[i]]) || + AreSamePoint( vPt[vPol[k]], vPt[vPol[vNext[i]]])) { + // Test previous segment with triangle diagonal + if ( TestIntersection( vPt[vPol[vPrev[k]]], vPt[vPol[k]], vPt[vPol[vPrev[i]]], vPt[vPol[vNext[i]]])) { + bIsEar = false ; + break ; + } + // Test next segment with triangle diagonal + if ( TestIntersection( vPt[vPol[k]], vPt[vPol[vNext[k]]], vPt[vPol[vPrev[i]]], vPt[vPol[vNext[i]]])) { + bIsEar = false ; + break ; + } + } + // If vertex k is inside the ear triangle, then this is not an ear + else if ( TestPointInTriangle( vPt[vPol[k]], vPt[vPol[vPrev[i]]], vPt[vPol[i]], vPt[vPol[vNext[i]]])) { + bIsEar = false ; + break ; + } + } } else { // The ‘ear’ triangle is clockwise so v[i] is not an ear @@ -279,36 +591,103 @@ Triangulate::TestTriangle( const PNTVECTOR& vPt, const INTVECTOR& vPol, } //---------------------------------------------------------------------------- +// Ratio between max side and opposite height +double +Triangulate::CalcTriangleAspectRatio( const Point3d& ptPa, const Point3d& ptPb, const Point3d& ptPc) +{ + double dSqDistA = SquareDist( ptPa, ptPb) ; + double dSqDistB = SquareDist( ptPb, ptPc) ; + double dSqDistC = SquareDist( ptPc, ptPa) ; + double dTwoArea = fabs( TwoArea( ptPa, ptPb, ptPc)) ; + if ( dTwoArea < EPS_SMALL * EPS_SMALL) + return INFINITO ; + else + return ( std::max( dSqDistA, std::max( dSqDistB, dSqDistC)) / dTwoArea) ; +} + +//---------------------------------------------------------------------------- +double +Triangulate::SquareDist( const Point3d& ptA, const Point3d& ptB) +{ + switch ( m_nPlane) { + default : // PL_XY + return (( ptB.x - ptA.x) * ( ptB.x - ptA.x) + ( ptB.y - ptA.y) * ( ptB.y - ptA.y)) ; + case PL_YZ : + return (( ptB.y - ptA.y) * ( ptB.y - ptA.y) + ( ptB.z - ptA.z) * ( ptB.z - ptA.z)) ; + case PL_ZX : + return (( ptB.z - ptA.z) * ( ptB.z - ptA.z) + ( ptB.x - ptA.x) * ( ptB.x - ptA.x)) ; + } +} + +//---------------------------------------------------------------------------- +// Vector product is double of the area (positive if CCW) +double +Triangulate::TwoArea( const Point3d& ptA, const Point3d& ptB, const Point3d& ptC) +{ + switch ( m_nPlane) { + default : // PL_XY + return ( ptB.x - ptA.x) * ( ptC.y - ptB.y) - ( ptB.y - ptA.y) * ( ptC.x - ptB.x) ; + case PL_YZ : + return ( ptB.y - ptA.y) * ( ptC.z - ptB.z) - ( ptB.z - ptA.z) * ( ptC.y - ptB.y) ; + case PL_ZX : + return ( ptB.z - ptA.z) * ( ptC.x - ptB.x) - ( ptB.x - ptA.x) * ( ptC.z - ptB.z) ; + } +} + +//---------------------------------------------------------------------------- +// Distance less than tolerance +bool +Triangulate::AreSamePoint( const Point3d& ptA, const Point3d& ptB, double dToler) +{ + return ( SquareDist( ptA, ptB) < dToler * dToler) ; +} + +//---------------------------------------------------------------------------- +// Collinear <--> Null Area +bool +Triangulate::Aligned( const Point3d& ptA, const Point3d& ptB, const Point3d& ptC) +{ + switch ( m_nPlane) { + default : // PL_XY + return (( ptB.x - ptA.x) * ( ptC.x - ptB.x) + ( ptB.y - ptA.y) * ( ptC.y - ptB.y)) > 0 ; + case PL_YZ : + return (( ptB.y - ptA.y) * ( ptC.y - ptB.y) + ( ptB.z - ptA.z) * ( ptC.z - ptB.z)) > 0 ; + case PL_ZX : + return (( ptB.z - ptA.z) * ( ptC.z - ptB.z) + ( ptB.x - ptA.x) * ( ptC.x - ptB.x)) > 0 ; + } +} + +//---------------------------------------------------------------------------- +// Collinear <--> Null Area +bool +Triangulate::Collinear( const Point3d& ptA, const Point3d& ptB, const Point3d& ptC, double dToler) +{ + return ( fabs( TwoArea( ptA, ptB, ptC)) < dToler * dToler) ; +} + +//---------------------------------------------------------------------------- +// Positive Area means A -> B -> C counterclockwise bool Triangulate::TriangleIsCCW( const Point3d& ptA, const Point3d& ptB, const Point3d& ptC, double dToler) { - double dV11 ; - double dV12 ; - double dV21 ; - double dV22 ; + return ( TwoArea( ptA, ptB, ptC) > dToler * dToler) ; +} - switch ( m_nPlane) { - default : // PL_XY - dV11 = ptB.x - ptA.x ; - dV12 = ptB.y - ptA.y ; - dV21 = ptC.x - ptB.x ; - dV22 = ptC.y - ptB.y ; - break ; - case PL_YZ : - dV11 = ptB.y - ptA.y ; - dV12 = ptB.z - ptA.z ; - dV21 = ptC.y - ptB.y ; - dV22 = ptC.z - ptB.z ; - break ; - case PL_ZX : - dV11 = ptB.z - ptA.z ; - dV12 = ptB.x - ptA.x ; - dV21 = ptC.z - ptB.z ; - dV22 = ptC.x - ptB.x ; - break ; - } - - return ( ( dV11 * dV22 - dV12 * dV21) > dToler * dToler) ; +//---------------------------------------------------------------------------- +// Proper intersection between line segments +bool +Triangulate::TestIntersection( const Point3d& ptA1, const Point3d& ptA2, + const Point3d& ptB1, const Point3d& ptB2) +{ + // Collinearity is not considered intersection + if ( Collinear( ptA1, ptA2, ptB1) || + Collinear( ptA1, ptA2, ptB2) || + Collinear( ptB1, ptB2, ptA1) || + Collinear( ptB1, ptB2, ptA2)) + return false ; + // To intersect, the endpoints of a line segment must be on opposite side of the other line segment + return ( TriangleIsCCW( ptA1, ptA2, ptB1) != TriangleIsCCW( ptA1, ptA2, ptB2) && + TriangleIsCCW( ptB1, ptB2, ptA1) != TriangleIsCCW( ptB1, ptB2, ptA2)) ; } //---------------------------------------------------------------------------- @@ -316,20 +695,279 @@ 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 P is on a vertex is considered outside - if ( AreSamePointApprox( ptP, ptA)) + // If P is on a vertex is considered outside + if ( AreSamePoint( ptP, ptA)) return false ; - if ( AreSamePointApprox( ptP, ptB)) + if ( AreSamePoint( ptP, ptB)) return false ; - if ( AreSamePointApprox( ptP, ptC)) + if ( AreSamePoint( ptP, ptC)) return false ; - // if P is on the right of at least one edge is outside + // 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 ; + // P is in triangle return true ; } +//---------------------------------------------------------------------------- +bool +Triangulate::SortInternalLoops( const POLYLINEVECTOR& vPL, INTVECTOR& vOrd) +{ + // riempio vettore con indice e massimo specifico per ogni loop interno + typedef std::pair INDMAX ; // coppia indice, massimo + typedef std::vector INDMAXVECTOR ; // vettore di coppie indice, massimo + INDMAXVECTOR vMax ; + vMax.reserve( vPL.size()) ; + vMax.emplace_back( 0, 0.0) ; + for ( int i = 1 ; i < int( vPL.size()) ; ++ i) { + BBox3d b3Loc ; + vPL[i].GetLocalBBox( b3Loc) ; + switch ( m_nPlane) { + default : // PL_XY + vMax.emplace_back( i, b3Loc.GetMax().x) ; + break ; + case PL_YZ : + vMax.emplace_back( i, b3Loc.GetMax().y) ; + break ; + case PL_ZX : + vMax.emplace_back( i, b3Loc.GetMax().z) ; + break ; + } + } + // ordino vettore in senso decrescente rispetto al massimo + sort( vMax.begin() + 1, vMax.end(), + []( const INDMAX& a, const INDMAX&b) { return a.second > b.second ; }) ; + // copio indice nel vettore di ordine + vOrd.reserve( vPL.size()) ; + for ( int i = 0 ; i < int( vPL.size()) ; ++ i) + vOrd.push_back( vMax[i].first) ; + + return true ; +} + +//---------------------------------------------------------------------------- +bool +Triangulate::GetPntVectorFromPolyline( const PolyLine& PL, bool bXmaxStart, PNTVECTOR& vPi) +{ + // copio i punti nel vettore (tranne il primo che coincide con l'ultimo) + Point3d ptP ; + if ( ! PL.GetFirstPoint( ptP)) + return false ; + while ( PL.GetNextPoint( ptP)) { + vPi.push_back( ptP) ; + } + // se non richiesto riordino per Xmax o Ymax o Zmax, ho finito + if ( ! bXmaxStart) + return true ; + // determino l'indice del punto con Xmax + int nI = - 1 ; + double dXmax = - INFINITO ; + for ( int i = 0 ; i < int( vPi.size()) ; ++ i) { + switch ( m_nPlane) { + default : // PL_XY + if ( vPi[i].x > dXmax) { + dXmax = vPi[i].x ; + nI = i ; + } + break ; + case PL_YZ : + if ( vPi[i].y > dXmax) { + dXmax = vPi[i].y ; + nI = i ; + } + break ; + case PL_ZX : + if ( vPi[i].z > dXmax) { + dXmax = vPi[i].z ; + nI = i ; + } + break ; + } + } + if ( nI == - 1) + return false ; + // riordino il vettore, per avere il punto con Xmax in prima posizione + return ChangeStartPntVector( nI, vPi) ; +} + +//---------------------------------------------------------------------------- +bool +Triangulate::GetOuterPntToJoin( const PNTVECTOR& vPt, const Point3d& ptP, int& nI) +{ + // cerco prima intersezione del raggio dal punto con direzione X+ al contorno esterno + nI = - 1 ; + double dMinDist = INFINITO ; + Point3d ptInt ; + int nNumPt = int( vPt.size()) ; + for ( int i = 0 ; i < nNumPt ; ++ i) { + // indice punto precedente + int h = ( i - 1 >= 0 ) ? i - 1 : nNumPt - 1 ; + // indice punto successivo + int j = ( i + 1 < nNumPt) ? i + 1 : 0 ; + // mi metto nel piano principale + switch (m_nPlane) { + default : // PL_XY + // se punto esattamente sul raggio e raggio interno al settore + if ( vPt[i].x > ptP.x && fabs( vPt[i].y - ptP.y) < EPS_SMALL && + PointInSector( ptP, vPt[h], vPt[i], vPt[j])) { + // se distanza minore della minima, nuovo minimo + if ( ( vPt[i].x - ptP.x) < dMinDist) { + dMinDist = vPt[i].x - ptP.x ; + nI = i ; + ptInt = vPt[i] ; + } + } + // se segmento al punto successivo che attraversa il raggio e raggio interno ovvero a sinistra ( segmento crescente in Y) + else if ( vPt[i].y < ptP.y && vPt[j].y > ptP.y) { + // calcolo l'ascissa di intersezione + double dCoeff = ( ptP.y - vPt[i].y) / ( vPt[j].y - vPt[i].y) ; + double dX = vPt[i].x + ( vPt[j].x - vPt[i].x) * dCoeff ; + // se sta sul raggio e distanza minore della minima + if ( dX > ptP.x && ( dX - ptP.x) < dMinDist) { + dMinDist = dX - ptP.x ; + nI = ( vPt[i].x >= vPt[j].x) ? i : j ; + double dZ = vPt[i].z + ( vPt[j].z - vPt[i].z) * dCoeff ; + ptInt.Set( dX, ptP.y, dZ) ; + } + } + break ; + case PL_YZ : + // se punto esattamente sul raggio e raggio interno al settore + if ( vPt[i].y > ptP.y && fabs( vPt[i].z - ptP.z) < EPS_SMALL && + PointInSector( ptP, vPt[h], vPt[i], vPt[j])) { + // se distanza minore della minima, nuovo minimo + if ( ( vPt[i].y - ptP.y) < dMinDist) { + dMinDist = vPt[i].y - ptP.y ; + nI = i ; + ptInt = vPt[i] ; + } + } + // se segmento al punto successivo che attraversa il raggio e raggio interno ovvero a sinistra ( segmento crescente in Y) + else if ( vPt[i].z < ptP.z && vPt[j].z > ptP.z) { + // calcolo l'ascissa di intersezione + double dCoeff = ( ptP.z - vPt[i].z) / ( vPt[j].z - vPt[i].z) ; + double dY = vPt[i].y + ( vPt[j].y - vPt[i].y) * dCoeff ; + // se sta sul raggio e distanza minore della minima + if ( dY > ptP.y && ( dY - ptP.y) < dMinDist) { + dMinDist = dY - ptP.y ; + nI = ( vPt[i].y >= vPt[j].y) ? i : j ; + double dX = vPt[i].x + ( vPt[j].x - vPt[i].x) * dCoeff ; + ptInt.Set( dX, dY, ptP.z) ; + } + } + break ; + case PL_ZX : + // se punto esattamente sul raggio e raggio interno al settore + if ( vPt[i].z > ptP.z && fabs( vPt[i].x - ptP.x) < EPS_SMALL && + PointInSector( ptP, vPt[h], vPt[i], vPt[j])) { + // se distanza minore della minima, nuovo minimo + if ( ( vPt[i].z - ptP.z) < dMinDist) { + dMinDist = vPt[i].z - ptP.z ; + nI = i ; + ptInt = vPt[i] ; + } + } + // se segmento al punto successivo che attraversa il raggio e raggio interno ovvero a sinistra ( segmento crescente in Y) + else if ( vPt[i].x < ptP.x && vPt[j].x > ptP.x) { + // calcolo l'ascissa di intersezione + double dCoeff = ( ptP.x - vPt[i].x) / ( vPt[j].x - vPt[i].x) ; + double dZ = vPt[i].z + ( vPt[j].z - vPt[i].z) * dCoeff ; + // se sta sul raggio e distanza minore della minima + if ( dZ > ptP.z && ( dZ - ptP.z) < dMinDist) { + dMinDist = dZ - ptP.z ; + nI = ( vPt[i].z >= vPt[j].z) ? i : j ; + double dY = vPt[i].y + ( vPt[j].y - vPt[i].y) * dCoeff ; + ptInt.Set( ptP.y, dY, dZ) ; + } + } + break ; + } + } + // non ho trovato alcunché, errore + if ( nI == - 1) + return false ; + // se ho trovato un punto esatto del contorno, non devo fare altri controlli + if ( AreSamePointApprox( ptInt, vPt[nI])) + return true ; + // devo ora verificare che il segmento che unisce i punti non intersechi altri lati del contorno esterno + // altrimenti tengo il punto con raggio più vicino a X_AX o Y_AX o Z_AX secondo m_nPlane + int nJ = nI ; + Point3d ptPa = ptP ; + Point3d ptPb = vPt[nI] ; + Point3d ptPc = ptInt ; + bool bSwap = false ; + switch ( m_nPlane) { + default : /* PL_XY */ bSwap = ( ptPb.y > ptPc.y) ; break ; + case PL_YZ : bSwap = ( ptPb.z > ptPc.z) ; break ; + case PL_ZX : bSwap = ( ptPb.x > ptPc.x) ; break ; + } + if ( bSwap) + swap( ptPb, ptPc) ; + double dMinTan = INFINITO ; + double dMinSqDist = INFINITO * INFINITO ; + for ( int i = 0 ; i < nNumPt ; ++ i) { + // salto il punto già trovato + if ( i == nJ) + continue ; + // verifico se sta nel triangolo + if ( TestPointInTriangle( vPt[i], ptPa, ptPb, ptPc)) { + double dTan = INFINITO ; + switch ( m_nPlane) { + default : /* PL_XY */ dTan = fabs(vPt[i].y - ptP.y) / (vPt[i].x - ptP.x) ; break ; + case PL_YZ : dTan = fabs(vPt[i].z - ptP.z) / (vPt[i].y - ptP.y) ; break ; + case PL_ZX : dTan = fabs(vPt[i].x - ptP.x) / (vPt[i].z - ptP.z) ; break ; + } + if ( dTan < dMinTan + EPS_ZERO) { + // indice punto precedente + int h = ( i - 1 >= 0) ? i - 1 : nNumPt - 1 ; + // indice punto successivo + int j = ( i + 1 < nNumPt) ? i + 1 : 0 ; + // verifico che il raggio che unisce i punti stia nel settore + if ( PointInSector( ptP, vPt[h], vPt[i], vPt[j])) { + double dSqDist = SqDist( vPt[i], ptP) ; + if ( dTan < dMinTan - EPS_ZERO || dSqDist < dMinSqDist) { + dMinTan = dTan ; + dMinSqDist = dSqDist ; + nI = i ; + } + } + } + } + } + + return true ; +} + +//---------------------------------------------------------------------------- +bool +Triangulate::PointInSector( const Point3d& ptTest, const Point3d& ptPrev, const Point3d& ptCorn, const Point3d& ptNext) +{ + // la parte valida del settore è a sinistra dei segmenti ptPrev --> ptCorn --> ptNext + // se corner convesso + if ( TriangleIsCCW( ptPrev, ptCorn, ptNext, 0)) + return ( TriangleIsCCW( ptPrev, ptCorn, ptTest) && + TriangleIsCCW( ptTest, ptCorn, ptNext)) ; + // altrimenti corner concavo ( reflex) + else + return ! ( TriangleIsCCW( ptTest, ptCorn, ptPrev, 0) && + TriangleIsCCW( ptNext, ptCorn, ptTest, 0)) ; +} + +//---------------------------------------------------------------------------- +bool +ChangeStartPntVector( int nNewStart, PNTVECTOR& vPi) +{ + // se il nuovo inizio coincide col vecchio, non devo fare alcunché + if ( nNewStart == 0) + return true ; + // se il nuovo indice è oltre la dimensione del vettore, errore + if ( nNewStart >= int( vPi.size())) + return false ; + // ciclo di aggiustamento + rotate( vPi.begin(), vPi.begin() + nNewStart, vPi.end()) ; + return true ; +} diff --git a/Triangulate.h b/Triangulate.h index 01cb650..a921f23 100644 --- a/Triangulate.h +++ b/Triangulate.h @@ -12,27 +12,41 @@ //---------------------------------------------------------------------------- #pragma once -#include "/EgtDev/Include/EGkGeoCollection.h" -class PolyLine ; - -//---------------------------------------------------------------------------- -enum { PL_XY = 1, PL_YZ = 2, PL_ZX = 3} ; +#include "/EgtDev/Include/EGkPolyLine.h" +#include "\EgtDev\Include\EGkPointGrid3d.h" //---------------------------------------------------------------------------- class Triangulate { public : bool Make( const PolyLine& PL, PNTVECTOR& vPt, INTVECTOR& vTr) ; + bool Make( const POLYLINEVECTOR& vPL, PNTVECTOR& vPt, INTVECTOR& vTr) ; private : + bool PrepareGrid( const PNTVECTOR& vPt, const INTVECTOR& vPol, + const INTVECTOR& vPrev, const INTVECTOR& vNext) ; bool MakeByEC( const PNTVECTOR& vPt, const INTVECTOR& vPol, INTVECTOR& vTr) ; bool MakeByEC2( const PNTVECTOR& vPt, const INTVECTOR& vPol, INTVECTOR& vTr) ; + bool MakeByEC3( const PNTVECTOR& vPt, const INTVECTOR& vPol, 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) ; + 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) ; + bool Aligned( const Point3d& ptA, const Point3d& ptB, const Point3d& ptC) ; + bool Collinear( const Point3d& ptA, const Point3d& ptB, const Point3d& ptC, double dToler = EPS_SMALL) ; bool TriangleIsCCW( const Point3d& ptA, const Point3d& ptB, const Point3d& ptC, double dToler = EPS_SMALL) ; + bool TestIntersection( const Point3d& ptA1, const Point3d& ptA2, const Point3d& ptB1, const Point3d& ptB2) ; bool TestPointInTriangle( const Point3d& ptP, const Point3d& ptA, const Point3d& ptB, const Point3d& ptC) ; + bool SortInternalLoops( const POLYLINEVECTOR& vPL, INTVECTOR& vOrd) ; + bool GetPntVectorFromPolyline( const PolyLine& PL, bool bXmaxStart, PNTVECTOR& vPi) ; + bool GetOuterPntToJoin( const PNTVECTOR& vPt, const Point3d& ptP, int& nI) ; + bool PointInSector( const Point3d& ptTest, const Point3d& ptPrev, const Point3d& ptCorn, const Point3d& ptNext) ; private : int m_nPlane ; -} ; \ No newline at end of file + PointGrid3d m_VertGrid ; + INTVECTOR m_vVert ; +} ;