EgtGeomKernel :

- correzioni a triangolazione superficie Zmap per congruenza scelta topologia tra voxel adiacenti.
This commit is contained in:
Dario Sassi
2018-05-22 11:16:47 +00:00
parent 1cf91116aa
commit 6d82a47214
2 changed files with 302 additions and 133 deletions
+275 -130
View File
@@ -1,7 +1,7 @@
//----------------------------------------------------------------------------
// EgalTech 2015-2018
//----------------------------------------------------------------------------
// File : VolZmap.cpp Data : 27.01.18 Versione : 1.8a3
// File : VolZmap.cpp Data : 22.05.18 Versione : 1.9e4
// Contenuto : Implementazione della classe Volume Zmap (tre griglie)
//
//
@@ -1220,12 +1220,18 @@ VolZmap::ExtMarchingCubes( int nBlock, TRIA3DEXLIST& lstTria, VoxelContainer& vV
FlatVoxelContainer VoxContXYSup( nDim) ;
FlatVoxelContainer VoxContYZInf( nDim) ;
FlatVoxelContainer VoxContYZSup( nDim) ;
// Unordered map per la coerenza topologica
typedef std::unordered_map <int, bool> InterVoxMatter ;
InterVoxMatter SliceXY( 200) ;
InterVoxMatter SliceXZ( 200) ;
InterVoxMatter SliceYZ( 200) ;
// Ciclo su tutti i voxel del blocco
for ( int i = nLimits[0] ; i < nLimits[1] ; ++ i) {
for ( int j = nLimits[2] ; j < nLimits[3] ; ++ j) {
for ( int k = nLimits[4] ; k < nLimits[5] ; ++ k) {
// Classificazione dei vertici: interni o esterni al materiale
int nIndex = CalcIndex( i, j, k) ;
@@ -1273,6 +1279,7 @@ VolZmap::ExtMarchingCubes( int nBlock, TRIA3DEXLIST& lstTria, VoxelContainer& vV
bReg = false ;
}
// Determino il numero di componenti connesse nel voxel in caso di configurazione standard
int nComponents = TriangleTableEn[nIndex][1][0] ;
@@ -1367,13 +1374,44 @@ VolZmap::ExtMarchingCubes( int nBlock, TRIA3DEXLIST& lstTria, VoxelContainer& vV
continue ;
}
}
// Test sulla topologia: richiede campo regolare perchè si fonda su angoli tra normali
if ( bReg) {
// Configurazione 3
if ( nAllConfig[nIndex] == 3) {
// Test sulla topologia
// Configurazione 3
if ( nAllConfig[nIndex] == 3) {
// Test sulla topologia
bool bDefTopology = false ;
bool bNewTopology = false ;
int nCount = 0 ;
while ( nIndexConfig3[nCount] != nIndex)
++ nCount ;
int nISl = ( nAdjVox3[nCount] != 1 ? i : i + 1) ;
int nJSl = ( nAdjVox3[nCount] != 2 ? j : j + 1) ;
int nKSl = ( nAdjVox3[nCount] != 3 ? k : k + 1) ;
int nNSl ;
if ( GetVoxNFromIJK( nISl, nJSl, nKSl, nNSl)) {
if ( abs( nAdjVox3[nCount]) == 1) {
auto it = SliceYZ.find( nNSl) ;
if ( it != SliceYZ.end()) {
bNewTopology = it->second ;
bDefTopology = true ;
}
}
else if ( abs( nAdjVox3[nCount]) == 2) {
auto it = SliceXZ.find( nNSl) ;
if ( it != SliceXZ.end()) {
bNewTopology = it->second ;
bDefTopology = true ;
}
}
else if ( abs( nAdjVox3[nCount]) == 3) {
auto it = SliceXY.find( nNSl) ;
if ( it != SliceXY.end()) {
bNewTopology = it->second ;
bDefTopology = true ;
}
}
}
if ( ! bDefTopology && bReg) {
double dDotSum = 0 ;
for ( int nFV = 0 ; nFV < 3 ; ++ nFV) {
for ( int nTV = 0 ; nTV < 2 ; ++ nTV) {
@@ -1390,42 +1428,85 @@ VolZmap::ExtMarchingCubes( int nBlock, TRIA3DEXLIST& lstTria, VoxelContainer& vV
dDotSum -= CompoVert[1][nTVI].vtVec * CompoVert[1][nTVJ].vtVec ;
}
}
bool bTestOnSum = dDotSum > - EPS_SMALL ;
// Si passa alla seconda topologia
if ( bTestOnSum) {
// Ricerca del caso corrispondente della nuova topologia
int nt = 0 ;
while ( nIndexVsIndex3[nt][0] != nIndex)
++ nt ;
int nRotCase = nIndexVsIndex3[nt][1] ;
// Aggiorno numero di componenti
nComponents = Cases3Plus[nRotCase][1][0] ;
// Riaggiorno gli offsets
nExtTabOff = nComponents ;
nStdTabOff = 0 ;
// Modifico le matrici
for ( int nC = 1 ; nC <= nComponents ; ++ nC) {
// Numero vertici per componenti
nVertComp[nC - 1] = Cases3Plus[nRotCase][1][nC] ;
// Matrice dei vertici della base del fan
for ( int nFanVert = 0 ; nFanVert < nVertComp[nC - 1] ; ++ nFanVert)
CompoVert[nC - 1][nFanVert] = VecField[Cases3Plus[nRotCase][1][nFanVert + nExtTabOff + 1]] ;
// Matrici dei vertici dei triangoli in assenza di sharp feature
for ( int nTriVert = 0 ; nTriVert < 3 * ( nVertComp[nC - 1] - 2) ; nTriVert += 3) {
CompoTriVert[nC - 1][nTriVert] = VecField[Cases3Plus[nRotCase][0][nStdTabOff + nTriVert+2]] ;
CompoTriVert[nC - 1][nTriVert+1] = VecField[Cases3Plus[nRotCase][0][nStdTabOff + nTriVert+1]] ;
CompoTriVert[nC - 1][nTriVert+2] = VecField[Cases3Plus[nRotCase][0][nStdTabOff + nTriVert]] ;
}
// Aggiorno gli offsets per raggiungere i vertici della componente successiva.
nExtTabOff += nVertComp[nC - 1] ;
nStdTabOff += 3 * ( nVertComp[nC - 1] - 2) ;
}
bNewTopology = dDotSum > - EPS_SMALL ;
if ( GetVoxNFromIJK( nISl, nJSl, nKSl, nNSl)) {
if ( abs( nAdjVox3[nCount]) == 1)
SliceYZ.emplace( nNSl, ! bNewTopology) ;
else if ( abs( nAdjVox3[nCount]) == 2)
SliceXZ.emplace( nNSl, ! bNewTopology) ;
else if ( abs( nAdjVox3[nCount]) == 3)
SliceXY.emplace( nNSl, ! bNewTopology) ;
}
}
// Si passa alla seconda topologia
if ( bNewTopology) {
// Ricerca del caso corrispondente della nuova topologia
int nt = 0 ;
while ( nIndexVsIndex3[nt][0] != nIndex)
++ nt ;
int nRotCase = nIndexVsIndex3[nt][1] ;
// Aggiorno numero di componenti
nComponents = Cases3Plus[nRotCase][1][0] ;
// Riaggiorno gli offsets
nExtTabOff = nComponents ;
nStdTabOff = 0 ;
// Modifico le matrici
for ( int nC = 1 ; nC <= nComponents ; ++ nC) {
// Numero vertici per componenti
nVertComp[nC - 1] = Cases3Plus[nRotCase][1][nC] ;
// Matrice dei vertici della base del fan
for ( int nFanVert = 0 ; nFanVert < nVertComp[nC - 1] ; ++ nFanVert)
CompoVert[nC - 1][nFanVert] = VecField[Cases3Plus[nRotCase][1][nFanVert + nExtTabOff + 1]] ;
// Matrici dei vertici dei triangoli in assenza di sharp feature
for ( int nTriVert = 0 ; nTriVert < 3 * ( nVertComp[nC - 1] - 2) ; nTriVert += 3) {
CompoTriVert[nC - 1][nTriVert] = VecField[Cases3Plus[nRotCase][0][nStdTabOff + nTriVert+2]] ;
CompoTriVert[nC - 1][nTriVert+1] = VecField[Cases3Plus[nRotCase][0][nStdTabOff + nTriVert+1]] ;
CompoTriVert[nC - 1][nTriVert+2] = VecField[Cases3Plus[nRotCase][0][nStdTabOff + nTriVert]] ;
}
// Aggiorno gli offsets per raggiungere i vertici della componente successiva.
nExtTabOff += nVertComp[nC - 1] ;
nStdTabOff += 3 * ( nVertComp[nC - 1] - 2) ;
}
}
}
// Configurazione 6
else if ( nAllConfig[nIndex] == 6) {
// Test sulla topologia
// Configurazione 6
else if ( nAllConfig[nIndex] == 6) {
// Test sulla topologia
bool bDefTopology = false ;
bool bNewTopology = false ;
int nCount = 0 ;
while ( nIndexConfig6[nCount] != nIndex)
++ nCount ;
int nISl = ( nAdjVox6[nCount] != 1 ? i : i + 1) ;
int nJSl = ( nAdjVox6[nCount] != 2 ? j : j + 1) ;
int nKSl = ( nAdjVox6[nCount] != 3 ? k : k + 1) ;
int nNSl ;
if ( GetVoxNFromIJK( nISl, nJSl, nKSl, nNSl)) {
if ( abs( nAdjVox6[nCount]) == 1) {
auto it = SliceYZ.find( nNSl) ;
if ( it != SliceYZ.end()) {
bNewTopology = it->second ;
bDefTopology = true ;
}
}
else if ( abs( nAdjVox6[nCount]) == 2) {
auto it = SliceXZ.find( nNSl) ;
if ( it != SliceXZ.end()) {
bNewTopology = it->second ;
bDefTopology = true ;
}
}
else if ( abs( nAdjVox6[nCount]) == 3) {
auto it = SliceXY.find( nNSl) ;
if ( it != SliceXY.end()) {
bNewTopology = it->second ;
bDefTopology = true ;
}
}
}
if ( ! bDefTopology && bReg) {
// Test sulla topologia
double dDotSum = 0 ;
for ( int nFV = 0 ; nFV < 4 ; ++ nFV) {
for ( int nTV = 0 ; nTV < 3 ; ++ nTV) {
@@ -1442,100 +1523,164 @@ VolZmap::ExtMarchingCubes( int nBlock, TRIA3DEXLIST& lstTria, VoxelContainer& vV
dDotSum -= CompoVert[1][nTVI].vtVec * CompoVert[1][nTVJ].vtVec ;
}
}
bool bTestOnSum = dDotSum > - 4 ;
// Si deve passare alla seconda topologia
if ( bTestOnSum) {
// Ricerca del caso corrispondente della nuova topologia
int nt = 0 ;
while ( nIndexVsIndex6[nt][0] != nIndex)
++ nt ;
int nRotCase = nIndexVsIndex6[nt][1] ;
// Costruzione dei triangoli
for ( int TriIndex = 0 ; TriIndex < 15 ; TriIndex += 3) {
// Indici vertici
int i0 = Cases6Plus[nRotCase][TriIndex + 2] ;
int i1 = Cases6Plus[nRotCase][TriIndex + 1] ;
int i2 = Cases6Plus[nRotCase][TriIndex] ;
// Costruzione triangolo
Triangle3dEx CurrentTriangle ;
CurrentTriangle.Set( VecField[i0].ptPApp, VecField[i1].ptPApp, VecField[i2].ptPApp) ;
CurrentTriangle.Validate( true) ;
CurrentTriangle.SetVertexNorm( 0, VecField[i0].vtVec) ;
CurrentTriangle.SetVertexNorm( 1, VecField[i1].vtVec) ;
CurrentTriangle.SetVertexNorm( 2, VecField[i2].vtVec) ;
// Setto il numero di utensile (conta solo positivo, nullo o negativo)
int nTool0 = Clamp( VecField[i0].nPropIndex, -1, 1) ;
int nTool1 = Clamp( VecField[i1].nPropIndex, -1, 1) ;
int nTool2 = Clamp( VecField[i2].nPropIndex, -1, 1) ;
if ( nTool0 == nTool1 || nTool0 == nTool2)
CurrentTriangle.SetGrade( nTool0) ;
else if ( nTool1 == nTool2)
CurrentTriangle.SetGrade( nTool1) ;
// Valido il triangolo e setto le normali del campo vettoriale ai corrispondenti vertici
if ( CurrentTriangle.Validate( true)) {
Vector3d vtVertNorm = VecField[i0].vtVec ;
if ( CurrentTriangle.GetN() * vtVertNorm > 0.6)
CurrentTriangle.SetVertexNorm( 0, vtVertNorm) ;
else
CurrentTriangle.SetVertexNorm( 0, CurrentTriangle.GetN()) ;
vtVertNorm = VecField[i1].vtVec ;
if ( CurrentTriangle.GetN() * vtVertNorm > 0.6)
CurrentTriangle.SetVertexNorm( 1, vtVertNorm) ;
else
CurrentTriangle.SetVertexNorm( 1, CurrentTriangle.GetN()) ;
vtVertNorm = VecField[i2].vtVec ;
if ( CurrentTriangle.GetN() * vtVertNorm > 0.6)
CurrentTriangle.SetVertexNorm( 2, vtVertNorm) ;
else
CurrentTriangle.SetVertexNorm( 2, CurrentTriangle.GetN()) ;
}
// Riporto il triangolo nel sistema in cui è immerso quello dello Zmap
CurrentTriangle.ToGlob( m_MapFrame) ;
// Aggiungo alla lista
lstTria.emplace_back( CurrentTriangle) ;
}
continue ;
bNewTopology = dDotSum > - 4 ;
if ( GetVoxNFromIJK( nISl, nJSl, nKSl, nNSl)) {
if ( abs( nAdjVox6[nCount]) == 1)
SliceYZ.emplace( nNSl, ! bNewTopology) ;
else if ( abs( nAdjVox6[nCount]) == 2)
SliceXZ.emplace( nNSl, ! bNewTopology) ;
else if ( abs( nAdjVox6[nCount]) == 3)
SliceXY.emplace( nNSl, ! bNewTopology) ;
}
}
// Configurazione 10
else if ( nAllConfig[nIndex] == 10) {
// Verifico concordanza tra i versori di una stessa componente
// (ogni coppia di vettori di una medesima componente deve avere prodotto scalare non inferiore a 0.0)
Vector3d vtCmpAvg0, vtCmpAvg1 ;
bool bTest0 = DotTest( CompoVert[0], 4, vtCmpAvg0, 0.0) ;
bool bTest1 = DotTest( CompoVert[1], 4, vtCmpAvg1, 0.0) ;
// Si passa alla seconda topologia
if ( ! bTest0 || ! bTest1) {
// Ricerca del caso corrispondente della nuova topologia
int nt = 0 ;
while ( nIndexVsIndex10[nt][0] != nIndex)
++ nt ;
// Riaggiorno gli offsets
nExtTabOff = 2 ;
nStdTabOff = 0 ;
// Modifico le matrici
int nRotCase = nIndexVsIndex10[nt][1] ;
for ( int nC = 1 ; nC <= 2 ; ++ nC) {
// Numero vertici per componenti
nVertComp[nC - 1] = Cases10Plus[nRotCase][1][nC] ;
// Matrice dei vertici della base del fan
for ( int nFanVert = 0 ; nFanVert < 4 ; ++ nFanVert)
CompoVert[nC - 1][nFanVert] = VecField[Cases10Plus[nRotCase][1][nFanVert + nExtTabOff + 1]] ;
// Matrici dei vertici dei triangoli in assenza di sharp feature
for ( int nTriVert = 0 ; nTriVert < 6 ; nTriVert += 3) {
CompoTriVert[nC - 1][nTriVert] = VecField[Cases10Plus[nRotCase][0][nStdTabOff + nTriVert+2]] ;
CompoTriVert[nC - 1][nTriVert+1] = VecField[Cases10Plus[nRotCase][0][nStdTabOff + nTriVert+1]] ;
CompoTriVert[nC - 1][nTriVert+2] = VecField[Cases10Plus[nRotCase][0][nStdTabOff + nTriVert]] ;
}
// Aggiorno gli offsets per raggiungere i vertici della componente successiva.
nExtTabOff += nVertComp[nC - 1] ;
nStdTabOff += 3 * ( nVertComp[nC - 1] - 2) ;
// Si deve passare alla seconda topologia
if ( bNewTopology) {
// Ricerca del caso corrispondente della nuova topologia
int nt = 0 ;
while ( nIndexVsIndex6[nt][0] != nIndex)
++ nt ;
int nRotCase = nIndexVsIndex6[nt][1] ;
// Costruzione dei triangoli
for ( int TriIndex = 0 ; TriIndex < 15 ; TriIndex += 3) {
// Indici vertici
int i0 = Cases6Plus[nRotCase][TriIndex + 2] ;
int i1 = Cases6Plus[nRotCase][TriIndex + 1] ;
int i2 = Cases6Plus[nRotCase][TriIndex] ;
// Costruzione triangolo
Triangle3dEx CurrentTriangle ;
CurrentTriangle.Set( VecField[i0].ptPApp, VecField[i1].ptPApp, VecField[i2].ptPApp) ;
CurrentTriangle.Validate( true) ;
CurrentTriangle.SetVertexNorm( 0, VecField[i0].vtVec) ;
CurrentTriangle.SetVertexNorm( 1, VecField[i1].vtVec) ;
CurrentTriangle.SetVertexNorm( 2, VecField[i2].vtVec) ;
// Setto il numero di utensile (conta solo positivo, nullo o negativo)
int nTool0 = Clamp( VecField[i0].nPropIndex, -1, 1) ;
int nTool1 = Clamp( VecField[i1].nPropIndex, -1, 1) ;
int nTool2 = Clamp( VecField[i2].nPropIndex, -1, 1) ;
if ( nTool0 == nTool1 || nTool0 == nTool2)
CurrentTriangle.SetGrade( nTool0) ;
else if ( nTool1 == nTool2)
CurrentTriangle.SetGrade( nTool1) ;
// Valido il triangolo e setto le normali del campo vettoriale ai corrispondenti vertici
if ( CurrentTriangle.Validate( true)) {
Vector3d vtVertNorm = VecField[i0].vtVec ;
if ( CurrentTriangle.GetN() * vtVertNorm > 0.6)
CurrentTriangle.SetVertexNorm( 0, vtVertNorm) ;
else
CurrentTriangle.SetVertexNorm( 0, CurrentTriangle.GetN()) ;
vtVertNorm = VecField[i1].vtVec ;
if ( CurrentTriangle.GetN() * vtVertNorm > 0.6)
CurrentTriangle.SetVertexNorm( 1, vtVertNorm) ;
else
CurrentTriangle.SetVertexNorm( 1, CurrentTriangle.GetN()) ;
vtVertNorm = VecField[i2].vtVec ;
if ( CurrentTriangle.GetN() * vtVertNorm > 0.6)
CurrentTriangle.SetVertexNorm( 2, vtVertNorm) ;
else
CurrentTriangle.SetVertexNorm( 2, CurrentTriangle.GetN()) ;
}
}
// Riporto il triangolo nel sistema in cui è immerso quello dello Zmap
CurrentTriangle.ToGlob( m_MapFrame) ;
// Aggiungo alla lista
lstTria.emplace_back( CurrentTriangle) ;
}
continue ;
}
}
// Configurazione 10
else if ( nAllConfig[nIndex] == 10) {
// Test sulla topologia
bool bDefTopology = false ;
bool bNewTopology = false ;
int nCount = 0 ;
while ( nIndexConfig10[nCount] != nIndex)
++ nCount ;
int nISlSt = i ;
int nJSlSt = j ;
int nKSlSt = k ;
int nISlEn = ( nAdjVox10[nCount] != 1 ? i : i + 1) ;
int nJSlEn = ( nAdjVox10[nCount] != 2 ? j : j + 1) ;
int nKSlEn = ( nAdjVox10[nCount] != 3 ? k : k + 1) ;
int nNSlSt, nNSlEn ;
if ( GetVoxNFromIJK( nISlSt, nJSlSt, nKSlSt, nNSlSt) &&
GetVoxNFromIJK( nISlEn, nJSlEn, nKSlEn, nNSlEn)) {
if ( abs( nAdjVox10[nCount]) == 1) {
auto itSt = SliceYZ.find( nNSlSt) ;
if ( itSt != SliceYZ.end()) {
bNewTopology = itSt->second ;
bDefTopology = true ;
}
}
else if ( abs( nAdjVox10[nCount]) == 2) {
auto itSt = SliceXZ.find( nNSlSt) ;
if ( itSt != SliceXZ.end()) {
bNewTopology = itSt->second ;
bDefTopology = true ;
}
}
else if ( abs( nAdjVox10[nCount]) == 3) {
auto itSt = SliceXY.find( nNSlSt) ;
if ( itSt != SliceXY.end()) {
bNewTopology = itSt->second ;
bDefTopology = true ;
}
}
}
if ( ! bDefTopology && bReg) {
// Verifico concordanza tra i versori di una stessa componente
// (ogni coppia di vettori di una medesima componente deve avere prodotto scalare non inferiore a 0.0)
Vector3d vtCmpAvg0, vtCmpAvg1 ;
bool bTest0 = DotTest( CompoVert[0], 4, vtCmpAvg0, 0.0) ;
bool bTest1 = DotTest( CompoVert[1], 4, vtCmpAvg1, 0.0) ;
bNewTopology = ( ! bTest0 || ! bTest1) ;
if ( GetVoxNFromIJK( nISlSt, nJSlSt, nKSlSt, nNSlSt) &&
GetVoxNFromIJK( nISlEn, nJSlEn, nKSlEn, nNSlEn)) {
if ( abs( nAdjVox6[nCount]) == 1) {
SliceYZ.emplace( nNSlSt, bNewTopology) ;
SliceYZ.emplace( nNSlEn, bNewTopology) ;
}
else if ( abs( nAdjVox6[nCount]) == 2) {
SliceXZ.emplace( nNSlSt, bNewTopology) ;
SliceXZ.emplace( nNSlEn, bNewTopology) ;
}
else if ( abs( nAdjVox6[nCount]) == 3) {
SliceXY.emplace( nNSlSt, bNewTopology) ;
SliceXY.emplace( nNSlEn, bNewTopology) ;
}
}
}
// Si passa alla seconda topologia
if ( bNewTopology) {
// Ricerca del caso corrispondente della nuova topologia
int nt = 0 ;
while ( nIndexVsIndex10[nt][0] != nIndex)
++ nt ;
// Riaggiorno gli offsets
nExtTabOff = 2 ;
nStdTabOff = 0 ;
// Modifico le matrici
int nRotCase = nIndexVsIndex10[nt][1] ;
for ( int nC = 1 ; nC <= 2 ; ++ nC) {
// Numero vertici per componenti
nVertComp[nC - 1] = Cases10Plus[nRotCase][1][nC] ;
// Matrice dei vertici della base del fan
for ( int nFanVert = 0 ; nFanVert < 4 ; ++ nFanVert)
CompoVert[nC - 1][nFanVert] = VecField[Cases10Plus[nRotCase][1][nFanVert + nExtTabOff + 1]] ;
// Matrici dei vertici dei triangoli in assenza di sharp feature
for ( int nTriVert = 0 ; nTriVert < 6 ; nTriVert += 3) {
CompoTriVert[nC - 1][nTriVert] = VecField[Cases10Plus[nRotCase][0][nStdTabOff + nTriVert+2]] ;
CompoTriVert[nC - 1][nTriVert+1] = VecField[Cases10Plus[nRotCase][0][nStdTabOff + nTriVert+1]] ;
CompoTriVert[nC - 1][nTriVert+2] = VecField[Cases10Plus[nRotCase][0][nStdTabOff + nTriVert]] ;
}
// Aggiorno gli offsets per raggiungere i vertici della componente successiva.
nExtTabOff += nVertComp[nC - 1] ;
nStdTabOff += 3 * ( nVertComp[nC - 1] - 2) ;
}
}
}
Voxel VoxConf ;
VoxConf.nNumComp = 0 ;