Egalware/EgtGeomKernel
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

EgtGeomKernel :

Browse Source 
- migliorie e correzioni ulteriori a calcolo pelle di Zmap.
This commit is contained in:
Dario Sassi
2017-09-18 11:24:00 +00:00
parent 38590bf381
commit 7a83e36819
2 changed files with 77 additions and 209 deletions
Download Patch File Download Diff File Expand all files Collapse all files
VolTriZmapGraphics.cpp
+77 -206
View File
@@ -124,7 +124,7 @@ CanonicPlaneTest( const VectorField CompoVert[], int nDir, double& dPos)
case Z_MINUS : vtN = - Z_AX ; break ;
}
for ( int i = 0 ; i < 4 ; ++ i) {
if ( CompoVert[i].vtNorm * vtN < 0.99)
if ( CompoVert[i].vtNorm * vtN < 0.999)
return false ;
}
// Superati tutti i test
@@ -695,16 +695,10 @@ VolZmap::MarchingCubes( int nBlock, TRIA3DLIST& lstTria) const
int i1 = TriangleTableEn[nIndex][0][TriIndex + 1] ;
int i2 = TriangleTableEn[nIndex][0][TriIndex] ;
Triangle3d CurrentTriangle ;
Vector3d vtN = ( ptIntPoint[i1] - ptIntPoint[i0]) ^ ( ptIntPoint[i2] - ptIntPoint[i1]) ;
vtN.Normalize() ;
vtN.ToGlob( m_MapFrame[0]) ;
// Il triangolo è pronto
CurrentTriangle.Set( ptIntPoint[i0], ptIntPoint[i1], ptIntPoint[i2], vtN) ;
Triangle3d CurrentTriangle ;
CurrentTriangle.Set( ptIntPoint[i0], ptIntPoint[i1], ptIntPoint[i2]) ;
bool bV = CurrentTriangle.Validate( true) ;
// Aggiungo triangolo
lstTria.emplace_back( CurrentTriangle) ;
@@ -754,10 +748,6 @@ VolZmap::ExtMarchingCubes( int nBlock, TRIA3DLIST& lstTria, TriHolder& triHold)
if ( EdgeTable[nIndex] == 0)
continue ;
// Riconoscimento dei voxel di frontiera
int nVoxIndexes[3] = { i, j, k} ;
bool bBoundary = IsAVoxelOnBoundary( nLimits, nVoxIndexes, true) ;
// Indici i,j,k dei vertici
int IndexCorner[8][3] = {
{ i, j, k},
@@ -778,8 +768,7 @@ VolZmap::ExtMarchingCubes( int nBlock, TRIA3DLIST& lstTria, TriHolder& triHold)
// Array di strutture punto di intersezione e normale alla superficie in esso.
VectorField VecField[12] ;
// Flag sulla regolatrità dei campi scalare e vettoriale:
// vero se i campi sono regolari
// Flag di regolarità dei campi scalare e vettoriale
bool bReg = true ;
// Ciclo sui segmenti
@@ -787,17 +776,15 @@ VolZmap::ExtMarchingCubes( int nBlock, TRIA3DLIST& lstTria, TriHolder& triHold)
// Se il segmento non attraversa la superficie passo al successivo
if ( ( EdgeTable[nIndex] & ( 1 << EdgeIndex)) == 0)
continue ;
// Indici per linee di griglia sui vertici
int n1 = intersections[EdgeIndex][0] ;
int n2 = intersections[EdgeIndex][1] ;
// Flag posizione corner
bool bN1 = ( ( nIndex & ( 1 << n1)) != 0) ;
// Determino con precisione il punto di intersezione sullo spigolo,
// se i campi scalare e vettoriale non sono regolari bReg diviene falso.
if ( ! IntersPos( IndexCorner[n1], IndexCorner[n2], bN1,
VecField[EdgeIndex].ptInt,
VecField[EdgeIndex].vtNorm))
VecField[EdgeIndex].ptInt, VecField[EdgeIndex].vtNorm))
bReg = false ;
}
@@ -1024,6 +1011,7 @@ VolZmap::ExtMarchingCubes( int nBlock, TRIA3DLIST& lstTria, TriHolder& triHold)
// Il triangolo è pronto
CurrentTriangle.Set( VecField[i0].ptInt, VecField[i1].ptInt, VecField[i2].ptInt) ;
bool bV = CurrentTriangle.Validate( true) ;
CurrentTriangle.ToGlob( m_MapFrame[0]) ;
// Aggiungo alla lista
lstTria.emplace_back( CurrentTriangle) ;
@@ -1792,6 +1780,10 @@ VolZmap::ExtMarchingCubes( int nBlock, TRIA3DLIST& lstTria, TriHolder& triHold)
// non piana confermo ExtMC
if ( ! ( bOutside || bPlane) && ! bDangerInversion) {
// Riconoscimento se voxel di frontiera
int nVoxIndexes[3] = { i, j, k} ;
bool bBoundary = IsAVoxelOnBoundary( nLimits, nVoxIndexes, true) ;
TRIA3DVECTOR vInnerTriaTemp, vBorderTriaTemp ;
for ( int ni = 0 ; ni < nContSize ; ++ ni) {
@@ -1805,7 +1797,7 @@ VolZmap::ExtMarchingCubes( int nBlock, TRIA3DLIST& lstTria, TriHolder& triHold)
Triangle3d trTemp = triContainer[ni] ;
if ( ! IsATriangleOnBorder( trTemp, ptTemp, nLimits, nVoxIJK, bBoundary))
if ( ! bBoundary || ! IsATriangleOnBorder( trTemp, ptTemp, nLimits, nVoxIJK))
vInnerTriaTemp.emplace_back( triContainer[ni]) ;
else
vBorderTriaTemp.emplace_back( triContainer[ni]) ;
@@ -1813,7 +1805,7 @@ VolZmap::ExtMarchingCubes( int nBlock, TRIA3DLIST& lstTria, TriHolder& triHold)
}
// Porto la soluzione nel sistema in cui è immerso lo Zmap
ptSol.ToGlob( m_MapFrame[0]) ;
ptSol.ToGlob( m_MapFrame[0]) ;
// Metto i triangoli negli appositi contenitori:
@@ -1869,7 +1861,7 @@ VolZmap::ExtMarchingCubes( int nBlock, TRIA3DLIST& lstTria, TriHolder& triHold)
// Se siamo in presenza della prima feature del
// voxel, ridimensiono il vettore che contiene
// la struttura dati del voxel.
if ( nBorderFeatureInVoxel == 1) {
if ( nBorderFeatureInVoxel == 1) {
m_InterBlockTria[nBlock].resize( m_InterBlockTria[nBlock].size() + 1) ;
@@ -1912,7 +1904,7 @@ VolZmap::ExtMarchingCubes( int nBlock, TRIA3DLIST& lstTria, TriHolder& triHold)
CompoTriVert[nCompCount - 1][TriIndex+1].ptInt,
CompoTriVert[nCompCount - 1][TriIndex+2].ptInt) ;
bool bV = CurrentTriangle.Validate( true) ;
// Riporto lew coordinate nel sistema in cui è immerso lo Zmap
// Riporto le coordinate nel sistema in cui è immerso lo Zmap
CurrentTriangle.ToGlob( m_MapFrame[0]) ;
// Aggiungo alla lista
lstTria.emplace_back( CurrentTriangle) ;
@@ -1930,8 +1922,8 @@ VolZmap::ExtMarchingCubes( int nBlock, TRIA3DLIST& lstTria, TriHolder& triHold)
CompoTriVert[nCompCount - 1][TriIndex+1].ptInt,
CompoTriVert[nCompCount - 1][TriIndex+2].ptInt) ;
bool bV = CurrentTriangle.Validate( true) ;
// Riporto lew coordinate nel sistema in cui è immerso lo Zmap
CurrentTriangle.ToGlob( m_MapFrame[0]) ;
// Riporto le coordinate nel sistema in cui è immerso lo Zmap
CurrentTriangle.ToGlob( m_MapFrame[0]) ;
// Aggiungo alla lista
lstTria.emplace_back( CurrentTriangle) ;
}
@@ -2361,7 +2353,7 @@ VolZmap::IntersPos( int nVec1[], int nVec2[], Point3d& ptInt) const
bool
VolZmap::IntersPos( int nVec1[], int nVec2[], bool bFirstCorner, Point3d& ptInt, Vector3d& vtNormal) const
{
double dEps = 2 * EPS_SMALL ;
const double dEps = 2 * EPS_SMALL ;
bool bFound = false ;
@@ -2377,53 +2369,41 @@ VolZmap::IntersPos( int nVec1[], int nVec2[], bool bFirstCorner, Point3d& ptInt,
ptInt.z = ( nVec1[2] + 0.5) * m_dStep ;
size_t nDexel = nVec1[2] * m_nNx[1] + nVec1[1] ;
size_t nSize = m_Values[1][nDexel].size() ;
int nSize = int( m_Values[1][nDexel].size()) ;
if ( bFirstCorner) {
int n = int( nSize) - 1 ;
int n = nSize - 1 ;
double dX = m_Values[1][nDexel][n].dMax ;
while ( n >= 0 && dX > dMinX - dEps) {
while ( n >= 0 && dX > dMinX - dEps) {
if ( dX < dMaxX + dEps) {
ptInt.x = dX ;
vtNormal = m_Values[1][nDexel][n].vtMaxN ;
bFound = true ;
break ;
}
n -= 1 ;
if ( n >= 0)
dX = m_Values[1][nDexel][n].dMax ;
}
}
else {
size_t n = 0 ;
int n = 0 ;
double dX = m_Values[1][nDexel][n].dMin ;
while ( n < nSize && dX < dMaxX + dEps) {
while ( n < nSize && dX < dMaxX + dEps) {
if ( dX > dMinX - dEps) {
ptInt.x = dX ;
vtNormal = m_Values[1][nDexel][n].vtMinN ;
bFound = true ;
break ;
}
n += 1 ;
if ( n < nSize)
dX = m_Values[1][nDexel][n].dMin ;
}
}
if ( ! bFound)
ptInt.x = 0.5 * ( dMinX + dMaxX) ;
}
@@ -2439,12 +2419,12 @@ VolZmap::IntersPos( int nVec1[], int nVec2[], bool bFirstCorner, Point3d& ptInt,
ptInt.z = ( nVec1[2] + 0.5) * m_dStep ;
size_t nDexel = nVec1[0] * m_nNx[2] + nVec1[2] ;
size_t nSize = m_Values[2][nDexel].size() ;
int nSize = int( m_Values[2][nDexel].size()) ;
if ( bFirstCorner) {
int n = int( nSize) - 1 ;
int n = nSize - 1 ;
double dY = m_Values[2][nDexel][n].dMax ;
while ( n >= 0 && dY > dMinY - dEps) {
while ( n >= 0 && dY > dMinY - dEps) {
if ( dY < dMaxY + dEps) {
ptInt.y = dY ;
vtNormal = m_Values[2][nDexel][n].vtMaxN ;
@@ -2458,9 +2438,9 @@ VolZmap::IntersPos( int nVec1[], int nVec2[], bool bFirstCorner, Point3d& ptInt,
}
else {
size_t n = 0 ;
int n = 0 ;
double dY = m_Values[2][nDexel][n].dMin ;
while ( n < nSize && dY < dMaxY + dEps) {
while ( n < nSize && dY < dMaxY + dEps) {
if ( dY > dMinY - dEps) {
ptInt.y = dY ;
vtNormal = m_Values[2][nDexel][n].vtMinN ;
@@ -2489,10 +2469,10 @@ VolZmap::IntersPos( int nVec1[], int nVec2[], bool bFirstCorner, Point3d& ptInt,
ptInt.y = ( nVec1[1] + 0.5) * m_dStep ;
size_t nDexel = nVec1[1] * m_nNx[0] + nVec1[0] ;
size_t nSize = m_Values[0][nDexel].size() ;
int nSize = int( m_Values[0][nDexel].size()) ;
if ( bFirstCorner) {
int n = int( nSize) - 1 ;
int n = nSize - 1 ;
double dZ = m_Values[0][nDexel][n].dMax ;
while ( n >= 0 && dZ > dMinZ - dEps) {
if ( dZ < dMaxZ + dEps) {
@@ -2508,9 +2488,9 @@ VolZmap::IntersPos( int nVec1[], int nVec2[], bool bFirstCorner, Point3d& ptInt,
}
else {
size_t n = 0 ;
int n = 0 ;
double dZ = m_Values[0][nDexel][n].dMin ;
while ( n < nSize && dZ < dMaxZ + dEps) {
while ( n < nSize && dZ < dMaxZ + dEps) {
if ( dZ > dMinZ - dEps) {
ptInt.z = dZ ;
vtNormal = m_Values[0][nDexel][n].vtMinN ;
@@ -2627,47 +2607,25 @@ VolZmap::GetBlockLimitsIJK( const int nIJK[], int nLimits[]) const
return true ;
}
//----------------------------------------------------------------------------
bool
VolZmap::IsPointInsideVoxel( int nI, int nJ, int nK, const Point3d& ptP) const
{
// Controllo sull'ammissibilità del voxel
if ( nI < - 1 || nI >= int( m_nNx[0]) ||
nJ < - 1 || nJ >= int( m_nNy[0]) ||
nK < - 1 || nK >= int( m_nNy[1]))
return false ;
int nPointI = int( floor( ( ptP.x - 0.5 * m_dStep) / m_dStep)) ;
int nPointJ = int( floor( ( ptP.y - 0.5 * m_dStep) / m_dStep)) ;
int nPointK = int( floor( ( ptP.z - 0.5 * m_dStep) / m_dStep)) ;
return ( nPointI == nI &&
nPointJ == nJ &&
nPointK == nK) ;
}
//----------------------------------------------------------------------------
bool
VolZmap::IsPointInsideVoxelApprox( int nI, int nJ, int nK, const Point3d& ptP, double dPrec) const
{
// Controllo sull'ammissibilità del voxel
if ( nI < - 1 || nI >= int( m_nNx[0]) ||
nJ < - 1 || nJ >= int( m_nNy[0]) ||
nK < - 1 || nK >= int( m_nNy[1]))
if ( nI < - 1 || nI >= int( m_nNx[0]) ||
nJ < - 1 || nJ >= int( m_nNy[0]) ||
nK < - 1 || nK >= int( m_nNy[1]))
return false ;
if ( dPrec < EPS_ZERO)
dPrec = 0 ;
Point3d ptPZmap = ptP ;
ptPZmap.ToLoc( m_MapFrame[0]) ;
bool bI = ptPZmap.x > ( nI + 0.5) * m_dStep - dPrec &&
ptPZmap.x < ( nI + 1.5) * m_dStep + dPrec ;
bool bJ = ptPZmap.y > ( nJ + 0.5) * m_dStep - dPrec &&
ptPZmap.y < ( nJ + 1.5) * m_dStep + dPrec ;
bool bK = ptPZmap.z > ( nK + 0.5) * m_dStep - dPrec &&
ptPZmap.z < ( nK + 1.5) * m_dStep + dPrec ;
bool bI = ptP.x > ( nI + 0.5) * m_dStep - dPrec &&
ptP.x < ( nI + 1.5) * m_dStep + dPrec ;
bool bJ = ptP.y > ( nJ + 0.5) * m_dStep - dPrec &&
ptP.y < ( nJ + 1.5) * m_dStep + dPrec ;
bool bK = ptP.z > ( nK + 0.5) * m_dStep - dPrec &&
ptP.z < ( nK + 1.5) * m_dStep + dPrec ;
return ( bI && bJ && bK) ;
}
@@ -2747,141 +2705,54 @@ VolZmap::IsAVoxelOnBoundary( const int nLimits[], const int nIJK[], bool bType)
nIJK[2] <= -2 || nIJK[2] >= int( m_nNy[1]))
return false ;
// Cerchiamo i voxel che confinano con quelli di altri blocchi
if ( bType) {
// Se cerchiamo i voxel che sono sulla frontiera del blocco
if ( ! bType) {
return ( nIJK[0] == nLimits[0] || nIJK[0] == nLimits[1] - 1 ||
nIJK[1] == nLimits[2] || nIJK[1] == nLimits[3] - 1 ||
nIJK[2] == nLimits[4] || nIJK[2] == nLimits[5] - 1) ;
}
// Altrimenti cerchiamo i voxel che confinano con quelli di altri blocchi
// Condizione necessaria è che il voxel sia sulla frontiera
if ( nIJK[0] == nLimits[0] || nIJK[0] == nLimits[1] - 1 ||
nIJK[1] == nLimits[2] || nIJK[1] == nLimits[3] - 1 ||
nIJK[2] == nLimits[4] || nIJK[2] == nLimits[5] - 1) {
// Indici del blocco
int nCurrentBlockIJK[3] ;
GetVoxelBlockIJK( nIJK, nCurrentBlockIJK) ;
// Condizione necessaria è che il voxel sia sulla frontiera
if ( IsAVoxelOnBoundary( nLimits, nIJK, false)) {
// Ciclo sulle posizioni possibili dei voxel adiacenti
for ( int nDeltaI = -1 ; nDeltaI <= 1 ; ++ nDeltaI) {
for ( int nDeltaJ = -1 ; nDeltaJ <= 1 ; ++ nDeltaJ) {
for ( int nDeltaK = -1 ; nDeltaK <= 1 ; ++ nDeltaK) {
if ( nDeltaI == 0 && nDeltaJ == 0 && nDeltaK == 0)
continue ;
// Indici dei voxel adiacenti
int nAdjIJK[3] = { nIJK[0] + nDeltaI, nIJK[1] + nDeltaJ, nIJK[2] + nDeltaK} ;
// 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]) &&
( nCurrentBlockIJK[0] != nAdjBlockIJK[0] ||
nCurrentBlockIJK[1] != nAdjBlockIJK[1] ||
nCurrentBlockIJK[2] != nAdjBlockIJK[2])) {
return true ;
}
// Ciclo sulle possibili posizioni dei voxel adiacenti
for ( int nDeltaI = -1 ; nDeltaI <= 1 ; ++ nDeltaI) {
for ( int nDeltaJ = -1 ; nDeltaJ <= 1 ; ++ nDeltaJ) {
for ( int nDeltaK = -1 ; nDeltaK <= 1 ; ++ nDeltaK) {
// Evito controllo con se stesso
if ( nDeltaI == 0 && nDeltaJ == 0 && nDeltaK == 0)
continue ;
// Indici del voxel adiacente
int nAdjIJK[3] = { nIJK[0] + nDeltaI, nIJK[1] + nDeltaJ, nIJK[2] + nDeltaK} ;
// 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]) &&
( nCurrentBlockIJK[0] != nAdjBlockIJK[0] ||
nCurrentBlockIJK[1] != nAdjBlockIJK[1] ||
nCurrentBlockIJK[2] != nAdjBlockIJK[2])) {
return true ;
}
}
}
}
}
}
// altrimenti cerchiamo i voxel che sono sulla frontiera del blocco
else {
if ( nIJK[0] == nLimits[0] || nIJK[0] == nLimits[1] - 1 ||
nIJK[1] == nLimits[2] || nIJK[1] == nLimits[3] - 1 ||
nIJK[2] == nLimits[4] || nIJK[2] == nLimits[5] - 1)
return true ;
}
return false ;
}
//----------------------------------------------------------------------------
bool
VolZmap::IsATriangleOnBorder( const Triangle3d& trTria, const Point3d& ptVert,
const int nBlockLimits[], const int nVoxIJK[]) const
{
// Se il triangolo ha area nulla non proseguiamo
if ( trTria.GetArea() < EPS_SMALL)
return false ;
// Se il voxel non è sulla frontiera esco
if ( ! IsAVoxelOnBoundary( nBlockLimits, nVoxIJK, true))
return false ;
Point3d ptFirstGrPt, ptSecondGrPt ;
int nCount = 0 ;
if ( SqDist( trTria.GetP( 0), ptVert) > SQ_EPS_SMALL) {
ptFirstGrPt = trTria.GetP( 0) ;
nCount ++ ;
}
if ( SqDist( trTria.GetP( 1), ptVert) > SQ_EPS_SMALL) {
if ( nCount == 0) {
ptFirstGrPt = trTria.GetP( 1) ;
nCount ++ ;
}
else if ( nCount == 1) {
ptSecondGrPt = trTria.GetP( 1) ;
nCount ++ ;
}
}
if ( SqDist( trTria.GetP( 2), ptVert) > SQ_EPS_SMALL) {
if ( nCount == 1) {
ptSecondGrPt = trTria.GetP( 2) ;
nCount ++ ;
}
}
// Se non ho trovato due punti, errore
if ( nCount != 2)
return false ;
// Verifico se stanno sulla griglia
if ( nVoxIJK[0] == nBlockLimits[0]) {
double dXGrid = ( nBlockLimits[0] + 0.5) * m_dStep ;
if ( abs( ptFirstGrPt.x - dXGrid) < EPS_SMALL &&
abs( ptSecondGrPt.x - dXGrid) < EPS_SMALL)
return true ;
}
if ( nVoxIJK[1] == nBlockLimits[2]) {
double dYGrid = ( nBlockLimits[2] + 0.5) * m_dStep ;
if ( abs( ptFirstGrPt.y - dYGrid) < EPS_SMALL &&
abs( ptSecondGrPt.y - dYGrid) < EPS_SMALL)
return true ;
}
if ( nVoxIJK[2] == nBlockLimits[4]) {
double dZGrid = ( nBlockLimits[4] + 0.5) * m_dStep ;
if ( abs( ptFirstGrPt.z - dZGrid) < EPS_SMALL &&
abs( ptSecondGrPt.z - dZGrid) < EPS_SMALL)
return true ;
}
if ( nVoxIJK[0] + 1 == nBlockLimits[1]) {
double dXGrid = ( nBlockLimits[1] + 0.5) * m_dStep ;
if ( abs( ptFirstGrPt.x - dXGrid) < EPS_SMALL &&
abs( ptSecondGrPt.x - dXGrid) < EPS_SMALL)
return true ;
}
if ( nVoxIJK[1] + 1 == nBlockLimits[3]) {
double dYGrid = ( nBlockLimits[3] + 0.5) * m_dStep ;
if ( abs( ptFirstGrPt.y - dYGrid) < EPS_SMALL &&
abs( ptSecondGrPt.y - dYGrid) < EPS_SMALL)
return true ;
}
if ( nVoxIJK[2] + 1 == nBlockLimits[5]) {
double dZGrid = ( nBlockLimits[5] + 0.5) * m_dStep ;
if ( abs( ptFirstGrPt.z - dZGrid) < EPS_SMALL &&
abs( ptSecondGrPt.z - dZGrid) < EPS_SMALL)
return true ;
}
return false ;
}
//----------------------------------------------------------------------------
bool
VolZmap::IsATriangleOnBorder( const Triangle3d& trTria, const Point3d& ptVert,
const int nBlockLimits[], const int nVoxIJK[], bool bBorderBox) const
const int nBlockLimits[], const int nVoxIJK[]) const
{
// Se il voxel non è di frontiera, esco
if ( ! bBorderBox)
return false ;
// Se il triangolo ha area nulla, esco
if ( trTria.GetArea() < EPS_SMALL)
return false ;