//---------------------------------------------------------------------------- // EgalTech 2024 //---------------------------------------------------------------------------- // File : IntersLineSurfBez.cpp Data : 06.02.24 Versione : 2.6b1 // Contenuto : Implementazione della intersezione linea/superficie bezier. // // // // Modifiche : 06.02.24 DB Creazione modulo. // // //---------------------------------------------------------------------------- //--------------------------- Include ---------------------------------------- #include "stdafx.h" #include "CurveLine.h" #include "CurveBezier.h" #include "CurveComposite.h" #include "SurfFlatRegion.h" #include "/EgtDev/Include/EGkDistPointLine.h" #include "/EgtDev/Include/EGkDistLineLine.h" #include "/EgtDev/Include/EGkDistPointSurfFr.h" #include "/EgtDev/Include/EGkIntersLineTria.h" #include "/EgtDev/Include/EGkIntersLineSurfTm.h" #include "/EgtDev/Include/EGkIntersLineSurfBez.h" #include "/EgtDev/Include/EGkSurfBezier.h" #include "/EgtDev/Include/ENkPolynomialRoots.h" using namespace std ; //---------------------------------------------------------------------------- static bool RefineIntersNewton( const Point3d& ptL, const Vector3d& vtL, double dLen, bool bFinite, const ISurfBezier* pSurfBz, Point3d& ptSP, Point3d& ptIBz) { // la funzione raffina la posisione del punto ptSP, minimizzando la distanza dalla retta e restituisce il punto di intersezione ptIBz // usando un algoritmo di newton cerco di avvicinarmi il più possibile alla retta DistPointLine dpl( ptIBz, ptL, vtL, dLen, bFinite) ; double dDistNew = 0, dDistPre = 0 ; dpl.GetDist(dDistNew) ; int nCount = 0 ; double dh = EPS_SMALL ; pSurfBz->GetPointD1D2( ptSP.x, ptSP.y, ISurfBezier::FROM_MINUS, ISurfBezier::FROM_MINUS, ptIBz) ; // metodo di newton in più dimensioni // vario sia il parametro U che il parametro V e verifico se la distanza dalla retta diminuisce per scostamenti positivi o negativi. while ( dDistNew > EPS_SMALL && nCount < 100) { dDistPre = dDistNew ; Point3d ptIBzNew1 ; pSurfBz->GetPointD1D2( ( ptSP.x + dh), ptSP.y, ISurfBezier::FROM_MINUS, ISurfBezier::FROM_MINUS, ptIBzNew1) ; DistPointLine dplNewU( ptIBzNew1, ptL, vtL, dLen, bFinite) ; dplNewU.GetDist( dDistNew) ; double dfdU = ( dDistNew - dDistPre) / dh ; Point3d ptIBzNew2 ; pSurfBz->GetPointD1D2( ptSP.x, ( ptSP.y + dh), ISurfBezier::FROM_MINUS, ISurfBezier::FROM_MINUS, ptIBzNew2) ; DistPointLine dplNewV( ptIBzNew2, ptL, vtL, dLen, bFinite) ; dplNewV.GetDist( dDistNew) ; double dfdV = ( dDistNew - dDistPre) / dh ; // mi avvicino cercando di annullare la distanza in un colpo solo double dr = - dDistPre / ( dfdU + dfdV) ; pSurfBz->GetPointD1D2(( ptSP.x + dr * dfdU), ( ptSP.y + dr * dfdV), ISurfBezier::FROM_MINUS, ISurfBezier::FROM_MINUS, ptIBz) ; DistPointLine dplNew( ptIBz, ptL, vtL, dLen, bFinite) ; dplNew.GetDist( dDistNew) ; ++ nCount ; } return ( nCount != 99) ; } //---------------------------------------------------------------------------- static void UpdateInfoIntersLineSurfBz( const Point3d& ptL, const Vector3d& vtDir, int nILT, int nT, const Point3d& ptSP, const Point3d& ptIBz, double dCos, const Point3d& ptSP2, const Point3d& ptIBz2, double dCos2, ILSBIVECTOR& vInfo) { int nType = LSBT_NONE ; if ( dCos > EPS_ZERO) nType = LSBT_IN ; else if ( dCos < EPS_ZERO) nType = LSBT_OUT ; else nType = LSBT_TOUCH ; if ( nILT == ILTA_IN || nILT == ILTA_EDGE || nILT == ILTA_VERT || nILT == ILTA_NO_TRIA) { double dU = ( ptIBz - ptL) * vtDir ; vInfo.emplace_back( nType, dU, nT, dCos, ptIBz, ptSP) ; } else if ( nILT == ILTA_SEGM || nILT == ILTA_SEGM_ON_EDGE) { double dU = ( ptIBz - ptL) * vtDir ; double dU2 = ( ptIBz2 - ptL) * vtDir ; int nType2 = LSBT_NONE ; if ( dCos2 > EPS_ZERO) nType2 = LSBT_IN ; else if ( dCos2 < EPS_ZERO) nType2 = LSBT_OUT ; vInfo.emplace_back( nType, dU, 0, nT, dCos, ptIBz, P_INVALID, ptSP, P_INVALID) ; vInfo.emplace_back( nType2, dU2, 0, nT, dCos2, ptIBz2, P_INVALID, ptSP2, P_INVALID) ; } } //---------------------------------------------------------------------------- static void OrderInfoIntersLineSurfBz( ILSBIVECTOR& vInfo) { // se non trovati, esco if ( vInfo.empty()) return ; // ordino il vettore delle intersezioni secondo il senso crescente del parametro di linea sort( vInfo.begin(), vInfo.end(), []( const IntLinSbzInfo& a, const IntLinSbzInfo& b) { return ( a.dU < b.dU) ; }) ; } //---------------------------------------------------------------------------- // Intersezione di una linea con una superficie di Bezier //---------------------------------------------------------------------------- bool IntersLineSurfBz( const Point3d& ptL, const Vector3d& vtL, double dLen, const ISurfBezier* pSurfBz, ILSBIVECTOR& vInfo, bool bFinite) { PtrOwner pCL( CreateCurveLine()) ; if ( bFinite) pCL->SetPVL(ptL, vtL, dLen) ; else pCL->SetPVL(ptL, vtL, 1e6) ; // verifico linea Vector3d vtDir = vtL ; if ( ! vtDir.Normalize( EPS_ZERO)) return false ; // verifico superficie if ( pSurfBz == nullptr) return false ; // verifico parametro di ritorno if ( &vInfo == nullptr) return false ; vInfo.clear() ; // trovo le intersezioni con la trimesh ausiliaria const ISurfTriMesh* pSurfTm = pSurfBz->GetAuxSurf() ; ILSIVECTOR vInfoTm ; if ( ! IntersLineSurfTm( ptL, vtL, dLen, *pSurfTm, vInfoTm, bFinite)) return false ; // ricavo le intersezioni con la superficie di Bezier for ( IntLinStmInfo InfoTm : vInfoTm ) { // devo raffinare i parametri lungo la curva, l'angolo e i punti di intersezione Point3d ptI, ptI2 ; // devo trovare le intersezioni Point3d ptSP, ptSP2 ; // coordinate parametriche delle soluzioni pSurfBz->UnprojectPointFromStm( InfoTm.nT, InfoTm.ptI, ptSP, InfoTm.nILTT) ; Point3d ptIBz, ptIBz2 ; if ( ! RefineIntersNewton( ptL, vtL, dLen, bFinite, pSurfBz, ptSP, ptIBz)) { /////// posso provare anche a rilanciare newton con un punto di partenza diverso oppure con una direzione di avvicinamento diversa/////////////////////////////////// // per restare nel triangolo mi sposto verso un vertice int nVert[3] ; pSurfTm->GetTriangle( InfoTm.nT, nVert) ; double dU0, dV0 ; pSurfTm->GetVertexParam( nVert[0], dU0, dV0) ; ptSP = ptSP + Point3d( dU0, dV0, 0) ; if ( ! RefineIntersNewton( ptL,vtL, dLen, bFinite, pSurfBz, ptSP, ptIBz)) return false ; } Vector3d vtN ; pSurfBz->GetPointNrmD1D2(ptSP.x, ptSP.y, ISurfBezier::FROM_MINUS, ISurfBezier::FROM_MINUS, ptIBz, vtN) ; double dCos = vtN * vtL ; double dCos2 = 0 ; // eventualmente ripeto tutto per ptI2 ( se ho un'intersezione con sovrapposizione) if ( InfoTm.nILTT == ILTA_SEGM || InfoTm.nILTT == ILTA_SEGM_ON_EDGE ) { pSurfBz->UnprojectPointFromStm( InfoTm.nT, InfoTm.ptI2, ptSP2, InfoTm.nILTT) ; if ( ! RefineIntersNewton(ptL, vtL, dLen, bFinite, pSurfBz, ptSP2, ptIBz2) ) { int nVert[3] ; pSurfTm->GetTriangle( InfoTm.nT, nVert) ; double dU0, dV0 ; pSurfTm->GetVertexParam( nVert[0], dU0, dV0) ; ptSP = ptSP + Point3d(dU0, dV0, 0) ; if ( ! RefineIntersNewton( ptL,vtL, dLen, bFinite, pSurfBz, ptSP, ptIBz)) return false ; } pSurfBz->GetPointNrmD1D2( ptSP2.x, ptSP2.y, ISurfBezier::FROM_MINUS, ISurfBezier::FROM_MINUS, ptIBz2, vtN) ; dCos2 = vtN * vtL ; } UpdateInfoIntersLineSurfBz( ptL, vtL, InfoTm.nILTT, InfoTm.nT, ptSP, ptIBz, dCos, ptSP2, ptIBz2, dCos2, vInfo) ; } OrderInfoIntersLineSurfBz( vInfo) ; return true ; } //---------------------------------------------------------------------------- bool FilterLineSurfBzInters( const ILSBIVECTOR& vInfo, INTDBLVECTOR& vInters) { // tengo per buone la classificazione delle intersezioni fatte sulla trimesh // ciclo sulle intersezioni for ( const auto& Info : vInfo) { // se intersezione puntuale vInters.emplace_back( Info.nILSB, Info.dU) ; // se intersezione sovrapposta // da sviluppare } // elimino intersezioni ripetute for ( size_t j = 1 ; j < vInters.size() ; ) { // intersezione precedente size_t i = j - 1 ; // se hanno lo stesso parametro if ( abs( vInters[i].second - vInters[j].second) < EPS_SMALL) { // se sono entrambe entranti o uscenti, elimino la seconda if ( ( vInters[i].first == LSBT_IN && vInters[j].first == LSBT_IN) || ( vInters[i].first == LSBT_OUT && vInters[j].first == LSBT_OUT)) { vInters.erase( vInters.begin() + j) ; continue ; } // se una entrante e l'altra uscente, cambio in touch ed elimino la seconda else if ( ( vInters[i].first == LSBT_IN && vInters[j].first == LSBT_OUT) || ( vInters[i].first == LSBT_OUT && vInters[j].first == LSBT_IN)) { vInters[i].first = LSBT_TOUCH ; vInters.erase( vInters.begin() + j) ; continue ; } // se una puntuale e l'altra inizio di coincidenza, elimino la prima else if ( ( vInters[i].first == LSBT_IN || vInters[i].first == LSBT_OUT || vInters[i].first == LSBT_TOUCH) && vInters[j].first == LSBT_TG_INI) { vInters.erase( vInters.begin() + i) ; continue ; } // se una fine di coincidenza e l'altra puntuale, elimino la seconda else if ( vInters[i].first == LSBT_TG_FIN && ( vInters[j].first == LSBT_IN || vInters[j].first == LSBT_OUT || vInters[j].first == LSBT_TOUCH)) { vInters.erase( vInters.begin() + j) ; continue ; } // se una fine di coincidenza e l'altra inizio di coincidenza, elimino entrambe else if ( i > 0 && vInters[i].first == LSBT_TG_FIN && vInters[j].first == LSBT_TG_INI) { vInters.erase( vInters.begin() + j) ; vInters.erase( vInters.begin() + i) ; -- j ; continue ; } } // passo alla successiva ++ j ; } return true ; } //---------------------------------------------------------------------------- // Intersezione di una linea con una superficie di Bezier di grado 3x1 monopatch //---------------------------------------------------------------------------- bool IntersLineSurfBzCubicLinear( const Point3d& ptL, const Vector3d& vtL, double dLen, const ISurfBezier* pSurfBz, ILSBIVECTOR& vInfo, bool bFinite) { int nDegU, nDegV, nSpanU, nSpanV ; bool bRat, bTrimmed ; pSurfBz->GetInfo( nDegU, nDegV, nSpanU, nSpanV, bRat, bTrimmed) ; // funzione pensata per funzionare solo con una monopatch di grado 3x1 if ( nDegU != 3 || nDegV != 1 || nSpanU > 1 || nSpanV > 1 || bRat) return false ; int nInters = int( vInfo.size()) ; Point3d r = ptL ; Vector3d q = vtL ; bool bNeedToRotX = AreSameVectorApprox( q, X_AX) ; bool bNeedToRotY = AreSameVectorApprox( q, Y_AX) ; bool bNeedToRot = bNeedToRotX || bNeedToRotY ; Frame3d frRot ; if ( bNeedToRotX) frRot.Set( ORIG, X_AX) ; if ( bNeedToRotY) frRot.Set( ORIG, Y_AX) ; if ( bNeedToRot) { r.ToLoc( frRot) ; q.ToLoc( frRot) ; } PNTVECTOR vPntCtrl = pSurfBz->GetAllControlPoints() ; if ( bNeedToRot) { for ( Point3d& pt: vPntCtrl) pt.ToLoc( frRot) ; } Vector3d A = vPntCtrl[4] - vPntCtrl[0] ; Vector3d B = vPntCtrl[5] - vPntCtrl[1] ; Vector3d C = vPntCtrl[6] - vPntCtrl[2] ; Vector3d D = vPntCtrl[7] - vPntCtrl[3] ; Vector3d E = vPntCtrl[0] - ORIG ; Vector3d F = vPntCtrl[1] - ORIG ; Vector3d G = vPntCtrl[2] - ORIG ; Vector3d H = vPntCtrl[3] - ORIG ; Vector3d a3 = -A + 3 * B - 3 * C + D ; Vector3d a2 = 3 * A - 6 * B + 3 * C ; Vector3d a1 = -3 * A + 3 * B ; Vector3d a0 = A ; Vector3d b3 = -E + 3 * F - 3 * G + H ; Vector3d b2 = 3 * E - 6 * F + 3 * G ; Vector3d b1 = -3 * E + 3 * F ; Vector3d b0 = E ; DBLVECTOR vdCoeff, vdRoots ; // coefficienti dal grado più basso al grado più alto vdCoeff = { // c0 q.x*q.z*a0.y*b0.z - q.x*q.y*a0.z*b0.z // 3 - r.z*q.x*q.z*a0.y + r.z*q.x*q.y*a0.z + // 3 q.y*q.z*a0.z*b0.z - q.z*q.z*a0.y*b0.x // 4 - r.x*q.y*q.z*a0.z + r.x*q.z*q.z*a0.y + // 4 q.z*q.z*a0.x*b0.y - q.y*q.z*a0.x*b0.z - q.x*q.z*a0.z*b0.y + q.x*q.y*a0.z*b0.z // 5 - r.y*q.z*q.z*a0.x + r.z*q.y*q.z*a0.x + r.y*q.x*q.z*a0.z - r.z*q.x*q.y*a0.z, // 5 // c1 q.x*q.z*(a1.y*b0.z + a0.y*b1.z) - q.x*q.y*(a1.z*b0.z + a0.z*b1.z) // 3 - r.z*q.x*q.z*a1.y + r.z*q.x*q.y*a1.z + // 3 q.y*q.z*(a1.z*b0.x + a0.z*b1.x) - q.z*q.z*(a1.y*b0.x + a0.y*b1.x) // 4 - r.x*q.y*q.z*a1.z + r.x*q.z*q.z*a1.y + // 4 q.z*q.z*(a1.x*b0.y + a0.x*b1.y) - q.y*q.z*(a1.x*b0.z + a0.x*b1.z) // 5 - q.x*q.z*(a1.z*b0.y + a0.z*b1.y) + q.x*q.y*(a1.z*b0.z + a0.z*b1.z) // 5 - r.y*q.z*q.z*a1.x + r.z*q.y*q.z*a1.x + r.y*q.x*q.z*a1.z - r.z*q.x*q.y*a1.z, // 5 // c2 q.x*q.z*(a2.y*b0.z + a1.y*b1.z + a0.y*b2.z) - q.x*q.y*(a2.z*b0.z + a1.z*b1.z + a0.z*b2.z) // 3 - r.z*q.x*q.z*a2.y + r.z*q.x*q.y*a2.z + // 3 q.y*q.z*(a2.z*b0.x + a1.z*b1.x + a0.z*b2.x) - q.z*q.z*(a2.y*b0.x + a1.y*b1.x + a0.y*b2.x) // 4 - r.x*q.y*q.z*a2.z + r.x*q.z*q.z*a2.y + // 4 q.z*q.z*(a2.x*b0.y + a1.x*b1.y + a0.x*b2.y) - q.y*q.z*(a2.x*b0.z + a1.x*b1.z + a0.x*b2.z) // 5 - q.x*q.z*(a2.z*b0.y + a1.z*b1.y + a0.z*b2.y) + q.x*q.y*(a2.z*b0.z + a1.z*b1.z + a0.z*b2.z)// 5 - r.y*q.z*q.z*a2.x + r.z*q.y*q.z*a2.x + r.y*q.x*q.z*a2.z - r.z*q.x*q.y*a2.z, // 5 // c3 q.x*q.z*(a3.y*b0.z + a2.y*b1.z + a1.y*b2.z + a0.y*b3.z) - q.x*q.y*(a3.z*b0.z + a2.z*b1.z + a1.z*b2.z + a0.z*b3.z) // 3 - r.z*q.x*q.z*a3.y + r.z*q.x*q.y*a3.z + // 3 q.y*q.z*(a3.z*b0.x + a2.z*b1.x + a1.z*b2.x + a0.z*b3.x) - q.z*q.z*(a3.y*b0.x + a2.y*b1.x + a1.y*b2.x + a0.y*b3.x) // 4 - r.x*q.y*q.z*a3.z + r.x*q.z*q.z*a3.y + // 4 q.z*q.z*(a3.x*b0.y + a2.x*b1.y + a1.x*b2.y + a0.x*b3.y) - q.y*q.z*(a3.x*b0.z + a2.x*b1.z + a1.x*b2.z + a0.x*b3.z) // 5 - q.x*q.z*(a3.z*b0.y + a2.z*b1.y + a1.z*b2.y + a0.z*b3.y) + q.x*q.y*(a3.z*b0.z + a2.z*b1.z + a1.z*b2.z + a0.z*b3.z)// 5 - r.y*q.z*q.z*a3.x + r.z*q.y*q.z*a3.x + r.y*q.x*q.z*a3.z - r.z*q.x*q.y*a3.z, // 5 // c4 q.x*q.z*(a3.y*b1.z + a2.y*b2.z + a1.y*b3.z) - q.x*q.y*(a3.z*b1.z + a2.z*b2.z + a1.z*b3.z) + // 3 q.y*q.z*(a3.z*b1.x + a2.z*b2.x + a1.z*b3.x) - q.z*q.z*(a3.y*b1.x + a2.y*b2.x + a1.y*b3.x) + // 4 q.z*q.z*(a3.x*b1.y + a2.x*b2.y + a1.x*b3.y) - q.y*q.z*(a3.x*b1.z + a2.x*b2.z + a1.x*b3.z) // 5 - q.x*q.z*(a3.z*b1.y + a2.z*b2.y + a1.z*b3.y) + q.x*q.y*(a3.z*b1.z + a2.z*b2.z + a1.z*b3.z), // 5 // c5 q.x*q.z*(a3.y*b2.z + a2.y*b3.z) - q.x*q.y*(a3.z*b2.z + a2.z*b3.z) + // 3 q.y*q.z*(a3.z*b2.x + a2.z*b3.x) - q.z*q.z*(a3.y*b2.x + a2.y*b3.x) + // 4 q.z*q.z*(a3.x*b2.y + a2.x*b3.y) - q.y*q.z*(a3.x*b2.z + a2.x*b3.z) // 5 - q.x*q.z*(a3.z*b2.y + a2.z*b3.y) + q.x*q.y*(a3.z*b2.z + a2.z*b3.z), // 5 // c6 q.x*q.z*a3.y*b3.z - q.x*q.y*a3.z*b3.z + // 3 q.y*q.z*a3.z*b3.x - q.z*q.z*a3.y*b3.x + // 4 q.z*q.z*a3.x*b3.y - q.y*q.z*a3.x*b3.z - q.x*q.z*a3.z*b3.y + q.x*q.y*a3.z*b3.z} ; // 5 int nRoots = PolynomialRoots( 6, vdCoeff, vdRoots) ; bool bFound = false ; for ( int w = 0 ; w < nRoots ; ++w) { double dU = 0, dV = 0 ; if ( vdRoots[w] > 0 - EPS_ZERO && vdRoots[w] < 1 + EPS_ZERO) { dU = vdRoots[w] ; // verifico che non sia una soluzione con molteplicità > 1 bool bAlreadyFound = false ; for ( int k = w - 1 ; k >= 0 && ! bAlreadyFound ; --k) bAlreadyFound = ( abs( dU - vdRoots[k]) < EPS_PARAM) ; if ( ! bAlreadyFound) { Vector3d vAlpha = a3 * pow(dU, 3) + a2 * pow( dU, 2) + a1 * dU + a0 ; Vector3d vBeta = b3 * pow(dU, 3) + b2 * pow( dU, 2) + b1 * dU + b0 ; double dDen = ( vAlpha.x * q.z - vAlpha.z * q.x) ; if ( abs( dDen) > EPS_ZERO) dV = ( ( vBeta.z - r.z) * q.x - ( vBeta.x - r.x ) * q.z) / dDen ; else { // se la prima equazione risulta un x/0 allora uso la seconda equazione per trovare il secondo parametro double dDen2 = ( vAlpha.y * q.z - vAlpha.z * q.y) ; dV = ( ( vBeta.z - r.z) * q.y - ( vBeta.y - r.y ) * q.z) / dDen2 ; } if ( dV > - EPS_ZERO && dV < 1 + EPS_ZERO) { Point3d ptIBez, ptIBez2 ; Vector3d vtN ; pSurfBz->GetPointNrmD1D2(dU, dV, ISurfBezier::Side::FROM_MINUS, ISurfBezier::Side::FROM_MINUS, ptIBez, vtN) ; Point3d ptSP( dU, dV, 0), ptSP2 ; double dCos = vtN * vtL, dCos2 = 0 ; int nType = ILTA_NO_TRIA ; UpdateInfoIntersLineSurfBz( ptL, vtL, nType, -1, ptSP, ptIBez, dCos, ptSP2, ptIBez2, dCos2, vInfo) ; bFound = true ; } } } } //// se tutti i coefficienti sono zero allora potrei avere una linea che giace sulla superficie //// per trovare i punti di inizio e fine sovrapposizione trovo i punti a minima distanza tra la linea e gli edge della superficie //if ( ! bFound && abs( vdCoeff[0]) < EPS_ZERO && abs( vdCoeff[1]) < EPS_ZERO && abs( vdCoeff[2]) < EPS_ZERO) { // ICRVCOMPOPOVECTOR vCrvEdge( 4) ; // vCrvEdge[0].Set(pSurfBz->GetCurveOnU( 0)) ; // vCrvEdge[1].Set(pSurfBz->GetCurveOnV( 1)) ; // vCrvEdge[2].Set(pSurfBz->GetCurveOnU( 1)) ; // vCrvEdge[3].Set(pSurfBz->GetCurveOnV( 0)) ; // double dAngTolDeg = 5 ; // for ( int i = 0 ; i < 4 ; ++i) { // PolyLine plApprox ; vCrvEdge[0]->ApproxWithLines( EPS_SMALL, dAngTolDeg, ICurve::ApprLineType::APL_STD, plApprox) ; // //CurveComposite cCC ; // //cCC.FromPolyLine( plApprox) ; // int nClosestLine = -1 ; // double dMinDist = INFINITO ; // Point3d pt ; plApprox.GetFirstPoint( pt) ; // Point3d ptClosest ; // int c = 0 ; // int nTot = plApprox.GetPointNbr() ; // for ( int j = 0 ; j < nTot ; ++j) { // DistPointLine dpl( pt, ptL, vtL, dLen, bFinite) ; // double dDist = INFINITO ; // dpl.GetDist( dDist) ; // if ( dDist < dMinDist) { // nClosestLine = c ; // dMinDist = dDist ; // } // plApprox.GetNextPoint( pt) ; // ++ c ; // } // Point3d ptInt1, ptInt2 ; // if ( nClosestLine < nTot - 1 && nClosestLine > 0) { // // tra i due tratti dell'approssimazione che arrivano al punto selezionato come più vicino, devo trovare quale si avvicina di più // Point3d ptStart ; plApprox.GetFirstPoint( ptStart) ; // Point3d ptEnd ; // for ( int z = 1 ; z < nClosestLine - 1 ; ++z) // plApprox.GetNextPoint( ptStart) ; // plApprox.GetNextPoint( ptEnd) ; // // linea precedente al punto // Vector3d vtLinePre = ptEnd - ptStart ; // double dLenPre = vtLinePre.Len() ; // DistLineLine dllPre( ptStart, vtLinePre, dLenPre, ptL, vtL,dLen) ; // double dDistPre = INFINITO ; // dllPre.GetDist( dDistPre) ; // // linea che inzia con quel punto // ptStart = ptEnd ; // plApprox.GetNextPoint( ptEnd) ; // Vector3d vtLineCurr = ptEnd - ptStart ; // double dLenCurr = vtLineCurr.Len() ; // DistLineLine dllCurr( ptStart, vtLineCurr, dLenCurr, ptL, vtL,dLen) ; // double dDistCurr = INFINITO ; // dllCurr.GetDist( dDistCurr) ; // if ( dDistPre < dDistCurr) // dllPre.GetMinDistPoints( ptInt1, ptInt2) ; // else // dllCurr.GetMinDistPoints( ptInt1, ptInt2) ; // } // else if ( nClosestLine == 0) { // // il punto più vicino è sulla prima linea // Point3d ptStart ; plApprox.GetFirstPoint( ptStart) ; // Point3d ptEnd ; plApprox.GetNextPoint( ptEnd) ; // Vector3d vtLineCurr = ptEnd - ptStart ; // double dLenCurr = vtLineCurr.Len() ; // DistLineLine dllCurr( ptStart, vtLineCurr, dLenCurr, ptL, vtL,dLen) ; // dllCurr.GetMinDistPoints( ptInt1, ptInt2) ; // } // else if ( nClosestLine == nTot- 1) { // // il punto più vicino è sull'ultima linea // Point3d ptStart ; plApprox.GetFirstPoint( ptStart) ; // Point3d ptEnd ; // for ( int z = 1 ; z < nClosestLine - 1 ; ++z) // plApprox.GetNextPoint( ptStart) ; // plApprox.GetNextPoint( ptEnd) ; // Vector3d vtLinePre = ptEnd - ptStart ; // double dLenPre = vtLinePre.Len() ; // DistLineLine dllCurr( ptStart, vtLinePre, dLenPre, ptL, vtL,dLen) ; // dllCurr.GetMinDistPoints( ptInt1, ptInt2) ; // } // // double dU1 = 0, dV1 = 0, dU2 = 0, dV2 = 0 ; // // se ho trovato due punti vuol dire che la linea coincide con un edge e ho trovato tutto quello che serve // if ( ! AreSamePointExact( ptInt2, ORIG)) { // if ( i == 0) { // //dV1 = 0 ; dV2 = 0 ; // vCrvEdge[0]->GetParamAtPoint( ptInt1, dU1) ; // vCrvEdge[0]->GetParamAtPoint( ptInt2, dU2) ; // } // else if ( i == 1) { // //dU1 = 1 ; dU2 = 1 ; // vCrvEdge[1]->GetParamAtPoint( ptInt1, dV1) ; // vCrvEdge[1]->GetParamAtPoint( ptInt2, dV2) ; // } // else if ( i == 2){ // //dV1 = 1 ; dV2 = 1 ; // vCrvEdge[2]->GetParamAtPoint( ptInt1, dU1) ; // vCrvEdge[2]->GetParamAtPoint( ptInt2, dU2) ; // } // else if ( i == 3){ // //dU1 = 0 ; dU2 = 0 ; // vCrvEdge[3]->GetParamAtPoint( ptInt1, dV1) ; // vCrvEdge[3]->GetParamAtPoint( ptInt2, dV2) ; // } // Point3d ptIBez1, ptIBez2 ; // Vector3d vtN1, vtN2 ; // pSurfBz->GetPointNrmD1D2(dU1, dV1, ISurfBezier::Side::FROM_MINUS, ISurfBezier::Side::FROM_MINUS, ptIBez1, vtN1) ; // pSurfBz->GetPointNrmD1D2(dU2, dV2, ISurfBezier::Side::FROM_MINUS, ISurfBezier::Side::FROM_MINUS, ptIBez2, vtN2) ; // Point3d ptSP1( dU1, dV1, 0) ; // double dCos1 = vtN1 * vtL ; // Point3d ptSP2( dU2, dV2, 0) ; // double dCos2 = vtN2 * vtL ; // // se avevo già trovato un punto singolo che coincide col primo punto di questa intersezione sovrapposta, allora cancello l'intersezione singola che // // avevo salvato e aggiungo quella sovrapposto che ho trovato ora // if ( bFound) { // int nNewTot = int(vInfo.size()) ; // int nNewInters = nNewTot - nInters ; // bool bAlreadyFound = false ; // for ( int i = 0 ; i < nNewInters ; ++i) { // bAlreadyFound = AreSamePointApprox(vInfo[nNewTot - i].ptUV, ptSP1) || AreSamePointApprox(vInfo[nNewTot - i].ptUV, ptSP2) ; // if ( bAlreadyFound) { // vInfo.erase( vInfo.begin() + nNewTot - i) ; // break ; // } // } // } // UpdateInfoIntersLineSurfBz( ptL, vtL, ILTA_NO_TRIA, -1, ptSP1, ptIBez1, dCos1, ptSP2, ptIBez2, dCos2, vInfo) ; // bFound = true ; // break ; // } // // se ho trovato un punto a distanza zero dalla linea allora ho trovato l'intersezione // else if ( dMinDist < EPS_SMALL) { // if ( i == 0) { // //dV1 = 0 ; // vCrvEdge[0]->GetParamAtPoint( ptInt1, dU1) ; // } // else if ( i == 1) { // //dU1 = 1 ; // vCrvEdge[1]->GetParamAtPoint( ptInt1, dV1) ; // } // else if ( i == 2) { // //dV1 = 1 ; // vCrvEdge[2]->GetParamAtPoint( ptInt1, dU1) ; // } // else if ( i == 3) { // //dU1 = 0 ; // vCrvEdge[3]->GetParamAtPoint( ptInt1, dV1) ; // } // Point3d ptSP1( dU1, dV1, 0), ptSP2 ; // // se avevo trovato già altri punti controllo di non essere esattamente su una diagonale ( e quindi avere un'intersezione con ogni edge, ma due sono doppie) // if ( bFound) { // int nNewTot = int(vInfo.size()) ; // int nNewInters = nNewTot - nInters ; // bool bAlreadyFound = false ; // for ( int i = 0 ; i < nNewInters ; ++i) // bAlreadyFound = AreSamePointApprox(vInfo[nNewTot - i].ptUV, ptSP1) ; // if ( bAlreadyFound) // continue ; // } // Point3d ptIBez1, ptIBez2 ; // Vector3d vtN1, vtN2 ; // pSurfBz->GetPointNrmD1D2(dU1, dV1, ISurfBezier::Side::FROM_MINUS, ISurfBezier::Side::FROM_MINUS, ptIBez1, vtN1) ; // double dCos1 = vtN1 * vtL, dCos2 = 0 ; // UpdateInfoIntersLineSurfBz( ptL, vtL, ILTA_NO_TRIA, -1, ptSP1, ptIBez1, dCos1, ptSP2, ptIBez2, dCos2, vInfo) ; // bFound = true ; // } // } //} // se la superficie è trimmed verifico che i punti trovati siano all'interno del parametrico trimmato if ( bTrimmed && bFound) { int nNewTot = int(vInfo.size()) ; int nNewInters = nNewTot - nInters ; const ISurfFlatRegion* pFRTrim = pSurfBz->GetTrimRegion() ; for ( int i = 0 ; i < nNewInters ; ++i) { Point3d ptTest = vInfo[nNewTot - i].ptUV * SBZ_TREG_COEFF ; bool bInside = false ; double dDist = INFINITO ; IsPointInsideSurfFr( ptTest, pFRTrim, dDist, bInside) ; if ( ! bInside) vInfo.erase( vInfo.begin() + nNewTot - i) ; } } return true ; } //---------------------------------------------------------------------------- // Intersezione di una linea con una superficie di Bezier bilineare monopatch //---------------------------------------------------------------------------- bool IntersLineSurfBzBilinear( const Point3d& ptL, const Vector3d& vtL, double dLen, const ISurfBezier* pSurfBz, ILSBIVECTOR& vInfo, bool bFinite) { int nDegU, nDegV, nSpanU, nSpanV ; bool bRat, bTrimmed ; pSurfBz->GetInfo( nDegU, nDegV, nSpanU, nSpanV, bRat, bTrimmed) ; // funzione pensata per funzionare solo con una monopatch bilineare if ( nDegU > 1 || nDegV > 1 || nSpanU > 1 || nSpanV > 1 || bRat) return false ; int nInters = int( vInfo.size()) ; PNTVECTOR vPntCtrl ; for ( int p = 0 ; p < 4 ; ++p) { bool bOk = false ; vPntCtrl.push_back( pSurfBz->GetControlPoint( p, &bOk)) ; } Vector3d a = vPntCtrl[3] - vPntCtrl[1] + ( vPntCtrl[0] - vPntCtrl[2]) ; Vector3d b = vPntCtrl[1] - vPntCtrl[0] ; Vector3d c = vPntCtrl[2] - vPntCtrl[0] ; Vector3d d = vPntCtrl[0] - ORIG ; double A1 = a.x * vtL.z - a.z * vtL.x ; double B1 = b.x * vtL.z - b.z * vtL.x ; double C1 = c.x * vtL.z - c.z * vtL.x ; double A2 = a.y * vtL.z - a.z * vtL.y ; double B2 = b.y * vtL.z - b.z * vtL.y ; double C2 = c.y * vtL.z - c.z * vtL.y ; double D1 = ( d.x - ptL.x) * vtL.z - ( d.z - ptL.z) * vtL.x ; double D2 = ( d.y - ptL.y) * vtL.z - ( d.z - ptL.z) * vtL.y ; DBLVECTOR vdCoeff, vdRoots ; vdCoeff = { (B2 * D1 - B1 * D2), ( A2 * D1 - A1 * D2 + B2 * C1 - B1 * C2), ( A2 * C1 - A1 * C2)} ; int nRoots = PolynomialRoots( 2, vdCoeff, vdRoots) ; bool bFound = false ; for ( int w = 0 ; w < nRoots ; ++w) { if ( vdRoots[w] > 0 - EPS_ZERO && vdRoots[w] < 1 + EPS_ZERO ) { double dU = 0, dV = vdRoots[w] ; // verifico che non sia una soluzione con molteplicità > 1 bool bAlreadyFound = false ; for ( int k = w - 1 ; k >= 0 && ! bAlreadyFound ; --k) bAlreadyFound = abs( dV - vdRoots[k]) < EPS_PARAM ; if ( ! bAlreadyFound) { dU = (dV * (C1 - C2) + ( D1 - D2)) / ( dV * ( A2 - A1) + ( B2 - B1)) ; if ( dU > - EPS_ZERO && dU < 1 + EPS_ZERO) { Point3d ptIBez, ptIBez2 ; Vector3d vtN ; pSurfBz->GetPointNrmD1D2(dU, dV, ISurfBezier::Side::FROM_MINUS, ISurfBezier::Side::FROM_MINUS, ptIBez, vtN) ; Point3d ptSP( dU, dV, 0), ptSP2 ; double dCos = vtN * vtL, dCos2 = 0 ; UpdateInfoIntersLineSurfBz( ptL, vtL, ILTA_NO_TRIA, -1, ptSP, ptIBez, dCos, ptSP2, ptIBez2, dCos2, vInfo) ; bFound = true ; } } } } // se tutti i coefficienti sono zero allora potrei avere una linea che giace sulla superficie // per trovare i punti di inizio e fine sovrapposizione trovo i punti a minima distanza tra la linea e gli edge della superficie if ( ! bFound && abs( vdCoeff[0]) < EPS_ZERO && abs( vdCoeff[1]) < EPS_ZERO && abs( vdCoeff[2]) < EPS_ZERO) { ICRVCOMPOPOVECTOR vCrvEdge( 4) ; vCrvEdge[0].Set(pSurfBz->GetCurveOnU( 0)) ; vCrvEdge[1].Set(pSurfBz->GetCurveOnV( 1)) ; vCrvEdge[2].Set(pSurfBz->GetCurveOnU( 1)) ; vCrvEdge[3].Set(pSurfBz->GetCurveOnV( 0)) ; double dAngTolDeg = 5 ; for ( int i = 0 ; i < 4 ; ++i) { PolyLine plApprox ; vCrvEdge[0]->ApproxWithLines( EPS_SMALL, dAngTolDeg, ICurve::ApprLineType::APL_STD, plApprox) ; //CurveComposite cCC ; //cCC.FromPolyLine( plApprox) ; int nClosestLine = -1 ; double dMinDist = INFINITO ; Point3d pt ; plApprox.GetFirstPoint( pt) ; Point3d ptClosest ; int c = 0 ; int nTot = plApprox.GetPointNbr() ; for ( int j = 0 ; j < nTot ; ++j) { DistPointLine dpl( pt, ptL, vtL, dLen, bFinite) ; double dDist = INFINITO ; dpl.GetDist( dDist) ; if ( dDist < dMinDist) { nClosestLine = c ; dMinDist = dDist ; } plApprox.GetNextPoint( pt) ; ++ c ; } Point3d ptInt1, ptInt2 ; if ( nClosestLine < nTot - 1 && nClosestLine > 0) { // tra i due tratti dell'approssimazione che arrivano al punto selezionato come più vicino, devo trovare quale si avvicina di più Point3d ptStart ; plApprox.GetFirstPoint( ptStart) ; Point3d ptEnd ; for ( int z = 1 ; z < nClosestLine - 1 ; ++z) plApprox.GetNextPoint( ptStart) ; plApprox.GetNextPoint( ptEnd) ; // linea precedente al punto Vector3d vtLinePre = ptEnd - ptStart ; double dLenPre = vtLinePre.Len() ; DistLineLine dllPre( ptStart, vtLinePre, dLenPre, ptL, vtL,dLen) ; double dDistPre = INFINITO ; dllPre.GetDist( dDistPre) ; // linea che inzia con quel punto ptStart = ptEnd ; plApprox.GetNextPoint( ptEnd) ; Vector3d vtLineCurr = ptEnd - ptStart ; double dLenCurr = vtLineCurr.Len() ; DistLineLine dllCurr( ptStart, vtLineCurr, dLenCurr, ptL, vtL,dLen) ; double dDistCurr = INFINITO ; dllCurr.GetDist( dDistCurr) ; if ( dDistPre < dDistCurr) dllPre.GetMinDistPoints( ptInt1, ptInt2) ; else dllCurr.GetMinDistPoints( ptInt1, ptInt2) ; } else if ( nClosestLine == 0) { // il punto più vicino è sulla prima linea Point3d ptStart ; plApprox.GetFirstPoint( ptStart) ; Point3d ptEnd ; plApprox.GetNextPoint( ptEnd) ; Vector3d vtLineCurr = ptEnd - ptStart ; double dLenCurr = vtLineCurr.Len() ; DistLineLine dllCurr( ptStart, vtLineCurr, dLenCurr, ptL, vtL,dLen) ; dllCurr.GetMinDistPoints( ptInt1, ptInt2) ; } else if ( nClosestLine == nTot- 1) { // il punto più vicino è sull'ultima linea Point3d ptStart ; plApprox.GetFirstPoint( ptStart) ; Point3d ptEnd ; for ( int z = 1 ; z < nClosestLine - 1 ; ++z) plApprox.GetNextPoint( ptStart) ; plApprox.GetNextPoint( ptEnd) ; Vector3d vtLinePre = ptEnd - ptStart ; double dLenPre = vtLinePre.Len() ; DistLineLine dllCurr( ptStart, vtLinePre, dLenPre, ptL, vtL,dLen) ; dllCurr.GetMinDistPoints( ptInt1, ptInt2) ; } double dU1 = 0, dV1 = 0, dU2 = 0, dV2 = 0 ; // se ho trovato due punti vuol dire che la linea coincide con un edge e ho trovato tutto quello che serve if ( ! AreSamePointExact( ptInt2, ORIG)) { if ( i == 0) { //dV1 = 0 ; dV2 = 0 ; vCrvEdge[0]->GetParamAtPoint( ptInt1, dU1) ; vCrvEdge[0]->GetParamAtPoint( ptInt2, dU2) ; } else if ( i == 1) { //dU1 = 1 ; dU2 = 1 ; vCrvEdge[1]->GetParamAtPoint( ptInt1, dV1) ; vCrvEdge[1]->GetParamAtPoint( ptInt2, dV2) ; } else if ( i == 2){ //dV1 = 1 ; dV2 = 1 ; vCrvEdge[2]->GetParamAtPoint( ptInt1, dU1) ; vCrvEdge[2]->GetParamAtPoint( ptInt2, dU2) ; } else if ( i == 3){ //dU1 = 0 ; dU2 = 0 ; vCrvEdge[3]->GetParamAtPoint( ptInt1, dV1) ; vCrvEdge[3]->GetParamAtPoint( ptInt2, dV2) ; } Point3d ptIBez1, ptIBez2 ; Vector3d vtN1, vtN2 ; pSurfBz->GetPointNrmD1D2(dU1, dV1, ISurfBezier::Side::FROM_MINUS, ISurfBezier::Side::FROM_MINUS, ptIBez1, vtN1) ; pSurfBz->GetPointNrmD1D2(dU2, dV2, ISurfBezier::Side::FROM_MINUS, ISurfBezier::Side::FROM_MINUS, ptIBez2, vtN2) ; Point3d ptSP1( dU1, dV1, 0) ; double dCos1 = vtN1 * vtL ; Point3d ptSP2( dU2, dV2, 0) ; double dCos2 = vtN2 * vtL ; // se avevo già trovato un punto singolo che coincide col primo punto di questa intersezione sovrapposta, allora cancello l'intersezione singola che // avevo salvato e aggiungo quella sovrapposto che ho trovato ora if ( bFound) { int nNewTot = int(vInfo.size()) ; int nNewInters = nNewTot - nInters ; bool bAlreadyFound = false ; for ( int i = 0 ; i < nNewInters ; ++i) { bAlreadyFound = AreSamePointApprox(vInfo[nNewTot - i].ptUV, ptSP1) || AreSamePointApprox(vInfo[nNewTot - i].ptUV, ptSP2) ; if ( bAlreadyFound) { vInfo.erase( vInfo.begin() + nNewTot - i) ; break ; } } } UpdateInfoIntersLineSurfBz( ptL, vtL, ILTA_NO_TRIA, -1, ptSP1, ptIBez1, dCos1, ptSP2, ptIBez2, dCos2, vInfo) ; bFound = true ; break ; } // se ho trovato un punto a distanza zero dalla linea allora ho trovato l'intersezione else if ( dMinDist < EPS_SMALL) { if ( i == 0) { //dV1 = 0 ; vCrvEdge[0]->GetParamAtPoint( ptInt1, dU1) ; } else if ( i == 1) { //dU1 = 1 ; vCrvEdge[1]->GetParamAtPoint( ptInt1, dV1) ; } else if ( i == 2) { //dV1 = 1 ; vCrvEdge[2]->GetParamAtPoint( ptInt1, dU1) ; } else if ( i == 3) { //dU1 = 0 ; vCrvEdge[3]->GetParamAtPoint( ptInt1, dV1) ; } Point3d ptSP1( dU1, dV1, 0), ptSP2 ; // se avevo trovato già altri punti controllo di non essere esattamente su una diagonale ( e quindi avere un'intersezione con ogni edge, ma due sono doppie) if ( bFound) { int nNewTot = int(vInfo.size()) ; int nNewInters = nNewTot - nInters ; bool bAlreadyFound = false ; for ( int i = 0 ; i < nNewInters ; ++i) bAlreadyFound = AreSamePointApprox(vInfo[nNewTot - i].ptUV, ptSP1) ; if ( bAlreadyFound) continue ; } Point3d ptIBez1, ptIBez2 ; Vector3d vtN1, vtN2 ; pSurfBz->GetPointNrmD1D2(dU1, dV1, ISurfBezier::Side::FROM_MINUS, ISurfBezier::Side::FROM_MINUS, ptIBez1, vtN1) ; double dCos1 = vtN1 * vtL, dCos2 = 0 ; UpdateInfoIntersLineSurfBz( ptL, vtL, ILTA_NO_TRIA, -1, ptSP1, ptIBez1, dCos1, ptSP2, ptIBez2, dCos2, vInfo) ; bFound = true ; } } } // se la superficie è trimmed verifico che i punti trovati siano all'interno del parametrico trimmato if ( bTrimmed && bFound) { int nNewTot = int(vInfo.size()) ; int nNewInters = nNewTot - nInters ; const ISurfFlatRegion* pFRTrim = pSurfBz->GetTrimRegion() ; for ( int i = 0 ; i < nNewInters ; ++i) { Point3d ptTest = vInfo[nNewTot - i].ptUV * SBZ_TREG_COEFF ; bool bInside = false ; double dDist = INFINITO ; IsPointInsideSurfFr( ptTest, pFRTrim, dDist, bInside) ; if ( ! bInside) vInfo.erase( vInfo.begin() + nNewTot - i) ; } } return true ; }