//---------------------------------------------------------------------------- // EgalTech 2015-2016 //---------------------------------------------------------------------------- // File : VolZmap.cpp Data : 22.01.15 Versione : 1.6a4 // Contenuto : Implementazione della classe Volume Zmap (tre griglie) // // // // Modifiche : 22.01.15 DS Creazione modulo. // // //---------------------------------------------------------------------------- //--------------------------- Include ---------------------------------------- #include "stdafx.h" #include "CurveLine.h" #include "VolZmap.h" #include "GeoConst.h" #include "/EgtDev/Include/EGkIntersLineTria.h" #include "/EgtDev/Include/EGkIntersLinePlane.h" #include "/EgtDev/Include/EGkIntersLineSphere.h" #include "/EgtDev/Include/EGkChainCurves.h" #include "/EgtDev/Include/EgtNumUtils.h" using namespace std ; //---------------------------------------------------------------------------- // Punti e vettore devono esse espressi nel medesimo sistema di riferimento. // Il box è allineato con gli assi di tale sistema di riferimento. // Viene restituito true in caso di intersezione, false altrimenti. bool VolZmap::IntersLineBox( const Point3d& ptP, const Vector3d& vtV, const Point3d& ptMin, const Point3d& ptMax) const { // Il box è allineato agli assi double dU1, dU2 ; return IntersLineBox( ptP, vtV, ptMin, ptMax, dU1, dU2) ; } //---------------------------------------------------------------------------- // Punti e vettore devono esse espressi nel medesimo sistema di riferimento. // Il box è allineato con gli assi di tale sistema di riferimento. // In dU1, dU2 vengono restituiti i parametri a cui si trovano le intersezioni. // Viene restituito true in caso di intersezione, false altrimenti. bool VolZmap::IntersLineBox( const Point3d& ptP, const Vector3d& vtV, const Point3d& ptMin, const Point3d& ptMax, double& dU1, double& dU2) const { // Il box è allineato agli assi dU1 = - INFINITO ; dU2 = INFINITO ; // confronto con piani YZ (perpendicolari ad asse X) if ( vtV.x > EPS_ZERO) { dU1 = max( dU1, ( ptMin.x - ptP.x) / vtV.x) ; dU2 = min( dU2, ( ptMax.x - ptP.x) / vtV.x) ; } else if ( vtV.x < - EPS_ZERO) { dU1 = max( dU1, ( ptMax.x - ptP.x) / vtV.x) ; dU2 = min( dU2, ( ptMin.x - ptP.x) / vtV.x) ; } else if ( ptP.x < ptMin.x - EPS_SMALL || ptP.x > ptMax.x + EPS_SMALL) return false ; // confronto con piani ZX (perpendicolari ad asse Y) if ( vtV.y > EPS_ZERO) { dU1 = max( dU1, ( ptMin.y - ptP.y) / vtV.y) ; dU2 = min( dU2, ( ptMax.y - ptP.y) / vtV.y) ; } else if ( vtV.y < - EPS_ZERO) { dU1 = max( dU1, ( ptMax.y - ptP.y) / vtV.y) ; dU2 = min( dU2, ( ptMin.y - ptP.y) / vtV.y) ; } else if ( ptP.y < ptMin.y - EPS_SMALL || ptP.y > ptMax.y + EPS_SMALL) return false ; // confronto con piani XY (perpendicolari ad asse Z) if ( vtV.z > EPS_ZERO) { dU1 = max( dU1, ( ptMin.z - ptP.z) / vtV.z) ; dU2 = min( dU2, ( ptMax.z - ptP.z) / vtV.z) ; } else if ( vtV.z < - EPS_ZERO) { dU1 = max( dU1, ( ptMax.z - ptP.z) / vtV.z) ; dU2 = min( dU2, ( ptMin.z - ptP.z) / vtV.z) ; } else if ( ptP.z < ptMin.z - EPS_SMALL || ptP.z > ptMax.z + EPS_SMALL) return false ; return ( dU2 >= dU1) ; } //---------------------------------------------------------------------------- // La retta è rappresentata da un punto e dal versore espressi nel riferimento dello Zmap. // Se non vi è intersezione, viene restituito falso bool VolZmap::IntersLineZMapLattice( const Point3d& ptP, const Vector3d& vtV, double& dU1, double& dU2) const { // Punti estremi della diagonale del box del reticolo Point3d ptMin( - 0.5 * m_dStep, - 0.5 * m_dStep, - 0.5 * m_dStep) ; Point3d ptMax( ( m_nNx[0] + 0.5) * m_dStep, ( m_nNy[0] + 0.5) * m_dStep, ( m_nNy[1] + 0.5) * m_dStep) ; return IntersLineBox( ptP, vtV, ptMin, ptMax, dU1, dU2) ; } //---------------------------------------------------------------------------- // Trova le intersezioni di una retta, definita tramite un suo punto e il vettore della direzione, // e il bounding box dello Zmap. Punto e vettore devono essere espressi nel sistema di riferimento // intrinseco dello Zmap. Se vi è intersezione viene restitutito true, false altrimenti. // Se vi è intersezione in dU1 e dU2 vengono restituiti i parametri a cui avvengono le intersezioni. bool VolZmap::IntersLineZMapBBox( const Point3d& ptP, const Vector3d& vtV, double& dU1, double& dU2) const { // Punti estremi del box dello Zmap Point3d ptMin = ORIG ; Point3d ptMax = ptMin + Vector3d( m_nNx[0] * m_dStep, m_nNy[0] * m_dStep, m_dMaxZ[0]) ; return IntersLineBox( ptP, vtV, ptMin, ptMax, dU1, dU2) ; } //---------------------------------------------------------------------------- // Trova le intersezioni fra una retta e un dexel parallelepipedale (gli intervalli sono // parallelepipdi). Punto e vettore devono essere espressi nel sistema grilgia che coincide // con quello intrinseco dello Zmap solo nel caso della prima griglia. Per le altre griglie // è necessario permutare ciclicamente le coordinate. bool VolZmap::IntersLineDexel( const Point3d& ptP, const Vector3d& vtV, int nGrid, int nI, int nJ, double& dU1, double& dU2) const { // Determino l'indice del dexel e il numero di suoi intervalli int nDexelPos = nJ * m_nNx[nGrid] + nI ; int nDexelSize = int( m_Values[nGrid][nDexelPos].size()) ; // Se non c'è materiale non devo fare alcunché if ( nDexelSize == 0) return false ; // Determino estremi nel piano XY intrinseco del dexel double dXmin = nI * m_dStep ; double dYmin = nJ * m_dStep ; double dXmax = ( nI + 1) * m_dStep ; double dYmax = ( nJ + 1) * m_dStep ; // ciclo sugli intervalli dU1 = INFINITO ; dU2 = - INFINITO ; bool bInters = false ; for ( int nIndex = 0 ; nIndex < nDexelSize ; nIndex += 1) { // estremi del box del singolo intervallo Point3d ptE1( dXmin, dYmin, m_Values[nGrid][nDexelPos][nIndex].dMin) ; Point3d ptE2( dXmax, dYmax, m_Values[nGrid][nDexelPos][nIndex].dMax) ; double dP1, dP2 ; if ( IntersLineBox( ptP, vtV, ptE1, ptE2, dP1, dP2)) { bInters = true ; dU1 = min( dU1, dP1) ; dU2 = max( dU2, dP2) ; } } return bInters ; } //---------------------------------------------------------------------------- // Trova le intersezioni fra una semi-retta e un dexel parallelepipedale (gli intervalli sono // parallelepipdi). Punto e vettore devono essere espressi nel sistema grilgia che coincide // con quello intrinseco dello Zmap solo nel caso della prima griglia. Per le altre griglie // è necessario permutare ciclicamente le coordinate. bool VolZmap::IntersRayDexel( const Point3d& ptP, const Vector3d& vtV, int nGrid, int nI, int nJ, double& dU1, double& dU2) const { // Determino l'indice del dexel e il numero di suoi intervalli int nDexelPos = nJ * m_nNx[nGrid] + nI ; int nDexelSize = int( m_Values[nGrid][nDexelPos].size()) ; // Se non c'è materiale non devo fare alcunché if ( nDexelSize == 0) return false ; // Determino estremi nel piano XY intrinseco del dexel double dXmin = nI * m_dStep ; double dYmin = nJ * m_dStep ; double dXmax = ( nI + 1) * m_dStep ; double dYmax = ( nJ + 1) * m_dStep ; // ciclo sugli intervalli dU1 = INFINITO ; dU2 = - INFINITO ; bool bInters = false ; for ( int nIndex = 0 ; nIndex < nDexelSize ; nIndex += 1) { // estremi del box del singolo intervallo Point3d ptE1( dXmin, dYmin, m_Values[nGrid][nDexelPos][nIndex].dMin) ; Point3d ptE2( dXmax, dYmax, m_Values[nGrid][nDexelPos][nIndex].dMax) ; double dP1, dP2 ; if ( IntersLineBox( ptP, vtV, ptE1, ptE2, dP1, dP2)) { if ( dP2 >= 0) { bInters = true ; if ( dP1 >= 0) dU1 = min( dU1, dP1) ; else dU1 = - 1 ; dU2 = max( dU2, dP2) ; } } } return bInters ; } //---------------------------------------------------------------------------- // Calcola la profondità del materiale lungo un raggio, identificato da un punto e un versoreore. // Punto e versore devono essere espressi nel sistema di riferimento locale, in cui è immerso quello dello Zmap. // InLength = distanza di ingresso (se -1 il punto è interno, se -2 il punto è esterno e il raggio non interseca lo Zmap) // OutLength = distanza di uscita bool VolZmap::GetDepth( const Point3d& ptPLoc, const Vector3d& vtDLoc, double& dInLength, double& dOutLength, bool bExact) const { Point3d ptP = ptPLoc ; Vector3d vtD = vtDLoc ; ptP.ToLoc( m_MapFrame) ; vtD.ToLoc( m_MapFrame) ; if ( bExact && m_nMapNum == 3) return GetDepthWithVoxel( ptP, vtD, dInLength, dOutLength, true) ; else return GetDepthWithDexel( ptP, vtD, dInLength, dOutLength) ; } //---------------------------------------------------------------------------- // Calcola la profondità del materiale usando i parallelepipedi corrispondenti ai segmenti di dexel di una griglia. // Punto e versore devono essere espressi nel sistema di riferimento Zmap. Il vettore deve essere normalizzato. // InLength = distanza di ingresso (se -1 il punto è interno, se -2 il punto è esterno e il raggio non interseca lo Zmap) // OutLength = distanza di uscita bool VolZmap::GetDepthWithDexel( const Point3d& ptP, const Vector3d& vtV, double& dInLength, double& dOutLength) const { // Intersezione fra semiretta e BBox dello Zmap double dU1, dU2 ; bool bTest = IntersLineZMapBBox( ptP, vtV, dU1, dU2) && ( dU1 > 0 || dU2 > 0) ; // Semiretta esterna al box dello Zmap quindi esterna anche allo Zmap if ( ! bTest) { dInLength = - 2 ; dOutLength = - 2 ; return true ; } double dInLen[3] ; double dOutLen[3] ; // Ciclo sulle griglie for ( int nGrid = 0 ; nGrid < int( m_nMapNum) ; ++ nGrid) { Point3d ptP0 = ptP ; Vector3d vtV0 = vtV ; Point3d ptI, ptF ; // Una sola intersezione valida ( punto interno, intersezione valida 2) if ( dU1 < 0 && dU2 > 0) { ptI = ptP ; ptF = ptP + dU2 * vtV ; } // due soluzioni valide ( punto esterno) else { ptI = ptP + dU1 * vtV ; ptF = ptP + dU2 * vtV ; } // Passo dal sistema intrinseco alla griglia if ( nGrid == 1) { swap( ptP0.y, ptP0.z) ; swap( ptP0.x, ptP0.z) ; swap( vtV0.y, vtV0.z) ; swap( vtV0.x, vtV0.z) ; swap( ptI.y, ptI.z) ; swap( ptI.x, ptI.z) ; swap( ptF.y, ptF.z) ; swap( ptF.x, ptF.z) ; } else if ( nGrid == 2) { swap( ptP0.x, ptP0.z) ; swap( ptP0.y, ptP0.z) ; swap( vtV0.x, vtV0.z) ; swap( vtV0.y, vtV0.z) ; swap( ptI.x, ptI.z) ; swap( ptI.y, ptI.z) ; swap( ptF.x, ptF.z) ; swap( ptF.y, ptF.z) ; } // Inizializzo distanze dInLen[nGrid] = INFINITO ; dOutLen[nGrid] = - INFINITO ; // Determino asse di spostamento maggiore bool bOnX = ( abs( ptF.y - ptI.y) <= abs( ptF.x - ptI.x)) ; // Movimento crescente su asse di movimento principale if ( ( bOnX && ptF.x < ptI.x) || ( ! bOnX && ptF.y < ptI.y) ) swap( ptI, ptF) ; // Pendenza double dDeltaX = ( ptF.x - ptI.x) ; double dDeltaY = ( ptF.y - ptI.y) ; double dM = 0 ; if ( bOnX && abs( dDeltaX) > EPS_SMALL) dM = dDeltaY / dDeltaX ; else if ( ! bOnX && abs( dDeltaY) > EPS_SMALL) dM = dDeltaX / dDeltaY ; // Determinazione degli indici i j dei punti ptI e ptF int nIi = Clamp( int( floor( ptI.x / m_dStep)), 0, m_nNx[nGrid] - 1) ; int nIj = Clamp( int( floor( ptI.y / m_dStep)), 0, m_nNy[nGrid] - 1) ; int nFi = Clamp( int( floor( ptF.x / m_dStep)), 0, m_nNx[nGrid] - 1) ; int nFj = Clamp( int( floor( ptF.y / m_dStep)), 0, m_nNy[nGrid] - 1) ; // mi muovo int i = nIi ; int j = nIj ; while ( ( bOnX && i <= nFi) || ( ! bOnX && j <= nFj)) { // Eseguo controllo double dU1, dU2 ; if ( IntersRayDexel( ptP0, vtV0, nGrid, i, j, dU1, dU2)) { dInLen[nGrid] = min( dInLen[nGrid], dU1) ; dOutLen[nGrid] = max( dOutLen[nGrid], dU2) ; } // Calcolo spostamento (a destra o sopra) if ( bOnX) { double dMoveX = ( ( i + 1) * m_dStep - ptI.x) ; double dMoveY = dMoveX * dM ; double dY = ptI.y + dMoveY ; int nNewJ = Clamp( int( floor( dY / m_dStep)), 0, m_nNy[nGrid] - 1) ; if ( nNewJ != j) j = nNewJ ; else ++ i ; } else { double dMoveY = ( ( j + 1) * m_dStep - ptI.y) ; double dMoveX = dMoveY * dM ; double dX = ptI.x + dMoveX ; int nNewI = Clamp( int( floor( dX / m_dStep)), 0, m_nNx[nGrid] - 1) ; if ( nNewI != i) i = nNewI ; else ++ j ; } } // Se non abbiamo incontrato materiale if ( dInLen[nGrid] > dOutLen[nGrid] - EPS_SMALL) { dInLen[nGrid] = - 2 ; dOutLen[nGrid] = - 2 ; } // Se parto dall'interno else if ( dInLen[nGrid] < - EPS_SMALL) dInLen[nGrid] = - 1 ; } if ( m_nMapNum == 1) { dInLength = dInLen[0] ; dOutLength = dOutLen[0] ; return true ; } if ( abs( dInLen[0] + 2) < EPS_SMALL && abs( dInLen[1] + 2) < EPS_SMALL && abs( dInLen[2] + 2) < EPS_SMALL) { dInLength = - 2 ; dOutLength = - 2 ; return true ; } else { dInLength = INFINITO ; dOutLength = - INFINITO ; for ( int nGr = 0 ; nGr < 3 ; ++ nGr) { if ( dInLen[nGr] > - 2 && dInLen[nGr] < dInLength) dInLength = dInLen[nGr] ; if ( dOutLen[nGr] > - 2 && dOutLen[nGr] > dOutLength) dOutLength = dOutLen[nGr] ; } } return true ; } //---------------------------------------------------------------------------- // Calcola la profondità del materiale usando i triangoli della rappresentazione grafica. // Punto e versore devono essere espressi nel sistema di riferimento Zmap. // Se si vuole che vengano utilizzati i triangoli con sharp-featue bEnh deve valere true, false altrimenti. // InLength = distanza di ingresso (se -1 il punto è interno, se -2 il punto è esterno e il raggio non interseca lo Zmap) // OutLength = distanza di uscita. bool VolZmap::GetDepthWithVoxel( const Point3d& ptP, const Vector3d& vtD, double& dInLength, double& dOutLength, bool bEnh) const { // Serve che punto e vettore siano espressi sia nel sistema intrinseco dello Zmap (m_MapFrame) sia in quello // in cui esso è immerso; questo perché i dexel sono espressi in quello intrinseco e i triangoli in quello // in cui esso è immerso. Point3d ptOutP = ptP ; Vector3d vtOutD = vtD ; ptOutP.ToGlob( m_MapFrame) ; vtOutD.ToGlob( m_MapFrame) ; // Intersezione fra semiretta e BBox dello Zmap double dU1, dU2 ; bool bLineBBoxInters = IntersLineZMapBBox( ptP, vtD, dU1, dU2) && ( dU1 > 0 || dU2 > 0) ; // Semiretta esterna al box dello Zmap quindi esterna anche allo Zmap if ( ! bLineBBoxInters) { dInLength = - 2 ; dOutLength = - 2 ; return true ; } // Se la grafica non è aggiornata la ricalcolo bool bUpGrade = false ; for ( int nBl = 0 ; nBl < int( m_nNumBlock) ; ++ nBl) { bUpGrade = bUpGrade || m_BlockToUpdate[nBl] ; } INTVECTOR nModifiedBlocks ; if ( bUpGrade) UpGradeGraphics( false, nModifiedBlocks) ; // Determino il voxel di partenza int nVoxI, nVoxJ, nVoxK ; if ( ! GetPointVoxel( ptP, nVoxI, nVoxJ, nVoxK)) { if ( ! GetPointVoxel( ptP + dU1 * vtD, nVoxI, nVoxJ, nVoxK)) return false ; } // Incrementi degli indici per sguire la retta int nDeltaI = ( vtD.x > 0 ? 1 : ( vtD.x < 0 ? - 1 : 0)) ; int nDeltaJ = ( vtD.y > 0 ? 1 : ( vtD.y < 0 ? - 1 : 0)) ; int nDeltaK = ( vtD.z > 0 ? 1 : ( vtD.z < 0 ? - 1 : 0)) ; // Scelgo i piani del voxel con cui intersecare la retta, per determinare il voxel successivo int nPlaneI = ( vtD.x >= 0 ? 1 : 0) ; int nPlaneJ = ( vtD.y >= 0 ? 1 : 0) ; int nPlaneK = ( vtD.z >= 0 ? 1 : 0) ; // Ciclo sui voxel vector vVox ; while ( IsValidVoxel( nVoxI, nVoxJ, nVoxK)) { // Estremi della diagonale del voxel corrente Point3d ptMin( ( nVoxI * N_DEXVOXRATIO + 0.5) * m_dStep, ( nVoxJ * N_DEXVOXRATIO + 0.5) * m_dStep, ( nVoxK * N_DEXVOXRATIO + 0.5) * m_dStep) ; Point3d ptMax( ( ( nVoxI + 1) * N_DEXVOXRATIO + 0.5) * m_dStep, ( ( nVoxJ + 1) * N_DEXVOXRATIO + 0.5) * m_dStep, ( ( nVoxK + 1) * N_DEXVOXRATIO + 0.5) * m_dStep) ; // Studio il voxel corrente if ( IntersLineBox( ptP, vtD, ptMin, ptMax)) { int nCurVoxIndex = CalcIndex( nVoxI, nVoxJ, nVoxK) ; if ( nCurVoxIndex != 0 && nCurVoxIndex != 255) { VoxelIndexes NewVox ; NewVox.nI = nVoxI ; NewVox.nJ = nVoxJ ; NewVox.nK = nVoxK ; vVox.emplace_back( NewVox) ; } } // Interseco la retta con i piani frontiera del voxel double dMaxTX = ( abs( vtD.x) > EPS_ZERO ? abs( ( ( ( nVoxI + nPlaneI) * N_DEXVOXRATIO + 0.5) * m_dStep - ptP.x) / vtD.x) : INFINITO) ; double dMaxTY = ( abs( vtD.y) > EPS_ZERO ? abs( ( ( ( nVoxJ + nPlaneJ) * N_DEXVOXRATIO + 0.5) * m_dStep - ptP.y) / vtD.y) : INFINITO) ; double dMaxTZ = ( abs( vtD.z) > EPS_ZERO ? abs( ( ( ( nVoxK + nPlaneK) * N_DEXVOXRATIO + 0.5) * m_dStep - ptP.z) / vtD.z) : INFINITO) ; // Determino gli incrementi if ( dMaxTX < dMaxTY) { if ( dMaxTX < dMaxTZ) nVoxI += nDeltaI ; else nVoxK += nDeltaK ; } else { if ( dMaxTY < dMaxTZ) nVoxJ += nDeltaJ ; else nVoxK += nDeltaK ; } } // Dati dell'intersezione struct LineTriaInt { int nNum ; double dPar1 ; double dPar2 ; double dDot ; LineTriaInt( void) : nNum( 0) {} LineTriaInt( double dP, double dD) : nNum( 1), dPar1( dP), dDot( dD) {} LineTriaInt( double dP1, double dP2, double dD) : nNum( 2), dPar1( dP1), dPar2( dP2), dDot( dD) {} } ; vector vInt ; int nPrevBlockN = - 1 ; // Ciclo sui voxel for ( int nVx = 0 ; nVx < int( vVox.size()) ; ++ nVx) { int nCurVoxIJK[3] = { vVox[nVx].nI, vVox[nVx].nJ, vVox[nVx].nK} ; int nCurBlockIJK[3] ; if ( GetVoxelBlockIJK( nCurVoxIJK, nCurBlockIJK)) { int nCurBlockN ; GetBlockNFromIJK( nCurBlockIJK, nCurBlockN) ; // Triangoli sharp fra blocchi for ( int nBlVx = 0 ; nBlVx < int( m_InterBlockToBeFlippedSharpTria[nCurBlockN].size()) ; ++ nBlVx) { if ( m_InterBlockToBeFlippedSharpTria[nCurBlockN][nBlVx].i == nCurVoxIJK[0] && m_InterBlockToBeFlippedSharpTria[nCurBlockN][nBlVx].j == nCurVoxIJK[1] && m_InterBlockToBeFlippedSharpTria[nCurBlockN][nBlVx].k == nCurVoxIJK[2]) { for ( int nBlCm = 0 ; nBlCm < int( m_InterBlockToBeFlippedSharpTria[nCurBlockN][nBlVx].vCompoTria.size()) ; ++ nBlCm) { for ( int nBlTr = 0 ; nBlTr < int( m_InterBlockToBeFlippedSharpTria[nCurBlockN][nBlVx].vCompoTria[nBlCm].size()) ; ++ nBlTr) { Triangle3d CurrTria = m_InterBlockToBeFlippedSharpTria[nCurBlockN][nBlVx].vCompoTria[nBlCm][nBlTr] ; Point3d ptLineTria1, ptLineTria2 ; // Studio dell'intersezione della retta con il triangolo corrente int nIntType = IntersLineTria( ptOutP, vtOutD, 1.5 * dU2, CurrTria, ptLineTria1, ptLineTria2) ; // Se non ci sono intersezioni passo al prossimo triangolo if ( nIntType == ILTT_NO) continue ; // se altrimenti c'è una sola intersezione else if ( nIntType == ILTT_VERT || nIntType == ILTT_EDGE || nIntType == ILTT_IN) { vInt.emplace_back( ( ptLineTria1 - ptOutP) * vtOutD, vtOutD * CurrTria.GetN()) ; } // altrimenti ci sono due intersezioni else { double dP1 = ( ptLineTria1 - ptOutP) * vtOutD ; double dP2 = ( ptLineTria2 - ptOutP) * vtOutD ; double dD = vtOutD * CurrTria.GetN() ; vInt.emplace_back( ( dP1 < dP2 ? dP1 : dP2), ( dP1 < dP2 ? dP2 : dP1), dD) ; } } } } } // Triangoli sharp interni for ( int nBlVx = 0 ; nBlVx < int( m_BlockSharpTria[nCurBlockN].size()) ; ++ nBlVx) { if ( m_BlockSharpTria[nCurBlockN][nBlVx].i == nCurVoxIJK[0] && m_BlockSharpTria[nCurBlockN][nBlVx].j == nCurVoxIJK[1] && m_BlockSharpTria[nCurBlockN][nBlVx].k == nCurVoxIJK[2]) { for ( int nBlCm = 0 ; nBlCm < int( m_BlockSharpTria[nCurBlockN][nBlVx].vCompoTria.size()) ; ++ nBlCm) { for ( int nBlTr = 0 ; nBlTr < int( m_BlockSharpTria[nCurBlockN][nBlVx].vCompoTria[nBlCm].size()) ; ++ nBlTr) { Triangle3d CurrTria = m_BlockSharpTria[nCurBlockN][nBlVx].vCompoTria[nBlCm][nBlTr] ; Point3d ptLineTria1, ptLineTria2 ; // Studio dell'intersezione della retta con il triangolo corrente int nIntType = IntersLineTria( ptOutP, vtOutD, 1.5 * dU2, CurrTria, ptLineTria1, ptLineTria2) ; // Se non ci sono intersezioni passo al prossimo triangolo if ( nIntType == ILTT_NO) continue ; // se altrimenti c'è una sola intersezione else if ( nIntType == ILTT_VERT || nIntType == ILTT_EDGE || nIntType == ILTT_IN) { vInt.emplace_back( ( ptLineTria1 - ptOutP) * vtOutD, vtOutD * CurrTria.GetN()) ; } // altrimenti ci sono due intersezioni else { double dP1 = ( ptLineTria1 - ptOutP) * vtOutD ; double dP2 = ( ptLineTria2 - ptOutP) * vtOutD ; double dD = vtOutD * CurrTria.GetN() ; vInt.emplace_back( ( dP1 < dP2 ? dP1 : dP2), ( dP1 < dP2 ? dP2 : dP1), dD) ; } } } } } // Triangoli smooth for ( int nBlVx = 0 ; nBlVx < int( m_BlockSmoothTria[nCurBlockN].size()) ; ++ nBlVx) { if ( m_BlockSmoothTria[nCurBlockN][nBlVx].i == nCurVoxIJK[0] && m_BlockSmoothTria[nCurBlockN][nBlVx].j == nCurVoxIJK[1] && m_BlockSmoothTria[nCurBlockN][nBlVx].k == nCurVoxIJK[2]) { for ( int nBlTr = 0 ; nBlTr < int( m_BlockSmoothTria[nCurBlockN][nBlVx].vTria.size()) ; ++ nBlTr) { Triangle3d CurrTria = m_BlockSmoothTria[nCurBlockN][nBlVx].vTria[nBlTr] ; Point3d ptLineTria1, ptLineTria2 ; // Studio dell'intersezione della retta con il triangolo corrente int nIntType = IntersLineTria( ptOutP, vtOutD, 1.5 * dU2, CurrTria, ptLineTria1, ptLineTria2) ; // Se non ci sono intersezioni passo al prossimo triangolo if ( nIntType == ILTT_NO) continue ; // se altrimenti c'è una sola intersezione else if ( nIntType == ILTT_VERT || nIntType == ILTT_EDGE || nIntType == ILTT_IN) { vInt.emplace_back( ( ptLineTria1 - ptOutP) * vtOutD, vtOutD * CurrTria.GetN()) ; } // altrimenti ci sono due intersezioni else { double dP1 = ( ptLineTria1 - ptOutP) * vtOutD ; double dP2 = ( ptLineTria2 - ptOutP) * vtOutD ; double dD = vtOutD * CurrTria.GetN() ; vInt.emplace_back( ( dP1 < dP2 ? dP1 : dP2), ( dP1 < dP2 ? dP2 : dP1), dD) ; } } } } // Triangoli grandi if ( nCurBlockN != nPrevBlockN) { for ( int nBlTr = 0 ; nBlTr < int( m_BlockBigTria[nCurBlockN].size()) ; ++ nBlTr) { Triangle3d CurrTria = m_BlockBigTria[nCurBlockN][nBlTr] ; Point3d ptLineTria1, ptLineTria2 ; // Studio dell'intersezione della retta con il triangolo corrente int nIntType = IntersLineTria( ptOutP, vtOutD, 1.5 * dU2, CurrTria, ptLineTria1, ptLineTria2) ; // Se non ci sono intersezioni passo al prossimo triangolo if ( nIntType == ILTT_NO) continue ; // se altrimenti c'è una sola intersezione else if ( nIntType == ILTT_VERT || nIntType == ILTT_EDGE || nIntType == ILTT_IN) { vInt.emplace_back( ( ptLineTria1 - ptOutP) * vtOutD, vtOutD * CurrTria.GetN()) ; } // altrimenti ci sono due intersezioni else { double dP1 = ( ptLineTria1 - ptOutP) * vtOutD ; double dP2 = ( ptLineTria2 - ptOutP) * vtOutD ; double dD = vtOutD * CurrTria.GetN() ; vInt.emplace_back( ( dP1 < dP2 ? dP1 : dP2), ( dP1 < dP2 ? dP2 : dP1), dD) ; } } nPrevBlockN = nCurBlockN ; } } } // Ciclo sui triangoli dei voxel //for ( auto it = lstTria.begin() ; it != lstTria.end() ; ++it ) { // // Triangolo corrente e suoi punti di intersezione con la retta // const Triangle3d& CurrTria = *it ; // Point3d ptLineTria1, ptLineTria2 ; // // Studio dell'intersezione della retta con il triangolo corrente // int nIntType = IntersLineTria( ptP, vtD, 1.5 * dU2, CurrTria, ptLineTria1, ptLineTria2) ; // // Se non ci sono intersezioni passo al prossimo triangolo // if ( nIntType == ILTT_NO) // continue ; // // se altrimenti c'è una sola intersezione // else if ( nIntType == ILTT_VERT || // nIntType == ILTT_EDGE || // nIntType == ILTT_IN) { // vInt.emplace_back( ( ptLineTria1 - ptP) * vtD, vtD * CurrTria.GetN()) ; // } // // altrimenti ci sono due intersezioni // else { // double dP1 = ( ptLineTria1 - ptP) * vtD ; // double dP2 = ( ptLineTria2 - ptP) * vtD ; // double dD = vtD * CurrTria.GetN() ; // vInt.emplace_back( ( dP1 < dP2 ? dP1 : dP2), ( dP1 < dP2 ? dP2 : dP1), dD) ; // } //} // Ordino le intersezioni in base al parametro distanza con segno da ptP sort( vInt.begin(), vInt.end(), []( const LineTriaInt& a, const LineTriaInt& b) { double dFP = ( a.nNum == 2 ? 0.5 * ( a.dPar1 + a.dPar2) : a.dPar1) ; double dLP = ( b.nNum == 2 ? 0.5 * ( b.dPar1 + b.dPar2) : b.dPar1) ; return ( dLP > dFP) ; }) ; // Inizializzo le distanze di ingresso e uscita: // dInLength diminuisce, dOutLength aumenta. dInLength = INFINITO ; dOutLength = - INFINITO ; int nFirstPosN ; int nN = 0 ; for ( ; nN < int( vInt.size()) ; ++ nN) { if ( vInt[nN].dPar1 > - EPS_SMALL) { nFirstPosN = nN ; break ; } } if ( nN == int( vInt.size())) { dInLength = - 2 ; dOutLength = - 2 ; return true ; } if ( nFirstPosN > 0) { if ( vInt[nFirstPosN - 1].dDot < EPS_ZERO) dInLength = -1 ; } else if ( nFirstPosN == 0) { if ( vInt[nFirstPosN].dDot > EPS_ZERO) dInLength = -1 ; else if ( GetPointVoxel( ptP, nVoxI, nVoxJ, nVoxK)) { int nCubeType = CalcIndex( nVoxI, nVoxJ, nVoxK) ; if ( nCubeType == 255) dInLength = -1 ; } } for ( int nN = nFirstPosN ; nN < int( vInt.size()) ; ++ nN) { if ( vInt[nN].dDot > - EPS_ZERO) { if ( vInt[nN].nNum == 2 && dOutLength < vInt[nN].dPar2) dOutLength = vInt[nN].dPar2 ; else if ( dOutLength < vInt[nN].dPar1) dOutLength = vInt[nN].dPar1 ; } if ( vInt[nN].dDot < EPS_ZERO && dInLength > vInt[nN].dPar1) dInLength = vInt[nN].dPar1 ; } return true ; } //---------------------------------------------------------------------------- bool VolZmap::AvoidBox( const Frame3d& frBox, const Vector3d& vtDiag) const { // BBox BBox3d b3Box( ORIG, ORIG + vtDiag) ; // lo porto nel riferimento intrinseco dello Zmap b3Box.LocToLoc( frBox, m_MapFrame) ; // BBox dello Zmap nel suo riferimento intrinseco BBox3d b3Zmap( ORIG, Point3d( m_nNx[0] * m_dStep, m_nNy[0] * m_dStep, m_dMaxZ[0])) ; // Se non interferiscono, posso uscire BBox3d b3Int ; if ( ! b3Zmap.FindIntersection( b3Box, b3Int)) return true ; // Limiti su indici int nStI = Clamp( int( b3Int.GetMin().x / m_dStep), 0, m_nNx[0] -1) ; int nEnI = Clamp( int( b3Int.GetMax().x / m_dStep), 0, m_nNx[0] -1) ; int nStJ = Clamp( int( b3Int.GetMin().y / m_dStep), 0, m_nNy[0] -1) ; int nEnJ = Clamp( int( b3Int.GetMax().y / m_dStep), 0, m_nNy[0] -1) ; // Vettore direzione dei dexel nel riferimento del Box Vector3d vtK = Z_AX ; vtK.LocToLoc( m_MapFrame, frBox) ; // Riferimento intrinseco dei dexel nel riferimento del box Point3d ptO = ORIG ; ptO.LocToLoc( m_MapFrame, frBox) ; Vector3d vtX = X_AX ; vtX.LocToLoc( m_MapFrame, frBox) ; Vector3d vtY = Y_AX ; vtY.LocToLoc( m_MapFrame, frBox) ; // Ciclo di intersezione dei dexel con il BBox for ( int i = nStI ; i <= nEnI ; ++ i) { for ( int j = nStJ ; j <= nEnJ ; ++ j) { int nPos = j * m_nNx[0] + i ; int nSize = int( m_Values[0][nPos].size()) ; if ( nSize == 0) continue ; for ( int k = 0 ; k < 5 ; ++ k) { Point3d ptT = ptO + ( i + 0.5) * m_dStep * vtX + ( j + 0.5) * m_dStep * vtY ; if ( k == 0) ; else if ( k == 1) ptT += - 0.4 * m_dStep * vtX - 0.4 * m_dStep * vtY ; else if ( k == 2) ptT += + 0.4 * m_dStep * vtX - 0.4 * m_dStep * vtY ; else if ( k == 3) ptT += + 0.4 * m_dStep * vtX + 0.4 * m_dStep * vtY ; else if ( k == 4) ptT += - 0.4 * m_dStep * vtX + 0.4 * m_dStep * vtY ; double dZmin, dZmax ; if ( IntersLineBox( ptT, vtK, ORIG, ORIG + vtDiag, dZmin, dZmax)) { for ( int nIndex = 0 ; nIndex < nSize ; nIndex += 1) { if ( ! ( dZmax < m_Values[0][nPos][nIndex].dMin - EPS_SMALL || dZmin > m_Values[0][nPos][nIndex].dMax + EPS_SMALL)) return false ; } } } } } return true ; } //---------------------------------------------------------------------------- bool VolZmap::AvoidSphere( const Point3d& ptCenter, double dRad) const { // Porto la sfera nel riferimento intrinseco dello Zmap Point3d ptC = ptCenter ; ptC.ToLoc( m_MapFrame) ; // BBox della sfera BBox3d b3Box( ptC) ; b3Box.Expand( dRad) ; // BBox dello Zmap nel suo riferimento intrinseco BBox3d b3Zmap( ORIG, Point3d( m_nNx[0] * m_dStep, m_nNy[0] * m_dStep, m_dMaxZ[0])) ; // Se non interferiscono, posso uscire BBox3d b3Int ; if ( ! b3Zmap.FindIntersection( b3Box, b3Int)) return true ; // Limiti su indici int nStI = Clamp( int( b3Int.GetMin().x / m_dStep), 0, m_nNx[0] -1) ; int nEnI = Clamp( int( b3Int.GetMax().x / m_dStep), 0, m_nNx[0] -1) ; int nStJ = Clamp( int( b3Int.GetMin().y / m_dStep), 0, m_nNy[0] -1) ; int nEnJ = Clamp( int( b3Int.GetMax().y / m_dStep), 0, m_nNy[0] -1) ; // Ciclo di intersezione dei dexel con la sfera (nel riferimento intrinseco) for ( int i = nStI ; i <= nEnI ; ++ i) { for ( int j = nStJ ; j <= nEnJ ; ++ j) { int nPos = j * m_nNx[0] + i ; int nSize = int( m_Values[0][nPos].size()) ; if ( nSize == 0) continue ; for ( int k = 0 ; k < 5 ; ++ k) { Point3d ptT = ORIG + ( i + 0.5) * m_dStep * X_AX + ( j + 0.5) * m_dStep * Y_AX ; if ( k == 0) ; else if ( k == 1) ptT += - 0.4 * m_dStep * X_AX - 0.4 * m_dStep * Y_AX ; else if ( k == 2) ptT += + 0.4 * m_dStep * X_AX - 0.4 * m_dStep * Y_AX ; else if ( k == 3) ptT += + 0.4 * m_dStep * X_AX + 0.4 * m_dStep * Y_AX ; else if ( k == 4) ptT += - 0.4 * m_dStep * X_AX + 0.4 * m_dStep * Y_AX ; Point3d ptI1, ptI2 ; if ( IntersLineSphere( ptT, Z_AX, ptC, dRad, ptI1, ptI2) != ILST_NO) { double dZmin = ptI1.z ; double dZmax = ptI2.z ; for ( int nIndex = 0 ; nIndex < nSize ; nIndex += 1) { if ( ! ( dZmax < m_Values[0][nPos][nIndex].dMin - EPS_SMALL || dZmin > m_Values[0][nPos][nIndex].dMax + EPS_SMALL)) return false ; } } } } } return true ; } //---------------------------------------------------------------------------- bool VolZmap::AvoidCylinder( const Frame3d& frCyl, double dL, double dR) const { return true ; } //---------------------------------------------------------------------------- // Riferimento con origine nel centro della base e asse di simmetria coincidente con l'asse Z. // La funzione restituisce true in caso di intersezione, false altrimenti. //---------------------------------------------------------------------------- bool VolZmap::IntersLineCylinder( const Point3d& ptLineSt, const Vector3d& vtLineDir, const Frame3d& CylFrame, double dH, double dRad, bool bTapLow, bool bTapUp, Point3d& ptInt1, Vector3d& vtN1, Point3d& ptInt2, Vector3d& vtN2) { // Porto la linea nel riferimento del cilindro Point3d ptP = ptLineSt ; ptP.ToLoc( CylFrame) ; Vector3d vtV = vtLineDir ; vtV.ToLoc( CylFrame) ; // Determino le eventuali intersezioni con le due basi a quota minima e massima (solo se linea non parallela ad esse) int nBasInt = 0 ; if ( abs( vtV.z) > EPS_ZERO) { // le linee tangenti al cilindro non sono considerate intersecanti double EpsRad = ( vtV.IsZeroXY() ? - EPS_SMALL : EPS_SMALL) ; ptInt1 = ptP + ( ( 0 - ptP.z) / vtV.z) * vtV ; if ( ptInt1.x * ptInt1.x + ptInt1.y * ptInt1.y < dRad * dRad + 2 * dRad * EpsRad) { nBasInt += 1 ; vtN1 = Z_AX ; } ptInt2 = ptP + ( ( dH - ptP.z) / vtV.z) * vtV ; if ( ptInt2.x * ptInt2.x + ptInt2.y * ptInt2.y < dRad * dRad + 2 * dRad * EpsRad) { nBasInt += 2 ; vtN2 = - Z_AX ; } } // Se la linea interseca entrambe le basi, si sono trovate le due intersezioni if ( nBasInt == 3) { // Porto i punti e i versori nel riferimento globale ptInt1.ToGlob( CylFrame) ; vtN1.ToGlob( CylFrame) ; ptInt2.ToGlob( CylFrame) ; vtN2.ToGlob( CylFrame) ; // Trovate intersezioni return true ; } // Determino le intersezioni con la superficie laterale del cilindro DBLVECTOR vdCoef(3) ; double dSqRad = dRad * dRad ; vdCoef[0] = ptP.x * ptP.x + ptP.y * ptP.y - dSqRad ; vdCoef[1] = 2 * ( ptP.x * vtV.x + ptP.y * vtV.y) ; vdCoef[2] = vtV.x * vtV.x + vtV.y * vtV.y ; DBLVECTOR vdRoots ; int nRoot = PolynomialRoots( 2, vdCoef, vdRoots) ; // Epsilon per piani di tappo double dEpsLow = ( bTapLow ? - EPS_SMALL : EPS_SMALL) ; double dEpsUp = ( bTapUp ? EPS_SMALL : - EPS_SMALL) ; // Elimino le soluzioni cha danno intersezioni fuori dai limiti in Z del cilindro if ( nRoot == 2) { double dIntZ2 = ptP.z + vdRoots[1] * vtV.z ; if ( dIntZ2 < 0 + dEpsLow || dIntZ2 > dH + dEpsUp) nRoot = 1 ; } if ( nRoot >= 1) { double dIntZ1 = ptP.z + vdRoots[0] * vtV.z ; if ( dIntZ1 < 0 + dEpsLow || dIntZ1 > dH + dEpsUp) { if ( nRoot == 2) vdRoots[0] = vdRoots[1] ; -- nRoot ; } } // Due soluzioni: la retta interseca due volte la superficie laterale if ( nRoot == 2) { // Punti di intersezione con la superficie del cilindro ptInt1 = ptP + vdRoots[0] * vtV ; ptInt2 = ptP + vdRoots[1] * vtV ; if ( ptInt1.z > ptInt2.z) swap( ptInt1, ptInt2) ; // Determino le normali vtN1.Set( -ptInt1.x, -ptInt1.y, 0) ; vtN1.Normalize() ; vtN2.Set( -ptInt2.x, -ptInt2.y, 0) ; vtN2.Normalize() ; // Porto i punti e i versori nel riferimento globale ptInt1.ToGlob( CylFrame) ; vtN1.ToGlob( CylFrame) ; ptInt2.ToGlob( CylFrame) ; vtN2.ToGlob( CylFrame) ; // Trovate intersezioni return true ; } // Una soluzione : la retta interseca la superficie laterale e un piano else if ( nRoot == 1) { // Se piano superiore if ( nBasInt == 2) { // Punto di intersezione ptInt1 = ptP + vdRoots[0] * vtV ; // Normale alla superficie del cilindro verso l'interno vtN1.Set( -ptInt1.x, -ptInt1.y, 0) ; vtN1.Normalize() ; } // altrimenti piano inferiore else if ( nBasInt == 1) { // Punto di intersezione ptInt2 = ptP + vdRoots[0] * vtV ; // Normale alla superficie del cilindro verso l'interno vtN2.Set( -ptInt2.x, -ptInt2.y, 0) ; vtN2.Normalize() ; } // altrimenti niente else return false ; // Porto i punti e i versori nel riferimento globale ptInt1.ToGlob( CylFrame) ; vtN1.ToGlob( CylFrame) ; ptInt2.ToGlob( CylFrame) ; vtN2.ToGlob( CylFrame) ; // Trovate intersezioni return true ; } // Nessuna soluzione : nessuna intersezione else return false ; } //---------------------------------------------------------------------------- // Riferimento con origine nel vertice del cono e asse di simmetria coincidente con l'asse Z. // La funzione restituisce true in caso di intersezione, false altrimenti. //---------------------------------------------------------------------------- bool VolZmap::IntersLineConus( const Point3d& ptLineSt, const Vector3d& vtLineDir, const Frame3d& ConusFrame, double dTan, double dMinH, double dMaxH, bool bTapLow, bool bTapUp, Point3d& ptInt1, Vector3d& vtN1, Point3d& ptInt2, Vector3d& vtN2) { // Porto la linea nel riferimento del cono Point3d ptP = ptLineSt ; ptP.ToLoc( ConusFrame) ; Vector3d vtV = vtLineDir ; vtV.ToLoc( ConusFrame) ; // Raggi delle due basi double dMinRad = dTan * dMinH ; double dMaxRad = dTan * dMaxH ; // Epsilon per piani di tappo double dEpsLow = ( bTapLow ? - EPS_SMALL : EPS_SMALL) ; double dEpsUp = ( bTapUp ? EPS_SMALL : - EPS_SMALL) ; // Determino le eventuali intersezioni con le due basi a quota minima e massima (solo se linea non parallela ad esse) int nBasInt = 0 ; if ( abs( vtV.z) > EPS_ZERO) { ptInt1 = ptP + ( ( dMinH - ptP.z) / vtV.z) * vtV ; if ( ptInt1.x * ptInt1.x + ptInt1.y * ptInt1.y < dMinRad * dMinRad + 2 * dMinRad * dTan * dEpsLow) { nBasInt += 1 ; vtN1 = Z_AX ; } ptInt2 = ptP + ( ( dMaxH - ptP.z) / vtV.z) * vtV ; if ( ptInt2.x * ptInt2.x + ptInt2.y * ptInt2.y < dMaxRad * dMaxRad + 2 * dMaxRad * dTan * dEpsUp) { nBasInt += 2 ; vtN2 = - Z_AX ; } } // Se la linea interseca entrambe le basi, si sono trovate le due intersezioni if ( nBasInt == 3) { // Porto i punti e i versori nel riferimento globale ptInt1.ToGlob( ConusFrame) ; vtN1.ToGlob( ConusFrame) ; ptInt2.ToGlob( ConusFrame) ; vtN2.ToGlob( ConusFrame) ; // Trovate intersezioni return true ; } // Determino le intersezioni con la superficie laterale del cono DBLVECTOR vdCoef( 3) ; double dSqTan = dTan * dTan ; vdCoef[0] = dSqTan * ptP.z * ptP.z - ptP.x * ptP.x - ptP.y * ptP.y ; vdCoef[1] = 2 * ( dSqTan * ptP.z * vtV.z - ptP.x * vtV.x - ptP.y * vtV.y) ; vdCoef[2] = dSqTan * vtV.z * vtV.z - vtV.x * vtV.x - vtV.y * vtV.y ; DBLVECTOR vdRoots ; int nRoot = PolynomialRoots( 2, vdCoef, vdRoots) ; // Elimino le soluzioni cha danno intersezioni fuori dai limiti in Z del tronco if ( nRoot == 2) { double dIntZ2 = ptP.z + vdRoots[1] * vtV.z ; if ( dIntZ2 < dMinH + dEpsLow || dIntZ2 > dMaxH + dEpsUp) nRoot = 1 ; } if ( nRoot >= 1) { double dIntZ1 = ptP.z + vdRoots[0] * vtV.z ; if ( dIntZ1 < dMinH + dEpsLow || dIntZ1 > dMaxH + dEpsUp) { if ( nRoot == 2) vdRoots[0] = vdRoots[1] ; -- nRoot ; } } // Due soluzioni: la retta interseca due volte la superficie laterale if ( nRoot == 2) { // Punti di intersezione con la superficie del cono ptInt1 = ptP + vdRoots[0] * vtV ; ptInt2 = ptP + vdRoots[1] * vtV ; if ( ptInt1.z > ptInt2.z) swap( ptInt1, ptInt2) ; // Determino le normali vtN1.Set( -ptInt1.x, -ptInt1.y, ( ptInt1.x * ptInt1.x + ptInt1.y * ptInt1.y) / ptInt1.z) ; vtN1.Normalize() ; vtN2.Set( -ptInt2.x, -ptInt2.y, ( ptInt2.x * ptInt2.x + ptInt2.y * ptInt2.y) / ptInt2.z) ; vtN2.Normalize() ; // Porto i punti e i versori nel riferimento globale ptInt1.ToGlob( ConusFrame) ; vtN1.ToGlob( ConusFrame) ; ptInt2.ToGlob( ConusFrame) ; vtN2.ToGlob( ConusFrame) ; // Trovate intersezioni return true ; } // Una soluzione : la retta interseca la superficie laterale e un piano else if ( nRoot == 1) { // Se piano superiore if ( nBasInt == 2) { // Punto di intersezione ptInt1 = ptP + vdRoots[0] * vtV ; // Normale alla superficie del cono verso l'interno vtN1.Set( -ptInt1.x, -ptInt1.y, ( ptInt1.x * ptInt1.x + ptInt1.y * ptInt1.y) / ptInt1.z) ; vtN1.Normalize() ; } // altrimenti piano inferiore else if( nBasInt == 1) { // Punto di intersezione ptInt2 = ptP + vdRoots[0] * vtV ; // Normale alla superficie del cono verso l'interno vtN2.Set( -ptInt2.x, -ptInt2.y, ( ptInt2.x * ptInt2.x + ptInt2.y * ptInt2.y) / ptInt2.z) ; vtN2.Normalize() ; } // altrimenti niente else return false ; // Porto i punti e i versori nel riferimento globale ptInt1.ToGlob( ConusFrame) ; vtN1.ToGlob( ConusFrame) ; ptInt2.ToGlob( ConusFrame) ; vtN2.ToGlob( ConusFrame) ; // Trovate intersezioni return true ; } // Nessuna soluzione : nessuna intersezione else return false ; } //---------------------------------------------------------------------------- // L'origine del sistema di riferimento deve essere // nel centro della circonferenza di base, la cui traslazione obliqua // genera il cilindro ellittico, e l'asse z deve essere l'asse // di simmetria di tale circonferenza. // La funzione restituisce true in caso di intersezione, false altrimenti. // NB: dSqRad è il quadrato del raggio della circonferenza la cui // traslazione obliqua genera il cilindro ellittico, dLongMvLen e // dOrtMvLen sono rispettivamente le lunghezze delle proiezioni del // movimento su z e x del sistema di riferimento CircFrame. bool VolZmap::IntersLineEllipticalCylinder( const Point3d& ptLineSt, const Vector3d& vtLineDir, const Frame3d& CircFrame, double dRad, double dLongMvLen, double dOrtMvLen, bool bTapLow, bool bTapUp, Point3d& ptInt1, Vector3d& vtN1, Point3d& ptInt2, Vector3d& vtN2) { // Se il cilindrico ellittico degenera in un piano, non bisogna tagliare if ( abs( dLongMvLen) < EPS_SMALL) return false ; // Porto la linea nel riferimento del cilindro Point3d ptP = ptLineSt ; ptP.ToLoc( CircFrame) ; Vector3d vtV = vtLineDir ; vtV.ToLoc( CircFrame) ; // Quadrato del raggio double dSqRad = dRad * dRad ; // Asse cilindro ellittico Vector3d vtAx( dOrtMvLen, 0, dLongMvLen) ; vtAx.Normalize() ; // Retta parallela all'asse del cilindro if ( AreSameOrOppositeVectorApprox( vtV, vtAx)) { // Interseco la retta con i piani delle circonferenze Point3d ptOLsCirc( dOrtMvLen, 0, dLongMvLen) ; ptInt1 = ptP - ( ptP.z / vtV.z) * vtV ; ptInt2 = ptP - ( ( ptP.z - dLongMvLen) / vtV.z) * vtV ; double dSafeSqRad = dSqRad - 2 * dRad * EPS_SMALL ; if ( ( ptInt1 - ORIG).SqLenXY() < dSafeSqRad && ( ptInt2 - ptOLsCirc).SqLenXY() < dSafeSqRad) { vtN1 = Z_AX ; vtN2 = - Z_AX ; ptInt1.ToGlob( CircFrame) ; ptInt2.ToGlob( CircFrame) ; vtN1.ToGlob( CircFrame) ; vtN2.ToGlob( CircFrame) ; return true ; } else return false ; } // Coefficiente angolare della retta di movimento nel // piano ZX del sistema di riferimento del movimento e suo quadrato double dObCoef = dOrtMvLen / dLongMvLen ; double dSqCoef = dObCoef * dObCoef ; // Setto i coeficienti dell'equazione DBLVECTOR vdCoef(3) ; vdCoef[0] = dSqCoef * ptP.z * ptP.z + ptP.x * ptP.x + ptP.y * ptP.y - 2 * dObCoef * ptP.z * ptP.x - dSqRad ; vdCoef[1] = 2 * ( dSqCoef * vtV.z * ptP.z + vtV.x * ptP.x + vtV.y * ptP.y - dObCoef * ( vtV.z * ptP.x + vtV.x * ptP.z)) ; vdCoef[2] = dSqCoef * vtV.z * vtV.z + vtV.x * vtV.x + vtV.y * vtV.y - 2 * dObCoef * vtV.z * vtV.x ; // Numero di soluzioni DBLVECTOR vdRoots ; int nRoot = PolynomialRoots( 2, vdCoef, vdRoots) ; // L'equazione ammette o due soluzioni (eventualmente // coincidenti) oppure nessuna o infinite se la la retta // appartiene alla superficie // Se numero di soluzioni diverso da due le eventuali intersezioni sono già state trovate if ( nRoot != 2) return false ; // Due soluzioni trovate // Vettore di movimento Vector3d vtMv( dOrtMvLen, 0, dLongMvLen) ; // Punti di intersezione ptInt1 = ptP + vdRoots[0] * vtV ; ptInt2 = ptP + vdRoots[1] * vtV ; // Simmetria del problema if ( ptInt1.z > ptInt2.z) swap( ptInt1, ptInt2) ; // Determino le normali alla superficie nei punti d'intersezione // Vettore 1 Vector3d vtTest1 = ( ptInt1 - ORIG) - ( ptInt1 - ORIG) * vtAx * vtAx ; double dX0_1 = ( vtTest1.x > 0 ? 1 : -1) * sqrt( max( dSqRad - ptInt1.y * ptInt1.y, 0.)) ; Vector3d vtCirc1( - dX0_1, - ptInt1.y, 0) ; Vector3d vtTan1( vtCirc1.y, - vtCirc1.x, 0) ; Vector3d vtCross1 = vtTan1 ^ vtMv ; vtN1 = ( vtCross1 * vtCirc1 > - EPS_ZERO ? vtCross1 : - vtCross1) ; vtN1.Normalize() ; // Vettore 2 Vector3d vtTest2 = ( ptInt2 - ORIG) - ( ptInt2 - ORIG) * vtAx * vtAx ; double dX0_2 = ( vtTest2.x > 0 ? 1 : -1) * sqrt( max( dSqRad - ptInt2.y * ptInt2.y, 0.)) ; Vector3d vtCirc2( - dX0_2, - ptInt2.y, 0) ; Vector3d vtTan2( vtCirc2.y, - vtCirc2.x, 0) ; Vector3d vtCross2 = vtTan2 ^ vtMv ; vtN2 = ( vtCross2 * vtCirc2 > - EPS_ZERO ? vtCross2 : - vtCross2) ; vtN2.Normalize() ; // Flag per i tappi double dEpsLow = ( bTapLow ? - EPS_SMALL : EPS_SMALL) ; double dEpsUp = ( bTapUp ? EPS_SMALL : - EPS_SMALL) ; // Studio le soluzioni: se una è fuori dalla regione // ammissibile, vuol dire che la retta esce da un tappo. if ( ptInt1.z < dLongMvLen + dEpsUp) { if ( ptInt1.z > + dEpsLow) { // ptInt2 è sul tappo if ( ptInt2.z > dLongMvLen + dEpsUp) { ptInt2 = ptP + ( ( dLongMvLen - ptP.z) / vtV.z) * vtV ; vtN2 = - Z_AX ; } } else { // Entrambe le soluzioni sono su un tappo if ( ptInt2.z > dLongMvLen + dEpsUp) { ptInt1 = ptP - ( ptP.z / vtV.z) * vtV ; ptInt2 = ptP + ( ( dLongMvLen - ptP.z) / vtV.z) * vtV ; vtN1.Set( 0, 0, 1) ; vtN2.Set( 0, 0, -1) ; if ( ptInt1.x * ptInt1.x + ptInt1.y * ptInt1.y > dSqRad && ptInt2.x * ptInt2.x + ptInt2.y * ptInt2.y > dSqRad) return false ; } // La prima soluzione è sul tappo else if ( ptInt2.z > dEpsLow) { ptInt1 = ptP - ( ptP.z / vtV.z) * vtV ; vtN1.Set( 0, 0, 1) ; } else return false ; } } else return false ; // Porto i punti e i versori nel riferimento globale ptInt1.ToGlob( CircFrame) ; vtN1.ToGlob( CircFrame) ; ptInt2.ToGlob( CircFrame) ; vtN2.ToGlob( CircFrame) ; return true ; } //---------------------------------------------------------------------------- bool VolZmap::IntersLineMyPolyhedron( const Point3d& ptLineSt, const Vector3d& vtLineDir, const Frame3d& PolyFrame, double dLenX, double dLenY, double dLenZ, double dDeltaZ, Point3d& ptInt1, Vector3d& vtN1, Point3d& ptInt2, Vector3d& vtN2) { // Controllo sulle dimensioni lineari affinché sia valido il poliedro if ( dLenX <= EPS_SMALL || dLenY <= EPS_SMALL || dLenZ <= EPS_SMALL) return false ; // Porto la linea nel riferimento del poliedro Point3d ptP = ptLineSt ; ptP.ToLoc( PolyFrame) ; Vector3d vtV = vtLineDir ; vtV.ToLoc( PolyFrame) ; // Facce 1 e 2 parallele a XY // Facce 3 e 4 parallele a XZ // Facce 5 e 6 oblique Point3d ptFacet135( 0, dLenY / 2, 0) ; Point3d ptFacet246( dLenX, - dLenY / 2, dLenZ + dDeltaZ) ; // Vettori per descrizione piani obliqui Vector3d vtFacet5 = ptFacet135 - ptP ; Vector3d vtFacet6 = ptFacet246 - ptP ; Vector3d vtOb( - dDeltaZ, 0, dLenX) ; vtOb.Normalize() ; // Controlli affinché non vengano tagliati dexel a filo // con il passaggio dell'utensile: // Controllo sulle facce 1 e 2 if ( abs( vtV.y) < EPS_ZERO && abs( ptP.y) > ptFacet135.y) return false ; // Controllo sulle facce 3 e 4 if ( abs( vtV.x) < EPS_ZERO && ( ptP.x < ptFacet135.x || ptP.x > ptFacet246.x)) return false ; // Controllo sulle facce 5 e 6 double dP1 = abs ( ( ptFacet135 - ptP) * vtOb) ; double dP2 = abs ( ( ptFacet246 - ptP) * vtOb) ; if ( abs( vtV * vtOb) < EPS_ZERO && ( dP1 < EPS_SMALL || dP2 < EPS_SMALL)) return false ; // Punti notevoli Point3d ptI1 = ptP + ( ( ptFacet135.y - ptP.y) / vtV.y) * vtV ; Point3d ptI2 = ptP + ( ( ptFacet246.y - ptP.y) / vtV.y) * vtV ; Point3d ptI3 = ptP + ( ( ptFacet135.x - ptP.x) / vtV.x) * vtV ; Point3d ptI4 = ptP + ( ( ptFacet246.x - ptP.x) / vtV.x) * vtV ; Point3d ptI5 = ptP + ( ( vtFacet5 * vtOb) / ( vtV * vtOb)) * vtV ; Point3d ptI6 = ptP + ( ( vtFacet6 * vtOb) / ( vtV * vtOb)) * vtV ; // Ricerca intersezioni con le facce int nIntNum = 0 ; // Intersezione con la prima faccia if ( ptI1.x >= 0 && ptI1.x <= dLenX && ptI1.z * dLenX >= dDeltaZ * ptI1.x && ( ptI1.z - dLenZ) * dLenX <= dDeltaZ * ptI1.x) { ptInt1 = ptI1 ; vtN1 = - Y_AX ; ++ nIntNum ; } // Intersezione con la seconda faccia if ( ptI2.x >= 0 && ptI2.x <= dLenX && ptI2.z * dLenX >= dDeltaZ * ptI2.x && ( ptI2.z - dLenZ) * dLenX <= dDeltaZ * ptI2.x) { if ( nIntNum == 0) { ptInt1 = ptI2 ; vtN1 = Y_AX ; ++ nIntNum ; } else if ( ( ptInt1 - ptI2).SqLen() > SQ_EPS_SMALL) { ptInt2 = ptI2 ; vtN2 = Y_AX ; ++ nIntNum ; } } // Intersezione con la terza faccia if ( nIntNum < 2 && ptI3.z >= 0 && ptI3.z <= dLenZ && ptI3.y >= - ptFacet135.y && ptI3.y <= ptFacet135.y) { if ( nIntNum == 0) { ptInt1 = ptI3 ; vtN1 = X_AX ; ++ nIntNum ; } else if ( ( ptInt1 - ptI3).SqLen() > SQ_EPS_SMALL) { ptInt2 = ptI3 ; vtN2 = X_AX ; ++ nIntNum ; } } // Intersezione con la quarta faccia if ( nIntNum < 2 && ptI4.z >= dDeltaZ && ptI4.z <= dLenZ + dDeltaZ && ptI4.y >= - ptFacet135.y && ptI4.y <= ptFacet135.y) { if ( nIntNum == 0) { ptInt1 = ptI4 ; vtN1 = - X_AX ; ++ nIntNum ; } else if ( ( ptInt1 - ptI4).SqLen() > SQ_EPS_SMALL) { ptInt2 = ptI4 ; vtN2 = - X_AX ; ++ nIntNum ; } } // Intersezione con la quinta faccia if ( nIntNum < 2 && ptI5.x >= 0 && ptI5.x <= dLenX && ptI5.y >= - ptFacet135.y && ptI5.y <= ptFacet135.y) { if ( nIntNum == 0) { ptInt1 = ptI5 ; vtN1 = vtOb ; ++ nIntNum ; } else if ( ( ptInt1 - ptI5).SqLen() > SQ_EPS_SMALL) { ptInt2 = ptI5 ; vtN2 = vtOb ; ++ nIntNum ; } } // Intersezione con la sesta faccia if ( nIntNum < 2 && ptI6.x >= 0 && ptI6.x <= dLenX && ptI6.y >= - ptFacet135.y && ptI6.y <= ptFacet135.y) { if ( nIntNum == 0) { ptInt1 = ptI6; vtN1 = - vtOb ; ++ nIntNum ; } else if ( ( ptInt1 - ptI6).SqLen() > SQ_EPS_SMALL) { ptInt2 = ptI6; vtN2 = - vtOb ; ++ nIntNum ; } } if ( nIntNum != 2) return false ; // Porto i punti e i versori nel riferimento globale ptInt1.ToGlob( PolyFrame) ; ptInt2.ToGlob( PolyFrame) ; vtN1.ToGlob( PolyFrame) ; vtN2.ToGlob( PolyFrame) ; return true ; } //---------------------------------------------------------------------------- bool VolZmap::GetVolume( double& dVol) const { // verifico lo stato if ( m_nStatus != OK) return false ; // Eseguo il calcolo della lunghezza totale degli spilloni double dLen = 0 ; for ( int nMap = 0 ; nMap < int( m_nMapNum) ; ++ nMap) { for ( int nDex = 0 ; nDex < int( m_Values[nMap].size()) ; ++ nDex) { for ( int nInt = 0 ; nInt < int( m_Values[nMap][nDex].size()) ; ++ nInt) { dLen += m_Values[nMap][nDex][nInt].dMax - m_Values[nMap][nDex][nInt].dMin ; } } } // Il volume si trova moltiplicando per l'area di ogni spillone e dividendo per il numero di mappe dVol = dLen * ( m_dStep * m_dStep) / m_nMapNum ; return true ; } //---------------------------------------------------------------------------- bool VolZmap::GetPartVolume( int nPart, double& dVol) const { // verifico lo stato if ( m_nStatus != OK) return false ; // Se una sola mappa o il numero di componenti è indefinito, errore if ( m_nMapNum == 1 || m_nConnectedCompoCount == -1) return false ; // Se la componente richiesta non esiste, errore if ( nPart < 0 || nPart > m_nConnectedCompoCount - 1) return false ; // Eseguo il calcolo della lunghezza totale degli spilloni double dLen = 0 ; for ( int nMap = 0 ; nMap < int( m_nMapNum) ; ++ nMap) { for ( int nDex = 0 ; nDex < int( m_Values[nMap].size()) ; ++ nDex) { for ( int nInt = 0 ; nInt < int( m_Values[nMap][nDex].size()) ; ++ nInt) { // se il pezzo di spillone appartiene alla parte, ne aggiungo la lunghezza if ( m_Values[nMap][nDex][nInt].nCompo == nPart + 1) { dLen += m_Values[nMap][nDex][nInt].dMax - m_Values[nMap][nDex][nInt].dMin ; } } } } // Il volume si trova moltiplicando per l'area di ogni spillone e dividendo per il numero di mappe dVol = dLen * ( m_dStep * m_dStep) / 3 ; return true ; } //---------------------------------------------------------------------------- // Il piano è espresso nel sistema locale, la funzione lo esprime nel sistema Zmap. // I loop della flat region ottenuta dall'intersezione vengono salvati nello // standard vector vLoop passato per riferimento. // Se il processo va a buon fine viene restituito true, false altrimenti. bool VolZmap::GetPlaneIntersection( const Plane3d& plPlane, ICURVEPOVECTOR& vpLoop) const { // Verifico validità parametri if ( ! plPlane.IsValid()) return false ; // Porto il piano nel sistema Zmap Plane3d plPlaneL = plPlane ; plPlaneL.ToLoc( m_MapFrame) ; // Vettore del piano Vector3d vtNL = plPlaneL.GetVersN() ; // Se il piano non interseca il bounding box del solido, abbiamo finito if ( ! TestIntersPlaneZmapBBox( plPlaneL)) return true ; // Cerco la mappa con i dexel più perpendicolari al piano e valuto // se il versore del piano è equiverso o controverso ai dexel. int nGrid = 0 ; int nSign = ( vtNL.z > 0 ? 1 : - 1) ; if ( abs( vtNL.x) > abs( vtNL.y) && abs( vtNL.x) > abs( vtNL.z)) { nGrid = 1 ; nSign = ( vtNL.x > 0 ? 1 : - 1) ; } else if ( abs( vtNL.y) > abs( vtNL.x) && abs( vtNL.y) > abs( vtNL.z)) { nGrid = 2 ; nSign = ( vtNL.y > 0 ? 1 : - 1) ; } // Ciclo sulle celle vector vLine ; for ( int ni = - 1 ; ni < int( m_nNx[nGrid]) ; ++ ni) { for ( int nj = - 1 ; nj < int( m_nNy[nGrid]) ; ++ nj) { ProcessCell( nGrid, ni, nj, plPlaneL, vLine) ; } } // Creo i loop ChainCurves LoopCreator ; LoopCreator.Init( false, EPS_SMALL, int( vLine.size())) ; // Carico le curve per concatenarle for ( int nCv = 0 ; nCv < int( vLine.size()) ; ++ nCv) { Point3d ptSt = vLine[nCv].GetStart() ; Point3d ptEn = vLine[nCv].GetEnd() ; Vector3d vtDir ; vLine[nCv].GetStartDir( vtDir) ; LoopCreator.AddCurve( nCv + 1, ptSt, vtDir, ptEn, vtDir) ; } // Recupero i concatenamenti INTVECTOR vIds ; Point3d ptNearStart ; while ( LoopCreator.GetChainFromNear( ptNearStart, false, vIds)) { PtrOwner pLoop( CreateBasicCurveComposite()) ; if ( IsNull( pLoop)) return false ; for ( auto i : vIds) { // Aggiungo la linea alla curva composta. if ( ! pLoop->AddCurve( vLine[i-1], true, EPS_SMALL)) return false ; } pLoop->SetExtrusion( vtNL) ; pLoop->ToGlob( m_MapFrame) ; // Se normali controverse, devo invertire il verso del loop. if ( nSign < 0) pLoop->Invert() ; // eseguo approssimazione PolyLine PL ; if ( ! pLoop->ApproxWithLines( m_dStep, ANG_TOL_STD_DEG, ICurve::APL_RIGHT, PL) || ! pLoop->Clear() || ! pLoop->FromPolyLine( PL)) return false ; pLoop->MergeCurves( m_dStep, ANG_TOL_STD_DEG) ; // Inserisco la curva composita nella raccolta da ritornare vpLoop.emplace_back( Release( pLoop)) ; } return true ; } //---------------------------------------------------------------------------- // Data una cella, identificata dal numero della sua griglia e dai suoi indici i,j // all'interno della griglia, genera gli eventuali segmenti che costituiscono la curva // corrispondente alla frontiera della regione di piano interna la solido. // Se la griglia o la cella non è valida, la funzione restituisce false, // altrimenti valuta la tipologia della cella e crea gli eventuali segmenti. // Se il processo va a buon fine la funzione restituisce true. bool VolZmap::ProcessCell( int nGrid, int nCellI, int nCellJ, const Plane3d& plPlane, vector& vLine) const { // Se la griglia non esiste vi è un errore if ( nGrid < 0 || nGrid > 2) return false ; // Se la cella non esiste vi è un errore if ( nCellI < - 1 || nCellI >= int( m_nNx[nGrid]) || nCellJ < - 1 || nCellJ >= int( m_nNy[nGrid])) return false ; // Determino la configurazione della cella int nIndex = CalcIndexForPlaneCells( plPlane, nGrid, nCellI, nCellJ) ; // Tabella segmenti static int nLineTable[16][5] = { { -1, -1, -1, -1, -1}, { 0, 3, -1, -1, -1}, { 1, 0, -1, -1, -1}, { 1, 3, -1, -1, -1}, { 2, 1, -1, -1, -1}, { 0, 1, 2, 3, -1}, { 2, 0, -1, -1, -1}, { 2, 3, -1, -1, -1}, { 3, 2, -1, -1, -1}, { 0, 2, -1, -1, -1}, { 1, 2, 3, 0, -1}, { 1, 2, -1, -1, -1}, { 3, 1, -1, -1, -1}, { 0, 1, -1, -1, -1}, { 3, 0, -1, -1, -1}, { -1, -1, -1, -1, -1} } ; // Tabella dei punti medi dei segmenti: le righe rappresentano // il numero di segmento (0, 1, 2, 3) le colonne rappresentano // i e j. static double dDeltaIJ[4][2] = { { 0.5, 0}, { 1, 0.5}, { 0.5, 1}, { 0, 0.5} } ; // Se la cella è di frontiera costruisco i segmenti della curva if ( nIndex != 0 && nIndex != 15) { // Ciclo su tutti gli spigoli della cella che vengono attraversati dalla curva. for ( int nEdge = 0 ; nLineTable[nIndex][nEdge] != -1 ; nEdge += 2) { int nStEdge = nLineTable[nIndex][nEdge] ; int nEnEdge = nLineTable[nIndex][nEdge + 1] ; double nStDeltaI = dDeltaIJ[nStEdge][0] ; double nStDeltaJ = dDeltaIJ[nStEdge][1] ; double nEnDeltaI = dDeltaIJ[nEnEdge][0] ; double nEnDeltaJ = dDeltaIJ[nEnEdge][1] ; if ( nGrid == 0) { // Punti sulla griglia Point3d ptGrSt( ( nCellI + nStDeltaI + 0.5) * m_dStep, ( nCellJ + nStDeltaJ + 0.5) * m_dStep, 0) ; Point3d ptGrEn( ( nCellI + nEnDeltaI + 0.5) * m_dStep, ( nCellJ + nEnDeltaJ + 0.5) * m_dStep, 0) ; // Corrispondenti punti sul piano Point3d ptSt, ptEn ; if ( IntersLinePlane( ptGrSt, Z_AX, 10, plPlane, ptSt, false) && IntersLinePlane( ptGrEn, Z_AX, 10, plPlane, ptEn, false)) { // Costruisco il tratto di curva CurveLine cvLine ; cvLine.Set( ptSt, ptEn) ; vLine.emplace_back( cvLine) ; } } else if ( nGrid == 1) { // Punti sulla griglia Point3d ptGrSt( 0, ( nCellI + nStDeltaI + 0.5) * m_dStep, ( nCellJ + nStDeltaJ + 0.5) * m_dStep) ; Point3d ptGrEn( 0, ( nCellI + nEnDeltaI + 0.5) * m_dStep, ( nCellJ + nEnDeltaJ + 0.5) * m_dStep) ; // Corrispondenti punti sul piano Point3d ptSt, ptEn ; if ( IntersLinePlane( ptGrSt, X_AX, 10, plPlane, ptSt, false) && IntersLinePlane( ptGrEn, X_AX, 10, plPlane, ptEn, false)) { // Costruisco il tratto di curva CurveLine cvLine ; cvLine.Set( ptSt, ptEn) ; vLine.emplace_back( cvLine) ; } } else { // Punti sulla griglia Point3d ptGrSt( ( nCellJ + nStDeltaJ + 0.5) * m_dStep, 0, ( nCellI + nStDeltaI + 0.5) * m_dStep) ; Point3d ptGrEn( ( nCellJ + nEnDeltaJ + 0.5) * m_dStep, 0, ( nCellI + nEnDeltaI + 0.5) * m_dStep) ; // Corrispondenti punti sul piano Point3d ptSt, ptEn ; if ( IntersLinePlane( ptGrSt, Y_AX, 10, plPlane, ptSt, false) && IntersLinePlane( ptGrEn, Y_AX, 10, plPlane, ptEn, false)) { // Costruisco il tratto di curva CurveLine cvLine ; cvLine.Set( ptSt, ptEn) ; vLine.emplace_back( cvLine) ; } } } } return true ; } //---------------------------------------------------------------------------- // Dato un dexel, identificato dal numero della sua griglia e dai suoi indici, // determina se un suo intervallo attraversa il piano dato. // Il numero di griglia e gli indici del dexel devono essere validi, se // non lo sono, la funzione restituisce false, altrimenti cerca un // tratto di dexel che intersechi il piano, se lo trova restituisce true, // false altrimenti. bool VolZmap::InOut( const Plane3d& plPlane, int nGrid, int nI, int nJ) const { // Se la griglia non esiste, vi è un errore. if ( nGrid < 0 || nGrid > 2) return false ; // Se gli indici sono di frontiera per lo // Zmap o non esistono, non sono interni. if ( nI <= - 1 || nI >= int( m_nNx[nGrid]) || nJ <= - 1 || nJ >= int( m_nNy[nGrid])) return false ; // Numero del dexel int nDex = nJ * m_nNx[nGrid] + nI ; // Ciclo sui segmenti del dexel bool bNotFound = true ; for ( int nK = 0 ; nK < int( m_Values[nGrid][nDex].size()) ; ++ nK) { if ( nGrid == 0) { // Punti estremi del segmento Point3d ptSt( ( nI + 0.5) * m_dStep, ( nJ + 0.5) * m_dStep, m_Values[nGrid][nDex][nK].dMin) ; Point3d ptEn( ( nI + 0.5) * m_dStep, ( nJ + 0.5) * m_dStep, m_Values[nGrid][nDex][nK].dMax) ; // Se il segmento interseca il piano abbiamo finito Point3d ptInt ; if ( IntersLinePlane( ptSt, ptEn, plPlane, ptInt)) return true ; } else if ( nGrid == 1) { // Punti estremi del segmento Point3d ptSt( m_Values[nGrid][nDex][nK].dMin, ( nI + 0.5) * m_dStep, ( nJ + 0.5) * m_dStep) ; Point3d ptEn( m_Values[nGrid][nDex][nK].dMax, ( nI + 0.5) * m_dStep, ( nJ + 0.5) * m_dStep) ; // Se il segmento interseca il piano abbiamo finito Point3d ptInt ; if ( IntersLinePlane( ptSt, ptEn, plPlane, ptInt)) return true ; } else { // Punti estremi del segmento Point3d ptSt( ( nJ + 0.5) * m_dStep, m_Values[nGrid][nDex][nK].dMin, ( nI + 0.5) * m_dStep) ; Point3d ptEn( ( nJ + 0.5) * m_dStep, m_Values[nGrid][nDex][nK].dMax, ( nI + 0.5) * m_dStep) ; // Se il segmento interseca il piano abbiamo finito Point3d ptInt ; if ( IntersLinePlane( ptSt, ptEn, plPlane, ptInt)) return true ; } } return false ; } //---------------------------------------------------------------------------- int VolZmap::CalcIndexForPlaneCells( const Plane3d& plPlane, int nGrid, int nCellI, int nCellJ) const { int nIndex = 0 ; if ( InOut( plPlane, nGrid, nCellI, nCellJ)) nIndex |= ( 1 << 0) ; if ( InOut( plPlane, nGrid, nCellI + 1, nCellJ)) nIndex |= ( 1 << 1) ; if ( InOut( plPlane, nGrid, nCellI + 1, nCellJ + 1)) nIndex |= ( 1 << 2) ; if ( InOut( plPlane, nGrid, nCellI, nCellJ + 1)) nIndex |= ( 1 << 3) ; return nIndex ; } //---------------------------------------------------------------------------- // dZ è una quota relativa alla mappa nGrid, se il dexel nDex esiste e // dZ è interna a un suo intervallo restituisce true, false altrimenti. bool VolZmap::IsZInsideInterval( int nGrid, int nDex, double dZ) const { // Se la griglia non esiste vi è un errore. if ( nGrid < 0 || nGrid > 2) return false ; // Se il dexel non esiste vi è un errore. if ( nDex < 0 && nDex > int( m_Values[nGrid].size() - 1)) return false ; // Valuto se dZ è interna a un intervallo. for ( int nk = 0 ; nk < int( m_Values[nGrid][nDex].size()) ; ++ nk) { double dZ1 = m_Values[nGrid][nDex][nk].dMin ; double dZ2 = m_Values[nGrid][nDex][nk].dMax ; // Se troviamo dZ in un intervallo abbiamo finito. if ( dZ > dZ1 && dZ < dZ2) return true ; } // dZ non sta in un intervallo. return false ; } //---------------------------------------------------------------------------- // Box e piano devono essere nello stesso riferimento. Il box deve essere axis aligned. bool TestIntersPlaneBox( const BBox3d& b3Box, const Plane3d& plPlane) { // Calcolo le distanze con segno dei punti dal piano Point3d ptE0 = b3Box.GetMin() ; double dDist0 = DistPointPlane( ptE0, plPlane) ; Point3d ptE1( b3Box.GetMax().x, b3Box.GetMin().y, b3Box.GetMin().z) ; double dDist1 = DistPointPlane( ptE1, plPlane) ; Point3d ptE2( b3Box.GetMax().x, b3Box.GetMax().y, b3Box.GetMin().z) ; double dDist2 = DistPointPlane( ptE2, plPlane) ; Point3d ptE3( b3Box.GetMin().x, b3Box.GetMax().y, b3Box.GetMin().z) ; double dDist3 = DistPointPlane( ptE3, plPlane) ; Point3d ptE4( b3Box.GetMin().x, b3Box.GetMin().y, b3Box.GetMax().z) ; double dDist4 = DistPointPlane( ptE4, plPlane) ; Point3d ptE5( b3Box.GetMax().x, b3Box.GetMin().y, b3Box.GetMax().z) ; double dDist5 = DistPointPlane( ptE5, plPlane) ; Point3d ptE6 = b3Box.GetMax() ; double dDist6 = DistPointPlane( ptE6, plPlane) ; Point3d ptE7( b3Box.GetMin().x, b3Box.GetMax().y, b3Box.GetMax().z) ; double dDist7 = DistPointPlane( ptE7, plPlane) ; // Distanze tutte positive if ( dDist0 > EPS_SMALL && dDist1 > EPS_SMALL && dDist2 > EPS_SMALL && dDist3 > EPS_SMALL && dDist4 > EPS_SMALL && dDist5 > EPS_SMALL && dDist6 > EPS_SMALL && dDist7 > EPS_SMALL) return false ; // Distanze tutte negative if ( dDist0 < - EPS_SMALL && dDist1 < - EPS_SMALL && dDist2 < - EPS_SMALL && dDist3 < - EPS_SMALL && dDist4 < - EPS_SMALL && dDist5 < - EPS_SMALL && dDist6 < - EPS_SMALL && dDist7 < - EPS_SMALL) return false ; // Il piano interseca il box return true ; } //---------------------------------------------------------------------------- // Il piano deve essere espresso nel sistema di riferimento intrinseco dello Zmap. bool VolZmap::TestIntersPlaneZmapBBox( const Plane3d& plPlane) const { BBox3d b3Zmap( ORIG, Point3d( m_nNx[0] * m_dStep, m_nNy[0] * m_dStep, m_dMaxZ[0])) ; return TestIntersPlaneBox( b3Zmap, plPlane) ; }