EgtGeomKernel 1.8h1 :

- inserite ultime migliorie per grafica Zmap
- corretta AvoidBox di Zmap.
This commit is contained in:
Dario Sassi
2017-08-10 17:16:30 +00:00
parent 99bbdda1d0
commit c17e44a3e1
4 changed files with 376 additions and 238 deletions
BIN
View File
Binary file not shown.
+1 -1
View File
@@ -480,7 +480,7 @@ VolZmap::AvoidBox( const Frame3d& frBox, const Vector3d& vtDiag)
for ( int nIndex = 0 ; nIndex < nSize ; nIndex += 1) {
if ( ! ( dZmax < m_Values[0][nPos][nIndex].dMin - EPS_SMALL ||
dZmin > m_Values[0][nPos][nIndex + 1].dMax + EPS_SMALL))
dZmin > m_Values[0][nPos][nIndex].dMax + EPS_SMALL))
return false ;
}
}
+369 -232
View File
@@ -713,7 +713,7 @@ VolZmap::ExtMarchingCubes( int nBlock, TRIA3DLIST& lstTria, TriHolder& triHold)
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) {
// Riconoscimento dei voxel di frontiera
int nVoxIndexes[3] = { i, j, k} ;
bool bBoundary = IsAVoxelOnBoundary( nLimits, nVoxIndexes, true) ;
@@ -851,169 +851,173 @@ VolZmap::ExtMarchingCubes( int nBlock, TRIA3DLIST& lstTria, TriHolder& triHold)
nStdTabOff += 3 * ( nVertComp[nCompCount - 1] - 2) ;
}
// Test sulla topologia
if ( nAllConfig[nIndex] == 3) {
// Test sulla topologia: dal momento che il nostro test
// si fonda sugli angoli compresi fra le normali, esso ha
// senso solo se il campo è regolare.
if ( bReg) {
if ( nAllConfig[nIndex] == 3) {
Vector3d vtCmpAvg0, vtCmpAvg1 ;
Vector3d vtCmpAvg0, vtCmpAvg1 ;
// Verifico se i versori delle componenti sono tutti
// più o meno concordi (per esserlo non devono esserci
// due vettori di una medesima componente con prodotto
// scalare inferiore a 0.7).
bool bTest0 = DotTest( CompoVert[0], 3, vtCmpAvg0, 0.7) ;
bool bTest1 = DotTest( CompoVert[1], 3, vtCmpAvg1, 0.7) ;
// Verifico se i versori delle componenti sono tutti
// più o meno concordi (per esserlo non devono esserci
// due vettori di una medesima componente con prodotto
// scalare inferiore a 0.7).
bool bTest0 = DotTest( CompoVert[0], 3, vtCmpAvg0, 0.7) ;
bool bTest1 = DotTest( CompoVert[1], 3, vtCmpAvg1, 0.7) ;
// Se i versori di entrambe le componenti sono concordi
// ha senso parlare di vettori medi, altrimenti non ha
// senso. Se non ha senso parlare di vettori medi non
// ha senso parlare di prodotti scalari fra loro,
// quindi pongo il loro prodotto a un valore assurdo -2
// (il prodotto scalare fra versori ha modulo non superiore
// a uno).
double dScProd = - 2 ;
// Se i versori di entrambe le componenti sono concordi
// ha senso parlare di vettori medi, altrimenti non ha
// senso. Se non ha senso parlare di vettori medi non
// ha senso parlare di prodotti scalari fra loro,
// quindi pongo il loro prodotto a un valore assurdo -2
// (il prodotto scalare fra versori ha modulo non superiore
// a uno).
double dScProd = - 2 ;
if ( bTest0 && bTest1)
dScProd = vtCmpAvg0 * vtCmpAvg1 ;
if ( bTest0 && bTest1)
dScProd = vtCmpAvg0 * vtCmpAvg1 ;
double dThreshold = 0.7 ;
double dThreshold = 0.7 ;
if ( ( ! ( bTest0 && bTest1)) || ( bTest0 && bTest1 && dScProd > dThreshold)) {
if ( ( ! ( bTest0 && bTest1)) || ( bTest0 && bTest1 && dScProd > dThreshold)) {
int nt = 0 ;
int nt = 0 ;
while ( nIndexVsIndex3[nt][0] != nIndex)
++ nt ;
while ( nIndexVsIndex3[nt][0] != nIndex)
++ nt ;
int nRotCase = nIndexVsIndex3[nt][1] ;
int nRotCase = nIndexVsIndex3[nt][1] ;
nComponents = Cases3Plus[nRotCase][1][0] ;
nComponents = Cases3Plus[nRotCase][1][0] ;
// Riaggiorno gli offsets
nExtTabOff = nComponents ;
nStdTabOff = 0 ;
// Riaggiorno gli offsets
nExtTabOff = nComponents ;
nStdTabOff = 0 ;
// Modifico le matrici
for ( int nC = 1 ; nC <= nComponents ; ++ nC) {
// Modifico le matrici
for ( int nC = 1 ; nC <= nComponents ; ++ nC) {
// Numero vertici per componenti
nVertComp[nC - 1] = Cases3Plus[nRotCase][1][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)
// 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]] ;
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) {
// 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]] ;
}
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) ;
}
// Aggiorno gli offsets per raggiungere i
// vertici della componente successiva.
nExtTabOff += nVertComp[nC - 1] ;
nStdTabOff += 3 * ( nVertComp[nC - 1] - 2) ;
}
}
}
else if ( nAllConfig[nIndex] == 6) {
// Procedura analoga a quella della configurazione 3
Vector3d vtCmpAvg0, vtCmpAvg1 ;
bool bTest0 = DotTest( CompoVert[0], 4, vtCmpAvg0, 0.7) ;
bool bTest1 = DotTest( CompoVert[1], 3, vtCmpAvg1, 0.7) ;
double dScProd = - 2 ;
if ( bTest0 && bTest1)
dScProd = vtCmpAvg0 * vtCmpAvg1 ;
double dThreshold = 0.7 ;
if ( ( ! ( bTest0 && bTest1)) || ( bTest0 && bTest1 && dScProd > dThreshold)) {
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) {
// Costruzione triangolo
int i0 = Cases6Plus[nRotCase][TriIndex + 2] ;
int i1 = Cases6Plus[nRotCase][TriIndex + 1] ;
int i2 = Cases6Plus[nRotCase][TriIndex] ;
Triangle3d CurrentTriangle ;
// Il triangolo è pronto
CurrentTriangle.Set( VecField[i0].ptInt, VecField[i1].ptInt, VecField[i2].ptInt) ;
bool bV = CurrentTriangle.Validate( true) ;
// Aggiungo alla lista
lstTria.emplace_back( CurrentTriangle) ;
}
continue ;
}
}
}
else if ( nAllConfig[nIndex] == 10) {
else if ( nAllConfig[nIndex] == 6) {
// Procedura analoga a quella della configurazione 3
Vector3d vtCmpAvg0, vtCmpAvg1 ;
Vector3d vtCmpAvg0, vtCmpAvg1 ;
bool bTest0 = DotTest( CompoVert[0], 4, vtCmpAvg0, 0.7) ;
bool bTest1 = DotTest( CompoVert[1], 3, vtCmpAvg1, 0.7) ;
double dScProd = - 2 ;
// Verifico se i versori delle componenti sono tutti
// più o meno concordi (per esserlo non devono esserci
// due vettori di una medesima componente con prodotto
// scalare inferiore a 0). decidere se 0.0 o 0.7
bool bTest0 = DotTest( CompoVert[0], 4, vtCmpAvg0) ;
bool bTest1 = DotTest( CompoVert[1], 4, vtCmpAvg1) ;
if ( bTest0 && bTest1)
dScProd = vtCmpAvg0 * vtCmpAvg1 ;
if ( ! ( bTest0 && bTest1)) {
double dThreshold = 0.7 ;
int nt = 0 ;
if ( ( ! ( bTest0 && bTest1)) || ( bTest0 && bTest1 && dScProd > dThreshold)) {
int nt = 0 ;
while ( nIndexVsIndex6[nt][0] != nIndex)
++ nt ;
int nRotCase = nIndexVsIndex6[nt][1] ;
while ( nIndexVsIndex10[nt][0] != nIndex)
++ nt ;
int nRotCase = nIndexVsIndex10[nt][1] ;
// Costruzione dei triangoli
for ( int TriIndex = 0 ; TriIndex < 15 ; TriIndex += 3) {
// Riaggiorno gli offsets
nExtTabOff = 2 ;
nStdTabOff = 0 ;
// Costruzione triangolo
int i0 = Cases6Plus[nRotCase][TriIndex + 2] ;
int i1 = Cases6Plus[nRotCase][TriIndex + 1] ;
int i2 = Cases6Plus[nRotCase][TriIndex] ;
// Modifico le matrici
for ( int nC = 1 ; nC <= 2 ; ++ nC) {
// Numero vertici per componenti
nVertComp[nC - 1] = Cases10Plus[nRotCase][1][nC] ;
Triangle3d CurrentTriangle ;
// Matrice dei vertici della base del fan
for ( int nFanVert = 0 ; nFanVert < 4 ; ++ nFanVert)
// Il triangolo è pronto
CurrentTriangle.Set( VecField[i0].ptInt, VecField[i1].ptInt, VecField[i2].ptInt) ;
bool bV = CurrentTriangle.Validate( true) ;
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]] ;
// Aggiungo alla lista
lstTria.emplace_back( CurrentTriangle) ;
}
continue ;
}
}
else if ( nAllConfig[nIndex] == 10) {
Vector3d vtCmpAvg0, vtCmpAvg1 ;
// Verifico se i versori delle componenti sono tutti
// più o meno concordi (per esserlo non devono esserci
// due vettori di una medesima componente con prodotto
// scalare inferiore a 0). decidere se 0.0 o 0.7
bool bTest0 = DotTest( CompoVert[0], 4, vtCmpAvg0) ;
bool bTest1 = DotTest( CompoVert[1], 4, vtCmpAvg1) ;
if ( ! ( bTest0 && bTest1)) {
int nt = 0 ;
while ( nIndexVsIndex10[nt][0] != nIndex)
++ nt ;
// Aggiorno gli offsets per raggiungere i
// vertici della componente successiva.
nExtTabOff += nVertComp[nC - 1] ;
nStdTabOff += 3 * ( nVertComp[nC - 1] - 2) ;
}
}
int nRotCase = nIndexVsIndex10[nt][1] ;
// Riaggiorno gli offsets
nExtTabOff = 2 ;
nStdTabOff = 0 ;
// Modifico le matrici
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) ;
}
}
}
}
@@ -1288,7 +1292,7 @@ VolZmap::ExtMarchingCubes( int nBlock, TRIA3DLIST& lstTria, TriHolder& triHold)
int nPosD ;
Vector3d vtD, vtE ;
if ( nVertComp[nCompCount - 1] == 4 && nFeatureType == 2) {
if ( nVertComp[nCompCount - 1] == 4 && nFeatureType == EDGE) {
int nPosEq ;
for ( int ni = 0 ; ni < 2 ; ++ ni) {
for ( int nj = ni + 1 ; nj < nVertComp[nCompCount - 1] ; ++ nj) {
@@ -1410,17 +1414,16 @@ VolZmap::ExtMarchingCubes( int nBlock, TRIA3DLIST& lstTria, TriHolder& triHold)
dUnknownVector( 2)) ;
double dDistFeature = vtFeature.Len() ;
const double dMaxDist = sqrt( 3) * m_dStep ;
const double MAX_DIST = sqrt( 3) * m_dStep ;
// Limito la feature a una distanza dal
// baricentro pari alla diagonale del voxel.
if ( dDistFeature > dMaxDist)
if ( dDistFeature > MAX_DIST)
bOutside = true ;
Triangle3d CurrentTriangle ;
TRIA3DVECTOR triContainer ;
// Costruisco triangoli di prova
for ( int ni = 0 ; ni < nVertComp[nCompCount - 1] ; ++ ni) {
@@ -1432,106 +1435,242 @@ VolZmap::ExtMarchingCubes( int nBlock, TRIA3DLIST& lstTria, TriHolder& triHold)
// Aggiungo triangolo al vettore temporaneo
triContainer.emplace_back( CurrentTriangle) ;
}
// Controllo delle inversioni dei triangoli
bool bInversione = false ;
// Controllo delle inversioni dei triangoli
bool bDangerInversion = false ;
// Caso ventaglio con tre vertici di base
if ( nVertComp[nCompCount - 1] == 3) {
// Controllo se esiste almeno un triangolo con normale avente prodotto scalare
// negativo con la normale di almeno uno dei vertici di base del ventaglio.
bool bInversione = false ;
for ( int ni = 0 ; ni < nVertComp[nCompCount - 1] ; ++ ni) {
for ( int ni = 0 ; ni < nVertComp[nCompCount - 1] ; ++ ni) {
int nj = ( ni + 1 < nVertComp[nCompCount - 1]) ? ni + 1 : 0 ;
int nj = ( ni + 1 < nVertComp[nCompCount - 1]) ? ni + 1 : 0 ;
if ( triContainer[ni].GetN() * CompoVert[nCompCount - 1][ni].vtNorm < - EPS_SMALL ||
triContainer[ni].GetN() * CompoVert[nCompCount - 1][nj].vtNorm < - EPS_SMALL)
double dDI = triContainer[ni].GetN() * CompoVert[nCompCount - 1][ni].vtNorm ;
double dDJ = triContainer[ni].GetN() * CompoVert[nCompCount - 1][nj].vtNorm ;
if ( dDI < - EPS_SMALL || dDJ < - EPS_SMALL)
bInversione = true ;
}
if ( bInversione &&
( nVertComp[nCompCount - 1] == 3 ||
nVertComp[nCompCount - 1] == 4)) {
int nNegDot = 0 ;
for ( int ni = 0 ; ni < nVertComp[nCompCount - 1] - 1 ; ++ ni) {
for ( int nj = ni + 1 ; nj < nVertComp[nCompCount - 1] ; ++ nj) {
if ( CompoVert[nCompCount - 1][ni].vtNorm * CompoVert[nCompCount - 1][nj].vtNorm < - EPS_SMALL)
nNegDot ++ ;
}
}
if ( nNegDot == nVertComp[nCompCount - 1] - 1) {
// Se tale triangolo esiste procedo
if ( bInversione) {
Point3d ptSolZMapFrame = ptSol ;
ptSolZMapFrame.ToLoc( m_MapFrame[0]) ;
if ( ! IsPointInsideVoxelApprox( i, j, k, ptSolZMapFrame)) {
Vector3d vtNullSpace( dMatrixV( 0, 2), dMatrixV( 1, 2), dMatrixV( 2, 2)) ;
vtNullSpace.ToLoc( m_MapFrame[0]) ;
double dParInt1, dParInt2 ;
Point3d ptVoxMin( ( i + 0.5) * m_dStep, ( j + 0.5) * m_dStep, ( k + 0.5) * m_dStep) ;
Point3d ptVoxMax( ( i + 1.5) * m_dStep, ( j + 1.5) * m_dStep, ( k + 1.5) * m_dStep) ;
if ( IntersLineBox( ptSolZMapFrame, vtNullSpace, ptVoxMin, ptVoxMax, dParInt1, dParInt2)) {
triContainer.resize( 0) ;
double dPar = abs( dParInt1) < abs( dParInt2) ? dParInt1 + ( dParInt2 - dParInt1) / 100 :
dParInt2 + ( dParInt1 - dParInt2) / 100 ;
Point3d ptNewSol = ptSolZMapFrame + dPar * vtNullSpace ;
ptNewSol.ToGlob( m_MapFrame[0]) ;
ptSol = ptNewSol ;
// Costruisco triangoli di prova
for ( int ni = 0 ; ni < nVertComp[nCompCount - 1] ; ++ ni) {
int nj = ( ni + 1 < nVertComp[nCompCount - 1]) ? ni + 1 : 0 ;
// Il triangolo è pronto
CurrentTriangle.Set( ptSol, CompoVert[nCompCount - 1][nj].ptInt, CompoVert[nCompCount - 1][ni].ptInt) ;
CurrentTriangle.Validate( true) ;
// Aggiungo triangolo al vettore temporaneo
triContainer.emplace_back( CurrentTriangle) ;
}
// Conto le coppie di normali con angolo compreso maggiore di 90°
int nNegDot = 0 ;
for ( int ni = 0 ; ni < nVertComp[nCompCount - 1] - 1 ; ++ ni) {
for ( int nj = ni + 1 ; nj < nVertComp[nCompCount - 1] ; ++ nj) {
if ( CompoVert[nCompCount - 1][ni].vtNorm * CompoVert[nCompCount - 1][nj].vtNorm < - EPS_SMALL)
nNegDot ++ ;
}
else {
}
if ( nNegDot == nVertComp[nCompCount - 1] - 1) {
// Cerco se esiste un punto in cui la normale ha prodotto scalre negativo
// con le normali di entrambi i triangoli che lo contengono
bool bInversione2 = false ;
for ( int ni = 0 ; ni < nVertComp[nCompCount - 1] ; ++ ni) {
int nj = ( ni == 0 ? nVertComp[nCompCount - 1] - 1 : ni - 1) ;
double dDLast = triContainer[ni].GetN() * CompoVert[nCompCount - 1][ni].vtNorm ;
double dDPrev = triContainer[nj].GetN() * CompoVert[nCompCount - 1][ni].vtNorm ;
if ( dDLast < EPS_SMALL && dDPrev < EPS_SMALL)
bInversione2 = true ;
}
// Se tale normale esiste
if ( bInversione2) {
// Soluzione del sistema nel sistema Zmap
Point3d ptSolZMapFrame = ptSol ;
ptSolZMapFrame.ToLoc( m_MapFrame[0]) ;
// Se la soluzione non cade nel voxel di appartenenza vedo se può
// essere riportata dentro muovendosi lungo la linea di feature.
if ( ! IsPointInsideVoxelApprox( i, j, k, ptSolZMapFrame)) {
Vector3d vtNullSpace( dMatrixV( 0, 2), dMatrixV( 1, 2), dMatrixV( 2, 2)) ;
vtNullSpace.ToLoc( m_MapFrame[0]) ;
double dParInt1, dParInt2 ;
Point3d ptVoxMin( ( i + 0.5) * m_dStep, ( j + 0.5) * m_dStep, ( k + 0.5) * m_dStep) ;
Point3d ptVoxMax( ( i + 1.5) * m_dStep, ( j + 1.5) * m_dStep, ( k + 1.5) * m_dStep) ;
// Caso in cui può essere riportata dentro: se il voxel in cui cade la feature è pieno
// la riporto muovendola lungo la sua linea
if ( IntersLineBox( ptSolZMapFrame, vtNullSpace, ptVoxMin, ptVoxMax, dParInt1, dParInt2)) {
triContainer.resize( 0) ;
double dPar = abs( dParInt1) < abs( dParInt2) ? dParInt1 + ( dParInt2 - dParInt1) / 100 :
dParInt2 + ( dParInt1 - dParInt2) / 100 ;
Point3d ptNewSol = ptSolZMapFrame + dPar * vtNullSpace ;
ptNewSol.ToGlob( m_MapFrame[0]) ;
ptSol = ptNewSol ;
// Costruisco triangoli di prova
for ( int ni = 0 ; ni < nVertComp[nCompCount - 1] ; ++ ni) {
int nj = ( ni + 1 < nVertComp[nCompCount - 1]) ? ni + 1 : 0 ;
// Il triangolo è pronto
CurrentTriangle.Set( ptSol, CompoVert[nCompCount - 1][nj].ptInt, CompoVert[nCompCount - 1][ni].ptInt) ;
CurrentTriangle.Validate( true) ;
// Aggiungo triangolo al vettore temporaneo
triContainer.emplace_back( CurrentTriangle) ;
}
}
// Caso in cui non può essere riportata dentro:
// si passa alla routine standard se il voxel
// in cui cade è pieno.
else {
int nAdjVoxI, nAdjVoxJ, nAdjVoxK ;
if ( GetPointVoxel( ptSolZMapFrame, nAdjVoxI, nAdjVoxJ, nAdjVoxK)) {
int nAdjVoxI, nAdjVoxJ, nAdjVoxK ;
if ( GetPointVoxel( ptSolZMapFrame, nAdjVoxI, nAdjVoxJ, nAdjVoxK)) {
// Classificazione del voxel adiacente
int nAdjIndex = 0 ;
if ( IsThereMat( nAdjVoxI, nAdjVoxJ, nAdjVoxK))
nAdjIndex |= ( 1 << 0) ;
if ( IsThereMat( nAdjVoxI + 1, nAdjVoxJ, nAdjVoxK))
nAdjIndex |= ( 1 << 1) ;
if ( IsThereMat( nAdjVoxI + 1, nAdjVoxJ + 1, nAdjVoxK))
nAdjIndex |= ( 1 << 2) ;
if ( IsThereMat( nAdjVoxI, nAdjVoxJ + 1, nAdjVoxK))
nAdjIndex |= ( 1 << 3) ;
if ( IsThereMat( nAdjVoxI, nAdjVoxJ, nAdjVoxK + 1))
nAdjIndex |= ( 1 << 4) ;
if ( IsThereMat( nAdjVoxI + 1, nAdjVoxJ, nAdjVoxK + 1))
nAdjIndex |= ( 1 << 5) ;
if ( IsThereMat( nAdjVoxI + 1, nAdjVoxJ + 1, nAdjVoxK + 1))
nAdjIndex |= ( 1 << 6) ;
if ( IsThereMat( nAdjVoxI, nAdjVoxJ + 1, nAdjVoxK + 1))
nAdjIndex |= ( 1 << 7) ;
// Classificazione del voxel adiacente
int nAdjIndex = 0 ;
if ( IsThereMat( nAdjVoxI, nAdjVoxJ, nAdjVoxK))
nAdjIndex |= ( 1 << 0) ;
if ( IsThereMat( nAdjVoxI + 1, nAdjVoxJ, nAdjVoxK))
nAdjIndex |= ( 1 << 1) ;
if ( IsThereMat( nAdjVoxI + 1, nAdjVoxJ + 1, nAdjVoxK))
nAdjIndex |= ( 1 << 2) ;
if ( IsThereMat( nAdjVoxI, nAdjVoxJ + 1, nAdjVoxK))
nAdjIndex |= ( 1 << 3) ;
if ( IsThereMat( nAdjVoxI, nAdjVoxJ, nAdjVoxK + 1))
nAdjIndex |= ( 1 << 4) ;
if ( IsThereMat( nAdjVoxI + 1, nAdjVoxJ, nAdjVoxK + 1))
nAdjIndex |= ( 1 << 5) ;
if ( IsThereMat( nAdjVoxI + 1, nAdjVoxJ + 1, nAdjVoxK + 1))
nAdjIndex |= ( 1 << 6) ;
if ( IsThereMat( nAdjVoxI, nAdjVoxJ + 1, nAdjVoxK + 1))
nAdjIndex |= ( 1 << 7) ;
// Se il voxel è pieno
if ( EdgeTable[nAdjIndex] != 0)
bDangerInversion = true ;
}
}
}
}
else {
// Se il voxel è pieno
if ( EdgeTable[nAdjIndex] != 0)
bDangerInversion = true ;
}
}
}
}
}
}
}
// Ventaglio con base a quattro vertici
else if ( nVertComp[nCompCount - 1] == 4) {
// Controllo se esiste almeno un triangolo con normale avente prodotto scalare
// negativo con la normale di almeno uno dei vertici di base del ventaglio.
bool bInversione = false ;
for ( int ni = 0 ; ni < nVertComp[nCompCount - 1] ; ++ ni) {
if ( nVertComp[nCompCount - 1] == 4) {
int nj = ( ni + 1 < nVertComp[nCompCount - 1]) ? ni + 1 : 0 ;
double dDI = triContainer[ni].GetN() * CompoVert[nCompCount - 1][ni].vtNorm ;
double dDJ = triContainer[ni].GetN() * CompoVert[nCompCount - 1][nj].vtNorm ;
if ( dDI < - EPS_SMALL || dDJ < - EPS_SMALL)
bInversione = true ;
}
// Se tale triangolo esiste continuo i test
if ( bInversione) {
// Conto il numero di coppie di normali con prodotto scalare negativo
int nNegDot = 0 ;
for ( int ni = 0 ; ni < nVertComp[nCompCount - 1] - 1 ; ++ ni) {
for ( int nj = ni + 1 ; nj < nVertComp[nCompCount - 1] ; ++ nj) {
double dDot = CompoVert[nCompCount - 1][ni].vtNorm * CompoVert[nCompCount - 1][nj].vtNorm ;
if ( dDot < - EPS_SMALL)
nNegDot ++ ;
}
}
// Caso in cui tale numero è 3
if ( nNegDot == nVertComp[nCompCount - 1] - 1) {
Point3d ptSolZMapFrame = ptSol ;
ptSolZMapFrame.ToLoc( m_MapFrame[0]) ;
// Se la feature non cade nel suo voxel
if ( ! IsPointInsideVoxelApprox( i, j, k, ptSolZMapFrame)) {
Vector3d vtNullSpace( dMatrixV( 0, 2), dMatrixV( 1, 2), dMatrixV( 2, 2)) ;
vtNullSpace.ToLoc( m_MapFrame[0]) ;
double dParInt1, dParInt2 ;
Point3d ptVoxMin( ( i + 0.5) * m_dStep, ( j + 0.5) * m_dStep, ( k + 0.5) * m_dStep) ;
Point3d ptVoxMax( ( i + 1.5) * m_dStep, ( j + 1.5) * m_dStep, ( k + 1.5) * m_dStep) ;
// Se è possibile riportarla dentro e il voxel in cui cade è pieno, la riporto nel suo voxel
// lungo la sua linea
if ( IntersLineBox( ptSolZMapFrame, vtNullSpace, ptVoxMin, ptVoxMax, dParInt1, dParInt2)) {
triContainer.resize( 0) ;
double dPar = abs( dParInt1) < abs( dParInt2) ? dParInt1 + ( dParInt2 - dParInt1) / 100 :
dParInt2 + ( dParInt1 - dParInt2) / 100 ;
Point3d ptNewSol = ptSolZMapFrame + dPar * vtNullSpace ;
ptNewSol.ToGlob( m_MapFrame[0]) ;
ptSol = ptNewSol ;
// Costruisco triangoli di prova
for ( int ni = 0 ; ni < nVertComp[nCompCount - 1] ; ++ ni) {
int nj = ( ni + 1 < nVertComp[nCompCount - 1]) ? ni + 1 : 0 ;
// Il triangolo è pronto
CurrentTriangle.Set( ptSol, CompoVert[nCompCount - 1][nj].ptInt, CompoVert[nCompCount - 1][ni].ptInt) ;
CurrentTriangle.Validate( true) ;
// Aggiungo triangolo al vettore temporaneo
triContainer.emplace_back( CurrentTriangle) ;
}
}
// Se non è possibile riportarla dentro e il voxel in
// cui cade è pieno passo alla routine standard
else {
int nAdjVoxI, nAdjVoxJ, nAdjVoxK ;
if ( GetPointVoxel( ptSolZMapFrame, nAdjVoxI, nAdjVoxJ, nAdjVoxK)) {
// Classificazione del voxel adiacente
int nAdjIndex = 0 ;
if ( IsThereMat( nAdjVoxI, nAdjVoxJ, nAdjVoxK))
nAdjIndex |= ( 1 << 0) ;
if ( IsThereMat( nAdjVoxI + 1, nAdjVoxJ, nAdjVoxK))
nAdjIndex |= ( 1 << 1) ;
if ( IsThereMat( nAdjVoxI + 1, nAdjVoxJ + 1, nAdjVoxK))
nAdjIndex |= ( 1 << 2) ;
if ( IsThereMat( nAdjVoxI, nAdjVoxJ + 1, nAdjVoxK))
nAdjIndex |= ( 1 << 3) ;
if ( IsThereMat( nAdjVoxI, nAdjVoxJ, nAdjVoxK + 1))
nAdjIndex |= ( 1 << 4) ;
if ( IsThereMat( nAdjVoxI + 1, nAdjVoxJ, nAdjVoxK + 1))
nAdjIndex |= ( 1 << 5) ;
if ( IsThereMat( nAdjVoxI + 1, nAdjVoxJ + 1, nAdjVoxK + 1))
nAdjIndex |= ( 1 << 6) ;
if ( IsThereMat( nAdjVoxI, nAdjVoxJ + 1, nAdjVoxK + 1))
nAdjIndex |= ( 1 << 7) ;
// Se il voxel è pieno
if ( EdgeTable[nAdjIndex] != 0)
bDangerInversion = true ;
}
}
}
}
// Caso in cui il numero di coppie di normali con prodotto
// scalare negativo non è 3
else {
Point3d ptSolZMapFrame = ptSol ;
ptSolZMapFrame.ToLoc( m_MapFrame[0]) ;
@@ -1571,10 +1710,10 @@ VolZmap::ExtMarchingCubes( int nBlock, TRIA3DLIST& lstTria, TriHolder& triHold)
}
}
}
}
}
}
}
}
// Valuto normali: questo è ancora un controllo
// sulle normali, se risultano in tutti i punti
// approssimativamente uguali passiamo alla
@@ -1598,7 +1737,7 @@ VolZmap::ExtMarchingCubes( int nBlock, TRIA3DLIST& lstTria, TriHolder& triHold)
if ( ! bPlane)
break ;
}
// Se la feature non è fuori dai limiti
// e i triangoli formano una superficie
// non piana confermo ExtMC
@@ -1717,7 +1856,7 @@ VolZmap::ExtMarchingCubes( int nBlock, TRIA3DLIST& lstTria, TriHolder& triHold)
}
// ExtMC non confermato, si passa a MC
else {
// Costruzione dei triangoli
for ( int TriIndex = 0; TriIndex < ( nVertComp[nCompCount - 1] - 2) * 3 ; TriIndex += 3) {
@@ -1736,7 +1875,7 @@ VolZmap::ExtMarchingCubes( int nBlock, TRIA3DLIST& lstTria, TriHolder& triHold)
}
// Standard MC
else {
// Costruzione dei triangoli
for ( int TriIndex = 0; TriIndex < ( nVertComp[nCompCount - 1] - 2) * 3 ; TriIndex += 3) {
@@ -2796,5 +2935,3 @@ VolZmap::IsATriangleOnBorder( const Triangle3d& trTria, const Point3d& ptVert,
else
return false ;
}
+6 -5
View File
@@ -3272,7 +3272,7 @@ VolZmap::CompCyl_ZDrilling( unsigned int nGrid, const Point3d & ptS, const Point
// Parametri geometrici dell'utensile
double dSqRad = dRad * dRad ;
double dSqfeSqRad = dSqRad - 2 * dRad * EPS_SMALL ;
double dSafeSqRad = dSqRad - 2 * dRad * EPS_SMALL ;
// Punte del gambo
Point3d ptTStemS = ptS - vtToolDir * dHei ;
@@ -3293,7 +3293,7 @@ VolZmap::CompCyl_ZDrilling( unsigned int nGrid, const Point3d & ptS, const Point
double dSqLen = vtC.SqLen() ;
// Se il punto si trova dentro il cerchio taglio
if ( dSqLen < dSqfeSqRad)
if ( dSqLen < dSafeSqRad)
SubtractIntervals( nGrid, i, j, dMinStemZ, dMaxStemZ, Z_AX, - Z_AX) ;
}
}
@@ -3318,8 +3318,9 @@ VolZmap::CompConus_ZDrilling( unsigned int nGrid, const Point3d & ptS, const Poi
double dAngC = dHei / ( dMaxRad - dMinRad) ;
double dSqMinRad = dMinRad * dMinRad ;
double dSqMaxRad = dMaxRad * dMaxRad ;
double dDeltaR = dMaxRad - dMinRad ;
double dSafeSqMaxRad = dSqMaxRad - 2 * dMaxRad * EPS_SMALL ; // Questa variabile è sperimentale: serve per evitare il taglio di un dexel dalla parte cilindrica del volume spazzato dalla traslazione del cono.
double dDeltaR = dMaxRad - dMinRad ; // Per tornare alla versione precedente basta sostituire dSafeSqMaxRad con dSqMaxRad. Per risolvere il problema in modo forse più sicuro, ma
// computazionalmente più pesante è sottrarre prima il cilindro con dSafeSqMaxRad e dopo il cono con dSqMaxRad.
// Studio delle simmetrie
if ( vtToolDir.z > 0) {
@@ -3354,7 +3355,7 @@ VolZmap::CompConus_ZDrilling( unsigned int nGrid, const Point3d & ptS, const Poi
SubtractIntervals( nGrid, i, j, dZMin, dZMax, Z_AX, -Z_AX) ;
else if ( dSqDist < dSqMaxRad) {
else if ( dSqDist < dSafeSqMaxRad) { // dSafeSqMaxRad è sperimentale
double dr = sqrt( dSqDist) ;