diff --git a/VolZmapCalculus.cpp b/VolZmapCalculus.cpp index bb66301..24762d9 100644 --- a/VolZmapCalculus.cpp +++ b/VolZmapCalculus.cpp @@ -248,15 +248,11 @@ VolZmap::IntersRayDexel( const Point3d& ptP, const Vector3d& vtV, int nGrid, int 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) ; - else - return GetDepthWithDexel( ptP, vtD, dInLength, dOutLength) ; + return GetDepthWithVoxel( ptPLoc, vtDLoc, dInLength, dOutLength) ; + else { + return GetDepthWithDexel( ptPLoc, vtDLoc, dInLength, dOutLength) ; + } } //---------------------------------------------------------------------------- @@ -265,38 +261,45 @@ VolZmap::GetDepth( const Point3d& ptPLoc, const Vector3d& vtDLoc, double& dInLen // 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 +VolZmap::GetDepthWithDexel( const Point3d& ptPLoc, const Vector3d& vtDLoc, double& dInLength, double& dOutLength) const { + // Punto e direzione nel riferimento intrinseco + Point3d ptP = ptPLoc ; + Vector3d vtD = vtDLoc ; + ptP.ToLoc( m_MapFrame) ; + vtD.ToLoc( m_MapFrame) ; + vtD.Normalize() ; + // Intersezione fra semiretta e BBox dello Zmap double dU1, dU2 ; - bool bTest = IntersLineZMapBBox( ptP, vtV, dU1, dU2) && ( dU1 > 0 || dU2 > 0) ; + bool bTest = IntersLineZMapBBox( ptP, vtD, dU1, dU2) && ( dU1 > 0 || dU2 > 0) ; // Semiretta esterna al box dello Zmap quindi esterna anche allo Zmap if ( ! bTest) { - dInLength = - 2 ; - dOutLength = - 2 ; + dInLength = -2 ; + dOutLength = -2 ; return true ; } - double dInLen[3] ; - double dOutLen[3] ; + double dInLen[N_MAPS] ; + double dOutLen[N_MAPS] ; // Ciclo sulle griglie - for ( int nGrid = 0 ; nGrid < int( m_nMapNum) ; ++ nGrid) { + for ( int nGrid = 0 ; nGrid < m_nMapNum ; ++ nGrid) { Point3d ptP0 = ptP ; - Vector3d vtV0 = vtV ; + Vector3d vtV0 = vtD ; Point3d ptI, ptF ; // Una sola intersezione valida ( punto interno, intersezione valida 2) if ( dU1 < 0 && dU2 > 0) { ptI = ptP ; - ptF = ptP + dU2 * vtV ; + ptF = ptP + dU2 * vtD ; } // due soluzioni valide ( punto esterno) else { - ptI = ptP + dU1 * vtV ; - ptF = ptP + dU2 * vtV ; + ptI = ptP + dU1 * vtD ; + ptF = ptP + dU2 * vtD ; } // Passo dal sistema intrinseco alla griglia if ( nGrid == 1) { @@ -382,12 +385,12 @@ VolZmap::GetDepthWithDexel( const Point3d& ptP, const Vector3d& vtV, double& dIn // Se non abbiamo incontrato materiale if ( dInLen[nGrid] > dOutLen[nGrid] - EPS_SMALL) { - dInLen[nGrid] = - 2 ; - dOutLen[nGrid] = - 2 ; + dInLen[nGrid] = -2 ; + dOutLen[nGrid] = -2 ; } // Se parto dall'interno else if ( dInLen[nGrid] < - EPS_SMALL) - dInLen[nGrid] = - 1 ; + dInLen[nGrid] = -1 ; } if ( m_nMapNum == 1) { @@ -399,17 +402,17 @@ VolZmap::GetDepthWithDexel( const Point3d& ptP, const Vector3d& vtV, double& dIn if ( abs( dInLen[0] + 2) < EPS_SMALL && abs( dInLen[1] + 2) < EPS_SMALL && abs( dInLen[2] + 2) < EPS_SMALL) { - dInLength = - 2 ; - dOutLength = - 2 ; + 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) + if ( dInLen[nGr] > -2 && dInLen[nGr] < dInLength) dInLength = dInLen[nGr] ; - if ( dOutLen[nGr] > - 2 && dOutLen[nGr] > dOutLength) + if ( dOutLen[nGr] > -2 && dOutLen[nGr] > dOutLength) dOutLength = dOutLen[nGr] ; } } @@ -430,14 +433,10 @@ VolZmap::GetDepthWithVoxel( const Point3d& ptP, const Vector3d& vtD, double& dIn return false ; if ( vIntersInfo.empty()) { - dInLength = - 2. ; - dOutLength = - 2. ; + dInLength = -2. ; + dOutLength = -2. ; return true ; } - - // Porto nel riferimento intrinseco dello Zmap - Vector3d vtOutD = vtD ; - vtOutD.ToGlob( m_MapFrame) ; // Inizializzo le distanze di ingresso e uscita: // dInLength diminuisce, dOutLength aumenta. @@ -445,46 +444,48 @@ VolZmap::GetDepthWithVoxel( const Point3d& ptP, const Vector3d& vtD, double& dIn dOutLength = -INFINITO ; int nFirstPosN ; - int nN = 0 ; - for ( ; nN < int( vIntersInfo.size()) ; ++ nN) { + for ( int nN = 0; nN < int( vIntersInfo.size()) ; ++ nN) { if ( vIntersInfo[nN].dU > - EPS_SMALL) { nFirstPosN = nN ; break ; } } - if ( nN == int( vIntersInfo.size())) { - dInLength = - 2 ; - dOutLength = - 2 ; + if ( nFirstPosN == int( vIntersInfo.size())) { + dInLength = -2 ; + dOutLength = -2 ; return true ; } - int nVoxIJK[3] ; - GetBlockIJKFromN( vIntersInfo[0].nVox, nVoxIJK) ; - if ( nFirstPosN > 0) { - if ( vIntersInfo[nFirstPosN - 1].trTria.GetN() * vtOutD < EPS_ZERO) - dInLength = - 1 ; + if ( vIntersInfo[nFirstPosN - 1].trTria.GetN() * vtD < EPS_ZERO) + dInLength = -1 ; } else if ( nFirstPosN == 0) { - if ( vIntersInfo[nFirstPosN].trTria.GetN() * vtOutD > EPS_ZERO) - dInLength = - 1 ; - else if ( GetPointVoxel( ptP, nVoxIJK[0], nVoxIJK[1], nVoxIJK[2])) { - int nCubeType = CalcIndex( nVoxIJK[0], nVoxIJK[1], nVoxIJK[2]) ; - if ( nCubeType == 255) - dInLength = - 1 ; + if ( vIntersInfo[nFirstPosN].trTria.GetN() * vtD > EPS_ZERO) + dInLength = -1 ; + else { + Point3d ptPi = ptP ; + ptPi.ToLoc( m_MapFrame) ; + int nVoxIJK[3] ; + GetBlockIJKFromN( vIntersInfo[0].nVox, nVoxIJK) ; + if ( GetPointVoxel( ptPi, nVoxIJK[0], nVoxIJK[1], nVoxIJK[2])) { + int nCubeType = CalcIndex( nVoxIJK[0], nVoxIJK[1], nVoxIJK[2]) ; + if ( nCubeType == 255) + dInLength = -1 ; + } } } for ( int nN = nFirstPosN ; nN < int( vIntersInfo.size()) ; ++ nN) { - if ( vIntersInfo[nN].trTria.GetN() * vtOutD > - EPS_ZERO) { + if ( vIntersInfo[nN].trTria.GetN() * vtD > - EPS_ZERO) { if ( ( vIntersInfo[nN].nILTT == ILTT_SEGM || vIntersInfo[nN].nILTT == ILTT_SEGM_ON_EDGE) && dOutLength < vIntersInfo[nN].dU2) dOutLength = vIntersInfo[nN].dU2 ; else if ( dOutLength < vIntersInfo[nN].dU) dOutLength = vIntersInfo[nN].dU ; } - if ( vIntersInfo[nN].trTria.GetN() * vtOutD < EPS_ZERO && + if ( vIntersInfo[nN].trTria.GetN() * vtD < EPS_ZERO && dInLength > vIntersInfo[nN].dU) dInLength = vIntersInfo[nN].dU ; } @@ -521,10 +522,10 @@ VolZmap::AvoidBox( const Frame3d& frBox, const Vector3d& vtDiag, double dSafeDis 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) ; + 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, frB) ; @@ -593,10 +594,10 @@ VolZmap::AvoidSphere( const Point3d& ptCenter, double dRad, double dSafeDist) co 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) ; + 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) { @@ -681,10 +682,10 @@ VolZmap::AvoidCylinder( const Frame3d& frCyl, double dH, double dR, double dSafe 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) ; + 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 il cilindro (nel riferimento intrinseco) for ( int i = nStI ; i <= nEnI ; ++ i) { @@ -1331,15 +1332,19 @@ VolZmap::GetPartVolume( int nPart, double& dVol) const bool VolZmap::GetLineIntersection( const Point3d& ptP, const Vector3d& vtD, ILZIVECTOR& vIntersInfo) const { - // Calcolo punto e vettore espressi nel riferimento intrinseco dello Zmap + // Direzione normalizzata + Vector3d vtDir = vtD ; + vtDir.Normalize() ; + + // Punto e vettore espressi nel riferimento intrinseco dello Zmap Point3d ptLocP = ptP ; - Vector3d vtLocD = vtD ; + Vector3d vtDirL = vtDir ; ptLocP.ToLoc( m_MapFrame) ; - vtLocD.ToLoc( m_MapFrame) ; + vtDirL.ToLoc( m_MapFrame) ; // Intersezione fra semiretta e BBox dello Zmap double dU1, dU2 ; - bool bLineBBoxInters = IntersLineZMapBBox( ptLocP, vtLocD, dU1, dU2) && ( dU1 > 0 || dU2 > 0) ; + bool bLineBBoxInters = IntersLineZMapBBox( ptLocP, vtDirL, dU1, dU2) && ( dU1 > 0 || dU2 > 0) ; // Semiretta esterna al box dello Zmap quindi esterna anche allo Zmap if ( ! bLineBBoxInters) @@ -1362,7 +1367,7 @@ VolZmap::GetLineIntersection( const Point3d& ptP, const Vector3d& vtD, ILZIVECTO if ( ! GetBlockBox( nBlockIJK, b3BlockBox)) return false ; // Se c'è intersezione valuto tutti i voxel interni - if ( IntersLineBox( ptLocP, vtLocD, b3BlockBox.GetMin(), b3BlockBox.GetMax())) { + if ( IntersLineBox( ptLocP, vtDirL, b3BlockBox.GetMin(), b3BlockBox.GetMax())) { // Ciclo sui voxel del blocco. // Triangoli smooth for ( int nV = 0 ; nV < int( m_BlockSmoothTria[nB].size()) ; ++ nV) { @@ -1373,7 +1378,7 @@ VolZmap::GetLineIntersection( const Point3d& ptP, const Vector3d& vtD, ILZIVECTO m_BlockSmoothTria[nB][nV].k } ; GetVoxelBox( nCurVoxIJK[0], nCurVoxIJK[1], nCurVoxIJK[2], b3Vox) ; // Se non c'è intersezione col voxel, passo al successivo. - if ( ! IntersLineBox( ptLocP, vtLocD, b3Vox.GetMin(), b3Vox.GetMax())) + if ( ! IntersLineBox( ptLocP, vtDirL, b3Vox.GetMin(), b3Vox.GetMax())) continue ; for ( int nT = 0 ; nT < int( m_BlockSmoothTria[nB][nV].vTria.size()) ; ++ nT) { Triangle3d trTria = m_BlockSmoothTria[nB][nV].vTria[nT] ; @@ -1381,7 +1386,7 @@ VolZmap::GetLineIntersection( const Point3d& ptP, const Vector3d& vtD, ILZIVECTO continue ; Point3d ptLineTria1, ptLineTria2 ; // Studio dell'intersezione della retta con il triangolo corrente - int nIntType = IntersLineTria( ptP, vtD, 1.5 * dU2, trTria, ptLineTria1, ptLineTria2) ; + int nIntType = IntersLineTria( ptP, vtDir, 1.5 * dU2, trTria, ptLineTria1, ptLineTria2) ; // Se non ci sono intersezioni passo al prossimo triangolo if ( nIntType == ILTT_NO) continue ; @@ -1389,15 +1394,15 @@ VolZmap::GetLineIntersection( const Point3d& ptP, const Vector3d& vtD, ILZIVECTO else if ( nIntType == ILTT_VERT || nIntType == ILTT_EDGE || nIntType == ILTT_IN) { int nNumVox ; GetVoxNFromIJK( nCurVoxIJK[0], nCurVoxIJK[1], nCurVoxIJK[2], nNumVox) ; - vIntersInfo.emplace_back( nIntType, ( ptLineTria1 - ptP) * vtD, + vIntersInfo.emplace_back( nIntType, ( ptLineTria1 - ptP) * vtDir, nNumVox, nB, ptLineTria1, trTria) ; } // altrimenti ci sono due intersezioni else { int nNumVox ; GetVoxNFromIJK( nCurVoxIJK[0], nCurVoxIJK[1], nCurVoxIJK[2], nNumVox) ; - double dP1 = ( ptLineTria1 - ptP) * vtD ; - double dP2 = ( ptLineTria2 - ptP) * vtD ; + double dP1 = ( ptLineTria1 - ptP) * vtDir ; + double dP2 = ( ptLineTria2 - ptP) * vtDir ; vIntersInfo.emplace_back( nIntType, ( dP1 < dP2 ? dP1 : dP2), ( dP1 < dP2 ? dP2 : dP1), nNumVox, nB, ptLineTria1, ptLineTria2, trTria) ; } @@ -1417,7 +1422,7 @@ VolZmap::GetLineIntersection( const Point3d& ptP, const Vector3d& vtD, ILZIVECTO continue ; Point3d ptLineTria1, ptLineTria2 ; // Studio dell'intersezione della retta con il triangolo corrente - int nIntType = IntersLineTria( ptP, vtD, 1.5 * dU2, trTria, ptLineTria1, ptLineTria2) ; + int nIntType = IntersLineTria( ptP, vtDir, 1.5 * dU2, trTria, ptLineTria1, ptLineTria2) ; // Se non ci sono intersezioni passo al prossimo triangolo if ( nIntType == ILTT_NO) continue ; @@ -1425,15 +1430,15 @@ VolZmap::GetLineIntersection( const Point3d& ptP, const Vector3d& vtD, ILZIVECTO else if ( nIntType == ILTT_VERT || nIntType == ILTT_EDGE || nIntType == ILTT_IN) { int nNumVox ; GetVoxNFromIJK( nCurVoxIJK[0], nCurVoxIJK[1], nCurVoxIJK[2], nNumVox) ; - vIntersInfo.emplace_back( nIntType, ( ptLineTria1 - ptP) * vtD, + vIntersInfo.emplace_back( nIntType, ( ptLineTria1 - ptP) * vtDir, nNumVox, nB, ptLineTria1, trTria) ; } // altrimenti ci sono due intersezioni else { int nNumVox ; GetVoxNFromIJK( nCurVoxIJK[0], nCurVoxIJK[1], nCurVoxIJK[2], nNumVox) ; - double dP1 = ( ptLineTria1 - ptP) * vtD ; - double dP2 = ( ptLineTria2 - ptP) * vtD ; + double dP1 = ( ptLineTria1 - ptP) * vtDir ; + double dP2 = ( ptLineTria2 - ptP) * vtDir ; vIntersInfo.emplace_back( nIntType, ( dP1 < dP2 ? dP1 : dP2), ( dP1 < dP2 ? dP2 : dP1), nNumVox, nB, ptLineTria1, ptLineTria2, trTria) ; } @@ -1445,19 +1450,19 @@ VolZmap::GetLineIntersection( const Point3d& ptP, const Vector3d& vtD, ILZIVECTO Triangle3d trTria = m_BlockBigTria[nB][nT] ; Point3d ptLineTria1, ptLineTria2 ; // Studio dell'intersezione della retta con il triangolo corrente - int nIntType = IntersLineTria( ptP, vtD, 1.5 * dU2, trTria, ptLineTria1, ptLineTria2) ; + int nIntType = IntersLineTria( ptP, vtDir, 1.5 * dU2, trTria, 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) { - vIntersInfo.emplace_back( nIntType, ( ptLineTria1 - ptP) * vtD, + vIntersInfo.emplace_back( nIntType, ( ptLineTria1 - ptP) * vtDir, -1, nB, ptLineTria1, trTria) ; } // altrimenti ci sono due intersezioni else { - double dP1 = ( ptLineTria1 - ptP) * vtD ; - double dP2 = ( ptLineTria2 - ptP) * vtD ; + double dP1 = ( ptLineTria1 - ptP) * vtDir ; + double dP2 = ( ptLineTria2 - ptP) * vtDir ; vIntersInfo.emplace_back( nIntType, ( dP1 < dP2 ? dP1 : dP2), ( dP1 < dP2 ? dP2 : dP1), -1, nB, ptLineTria1, ptLineTria2, trTria) ; } @@ -1476,7 +1481,7 @@ VolZmap::GetLineIntersection( const Point3d& ptP, const Vector3d& vtD, ILZIVECTO continue ; Point3d ptLineTria1, ptLineTria2 ; // Studio dell'intersezione della retta con il triangolo corrente - int nIntType = IntersLineTria( ptP, vtD, 1.5 * dU2, trTria, ptLineTria1, ptLineTria2) ; + int nIntType = IntersLineTria( ptP, vtDir, 1.5 * dU2, trTria, ptLineTria1, ptLineTria2) ; // Se non ci sono intersezioni passo al prossimo triangolo if ( nIntType == ILTT_NO) continue ; @@ -1484,15 +1489,15 @@ VolZmap::GetLineIntersection( const Point3d& ptP, const Vector3d& vtD, ILZIVECTO else if ( nIntType == ILTT_VERT || nIntType == ILTT_EDGE || nIntType == ILTT_IN) { int nNumVox ; GetVoxNFromIJK( nCurVoxIJK[0], nCurVoxIJK[1], nCurVoxIJK[2], nNumVox) ; - vIntersInfo.emplace_back( nIntType, ( ptLineTria1 - ptP) * vtD, + vIntersInfo.emplace_back( nIntType, ( ptLineTria1 - ptP) * vtDir, nNumVox, nB, ptLineTria1, trTria) ; } // altrimenti ci sono due intersezioni else { int nNumVox ; GetVoxNFromIJK( nCurVoxIJK[0], nCurVoxIJK[1], nCurVoxIJK[2], nNumVox) ; - double dP1 = ( ptLineTria1 - ptP) * vtD ; - double dP2 = ( ptLineTria2 - ptP) * vtD ; + double dP1 = ( ptLineTria1 - ptP) * vtDir ; + double dP2 = ( ptLineTria2 - ptP) * vtDir ; vIntersInfo.emplace_back( nIntType, ( dP1 < dP2 ? dP1 : dP2), ( dP1 < dP2 ? dP2 : dP1), nNumVox, nB, ptLineTria1, ptLineTria2, trTria) ; }