diff --git a/IntersPlanePlane.cpp b/IntersPlanePlane.cpp index a4259a5..0b83506 100644 --- a/IntersPlanePlane.cpp +++ b/IntersPlanePlane.cpp @@ -42,5 +42,5 @@ Inters3Planes( const Plane3d& plPlane1, const Plane3d& plPlane2, const Plane3d& if ( IntersPlanePlane( plPlane1, plPlane2, ptL, vtL) != IPPT_YES) return IPPT_NO ; // intersezione della linea con il terzo piano - return ( IntersLinePlane( ptL, vtL, 1, plPlane3, ptInt) == ILPT_YES ? IPPT_YES : IPPT_NO) ; + return ( IntersLinePlane( ptL, vtL, 1, plPlane3, ptInt, false) == ILPT_YES ? IPPT_YES : IPPT_NO) ; } diff --git a/ProjectCurveSurf.cpp b/ProjectCurveSurf.cpp index 9da5bb9..d37054d 100644 --- a/ProjectCurveSurf.cpp +++ b/ProjectCurveSurf.cpp @@ -19,19 +19,187 @@ #include "/EgtDev/Include/EGkDistPointLine.h" #include "/EgtDev/Include/EGkDistPointCurve.h" #include "/EgtDev/Include/EGkDistPointSurfTm.h" +#include "/EgtDev/Include/EGkIntersPlanePlane.h" +#include "/EgtDev/Include/EGkIntersLinePlane.h" #include "/EgtDev/Include/EGkIntersLineSurfTm.h" #include "/EgtDev/Include/EGkProjectCurveSurf.h" using namespace std ; //---------------------------------------------------------------------------- -// Angolo limite tra normale al triangolo e direzione di proiezione 89° +// Angolo limite tra normale al triangolo e direzione di proiezione 89° const double COS_ANG_LIM = 0.0175 ; +// Massimo numero di triangoli per raffinare ricerca su spigoli +const int MAX_FACET_FOR_CORNER = 1000 ; +// Angolo massimo tra normali per effettuare bisezione su spigolo +const double COS_ANG_MAX_CORNER = 0.8660 ; +// Tipologia di punto +const int P5AX_TO_DELETE = -1 ; // da cancellare +const int P5AX_OUT = 0 ; // aggiunto prima di inizio o dopo fine +const int P5AX_STD = 1 ; // standard +const int P5AX_CVEX = 2 ; // su angolo convesso +const int P5AX_CONC = 3 ; // in angolo concavo +const int P5AX_BEFORE_CONC = 4 ; // adiacente ad angolo concavo +const int P5AX_AFTER_CONC = 5 ; // adiacente ad angolo concavo + +//---------------------------------------------------------------------------- +static bool +PointsInTolerance( const PNT5AXVECTOR& vPt5ax, int nPrec, int nCurr, int nNext, double dSqTol) +{ + for ( int i = nPrec + 1 ; i < nCurr ; ++ i) { + double dSqDist ; + if ( ! DistPointLine( vPt5ax[i].ptP, vPt5ax[nPrec].ptP, vPt5ax[nNext].ptP).GetSqDist( dSqDist) || dSqDist > dSqTol) + return false ; + } + return true ; +} + +//---------------------------------------------------------------------------- +static bool +RemovePointsInExcess( PNT5AXVECTOR& vPt5ax, double dLinTol, double dMaxSegmLen, bool bTestDir) +{ + // Parametri di riferimento + double dSqMaxLen = dMaxSegmLen * dMaxSegmLen ; + double dSqTol = dLinTol * dLinTol ; + const double LENREF = 200 ; + double dCosAngLim = 1 - dSqTol / ( 2 * LENREF * LENREF) ; + // Cerco gli angoli interni e marco opportunamente i punti nelle vicinanze fino ai limiti prima e dopo + int nInd = 0 ; + while ( nInd < int( vPt5ax.size())) { + if ( vPt5ax[nInd].nFlag == P5AX_CONC) { + // analizzo i punti appena precedenti + int nIpv = nInd - 1 ; + while ( nIpv >= 0) { + double dSqLen = SqDist( vPt5ax[nInd].ptP, vPt5ax[nIpv].ptP) ; + if ( dSqLen < dSqMaxLen) + vPt5ax[nIpv].nFlag = P5AX_TO_DELETE ; + else { + vPt5ax[nIpv].nFlag = P5AX_BEFORE_CONC ; + break ; + } + -- nIpv ; + } + // analizzo i punti appena successivi + int nInx = nInd + 1 ; + while ( nInx < int( vPt5ax.size())) { + double dSqLen = SqDist( vPt5ax[nInd].ptP, vPt5ax[nInx].ptP) ; + if ( dSqLen < dSqMaxLen) + vPt5ax[nInx].nFlag = P5AX_TO_DELETE ; + else { + vPt5ax[nInx].nFlag = P5AX_AFTER_CONC ; + break ; + } + ++ nInx ; + } + } + ++ nInd ; + } + // Rimuovo i punti allineati entro la tolleranza e non più lontani tra loro del massimo + int nCnt = 0 ; + int nPrec = 0 ; + int nCurr = 1 ; + int nNext = 2 ; + while ( nNext < int( vPt5ax.size())) { + bool bRemove = false ; + // lunghezza del segmento che unisce gli adiacenti + double dSqLen = SqDist( vPt5ax[nPrec].ptP, vPt5ax[nNext].ptP) ; + // se rimovibile (Flag standard) e lunghezza inferiore al massimo, passo agli altri controlli + if ( vPt5ax[nCurr].nFlag == P5AX_STD && dSqLen <= dSqMaxLen) { + // distanza del punto corrente dal segmento che unisce gli adiacenti + DistPointLine dPL( vPt5ax[nCurr].ptP, vPt5ax[nPrec].ptP, vPt5ax[nNext].ptP) ; + double dSqDist ; + // se distanza inferiore a tolleranza lineare + if ( dPL.GetSqDist( dSqDist) && dSqDist < dSqTol && PointsInTolerance( vPt5ax, nPrec, nCurr, nNext, dSqTol)) { + // verifico se errore angolare inferiore a limite + double dPar ; dPL.GetParamAtMinDistPoint( dPar) ; + if ( bTestDir) { + Vector3d vtNew = Media( vPt5ax[nPrec].vtDir1, vPt5ax[nNext].vtDir1, dPar) ; + if ( vtNew.Normalize() && vtNew * vPt5ax[nCurr].vtDir1 > dCosAngLim) + bRemove = true ; + } + else { + Vector3d vtNew = Media( vPt5ax[nPrec].vtDir2, vPt5ax[nNext].vtDir2, dPar) ; + if ( vtNew.Normalize() && vtNew * vPt5ax[nCurr].vtDir2 > dCosAngLim) + bRemove = true ; + } + } + } + // se da eliminare + if ( bRemove) { + // dichiaro da eliminare il punto + vPt5ax[nCurr].nFlag = P5AX_TO_DELETE ; + // avanzo con corrente e successivo + nCurr = nNext ; + ++ nNext ; + } + // altrimenti da tenere + else { + // avanzo il terzetto di uno step + nPrec = nCurr ; + nCurr = nNext ; + ++ nNext ; + // incremento contatore dei punti conservati + ++ nCnt ; + } + } + + // Copio i punti da conservare in un vettore temporaneo + PNT5AXVECTOR vMyPt5ax ; + vMyPt5ax.reserve( nCnt) ; + for ( const auto& Pt5ax : vPt5ax) { + if ( Pt5ax.nFlag != P5AX_TO_DELETE) + vMyPt5ax.emplace_back( Pt5ax) ; + } + // scambio i due vettori + vPt5ax.swap( vMyPt5ax) ; + + return true ; +} + +//---------------------------------------------------------------------------- +static bool +ProjectPointOnSurf( const Point3d& ptP, const SurfTriMesh& stmSurf, const Frame3d& frRefLine, const IntersParLinesSurfTm& intPLSTM, + double dPar, Point5ax& Pt5ax) +{ + // intersezione retta di proiezione con superficie + Point3d ptL = GetToLoc( ptP, frRefLine) ; + ILSIVECTOR vIntRes ; + intPLSTM.GetInters( ptL, 1, vIntRes, false) ; + // cerco la prima intersezione valida a partire dall'ultima (è la più alta) + int nI = int( vIntRes.size()) - 1 ; + while ( nI >= 0 && abs( vIntRes[nI].dCosDN) < COS_ANG_LIM) + --nI ; + // se trovata + if ( nI >= 0) { + // calcolo il punto + Point3d ptInt ; + if ( vIntRes[nI].nILTT == ILTT_SEGM || vIntRes[nI].nILTT == ILTT_SEGM_ON_EDGE) + ptInt = vIntRes[nI].ptI2 ; + else + ptInt = vIntRes[nI].ptI ; + // calcolo la normale (si calcola smooth, in caso di errore si prende quella del triangolo) + Triangle3dEx trTria ; + if ( ! stmSurf.GetTriangle( vIntRes[nI].nT, trTria)) + return false ; + Vector3d vtN ; + if ( ! CalcNormal( ptInt, trTria, vtN)) + vtN = trTria.GetN() ; + // assegno valori al punto 5assi + Pt5ax.ptP = ptInt ; + Pt5ax.vtDir1 = vtN ; + Pt5ax.vtDir2 = frRefLine.VersZ() ; + Pt5ax.dPar = dPar ; + Pt5ax.nFlag = P5AX_STD ; + // ritorno con successo + return true ; + } + return false ; +} //---------------------------------------------------------------------------- bool ProjectCurveOnSurf( const ICurve& crCrv, const ISurf& sfSurf, const Vector3d& vtDir, - double dLinTol, double dMaxSegmLen, PNT5AXVECTOR& vPt5ax) + double dLinTol, double dMaxSegmLen, bool bSharpEdges, PNT5AXVECTOR& vPt5ax) { // sistemazioni per tipo di superficie const SurfTriMesh* pSurfTm = nullptr ; @@ -45,6 +213,8 @@ ProjectCurveOnSurf( const ICurve& crCrv, const ISurf& sfSurf, const Vector3d& vt case SRF_FLATRGN : pSurfTm = GetBasicSurfFlatRegion( &sfSurf)->GetAuxSurf() ; break ; + default : + break ; } if ( pSurfTm == nullptr) return false ; @@ -66,60 +236,143 @@ ProjectCurveOnSurf( const ICurve& crCrv, const ISurf& sfSurf, const Vector3d& vt return false ; IntersParLinesSurfTm intPLSTM( frRefLine, *pSurfTm) ; - // Vettore locale dei punti risultanti - PNT5AXVECTOR vMyPt5ax ; - vMyPt5ax.reserve( PL.GetPointNbr()) ; + // Pulisco e riservo spazio nel vettore dei punti risultanti + vPt5ax.clear() ; + vPt5ax.reserve( PL.GetPointNbr()) ; // proietto i punti della polilinea sulla superficie double dPar ; Point3d ptP ; bool bFound = PL.GetFirstUPoint( &dPar, &ptP) ; while ( bFound) { - // intersezione retta di proiezione con superficie - Point3d ptL = GetToLoc( ptP, frRefLine) ; - ILSIVECTOR vIntRes ; - intPLSTM.GetInters( ptL, 1, vIntRes, false) ; - // cerco la prima intersezione valida a partire dall'ultima (è la più alta) - int nI = int( vIntRes.size()) - 1 ; - while ( nI >= 0 && abs( vIntRes[nI].dCosDN) < COS_ANG_LIM) - --nI ; - // se trovata - if ( nI >= 0) { - // calcolo il punto - Point3d ptInt ; - if ( vIntRes[nI].nILTT == ILTT_SEGM || vIntRes[nI].nILTT == ILTT_SEGM_ON_EDGE) - ptInt = vIntRes[nI].ptI2 ; - else - ptInt = vIntRes[nI].ptI ; - // calcolo la normale (si calcola smooth, in caso di errore si prende quella del triangolo) - Triangle3dEx trTria ; - if ( ! pSurfTm->GetTriangle( vIntRes[nI].nT, trTria)) - return false ; - Vector3d vtN ; - if ( ! CalcNormal( ptInt, trTria, vtN)) - vtN = trTria.GetN() ; - // aggiungo al vettore dei proiettati - vMyPt5ax.emplace_back( ptInt, vtN, vtDir, dPar, 1) ; - } + // se trovo proiezione, la salvo + Point5ax Pt5ax ; + if ( ProjectPointOnSurf( ptP, *pSurfTm, frRefLine, intPLSTM, dPar, Pt5ax)) + vPt5ax.emplace_back( Pt5ax) ; + // passo al successivo bFound = PL.GetNextUPoint( &dPar, &ptP) ; } - // eventuale rimozione punti in eccesso rispetto alle tolleranze lasciata alla funzione chiamante - - // copio i punti nel vettore di ritorno - vPt5ax.clear() ; - for ( const auto& Pt5ax : vMyPt5ax) { - if ( Pt5ax.nFlag != -1) - vPt5ax.emplace_back( Pt5ax) ; + // se richiesto, inserimento punti intermedi in presenza di spigoli + if ( bSharpEdges) { + for ( int i = 1 ; i < int( vPt5ax.size()) ; ++ i) { + // precedente + int j = i - 1 ; + // se normali tra corrente e precedente oltre limite e punti abbastanza lontani + if ( vPt5ax[i].vtDir1 * vPt5ax[j].vtDir1 < COS_ANG_MAX_CORNER && + SqDist( vPt5ax[i].ptP, vPt5ax[j].ptP) > 25 * SQ_EPS_SMALL) { + // intersezione tra le due facce + Plane3d plPlane1 ; plPlane1.Set( vPt5ax[j].ptP, vPt5ax[j].vtDir1) ; + Plane3d plPlane2 ; plPlane2.Set( vPt5ax[i].ptP, vPt5ax[i].vtDir1) ; + Point3d ptEdge ; Vector3d vtEdge ; + if ( IntersPlanePlane( plPlane1, plPlane2, ptEdge, vtEdge) == IPPT_YES) { + Plane3d plPlane3 ; plPlane3.Set( vPt5ax[i].ptP, ( vPt5ax[i].ptP - vPt5ax[j].ptP) ^ vtDir) ; + Point3d ptInt ; + if ( IntersLinePlane( ptEdge, vtEdge, 1, plPlane3, ptInt, false) == ILPT_YES) { + // verifico se spigolo convesso o concavo + bool bConvex ; + if ( ! AreSamePointApprox( ptInt, vPt5ax[j].ptP)) + bConvex = ( ( vPt5ax[j].vtDir1 ^ ( ptInt - vPt5ax[j].ptP)) * vtEdge > 0) ; + else + bConvex = (( vPt5ax[i].vtDir1 ^ ( ptInt - vPt5ax[i].ptP)) * vtEdge < 0) ; + // se convesso, metto due punti con direzione appena prima e appena dopo + if ( bConvex) { + Vector3d vtLine1 = ptInt - vPt5ax[j].ptP ; double dLen1 = vtLine1.Len() ; + Vector3d vtLine2 = vPt5ax[i].ptP - ptInt ; double dLen2 = vtLine2.Len() ; + if ( dLen1 > 2 * EPS_SMALL) { + Point5ax Pt5ax ; + Pt5ax.ptP = ptInt - vtLine1 / dLen1 * 2 * EPS_SMALL ; + Pt5ax.vtDir1 = vPt5ax[j].vtDir1 ; + Pt5ax.vtDir2 = vtDir ; + Pt5ax.dPar = ( vPt5ax[i].dPar + vPt5ax[j].dPar) / 2 ; + Pt5ax.nFlag = P5AX_CVEX ; + vPt5ax.insert( vPt5ax.begin() + i, Pt5ax) ; + ++ i ; + } + if ( dLen2 > 2 * EPS_SMALL) { + Point5ax Pt5ax ; + Pt5ax.ptP = ptInt + vtLine2 / dLen2 * 2 * EPS_SMALL ; + Pt5ax.vtDir1 = vPt5ax[i].vtDir1 ; + Pt5ax.vtDir2 = vtDir ; + Pt5ax.dPar = ( vPt5ax[i].dPar + vPt5ax[j].dPar) / 2 ; + Pt5ax.nFlag = P5AX_CVEX ; + vPt5ax.insert( vPt5ax.begin() + i, Pt5ax) ; + ++ i ; + } + } + // altrimenti concavo, aggiungo un solo punto con la direzione media + else { + Point5ax Pt5ax ; + Pt5ax.ptP = ptInt ; + Pt5ax.vtDir1 = Media( vPt5ax[i].vtDir1, vPt5ax[j].vtDir1) ; + Pt5ax.vtDir2 = vtDir ; + Pt5ax.dPar = ( vPt5ax[i].dPar + vPt5ax[j].dPar) / 2 ; + Pt5ax.nFlag = P5AX_CONC ; + vPt5ax.insert( vPt5ax.begin() + i, Pt5ax) ; + ++ i ; + } + } + } + } + } } + // rimozione punti in eccesso rispetto alle tolleranze + RemovePointsInExcess( vPt5ax, dLinTol, dMaxSegmLen, bSharpEdges) ; + return true ; } +//---------------------------------------------------------------------------- +static bool +ProjectPointOnSurf( const Point3d& ptP, const SurfTriMesh& stmSurf, const IGeoPoint3d& gpRef, double dPar, Point5ax& Pt5ax) +{ + // punto di riferimento + Point3d ptMin = gpRef.GetPoint() ; + // intersezione della retta di minima distanza con la superficie + Vector3d vtLine = ptP - ptMin ; + double dLineLen = vtLine.Len() ; + if ( dLineLen > EPS_SMALL) { + vtLine /= dLineLen ; + ILSIVECTOR vIntRes ; + if ( IntersLineSurfTm( ptP, vtLine, dLineLen, stmSurf, vIntRes, false)) { + // cerco la prima intersezione valida a partire dall'ultima (è la più alta) + int nI = int( vIntRes.size()) - 1 ; + while ( nI >= 0 && abs( vIntRes[nI].dCosDN) < COS_ANG_LIM) + --nI ; + // se trovata + if ( nI >= 0) { + // calcolo il punto + Point3d ptInt ; + if ( vIntRes[nI].nILTT == ILTT_SEGM || vIntRes[nI].nILTT == ILTT_SEGM_ON_EDGE) + ptInt = vIntRes[nI].ptI2 ; + else + ptInt = vIntRes[nI].ptI ; + // calcolo la normale (si calcola smooth, in caso di errore si prende quella del triangolo) + Triangle3dEx trTria ; + if ( ! stmSurf.GetTriangle( vIntRes[nI].nT, trTria)) + return false ; + Vector3d vtN ; + if ( ! CalcNormal( ptInt, trTria, vtN)) + vtN = trTria.GetN() ; + // assegno valori al punto 5assi + Pt5ax.ptP = ptInt ; + Pt5ax.vtDir1 = vtN ; + Pt5ax.vtDir2 = vtLine ; + Pt5ax.dPar = dPar ; + Pt5ax.nFlag = P5AX_STD ; + // ritorno con successo + return true ; + } + } + } + return false ; +} + //---------------------------------------------------------------------------- bool ProjectCurveOnSurf( const ICurve& crCrv, const ISurf& sfSurf, const IGeoPoint3d& gpRef, - double dLinTol, double dMaxSegmLen, PNT5AXVECTOR& vPt5ax) + double dLinTol, double dMaxSegmLen, bool bSharpEdges, PNT5AXVECTOR& vPt5ax) { // sistemazioni per tipo di superficie const SurfTriMesh* pSurfTm = nullptr ; @@ -133,6 +386,8 @@ ProjectCurveOnSurf( const ICurve& crCrv, const ISurf& sfSurf, const IGeoPoint3d& case SRF_FLATRGN : pSurfTm = GetBasicSurfFlatRegion( &sfSurf)->GetAuxSurf() ; break ; + default : + break ; } if ( pSurfTm == nullptr) return false ; @@ -141,7 +396,7 @@ ProjectCurveOnSurf( const ICurve& crCrv, const ISurf& sfSurf, const IGeoPoint3d& dLinTol = max( dLinTol, LIN_TOL_MIN) ; dMaxSegmLen = max( dMaxSegmLen, 10 * EPS_SMALL) ; - // approssimo la curva con una polilinea entro la metà della tolleranza + // approssimo la curva con una polilinea entro la metà della tolleranza PolyLine PL ; if ( ! crCrv.ApproxWithLines( dLinTol, ANG_TOL_STD_DEG, ICurve::APL_STD, PL)) return false ; @@ -149,25 +404,67 @@ ProjectCurveOnSurf( const ICurve& crCrv, const ISurf& sfSurf, const IGeoPoint3d& if ( ! PL.AdjustForMaxSegmentLen( MAX_SEG_LEN)) return false ; - // Vettore locale dei punti risultanti - PNT5AXVECTOR vMyPt5ax ; - vMyPt5ax.reserve( PL.GetPointNbr()) ; + // Pulisco e riservo spazio nel vettore dei punti risultanti + vPt5ax.clear() ; + vPt5ax.reserve( PL.GetPointNbr()) ; // proietto i punti della polilinea sulla superficie con direzione data dal punto di riferimento double dPar ; Point3d ptP ; bool bFound = PL.GetFirstUPoint( &dPar, &ptP) ; while ( bFound) { - // punto di riferimento - Point3d ptMin = gpRef.GetPoint() ; + // se trovo proiezione, la salvo + Point5ax Pt5ax ; + if ( ProjectPointOnSurf( ptP, *pSurfTm, gpRef, dPar, Pt5ax)) + vPt5ax.emplace_back( Pt5ax) ; + // passo al successivo + bFound = PL.GetNextUPoint( &dPar, &ptP) ; + } + + // se superfici con non troppi triangoli, inserimento punti intermedi in presenza di spigoli + if ( pSurfTm->GetFacetCount() < MAX_FACET_FOR_CORNER) { + for ( int i = 1 ; i < int( vPt5ax.size()) ; ++ i) { + // precedente + int j = i - 1 ; + // se normali tra corrente e precedente oltre limite e punti abbastanza lontani + if ( vPt5ax[i].vtDir1 * vPt5ax[j].vtDir1 < COS_ANG_MAX_CORNER && + SqDist( vPt5ax[i].ptP, vPt5ax[j].ptP) > 25 * SQ_EPS_SMALL) { + // punto medio + Point3d ptMid = Media( vPt5ax[i].ptP, vPt5ax[j].ptP) ; + double dMid = ( vPt5ax[i].dPar + vPt5ax[j].dPar) / 2 ; + // se trovo proiezione, la salvo + Point5ax Pt5ax ; + if ( ProjectPointOnSurf( ptMid, *pSurfTm, gpRef, dMid, Pt5ax)) { + vPt5ax.insert( vPt5ax.begin() + i, Pt5ax) ; + -- i ; + } + } + } + } + + // rimozione punti in eccesso rispetto alle tolleranze + RemovePointsInExcess( vPt5ax, dLinTol, dMaxSegmLen, bSharpEdges) ; + + return true ; +} + +//---------------------------------------------------------------------------- +static bool +ProjectPointOnSurf( const Point3d& ptP, const SurfTriMesh& stmSurf, const ICurve& crRef, double dPar, Point5ax& Pt5ax) +{ + // punto a minima distanza + DistPointCurve dPC( ptP, crRef) ; + Point3d ptMin ; + int nFlag ; + if ( dPC.GetMinDistPoint( 0, ptMin, nFlag)) { // intersezione della retta di minima distanza con la superficie Vector3d vtLine = ptP - ptMin ; double dLineLen = vtLine.Len() ; if ( dLineLen > EPS_SMALL) { vtLine /= dLineLen ; ILSIVECTOR vIntRes ; - if ( IntersLineSurfTm( ptP, vtLine, dLineLen, *pSurfTm, vIntRes, false)) { - // cerco la prima intersezione valida a partire dall'ultima (è la più alta) + if ( IntersLineSurfTm( ptP, vtLine, dLineLen, stmSurf, vIntRes, false)) { + // cerco la prima intersezione valida a partire dall'ultima (è la più alta) int nI = int( vIntRes.size()) - 1 ; while ( nI >= 0 && abs( vIntRes[nI].dCosDN) < COS_ANG_LIM) --nI ; @@ -181,37 +478,32 @@ ProjectCurveOnSurf( const ICurve& crCrv, const ISurf& sfSurf, const IGeoPoint3d& ptInt = vIntRes[nI].ptI ; // calcolo la normale (si calcola smooth, in caso di errore si prende quella del triangolo) Triangle3dEx trTria ; - if ( ! pSurfTm->GetTriangle( vIntRes[nI].nT, trTria)) + if ( ! stmSurf.GetTriangle( vIntRes[nI].nT, trTria)) return false ; Vector3d vtN ; if ( ! CalcNormal( ptInt, trTria, vtN)) vtN = trTria.GetN() ; - // aggiungo al vettore dei proiettati - vMyPt5ax.emplace_back( ptInt, vtN, vtLine, dPar, 1) ; + // assegno valori al punto 5assi + Pt5ax.ptP = ptInt ; + Pt5ax.vtDir1 = vtN ; + Pt5ax.vtDir2 = vtLine ; + Pt5ax.dPar = dPar ; + Pt5ax.nFlag = P5AX_STD ; + // ritorno con successo + return true ; } } } - bFound = PL.GetNextUPoint( &dPar, &ptP) ; } - - // eventuale rimozione punti in eccesso rispetto alle tolleranze lasciata alla funzione chiamante - - // copio i punti nel vettore di ritorno - vPt5ax.clear() ; - for ( const auto& Pt5ax : vMyPt5ax) { - if ( Pt5ax.nFlag != -1) - vPt5ax.emplace_back( Pt5ax) ; - } - - return true ; + return false ; } //---------------------------------------------------------------------------- bool ProjectCurveOnSurf( const ICurve& crCrv, const ISurf& sfSurf, const ICurve& crRef, - double dLinTol, double dMaxSegmLen, PNT5AXVECTOR& vPt5ax) + double dLinTol, double dMaxSegmLen, bool bSharpEdges, PNT5AXVECTOR& vPt5ax) { - // sistemazioni per tipo di superficie + // Sistemazioni per tipo di superficie const SurfTriMesh* pSurfTm = nullptr ; switch ( sfSurf.GetType()) { case SRF_TRIMESH : @@ -223,15 +515,17 @@ ProjectCurveOnSurf( const ICurve& crCrv, const ISurf& sfSurf, const ICurve& crRe case SRF_FLATRGN : pSurfTm = GetBasicSurfFlatRegion( &sfSurf)->GetAuxSurf() ; break ; + default : + break ; } if ( pSurfTm == nullptr) return false ; - // controllo le tolleranze + // Controllo le tolleranze dLinTol = max( dLinTol, LIN_TOL_MIN) ; dMaxSegmLen = max( dMaxSegmLen, 10 * EPS_SMALL) ; - // approssimo la curva con una polilinea alla massima risoluzione + // Approssimo la curva con una polilinea alla massima risoluzione PolyLine PL ; if ( ! crCrv.ApproxWithLines( EPS_SMALL, ANG_TOL_STD_DEG, ICurve::APL_STD, PL)) return false ; @@ -239,71 +533,163 @@ ProjectCurveOnSurf( const ICurve& crCrv, const ISurf& sfSurf, const ICurve& crRe if ( ! PL.AdjustForMaxSegmentLen( MAX_SEG_LEN)) return false ; - // Vettore locale dei punti risultanti - PNT5AXVECTOR vMyPt5ax ; - vMyPt5ax.reserve( PL.GetPointNbr()) ; + // Pulisco e riservo spazio nel vettore dei punti risultanti + vPt5ax.clear() ; + vPt5ax.reserve( PL.GetPointNbr()) ; // proietto i punti della polilinea sulla superficie con direzione normale alla curva di riferimento double dPar ; Point3d ptP ; bool bFound = PL.GetFirstUPoint( &dPar, &ptP) ; while ( bFound) { - // punto sulla curva a minima distanza - DistPointCurve dPC( ptP, crRef) ; - Point3d ptMin ; - int nFlag ; - if ( dPC.GetMinDistPoint( 0, ptMin, nFlag)) { - // intersezione della retta di minima distanza con la superficie - Vector3d vtLine = ptP - ptMin ; - double dLineLen = vtLine.Len() ; - if ( dLineLen > EPS_SMALL) { - vtLine /= dLineLen ; - ILSIVECTOR vIntRes ; - if ( IntersLineSurfTm( ptP, vtLine, dLineLen, *pSurfTm, vIntRes, false)) { - // cerco la prima intersezione valida a partire dall'ultima (è la più alta) - int nI = int( vIntRes.size()) - 1 ; - while ( nI >= 0 && abs( vIntRes[nI].dCosDN) < COS_ANG_LIM) - --nI ; - // se trovata - if ( nI >= 0) { - // calcolo il punto - Point3d ptInt ; - if ( vIntRes[nI].nILTT == ILTT_SEGM || vIntRes[nI].nILTT == ILTT_SEGM_ON_EDGE) - ptInt = vIntRes[nI].ptI2 ; + // se trovo proiezione, la salvo + Point5ax Pt5ax ; + if ( ProjectPointOnSurf( ptP, *pSurfTm, crRef, dPar, Pt5ax)) + vPt5ax.emplace_back( Pt5ax) ; + // passo al successivo + bFound = PL.GetNextUPoint( &dPar, &ptP) ; + } + + // se richiesto, inserimento punti intermedi in presenza di spigoli + if ( bSharpEdges) { + for ( int i = 1 ; i < int( vPt5ax.size()) ; ++ i) { + // precedente + int j = i - 1 ; + // se normali tra corrente e precedente oltre limite e punti abbastanza lontani + if ( vPt5ax[i].vtDir1 * vPt5ax[j].vtDir1 < COS_ANG_MAX_CORNER && + SqDist( vPt5ax[i].ptP, vPt5ax[j].ptP) > 25 * SQ_EPS_SMALL) { + // intersezione tra le due facce + Plane3d plPlane1 ; plPlane1.Set( vPt5ax[j].ptP, vPt5ax[j].vtDir1) ; + Plane3d plPlane2 ; plPlane2.Set( vPt5ax[i].ptP, vPt5ax[i].vtDir1) ; + Point3d ptEdge ; Vector3d vtEdge ; + if ( IntersPlanePlane( plPlane1, plPlane2, ptEdge, vtEdge) == IPPT_YES) { + Plane3d plPlane3 ; plPlane3.Set( vPt5ax[i].ptP, vPt5ax[i].vtDir2 ^ vPt5ax[j].vtDir2) ; + Point3d ptInt ; + if ( IntersLinePlane( ptEdge, vtEdge, 1, plPlane3, ptInt, false) == ILPT_YES) { + // verifico se spigolo convesso o concavo + bool bConvex ; + if ( ! AreSamePointApprox( ptInt, vPt5ax[j].ptP)) + bConvex = ( ( vPt5ax[j].vtDir1 ^ ( ptInt - vPt5ax[j].ptP)) * vtEdge > 0) ; else - ptInt = vIntRes[nI].ptI ; - // calcolo la normale (si calcola smooth, in caso di errore si prende quella del triangolo) - Triangle3dEx trTria ; - if ( ! pSurfTm->GetTriangle( vIntRes[nI].nT, trTria)) - return false ; - Vector3d vtN ; - if ( ! CalcNormal( ptInt, trTria, vtN)) - vtN = trTria.GetN() ; - // aggiungo al vettore dei proiettati - vMyPt5ax.emplace_back( ptInt, vtN, vtLine, dPar, 1) ; + bConvex = (( vPt5ax[i].vtDir1 ^ ( ptInt - vPt5ax[i].ptP)) * vtEdge < 0) ; + // se convesso, metto due punti con direzione appena prima e appena dopo + if ( bConvex) { + Vector3d vtLine1 = ptInt - vPt5ax[j].ptP ; double dLen1 = vtLine1.Len() ; + Vector3d vtLine2 = vPt5ax[i].ptP - ptInt ; double dLen2 = vtLine2.Len() ; + if ( dLen1 > 2 * EPS_SMALL) { + Point5ax Pt5ax ; + Pt5ax.ptP = ptInt - vtLine1 / dLen1 * 2 * EPS_SMALL ; + Pt5ax.vtDir1 = vPt5ax[j].vtDir1 ; + Pt5ax.vtDir2 = Media( vPt5ax[i].vtDir2, vPt5ax[j].vtDir2) ; + Pt5ax.dPar = ( vPt5ax[i].dPar + vPt5ax[j].dPar) / 2 ; + Pt5ax.nFlag = P5AX_CVEX ; + vPt5ax.insert( vPt5ax.begin() + i, Pt5ax) ; + ++ i ; + } + if ( dLen2 > 2 * EPS_SMALL) { + Point5ax Pt5ax ; + Pt5ax.ptP = ptInt + vtLine2 / dLen2 * 2 * EPS_SMALL ; + Pt5ax.vtDir1 = vPt5ax[i].vtDir1 ; + Pt5ax.vtDir2 = Media( vPt5ax[i].vtDir2, vPt5ax[j].vtDir2) ; + Pt5ax.dPar = ( vPt5ax[i].dPar + vPt5ax[j].dPar) / 2 ; + Pt5ax.nFlag = P5AX_CVEX ; + vPt5ax.insert( vPt5ax.begin() + i, Pt5ax) ; + ++ i ; + } + } + // altrimenti concavo, aggiungo un solo punto con la direzione media + else { + Point5ax Pt5ax ; + Pt5ax.ptP = ptInt ; + Pt5ax.vtDir1 = Media( vPt5ax[i].vtDir1, vPt5ax[j].vtDir1) ; + Pt5ax.vtDir2 = Media( vPt5ax[i].vtDir2, vPt5ax[j].vtDir2) ; + Pt5ax.dPar = ( vPt5ax[i].dPar + vPt5ax[j].dPar) / 2 ; + Pt5ax.nFlag = P5AX_CONC ; + vPt5ax.insert( vPt5ax.begin() + i, Pt5ax) ; + ++ i ; + } } } } } - bFound = PL.GetNextUPoint( &dPar, &ptP) ; } - // eventuale rimozione punti in eccesso rispetto alle tolleranze lasciata alla funzione chiamante - - // copio i punti nel vettore di ritorno - vPt5ax.clear() ; - for ( const auto& Pt5ax : vMyPt5ax) { - if ( Pt5ax.nFlag != -1) - vPt5ax.emplace_back( Pt5ax) ; - } + // rimozione punti in eccesso rispetto alle tolleranze + RemovePointsInExcess( vPt5ax, dLinTol, dMaxSegmLen, bSharpEdges) ; return true ; } +//---------------------------------------------------------------------------- +static bool +ProjectPointOnSurf( const Point3d& ptP, const SurfTriMesh& stmSurf, const SurfTriMesh& stmRef, double dPar, Point5ax& Pt5ax) +{ + // punto sulla superficie guida a minima distanza + DistPointSurfTm dPS( ptP, stmRef) ; + Point3d ptMin ; + int nTriaMin ; + if ( dPS.GetMinDistPoint( ptMin) && dPS.GetMinDistTriaIndex ( nTriaMin)) { + // recupero direzione della retta di minima distanza, altrimenti normale alla superficie + Vector3d vtLine = ptP - ptMin ; + double dLineLen = vtLine.Len() ; + if ( dLineLen > EPS_SMALL) + vtLine /= dLineLen ; + else { + // calcolo la normale della superficie guida + Triangle3dEx trGuide ; + if ( ! stmRef.GetTriangle( nTriaMin, trGuide)) + return false ; + if ( ! CalcNormal( ptMin, trGuide, vtLine)) + vtLine = trGuide.GetN() ; + dLineLen = 100 ; + } + // intersezione della retta con la superficie + ILSIVECTOR vIntRes ; + if ( IntersLineSurfTm( ptP, vtLine, dLineLen, stmSurf, vIntRes, false)) { + // cerco la prima intersezione valida a partire dall'ultima (è la più alta) + int nI = int( vIntRes.size()) - 1 ; + while ( nI >= 0 && abs( vIntRes[nI].dCosDN) < COS_ANG_LIM) + --nI ; + // se trovata + if ( nI >= 0) { + // calcolo il punto + Point3d ptInt ; + if ( vIntRes[nI].nILTT == ILTT_SEGM || vIntRes[nI].nILTT == ILTT_SEGM_ON_EDGE) + ptInt = vIntRes[nI].ptI2 ; + else + ptInt = vIntRes[nI].ptI ; + // calcolo la normale (si calcola smooth, in caso di errore si prende quella del triangolo) + Triangle3dEx trTria ; + if ( ! stmSurf.GetTriangle( vIntRes[nI].nT, trTria)) + return false ; + Vector3d vtN ; + if ( ! CalcNormal( ptMin, trTria, vtN)) + vtN = trTria.GetN() ; + // calcolo la normale della superficie guida + Triangle3dEx trGuide ; + if ( ! stmRef.GetTriangle( nTriaMin, trGuide)) + return false ; + Vector3d vtN2 ; + if ( ! CalcNormal( ptMin, trGuide, vtN2)) + vtN2 = trGuide.GetN() ; + // assegno valori al punto 5assi + Pt5ax.ptP = ptInt ; + Pt5ax.vtDir1 = vtN ; + Pt5ax.vtDir2 = vtN2 ; + Pt5ax.dPar = dPar ; + Pt5ax.nFlag = P5AX_STD ; + // ritorno con successo + return true ; + } + } + } + return false ; +} + //---------------------------------------------------------------------------- bool ProjectCurveOnSurf( const ICurve& crCrv, const ISurf& sfSurf, const ISurf& sfRef, - double dLinTol, double dMaxSegmLen, PNT5AXVECTOR& vPt5ax) + double dLinTol, double dMaxSegmLen, bool bSharpEdges, PNT5AXVECTOR& vPt5ax) { // sistemazioni per tipo di superficie const SurfTriMesh* pSurfTm = nullptr ; @@ -317,6 +703,8 @@ ProjectCurveOnSurf( const ICurve& crCrv, const ISurf& sfSurf, const ISurf& sfRef case SRF_FLATRGN : pSurfTm = GetBasicSurfFlatRegion( &sfSurf)->GetAuxSurf() ; break ; + default : + break ; } if ( pSurfTm == nullptr) return false ; @@ -333,6 +721,8 @@ ProjectCurveOnSurf( const ICurve& crCrv, const ISurf& sfSurf, const ISurf& sfRef case SRF_FLATRGN : pRefTm = GetBasicSurfFlatRegion( &sfRef)->GetAuxSurf() ; break ; + default : + break ; } if ( pRefTm == nullptr) return false ; @@ -341,7 +731,7 @@ ProjectCurveOnSurf( const ICurve& crCrv, const ISurf& sfSurf, const ISurf& sfRef dLinTol = max( dLinTol, LIN_TOL_MIN) ; dMaxSegmLen = max( dMaxSegmLen, 10 * EPS_SMALL) ; - // approssimo la curva con una polilinea entro la metà della tolleranza + // approssimo la curva con una polilinea entro la metà della tolleranza PolyLine PL ; if ( ! crCrv.ApproxWithLines( dLinTol, ANG_TOL_STD_DEG, ICurve::APL_STD, PL)) return false ; @@ -349,147 +739,46 @@ ProjectCurveOnSurf( const ICurve& crCrv, const ISurf& sfSurf, const ISurf& sfRef if ( ! PL.AdjustForMaxSegmentLen( MAX_SEG_LEN)) return false ; - // Vettore locale dei punti risultanti - PNT5AXVECTOR vMyPt5ax ; - vMyPt5ax.reserve( PL.GetPointNbr()) ; + // Pulisco e riservo spazio nel vettore dei punti risultanti + vPt5ax.clear() ; + vPt5ax.reserve( PL.GetPointNbr()) ; // proietto i punti della polilinea sulla superficie con direzione normale alla curva di riferimento double dPar ; Point3d ptP ; bool bFound = PL.GetFirstUPoint( &dPar, &ptP) ; while ( bFound) { - // punto sulla superficie guida a minima distanza - DistPointSurfTm dPS( ptP, *pRefTm) ; - Point3d ptMin ; - int nTriaMin ; - if ( dPS.GetMinDistPoint( ptMin) && dPS.GetMinDistTriaIndex ( nTriaMin)) { - // recupero direzione della retta di minima distanza, altrimenti normale alla superficie - Vector3d vtLine = ptP - ptMin ; - double dLineLen = vtLine.Len() ; - if ( dLineLen > EPS_SMALL) - vtLine /= dLineLen ; - else { - // calcolo la normale della superficie guida - Triangle3dEx trGuide ; - if ( ! pRefTm->GetTriangle( nTriaMin, trGuide)) - return false ; - if ( ! CalcNormal( ptMin, trGuide, vtLine)) - vtLine = trGuide.GetN() ; - dLineLen = 100 ; - } - // intersezione della retta con la superficie - ILSIVECTOR vIntRes ; - if ( IntersLineSurfTm( ptP, vtLine, dLineLen, *pSurfTm, vIntRes, false)) { - // cerco la prima intersezione valida a partire dall'ultima (è la più alta) - int nI = int( vIntRes.size()) - 1 ; - while ( nI >= 0 && abs( vIntRes[nI].dCosDN) < COS_ANG_LIM) - --nI ; - // se trovata - if ( nI >= 0) { - // calcolo il punto - Point3d ptInt ; - if ( vIntRes[nI].nILTT == ILTT_SEGM || vIntRes[nI].nILTT == ILTT_SEGM_ON_EDGE) - ptInt = vIntRes[nI].ptI2 ; - else - ptInt = vIntRes[nI].ptI ; - // calcolo la normale (si calcola smooth, in caso di errore si prende quella del triangolo) - Triangle3dEx trTria ; - if ( ! pSurfTm->GetTriangle( vIntRes[nI].nT, trTria)) - return false ; - Vector3d vtN ; - if ( ! CalcNormal( ptMin, trTria, vtN)) - vtN = trTria.GetN() ; - // calcolo la normale della superficie guida - Triangle3dEx trGuide ; - if ( ! pRefTm->GetTriangle( nTriaMin, trGuide)) - return false ; - Vector3d vtN2 ; - if ( ! CalcNormal( ptMin, trGuide, vtN2)) - vtN2 = trGuide.GetN() ; - // aggiungo al vettore dei proiettati - vMyPt5ax.emplace_back( ptInt, vtN, vtN2, dPar, 1) ; - } - } - } + // se trovo proiezione, la salvo + Point5ax Pt5ax ; + if ( ProjectPointOnSurf( ptP, *pSurfTm, *pRefTm, dPar, Pt5ax)) + vPt5ax.emplace_back( Pt5ax) ; + // passo al successivo bFound = PL.GetNextUPoint( &dPar, &ptP) ; } - // eventuale rimozione punti in eccesso rispetto alle tolleranze lasciata alla funzione chiamante - - // copio i punti nel vettore di ritorno - vPt5ax.clear() ; - for ( const auto& Pt5ax : vMyPt5ax) { - if ( Pt5ax.nFlag != -1) - vPt5ax.emplace_back( Pt5ax) ; - } - - return true ; -} - -//---------------------------------------------------------------------------- -static bool -PointsInTolerance( const PNT5AXVECTOR& vPt5ax, int nPrec, int nCurr, int nNext, double dSqTol) -{ - for ( int i = nPrec + 1 ; i < nCurr ; ++ i) { - double dSqDist ; - if ( ! DistPointLine( vPt5ax[i].ptP, vPt5ax[nPrec].ptP, vPt5ax[nNext].ptP).GetSqDist( dSqDist) || dSqDist > dSqTol) - return false ; - } - return true ; -} - -//---------------------------------------------------------------------------- -bool -RemovePointsInExcess( PNT5AXVECTOR& vMyPt5ax, double dLinTol, double dMaxSegmLen, bool bTestDir) -{ - // rimuovo i punti allineati entro la tolleranza e non più lontani tra loro del massimo - double dSqMaxLen = dMaxSegmLen * dMaxSegmLen ; - double dSqTol = dLinTol * dLinTol ; - const double LENREF = 100 ; - double dCosAngLim = 1 - dSqTol / ( 2 * LENREF * LENREF) ; - int nPrec = 0 ; - int nCurr = 1 ; - int nNext = 2 ; - while ( nNext < int( vMyPt5ax.size())) { - bool bRemove = false ; - // lunghezza del segmento che unisce gli adiacenti - double dSqLen = SqDist( vMyPt5ax[nPrec].ptP, vMyPt5ax[nNext].ptP) ; - // se lunghezza inferiore al massimo, passo agli altri controlli - if ( dSqLen <= dSqMaxLen) { - // distanza del punto corrente dal segmento che unisce gli adiacenti - DistPointLine dPL( vMyPt5ax[nCurr].ptP, vMyPt5ax[nPrec].ptP, vMyPt5ax[nNext].ptP) ; - double dSqDist ; - // se distanza inferiore a tolleranza lineare - if ( dPL.GetSqDist( dSqDist) && dSqDist < dSqTol && PointsInTolerance( vMyPt5ax, nPrec, nCurr, nNext, dSqTol)) { - // verifico se errore angolare inferiore a limite - double dPar ; dPL.GetParamAtMinDistPoint( dPar) ; - if ( bTestDir) { - Vector3d vtNew = Media( vMyPt5ax[nPrec].vtDir, vMyPt5ax[nNext].vtDir, dPar) ; - if ( vtNew.Normalize() && vtNew * vMyPt5ax[nCurr].vtDir > dCosAngLim) - bRemove = true ; - } - else { - Vector3d vtNew = Media( vMyPt5ax[nPrec].vtDir2, vMyPt5ax[nNext].vtDir2, dPar) ; - if ( vtNew.Normalize() && vtNew * vMyPt5ax[nCurr].vtDir2 > dCosAngLim) - bRemove = true ; + // se superfici con non troppi triangoli, inserimento punti intermedi in presenza di spigoli + if ( pSurfTm->GetFacetCount() < MAX_FACET_FOR_CORNER) { + for ( int i = 1 ; i < int( vPt5ax.size()) ; ++ i) { + // precedente + int j = i - 1 ; + // se normali tra corrente e precedente oltre limite e punti abbastanza lontani + if ( vPt5ax[i].vtDir1 * vPt5ax[j].vtDir1 < COS_ANG_MAX_CORNER && + SqDist( vPt5ax[i].ptP, vPt5ax[j].ptP) > 25 * SQ_EPS_SMALL) { + // punto medio + Point3d ptMid = Media( vPt5ax[i].ptP, vPt5ax[j].ptP) ; + double dMid = ( vPt5ax[i].dPar + vPt5ax[j].dPar) / 2 ; + // se trovo proiezione, la salvo + Point5ax Pt5ax ; + if ( ProjectPointOnSurf( ptMid, *pSurfTm, *pRefTm, dMid, Pt5ax)) { + vPt5ax.insert( vPt5ax.begin() + i, Pt5ax) ; + -- i ; } } } - // se da eliminare - if ( bRemove) { - // dichiaro da eliminare il punto - vMyPt5ax[nCurr].nFlag = -1 ; - // avanzo con corrente e successivo - nCurr = nNext ; - ++ nNext ; - } - // altrimenti da tenere - else { - // avanzo il terzetto di uno step - nPrec = nCurr ; - nCurr = nNext ; - ++ nNext ; - } } + + // rimozione punti in eccesso rispetto alle tolleranze + RemovePointsInExcess( vPt5ax, dLinTol, dMaxSegmLen, bSharpEdges) ; + return true ; }