EgtGeomKernel 1.8j6 :

- correzioni e migliorie a Triangle3d
- migliorie a Zmap.
This commit is contained in:
Dario Sassi
2017-11-03 12:29:26 +00:00
parent 5396990de4
commit 7cf76d3175
5 changed files with 568 additions and 751 deletions
+308 -122
View File
@@ -24,6 +24,7 @@
#include "/EgtDev/Include/EGkStringUtils3d.h"
#include "/EgtDev/Extern/Eigen\Core"
#include "/EgtDev/Extern/Eigen\SVD"
#include "PolygonPlane.h"
using namespace std ;
@@ -1665,7 +1666,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) {
// Classificazione dei vertici: interni o esterni al materiale
int nIndex = CalcIndex( i, j , k) ;
@@ -2597,6 +2598,224 @@ VolZmap::ExtMarchingCubes( int nBlock, TRIA3DLIST& lstTria, TriHolder& triHold)
}
}
// 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) {
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 ;
// 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)) ;
double dParInt1, dParInt2 ;
Point3d ptVoxMin( ( i * N_DEXVOXRATIO + 0.5) * m_dStep,
( j * N_DEXVOXRATIO + 0.5) * m_dStep, ( k * N_DEXVOXRATIO + 0.5) * m_dStep) ;
Point3d ptVoxMax( ( ( i + 1) * N_DEXVOXRATIO + 0.5) * m_dStep,
( ( j + 1) * N_DEXVOXRATIO + 0.5) * m_dStep, ( ( k + 1) * N_DEXVOXRATIO + 0.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 ;
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
Triangle3d CurrentTriangle ;
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 = CalcIndex( nAdjVoxI, nAdjVoxJ, nAdjVoxK) ;
// 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 ;
if ( ! IsPointInsideVoxelApprox( i, j, k, ptSolZMapFrame)) {
Vector3d vtNullSpace( dMatrixV( 0, 2), dMatrixV( 1, 2), dMatrixV( 2, 2)) ;
double dParInt1, dParInt2 ;
Point3d ptVoxMin( ( i * N_DEXVOXRATIO + 0.5) * m_dStep,
( j * N_DEXVOXRATIO + 0.5) * m_dStep, ( k * N_DEXVOXRATIO + 0.5) * m_dStep) ;
Point3d ptVoxMax( ( ( i + 1) * N_DEXVOXRATIO + 0.5) * m_dStep,
( ( j + 1) * N_DEXVOXRATIO + 0.5) * m_dStep, ( ( k + 1) * N_DEXVOXRATIO + 0.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 ;
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
Triangle3d CurrentTriangle ;
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) ;
}
}
else {
int nCouple = 0 ;
int nCoupleIndex[4] ;
// Valuto il numero di coppie di vettori
// quasi coincidenti e per ogni coppia salvo gli
// indici dei vettori
for ( int nNI = 0 ; nNI < 3 ; ++ nNI) {
for ( int nNJ = nNI + 1 ; nNJ < 4 ; ++ nNJ) {
if ( AreSameVectorApprox( CompoVert[nCompCount - 1][nNI].vtNorm,
CompoVert[nCompCount - 1][nNJ].vtNorm)) {
++ nCouple ;
if ( nCouple == 1) {
nCoupleIndex[0] = nNI ;
nCoupleIndex[1] = nNJ ;
}
else if ( nCouple == 2) {
nCoupleIndex[2] = nNI ;
nCoupleIndex[3] = nNJ ;
}
}
}
}
// caso due coppie
if ( nCouple == 2) {
// vedo se c'è un triangolo con normale invertita rispetto a quelle
// si entrambi i vertici, se esiste si passerà a std MC
for ( int ni = 0 ; ni < 4 ; ++ ni) {
int nj = ( ni == 3 ? 0 : ni + 1) ;
if ( triContainer[ni].GetN() * CompoVert[nCompCount - 1][ni].vtNorm < - 0.95 &&
triContainer[ni].GetN() * CompoVert[nCompCount - 1][nj].vtNorm < - 0.95) {
bDangerInversion = true ;
break ;
}
}
}
// caso una coppia
else if ( nCouple == 1) {
// cerco gli indici dei vettori non appartenenti alla coppia
for ( int ni = 0 ; ni < 4 ; ++ ni) {
if ( ni != nCoupleIndex[0] && ni != nCoupleIndex[1]) {
nCoupleIndex[2] = ni ;
break ;
}
}
for ( int ni = 0 ; ni < 4 ; ++ ni) {
if ( ni != nCoupleIndex[0] && ni != nCoupleIndex[1] && ni != nCoupleIndex[2]) {
nCoupleIndex[3] = ni ;
break ;
}
}
// Media dei vettori coppia
Vector3d vtAv01 = 0.5 * ( CompoVert[nCompCount - 1][nCoupleIndex[0]].vtNorm +
CompoVert[nCompCount - 1][nCoupleIndex[1]].vtNorm) ;
// vettore nello spazio genenrato dai due non appartenenti alla coppia
Vector3d vtAv23 = 0.5 * ( CompoVert[nCompCount - 1][nCoupleIndex[2]].vtNorm +
CompoVert[nCompCount - 1][nCoupleIndex[3]].vtNorm) ;
vtAv01.Normalize() ;
vtAv23.Normalize() ;
double dDAvAV = vtAv01 * vtAv23 ;
// se angolo grande si esegue std MC
if ( abs( vtAv01 * vtAv23) < EPS_SMALL) {
for ( int ni = 0 ; ni < 4 ; ++ ni) {
int nj = ni == 3 ? 0 : ni + 1 ;
if ( triContainer[ni].GetN() * CompoVert[nCompCount - 1][ni].vtNorm < - 0.95 &&
triContainer[ni].GetN() * CompoVert[nCompCount - 1][nj].vtNorm < - 0.95) {
bDangerInversion = true ;
break ;
}
}
}
else {
double dD23 = CompoVert[nCompCount - 1][nCoupleIndex[2]].vtNorm *
CompoVert[nCompCount - 1][nCoupleIndex[3]].vtNorm ;
if ( dD23 > 0.7 && dDAvAV < 0.7) {
for ( int ni = 0 ; ni < 4 ; ++ ni) {
int nj = ni == 3 ? 0 : ni + 1 ;
if ( triContainer[ni].GetN() * CompoVert[nCompCount - 1][ni].vtNorm < - 0.9 &&
triContainer[ni].GetN() * CompoVert[nCompCount - 1][nj].vtNorm < - 0.9) {
bDangerInversion = true ;
break ;
}
}
}
}
}
}
}
}
}
}
// Controllo sulle normali, se risultano in tutti i punti
// approssimativamente uguali passiamo alla routine standard
int nContSize = int( triContainer.size()) ;
@@ -2934,20 +3153,55 @@ VolZmap::FlipEdgesII( TriHolder& TriHold, bool bGraph) const
// Si deve operare la modifica dei triangoli
if ( SharedIndex.size() > 2) {
// verifico che i due lati adiacenti siano controversi
// Verifico che i due lati adiacenti siano controversi
int nProd = ( SharedIndex[2] - SharedIndex[0]) * ( SharedIndex[3] - SharedIndex[1]) ;
if ( nProd != 1 && nProd != - 2 && nProd != 4 && nProd != 0) {
TriHold[n1].vCompoTria[nCompo1][nTri1].SetP( SharedIndex[0],
// Assegno l'array dei punti di contorno
Point3d vPnt[4] ;
vPnt[0] = TriHold[n1].vCompoTria[nCompo1][nTri1].GetP( SharedIndex[0]) ;
vPnt[1] = TriHold[n1].ptCompoVert[nCompo1] ;
vPnt[2] = TriHold[n1].vCompoTria[nCompo1][nTri1].GetP( SharedIndex[2]) ;
vPnt[3] = TriHold[n2].ptCompoVert[nCompo2] ;
// Valuto se i triangoli giacciono su un piano
PolygonPlane Polygon ;
for ( int i = 0 ; i < 4 ; ++ i)
Polygon.AddPoint( vPnt[i]) ;
Plane3d plPlane ;
bool bOnPlane = Polygon.GetPlane( plPlane) ;
for ( int i = 0 ; i < 4 && bOnPlane ; ++ i)
bOnPlane = PointInPlaneApprox( vPnt[i], plPlane) ;
// Se sono su un piano controllo se avviene inversione
bool bInv = false ;
if ( bOnPlane) {
Triangle3d trT1 = TriHold[n1].vCompoTria[nCompo1][nTri1] ;
Triangle3d trT2 = TriHold[n2].vCompoTria[nCompo2][nTri2] ;
int nVert1, nVert2 ;
// Determino gli indici dei punti sharp-feature
for ( int nP = 0 ; nP < 3 ; ++ nP) {
if ( nP != SharedIndex[0] && nP != SharedIndex[2])
nVert1 = nP ;
if ( nP != SharedIndex[1] && nP != SharedIndex[3])
nVert2 = nP ;
}
trT1.SetP( SharedIndex[0], trT2.GetP( nVert2)) ;
trT2.SetP( SharedIndex[3], trT1.GetP( nVert1)) ;
trT1.Validate( true) ;
trT2.Validate( true) ;
bInv = ( trT1.GetN() * trT2.GetN() < 0) ;
}
// Se non vi è inversione eseguo il flipping
if ( ! bInv) {
TriHold[n1].vCompoTria[nCompo1][nTri1].SetP( SharedIndex[0],
TriHold[n2].ptCompoVert[nCompo2]) ;
TriHold[n2].vCompoTria[nCompo2][nTri2].SetP( SharedIndex[3],
TriHold[n2].vCompoTria[nCompo2][nTri2].SetP( SharedIndex[3],
TriHold[n1].ptCompoVert[nCompo1]) ;
TriHold[n1].vCompoTria[nCompo1][nTri1].Validate( true) ;
TriHold[n2].vCompoTria[nCompo2][nTri2].Validate( true) ;
TriHold[n1].vCompoTria[nCompo1][nTri1].Validate( true) ;
TriHold[n2].vCompoTria[nCompo2][nTri2].Validate( true) ;
bModified = true ;
break ;
bModified = true ;
break ;
}
}
}
}
@@ -2963,108 +3217,6 @@ VolZmap::FlipEdgesII( TriHolder& TriHold, bool bGraph) const
return true ;
}
//////// SPERIMENTALE ///////////////////
//bool
//VolZmap::FlipEdgesII( TriHolder& TriHold) const
//{
// // Numero di voxel in cui si presentano sharp feature
// int nVoxelNum = int( TriHold.size()) ;
//
// // Ciclo su tali voxel
// for ( int n1 = 0 ; n1 < nVoxelNum ; ++ n1) {
//
// for ( int n2 = n1 ; n2 < nVoxelNum ; ++ n2) {
//
// // Se i voxel sono adiacenti proseguo
// if ( abs( TriHold[n2].i - TriHold[n1].i) == 1 ||
// abs( TriHold[n2].j - TriHold[n1].j) == 1 ||
// abs( TriHold[n2].k - TriHold[n1].k) == 1) {
//
// // Numero delle componenti connesse nei due voxel
// int nNumCompo1 = int( TriHold[n1].ptCompoVert.size()) ;
// int nNumCompo2 = int( TriHold[n2].ptCompoVert.size()) ;
//
// int nCompo1 = 0 ;
//
// // Ciclo sulle componenti
// for ( ; nCompo1 < nNumCompo1 ; ++ nCompo1) {
//
// int nCompo2 = ( n1 == n2 ? nCompo1 + 1 : 0) ;
//
// for ( ; nCompo2 < nNumCompo2 ; ++ nCompo2) {
//
// // Numero di triangoli per le componenti connesse
// int nTriNum1 = int( TriHold[n1].vCompoTria[nCompo1].size()) ;
// int nTriNum2 = int( TriHold[n2].vCompoTria[nCompo2].size()) ;
//
// for ( int nTri1 = 0 ; nTri1 < nTriNum1 ; ++ nTri1) {
//
// bool bModified = false ;
//
// for ( int nTri2 = 0 ; nTri2 < nTriNum2 ; ++ nTri2) {
//
// INTVECTOR SharedIndex ;
//
// for ( int nVert1 = 0 ; nVert1 < 3 ; ++ nVert1) {
//
// for ( int nVert2 = 0 ; nVert2 < 3 ; ++ nVert2) {
//
// Point3d ptP1 = TriHold[n1].vCompoTria[nCompo1][nTri1].GetP( nVert1) ;
// Point3d ptP2 = TriHold[n2].vCompoTria[nCompo2][nTri2].GetP( nVert2) ;
//
// if ( AreSamePointEpsilon( ptP1, ptP2, EPS_ZERO)) {
//
// Point3d ptVert1 = TriHold[n1].ptCompoVert[nCompo1] ;
// Point3d ptVert2 = TriHold[n2].ptCompoVert[nCompo2] ;
//
// if ( ! ( AreSamePointEpsilon( ptP1, ptVert1, EPS_ZERO) ||
// AreSamePointEpsilon( ptP2, ptVert2, EPS_ZERO))) {
//
// SharedIndex.emplace_back( nVert1) ;
// SharedIndex.emplace_back( nVert2) ;
// }
// }
//
// if ( SharedIndex.size() > 2)
// break ;
// }
//
// if ( SharedIndex.size() > 2)
// break ;
// }
//
// // Si deve operare la modifica dei triangoli
// if ( SharedIndex.size() > 2) {
//
// // verifico che i due lati adiacenti siano controversi
// int nProd = ( SharedIndex[2] - SharedIndex[0]) * ( SharedIndex[3] - SharedIndex[1]) ;
// if ( nProd != 1 && nProd != - 2 && nProd != 4 && nProd != 0) {
//
// TriHold[n1].vCompoTria[nCompo1][nTri1].SetP( SharedIndex[0],
// TriHold[n2].ptCompoVert[nCompo2]) ;
// TriHold[n2].vCompoTria[nCompo2][nTri2].SetP( SharedIndex[3],
// TriHold[n1].ptCompoVert[nCompo1]) ;
//
// TriHold[n1].vCompoTria[nCompo1][nTri1].Validate( true) ;
// TriHold[n2].vCompoTria[nCompo2][nTri2].Validate( true) ;
//
// bModified = true ;
// break ;
// }
// }
// }
//
// if ( bModified)
// break ;
// }
// }
// }
// }
// }
// }
//
// return true ;
//}
//----------------------------------------------------------------------------
bool
@@ -3154,21 +3306,55 @@ VolZmap::FlipEdgesBB( TriaMatrix& InterTria) const
// Si deve operare la modifica dei triangoli
if ( SharedIndex.size() > 2) {
// verifico che i due lai adiacenti siano controversi
// Verifico che i due lati adiacenti siano controversi
int nProd = ( SharedIndex[2] - SharedIndex[0]) * ( SharedIndex[3] - SharedIndex[1]) ;
if ( nProd != 1 && nProd != - 2 && nProd != 4 && nProd != 0) {
// Assegno l'array dei punti di contorno
Point3d vPnt[4] ;
vPnt[0] = InterTria[tFB][tVFB].vCompoTria[tCmpF][tTriFB].GetP( SharedIndex[0]) ;
vPnt[1] = InterTria[tFB][tVFB].ptCompoVert[tCmpF] ;
vPnt[2] = InterTria[tFB][tVFB].vCompoTria[tCmpF][tTriFB].GetP( SharedIndex[2]) ;
vPnt[3] = InterTria[tLB][tVLB].ptCompoVert[tCmpL] ;
// Valuto se i triangoli giacciono su un piano
PolygonPlane Polygon ;
for ( int i = 0 ; i < 4 ; ++ i)
Polygon.AddPoint( vPnt[i]) ;
Plane3d plPlane ;
bool bOnPlane = Polygon.GetPlane( plPlane) ;
for ( int i = 0 ; i < 4 && bOnPlane ; ++ i)
bOnPlane = PointInPlaneApprox( vPnt[i], plPlane) ;
// Se sono su un piano controllo se avviene inversione
bool bInv = false ;
if ( bOnPlane) {
Triangle3d trTF = InterTria[tFB][tVFB].vCompoTria[tCmpF][tTriFB] ;
Triangle3d trTL = InterTria[tLB][tVLB].vCompoTria[tCmpL][tTriLB] ;
int nVertF, nVertL ;
// Determino gli indici dei punti sharp-feature
for ( int nP = 0 ; nP < 3 ; ++ nP) {
if ( nP != SharedIndex[0] && nP != SharedIndex[2])
nVertF = nP ;
if ( nP != SharedIndex[1] && nP != SharedIndex[3])
nVertL = nP ;
}
trTF.SetP( SharedIndex[0], trTL.GetP( nVertL)) ;
trTF.Validate( true) ;
trTL.SetP( SharedIndex[3], trTF.GetP( nVertF)) ;
trTL.Validate( true) ;
bInv = ( trTF.GetN() * trTL.GetN() < 0) ;
}
// Se non vi è inversione eseguo il flipping
if ( ! bInv) {
InterTria[tFB][tVFB].vCompoTria[tCmpF][tTriFB].SetP( SharedIndex[0],
InterTria[tLB][tVLB].ptCompoVert[tCmpL]) ;
InterTria[tLB][tVLB].vCompoTria[tCmpL][tTriLB].SetP( SharedIndex[3],
InterTria[tFB][tVFB].ptCompoVert[tCmpF]) ;
InterTria[tFB][tVFB].vCompoTria[tCmpF][tTriFB].SetP( SharedIndex[0],
InterTria[tLB][tVLB].ptCompoVert[tCmpL]) ;
InterTria[tLB][tVLB].vCompoTria[tCmpL][tTriLB].SetP( SharedIndex[3],
InterTria[tFB][tVFB].ptCompoVert[tCmpF]) ;
InterTria[tFB][tVFB].vCompoTria[tCmpF][tTriFB].Validate( true) ;
InterTria[tLB][tVLB].vCompoTria[tCmpL][tTriLB].Validate( true) ;
InterTria[tFB][tVFB].vCompoTria[tCmpF][tTriFB].Validate( true) ;
InterTria[tLB][tVLB].vCompoTria[tCmpL][tTriLB].Validate( true) ;
bModified = true ;
break ;
bModified = true ;
break ;
}
}
}
}