diff --git a/EgtGeomKernel.rc b/EgtGeomKernel.rc index 02c78a1..85102c0 100644 Binary files a/EgtGeomKernel.rc and b/EgtGeomKernel.rc differ diff --git a/VolTriZmapCreation.cpp b/VolTriZmapCreation.cpp index 970e75b..0bf1aef 100644 --- a/VolTriZmapCreation.cpp +++ b/VolTriZmapCreation.cpp @@ -75,7 +75,7 @@ VolZmap::Create( const Point3d& ptO, double dLengthX, double dLengthY, double dL // Definisco il numero di blocchi lungo z m_nFracLin[2] = max( nMinBlockNum, m_nNy[1] / m_nDexNumPBlock + - ( m_nNy[1] % m_nDexNumPBlock >= m_nDexNumPBlock / 2 ? 1 : 0)) ; + ( m_nNy[1] % m_nDexNumPBlock >= m_nDexNumPBlock / 2 ? 1 : 0)) ; } else { @@ -153,6 +153,8 @@ VolZmap::Create( const Point3d& ptO, double dLengthX, double dLengthY, double dL for ( unsigned int nCount = 0 ; nCount < m_nNumBlock ; ++ nCount) m_BlockToUpdate[nCount] = true ; + + m_InterBlockTria.resize( m_nNumBlock) ; // Aggiornamento dello stato m_nStatus = OK ; @@ -479,7 +481,9 @@ VolZmap::CreateFromFlatRegion( const ISurfFlatRegion& Surf, double dDimZ, double for ( unsigned int nCount = 0 ; nCount < m_nNumBlock ; ++ nCount) m_BlockToUpdate[nCount] = true ; - + + m_InterBlockTria.resize( m_nNumBlock) ; + // Aggiornamento dello stato m_nStatus = OK ; @@ -815,8 +819,10 @@ VolZmap::CreateFromTriMesh( const ISurfTriMesh& Surf, double dPrec, bool bFlag) for ( unsigned int nCount = 0 ; nCount < m_nNumBlock ; ++ nCount) - m_BlockToUpdate[nCount] = true ; + m_BlockToUpdate[nCount] = true ; + m_InterBlockTria.resize( m_nNumBlock) ; + m_nStatus = OK ; return true ; diff --git a/VolTriZmapGraphics.cpp b/VolTriZmapGraphics.cpp index bd4a5eb..5bb27e9 100644 --- a/VolTriZmapGraphics.cpp +++ b/VolTriZmapGraphics.cpp @@ -332,35 +332,107 @@ VolZmap::GetTriangles( bool bAllBlocks, INTVECTOR& nModifiedBlocks, TRIA3DLISTVE // Caso con tre mappe else { - nModifiedBlocks.resize( m_nNumBlock) ; - vLstTria.reserve( m_nNumBlock) ; + nModifiedBlocks.resize( m_nNumBlock + 1) ; + vLstTria.reserve( m_nNumBlock + 1) ; - std::vector VecTriHold ; + TriaMatrix VecTriHold ; VecTriHold.resize( m_nNumBlock) ; - + // Ciclo sui blocchi for ( size_t t = 0 ; t < m_nNumBlock ; ++ t) { + // Calcolo i limiti sugli indici dei voxel del blocco int nIJK[3] ; // Vettore indici i,j,k del blocco - GetBlockIJK( int( t), nIJK) ; + GetBlockIJKFromN( int( t), nIJK) ; // Vettore limiti sugli indici dei voxel del blocco int LimitIndexes[6] ; GetBlockLimitsIJK( nIJK, LimitIndexes) ; + // Se il blocco deve essere processato, eseguo if ( bAllBlocks || m_BlockToUpdate[t]) { + vLstTria.emplace_back() ; nModifiedBlocks[t] = int( vLstTria.size()) - 1 ; - ExtMarchingCubes( LimitIndexes, vLstTria.back(), VecTriHold[t]) ; + m_InterBlockTria[t].clear() ; + ExtMarchingCubes( LimitIndexes, t, vLstTria.back(), VecTriHold[t]) ; + + // Flipping fra voxel interni + FlipEdgesII( VecTriHold[t]) ; + + //for ( int nDeltaI = -1 ; nDeltaI <= 1 ; ++ nDeltaI) { + // for ( int nDeltaJ = -1 ; nDeltaJ <= 1 ; ++ nDeltaJ) { + // for ( int nDeltaK = -1 ; nDeltaK <= 1 ; ++ nDeltaK) { + // + // // Flipping fra voxel dello stesso blocco + // if ( nDeltaI == 0 && nDeltaJ == 0 && nDeltaK == 0) { + // + // // Flipping fra voxel interni + // FlipEdgesII( VecTriHold[t]) ; + // // Flipping fra voxel interni e di frontiera + // //FlipEdgesIB( VecTriHold[t], t) ; + // } + // // Flipping fra voxel adiacenti + // else { + + // int nAdjN ; + // + // if ( GetAdjBlockToBlock( int( t), nDeltaI, nDeltaJ, nDeltaK, nAdjN)) { + // + // // Qui bisogna fare il flipping fra i blocchi: t indica il blocco + // // corrente appena aggiornato. Guardiamo il suo vettore di voxel bislacchi + // // e cicliamo su tali voxel. Per ogni voxel bislacco guardiamo i voxel adiacenti. + // // Se un voxel adiacente è nel blocco adiacente corrente ricalcoliamo i triangoli e flippiamo. + // size_t nBoundaryVoxelNum = m_InterBlockTria[t].size() ; + + // for ( size_t a = 0 ; a < nBoundaryVoxelNum ; ++ a) { + // + // // vedere se il voxel appartiene al blocco adiacente + // int nBVoxI = m_InterBlockTria[t][a].i ; + // int nBVoxJ = m_InterBlockTria[t][a].j ; + // int nBVoxK = m_InterBlockTria[t][a].k ; + + // // Indici blocco adiacente + // int nAdjBlockIJK[3] ; + // GetBlockIJKFromN( nAdjN, nAdjBlockIJK) ; + // // Limiti indici di voxel per il blocco adiacente + // int nLimitAdjBlock[6] ; + + // GetBlockLimitsIJK( nAdjBlockIJK, nLimitAdjBlock) ; + // + // // Se il blocco adiacente esiste eseguo il flipping + // if ( nBVoxI >= nLimitAdjBlock[0] && nBVoxI < nLimitAdjBlock[1] && + // nBVoxJ >= nLimitAdjBlock[2] && nBVoxJ < nLimitAdjBlock[3] && + // nBVoxK >= nLimitAdjBlock[4] && nBVoxK < nLimitAdjBlock[5]) { + // + // m_InterBlockTria[nAdjN].clear() ; + + + + // } + // } + // } + // } + // } + // } + //} m_BlockToUpdate[t] = false ; } else nModifiedBlocks[t] = -1 ; } - // Eseguo il flipping - FlipEdges( VecTriHold) ; + // Copio i triangoli di frontiera in una matrice gemella + // di m_InterBlockTria per avere sempre a disposizione + // i triangoli non flippati. + TriaMatrix InterBlockTria ; + InterBlockTria = m_InterBlockTria ; + + FlipEdgesBB( InterBlockTria) ; + + + // ciclo sui blocchi for ( size_t t = 0 ; t < m_nNumBlock ; ++ t) { // ciclo sui voxel del blocco @@ -375,7 +447,28 @@ VolZmap::GetTriangles( bool bAllBlocks, INTVECTOR& nModifiedBlocks, TRIA3DLISTVE } } } - } + } + + + vLstTria.resize( vLstTria.size() + 1) ; + size_t nPos = size_t( vLstTria.size() - 1) ; + + for ( size_t t = 0 ; t < m_InterBlockTria.size() ; ++ t) { + for ( size_t t1 = 0 ; t1 < m_InterBlockTria[t].size() ; ++ t1) { + for ( size_t t2 = 0 ; t2 < m_InterBlockTria[t][t1].vCompoTria.size() ; ++ t2) { + for ( size_t t3 = 0 ; t3 < m_InterBlockTria[t][t1].vCompoTria[t2].size() ; ++ t3) { + + if ( m_InterBlockTria[t][t1].vCompoTria[t2][t3].GetArea() > EPS_SMALL) { + + + vLstTria[nPos].emplace_back( InterBlockTria[t][t1].vCompoTria[t2][t3]) ; + } + } + } + } + } + + nModifiedBlocks[nModifiedBlocks.size()-1] = int( vLstTria.size() - 1) ; } return true ; @@ -385,7 +478,7 @@ VolZmap::GetTriangles( bool bAllBlocks, INTVECTOR& nModifiedBlocks, TRIA3DLISTVE int VolZmap::GetBlockCount( void) const { - return m_nNumBlock ; + return m_nNumBlock + 1 ; } //---------------------------------------------------------------------------- @@ -792,7 +885,7 @@ VolZmap::MarchingCubes( int nBlock, TRIA3DLIST& lstTria) const //---------------------------------------------------------------------------- bool -VolZmap::ExtMarchingCubes( const int nLimits[], TRIA3DLIST& lstTria, TriHolder& triHold) const +VolZmap::ExtMarchingCubes( const int nLimits[], size_t nBlockNumber, TRIA3DLIST& lstTria, TriHolder& triHold) const { Point3d ptMapOrig = m_MapFrame[0].Orig() ; @@ -801,7 +894,17 @@ VolZmap::ExtMarchingCubes( const int nLimits[], TRIA3DLIST& lstTria, TriHolder& for ( int j = nLimits[2] ; j < nLimits[3] ; ++ j) { for ( int k = nLimits[4] ; k < nLimits[5] ; ++ k) { - + if ( !(i == 31 && j == -1 && k == 49)) + ;// continue ; + + // Riconoscimento dei voxel di frontiera + bool bBoundary = false ; + int nVoxIndexes[3] = { i, j, k} ; + + if ( IsAVoxelOnBoundary( nLimits, nVoxIndexes, true)) + bBoundary = true ; + + // Indici i,j,k dei vertici int IndexCorner[8][3] = { { i, j, k}, @@ -891,7 +994,8 @@ VolZmap::ExtMarchingCubes( const int nLimits[], TRIA3DLIST& lstTria, TriHolder& // Numero di feature nel voxel: al più vi è // una feature per componente connessa. - int nFeatureInVoxel = 0 ; + int nInnerFeatureInVoxel = 0 ; + int nBorderFeatureInVoxel = 0 ; // Ciclo sulle componenti for ( int nCompCount = 1 ; nCompCount <= nComponents ; ++ nCompCount) { @@ -907,8 +1011,166 @@ VolZmap::ExtMarchingCubes( const int nLimits[], TRIA3DLIST& lstTria, TriHolder& // Nota che il primo elemento dell'array (0-esimo) non viene inizializzato CompoVert[nVertCount] = VecField[TriangleTableEn[nIndex][1][nVertCount + nTableOffset + 1]] ; - int nFeatureType ; + // Controllo per il caso piano su una griglia + // con versori normali a due a due paralleli. + bool bGridControl = true ; + if ( nVertComp == 4) { + + int nIndArrey[4] ; + // Riempio un vettore con gli indici dei punti + for ( int nCpInd = 0 ; nCpInd < nVertComp ; ++ nCpInd) + nIndArrey[nCpInd] = TriangleTableEn[nIndex][1][nCpInd + nTableOffset + 1] ; + // Ordino i 4 indici in senso crescente + for ( int nSrtInd1 = 0 ; nSrtInd1 < nVertComp - 1 ; ++ nSrtInd1) { + for ( int nSrtInd2 = nSrtInd1 + 1 ; nSrtInd2 < nVertComp ; ++ nSrtInd2) { + if ( nIndArrey[nSrtInd1] > nIndArrey[nSrtInd2]) + swap( nIndArrey[nSrtInd1], nIndArrey[nSrtInd2]) ; + } + } + + if ( ( nIndArrey[0] == 0 && nIndArrey[1] == 2 && nIndArrey[2] == 9 && nIndArrey[3] == 10) || + ( nIndArrey[0] == 4 && nIndArrey[1] == 6 && nIndArrey[2] == 9 && nIndArrey[3] == 10) || + ( nIndArrey[0] == 4 && nIndArrey[1] == 6 && nIndArrey[2] == 8 && nIndArrey[3] == 11) || + ( nIndArrey[0] == 0 && nIndArrey[1] == 2 && nIndArrey[2] == 8 && nIndArrey[3] == 11) || + ( nIndArrey[0] == 1 && nIndArrey[1] == 3 && nIndArrey[2] == 8 && nIndArrey[3] == 9 ) || + ( nIndArrey[0] == 1 && nIndArrey[1] == 3 && nIndArrey[2] == 10 && nIndArrey[3] == 11) || + ( nIndArrey[0] == 5 && nIndArrey[1] == 7 && nIndArrey[2] == 10 && nIndArrey[3] == 11) || + ( nIndArrey[0] == 5 && nIndArrey[1] == 7 && nIndArrey[2] == 8 && nIndArrey[3] == 9 )) { + + if ( AreSameVectorApprox( VecField[nIndArrey[0]].vtNorm, VecField[nIndArrey[1]].vtNorm) && + AreSameVectorApprox( VecField[nIndArrey[2]].vtNorm, VecField[nIndArrey[3]].vtNorm) && + abs( VecField[nIndArrey[0]].vtNorm * VecField[nIndArrey[2]].vtNorm) < EPS_SMALL) { + + Point3d ptBarycenter = ( CompoVert[0].ptInt + CompoVert[1].ptInt + + CompoVert[2].ptInt + CompoVert[3].ptInt) / 4 ; + + if ( abs( CompoVert[0].ptInt.x - ptBarycenter.x) < EPS_SMALL && + abs( CompoVert[1].ptInt.x - ptBarycenter.x) < EPS_SMALL && + abs( CompoVert[2].ptInt.x - ptBarycenter.x) < EPS_SMALL && + abs( CompoVert[3].ptInt.x - ptBarycenter.x) < EPS_SMALL) { + + double dXBar = ptBarycenter.x / m_dStep - 0.5 ; + int nBarLimSup = int( m_nNx[0]) ; + int nBarInd = 0 ; + + while ( nBarInd < nBarLimSup) { + + double dXInd = double( nBarInd) ; + + if ( abs( dXInd - dXBar) < EPS_SMALL) { + + bGridControl = false ; + break ; + } + ++ nBarInd ; + } + } + else if ( abs( CompoVert[0].ptInt.y - ptBarycenter.y) < EPS_SMALL && + abs( CompoVert[1].ptInt.y - ptBarycenter.y) < EPS_SMALL && + abs( CompoVert[2].ptInt.y - ptBarycenter.y) < EPS_SMALL && + abs( CompoVert[3].ptInt.y - ptBarycenter.y) < EPS_SMALL) { + + double dYBar = ptBarycenter.y / m_dStep - 0.5 ; + int nBarLimSup = int( m_nNy[0]) ; + int nBarInd = 0 ; + + while ( nBarInd < nBarLimSup) { + + double dYInd = double( nBarInd) ; + + if ( abs( dYInd - dYBar) < EPS_SMALL) { + + bGridControl = false ; + break ; + } + ++ nBarInd ; + } + } + else if ( abs( CompoVert[0].ptInt.z - ptBarycenter.z) < EPS_SMALL && + abs( CompoVert[1].ptInt.z - ptBarycenter.z) < EPS_SMALL && + abs( CompoVert[2].ptInt.z - ptBarycenter.z) < EPS_SMALL && + abs( CompoVert[3].ptInt.z - ptBarycenter.z) < EPS_SMALL) { + + double dZBar = ptBarycenter.z / m_dStep - 0.5 ; + int nBarLimSup = int( m_nNy[1]) ; + int nBarInd = 0 ; + + while ( nBarInd < nBarLimSup) { + + double dZInd = double( nBarInd) ; + + if ( abs( dZInd - dZBar) < EPS_SMALL) { + + bGridControl = false ; + break ; + } + ++ nBarInd ; + } + } + } + } + else if ( ( nIndArrey[0] == 0 && nIndArrey[1] == 1 && nIndArrey[2] == 4 && nIndArrey[3] == 5) || + ( nIndArrey[0] == 1 && nIndArrey[1] == 2 && nIndArrey[2] == 5 && nIndArrey[3] == 6) || + ( nIndArrey[0] == 2 && nIndArrey[1] == 3 && nIndArrey[2] == 6 && nIndArrey[3] == 7) || + ( nIndArrey[0] == 0 && nIndArrey[1] == 3 && nIndArrey[2] == 4 && nIndArrey[3] == 7)) { + + if ( AreSameVectorApprox( VecField[nIndArrey[0]].vtNorm, VecField[nIndArrey[2]].vtNorm) && + AreSameVectorApprox( VecField[nIndArrey[1]].vtNorm, VecField[nIndArrey[3]].vtNorm) && + abs( VecField[nIndArrey[0]].vtNorm * VecField[nIndArrey[1]].vtNorm) < EPS_SMALL) { + + Point3d ptBarycenter = ( CompoVert[0].ptInt + CompoVert[1].ptInt + + CompoVert[2].ptInt + CompoVert[3].ptInt) / 4 ; + + if ( abs( CompoVert[0].ptInt.x - ptBarycenter.x) < EPS_SMALL && + abs( CompoVert[1].ptInt.x - ptBarycenter.x) < EPS_SMALL && + abs( CompoVert[2].ptInt.x - ptBarycenter.x) < EPS_SMALL && + abs( CompoVert[3].ptInt.x - ptBarycenter.x) < EPS_SMALL) { + + double dXBar = ptBarycenter.x / m_dStep - 0.5 ; + int nBarLimSup = int( m_nNx[0]) ; + int nBarInd = 0 ; + + while ( nBarInd < nBarLimSup) { + + double dXInd = double( nBarInd) ; + + if ( abs( dXInd - dXBar) < EPS_SMALL) { + + bGridControl = false ; + break ; + } + ++ nBarInd ; + } + } + else if ( abs( CompoVert[0].ptInt.y - ptBarycenter.y) < EPS_SMALL && + abs( CompoVert[1].ptInt.y - ptBarycenter.y) < EPS_SMALL && + abs( CompoVert[2].ptInt.y - ptBarycenter.y) < EPS_SMALL && + abs( CompoVert[3].ptInt.y - ptBarycenter.y) < EPS_SMALL) { + + double dYBar = ptBarycenter.y / m_dStep - 0.5 ; + int nBarLimSup = int( m_nNy[0]) ; + int nBarInd = 0 ; + + while ( nBarInd < nBarLimSup) { + + double dYInd = double( nBarInd) ; + + if ( abs( dYInd - dYBar) < EPS_SMALL) { + + bGridControl = false ; + break ; + } + ++ nBarInd ; + } + } + } + } + } + + + int nFeatureType ; + // Valuto le relazioni reciproche fra le normali e // se i punti sono su un piano di griglia. Vector3d vtFeatureAxis ; @@ -926,10 +1188,10 @@ VolZmap::ExtMarchingCubes( const int nLimits[], TRIA3DLIST& lstTria, TriHolder& bNormal = TestOnNormal( CompoVert, nVertComp, nFeatureType, vtFeatureAxis, nMin1, nMin2) ; // Flag ExtMC - bool bExtMC = bNormal ; + bool bExtMC = bNormal && bGridControl ; // Extended MC - if ( bExtMC) { + if ( bExtMC) { // Passo al sistema di riferimento del baricentro Point3d ptGravityCenter( 0, 0, 0) ; @@ -1297,55 +1559,158 @@ VolZmap::ExtMarchingCubes( const int nLimits[], TRIA3DLIST& lstTria, TriHolder& break ; } - // Se la feaure non è fuori dai limiti + // Se la feature non è fuori dai limiti // e i triangoli formano una superficie // non piana confermo ExtMC if ( ! ( bOutside || bPlane)) { // Aggiorno il numero di feature. - ++ nFeatureInVoxel ; + // ++ nInnerFeatureInVoxel ; - // Se siamo in presenza della prima feature del - // voxel, ridimensiono il vettore che contiene - // la struttura dati del voxel. - if ( nFeatureInVoxel == 1) { + //// Se siamo in presenza della prima feature del + //// voxel, ridimensiono il vettore che contiene + //// la struttura dati del voxel. + // if ( nInnerFeatureInVoxel == 1) { - triHold.resize( triHold.size() + 1) ; + // triHold.resize( triHold.size() + 1) ; - // Questi dati dipendono solo dal voxel, - // quindi sono validi per tutte le - // componenti che vi appartengono. - int nCurrent = int( triHold.size()) - 1 ; + // // Questi dati dipendono solo dal voxel, + // // quindi sono validi per tutte le + // // componenti che vi appartengono. + // int nCurrent = int( triHold.size()) - 1 ; - triHold[nCurrent].i = i ; - triHold[nCurrent].j = j ; - triHold[nCurrent].k = k ; - } + // triHold[nCurrent].i = i ; + // triHold[nCurrent].j = j ; + // triHold[nCurrent].k = k ; + // } - // Indice che identifica la posizione del voxel - // nel vector. - int nCurrent = int( triHold.size()) - 1 ; + //// Indice che identifica la posizione del voxel + //// nel vector. + // int nCurrent = int( triHold.size()) - 1 ; - // Aggiungo vertice della componente - // connessa alla lista dei vertici. - triHold[nCurrent].ptCompoVert.emplace_back( ptSol) ; + //// Aggiungo vertice della componente + //// connessa alla lista dei vertici. + // triHold[nCurrent].ptCompoVert.emplace_back( ptSol) ; - int nOldFeatureNum = int( triHold[nCurrent].vCompoTria.size()) ; - int nNewFeatureNum = nOldFeatureNum + 1 ; + // int nOldFeatureNum = int( triHold[nCurrent].vCompoTria.size()) ; + // int nNewFeatureNum = nOldFeatureNum + 1 ; - // Aggiungo una componente per il vettore - // dei triangoli della componente connessa. - triHold[nCurrent].vCompoTria.resize( nNewFeatureNum) ; + //// Aggiungo una componente per il vettore + //// dei triangoli della componente connessa. + // triHold[nCurrent].vCompoTria.resize( nNewFeatureNum) ; + + TRIA3DVECTOR vInnerTriaTemp, vBorderTriaTemp ; for ( int ni = 0 ; ni < nContSize ; ++ ni) { // Se l'area è non nulla aggiungo il triangolo + // a uno dei due vettori. if ( triContainer[ni].GetArea() > EPS_SMALL) { - - triHold[nCurrent].vCompoTria[nOldFeatureNum].emplace_back( triContainer[ni]) ; + + int nVoxIJK[3] = { i, j, k} ; + + Point3d ptTemp = ptSol ; + ptTemp.ToLoc( m_MapFrame[0]) ; + Triangle3d trTemp = triContainer[ni] ; + trTemp.ToLoc( m_MapFrame[0]) ; + + if ( ! IsATriangleOnBorder( trTemp, ptTemp, nLimits, nVoxIJK, bBoundary)) + + vInnerTriaTemp.emplace_back( triContainer[ni]) ; + else + vBorderTriaTemp.emplace_back( triContainer[ni]) ; } - } + } + + // Metto i triangoli negli appositi contenitori: + + // Triangoli interni + if ( vInnerTriaTemp.size() > 0) { + + // Aggiorno il numero di feature. + ++ nInnerFeatureInVoxel ; + + // Se siamo in presenza della prima feature del + // voxel, ridimensiono il vettore che contiene + // la struttura dati del voxel. + if ( nInnerFeatureInVoxel == 1) { + + triHold.resize( triHold.size() + 1) ; + + // Questi dati dipendono solo dal voxel, + // quindi sono validi per tutte le + // componenti che vi appartengono. + int nCurrent = int( triHold.size()) - 1 ; + + triHold[nCurrent].i = i ; + triHold[nCurrent].j = j ; + triHold[nCurrent].k = k ; + } + + // Indice che identifica la posizione del voxel + // nel vector. + int nCurrent = int( triHold.size()) - 1 ; + + + // Aggiungo vertice della componente + // connessa alla lista dei vertici. + triHold[nCurrent].ptCompoVert.emplace_back( ptSol) ; + + int nOldFeatureNum = int( triHold[nCurrent].vCompoTria.size()) ; + int nNewFeatureNum = nOldFeatureNum + 1 ; + + // 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]) ; + } + + // Triangoli di frontiera + if ( vBorderTriaTemp.size() > 0) { + + // Aggiorno il numero di feature. + ++ nBorderFeatureInVoxel ; + + // Se siamo in presenza della prima feature del + // voxel, ridimensiono il vettore che contiene + // la struttura dati del voxel. + if ( nBorderFeatureInVoxel == 1) { + + m_InterBlockTria[nBlockNumber].resize( m_InterBlockTria[nBlockNumber].size() + 1) ; + + // Questi dati dipendono solo dal voxel, + // quindi sono validi per tutte le + // componenti che vi appartengono. + int nCurrent = int( m_InterBlockTria[nBlockNumber].size()) - 1 ; + + m_InterBlockTria[nBlockNumber][nCurrent].i = i ; + m_InterBlockTria[nBlockNumber][nCurrent].j = j ; + m_InterBlockTria[nBlockNumber][nCurrent].k = k ; + } + + // Indice che identifica la posizione del voxel + // nel vector. + int nCurrent = int( m_InterBlockTria[nBlockNumber].size()) - 1 ; + + // Aggiungo vertice della componente + // connessa alla lista dei vertici. + m_InterBlockTria[nBlockNumber][nCurrent].ptCompoVert.emplace_back( ptSol) ; + + int nOldFeatureNum = int( m_InterBlockTria[nBlockNumber][nCurrent].vCompoTria.size()) ; + int nNewFeatureNum = nOldFeatureNum + 1 ; + + // Aggiungo una componente per il vettore + // dei triangoli della componente connessa. + m_InterBlockTria[nBlockNumber][nCurrent].vCompoTria.resize( nNewFeatureNum) ; + + for ( int ni = 0 ; ni < int( vBorderTriaTemp.size()) ; ++ ni) + + m_InterBlockTria[nBlockNumber][nCurrent].vCompoTria[nOldFeatureNum].emplace_back( vBorderTriaTemp[ni]) ; + } } // ExtMC non confermato, si passa a MC else { @@ -1391,9 +1756,8 @@ VolZmap::ExtMarchingCubes( const int nLimits[], TRIA3DLIST& lstTria, TriHolder& lstTria.emplace_back( CurrentTriangle) ; } } - nTableOffset += nVertComp ; - } + } } } } @@ -1401,9 +1765,148 @@ VolZmap::ExtMarchingCubes( const int nLimits[], TRIA3DLIST& lstTria, TriHolder& return true ; } +//---------------------------------------------------------------------------- +//bool +//VolZmap::FlipEdges( std::vector& VecTriHold) const +//{ +// double dSqEps = EPS_SMALL * EPS_SMALL ; +// +// // Ciclo sui blocchi +// for ( size_t t1 = 0 ; t1 < m_nNumBlock ; ++ t1) { +// +// // Determino i,j,k del primo blocco +// int nIJK1[3] ; +// GetBlockIJK( int( t1), nIJK1) ; +// +// for ( size_t t2 = t1 ; t2 < m_nNumBlock ; ++ t2) { +// +// // Determino i,j,k del secondo blocco +// int nIJK2[3] ; +// GetBlockIJK( int( t2), nIJK2) ; +// +// // Se i blocchi sono adiacenti o coincidenti proseguo +// if ( ( nIJK2[0] >= nIJK1[0] - 1 && nIJK2[0] <= nIJK1[0] + 1) || +// ( nIJK2[1] >= nIJK1[1] - 1 && nIJK2[1] <= nIJK1[1] + 1) || +// ( nIJK2[2] >= nIJK1[2] - 1 && nIJK2[2] <= nIJK1[2] + 1)) { +// +// // Numero di voxel in cui si presentano sharp feature +// int nVoxelNum1 = int( VecTriHold[t1].size()) ; +// int nVoxelNum2 = int( VecTriHold[t2].size()) ; +// +// // Determino estremi del ciclo sui voxel esterno +// int nSt1 = 0 ; +// int nEn1 = nVoxelNum1 - ( t1 == t2 ? 1 : 0) ; +// +// // Ciclo su tali voxel +// for ( int n1 = nSt1 ; n1 < nEn1 ; ++ n1) { +// +// // Determino estremi del ciclo sui voxel interno +// int nSt2 = ( t1 == t2 ? nSt1 : 0) ; +// int nEn2 = nVoxelNum2 ; +// +// for ( int n2 = n1 ; n2 < nVoxelNum2 ; ++ n2) { +// +// // Se i voxel sono adiacenti proseguo +// if ( VecTriHold[t2][n2].i == VecTriHold[t1][n1].i - 1 || +// VecTriHold[t2][n2].i == VecTriHold[t1][n1].i + 1 || +// VecTriHold[t2][n2].j == VecTriHold[t1][n1].j - 1 || +// VecTriHold[t2][n2].j == VecTriHold[t1][n1].j + 1 || +// VecTriHold[t2][n2].k == VecTriHold[t1][n1].k - 1 || +// VecTriHold[t2][n2].k == VecTriHold[t1][n1].k + 1) { +// +// // Numero delle componenti connesse nei due voxel +// int nNumCompo1 = int( VecTriHold[t1][n1].ptCompoVert.size()) ; +// int nNumCompo2 = int( VecTriHold[t2][n2].ptCompoVert.size()) ; +// +// int nCompo1 = 0 ; +// +// // Ciclo sulle componenti +// for ( ; nCompo1 < nNumCompo1 ; ++ nCompo1) { +// +// int nCompo2 = ( t1 == t2 && n1 == n2 ? nCompo1 + 1 : 0) ; +// +// for ( ; nCompo2 < nNumCompo2 ; ++ nCompo2) { +// +// // Numero di triangoli per le componenti connesse +// int nTriNum1 = int( VecTriHold[t1][n1].vCompoTria[nCompo1].size()) ; +// int nTriNum2 = int( VecTriHold[t2][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 = VecTriHold[t1][n1].vCompoTria[nCompo1][nTri1].GetP( nVert1) ; +// Point3d ptP2 = VecTriHold[t2][n2].vCompoTria[nCompo2][nTri2].GetP( nVert2) ; +// +// if ( SqDist( ptP1, ptP2) < dSqEps) { +// +// Point3d ptVert1 = VecTriHold[t1][n1].ptCompoVert[nCompo1] ; +// Point3d ptVert2 = VecTriHold[t2][n2].ptCompoVert[nCompo2] ; +// +// if ( ! ( AreSamePointApprox( ptP1, ptVert1) || +// AreSamePointApprox( ptP2, ptVert2))) { +// +// 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) { +// +// int nProd = ( SharedIndex[2] - SharedIndex[0]) * ( SharedIndex[3] - SharedIndex[1]) ; +// +// // --- +// if ( nProd != 1 && nProd != - 2 && nProd != 4) { +// +// VecTriHold[t1][n1].vCompoTria[nCompo1][nTri1].SetP( SharedIndex[0], +// VecTriHold[t2][n2].ptCompoVert[nCompo2]) ; +// VecTriHold[t2][n2].vCompoTria[nCompo2][nTri2].SetP( SharedIndex[3], +// VecTriHold[t1][n1].ptCompoVert[nCompo1]) ; +// +// VecTriHold[t1][n1].vCompoTria[nCompo1][nTri1].Validate( true) ; +// VecTriHold[t2][n2].vCompoTria[nCompo2][nTri2].Validate( true) ; +// +// bModified = true ; +// break ; +// } +// } +// } +// +// if ( bModified) +// break ; +// } +// } +// } +// } +// } +// } +// } +// } +// } +// +// return true ; +//} + //---------------------------------------------------------------------------- bool -VolZmap::FlipEdges( std::vector& VecTriHold) const +VolZmap::FlipEdges( TriaMatrix& VecTriHold) const { double dSqEps = EPS_SMALL * EPS_SMALL ; @@ -1412,13 +1915,13 @@ VolZmap::FlipEdges( std::vector& VecTriHold) const // Determino i,j,k del primo blocco int nIJK1[3] ; - GetBlockIJK( int( t1), nIJK1) ; + GetBlockIJKFromN( int( t1), nIJK1) ; for ( size_t t2 = t1 ; t2 < m_nNumBlock ; ++ t2) { // Determino i,j,k del secondo blocco int nIJK2[3] ; - GetBlockIJK( int( t2), nIJK2) ; + GetBlockIJKFromN( int( t2), nIJK2) ; // Se i blocchi sono adiacenti o coincidenti proseguo if ( ( nIJK2[0] >= nIJK1[0] - 1 && nIJK2[0] <= nIJK1[0] + 1) || @@ -1437,10 +1940,10 @@ VolZmap::FlipEdges( std::vector& VecTriHold) const for ( int n1 = nSt1 ; n1 < nEn1 ; ++ n1) { // Determino estremi del ciclo sui voxel interno - int nSt2 = ( t1 == t2 ? nSt1 : 0) ; + int nSt2 = ( t1 == t2 ? n1 + 1 : 0) ; int nEn2 = nVoxelNum2 ; - for ( int n2 = n1 ; n2 < nVoxelNum2 ; ++ n2) { + for ( int n2 = nSt2 ; n2 < nVoxelNum2 ; ++ n2) { // Se i voxel sono adiacenti proseguo if ( VecTriHold[t2][n2].i == VecTriHold[t1][n1].i - 1 || @@ -1518,7 +2021,9 @@ VolZmap::FlipEdges( std::vector& VecTriHold) const VecTriHold[t1][n1].vCompoTria[nCompo1][nTri1].Validate( true) ; VecTriHold[t2][n2].vCompoTria[nCompo2][nTri2].Validate( true) ; - + if ( t1 != t2) { + int bau = 1 ; + } bModified = true ; break ; } @@ -1540,6 +2045,684 @@ VolZmap::FlipEdges( std::vector& VecTriHold) const return true ; } +//---------------------------------------------------------------------------- +//bool +//VolZmap::FlipEdges( std::vector& VecTriHold) const +//{ +// double dSqEps = EPS_SMALL * EPS_SMALL ; +// size_t nEnd = VecTriHold.size() ; //m_nNumBlock ; +// // Ciclo sui blocchi +// for ( size_t t = 0 ; t < nEnd ; ++ t) { +// +// // Numero di voxel in cui si presentano sharp feature +// int nVoxelNum = int( VecTriHold[t].size()) ; +// +// // Determino estremi del ciclo sui voxel esterno +// int nSt1 = 0 ; +// int nEn1 = nVoxelNum - 1 ; +// +// // Ciclo su tali voxel +// for ( int n1 = nSt1 ; n1 < nEn1 ; ++ n1) { +// +// // Determino estremi del ciclo sui voxel interno +// int nSt2 = nSt1 + 1 ; +// int nEn2 = nVoxelNum ; +// +// for ( int n2 = nSt2 ; n2 < nEn2 ; ++ n2) { +// +// // Se i voxel sono adiacenti proseguo +// if ( VecTriHold[t][n2].i == VecTriHold[t][n1].i + 1 || +// VecTriHold[t][n2].j == VecTriHold[t][n1].j + 1 || +// VecTriHold[t][n2].k == VecTriHold[t][n1].k + 1) { +// +// // Numero delle componenti connesse nei due voxel +// int nNumCompo1 = int( VecTriHold[t][n1].ptCompoVert.size()) ; +// int nNumCompo2 = int( VecTriHold[t][n2].ptCompoVert.size()) ; +// +// int nCompo1 = 0 ; +// +// // Ciclo sulle componenti +// for ( ; nCompo1 < nNumCompo1 ; ++ nCompo1) { +// +// int nCompo2 = 0 ; +// +// for ( ; nCompo2 < nNumCompo2 ; ++ nCompo2) { +// +// // Numero di triangoli per le componenti connesse +// int nTriNum1 = int( VecTriHold[t][n1].vCompoTria[nCompo1].size()) ; +// int nTriNum2 = int( VecTriHold[t][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 = VecTriHold[t][n1].vCompoTria[nCompo1][nTri1].GetP( nVert1) ; +// Point3d ptP2 = VecTriHold[t][n2].vCompoTria[nCompo2][nTri2].GetP( nVert2) ; +// +// if ( SqDist( ptP1, ptP2) < dSqEps) { +// +// Point3d ptVert1 = VecTriHold[t][n1].ptCompoVert[nCompo1] ; +// Point3d ptVert2 = VecTriHold[t][n2].ptCompoVert[nCompo2] ; +// +// if ( ! ( AreSamePointApprox( ptP1, ptVert1) || +// AreSamePointApprox( ptP2, ptVert2))) { +// +// 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) { +// +// int nProd = ( SharedIndex[2] - SharedIndex[0]) * ( SharedIndex[3] - SharedIndex[1]) ; +// +// // --- +// if ( nProd != 1 && nProd != - 2 && nProd != 4) { +// +// Triangle3d TriVox1 = VecTriHold[t][n1].vCompoTria[nCompo1][nTri1] ; +// Triangle3d TriVox2 = VecTriHold[t][n2].vCompoTria[nCompo2][nTri2] ; +// +// TriVox1.SetP( SharedIndex[0], VecTriHold[t][n2].ptCompoVert[nCompo2]) ; +// TriVox2.SetP( SharedIndex[3], VecTriHold[t][n1].ptCompoVert[nCompo1]) ; +// +// TriVox1.Validate( true) ; +// TriVox2.Validate( true) ; +// // questo dentro il commento è un'assurdità +// /*bool bContinue = true ; +// +// for ( int nControl = 0 ; nControl < nTriNum1 ; ++ nControl) { +// +// Vector3d vtControl = VecTriHold[t][n1].vCompoTria[nCompo1][nControl].GetN() ; +// Vector3d vtVox = TriVox1.GetN() ; +// +// if ( vtControl * vtVox < - 0.9) { +// bContinue = false ; +// break ; +// } +// } +// +// if ( bContinue) { +// +// for ( int nControl = 0 ; nControl < nTriNum2 ; ++ nControl) { +// +// Vector3d vtControl = VecTriHold[t][n2].vCompoTria[nCompo2][nControl].GetN() ; +// Vector3d vtVox = TriVox2.GetN() ; +// +// if ( vtControl * vtVox < - 0.9) { +// bContinue = false ; +// break ; +// } +// } +// } +// +// if ( bContinue) {*/ +// +// VecTriHold[t][n1].vCompoTria[nCompo1][nTri1].SetP( SharedIndex[0], +// VecTriHold[t][n2].ptCompoVert[nCompo2]) ; +// VecTriHold[t][n2].vCompoTria[nCompo2][nTri2].SetP( SharedIndex[3], +// VecTriHold[t][n1].ptCompoVert[nCompo1]) ; +// +// VecTriHold[t][n1].vCompoTria[nCompo1][nTri1].Validate( true) ; +// VecTriHold[t][n2].vCompoTria[nCompo2][nTri2].Validate( true) ; +// +// bModified = true ; +// break ; +// /*}*/ +// } +// } +// } +// +// if ( bModified) +// break ; +// } +// } +// } +// } +// } +// } +// } +// return true ; +//} + +//---------------------------------------------------------------------------- +bool +VolZmap::FlipEdgesII( TriHolder& TriHold) const +{ + double dSqEps = EPS_SMALL * EPS_SMALL ; + + // Numero di voxel in cui si presentano sharp feature + int nVoxelNum = int( TriHold.size()) ; + + // Determino estremi del ciclo sui voxel esterno + int nSt1 = 0 ; + int nEn1 = nVoxelNum - 1 ; + + // Ciclo su tali voxel + for ( int n1 = nSt1 ; n1 < nEn1 ; ++ n1) { + + // Determino estremi del ciclo sui voxel interno + int nSt2 = nSt1 + 1 ; + int nEn2 = nVoxelNum ; + + for ( int n2 = nSt2 ; n2 < nEn2 ; ++ n2) { + + // Se i voxel sono adiacenti proseguo + if ( TriHold[n2].i == TriHold[n1].i + 1 || + TriHold[n2].j == TriHold[n1].j + 1 || + 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 = 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 ( SqDist( ptP1, ptP2) < dSqEps) { + + Point3d ptVert1 = TriHold[n1].ptCompoVert[nCompo1] ; + Point3d ptVert2 = TriHold[n2].ptCompoVert[nCompo2] ; + + if ( ! ( AreSamePointApprox( ptP1, ptVert1) || + AreSamePointApprox( ptP2, ptVert2))) { + + 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) { + + int nProd = ( SharedIndex[2] - SharedIndex[0]) * ( SharedIndex[3] - SharedIndex[1]) ; + + // --- + if ( nProd != 1 && nProd != - 2 && nProd != 4) { + + Triangle3d TriVox1 = TriHold[n1].vCompoTria[nCompo1][nTri1] ; + Triangle3d TriVox2 = TriHold[n2].vCompoTria[nCompo2][nTri2] ; + + TriVox1.SetP( SharedIndex[0], TriHold[n2].ptCompoVert[nCompo2]) ; + TriVox2.SetP( SharedIndex[3], TriHold[n1].ptCompoVert[nCompo1]) ; + + TriVox1.Validate( true) ; + TriVox2.Validate( true) ; + // questo dentro il commento è un'assurdità + /*bool bContinue = true ; + + for ( int nControl = 0 ; nControl < nTriNum1 ; ++ nControl) { + + Vector3d vtControl = VecTriHold[t][n1].vCompoTria[nCompo1][nControl].GetN() ; + Vector3d vtVox = TriVox1.GetN() ; + + if ( vtControl * vtVox < - 0.9) { + bContinue = false ; + break ; + } + } + + if ( bContinue) { + + for ( int nControl = 0 ; nControl < nTriNum2 ; ++ nControl) { + + Vector3d vtControl = VecTriHold[t][n2].vCompoTria[nCompo2][nControl].GetN() ; + Vector3d vtVox = TriVox2.GetN() ; + + if ( vtControl * vtVox < - 0.9) { + bContinue = false ; + break ; + } + } + } + + if ( bContinue) {*/ + + 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 +VolZmap::FlipEdgesBB( TriaMatrix& InterTria) const { + + double dSqEps = EPS_SMALL * EPS_SMALL ; + + // Numero di blocchi + size_t nBlocksNum = InterTria.size() ; + + // ciclo sui blocchi + for ( size_t tFB = 0 ; tFB < nBlocksNum - 1 ; ++ tFB) { + for ( size_t tLB = tFB + 1 ; tLB < nBlocksNum ; ++ tLB) { + + int nFBijk[3] ; + int nLBijk[3] ; + + GetBlockIJKFromN( int( tFB), nFBijk) ; + GetBlockIJKFromN( int( tLB), nLBijk) ; + + // Se i blocchi non sono adiacenti salto l'iterazione + if ( ! ( abs( nFBijk[0] - nLBijk[0]) <= 1 && + abs( nFBijk[1] - nLBijk[1]) <= 1 && + abs( nFBijk[2] - nLBijk[2]) <= 1)) + + continue ; + + // Numero di voxel nei blocchi correnti + size_t nVoxelNumFB = InterTria[tFB].size() ; + size_t nVoxelNumLB = InterTria[tLB].size() ; + + // Ciclo sui voxel dei due blocchi + for ( size_t tVFB = 0 ; tVFB < nVoxelNumFB ; ++ tVFB) { + for ( size_t tVLB = 0 ; tVLB < nVoxelNumLB ; ++ tVLB) { + + // Se i voxel non sono adiacenti salto l'iterazione + if ( ! ( abs( InterTria[tFB][tVFB].i - InterTria[tLB][tVLB].i) <= 1 && + abs( InterTria[tFB][tVFB].j - InterTria[tLB][tVLB].j) <= 1 && + abs( InterTria[tFB][tVFB].k - InterTria[tLB][tVLB].k) <= 1)) + + continue ; + + // Numero di componenti connesse dei voxel + size_t nCompoVFBNum = InterTria[tFB][tVFB].ptCompoVert.size() ; + size_t nCompoVLBNum = InterTria[tLB][tVLB].ptCompoVert.size() ; + + // Ciclo sulle componenti connesse + for ( size_t tCmpF = 0 ; tCmpF < nCompoVFBNum ; ++ tCmpF) { + for ( size_t tCmpL = 0 ; tCmpL < nCompoVLBNum ; ++ tCmpL) { + + // Numero di triangoli delle componenti connesse + size_t nTriFBNum = InterTria[tFB][tVFB].vCompoTria[tCmpF].size() ; + size_t nTriLBNum = InterTria[tLB][tVLB].vCompoTria[tCmpL].size() ; + + for ( size_t tTriFB = 0 ; tTriFB < nTriFBNum ; ++ tTriFB) { + + bool bModified = false ; + + for ( size_t tTriLB = 0 ; tTriLB < nTriLBNum ; ++ tTriLB) { + + INTVECTOR SharedIndex ; + + for ( int nVertF = 0 ; nVertF < 3 ; ++ nVertF) { + for ( int nVertL = 0 ; nVertL < 3 ; ++ nVertL) { + + Point3d ptPF = InterTria[tFB][tVFB].vCompoTria[tCmpF][tTriFB].GetP( nVertF) ; + Point3d ptPL = InterTria[tLB][tVLB].vCompoTria[tCmpL][tTriLB].GetP( nVertL) ; + + if ( SqDist( ptPF, ptPL) < dSqEps) { + + Point3d ptVertF = InterTria[tFB][tVFB].ptCompoVert[tCmpF] ; + Point3d ptVertL = InterTria[tLB][tVLB].ptCompoVert[tCmpL] ; + + if ( ! ( AreSamePointApprox( ptPF, ptVertF) || + AreSamePointApprox( ptPL, ptVertL))) { + + SharedIndex.emplace_back( nVertF) ; + SharedIndex.emplace_back( nVertL) ; + } + } + + if ( SharedIndex.size() > 2) + break ; + } + + if ( SharedIndex.size() > 2) + break ; + } + + // Si deve operare la modifica dei triangoli + if ( SharedIndex.size() > 2) { + + int nProd = ( SharedIndex[2] - SharedIndex[0]) * ( SharedIndex[3] - SharedIndex[1]) ; + + + if ( nProd != 1 && nProd != - 2 && nProd != 4) { + + Triangle3d TriVoxF = InterTria[tFB][tVFB].vCompoTria[tCmpF][tTriFB] ; + Triangle3d TriVoxL = InterTria[tLB][tVLB].vCompoTria[tCmpL][tTriLB] ; + + TriVoxF.SetP( SharedIndex[0], InterTria[tLB][tVLB].ptCompoVert[tCmpL]) ; + TriVoxL.SetP( SharedIndex[3], InterTria[tFB][tVFB].ptCompoVert[tCmpF]) ; + + TriVoxF.Validate( true) ; + TriVoxL.Validate( true) ; + + + 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) ; + + bModified = true ; + break ; + + } + } + } + if ( bModified) + break ; + } + } + } + } + } + } + } + return true ; +} + +//---------------------------------------------------------------------------- +bool +VolZmap::FlipEdgesIB( TriHolder& InnerVox, size_t nBlock) const +{ + double dSqEps = EPS_SMALL * EPS_SMALL ; + + size_t nNumI = InnerVox.size() ; + size_t nNumB = m_InterBlockTria[nBlock].size() ; + + for ( size_t tB = 0 ; tB < nNumB ; ++ tB) { + for ( size_t tI = 0 ; tI < nNumI ; ++ tI) { + + if ( abs( InnerVox[tI].i - m_InterBlockTria[nBlock][tB].i) <= 1 && + abs( InnerVox[tI].j - m_InterBlockTria[nBlock][tB].j) <= 1 && + abs( InnerVox[tI].k - m_InterBlockTria[nBlock][tB].k) <= 1) { + + // Numero delle componenti connessVecTre nei due voxel + int nNumCompoB = int( m_InterBlockTria[nBlock][tB].ptCompoVert.size()) ; + int nNumCompoI = int( InnerVox[tI].ptCompoVert.size()) ; + + int nCompoB = 0 ; + + // Ciclo sulle componenti + for ( ; nCompoB < nNumCompoB ; ++ nCompoB) { + + int nCompoI = 0 ; + + for ( ; nCompoI < nNumCompoI ; ++ nCompoI) { + + // Numero di triangoli per le componenti connesse + int nTriNumB = int( m_InterBlockTria[nBlock][tB].vCompoTria[nCompoB].size()) ; + int nTriNumI = int( InnerVox[tI].vCompoTria[nCompoI].size()) ; + + for ( int nTriB = 0 ; nTriB < nTriNumB ; ++ nTriB) { + + bool bModified = false ; + + for ( int nTriI = 0 ; nTriI < nTriNumI ; ++ nTriI) { + + INTVECTOR SharedIndex ; + + for ( int nVertB = 0 ; nVertB < 3 ; ++ nVertB) { + + for ( int nVertI = 0 ; nVertI < 3 ; ++ nVertI) { + + Point3d ptPB = m_InterBlockTria[nBlock][tB].vCompoTria[nCompoB][nTriB].GetP( nVertB) ; + Point3d ptPI = InnerVox[tI].vCompoTria[nCompoI][nTriI].GetP( nVertI) ; + + if ( SqDist( ptPB, ptPI) < dSqEps) { + + Point3d ptVertB = m_InterBlockTria[nBlock][tB].ptCompoVert[nCompoB] ; + Point3d ptVertI = InnerVox[tI].ptCompoVert[nCompoI] ; + + if ( ! ( AreSamePointApprox( ptPB, ptVertB) || + AreSamePointApprox( ptPI, ptVertI))) { + + SharedIndex.emplace_back( nVertB) ; + SharedIndex.emplace_back( nVertI) ; + } + } + + if ( SharedIndex.size() > 2) + break ; + } + + if ( SharedIndex.size() > 2) + break ; + } + + // Si deve operare la modifica dei triangoli + if ( SharedIndex.size() > 2) { + + int nProd = ( SharedIndex[2] - SharedIndex[0]) * ( SharedIndex[3] - SharedIndex[1]) ; + + // --- + if ( nProd != 1 && nProd != - 2 && nProd != 4) { + + Triangle3d TriVoxB = m_InterBlockTria[nBlock][tB].vCompoTria[nCompoB][nTriB] ; + Triangle3d TriVoxI = InnerVox[tI].vCompoTria[nCompoI][nTriI] ; + + TriVoxB.SetP( SharedIndex[0], InnerVox[tI].ptCompoVert[nCompoI]) ; + TriVoxI.SetP( SharedIndex[3], m_InterBlockTria[nBlock][tB].ptCompoVert[nCompoB]) ; + + TriVoxB.Validate( true) ; + TriVoxI.Validate( true) ; + + + m_InterBlockTria[nBlock][tB].vCompoTria[nCompoB][nTriB].SetP( SharedIndex[0], + InnerVox[tI].ptCompoVert[nCompoI]) ; + InnerVox[tI].vCompoTria[nCompoI][nTriI].SetP( SharedIndex[3], + m_InterBlockTria[nBlock][tB].ptCompoVert[nCompoB]) ; + + m_InterBlockTria[nBlock][tB].vCompoTria[nCompoB][nTriB].Validate( true) ; + InnerVox[tI].vCompoTria[nCompoI][nTriI].Validate( true) ; + + bModified = true ; + break ; + } + } + } + + if ( bModified) + break ; + } + } + } + + } + } + } + return true ; + +} + +//---------------------------------------------------------------------------- +bool +VolZmap::FlipEdgesBB( size_t nCurBlock, size_t nAdjBlock) const +{ + double dSqEps = EPS_SMALL * EPS_SMALL ; + + return true ; + +} + +//---------------------------------------------------------------------------- +bool +VolZmap::FlipEdgesIB( TriHolder& InnerVox, TriHolder& BoundaryVox) const +{ + double dSqEps = EPS_SMALL * EPS_SMALL ; + + size_t nNumI = InnerVox.size() ; + size_t nNumB = BoundaryVox.size() ; + + for ( size_t tB = 0 ; tB < nNumB ; ++ tB) { + for ( size_t tI = 0 ; tI < nNumI ; ++ tI) { + + if ( abs( InnerVox[tI].i - BoundaryVox[tB].i) <= 1 && + abs( InnerVox[tI].j - BoundaryVox[tB].j) <= 1 && + abs( InnerVox[tI].k - BoundaryVox[tB].k) <= 1) { + + // Numero delle componenti connessVecTre nei due voxel + int nNumCompoB = int( BoundaryVox[tB].ptCompoVert.size()) ; + int nNumCompoI = int( InnerVox[tI].ptCompoVert.size()) ; + + int nCompoB = 0 ; + + // Ciclo sulle componenti + for ( ; nCompoB < nNumCompoB ; ++ nCompoB) { + + int nCompoI = 0 ; + + for ( ; nCompoI < nNumCompoI ; ++ nCompoI) { + + // Numero di triangoli per le componenti connesse + int nTriNumB = int( BoundaryVox[tB].vCompoTria[nCompoB].size()) ; + int nTriNumI = int( InnerVox[tI].vCompoTria[nCompoI].size()) ; + + for ( int nTriB = 0 ; nTriB < nTriNumB ; ++ nTriB) { + + bool bModified = false ; + + for ( int nTriI = 0 ; nTriI < nTriNumI ; ++ nTriI) { + + INTVECTOR SharedIndex ; + + for ( int nVertB = 0 ; nVertB < 3 ; ++ nVertB) { + + for ( int nVertI = 0 ; nVertI < 3 ; ++ nVertI) { + + Point3d ptPB = BoundaryVox[tB].vCompoTria[nCompoB][nTriB].GetP( nVertB) ; + Point3d ptPI = InnerVox[tI].vCompoTria[nCompoI][nTriI].GetP( nVertI) ; + + if ( SqDist( ptPB, ptPI) < dSqEps) { + + Point3d ptVertB = BoundaryVox[tB].ptCompoVert[nCompoB] ; + Point3d ptVertI = InnerVox[tI].ptCompoVert[nCompoI] ; + + if ( ! ( AreSamePointApprox( ptPB, ptVertB) || + AreSamePointApprox( ptPI, ptVertI))) { + + SharedIndex.emplace_back( nVertB) ; + SharedIndex.emplace_back( nVertI) ; + } + } + + if ( SharedIndex.size() > 2) + break ; + } + + if ( SharedIndex.size() > 2) + break ; + } + + // Si deve operare la modifica dei triangoli + if ( SharedIndex.size() > 2) { + + int nProd = ( SharedIndex[2] - SharedIndex[0]) * ( SharedIndex[3] - SharedIndex[1]) ; + + // --- + if ( nProd != 1 && nProd != - 2 && nProd != 4) { + + Triangle3d TriVoxB = BoundaryVox[tB].vCompoTria[nCompoB][nTriB] ; + Triangle3d TriVoxI = InnerVox[tI].vCompoTria[nCompoI][nTriI] ; + + TriVoxB.SetP( SharedIndex[0], InnerVox[tI].ptCompoVert[nCompoI]) ; + TriVoxI.SetP( SharedIndex[3], BoundaryVox[tB].ptCompoVert[nCompoB]) ; + + TriVoxB.Validate( true) ; + TriVoxI.Validate( true) ; + + + BoundaryVox[tB].vCompoTria[nCompoB][nTriB].SetP( SharedIndex[0], + InnerVox[tI].ptCompoVert[nCompoI]) ; + InnerVox[tI].vCompoTria[nCompoI][nTriI].SetP( SharedIndex[3], + BoundaryVox[tB].ptCompoVert[nCompoB]) ; + + BoundaryVox[tB].vCompoTria[nCompoB][nTriB].Validate( true) ; + InnerVox[tI].vCompoTria[nCompoI][nTriI].Validate( true) ; + + bModified = true ; + break ; + } + } + } + + if ( bModified) + break ; + } + } + } + + } + } + } + return true ; +} + //---------------------------------------------------------------------------- bool VolZmap::IsThereMat( int nI, int nJ, int nK) const @@ -1911,7 +3094,7 @@ VolZmap::IntersPos( int nVec1[], int nVec2[], bool bFirstCorner, Point3d& ptInt, //---------------------------------------------------------------------------- bool -VolZmap::GetBlockIJK( int nBlock, int nIJK[]) const +VolZmap::GetBlockIJKFromN( int nBlock, int nIJK[]) const { // Controllo sulla validità del blocco if ( nBlock < 0 || nBlock >= int( m_nNumBlock)) @@ -1925,6 +3108,22 @@ VolZmap::GetBlockIJK( int nBlock, int nIJK[]) const return true ; } +//---------------------------------------------------------------------------- +bool +VolZmap::GetBlockNFromIJK( int nIJK[], int& nBlock) const +{ + // Controllo sulla validità degli indici i, j, k del blocco + if ( nIJK[0] < 0 || nIJK[0] >= int( m_nFracLin[0]) || + nIJK[1] < 0 || nIJK[1] >= int( m_nFracLin[1]) || + nIJK[2] < 0 || nIJK[2] >= int( m_nFracLin[2])) + return false ; + + // Determino il blocco + nBlock = m_nFracLin[0] * m_nFracLin[1] * nIJK[2] + m_nFracLin[0] * nIJK[1] + nIJK[0] ; + + return true ; +} + //---------------------------------------------------------------------------- bool VolZmap::GetBlockLimitsIJK( const int nIJK[], int nLimits[]) const @@ -2008,6 +3207,339 @@ VolZmap::GetPointVoxel( Point3d& ptP, int& nVoxI, int& nVoxJ, int& nVoxK) const ( nVoxK >= -1 && nVoxK < int( m_nNy[1])) ; } +//---------------------------------------------------------------------------- +bool +VolZmap::GetVoxelBlockIJK( const int nVoxIJK[], int nBlockIJK[]) const +{ + // Controllo sull'ammissibilità del voxel + if ( nVoxIJK[0] <= -2 || nVoxIJK[0] >= int( m_nNx[0]) || + nVoxIJK[1] <= -2 || nVoxIJK[1] >= int( m_nNy[0]) || + nVoxIJK[2] <= -2 || nVoxIJK[2] >= int( m_nNy[1])) + return false ; + + // Divisioni intere + int nIntRatio0 = nVoxIJK[0] / m_nDexNumPBlock ; + int nIntRatio1 = nVoxIJK[1] / m_nDexNumPBlock ; + int nIntRatio2 = nVoxIJK[2] / m_nDexNumPBlock ; + + // Resti delle divisioni intere + /*int nMod0 = nVoxIJK[0] % m_nDexNumPBlock ; + int nMod1 = nVoxIJK[1] % m_nDexNumPBlock ; + int nMod2 = nVoxIJK[2] % m_nDexNumPBlock ;*/ + + // Calcolo indici del blocco + nBlockIJK[0] = ( nVoxIJK[0] == -1 ? 0 : max( 0, nIntRatio0 - ( nIntRatio0 == m_nFracLin[0] ? 1 : 0))) ; + nBlockIJK[1] = ( nVoxIJK[1] == -1 ? 0 : max( 0, nIntRatio1 - ( nIntRatio1 == m_nFracLin[1] ? 1 : 0))) ; + nBlockIJK[2] = ( nVoxIJK[2] == -1 ? 0 : max( 0, nIntRatio2 - ( nIntRatio2 == m_nFracLin[2] ? 1 : 0))) ; + + return true ; +} + +//---------------------------------------------------------------------------- +bool +VolZmap::GetAdjBlockToBlock( int nBlockN, int nDeltaI, int nDeltaJ, int nDeltaK, int& nAdjBlockN) const +{ + // Test sulla validità degli incrementi su i,j,k + if ( nDeltaI < - 1 || nDeltaI > 1 || + nDeltaJ < - 1 || nDeltaJ > 1 || + nDeltaK < - 1 || nDeltaK > 1) + return false ; + + // Determino blocco adiacente + nAdjBlockN = nBlockN ; + nAdjBlockN += nDeltaI ; + nAdjBlockN += nDeltaJ * m_nFracLin[0] ; + nAdjBlockN += nDeltaK * m_nFracLin[0] * m_nFracLin[1] ; + + // Se il blocco adiacente esiste restituisco vero, altrimenti falso. + return ( nAdjBlockN > -1 && nAdjBlockN < int( m_nNumBlock)) ; +} + +//---------------------------------------------------------------------------- +bool +VolZmap::IsAVoxelOnBoundary( const int nLimits[], const int nIJK[], bool bType) const +{ + // Test sulla validità dei limiti + if ( nLimits[0] < -1 || nLimits[0] > int( m_nNx[0]) || + nLimits[1] < -1 || nLimits[1] > int( m_nNx[0]) || + nLimits[2] < -1 || nLimits[2] > int( m_nNy[0]) || + nLimits[3] < -1 || nLimits[3] > int( m_nNy[0]) || + nLimits[4] < -1 || nLimits[4] > int( m_nNy[1]) || + nLimits[5] < -1 || nLimits[5] > int( m_nNy[1])) + return false ; + + // Controllo sull'ammissibilità del voxel + if ( nIJK[0] <= -2 || nIJK[0] >= int( m_nNx[0]) || + nIJK[1] <= -2 || nIJK[1] >= int( m_nNx[1]) || + nIJK[2] <= -2 || nIJK[2] >= int( m_nNy[1])) + return false ; + + // Se bType è vero cerchiamo i voxel che + // confinano con quelli di atri blocchi, + if ( bType) { + + 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} ; + + int nAdjBlockIJK[3] ; + + // Determino (se esiste) il blocco in cui cade il voxel adiacente. + 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()) + return false ; + + // Se il voxel sta sulla frontiera contino, + if ( IsAVoxelOnBoundary( nBlockLimits, nVoxIJK, true)) { + + double dSqEps = EPS_SMALL * EPS_SMALL ; + + Point3d ptFirstGrPt, ptSecondGrPt ; + + int nCount = 0 ; + + if ( SqDist( trTria.GetP( 0), ptVert) > dSqEps) { + + ptFirstGrPt = trTria.GetP( 0) ; + nCount ++ ; + } + + if ( SqDist( trTria.GetP( 1), ptVert) > dSqEps) { + + if ( nCount == 0) { + + ptFirstGrPt = trTria.GetP( 1) ; + nCount ++ ; + } + else if ( nCount == 1) { + + ptSecondGrPt = trTria.GetP( 1) ; + nCount ++ ; + } + } + + if ( SqDist( trTria.GetP( 2), ptVert) > dSqEps) { + + if ( nCount == 1) { + + ptSecondGrPt = trTria.GetP( 2) ; + nCount ++ ; + } + } + + /* Altro modo: + std::vector vPt ; + for ( int nC = 0 ; nC < 3 ; ++ nC) + if ( SqDist( trTria.GetP( 0), ptVert) > dSqEps) + vPt.emplace_back( trTria.GetP( 0)) ; + if ( vPt.size() == 2) + ......*/ + + + 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 true ; + + /*Altro modo + for ( int nC = 0 ; nC < 3 ; ++ nC) { + + if ( nVoxIJK[nC] == nBlockLimits[2*nC]) { + + double dGrid = ( nBlockLimits[2*nC] + 0.5) * m_dStep ; + + if ( abs( ptFirstGrPt.v[nC] - dGrid) < EPS_SMALL && + abs( ptSecondGrPt.v[nC] - dGrid) < EPS_SMALL) + + return true ; + } + + if ( nVoxIJK[nC] + 1 == nBlockLimits[2*nC+1]) { + + double dGrid = ( nBlockLimits[2*nC+1] + 0.5) * m_dStep ; + + if ( abs( ptFirstGrPt.v[nC] - dGrid) < EPS_SMALL && + abs( ptSecondGrPt.v[nC] - dGrid) < EPS_SMALL) + + return true ; + } + }*/ + } + // altrimenti non può starci il triangolo + else + return false ; +} + +//---------------------------------------------------------------------------- +bool +VolZmap::IsATriangleOnBorder( const Triangle3d& trTria, const Point3d& ptVert, + const int nBlockLimits[], const int nVoxIJK[], bool bBorderBox) const +{ + // Se il voxel è di frontiera continuo, + if ( bBorderBox) { + + // Se il triangolo ha area nulla non proseguiamo + if ( trTria.GetArea() < EPS_SMALL) + return false ; + + // Determino i punti del triangolo sulla griglia + double dSqEps = EPS_SMALL * EPS_SMALL ; + + int nNotVert[2] ; + int nCount = 0 ; + + for ( int nV = 0 ; nV < 3 ; ++ nV) { + + if ( SqDist( trTria.GetP( nV), ptVert) > dSqEps) { + + if ( nCount == 0) { + + nNotVert[nCount] = nV ; + } + else if ( nCount == 1) { + + nNotVert[nCount] = nV ; + } + + nCount ++ ; + } + } + + // Se nCount != 2 deve esserci un errore + if ( nCount != 2) + return false ; + + // Punti del triangolo sulla griglia + Point3d ptFirstGrPt = trTria.GetP( nNotVert[0]) ; + Point3d ptSecondGrPt = trTria.GetP( nNotVert[1]) ; + + // Verifico se tali punti sono sula griglia + for ( int nC = 0 ; nC < 3 ; ++ nC) { + + if ( nVoxIJK[nC] == nBlockLimits[2*nC]) { + + double dGrid = ( nBlockLimits[2*nC] + 0.5) * m_dStep ; + + if ( abs( ptFirstGrPt.v[nC] - dGrid) < EPS_SMALL && + abs( ptSecondGrPt.v[nC] - dGrid) < EPS_SMALL) + + return true ; + } + + if ( nVoxIJK[nC] + 1 == nBlockLimits[2*nC+1]) { + + double dGrid = ( nBlockLimits[2*nC+1] + 0.5) * m_dStep ; + + if ( abs( ptFirstGrPt.v[nC] - dGrid) < EPS_SMALL && + abs( ptSecondGrPt.v[nC] - dGrid) < EPS_SMALL) + + return true ; + } + } + return false ; + } + // altrimenti non ha senso continuare + else + return false ; +} + diff --git a/VolTriZmapVolume.cpp b/VolTriZmapVolume.cpp index 5a8358f..e223569 100644 --- a/VolTriZmapVolume.cpp +++ b/VolTriZmapVolume.cpp @@ -157,17 +157,17 @@ VolZmap::SubtractIntervals( unsigned int nGrid, unsigned int nI, unsigned int nJ int nMaxZBlock = min( int( m_nFracLin[2] - 1), nMaxK / int( m_nDexNumPBlock)) ; for ( int k = nMinZBlock ; k <= nMaxZBlock ; ++ k) { - for ( int nIndX = nXBlock - 1 ; nIndX <= nXBlock + 1 ; nIndX ++) { - for ( int nIndY = nYBlock - 1 ; nIndY <= nYBlock + 1 ; nIndY ++) { + /* for ( int nIndX = nXBlock - 1 ; nIndX <= nXBlock + 1 ; nIndX ++) { + for ( int nIndY = nYBlock - 1 ; nIndY <= nYBlock + 1 ; nIndY ++) {*/ - if ( nIndX >= 0 && nIndX < int( m_nFracLin[0]) && - nIndY >= 0 && nIndY < int( m_nFracLin[1])) { + /*if ( nIndX >= 0 && nIndX < int( m_nFracLin[0]) && + nIndY >= 0 && nIndY < int( m_nFracLin[1])) {*/ - int nBlockNum = k * nLayerBlock + nIndY * m_nFracLin[0] + nIndX ; - m_BlockToUpdate[nBlockNum] = true ; - } + int nBlockNum = k * nLayerBlock + nYBlock * m_nFracLin[0] + nXBlock ; + m_BlockToUpdate[nBlockNum] = true ; + /* } } - } + }*/ } } else if ( nGrid == 1) { @@ -182,17 +182,17 @@ VolZmap::SubtractIntervals( unsigned int nGrid, unsigned int nI, unsigned int nJ int nMaxXBlock = min( int( m_nFracLin[0] - 1), nMaxI / int( m_nDexNumPBlock)) ; for ( int k = nMinXBlock ; k <= nMaxXBlock ; ++ k) { - for( int nIndY = nYBlock - 1 ; nIndY <= nYBlock + 1 ; ++ nIndY) { + /* for( int nIndY = nYBlock - 1 ; nIndY <= nYBlock + 1 ; ++ nIndY) { for ( int nIndZ = nZBlock - 1 ; nIndZ <= nZBlock + 1 ; ++ nIndZ) { if ( nIndY >= 0 && nIndY < int( m_nFracLin[1]) && - nIndZ >= 0 && nIndZ < int( m_nFracLin[2])) { + nIndZ >= 0 && nIndZ < int( m_nFracLin[2])) {*/ - int nBlockNum = nIndZ * nLayerBlock + nIndY * m_nFracLin[0] + k ; + int nBlockNum = nZBlock * nLayerBlock + nYBlock * m_nFracLin[0] + k ; m_BlockToUpdate[nBlockNum] = true ; - } + /* } } - } + }*/ } } else if ( nGrid == 2) { @@ -207,17 +207,17 @@ VolZmap::SubtractIntervals( unsigned int nGrid, unsigned int nI, unsigned int nJ int nMaxYBlock = min( int( m_nFracLin[1] - 1), nMaxJ / int( m_nDexNumPBlock)) ; for ( int k = nMinYBlock ; k <= nMaxYBlock ; ++ k) { - for ( int nIndX = nXBlock - 1 ; nIndX <= nXBlock + 1 ; ++ nIndX) { + /* for ( int nIndX = nXBlock - 1 ; nIndX <= nXBlock + 1 ; ++ nIndX) { for ( int nIndZ = nZBlock - 1 ; nIndZ <= nZBlock + 1 ; ++ nIndZ) { if ( nIndX >= 0 && nIndX < int( m_nFracLin[0]) && nIndZ >= 0 && nIndZ < int( m_nFracLin[2])) { - - int nBlockNum = nIndZ * nLayerBlock + k * m_nFracLin[0] + nIndX ; + */ + int nBlockNum = nZBlock * nLayerBlock + k * m_nFracLin[0] + nXBlock ; m_BlockToUpdate[nBlockNum] = true ; - } + /* } } - } + }*/ } } diff --git a/VolZmap.h b/VolZmap.h index 45d103c..25d8158 100644 --- a/VolZmap.h +++ b/VolZmap.h @@ -29,7 +29,29 @@ struct TriaStruct { std::vector ptCompoVert ; std::vector vCompoTria ; } ; -typedef std::vector TriHolder ; + +// Vettore di TriaStruct con sharp fature interni a un blocco +typedef std::vector TriHolder ; + +// Vettori di TriHolder con sharp feature di frontiera: +// il primo indice individua il blocco, il secondo il voxel +typedef std::vector TriaMatrix ; + +struct IndexStruct { + + int i, j, k ; +} ; + +typedef std::vector> IndexMatrix ; +//struct BoundaryTriangle { +// +// int i, j, k ; +// int nCC ; +// +// Point3d ptVertex ; +// +// bool bInterBlockFlipping ; +//} ; //---------------------------------------------------------------------------- class VolZmap : public IVolZmap, public IGeoObjRW @@ -116,8 +138,13 @@ class VolZmap : public IVolZmap, public IGeoObjRW const Vector3d& vtZ, const Vector3d& vtNorm, TRIA3DLIST& lstTria) const ; bool MarchingCubes( TRIA3DLIST& lstTria) const ; bool MarchingCubes( int nBlock, TRIA3DLIST& lstTria) const ; - bool ExtMarchingCubes( const int nLimits[], TRIA3DLIST& lstTria, TriHolder& triHold) const ; - bool FlipEdges( std::vector& VecTriHold) const ; + bool ExtMarchingCubes( const int nLimits[], size_t nBlockNumber, TRIA3DLIST& lstTria, TriHolder& triHold) const ; + bool FlipEdges( TriaMatrix& VecTriHold) const ; + bool FlipEdgesII( TriHolder& TriHold) const ; + bool FlipEdgesBB( TriaMatrix& InterTria) const ; + bool FlipEdgesIB( TriHolder& InnerVox, size_t nBlock) const ; + bool FlipEdgesBB( size_t nCurBlock, size_t nAdjBlock) const ; + bool FlipEdgesIB( TriHolder& InnerVox, TriHolder& BoundaryVox) const ; bool IsThereMat( int nI, int nJ, int nK) const ; bool IsThereMat( const int nMatr[][3], int nNum, double & dHx, double & dHy, double & dHz) const ; bool IntersPos( int nVec1[], int nVec2[], Point3d & ptInt) const ; @@ -245,13 +272,21 @@ class VolZmap : public IVolZmap, public IGeoObjRW Point3d& ptInt1, Point3d& ptInt2, Vector3d& vtN1, Vector3d& vtN2) ; // Funzioni di gestione dei blocchi - bool GetBlockIJK( int nBlock, int nIJK[]) const ; + bool GetBlockIJKFromN( int nBlock, int nIJK[]) const ; + bool GetBlockNFromIJK( int nIJK[], int& nBlock) const ; bool GetBlockLimitsIJK( const int nIJK[], int nLimits[]) const ; + bool GetVoxelBlockIJK( const int nVoxIJK[], int nBlockIJK[]) const ; + bool GetAdjBlockToBlock( int nBlockN, int nDeltaI, int nDeltaJ, int nDeltaK, int& nAdjBlockN) const ; + bool IsAVoxelOnBoundary( const int nLimits[], const int nIJK[], bool bType) const ; + bool IsATriangleOnBorder( const Triangle3d& trTria, const Point3d& ptVert, + const int nBlockLimits[], const int nVoxIJK[]) const ; + bool IsATriangleOnBorder( const Triangle3d& trTria, const Point3d& ptVert, + const int nBlockLimits[], const int nVoxIJK[], bool bBorderBox) const ; private : enum Status { ERR = 0, OK = 1, TO_VERIFY = 2} ; static const int N_MAPS = 3 ; - static const int N_DEXBLOCK = 32 ; // 10000 ;//20 ;//32 ; + static const int N_DEXBLOCK = 32 ; // 10000 ; // 20 ; // 32 ; private : ObjGraphicsMgr m_OGrMgr ; // gestore grafica dell'oggetto @@ -290,6 +325,8 @@ class VolZmap : public IVolZmap, public IGeoObjRW double m_dRadius ; double m_dRCorner ; double m_dTipRadius ; + + mutable TriaMatrix m_InterBlockTria ; } ;