EgtGeomKernel 1.8h3 :

- migliorie e correzioni su Zmap
- aggiunta a PointGrid3d nuova FindNearest
- migliorie a TriMesh.
This commit is contained in:
Dario Sassi
2017-08-31 07:43:41 +00:00
parent 9e400bee15
commit c8a38f5aae
6 changed files with 251 additions and 196 deletions
+194 -188
View File
@@ -712,8 +712,8 @@ VolZmap::ExtMarchingCubes( int nBlock, TRIA3DLIST& lstTria, TriHolder& triHold)
// Ciclo su tutti i voxel dello Zmap
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) {
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) ;
@@ -824,22 +824,17 @@ VolZmap::ExtMarchingCubes( int nBlock, TRIA3DLIST& lstTria, TriHolder& triHold)
// Riempio il nCompCount-esimo vettore di vertici della base del fan
for ( int nVertCount = 0 ; nVertCount < nVertComp[nCompCount - 1] ; ++ nVertCount)
CompoVert[nCompCount - 1][nVertCount] = VecField[TriangleTableEn[nIndex][1][nVertCount + nExtTabOff + 1]] ;
// Serve per la gestione del caso ...
if ( nVertComp[nCompCount - 1] == 4) {
for ( int nVertCount = 0 ; nVertCount < nVertComp[nCompCount - 1] ; ++ nVertCount)
nIndArrey[nCompCount - 1][nVertCount] = TriangleTableEn[nIndex][1][nVertCount + nExtTabOff + 1] ;
}
// Riempio il nCompCount-esimo vettore di vertici dei triangoli in assenza di
// sharp feature: in una mesh di triangoli con n vertici vi sono n - 2 triangoli.
for ( int nVertCount = 0 ; nVertCount < 3 * ( nVertComp[nCompCount - 1] - 2) ; nVertCount += 3) {
CompoTriVert[nCompCount - 1][nVertCount] = VecField[TriangleTableEn[nIndex][0][nStdTabOff + nVertCount+2]] ;
CompoTriVert[nCompCount - 1][nVertCount+1] = VecField[TriangleTableEn[nIndex][0][nStdTabOff + nVertCount+1]] ;
CompoTriVert[nCompCount - 1][nVertCount+2] = VecField[TriangleTableEn[nIndex][0][nStdTabOff + nVertCount]] ;
@@ -875,12 +870,12 @@ VolZmap::ExtMarchingCubes( int nBlock, TRIA3DLIST& lstTria, TriHolder& triHold)
// a uno).
double dScProd = - 2 ;
if ( bTest0 && bTest1)
if ( bTest0 && bTest1)
dScProd = vtCmpAvg0 * vtCmpAvg1 ;
double dThreshold = 0.7 ;
if ( ( ! ( bTest0 && bTest1)) || ( bTest0 && bTest1 && dScProd > dThreshold)) {
if ( ! ( bTest0 && bTest1) || ( bTest0 && bTest1 && dScProd > dThreshold)) {
int nt = 0 ;
@@ -933,12 +928,12 @@ VolZmap::ExtMarchingCubes( int nBlock, TRIA3DLIST& lstTria, TriHolder& triHold)
double dScProd = - 2 ;
if ( bTest0 && bTest1)
if ( bTest0 && bTest1)
dScProd = vtCmpAvg0 * vtCmpAvg1 ;
double dThreshold = 0.7 ;
if ( ( ! ( bTest0 && bTest1)) || ( bTest0 && bTest1 && dScProd > dThreshold)) {
if ( ! ( bTest0 && bTest1) || ( bTest0 && bTest1 && dScProd > dThreshold)) {
int nt = 0 ;
@@ -1049,22 +1044,22 @@ VolZmap::ExtMarchingCubes( int nBlock, TRIA3DLIST& lstTria, TriHolder& triHold)
}
}
if ( ( nIndArrey[nCompCount - 1][0] == 0 && nIndArrey[nCompCount - 1][1] == 2 &&
nIndArrey[nCompCount - 1][2] == 9 && nIndArrey[nCompCount - 1][3] == 10) ||
( nIndArrey[nCompCount - 1][0] == 4 && nIndArrey[nCompCount - 1][1] == 6 &&
nIndArrey[nCompCount - 1][2] == 9 && nIndArrey[nCompCount - 1][3] == 10) ||
( nIndArrey[nCompCount - 1][0] == 4 && nIndArrey[nCompCount - 1][1] == 6 &&
nIndArrey[nCompCount - 1][2] == 8 && nIndArrey[nCompCount - 1][3] == 11) ||
( nIndArrey[nCompCount - 1][0] == 0 && nIndArrey[nCompCount - 1][1] == 2 &&
nIndArrey[nCompCount - 1][2] == 8 && nIndArrey[nCompCount - 1][3] == 11) ||
( nIndArrey[nCompCount - 1][0] == 1 && nIndArrey[nCompCount - 1][1] == 3 &&
nIndArrey[nCompCount - 1][2] == 8 && nIndArrey[nCompCount - 1][3] == 9 ) ||
( nIndArrey[nCompCount - 1][0] == 1 && nIndArrey[nCompCount - 1][1] == 3 &&
nIndArrey[nCompCount - 1][2] == 10 && nIndArrey[nCompCount - 1][3] == 11) ||
( nIndArrey[nCompCount - 1][0] == 5 && nIndArrey[nCompCount - 1][1] == 7 &&
nIndArrey[nCompCount - 1][2] == 10 && nIndArrey[nCompCount - 1][3] == 11) ||
( nIndArrey[nCompCount - 1][0] == 5 && nIndArrey[nCompCount - 1][1] == 7 &&
nIndArrey[nCompCount - 1][2] == 8 && nIndArrey[nCompCount - 1][3] == 9 )) {
if ( ( nIndArrey[nCompCount - 1][0] == 0 && nIndArrey[nCompCount - 1][1] == 2 &&
nIndArrey[nCompCount - 1][2] == 9 && nIndArrey[nCompCount - 1][3] == 10) ||
( nIndArrey[nCompCount - 1][0] == 4 && nIndArrey[nCompCount - 1][1] == 6 &&
nIndArrey[nCompCount - 1][2] == 9 && nIndArrey[nCompCount - 1][3] == 10) ||
( nIndArrey[nCompCount - 1][0] == 4 && nIndArrey[nCompCount - 1][1] == 6 &&
nIndArrey[nCompCount - 1][2] == 8 && nIndArrey[nCompCount - 1][3] == 11) ||
( nIndArrey[nCompCount - 1][0] == 0 && nIndArrey[nCompCount - 1][1] == 2 &&
nIndArrey[nCompCount - 1][2] == 8 && nIndArrey[nCompCount - 1][3] == 11) ||
( nIndArrey[nCompCount - 1][0] == 1 && nIndArrey[nCompCount - 1][1] == 3 &&
nIndArrey[nCompCount - 1][2] == 8 && nIndArrey[nCompCount - 1][3] == 9 ) ||
( nIndArrey[nCompCount - 1][0] == 1 && nIndArrey[nCompCount - 1][1] == 3 &&
nIndArrey[nCompCount - 1][2] == 10 && nIndArrey[nCompCount - 1][3] == 11) ||
( nIndArrey[nCompCount - 1][0] == 5 && nIndArrey[nCompCount - 1][1] == 7 &&
nIndArrey[nCompCount - 1][2] == 10 && nIndArrey[nCompCount - 1][3] == 11) ||
( nIndArrey[nCompCount - 1][0] == 5 && nIndArrey[nCompCount - 1][1] == 7 &&
nIndArrey[nCompCount - 1][2] == 8 && nIndArrey[nCompCount - 1][3] == 9 )) {
VectorField LocVecF[12], LocCompV[12] ;
@@ -1079,19 +1074,19 @@ VolZmap::ExtMarchingCubes( int nBlock, TRIA3DLIST& lstTria, TriHolder& triHold)
LocCompV[LocInd].vtNorm.ToLoc( m_MapFrame[0]) ;
}
if ( ( AreSameVectorApprox( LocVecF[nIndArrey[nCompCount - 1][0]].vtNorm, LocVecF[nIndArrey[nCompCount - 1][1]].vtNorm) &&
abs( LocVecF[nIndArrey[nCompCount - 1][0]].vtNorm *LocVecF[nIndArrey[nCompCount - 1][2]].vtNorm) < EPS_SMALL &&
abs( LocVecF[nIndArrey[nCompCount - 1][0]].vtNorm * LocVecF[nIndArrey[nCompCount - 1][3]].vtNorm) < EPS_SMALL) ||
( AreSameVectorApprox( LocVecF[nIndArrey[nCompCount - 1][2]].vtNorm, LocVecF[nIndArrey[nCompCount - 1][3]].vtNorm) &&
abs( LocVecF[nIndArrey[nCompCount - 1][2]].vtNorm * LocVecF[nIndArrey[nCompCount - 1][0]].vtNorm) < EPS_SMALL &&
if ( ( AreSameVectorApprox( LocVecF[nIndArrey[nCompCount - 1][0]].vtNorm, LocVecF[nIndArrey[nCompCount - 1][1]].vtNorm) &&
abs( LocVecF[nIndArrey[nCompCount - 1][0]].vtNorm *LocVecF[nIndArrey[nCompCount - 1][2]].vtNorm) < EPS_SMALL &&
abs( LocVecF[nIndArrey[nCompCount - 1][0]].vtNorm * LocVecF[nIndArrey[nCompCount - 1][3]].vtNorm) < EPS_SMALL) ||
( AreSameVectorApprox( LocVecF[nIndArrey[nCompCount - 1][2]].vtNorm, LocVecF[nIndArrey[nCompCount - 1][3]].vtNorm) &&
abs( LocVecF[nIndArrey[nCompCount - 1][2]].vtNorm * LocVecF[nIndArrey[nCompCount - 1][0]].vtNorm) < EPS_SMALL &&
abs( LocVecF[nIndArrey[nCompCount - 1][2]].vtNorm * LocVecF[nIndArrey[nCompCount - 1][1]].vtNorm) < EPS_SMALL)) {
Point3d ptBarycenter = ( LocCompV[0].ptInt + LocCompV[1].ptInt +
LocCompV[2].ptInt + LocCompV[3].ptInt) / 4 ;
if ( abs( LocCompV[0].ptInt.x - ptBarycenter.x) < EPS_SMALL &&
abs( LocCompV[1].ptInt.x - ptBarycenter.x) < EPS_SMALL &&
abs( LocCompV[2].ptInt.x - ptBarycenter.x) < EPS_SMALL &&
if ( abs( LocCompV[0].ptInt.x - ptBarycenter.x) < EPS_SMALL &&
abs( LocCompV[1].ptInt.x - ptBarycenter.x) < EPS_SMALL &&
abs( LocCompV[2].ptInt.x - ptBarycenter.x) < EPS_SMALL &&
abs( LocCompV[3].ptInt.x - ptBarycenter.x) < EPS_SMALL) {
double dXBar = ptBarycenter.x / m_dStep - 0.5 ;
@@ -1110,9 +1105,9 @@ VolZmap::ExtMarchingCubes( int nBlock, TRIA3DLIST& lstTria, TriHolder& triHold)
++ nBarInd ;
}
}
else if ( abs( LocCompV[0].ptInt.y - ptBarycenter.y) < EPS_SMALL &&
abs( LocCompV[1].ptInt.y - ptBarycenter.y) < EPS_SMALL &&
abs( LocCompV[2].ptInt.y - ptBarycenter.y) < EPS_SMALL &&
else if ( abs( LocCompV[0].ptInt.y - ptBarycenter.y) < EPS_SMALL &&
abs( LocCompV[1].ptInt.y - ptBarycenter.y) < EPS_SMALL &&
abs( LocCompV[2].ptInt.y - ptBarycenter.y) < EPS_SMALL &&
abs( LocCompV[3].ptInt.y - ptBarycenter.y) < EPS_SMALL) {
double dYBar = ptBarycenter.y / m_dStep - 0.5 ;
@@ -1131,9 +1126,9 @@ VolZmap::ExtMarchingCubes( int nBlock, TRIA3DLIST& lstTria, TriHolder& triHold)
++ nBarInd ;
}
}
else if ( abs( LocCompV[0].ptInt.z - ptBarycenter.z) < EPS_SMALL &&
abs( LocCompV[1].ptInt.z - ptBarycenter.z) < EPS_SMALL &&
abs( LocCompV[2].ptInt.z - ptBarycenter.z) < EPS_SMALL &&
else if ( abs( LocCompV[0].ptInt.z - ptBarycenter.z) < EPS_SMALL &&
abs( LocCompV[1].ptInt.z - ptBarycenter.z) < EPS_SMALL &&
abs( LocCompV[2].ptInt.z - ptBarycenter.z) < EPS_SMALL &&
abs( LocCompV[3].ptInt.z - ptBarycenter.z) < EPS_SMALL) {
double dZBar = ptBarycenter.z / m_dStep - 0.5 ;
@@ -1154,14 +1149,14 @@ VolZmap::ExtMarchingCubes( int nBlock, TRIA3DLIST& lstTria, TriHolder& triHold)
}
}
}
else if ( ( nIndArrey[nCompCount - 1][0] == 0 && nIndArrey[nCompCount - 1][1] == 1 &&
nIndArrey[nCompCount - 1][2] == 4 && nIndArrey[nCompCount - 1][3] == 5) ||
( nIndArrey[nCompCount - 1][0] == 1 && nIndArrey[nCompCount - 1][1] == 2 &&
nIndArrey[nCompCount - 1][2] == 5 && nIndArrey[nCompCount - 1][3] == 6) ||
( nIndArrey[nCompCount - 1][0] == 2 && nIndArrey[nCompCount - 1][1] == 3 &&
nIndArrey[nCompCount - 1][2] == 6 && nIndArrey[nCompCount - 1][3] == 7) ||
( nIndArrey[nCompCount - 1][0] == 0 && nIndArrey[nCompCount - 1][1] == 3 &&
nIndArrey[nCompCount - 1][2] == 4 && nIndArrey[nCompCount - 1][3] == 7)) {
else if ( ( nIndArrey[nCompCount - 1][0] == 0 && nIndArrey[nCompCount - 1][1] == 1 &&
nIndArrey[nCompCount - 1][2] == 4 && nIndArrey[nCompCount - 1][3] == 5) ||
( nIndArrey[nCompCount - 1][0] == 1 && nIndArrey[nCompCount - 1][1] == 2 &&
nIndArrey[nCompCount - 1][2] == 5 && nIndArrey[nCompCount - 1][3] == 6) ||
( nIndArrey[nCompCount - 1][0] == 2 && nIndArrey[nCompCount - 1][1] == 3 &&
nIndArrey[nCompCount - 1][2] == 6 && nIndArrey[nCompCount - 1][3] == 7) ||
( nIndArrey[nCompCount - 1][0] == 0 && nIndArrey[nCompCount - 1][1] == 3 &&
nIndArrey[nCompCount - 1][2] == 4 && nIndArrey[nCompCount - 1][3] == 7)) {
VectorField LocVecF[12], LocCompV[12] ;
@@ -1176,19 +1171,19 @@ VolZmap::ExtMarchingCubes( int nBlock, TRIA3DLIST& lstTria, TriHolder& triHold)
LocCompV[LocInd].vtNorm.ToLoc( m_MapFrame[0]) ;
}
if ( ( AreSameVectorApprox( LocVecF[nIndArrey[nCompCount - 1][0]].vtNorm, LocVecF[nIndArrey[nCompCount - 1][2]].vtNorm) &&
abs( LocVecF[nIndArrey[nCompCount - 1][0]].vtNorm * LocVecF[nIndArrey[nCompCount - 1][1]].vtNorm) < EPS_SMALL &&
abs( LocVecF[nIndArrey[nCompCount - 1][0]].vtNorm * LocVecF[nIndArrey[nCompCount - 1][3]].vtNorm) < EPS_SMALL) ||
( AreSameVectorApprox( LocVecF[nIndArrey[nCompCount - 1][1]].vtNorm, LocVecF[nIndArrey[nCompCount - 1][3]].vtNorm) &&
abs( LocVecF[nIndArrey[nCompCount - 1][1]].vtNorm * LocVecF[nIndArrey[nCompCount - 1][0]].vtNorm) < EPS_SMALL &&
if ( ( AreSameVectorApprox( LocVecF[nIndArrey[nCompCount - 1][0]].vtNorm, LocVecF[nIndArrey[nCompCount - 1][2]].vtNorm) &&
abs( LocVecF[nIndArrey[nCompCount - 1][0]].vtNorm * LocVecF[nIndArrey[nCompCount - 1][1]].vtNorm) < EPS_SMALL &&
abs( LocVecF[nIndArrey[nCompCount - 1][0]].vtNorm * LocVecF[nIndArrey[nCompCount - 1][3]].vtNorm) < EPS_SMALL) ||
( AreSameVectorApprox( LocVecF[nIndArrey[nCompCount - 1][1]].vtNorm, LocVecF[nIndArrey[nCompCount - 1][3]].vtNorm) &&
abs( LocVecF[nIndArrey[nCompCount - 1][1]].vtNorm * LocVecF[nIndArrey[nCompCount - 1][0]].vtNorm) < EPS_SMALL &&
abs( LocVecF[nIndArrey[nCompCount - 1][1]].vtNorm * LocVecF[nIndArrey[nCompCount - 1][2]].vtNorm) < EPS_SMALL)) {
Point3d ptBarycenter = ( LocCompV[0].ptInt + LocCompV[1].ptInt +
LocCompV[2].ptInt + LocCompV[3].ptInt) / 4 ;
if ( abs( LocCompV[0].ptInt.x - ptBarycenter.x) < EPS_SMALL &&
abs( LocCompV[1].ptInt.x - ptBarycenter.x) < EPS_SMALL &&
abs( LocCompV[2].ptInt.x - ptBarycenter.x) < EPS_SMALL &&
if ( abs( LocCompV[0].ptInt.x - ptBarycenter.x) < EPS_SMALL &&
abs( LocCompV[1].ptInt.x - ptBarycenter.x) < EPS_SMALL &&
abs( LocCompV[2].ptInt.x - ptBarycenter.x) < EPS_SMALL &&
abs( LocCompV[3].ptInt.x - ptBarycenter.x) < EPS_SMALL) {
double dXBar = ptBarycenter.x / m_dStep - 0.5 ;
@@ -1207,9 +1202,9 @@ VolZmap::ExtMarchingCubes( int nBlock, TRIA3DLIST& lstTria, TriHolder& triHold)
++ nBarInd ;
}
}
else if ( abs( LocCompV[0].ptInt.y - ptBarycenter.y) < EPS_SMALL &&
abs( LocCompV[1].ptInt.y - ptBarycenter.y) < EPS_SMALL &&
abs( LocCompV[2].ptInt.y - ptBarycenter.y) < EPS_SMALL &&
else if ( abs( LocCompV[0].ptInt.y - ptBarycenter.y) < EPS_SMALL &&
abs( LocCompV[1].ptInt.y - ptBarycenter.y) < EPS_SMALL &&
abs( LocCompV[2].ptInt.y - ptBarycenter.y) < EPS_SMALL &&
abs( LocCompV[3].ptInt.y - ptBarycenter.y) < EPS_SMALL) {
double dYBar = ptBarycenter.y / m_dStep - 0.5 ;
@@ -1228,9 +1223,9 @@ VolZmap::ExtMarchingCubes( int nBlock, TRIA3DLIST& lstTria, TriHolder& triHold)
++ nBarInd ;
}
}
else if ( abs( LocCompV[0].ptInt.z - ptBarycenter.z) < EPS_SMALL &&
abs( LocCompV[1].ptInt.z - ptBarycenter.z) < EPS_SMALL &&
abs( LocCompV[2].ptInt.z - ptBarycenter.z) < EPS_SMALL &&
else if ( abs( LocCompV[0].ptInt.z - ptBarycenter.z) < EPS_SMALL &&
abs( LocCompV[1].ptInt.z - ptBarycenter.z) < EPS_SMALL &&
abs( LocCompV[2].ptInt.z - ptBarycenter.z) < EPS_SMALL &&
abs( LocCompV[3].ptInt.z - ptBarycenter.z) < EPS_SMALL) {
double dZBar = ptBarycenter.z / m_dStep - 0.5 ;
@@ -1262,37 +1257,29 @@ VolZmap::ExtMarchingCubes( int nBlock, TRIA3DLIST& lstTria, TriHolder& triHold)
// Passo al sistema di riferimento del baricentro
Point3d ptGravityCenter( 0, 0, 0) ;
for ( int ni = 0 ; ni < nVertComp[nCompCount - 1] ; ++ ni)
ptGravityCenter = ptGravityCenter + CompoVert[nCompCount - 1][ni].ptInt ;
ptGravityCenter = ptGravityCenter / nVertComp[nCompCount - 1] ;
Vector3d vtO = ptGravityCenter - ORIG ;
Point3d ptTrasf[12] ;
Vector3d vtTrasf[12] ;
for ( int ni = 0 ; ni < nVertComp[nCompCount - 1] ; ++ ni)
ptTrasf[ni] = CompoVert[nCompCount - 1][ni].ptInt - vtO ;
vtTrasf[ni] = CompoVert[nCompCount - 1][ni].ptInt - ptGravityCenter ;
// Preparo le matrici per il sistema
typedef Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic> dSystemMatrix ;
typedef Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic> dSystemVector ;
typedef Eigen::JacobiSVD<dSystemMatrix> DecomposerSVD ;
dSystemMatrix dMatrixN, dMatrixU, dMatrixV ;
dSystemVector dKnownVector, dUnknownVector, dSingularValue ;
dMatrixN.resize( nVertComp[nCompCount - 1], 3) ;
dKnownVector.resize( nVertComp[nCompCount - 1], 1) ;
dUnknownVector.resize( 3, 1) ;
dSystemMatrix dMatrixN( nVertComp[nCompCount - 1], 3) ;
dSystemVector dKnownVector( nVertComp[nCompCount - 1], 1) ;
#if 0
// Studio del caso 4 punti su un piano
int nEqual = 0 ;
int nPosD ;
Vector3d vtD, vtE ;
if ( nVertComp[nCompCount - 1] == 4 && nFeatureType == EDGE) {
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) {
@@ -1326,20 +1313,20 @@ VolZmap::ExtMarchingCubes( int nBlock, TRIA3DLIST& lstTria, TriHolder& triHold)
( CompoVert[nCompCount - 1][3].ptInt - CompoVert[nCompCount - 1][2].ptInt))) ;
// Caso superficie piana
if ( nVertComp[nCompCount - 1] == 4 && nEqual == 2 && dDot < EPS_SMALL) {
if ( false && nVertComp[nCompCount - 1] == 4 && nEqual == 2 && dDot < EPS_SMALL) {
for ( int ni = 0 ; ni < nVertComp[nCompCount - 1] ; ++ ni) {
if ( ni != nPosD) {
dMatrixN( ni, 0) = CompoVert[nCompCount - 1][ni].vtNorm.x ;
dMatrixN( ni, 1) = CompoVert[nCompCount - 1][ni].vtNorm.y ;
dMatrixN( ni, 2) = CompoVert[nCompCount - 1][ni].vtNorm.z ;
dKnownVector( ni) = CompoVert[nCompCount - 1][ni].vtNorm * ( ptTrasf[ni] - ORIG) ;
dKnownVector( ni) = CompoVert[nCompCount - 1][ni].vtNorm * vtTrasf[ni] ;
}
else {
dMatrixN( ni, 0) = vtE.x ;
dMatrixN( ni, 1) = vtE.y ;
dMatrixN( ni, 2) = vtE.z ;
dKnownVector( ni) = vtE * ( ptTrasf[ni] - ORIG) ;
dKnownVector( ni) = vtE *vtTrasf[ni] ;
}
}
}
@@ -1349,89 +1336,54 @@ VolZmap::ExtMarchingCubes( int nBlock, TRIA3DLIST& lstTria, TriHolder& triHold)
dMatrixN( ni, 0) = CompoVert[nCompCount - 1][ni].vtNorm.x ;
dMatrixN( ni, 1) = CompoVert[nCompCount - 1][ni].vtNorm.y ;
dMatrixN( ni, 2) = CompoVert[nCompCount - 1][ni].vtNorm.z ;
dKnownVector( ni) = CompoVert[nCompCount - 1][ni].vtNorm * ( ptTrasf[ni] - ORIG) ;
dKnownVector( ni) = CompoVert[nCompCount - 1][ni].vtNorm * vtTrasf[ni] ;
}
}
#else
// Caso generale
for ( int ni = 0 ; ni < nVertComp[nCompCount - 1] ; ++ ni) {
dMatrixN( ni, 0) = CompoVert[nCompCount - 1][ni].vtNorm.x ;
dMatrixN( ni, 1) = CompoVert[nCompCount - 1][ni].vtNorm.y ;
dMatrixN( ni, 2) = CompoVert[nCompCount - 1][ni].vtNorm.z ;
dKnownVector( ni) = CompoVert[nCompCount - 1][ni].vtNorm * vtTrasf[ni] ;
}
#endif
// calcolo SVD
DecomposerSVD svd( dMatrixN, Eigen::ComputeThinU | Eigen::ComputeThinV) ;
dMatrixU = svd.matrixU() ;
dMatrixV = svd.matrixV() ;
dSingularValue = svd.singularValues() ;
double dS0 = dSingularValue( 0) ;
double dS1 = dSingularValue( 1) ;
double dS2 = dSingularValue( 2) ;
dSystemMatrix dMatrixV = svd.matrixV() ;
dSystemVector dSingularValue = svd.singularValues() ;
// Se la feature è un edge annullo il valore singolare minore.
if ( nFeatureType == EDGE) {
dSingularValue( 2) = 0 ;
}
// Back substitution: risolvo il sistema USV*x = b
// Calcolo U^T b
double vTemp[3] ;
double dThres = 0.5 * ( dSingularValue( 1) + dSingularValue( 2)) / dSingularValue( 0) ;
svd.setThreshold( dThres) ;
}
for ( int ni = 0 ; ni < 3 ; ++ ni) {
// risolvo il sistema con SVD, quindi ai minimi quadrati
dSystemVector dUnknownVector( 3, 1) ;
dUnknownVector = svd.solve( dKnownVector) ;
double s = 0 ;
if ( dSingularValue( ni) > 0) {
for ( int nj = 0 ; nj < nVertComp[nCompCount - 1] ; ++ nj)
s += dMatrixU( nj, ni) * dKnownVector( nj) ;
s /= dSingularValue( ni) ;
}
vTemp[ni] = s ;
}
// Moltiplico per V
for ( int ni = 0 ; ni < 3 ; ++ ni) {
double s = 0 ;
for ( int nj = 0 ; nj < 3 ; ++ nj)
s += dMatrixV( ni, nj) * vTemp[nj] ;
dUnknownVector( ni) = s ;
}
// Vettore Baricentro-Feature
Vector3d vtFeature( dUnknownVector( 0), dUnknownVector( 1), dUnknownVector( 2)) ;
// Esprimo la soluzione nel sistema di riferimento
// in cui è immerso quello dello dello z-map.
Point3d ptSol( dUnknownVector( 0) + vtO.x,
dUnknownVector( 1) + vtO.y,
dUnknownVector( 2) + vtO.z) ;
// Flag sulla distanza del vertice dal
// baricentro del sistema di punti
bool bOutside = false ;
// Vettore Baricentro-Feature.
Vector3d vtFeature( dUnknownVector( 0),
dUnknownVector( 1),
dUnknownVector( 2)) ;
// Esprimo la soluzione nel sistema di riferimento in cui è immerso quello dello z-map
Point3d ptSol = ptGravityCenter + vtFeature ;
// Limito la feature a una distanza dal baricentro pari alla diagonale del voxel
double dDistFeature = vtFeature.Len() ;
const double MAX_DIST = sqrt( 3) * m_dStep ;
// Limito la feature a una distanza dal
// baricentro pari alla diagonale del voxel.
if ( dDistFeature > MAX_DIST)
bOutside = true ;
bool bOutside = ( dDistFeature > MAX_DIST) ;
Triangle3d CurrentTriangle ;
TRIA3DVECTOR triContainer ;
// 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
// 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) ;
}
@@ -1440,7 +1392,7 @@ VolZmap::ExtMarchingCubes( int nBlock, TRIA3DLIST& lstTria, TriHolder& triHold)
bool bDangerInversion = false ;
// Caso ventaglio con tre vertici di base
if ( nVertComp[nCompCount - 1] == 3) {
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.
@@ -1452,11 +1404,12 @@ VolZmap::ExtMarchingCubes( int nBlock, TRIA3DLIST& lstTria, TriHolder& triHold)
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 ( dDI < - EPS_SMALL || dDJ < - EPS_SMALL) {
bInversione = true ;
break ;
}
}
// Se tale triangolo esiste procedo
if ( bInversione) {
@@ -1481,8 +1434,11 @@ VolZmap::ExtMarchingCubes( int nBlock, TRIA3DLIST& lstTria, TriHolder& triHold)
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)
if ( ( dDLast < EPS_SMALL && dDPrev < EPS_SMALL) ||
( dDLast < - 0.75 || dDPrev < - 0.75)) {
bInversione2 = true ;
break ;
}
}
// Se tale normale esiste
@@ -1518,16 +1474,16 @@ VolZmap::ExtMarchingCubes( int nBlock, TRIA3DLIST& lstTria, TriHolder& triHold)
// 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) ;
}
}
// Caso in cui non può essere riportata dentro:
// si passa alla routine standard se il voxel
// in cui cade è pieno.
@@ -1564,7 +1520,8 @@ VolZmap::ExtMarchingCubes( int nBlock, TRIA3DLIST& lstTria, TriHolder& triHold)
}
}
}
}
}
// Ventaglio con base a quattro vertici
else if ( nVertComp[nCompCount - 1] == 4) {
@@ -1625,23 +1582,21 @@ VolZmap::ExtMarchingCubes( int nBlock, TRIA3DLIST& lstTria, TriHolder& triHold)
// 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
// 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))
@@ -1660,7 +1615,6 @@ VolZmap::ExtMarchingCubes( int nBlock, TRIA3DLIST& lstTria, TriHolder& triHold)
nAdjIndex |= ( 1 << 6) ;
if ( IsThereMat( nAdjVoxI, nAdjVoxJ + 1, nAdjVoxK + 1))
nAdjIndex |= ( 1 << 7) ;
// Se il voxel è pieno
if ( EdgeTable[nAdjIndex] != 0)
bDangerInversion = true ;
@@ -1674,7 +1628,6 @@ VolZmap::ExtMarchingCubes( int nBlock, TRIA3DLIST& lstTria, TriHolder& triHold)
Point3d ptSolZMapFrame = ptSol ;
ptSolZMapFrame.ToLoc( m_MapFrame[0]) ;
if ( ! IsPointInsideVoxelApprox( i, j, k, ptSolZMapFrame)) {
@@ -1699,16 +1652,87 @@ VolZmap::ExtMarchingCubes( int nBlock, TRIA3DLIST& lstTria, TriHolder& triHold)
// 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) ;
// 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 ;
}
}
}
}
}
}
}
}
@@ -1741,14 +1765,13 @@ VolZmap::ExtMarchingCubes( int nBlock, TRIA3DLIST& lstTria, TriHolder& triHold)
// Se la feature non è fuori dai limiti
// e i triangoli formano una superficie
// non piana confermo ExtMC
if ( ( ! ( bOutside || bPlane)) && ( ! bDangerInversion)) {
if ( ! ( bOutside || bPlane) && ! bDangerInversion) {
TRIA3DVECTOR vInnerTriaTemp, vBorderTriaTemp ;
for ( int ni = 0 ; ni < nContSize ; ++ ni) {
// Se l'area è non nulla aggiungo il triangolo
// a uno dei due vettori.
// Se l'area è non nulla aggiungo il triangolo a uno dei due vettori.
if ( triContainer[ni].GetArea() > EPS_SMALL * EPS_SMALL) {
int nVoxIJK[3] = { i, j, k} ;
@@ -1759,7 +1782,6 @@ VolZmap::ExtMarchingCubes( int nBlock, TRIA3DLIST& lstTria, TriHolder& triHold)
trTemp.ToLoc( m_MapFrame[0]) ;
if ( ! IsATriangleOnBorder( trTemp, ptTemp, nLimits, nVoxIJK, bBoundary))
vInnerTriaTemp.emplace_back( triContainer[ni]) ;
else
vBorderTriaTemp.emplace_back( triContainer[ni]) ;
@@ -1791,8 +1813,7 @@ VolZmap::ExtMarchingCubes( int nBlock, TRIA3DLIST& lstTria, TriHolder& triHold)
triHold[nCurrent].k = k ;
}
// Indice che identifica la posizione del voxel
// nel vector.
// Indice che identifica la posizione del voxel nel vector.
int nCurrent = int( triHold.size()) - 1 ;
// Aggiungo vertice della componente
@@ -1805,9 +1826,7 @@ VolZmap::ExtMarchingCubes( int nBlock, TRIA3DLIST& lstTria, TriHolder& triHold)
// Aggiungo una componente per il vettore
// dei triangoli della componente connessa.
triHold[nCurrent].vCompoTria.resize( nNewFeatureNum) ;
for ( int ni = 0 ; ni < int( vInnerTriaTemp.size()) ; ++ ni)
triHold[nCurrent].vCompoTria[nOldFeatureNum].emplace_back( vInnerTriaTemp[ni]) ;
}
@@ -1834,60 +1853,48 @@ VolZmap::ExtMarchingCubes( int nBlock, TRIA3DLIST& lstTria, TriHolder& triHold)
m_InterBlockTria[nBlock][nCurrent].k = k ;
}
// Indice che identifica la posizione del voxel
// nel vector.
// Indice che identifica la posizione del voxel nel vector
int nCurrent = int( m_InterBlockTria[nBlock].size()) - 1 ;
// Aggiungo vertice della componente
// connessa alla lista dei vertici.
// Aggiungo vertice della componente connessa alla lista dei vertici
m_InterBlockTria[nBlock][nCurrent].ptCompoVert.emplace_back( ptSol) ;
int nOldFeatureNum = int( m_InterBlockTria[nBlock][nCurrent].vCompoTria.size()) ;
int nNewFeatureNum = nOldFeatureNum + 1 ;
// Aggiungo una componente per il vettore
// dei triangoli della componente connessa.
// Aggiungo una componente per il vettore dei triangoli della componente connessa
m_InterBlockTria[nBlock][nCurrent].vCompoTria.resize( nNewFeatureNum) ;
for ( int ni = 0 ; ni < int( vBorderTriaTemp.size()) ; ++ ni)
m_InterBlockTria[nBlock][nCurrent].vCompoTria[nOldFeatureNum].emplace_back( vBorderTriaTemp[ni]) ;
}
}
// ExtMC non confermato, si passa a MC
else {
// Costruzione dei triangoli
for ( int TriIndex = 0; TriIndex < ( nVertComp[nCompCount - 1] - 2) * 3 ; TriIndex += 3) {
Triangle3d CurrentTriangle ;
// Il triangolo è pronto
Triangle3d CurrentTriangle ;
CurrentTriangle.Set( CompoTriVert[nCompCount - 1][TriIndex].ptInt,
CompoTriVert[nCompCount - 1][TriIndex+1].ptInt,
CompoTriVert[nCompCount - 1][TriIndex+2].ptInt) ;
bool bV = CurrentTriangle.Validate( true) ;
// Aggiungo alla lista
lstTria.emplace_back( CurrentTriangle) ;
}
}
}
// Standard MC
else {
// Costruzione dei triangoli
for ( int TriIndex = 0; TriIndex < ( nVertComp[nCompCount - 1] - 2) * 3 ; TriIndex += 3) {
Triangle3d CurrentTriangle ;
// Il triangolo è pronto
Triangle3d CurrentTriangle ;
CurrentTriangle.Set( CompoTriVert[nCompCount - 1][TriIndex].ptInt,
CompoTriVert[nCompCount - 1][TriIndex+1].ptInt,
CompoTriVert[nCompCount - 1][TriIndex+2].ptInt) ;
bool bV = CurrentTriangle.Validate( true) ;
// Aggiungo alla lista
lstTria.emplace_back( CurrentTriangle) ;
}
@@ -1899,7 +1906,6 @@ VolZmap::ExtMarchingCubes( int nBlock, TRIA3DLIST& lstTria, TriHolder& triHold)
return true ;
}
//
//----------------------------------------------------------------------------
bool