EgtGeomKernel 2.1e3 :

- modifiche a Zmap per intersezioni con linee e piani.
This commit is contained in:
Dario Sassi
2019-05-13 06:36:58 +00:00
parent 4e8508f731
commit 90df6d3a7a
7 changed files with 758 additions and 768 deletions
+149 -55
View File
@@ -1971,7 +1971,7 @@ VolZmap::ExtMarchingCubes( int nBlock, VoxelContainer& vVox) const
Vector3d vtNlm = ( CompoVert[nComp][l].ptPApp - ptSol) ^ ( CompoVert[nComp][m].ptPApp - ptSol) ;
vtNlm.Normalize() ;
double dDotl = vtNlm * CompoVert[nComp][l].vtVec ;
double dDotm = vtNlm * CompoVert[nComp][m].vtVec;
double dDotm = vtNlm * CompoVert[nComp][m].vtVec ;
if ( dDotl < - 0.85 || dDotm < - 0.85) {
bOverTurning = true ;
break ;
@@ -2000,6 +2000,44 @@ VolZmap::ExtMarchingCubes( int nBlock, VoxelContainer& vVox) const
}
}
// Gestione casi speciali: tutti i vettori quasi equiversi tranne uno
if ( bExtConfirmed && ( nVertComp[nComp] == 4 || nVertComp[nComp] == 5)) {
int nParNum = 0 ;
for ( int m = 0 ; m < nVertComp[nComp] - 1 ; ++ m) {
for ( int l = m + 1 ; l < nVertComp[nComp] ; ++ l) {
if ( CompoVert[nComp][m].vtVec * CompoVert[nComp][l].vtVec > 0.9) {
++ nParNum ;
}
}
}
if ( nParNum == ( nVertComp[nComp] - 1) * ( nVertComp[nComp] - 2) / 2) {
int nNotParIndex = - 1 ;
for ( int m = 0 ; m < nVertComp[nComp] ; ++ m) {
int nNumDiff = 0 ;
for ( int l = 0 ; l < nVertComp[nComp] ; ++ l) {
if ( CompoVert[nComp][m].vtVec * CompoVert[nComp][l].vtVec <= 0.9) {
++ nNumDiff ;
}
}
if ( nNumDiff == nVertComp[nComp] - 1) {
nNotParIndex = m ;
break ;
}
}
int nPrevL = nNotParIndex > 0 ? nNotParIndex - 1 : nVertComp[nComp] - 1 ;
int nNextL = nNotParIndex < nVertComp[nComp] - 1 ? nNotParIndex + 1 : 0 ;
Vector3d vtPrevN = ( CompoVert[nComp][nNotParIndex].ptPApp - ptSol) ^ ( CompoVert[nComp][nPrevL].ptPApp - ptSol) ;
Vector3d vtNextN = ( CompoVert[nComp][nNextL].ptPApp - ptSol) ^ ( CompoVert[nComp][nNotParIndex].ptPApp - ptSol) ;
vtPrevN.Normalize() ;
vtNextN.Normalize() ;
if ( nNotParIndex >= 0 &&
CompoVert[nComp][nNotParIndex].vtVec * vtPrevN < - EPS_SMALL &&
CompoVert[nComp][nNotParIndex].vtVec * vtNextN < - EPS_SMALL) {
bExtConfirmed = false ;
}
}
}
// ExtMC confermato
if ( bExtConfirmed) {
int tOldCompo = VoxConf.nNumComp ;
@@ -3450,9 +3488,9 @@ VolZmap::IsValidVoxel( int nI, int nJ, int nK) const
int nVoxNumY = int( m_nNy[0] / N_DEXVOXRATIO + ( m_nNy[0] % N_DEXVOXRATIO == 0 ? 1 : 2)) ;
int nVoxNumZ = int( m_nNy[1] / N_DEXVOXRATIO + ( m_nNy[1] % N_DEXVOXRATIO == 0 ? 1 : 2)) ;
// Verifico la validità del voxel
return ( nI >= - 1 && nI < nVoxNumX - 1 &&
nJ >= - 1 && nJ < nVoxNumY - 1 &&
nK >= - 1 && nK < nVoxNumZ - 1 ) ;
return ( nI >= -1 && nI <= nVoxNumX - 2 &&
nJ >= -1 && nJ <= nVoxNumY - 2 &&
nK >= -1 && nK <= nVoxNumZ - 2) ;
}
//----------------------------------------------------------------------------
@@ -3474,9 +3512,9 @@ VolZmap::GetVoxIJKFromN( int nN, int& nI, int& nJ, int& nK) const
nK = ( nN / ( nVoxNumX * nVoxNumY)) - 1 ;
// Controllo la sensatezza del risultato ottenuto
return ( nI >= - 1 && nI < nVoxNumX - 1 &&
nJ >= - 1 && nJ < nVoxNumY - 1 &&
nK >= - 1 && nK < nVoxNumZ - 1 ) ;
return ( nI >= -1 && nI <= nVoxNumX - 2 &&
nJ >= -1 && nJ <= nVoxNumY - 2 &&
nK >= -1 && nK <= nVoxNumZ - 2) ;
}
//----------------------------------------------------------------------------
@@ -3489,9 +3527,9 @@ VolZmap::GetVoxNFromIJK( int nI, int nJ, int nK, int& nN) const
int nVoxNumZ = int( m_nNy[1] / N_DEXVOXRATIO + ( m_nNy[1] % N_DEXVOXRATIO == 0 ? 1 : 2)) ;
// Controllo la validità del voxel
if ( nI <= - 2 || nI >= nVoxNumX - 1 ||
nJ <= - 2 || nJ >= nVoxNumY - 1 ||
nK <= - 2 || nK >= nVoxNumZ - 1 )
if ( nI < -1 || nI > nVoxNumX - 2 ||
nJ < -1 || nJ > nVoxNumY - 2 ||
nK < -1 || nK > nVoxNumZ - 2)
return false ;
// Calcolo il numero di Voxel
@@ -3501,6 +3539,32 @@ VolZmap::GetVoxNFromIJK( int nI, int nJ, int nK, int& nN) const
return ( nN >= 0 && nN < nVoxNumX * nVoxNumY * nVoxNumZ) ;
}
//----------------------------------------------------------------------------
bool
VolZmap::GetVoxelBox( int nI, int nJ, int nK, BBox3d& b3VoxBox) const
{
// Calcolo il numero di voxel lungo X,Y e Z
int nVoxNumX = int( m_nNx[0] / N_DEXVOXRATIO + ( m_nNx[0] % N_DEXVOXRATIO == 0 ? 1 : 2)) ;
int nVoxNumY = int( m_nNy[0] / N_DEXVOXRATIO + ( m_nNy[0] % N_DEXVOXRATIO == 0 ? 1 : 2)) ;
int nVoxNumZ = int( m_nNy[1] / N_DEXVOXRATIO + ( m_nNy[1] % N_DEXVOXRATIO == 0 ? 1 : 2)) ;
// Controllo la validità del voxel
if ( nI < -1 || nI > nVoxNumX - 2 ||
nJ < -1 || nJ > nVoxNumY - 2 ||
nK < -1 || nK > nVoxNumZ - 2)
return false ;
// Punti della diagonale del voxel
Point3d ptCubeInf( ( nI * N_DEXVOXRATIO + 0.5) * m_dStep,
( nJ * N_DEXVOXRATIO + 0.5) * m_dStep,
( nK * N_DEXVOXRATIO + 0.5) * m_dStep) ;
Point3d ptCubeSup( ( ( nI + 1) * N_DEXVOXRATIO + 0.5) * m_dStep,
( ( nJ + 1) * N_DEXVOXRATIO + 0.5) * m_dStep,
( ( nK + 1) * N_DEXVOXRATIO + 0.5) * m_dStep) ;
b3VoxBox.Set( ptCubeInf, ptCubeSup) ;
return true ;
}
//----------------------------------------------------------------------------
bool
VolZmap::GetBlockIJKFromN( int nBlock, int nIJK[]) const
@@ -3538,9 +3602,9 @@ bool
VolZmap::GetBlockLimitsIJK( const int nIJK[], int nLimits[]) const
{
// Controllo sulla validità degli indici i, j, k del blocco
if ( nIJK[0] < 0 || nIJK[0] >= int( m_nFracLin[0]) ||
nIJK[1] < 0 || nIJK[1] >= int( m_nFracLin[1]) ||
nIJK[2] < 0 || nIJK[2] >= int( m_nFracLin[2]))
if ( nIJK[0] < 0 || nIJK[0] >= m_nFracLin[0] ||
nIJK[1] < 0 || nIJK[1] >= m_nFracLin[1] ||
nIJK[2] < 0 || nIJK[2] >= m_nFracLin[2])
return false ;
// Calcolo il numero di voxel lungo X,Y e Z
@@ -3549,23 +3613,42 @@ VolZmap::GetBlockLimitsIJK( const int nIJK[], int nLimits[]) const
int nVoxNumZ = int( m_nNy[1] / N_DEXVOXRATIO + ( m_nNy[1] % N_DEXVOXRATIO == 0 ? 1 : 2)) ;
// Calcolo limiti per l'indice i
nLimits[0] = ( nIJK[0] == 0 ? - 1 : nIJK[0] * int( m_nVoxNumPerBlock)) ;
nLimits[1] = ( nIJK[0] + 1 == int( m_nFracLin[0]) ?
nVoxNumX - 1 : ( nIJK[0] + 1) * int( m_nVoxNumPerBlock)) ;
nLimits[0] = ( nIJK[0] == 0 ? - 1 : nIJK[0] * m_nVoxNumPerBlock) ;
nLimits[1] = ( nIJK[0] + 1 == m_nFracLin[0] ?
nVoxNumX - 1 : ( nIJK[0] + 1) * m_nVoxNumPerBlock) ;
// Calcolo limiti per l'indice j
nLimits[2] = ( nIJK[1] == 0 ? - 1 : nIJK[1] * int( m_nVoxNumPerBlock)) ;
nLimits[3] = ( nIJK[1] + 1 == int( m_nFracLin[1]) ?
nVoxNumY - 1 : ( nIJK[1] + 1) * int( m_nVoxNumPerBlock)) ;
nLimits[2] = ( nIJK[1] == 0 ? - 1 : nIJK[1] * m_nVoxNumPerBlock) ;
nLimits[3] = ( nIJK[1] + 1 == m_nFracLin[1] ?
nVoxNumY - 1 : ( nIJK[1] + 1) * m_nVoxNumPerBlock) ;
// Calcolo limiti per l'indice k
nLimits[4] = ( nIJK[2] == 0 ? - 1 : nIJK[2] * int( m_nVoxNumPerBlock)) ;
nLimits[5] = ( nIJK[2] + 1 == int( m_nFracLin[2]) ?
nVoxNumZ - 1 : ( nIJK[2] + 1) * int( m_nVoxNumPerBlock)) ;
nLimits[4] = ( nIJK[2] == 0 ? - 1 : nIJK[2] * m_nVoxNumPerBlock) ;
nLimits[5] = ( nIJK[2] + 1 == m_nFracLin[2] ?
nVoxNumZ - 1 : ( nIJK[2] + 1) * m_nVoxNumPerBlock) ;
return true ;
}
//----------------------------------------------------------------------------
bool
VolZmap::GetBlockBox( const int nIJK[], BBox3d& b3VoxBox) const
{
// Calcolo limiti sugli indici del blocco; se il blocco non è valido, errore
int nLimits[6] ;
if ( ! GetBlockLimitsIJK( nIJK, nLimits))
return false ;
// Costruisco il bounding-bx del blocco
Point3d ptMinDiag( ( nLimits[0] * N_DEXVOXRATIO + 0.5) * m_dStep,
( nLimits[2] * N_DEXVOXRATIO + 0.5) * m_dStep,
( nLimits[4] * N_DEXVOXRATIO + 0.5) * m_dStep) ;
Point3d ptMaxDiag( ( nLimits[1] * N_DEXVOXRATIO + 0.5) * m_dStep,
( nLimits[3] * N_DEXVOXRATIO + 0.5) * m_dStep,
( nLimits[5] * N_DEXVOXRATIO + 0.5) * m_dStep) ;
b3VoxBox.Set( ptMinDiag, ptMaxDiag) ;
return true ;
}
//----------------------------------------------------------------------------
bool
VolZmap::IsPointInsideVoxelApprox( int nI, int nJ, int nK, const Point3d& ptP, double dPrec) const
@@ -3575,9 +3658,9 @@ VolZmap::IsPointInsideVoxelApprox( int nI, int nJ, int nK, const Point3d& ptP, d
int nVoxNumY = int( m_nNy[0] / N_DEXVOXRATIO + ( m_nNy[0] % N_DEXVOXRATIO == 0 ? 1 : 2)) ;
int nVoxNumZ = int( m_nNy[1] / N_DEXVOXRATIO + ( m_nNy[1] % N_DEXVOXRATIO == 0 ? 1 : 2)) ;
// Controllo sulla validità del voxel
if ( nI <= - 2 || nI >= nVoxNumX - 1 ||
nJ <= - 2 || nJ >= nVoxNumY - 1 ||
nK <= - 2 || nK >= nVoxNumZ - 1 )
if ( nI < -1 || nI > nVoxNumX - 2 ||
nJ < -1 || nJ > nVoxNumY - 2 ||
nK < -1 || nK > nVoxNumZ - 2)
return false ;
// Se la tolleranza è minore di EPS_SMALL, la considero nulla
if ( dPrec < EPS_ZERO)
@@ -3606,9 +3689,9 @@ VolZmap::GetPointVoxel( const Point3d& ptP, int& nVoxI, int& nVoxJ, int& nVoxK)
nVoxJ = int( floor( ( ptP.y - 0.5 * m_dStep) / ( m_dStep * N_DEXVOXRATIO))) ;
nVoxK = int( floor( ( ptP.z - 0.5 * m_dStep) / ( m_dStep * N_DEXVOXRATIO))) ;
// Controllo la validità del voxel
return ( nVoxI >= - 1 && nVoxI < nVoxNumX - 1) &&
( nVoxJ >= - 1 && nVoxJ < nVoxNumY - 1) &&
( nVoxK >= - 1 && nVoxK < nVoxNumZ - 1) ;
return ( nVoxI >= - 1 && nVoxI <= nVoxNumX - 2 &&
nVoxJ >= - 1 && nVoxJ <= nVoxNumY - 2 &&
nVoxK >= - 1 && nVoxK <= nVoxNumZ - 2) ;
}
//----------------------------------------------------------------------------
@@ -3620,9 +3703,9 @@ VolZmap::GetVoxelBlockIJK( const int nVoxIJK[], int nBlockIJK[]) const
int nVoxNumY = int( m_nNy[0] / N_DEXVOXRATIO + ( m_nNy[0] % N_DEXVOXRATIO == 0 ? 1 : 2)) ;
int nVoxNumZ = int( m_nNy[1] / N_DEXVOXRATIO + ( m_nNy[1] % N_DEXVOXRATIO == 0 ? 1 : 2)) ;
// Controllo sull'ammissibilità del voxel
if ( nVoxIJK[0] <= - 2 || nVoxIJK[0] >= nVoxNumX - 1 ||
nVoxIJK[1] <= - 2 || nVoxIJK[1] >= nVoxNumY - 1 ||
nVoxIJK[2] <= - 2 || nVoxIJK[2] >= nVoxNumZ - 1 )
if ( nVoxIJK[0] < -1 || nVoxIJK[0] > nVoxNumX - 2 ||
nVoxIJK[1] < -1 || nVoxIJK[1] > nVoxNumY - 2 ||
nVoxIJK[2] < -1 || nVoxIJK[2] > nVoxNumZ - 2)
return false ;
// Divisioni intere
int nIntRatio0 = nVoxIJK[0] / m_nVoxNumPerBlock ;
@@ -3636,14 +3719,25 @@ VolZmap::GetVoxelBlockIJK( const int nVoxIJK[], int nBlockIJK[]) const
return true ;
}
//----------------------------------------------------------------------------
bool
VolZmap::GetVoxelBlockN( int nVox, int& nBlock) const
{
int nVoxIJK[3] ;
int nBlockIJK[3] ;
return ( GetVoxIJKFromN( nVox, nVoxIJK[0], nVoxIJK[1], nVoxIJK[2]) &&
GetVoxelBlockIJK( nVoxIJK, nBlockIJK) &&
GetBlockNFromIJK( nBlockIJK, nBlock)) ;
}
//----------------------------------------------------------------------------
bool
VolZmap::GetAdjBlockToBlock( int nBlockN, int nDeltaI, int nDeltaJ, int nDeltaK, int& nAdjBlockN) const
{
// Test sulla validità degli incrementi su i,j,k
if ( nDeltaI < - 1 || nDeltaI > 1 ||
nDeltaJ < - 1 || nDeltaJ > 1 ||
nDeltaK < - 1 || nDeltaK > 1)
if ( nDeltaI < -1 || nDeltaI > 1 ||
nDeltaJ < -1 || nDeltaJ > 1 ||
nDeltaK < -1 || nDeltaK > 1)
return false ;
// Determino blocco adiacente
@@ -3666,18 +3760,18 @@ VolZmap::IsAVoxelOnBoundary( const int nLimits[], const int nIJK[], bool bType)
int nVoxNumZ = int( m_nNy[1] / N_DEXVOXRATIO + ( m_nNy[1] % N_DEXVOXRATIO == 0 ? 1 : 2)) ;
// Test sulla validità dei limiti
if ( nLimits[0] < - 1 || nLimits[0] > nVoxNumX - 1 ||
nLimits[1] < - 1 || nLimits[1] > nVoxNumX - 1 ||
nLimits[2] < - 1 || nLimits[2] > nVoxNumY - 1 ||
nLimits[3] < - 1 || nLimits[3] > nVoxNumY - 1 ||
nLimits[4] < - 1 || nLimits[4] > nVoxNumZ - 1 ||
nLimits[5] < - 1 || nLimits[5] > nVoxNumZ - 1 )
if ( nLimits[0] < -1 || nLimits[0] > nVoxNumX - 1 ||
nLimits[1] < -1 || nLimits[1] > nVoxNumX - 1 ||
nLimits[2] < -1 || nLimits[2] > nVoxNumY - 1 ||
nLimits[3] < -1 || nLimits[3] > nVoxNumY - 1 ||
nLimits[4] < -1 || nLimits[4] > nVoxNumZ - 1 ||
nLimits[5] < -1 || nLimits[5] > nVoxNumZ - 1 )
return false ;
// Controllo sull'ammissibilità del voxel
if ( nIJK[0] <= -2 || nIJK[0] > nVoxNumX - 2 ||
nIJK[1] <= -2 || nIJK[1] > nVoxNumY - 2 ||
nIJK[2] <= -2 || nIJK[2] > nVoxNumZ - 2)
if ( nIJK[0] < -1 || nIJK[0] > nVoxNumX - 2 ||
nIJK[1] < -1 || nIJK[1] > nVoxNumY - 2 ||
nIJK[2] < -1 || nIJK[2] > nVoxNumZ - 2)
return false ;
// Se cerchiamo i voxel che sono sulla frontiera del blocco
@@ -3707,9 +3801,9 @@ VolZmap::IsAVoxelOnBoundary( const int nLimits[], const int nIJK[], bool bType)
// Determino (se esiste) il blocco in cui cade il voxel adiacente.
int nAdjBlockIJK[3] ;
if ( GetVoxelBlockIJK( nAdjIJK, nAdjBlockIJK)) {
if ( nAdjBlockIJK[0] > -1 && nAdjBlockIJK[0] < int( m_nFracLin[0]) &&
nAdjBlockIJK[1] > -1 && nAdjBlockIJK[1] < int( m_nFracLin[1]) &&
nAdjBlockIJK[2] > -1 && nAdjBlockIJK[2] < int( m_nFracLin[2]) &&
if ( nAdjBlockIJK[0] > -1 && nAdjBlockIJK[0] < m_nFracLin[0] &&
nAdjBlockIJK[1] > -1 && nAdjBlockIJK[1] < m_nFracLin[1] &&
nAdjBlockIJK[2] > -1 && nAdjBlockIJK[2] < m_nFracLin[2] &&
( nCurrentBlockIJK[0] != nAdjBlockIJK[0] ||
nCurrentBlockIJK[1] != nAdjBlockIJK[1] ||
nCurrentBlockIJK[2] != nAdjBlockIJK[2])) {
@@ -3732,17 +3826,17 @@ VolZmap::IsAVoxelOnBoundary( const int nLimits[], const int nIJK[], int nDeltaIn
int nVoxNumY = int( m_nNy[0] / N_DEXVOXRATIO + ( m_nNy[0] % N_DEXVOXRATIO == 0 ? 1 : 2)) ;
int nVoxNumZ = int( m_nNy[1] / N_DEXVOXRATIO + ( m_nNy[1] % N_DEXVOXRATIO == 0 ? 1 : 2)) ;
// Test sulla validità dei limiti
if ( nLimits[0] < - 1 || nLimits[0] > nVoxNumX - 1 ||
nLimits[1] < - 1 || nLimits[1] > nVoxNumX - 1 ||
nLimits[2] < - 1 || nLimits[2] > nVoxNumY - 1 ||
nLimits[3] < - 1 || nLimits[3] > nVoxNumY - 1 ||
nLimits[4] < - 1 || nLimits[4] > nVoxNumZ - 1 ||
nLimits[5] < - 1 || nLimits[5] > nVoxNumZ - 1 )
if ( nLimits[0] < -1 || nLimits[0] > nVoxNumX - 1 ||
nLimits[1] < -1 || nLimits[1] > nVoxNumX - 1 ||
nLimits[2] < -1 || nLimits[2] > nVoxNumY - 1 ||
nLimits[3] < -1 || nLimits[3] > nVoxNumY - 1 ||
nLimits[4] < -1 || nLimits[4] > nVoxNumZ - 1 ||
nLimits[5] < -1 || nLimits[5] > nVoxNumZ - 1 )
return false ;
// Controllo sull'ammissibilità del voxel
if ( nIJK[0] <= -2 || nIJK[0] > nVoxNumX - 2 ||
nIJK[1] <= -2 || nIJK[1] > nVoxNumY - 2 ||
nIJK[2] <= -2 || nIJK[2] > nVoxNumZ - 2)
if ( nIJK[0] < -1 || nIJK[0] > nVoxNumX - 2 ||
nIJK[1] < -1 || nIJK[1] > nVoxNumY - 2 ||
nIJK[2] < -1 || nIJK[2] > nVoxNumZ - 2)
return false ;
nDeltaIndex[0] = 0 ;
nDeltaIndex[1] = 0 ;