//---------------------------------------------------------------------------- // 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 "IntersLineSurfTm.h" #include "/EgtDev/Include/EgtNumUtils.h" #include "/EgtDev/Include/EGkIntersLineTria.h" #include "/EgtDev/Include/EGkIntersLinePlane.h" #include "/EgtDev/Include/EGkChainCurves.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, unsigned int nGrid, unsigned int nI, unsigned int nJ, double& dU1, double& dU2) const { // Determino l'indice del dexel e il numero di suoi intervalli unsigned int nDexelPos = nJ * m_nNx[nGrid] + nI ; unsigned int nDexelSize = unsigned 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 ( unsigned 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, unsigned int nGrid, unsigned int nI, unsigned int nJ, double& dU1, double& dU2) const { // Determino l'indice del dexel e il numero di suoi intervalli unsigned int nDexelPos = nJ * m_nNx[nGrid] + nI ; unsigned int nDexelSize = unsigned 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 ( unsigned 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 { // 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 ; } // 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)) { Voxel 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 ; } } // Triangoli di frontiera dei voxel TRIA3DEXLIST lstTria ; ExtMarchingCubes( vVox, lstTria, bEnh) ; // 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 ; // 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 < 4 ; ++ k) { Point3d ptC = ptO + i * m_dStep * vtX + j * m_dStep * vtY ; if ( k == 0) ptC += 0.1 * m_dStep * vtX + 0.1 * m_dStep * vtY ; else if ( k == 1) ptC += 0.9 * m_dStep * vtX + 0.1 * m_dStep * vtY ; else if ( k == 2) ptC += 0.9 * m_dStep * vtX + 0.9 * m_dStep * vtY ; else if ( k == 3) ptC += 0.1 * m_dStep * vtX + 0.9 * m_dStep * vtY ; double dZmin, dZmax ; if ( IntersLineBox( ptC, 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::IntersLineCylinder( const Point3d& ptLineSt, const Vector3d& vtLineDir, const Frame3d& CylFrame, double dL, double dR, Point3d& ptInt1, Point3d& ptInt2, Vector3d& vtN1, Vector3d& vtN2, bool bTapO, bool bTapL) { // NB: L'origine del sistema di riferimento deve essere // nel centro della circonferenza di base e l'asse di simmetria // deve coincidere con l'asse x. // La funzione restituisce true in caso di intersezione, // false altrimenti. Point3d ptP = ptLineSt ; Vector3d vtV = vtLineDir ; // Trasformazione delle coordinate: // l'asse del cilindro corrisponde con // l'asse x del sistema di riferimento ptP.ToLoc( CylFrame) ; vtV.ToLoc( CylFrame) ; DBLVECTOR vdCoef(3) ; DBLVECTOR vdRoots ; double dSqRad = dR * dR - 2 * dR * EPS_SMALL ; vdCoef[0] = ptP.y * ptP.y + ptP.z * ptP.z - dSqRad ; vdCoef[1] = 2 * ( ptP.y * vtV.y + ptP.z * vtV.z) ; vdCoef[2] = vtV.y * vtV.y + vtV.z * vtV.z ; // Computo radici int nRoot = PolynomialRoots( 2, vdCoef, vdRoots) ; // Nessuna soluzione if ( nRoot == 0 || nRoot == 1) { if ( abs( vtV.x) > EPS_ZERO) { ptInt1 = ptP - ( ptP.x / vtV.x) * vtV ; ptInt2 = ptP + ( ( dL - ptP.x) / vtV.x) * vtV ; vtN1 = X_AX ; vtN2 = - X_AX ; if ( ptInt1.y * ptInt1.y + ptInt1.z * ptInt1.z <= dSqRad && ptInt2.y * ptInt2.y + ptInt2.z * ptInt2.z <= dSqRad) { ptInt1.ToGlob( CylFrame) ; ptInt2.ToGlob( CylFrame) ; vtN1.ToGlob( CylFrame) ; vtN2.ToGlob( CylFrame) ; return true ; } // Nessuna intersezione else return false ; } // Nessuna intersezione else return false ; } // L'equazione ammette o due soluzioni (eventualmente // coincidenti) oppure nessuna o infinite se la la retta // appartiene alla superficie if ( nRoot == 2) { double dEpsO = ( bTapO ? - EPS_SMALL : EPS_SMALL) ; double dEpsL = ( bTapL ? EPS_SMALL : - EPS_SMALL) ; ptInt1 = ptP + vdRoots[0] * vtV ; ptInt2 = ptP + vdRoots[1] * vtV ; if ( ptInt1.x > ptInt2.x) swap( ptInt1, ptInt2) ; vtN1.Set( 0, ( ORIG - ptInt1).y, ( ORIG - ptInt1).z) ; vtN2.Set( 0, ( ORIG - ptInt2).y, ( ORIG - ptInt2).z) ; if ( ptInt1.x < dL + dEpsL) { if ( ptInt1.x > dEpsO) { if ( ptInt2.x > dL + dEpsL) { ptInt2 = ptP + ( ( dL - ptP.x) / vtV.x) * vtV ; vtN2.Set( -1, 0, 0) ; } } else { if ( ptInt2.x > dL + dEpsL) { ptInt1 = ptP - ( ptP.x / vtV.x) * vtV ; ptInt2 = ptP + ( ( dL - ptP.x) / vtV.x) * vtV ; vtN1.Set( 1, 0, 0) ; vtN2.Set( -1, 0, 0) ; } else if ( ptInt2.x > dEpsO) { ptInt1 = ptP - ( ptP.x / vtV.x) * vtV ; vtN1.Set( 1, 0, 0) ; } else return false ; } } else return false ; // Riporto le coordinate nel sistema di riferimento griglia ptInt1.ToGlob( CylFrame) ; ptInt2.ToGlob( CylFrame) ; vtN1.ToGlob( CylFrame) ; vtN2.ToGlob( CylFrame) ; vtN1.Normalize() ; vtN2.Normalize() ; } return true ; } //---------------------------------------------------------------------------- bool VolZmap::IntersZLineCylinder( const Point3d& ptLine, const Point3d& ptBase, const Point3d& ptTop, double dCylR, double& dInfZ, double& dSupZ) { // NB: Le coordinate sono espresse nel sistema griglia // La funzione restituisce true in caso di intersezione, // false altrimenti. double dSqRad = dCylR * dCylR ; // Cilindro verticale if ( AreSamePointXYApprox( ptBase, ptTop)) { // Intersezione if ( SqDistXY( ptLine, ptBase) < dSqRad) { dInfZ = min( ptBase.z, ptTop.z) ; dSupZ = max( ptBase.z, ptTop.z) ; return true ; } // Non vi è intersezione else return false ; } // Cilindro non verticale else { // Studio delle simmetrie Point3d ptS = ( ptBase.z < ptTop.z ? ptBase : ptTop) ; Point3d ptE = ( ptBase.z < ptTop.z ? ptTop : ptBase) ; Vector3d vtAx = ptE - ptS ; Vector3d vtV1( vtAx.x, vtAx.y, 0) ; double dLenXY = vtV1.LenXY() ; double dSZ = ptS.z ; double dEZ = ptE.z ; double dDeltaZ = dEZ - dSZ ; Vector3d vtL( ptLine.x - ptS.x, ptLine.y - ptS.y, 0) ; // vtV1 e vtV2 formano un sistema ortonormale // sul piano e insieme a ptSxy formano un sistema // di riferimento bidimensionale vtV1.Normalize() ; Vector3d vtV2 = vtV1 ; vtV2.Rotate( Z_AX, 90) ; double dLen = vtAx.Len() ; // Sono seno e coseno dell'angolo complementare // rispetto a quello formato dal vettore movimento // con il piano, per questo motivo si ha dCos con // dDeltaZ e dSin con dLenXY double dCos = dDeltaZ / dLen ; double dSin = dLenXY / dLen ; // Nuove coordinate piane del punto double dLocX1 = vtL * vtV1 ; double dLocX2 = vtL * vtV2 ; double dSqRoot = sqrt( dSqRad - dLocX2 * dLocX2) ; double dX1_0 = dCos * dSqRoot ; if ( dLocX1 >= - dX1_0 && dLocX1 <= dLenXY + dX1_0 && abs( dLocX2) < dCylR) { // Minimi if ( dLocX1 < dX1_0) { double dDotS = vtAx * ( ptS - ORIG) ; // Qui usiamo ptLine perché servono coordinate griglia dInfZ = ( dDotS - vtAx.x * ptLine.x - vtAx.y * ptLine.y) / vtAx.z ; } else { double dZ0 = - dSin * dSqRoot ; dInfZ = dSZ + dZ0 + ( dLocX1 - dX1_0) * dDeltaZ / dLenXY ; } // Massimi if ( dLocX1 <= dLenXY - dX1_0) { double dZ0 = dSin * dSqRoot ; dSupZ = dSZ + dZ0 + ( dLocX1 + dX1_0) * dDeltaZ / dLenXY ; } else { double dDotE = vtAx * ( ptE - ORIG) ; // Qui usiamo ptLine perché servono coordinate griglia dSupZ = ( dDotE - vtAx.x * ptLine.x - vtAx.y * ptLine.y) / vtAx.z ; } return true ; } return false ; } } //---------------------------------------------------------------------------- bool VolZmap::IntersLineConus( const Point3d& ptLineSt, const Vector3d& vtLineDir, const Frame3d& ConusFrame, double dTan, double dl, double dL, Point3d& ptInt1, Point3d& ptInt2, Vector3d& vtN1, Vector3d& vtN2, bool bTapLow, bool bTapUp) { // NB: L'origine del sistema di riferimento deve essere // nel vertice del cono e l'asse di simmetria deve coincidere // con l'asse x. // La funzione restituisce true in caso di intersezione, // false altrimenti. Point3d ptP = ptLineSt ; Vector3d vtV = vtLineDir ; // Trasformazione delle coordinate ptP.ToLoc( ConusFrame) ; vtV.ToLoc( ConusFrame) ; DBLVECTOR vdCoef(3) ; DBLVECTOR vdRoots ; double dSqTan = dTan * dTan ; double dMinRad = dTan * dl ; double dMaxRad = dTan * dL ; double dDeltaR = dMaxRad - dMinRad ; double dHei = dL - dl ; vdCoef[0] = dSqTan * ptP.x * ptP.x - ptP.y * ptP.y - ptP.z * ptP.z ; vdCoef[1] = 2 * ( dSqTan * ptP.x * vtV.x - ptP.y * vtV.y - ptP.z * vtV.z) ; vdCoef[2] = dSqTan * vtV.x * vtV.x - vtV.y * vtV.y - vtV.z * vtV.z ; // Computo radici int nRoot = PolynomialRoots( 2, vdCoef, vdRoots) ; // Nessuna soluzione if ( nRoot == 0) return false ; double dEpsLow = ( bTapLow ? - EPS_SMALL : EPS_SMALL) ; double dEpsUp = ( bTapUp ? EPS_SMALL : - EPS_SMALL) ; // Una soluzione: la retta iterseca superficie // laterale e un piano if ( nRoot == 1) { ptInt1 = ptP + vdRoots[0] * vtV ; Vector3d vtU = ( ptInt1 - ORIG) - ( ptInt1 - ORIG).x * X_AX ; vtU.Normalize() ; vtN1 = dDeltaR * X_AX - dHei * vtU ; vtN1.Normalize() ; if ( ptInt1.x < dL + dEpsUp) { if ( ptInt1.x > dl + dEpsLow) { ptInt2 = ptP + ( ( dL - ptP.x) / vtV.x) * vtV ; vtN2 = - X_AX ; } else if ( ptInt1.x > - EPS_SMALL) { ptInt1 = ptP + ( ( dl - ptP.x) / vtV.x) * vtV ; ptInt2 = ptP + ( ( dL - ptP.x) / vtV.x) * vtV ; vtN1 = X_AX ; vtN2 = - X_AX ; if ( ptInt2.y * ptInt2.y + ptInt2.z * ptInt2.z > dMaxRad * dMaxRad) return false ; } else return false ; ptInt1.ToGlob( ConusFrame) ; ptInt2.ToGlob( ConusFrame) ; vtN1.ToGlob( ConusFrame) ; vtN2.ToGlob( ConusFrame) ; return true ; } else return false ; } // Due soluzioni: la retta interseca due volte la // superficie laterale else if ( nRoot == 2) { ptInt1 = ptP + vdRoots[0] * vtV ; ptInt2 = ptP + vdRoots[1] * vtV ; if ( ptInt1.x > ptInt2.x) swap( ptInt1, ptInt2) ; Vector3d vtU1 = ( ptInt1 - ORIG) - ( ptInt1 - ORIG).x * X_AX ; Vector3d vtU2 = ( ptInt2 - ORIG) - ( ptInt2 - ORIG).x * X_AX ; vtU1.Normalize() ; vtU2.Normalize() ; vtN1 = dDeltaR * X_AX - dHei * vtU1 ; vtN2 = dDeltaR * X_AX - dHei * vtU2 ; vtN1.Normalize() ; vtN2.Normalize() ; if ( abs( vtV.x) < EPS_ZERO) { if ( ptInt1.x > dl + dEpsLow && ptInt1.x < dL + dEpsUp) { ptInt1.ToGlob( ConusFrame) ; ptInt2.ToGlob( ConusFrame) ; vtN1.ToGlob( ConusFrame) ; vtN2.ToGlob( ConusFrame) ; vtN1.Normalize() ; vtN2.Normalize() ; return true ; } else return false ; } if ( ptInt1.x < dL + dEpsUp) { if ( ptInt1.x > dl + dEpsLow) { if ( ptInt2.x > dL + dEpsUp) { ptInt2 = ptP + ( ( dL - ptP.x) / vtV.x) * vtV ; vtN2 = - X_AX ; } } else if ( ptInt1.x > - EPS_SMALL) { if ( ptInt2.x > dL + dEpsUp) { ptInt1 = ptP + ( ( dl - ptP.x) / vtV.x) * vtV ; ptInt2 = ptP + ( ( dL - ptP.x) / vtV.x) * vtV ; vtN1 = X_AX ; vtN2 = - X_AX ; } else if ( ptInt2.x > dl + dEpsLow) { ptInt1 = ptP + ( ( dl - ptP.x) / vtV.x) * vtV ; vtN1 = X_AX ; } else return false ; } else { if ( ptInt2.x < 0) return false ; else if ( ptInt2.x < dl + dEpsLow) { ptInt1 = ptP + ( ( dl - ptP.x) / vtV.x) * vtV ; ptInt2 = ptP + ( ( dL - ptP.x) / vtV.x) * vtV ; vtN1 = X_AX ; vtN2 = - X_AX ; } else if ( ptInt2.x < dL + dEpsUp) { ptInt1 = ptP + ( ( dL - ptP.x) / vtV.x) * vtV ; vtN1 = - X_AX ; } else return false ; } ptInt1.ToGlob( ConusFrame) ; ptInt2.ToGlob( ConusFrame) ; vtN1.ToGlob( ConusFrame) ; vtN2.ToGlob( ConusFrame) ; vtN1.Normalize() ; vtN2.Normalize() ; return true ; } else return false ; } return false ; } //---------------------------------------------------------------------------- bool VolZmap::IntersLineEllipticalCylinder( const Vector3d& vtLineDir, const Point3d& ptLineSt, const Frame3d& CircFrame, double dSqRad, double dLongMvLen, double dOrtMvLen, Point3d& ptInt1, Point3d& ptInt2, Vector3d& vtN1, Vector3d& vtN2, bool bTapLow, bool bTapUp) { // NB: L'origine del sistema di riferimento deve essere // nel centro della circonferenza di base, la cui tralsazione obliqua // genera il cilindro ellittico, e l'asse x 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 x e y del sistema di riferimento CircFrame. double dObCoef = dOrtMvLen / dLongMvLen ; double dSqCoef = dObCoef * dObCoef ; Point3d ptP = ptLineSt ; Vector3d vtV = vtLineDir ; // Asse cilindro ellittico Vector3d vtAx( dLongMvLen, dOrtMvLen, 0) ; vtAx.Normalize() ; // Trasformazione delle coordinate ptP.ToLoc( CircFrame) ; vtV.ToLoc( CircFrame) ; vector vdCoef(3) ; vector vdRoots ; vdCoef[0] = dSqCoef * ptP.x * ptP.x + ptP.y * ptP.y + ptP.z * ptP.z - 2 * dObCoef * ptP.x * ptP.y - dSqRad ; vdCoef[1] = 2 * ( dSqCoef * vtV.x * ptP.x + vtV.y * ptP.y + vtV.z * ptP.z - dObCoef * ( vtV.x * ptP.y + vtV.y * ptP.x)) ; vdCoef[2] = dSqCoef * vtV.x * vtV.x + vtV.y * vtV.y + vtV.z * vtV.z - 2 * dObCoef * vtV.x * vtV.y ; int nRoot = PolynomialRoots( 2, vdCoef, vdRoots) ; // Nessuna soluzione if ( nRoot == 0 || nRoot == 1) { if ( abs( vtV.x) > EPS_ZERO) { ptInt1 = ptP - ( ptP.x / vtV.x) * vtV ; ptInt2 = ptP + ( ( dLongMvLen - ptP.x) / vtV.x) * vtV ; if ( ptInt1.y * ptInt1.y + ptInt1.z * ptInt1.z < dSqRad && ( ptInt2.y - dOrtMvLen) * ( ptInt2.y - dOrtMvLen) + ptInt2.z * ptInt2.z < dSqRad) { ptInt1.ToGlob( CircFrame) ; ptInt2.ToGlob( CircFrame) ; vtN1 = X_AX ; vtN2 = - X_AX ; vtN1.ToGlob( CircFrame) ; vtN2.ToGlob( CircFrame) ; return true ; } // Nessuna intersezione else return false ; } // Nessuna intersezione else return false ; } double dEpsLow = ( bTapLow ? - EPS_SMALL : EPS_SMALL) ; double dEpsUp = ( bTapUp ? EPS_SMALL : - EPS_SMALL) ; // L'equazione ammette o due soluzioni (eventualmente // coincidenti) oppure nessuna o infinite se la la retta // appartiene alla superficie Vector3d vtMv( dLongMvLen, dOrtMvLen, 0) ; if ( nRoot == 2) { ptInt1 = ptP + vdRoots[0] * vtV ; ptInt2 = ptP + vdRoots[1] * vtV ; if ( ptInt1.x > ptInt2.x) swap( ptInt1, ptInt2) ; Vector3d vtTest1 = ( ptInt1 - ORIG) - ( ptInt1 - ORIG) * vtAx * vtAx ; Vector3d vtTest2 = ( ptInt2 - ORIG) - ( ptInt2 - ORIG) * vtAx * vtAx ; double dY0_1, dY0_2 ; if ( vtTest1.y > 0) { dY0_1 = ( dSqRad - ptInt1.z * ptInt1.z > 0 ? sqrt( dSqRad - ptInt1.z * ptInt1.z) : 0) ; } else { dY0_1 = ( dSqRad - ptInt1.z * ptInt1.z > 0 ? - sqrt( dSqRad - ptInt1.z * ptInt1.z) : 0) ; } Vector3d vtCirc1( 0, - dY0_1, - ptInt1.z) ; Vector3d vtTan1( 0, - vtCirc1.z, vtCirc1.y) ; Vector3d vtCross1 = vtTan1 ^ vtMv ; vtN1 = ( vtCross1 * vtCirc1 > - EPS_ZERO ? vtCross1 : - vtCross1) ; if ( vtTest2.y > 0) { dY0_2 = ( dSqRad - ptInt2.z * ptInt2.z > 0 ? sqrt( dSqRad - ptInt2.z * ptInt2.z) : 0) ; } else { dY0_2 = ( dSqRad - ptInt2.z * ptInt2.z > 0 ? - sqrt( dSqRad - ptInt2.z * ptInt2.z) : 0) ; } Vector3d vtCirc2( 0, - dY0_2, - ptInt2.z) ; Vector3d vtTan2( 0, - vtCirc2.z, vtCirc2.y) ; Vector3d vtCross2 = vtTan2 ^ vtMv ; vtN2 = ( vtCross2 * vtCirc2 > - EPS_ZERO ? vtCross2 : - vtCross2) ; if ( ptInt1.x < dLongMvLen + dEpsUp) { if ( ptInt1.x > + dEpsLow) { if ( ptInt2.x > dLongMvLen + dEpsUp) { ptInt2 = ptP + ( ( dLongMvLen - ptP.x) / vtV.x) * vtV ; vtN2 = - X_AX ; } } else { if ( ptInt2.x > dLongMvLen + dEpsUp) { ptInt1 = ptP - ( ptP.x / vtV.x) * vtV ; ptInt2 = ptP + ( ( dLongMvLen - ptP.x) / vtV.x) * vtV ; vtN1.Set( 1, 0, 0) ; vtN2.Set( -1, 0, 0) ; if ( ptInt1.y * ptInt1.y + ptInt1.z * ptInt1.z > dSqRad && ptInt2.y * ptInt2.y + ptInt2.z * ptInt2.z > dSqRad) return false ; } else if ( ptInt2.x > dEpsLow) { ptInt1 = ptP - ( ptP.x / vtV.x) * vtV ; vtN1.Set( 1, 0, 0) ; } else return false ; } } else return false ; // Riporto le coordinate nel sistema di riferimento // griglia ptInt1.ToGlob( CircFrame) ; ptInt2.ToGlob( CircFrame) ; vtN1.ToGlob( CircFrame) ; vtN2.ToGlob( CircFrame) ; vtN1.Normalize() ; vtN2.Normalize() ; } return true ; } //---------------------------------------------------------------------------- bool VolZmap::IntersLineMyPolyhedron( const Point3d& ptLineSt, const Vector3d& vtLineDir, const Frame3d& PolyFrame, double dLenX, double dLenY, double dLenZ, double dDeltaX, Point3d& ptInt1, Point3d& ptInt2, Vector3d& vtN1, Vector3d& vtN2) { Point3d ptP = ptLineSt ; Vector3d vtV = vtLineDir ; // Trasformazione delle coordinate ptP.ToLoc( PolyFrame) ; 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, 0, dLenZ /2) ; Point3d ptFacet246( dLenX + dDeltaX, dLenY, - dLenZ / 2) ; // Servono per descrivere i piani obliqui Vector3d vtFacet5 = ptFacet135 - ptP ; Vector3d vtFacet6 = ptFacet246 - ptP ; Vector3d vtOb( dLenY, - dDeltaX, 0) ; vtOb.Normalize() ; Point3d ptI1 = ptP + ( ( ptFacet135.z - ptP.z) / vtV.z) * vtV ; Point3d ptI2 = ptP + ( ( ptFacet246.z - ptP.z) / vtV.z) * vtV ; Point3d ptI3 = ptP + ( ( ptFacet135.y - ptP.y) / vtV.y) * vtV ; Point3d ptI4 = ptP + ( ( ptFacet246.y - ptP.y) / vtV.y) * vtV ; Point3d ptI5 = ptP + ( ( vtFacet5 * vtOb) / ( vtV * vtOb)) * vtV ; Point3d ptI6 = ptP + ( ( vtFacet6 * vtOb) / ( vtV * vtOb)) * vtV ; // Controlli affinché non vengano tagliati dexel a filo // con il passaggio dell'utensile: // Controllo sulle facce 1 e 2 if ( abs( vtV.z) < EPS_ZERO && abs( ptP.z) > dLenZ / 2 - EPS_SMALL) return false ; // Controllo sulle facce 3 e 4 if ( abs( vtV.y) < EPS_ZERO && ( ptP.y < EPS_SMALL || ptP.y > dLenY - EPS_SMALL)) return false ; // Controllo sulle facce 5 e 6 Vector3d vtW( dDeltaX, dLenY, 0) ; vtW.Normalize() ; Vector3d vtU = vtV - vtV.z * Z_AX - vtV * vtW * vtW ; if ( vtU.Len() < EPS_ZERO && ( ptP.x * dLenY < dDeltaX * ptP.y + dLenY * EPS_SMALL || ptP.x * dLenY > dDeltaX * ptP.y + dLenY * ( dLenX - EPS_SMALL))) return false ; // Ricerca intersezioni con le facce int nIntNum = 0 ; // Intersezione con la prima faccia if ( ptI1.y >= 0 && ptI1.y <= dLenY && ptI1.x * dLenY >= dDeltaX * ptI1.y && ( ptI1.x - dLenX) * dLenY <= dDeltaX * ptI1.y) { ptInt1 = ptI1 ; vtN1 = - Z_AX ; ++ nIntNum ; } // Intersezione con la seconda faccia if ( ptI2.y >= 0 && ptI2.y <= dLenY && ptI2.x * dLenY >= dDeltaX * ptI2.y && ( ptI2.x - dLenX) * dLenY <= dDeltaX * ptI2.y) { if ( nIntNum == 0) { ptInt1 = ptI2 ; vtN1 = Z_AX ; ++ nIntNum ; } else if ( ( ptInt1 - ptI2).SqLen() > SQ_EPS_SMALL) { ptInt2 = ptI2 ; vtN2 = Z_AX ; ++ nIntNum ; } } // Intersezione con la terza faccia if ( nIntNum < 2 && ptI3.x >= 0 && ptI3.x <= dLenX && ptI3.z >= - ptFacet135.z && ptI3.z <= ptFacet135.z) { if ( nIntNum == 0) { ptInt1 = ptI3 ; vtN1 = Y_AX ; ++ nIntNum ; } else if ( ( ptInt1 - ptI3).SqLen() > SQ_EPS_SMALL) { ptInt2 = ptI3 ; vtN2 = Y_AX ; ++ nIntNum ; } } // Intersezione con la quarta faccia if ( nIntNum < 2 && ptI4.x >= dDeltaX && ptI4.x <= dLenX + dDeltaX && ptI4.z >= - ptFacet135.z && ptI4.z <= ptFacet135.z) { if ( nIntNum == 0) { ptInt1 = ptI4 ; vtN1 = - Y_AX ; ++ nIntNum ; } else if ( ( ptInt1 - ptI4).SqLen() > SQ_EPS_SMALL) { ptInt2 = ptI4 ; vtN2 = - Y_AX ; ++ nIntNum ; } } // Intersezione con la quinta faccia if ( nIntNum < 2 && ptI5.y >= 0 && ptI5.y <= dLenY && ptI5.z >= - ptFacet135.z && ptI5.z <= ptFacet135.z) { 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.y >= 0 && ptI6.y <= dLenY && ptI6.z >= - ptFacet135.z && ptI6.z <= ptFacet135.z) { if ( nIntNum == 0) { ptInt1 = ptI6; vtN1 = - vtOb ; ++ nIntNum ; } else if ( ( ptInt1 - ptI6).SqLen() > SQ_EPS_SMALL) { ptInt2 = ptI6; vtN2 = - vtOb ; ++ nIntNum ; } } if ( nIntNum == 2) { ptInt1.ToGlob( PolyFrame) ; ptInt2.ToGlob( PolyFrame) ; vtN1.ToGlob( PolyFrame) ; vtN2.ToGlob( PolyFrame) ; return true ; } else return false ; } //---------------------------------------------------------------------------- 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, POCRVVECTOR& 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) ; }