From d9d8b40a02555603bd8d5b0dccb47e4405f6f761 Mon Sep 17 00:00:00 2001 From: Dario Sassi Date: Sat, 10 Jun 2017 15:43:14 +0000 Subject: [PATCH] EgtGeomKernel 1.8f2 : - migliorie e correzioni a Zmap. --- EgtGeomKernel.rc | Bin 11710 -> 11710 bytes MC_Tables.h | 4 +- VolTriZmapCreation.cpp | 212 +++++++++++++-------- VolTriZmapGraphics.cpp | 418 ++++++++++++++--------------------------- VolTriZmapVolume.cpp | 8 +- VolZmap.cpp | 9 - VolZmap.h | 3 +- 7 files changed, 287 insertions(+), 367 deletions(-) diff --git a/EgtGeomKernel.rc b/EgtGeomKernel.rc index ba4f43363547fd0edd68f7be9407bc50b980d11f..3503d7ba569ba8b7791e7d6e921fe64dd748d6f3 100644 GIT binary patch delta 94 zcmdlNy)SyhFE&P_&A-_cnHh~HD{|{@_Trkr0u;H;XNwSVW8B;$>;>dw2zN+>g;Df- LFmBFL4&ed-Vf-5U delta 94 zcmdlNy)SyhFE&QQ&A-_cnHdcyD{|{@_Trkr0u;H;XNwSVW8B;$>;>dw2zN+>g;Df- LFmBFL4&ed-U|t&X diff --git a/MC_Tables.h b/MC_Tables.h index d9d6d68..2980262 100644 --- a/MC_Tables.h +++ b/MC_Tables.h @@ -72,10 +72,10 @@ static int TriangleTableEn[256][2][17] = { // 10 ----------------------------------------------------------- {{1, 9, 0, 2, 3, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, - {2, 3, 3, 1, 9, 0, 2, 3,11, -1, -1, -1, -1, -1, -1, -1, -1}}, + {2, 3, 3, 1, 9, 0, 2, 3, 11, -1, -1, -1, -1, -1, -1, -1, -1}}, {{1, 11, 2, 1, 9, 11, 9, 8, 11, -1, -1, -1, -1, -1, -1, -1, -1}, - {1, 5, 9, 8, 11, 2, 1,-1, -1, -1, -1, -1, -1, -1, -1, -1, -1}}, + {1, 5, 9, 8, 11, 2, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}}, {{3, 10, 1, 11, 10, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, {1, 4, 3, 11,10, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}}, diff --git a/VolTriZmapCreation.cpp b/VolTriZmapCreation.cpp index e6c113e..26a925f 100644 --- a/VolTriZmapCreation.cpp +++ b/VolTriZmapCreation.cpp @@ -27,7 +27,7 @@ using namespace std ; //---------------------------------------------------------------------------- bool -VolZmap::Create( const Point3d& ptO, double dLengthX, double dLengthY, double dLengthZ, double dPrec, bool bFlag) +VolZmap::Create( const Point3d& ptO, double dLengthX, double dLengthY, double dLengthZ, double dPrec, bool bTriDex) { // Controlli sull'ammissibilità delle dimensioni lineari del grezzo e del passo if ( dPrec < EPS_SMALL || dLengthX < EPS_SMALL || dLengthY < EPS_SMALL || dLengthZ < EPS_SMALL) @@ -37,7 +37,7 @@ VolZmap::Create( const Point3d& ptO, double dLengthX, double dLengthY, double dL m_dStep = max( dPrec, 100 * EPS_SMALL) ; // Aggiorno la dimensione della mappa 1 o 3 - m_nMapNum = ( bFlag ? 3 : 1) ; + m_nMapNum = ( bTriDex ? 3 : 1) ; // Disponendo i sistemi di riferimento in una successione, le coordinate x,y,z // di uno si ottengono da una permutazione ciclica di quelle del precedente sistema. @@ -64,7 +64,7 @@ VolZmap::Create( const Point3d& ptO, double dLengthX, double dLengthY, double dL m_nConnectedCompoCount = 1 ; // Se tridexel - if ( bFlag) { + if ( bTriDex) { m_MapFrame[1].Set( ptO, Y_AX, Z_AX, X_AX) ; m_MapFrame[2].Set( ptO, Z_AX, X_AX, Y_AX) ; @@ -78,7 +78,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)) ; } // altrimenti mono dexel else { @@ -152,7 +152,7 @@ VolZmap::Create( const Point3d& ptO, double dLengthX, double dLengthY, double dL } // Ridimensiono e setto il vettore dei blocchi a falso - m_nNumBlock = ( bFlag ? m_nFracLin[0] * m_nFracLin[1] * m_nFracLin[2] : m_nFracLin[0] * m_nFracLin[1]) ; + m_nNumBlock = m_nFracLin[0] * m_nFracLin[1] * m_nFracLin[2] ; m_BlockToUpdate.resize( m_nNumBlock) ; for ( unsigned int nCount = 0 ; nCount < m_nNumBlock ; ++ nCount) @@ -166,12 +166,12 @@ VolZmap::Create( const Point3d& ptO, double dLengthX, double dLengthY, double dL return true ; } -//---------------------------------------------------------------------------------------------------- +//---------------------------------------------------------------------------- bool -VolZmap::CreateFromFlatRegion( const ISurfFlatRegion& Surf, double dDimZ, double dPrec, bool bFlag) +VolZmap::CreateFromFlatRegion( const ISurfFlatRegion& Surf, double dDimZ, double dPrec, bool bTriDex) { // Aggiorno la dimensione della mappa 1 o 3 - m_nMapNum = ( bFlag ? 3 : 1) ; + m_nMapNum = ( bTriDex ? 3 : 1) ; // Il passo di discretizzazione non può essere inferiore a 100 * EPS_SMALL m_dStep = max( dPrec, 100 * EPS_SMALL) ; @@ -212,10 +212,10 @@ VolZmap::CreateFromFlatRegion( const ISurfFlatRegion& Surf, double dDimZ, double ( m_nNy[0] % m_nDexNumPBlock >= m_nDexNumPBlock / 2 ? 1 : 0)) ; // Numero di componenti connesse - m_nConnectedCompoCount = - 1 ; + m_nConnectedCompoCount = Surf.GetChunkCount() ; // Se Tridexel ridimensiono anche gli altri vettori - if ( bFlag) { + if ( bTriDex) { m_MapFrame[1].Set( ptMapOrig, Y_AX, Z_AX, X_AX) ; m_MapFrame[2].Set( ptMapOrig, Z_AX, X_AX, Y_AX) ; @@ -252,6 +252,9 @@ VolZmap::CreateFromFlatRegion( const ISurfFlatRegion& Surf, double dDimZ, double m_nNy[2] = 0 ; m_nDim[2] = 0 ; + + // Definisco il numero di blocchi lungo z + m_nFracLin[2] = 1 ; } // Determinazione e ridimensionamento dei dexel @@ -289,6 +292,72 @@ VolZmap::CreateFromFlatRegion( const ISurfFlatRegion& Surf, double dDimZ, double int nStartJ = Clamp( int( floor( dt1 * dLengthY / m_dStep - EPS_SMALL + 0.5)), 0, m_nNy[0] - 1) ; int nEndJ = Clamp( int( floor( dt2 * dLengthY / m_dStep + EPS_SMALL - 0.5)), 0, m_nNy[0] - 1) ; + // Punti estremi della parte di retta corrente + Point3d ptP1 = ptP0 + dt1 * dLengthY * Y_AX ; + Point3d ptP2 = ptP0 + dt2 * dLengthY * Y_AX ; + + // Determino il numero della componente connessa e le normali + int nCompo ; + Vector3d vtN1, vtN2 ; + + // Numero di componenti connesse + int nChunkNum = m_nConnectedCompoCount ; + + // Ciclo sulle componenti connesse + for ( int nChunk = 0 ; nChunk < nChunkNum ; ++ nChunk) { + + int nFind = 0 ; + + // Numero di looop della componente connessa + int nLoopNum = Surf.GetLoopCount( nChunk) ; + + // Ciclo sui loop della componente + for ( int nLoop = 0 ; nLoop < nLoopNum ; ++ nLoop) { + + ICurve * pCurve = Surf.GetLoop( nChunk, nLoop) ; + + if ( pCurve -> IsPointOn( ptP1, 10 * EPS_SMALL)) { + + double dP1 ; + Point3d ptTemp1 ; + Vector3d vtT1 ; + + pCurve -> GetParamAtPoint(ptP1, dP1) ; + pCurve -> GetPointTang( dP1, ICurve::FROM_MINUS, ptTemp1, vtT1) ; + + vtN1 = vtT1 ^ Z_AX ; + + nFind ++ ; + } + + if ( pCurve -> IsPointOn( ptP2, 10 * EPS_SMALL)) { + + double dP2 ; + Point3d ptTemp2 ; + Vector3d vtT2 ; + + pCurve -> GetParamAtPoint(ptP2, dP2) ; + pCurve -> GetPointTang( dP2, ICurve::FROM_MINUS, ptTemp2, vtT2) ; + + vtN2 = vtT2 ^ Z_AX ; + + nFind ++ ; + } + + // Se abbiamo trovato la componente + // connessa, usciamo dal ciclo. + if ( nFind == 2) { + + nCompo = nChunk + 1 ; + break ; + } + } + // Se trovata componente esco dal cilco + if ( nFind == 2) + break ; + } + + // Ridimensiono e riempio i dexel for ( int j = nStartJ ; j <= nEndJ ; ++ j) { // Determino il dexel @@ -304,12 +373,12 @@ VolZmap::CreateFromFlatRegion( const ISurfFlatRegion& Surf, double dDimZ, double m_Values[0][nPos0][0].vtMinN = - Z_AX ; m_Values[0][nPos0][0].vtMaxN = Z_AX ; m_Values[0][nPos0][0].nToolMax = 0 ; - m_Values[0][nPos0][0].nCompo = 0 ; + m_Values[0][nPos0][0].nCompo = nCompo ; } // Se tridexel riempio i singoli dexel della // griglia 2 con gli intervalli - if ( bFlag) { + if ( bTriDex) { for ( size_t a = 0 ; a < m_nNx[2] ; ++ a) { @@ -325,63 +394,19 @@ VolZmap::CreateFromFlatRegion( const ISurfFlatRegion& Surf, double dDimZ, double m_Values[2][nPos2][nCurrentSize].nToolMin = 0 ; m_Values[2][nPos2][nCurrentSize].dMax = dt2 * dLengthY ; m_Values[2][nPos2][nCurrentSize].nToolMax = 0 ; - - Point3d ptP1 = ptP0 + dt1 * dLengthY * Y_AX ;//+ ( a + 0.5) * Z_AX ; - Point3d ptP2 = ptP0 + dt2 * dLengthY * Y_AX ;//+ ( a + 0.5) * Z_AX ; - - int nChunkNum = Surf.GetChunkCount() ; - - for ( int nChunk = 0 ; nChunk < nChunkNum ; ++ nChunk) { - - int nLoopNum = Surf.GetLoopCount( nChunk) ; - - for ( int nLoop = 0 ; nLoop < nLoopNum ; ++ nLoop) { - - ICurve * pCurve = Surf.GetLoop( nChunk, nLoop) ; - - if ( pCurve -> IsPointOn( ptP1)) { - - double dP1 ; - Vector3d vtT1, vtN1 ; - - pCurve -> GetParamAtPoint(ptP1, dP1) ; - pCurve -> GetPointTang( dP1, ICurve::FROM_MINUS, ptP1, vtT1) ; - - vtN1 = vtT1 ^ Z_AX ; - - // Aggiorno normale inferiore al tratto di dexel - m_Values[2][nPos2][nCurrentSize].vtMinN = vtN1 ; - - } - - if ( pCurve -> IsPointOn( ptP2)) { - - double dP2 ; - Vector3d vtT2, vtN2 ; - - pCurve -> GetParamAtPoint(ptP2, dP2) ; - pCurve -> GetPointTang( dP2, ICurve::FROM_MINUS, ptP1, vtT2) ; - - vtN2 = vtT2 ^ Z_AX ; - - // Aggiorno normale superiore al tratto di dexel - m_Values[2][nPos2][nCurrentSize].vtMaxN = vtN2 ; - } - } - } - + m_Values[2][nPos2][nCurrentSize].vtMinN = vtN1 ; + m_Values[2][nPos2][nCurrentSize].vtMaxN = vtN2 ; + // Aggiorno il numero della componente connessa - m_Values[2][nPos2][nCurrentSize].nCompo = 0 ; + m_Values[2][nPos2][nCurrentSize].nCompo = nCompo ; } - } - - - } + } + } ///// Fine parte relativa a un segmento } } // Se tridexel resta la griglia 1 - if ( bFlag) { + if ( bTriDex) { for ( unsigned int i = 0 ; i < m_nNx[1] ; ++ i) { @@ -424,8 +449,7 @@ VolZmap::CreateFromFlatRegion( const ISurfFlatRegion& Surf, double dDimZ, double m_Values[1][nPos1][nCurrentSize].nToolMin = 0 ; m_Values[1][nPos1][nCurrentSize].dMax = dt2 * dLengthX ; m_Values[1][nPos1][nCurrentSize].nToolMax = 0 ; - // Aggiorno il numero della componente connessa - m_Values[1][nPos1][nCurrentSize].nCompo = 0 ; + Point3d ptP1 = ptP0 + dt1 * dLengthX * X_AX ; Point3d ptP2 = ptP0 + dt2 * dLengthX * X_AX ; @@ -434,39 +458,56 @@ VolZmap::CreateFromFlatRegion( const ISurfFlatRegion& Surf, double dDimZ, double for ( int nChunk = 0 ; nChunk < nChunkNum ; ++ nChunk) { + // Indice per ricerca della componente connessa + int nFind = 0 ; + int nLoopNum = Surf.GetLoopCount( nChunk) ; for ( int nLoop = 0 ; nLoop < nLoopNum ; ++ nLoop) { ICurve * pCurve = Surf.GetLoop( nChunk, nLoop) ; - if ( pCurve -> IsPointOn( ptP1)) { + if ( pCurve -> IsPointOn( ptP1, 10 * EPS_SMALL)) { double dP1 ; + Point3d ptTemp1 ; Vector3d vtT1, vtN1 ; pCurve -> GetParamAtPoint(ptP1, dP1) ; - pCurve -> GetPointTang( dP1, ICurve::FROM_MINUS, ptP1, vtT1) ; + pCurve -> GetPointTang( dP1, ICurve::FROM_MINUS, ptTemp1, vtT1) ; vtN1 = vtT1 ^ Z_AX ; // Aggiorno normale inferiore al tratto di dexel m_Values[1][nPos1][nCurrentSize].vtMinN = vtN1 ; + nFind ++ ; } - if ( pCurve -> IsPointOn( ptP2)) { + if ( pCurve -> IsPointOn( ptP2, 10 * EPS_SMALL)) { double dP2 ; + Point3d ptTemp2 ; Vector3d vtT2, vtN2 ; pCurve -> GetParamAtPoint(ptP2, dP2) ; - pCurve -> GetPointTang( dP2, ICurve::FROM_MINUS, ptP1, vtT2) ; + pCurve -> GetPointTang( dP2, ICurve::FROM_MINUS, ptTemp2, vtT2) ; vtN2 = vtT2 ^ Z_AX ; // Aggiorno normale superiore al tratto di dexel m_Values[1][nPos1][nCurrentSize].vtMaxN = vtN2 ; - } + + nFind ++ ; + } + // Se trovata componente esco dal cilco + if ( nFind == 2) { + + // Aggiorno il numero della componente connessa + m_Values[1][nPos1][nCurrentSize].nCompo = nChunk + 1 ; + } } + // Se trovata componente esco dal cilco + if ( nFind == 2) + break ; } } } @@ -477,7 +518,7 @@ VolZmap::CreateFromFlatRegion( const ISurfFlatRegion& Surf, double dDimZ, double m_dMinZ[0] = 0 ; m_dMaxZ[0] = dDimZ ; - if ( bFlag) { + if ( bTriDex) { m_dMinZ[1] = 0 ; m_dMaxZ[1] = dLengthX ; @@ -493,7 +534,7 @@ VolZmap::CreateFromFlatRegion( const ISurfFlatRegion& Surf, double dDimZ, double } // Ridimensiono e setto il vettore dei blocchi a falso - m_nNumBlock = ( bFlag ? m_nFracLin[0] * m_nFracLin[1] * m_nFracLin[2] : m_nFracLin[0] * m_nFracLin[1]) ; + m_nNumBlock = m_nFracLin[0] * m_nFracLin[1] * m_nFracLin[2] ; m_BlockToUpdate.resize( m_nNumBlock) ; for ( unsigned int nCount = 0 ; nCount < m_nNumBlock ; ++ nCount) @@ -509,14 +550,14 @@ VolZmap::CreateFromFlatRegion( const ISurfFlatRegion& Surf, double dDimZ, double //---------------------------------------------------------------------------- bool -VolZmap::CreateFromTriMesh( const ISurfTriMesh& Surf, double dPrec, bool bFlag) +VolZmap::CreateFromTriMesh( const ISurfTriMesh& Surf, double dPrec, bool bTriDex) { // Se la superficie non è chiusa non ha senso continuare if ( ! Surf.IsClosed()) return false ; // Aggiorno la dimensione della mappa 1 o 3 - m_nMapNum = ( bFlag ? 3 : 1) ; + m_nMapNum = ( bTriDex ? 3 : 1) ; // Determino il bounding box della TriMesh BBox3d SurfBBox ; @@ -562,7 +603,7 @@ VolZmap::CreateFromTriMesh( const ISurfTriMesh& Surf, double dPrec, bool bFlag) m_nConnectedCompoCount = - 1 ; // Se Tridexel ridimensiono anche gli altri vettori - if ( bFlag) { + if ( bTriDex) { m_MapFrame[1].Set( ptMapOrig, Y_AX, Z_AX, X_AX) ; // Sarà Front Left m_MapFrame[2].Set( ptMapOrig, Z_AX, X_AX, Y_AX) ; @@ -585,7 +626,22 @@ VolZmap::CreateFromTriMesh( const ISurfTriMesh& Surf, double dPrec, bool bFlag) m_nNy[1] / m_nDexNumPBlock + ( m_nNy[1] % m_nDexNumPBlock >= m_nDexNumPBlock / 2 ? 1 : 0)) ; } + else { + m_MapFrame[1].Set( ptMapOrig, Y_AX, Z_AX, X_AX) ; + m_MapFrame[2].Set( ptMapOrig, Z_AX, X_AX, Y_AX) ; + + m_nNx[1] = 0 ; + m_nNy[1] = 0 ; + m_nDim[1] = 0 ; + + m_nNx[2] = 0 ; + m_nNy[2] = 0 ; + m_nDim[2] = 0 ; + + // Definisco il numero di blocchi lungo z + m_nFracLin[2] = 1 ; + } // Oggetto per calcolo massivo intersezioni IntersParLinesSurfTm intPLSTM( m_MapFrame[0], Surf) ; @@ -667,7 +723,7 @@ VolZmap::CreateFromTriMesh( const ISurfTriMesh& Surf, double dPrec, bool bFlag) } } - if ( bFlag) { + if ( bTriDex) { IntersParLinesSurfTm intPLSTM1( m_MapFrame[1], Surf) ; @@ -830,7 +886,7 @@ VolZmap::CreateFromTriMesh( const ISurfTriMesh& Surf, double dPrec, bool bFlag) m_dMinZ[0] = 0 ; m_dMaxZ[0] = dLengthZ ; - if ( bFlag) { + if ( bTriDex) { m_dMinZ[1] = 0 ; m_dMaxZ[1] = dLengthX ; m_dMinZ[2] = 0 ; @@ -844,7 +900,7 @@ VolZmap::CreateFromTriMesh( const ISurfTriMesh& Surf, double dPrec, bool bFlag) } // Ridimensiono e setto il vettore dei blocchi a falso - m_nNumBlock = ( bFlag ? m_nFracLin[0] * m_nFracLin[1] * m_nFracLin[2] : m_nFracLin[0] * m_nFracLin[1]) ; + m_nNumBlock = m_nFracLin[0] * m_nFracLin[1] * m_nFracLin[2] ; m_BlockToUpdate.resize( m_nNumBlock) ; for ( unsigned int nCount = 0 ; nCount < m_nNumBlock ; ++ nCount) diff --git a/VolTriZmapGraphics.cpp b/VolTriZmapGraphics.cpp index 6621fb1..eae0406 100644 --- a/VolTriZmapGraphics.cpp +++ b/VolTriZmapGraphics.cpp @@ -205,18 +205,14 @@ VolZmap::GetDexelLines( int nDir, int nPos1, int nPos2, POLYLINELIST& lstPL) con bool VolZmap::GetAllTriangles( TRIA3DLIST& lstTria) const { - if ( m_nMapNum == 1) { - const int MAX_DIM_CHUNK = 128 ; - for ( int i = 0 ; i < int( m_nNx[0]) ; i += MAX_DIM_CHUNK) { - int nDimChunkX = min( MAX_DIM_CHUNK, int( m_nNx[0]) - i) ; - for ( int j = 0 ; j < int( m_nNy[0]) ; j += MAX_DIM_CHUNK) { - int nDimChunkY = min( MAX_DIM_CHUNK, int( m_nNy[0]) - j) ; - GetChunkPrisms( i, j, nDimChunkX, nDimChunkY, MAX_DIM_CHUNK, lstTria) ; - } - } + INTVECTOR nModifiedBlocks ; + TRIA3DLISTVECTOR vLstTria ; + if ( ! GetTriangles( true, nModifiedBlocks, vLstTria)) + return false ; + lstTria.clear() ; + for ( size_t i = 0 ; i < vLstTria.size() ; ++ i) { + lstTria.splice( lstTria.end(), vLstTria[i]) ; } - else - MarchingCubes( lstTria) ; return true ; } @@ -225,6 +221,17 @@ VolZmap::GetAllTriangles( TRIA3DLIST& lstTria) const bool VolZmap::GetTriangles( bool bAllBlocks, INTVECTOR& nModifiedBlocks, TRIA3DLISTVECTOR& vLstTria) const { + // Se nessun blocco modificato, è richiesta esterna e li considero tutti modificati + bool bSomeModif = false ; + for ( size_t i = 0 ; i < m_nNumBlock ; ++ i) { + if ( m_BlockToUpdate[i]) { + bSomeModif = true ; + break ; + } + } + if ( ! bSomeModif) + bAllBlocks = true ; + // Caso di singola mappa if ( m_nMapNum == 1) { @@ -281,82 +288,82 @@ VolZmap::GetTriangles( bool bAllBlocks, INTVECTOR& nModifiedBlocks, TRIA3DLISTVE TriaMatrix VecTriHold ; VecTriHold.resize( m_nNumBlock) ; + + bool bCalcInterBlock = false ; - // Ciclo sui blocchi + // Calcolo i triangoli 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 - 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 + // Se il blocco deve essere processato if ( bAllBlocks || m_BlockToUpdate[t]) { + // processo ... vLstTria.emplace_back() ; nModifiedBlocks[t] = int( vLstTria.size()) - 1 ; m_InterBlockTria[t].clear() ; - ExtMarchingCubes( LimitIndexes, t, vLstTria.back(), VecTriHold[t]) ; - - // Flipping fra voxel interni - FlipEdgesII( VecTriHold[t]) ; - + #if 1 + ExtMarchingCubes( int( t), vLstTria.back(), VecTriHold[t]) ; + // Flipping fra voxel interni + FlipEdgesII( VecTriHold[t]) ; + bCalcInterBlock = true ; + #else + MarchingCubes( int( t), vLstTria.back()) ; + #endif m_BlockToUpdate[t] = false ; } else nModifiedBlocks[t] = -1 ; } - // Copio i triangoli di frontiera in una matrice gemella - // di m_InterBlockTria per avere sempre a disposizione - // i triangoli non flippati. + // Calcolo i triangoli di frontiera tra feature di blocchi diversi + // 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) ; - - + if ( bCalcInterBlock) { + InterBlockTria = m_InterBlockTria ; + FlipEdgesBB( InterBlockTria) ; + } - // ciclo sui blocchi + // Inserisco in lista i triangoli di feature derivanti dai blocchi for ( size_t t = 0 ; t < m_nNumBlock ; ++ t) { - // ciclo sui voxel del blocco - for ( size_t t1 = 0 ; t1 < VecTriHold[t].size() ; ++ t1) { - // ciclo sulle componenti connesse del voxel - for ( size_t t2 = 0 ; t2 < VecTriHold[t][t1].vCompoTria.size() ; ++ t2) { - // ciclo sui triangoli delle componenti connesse - for ( size_t t3 = 0 ; t3 < VecTriHold[t][t1].vCompoTria[t2].size() ; ++ t3) { - // aggiungo un singolo triangolo alla lista - if ( nModifiedBlocks[t] >= 0) + if ( nModifiedBlocks[t] >= 0) { + // ciclo sui voxel del blocco + for ( size_t t1 = 0 ; t1 < VecTriHold[t].size() ; ++ t1) { + // ciclo sulle componenti connesse del voxel + for ( size_t t2 = 0 ; t2 < VecTriHold[t][t1].vCompoTria.size() ; ++ t2) { + // ciclo sui triangoli delle componenti connesse + for ( size_t t3 = 0 ; t3 < VecTriHold[t][t1].vCompoTria[t2].size() ; ++ t3) { + // aggiungo triangolo alla lista vLstTria[nModifiedBlocks[t]].emplace_back( VecTriHold[t][t1].vCompoTria[t2][t3]) ; - } - } - } - } - - - 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]) ; } } } } + } + + // Inserisco in lista i triangoli di frontiera tra feature di blocchi diversi + if ( bCalcInterBlock) { + 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.back() = int( nPos) ; } - nModifiedBlocks[nModifiedBlocks.size()-1] = int( vLstTria.size() - 1) ; + else + nModifiedBlocks.back() = - 1 ; } return true ; @@ -560,106 +567,13 @@ VolZmap::AddDexelSideFace( int nPos, int nPosAdj, const Point3d& ptP, const Poin return true ; } -//---------------------------------------------------------------------------- -bool -VolZmap::MarchingCubes( TRIA3DLIST& lstTria) const -{ - // Ciclo su tutti i voxel dello Zmap - for ( int i = - 1 ; i < int( m_nNx[0]) ; ++ i) { - for ( int j = - 1 ; j < int( m_nNy[0]) ; ++ j) { - for ( int k = - 1 ; k < int( m_nNy[1]) ; ++ k) { - - // Indici i,j,k dei vertici - int IndexCorner[8][3] = { - { i, j, k}, - { i + 1, j, k}, - { i + 1, j + 1, k}, - { i, j + 1, k}, - { i, j, k + 1}, - { i + 1, j, k + 1}, - { i + 1, j + 1, k + 1}, - { i, j + 1, k + 1} - } ; - - // Classificazione dei vertici: interni o esterni al materiale - int nIndex = 0 ; - if ( IsThereMat( i, j, k)) - nIndex |= ( 1 << 0) ; - if ( IsThereMat( i + 1, j, k)) - nIndex |= ( 1 << 1) ; - if ( IsThereMat( i + 1, j + 1, k)) - nIndex |= ( 1 << 2) ; - if ( IsThereMat( i, j + 1, k)) - nIndex |= ( 1 << 3) ; - if ( IsThereMat( i, j, k + 1)) - nIndex |= ( 1 << 4) ; - if ( IsThereMat( i + 1, j, k + 1)) - nIndex |= ( 1 << 5) ; - if ( IsThereMat( i + 1, j + 1, k + 1)) - nIndex |= ( 1 << 6) ; - if ( IsThereMat( i, j + 1, k + 1)) - nIndex |= ( 1 << 7) ; - - // Se vi è qualche intersezione fra segmenti e superficie - // continuo altrimenti passo al prossimo voxel - if ( EdgeTable[nIndex] == 0) - continue ; - - static int intersections[12][2] = { - { 0, 1 }, { 1, 2 }, { 2, 3 }, { 3, 0 }, { 4, 5 }, { 5, 6 }, - { 6, 7 }, { 7, 4 }, { 0, 4 }, { 1, 5 }, { 2, 6 }, { 3, 7 } - } ; - - Point3d ptIntPoint[12] ; - - // Ciclo sui segmenti - for ( int EdgeIndex = 0 ; EdgeIndex < 12 ; ++ EdgeIndex) { - // Se il segmento non attraversa la superficie - // passo al successivo - if ( ! ( EdgeTable[nIndex] & ( 1 << EdgeIndex))) - continue ; - - int n1 = intersections[EdgeIndex][0] ; - int n2 = intersections[EdgeIndex][1] ; - - // Determino con precisione il punto di intersezione sullo spigolo - IntersPos( IndexCorner[n1], IndexCorner[n2], ptIntPoint[EdgeIndex]) ; - - ptIntPoint[EdgeIndex].ToGlob( m_MapFrame[0]) ; - } - - // Costruzione dei triangoli - for ( int TriIndex = 0 ; TriangleTableEn[nIndex][0][TriIndex] != - 1 ; TriIndex += 3) { - - // Costruzione triangolo - int i0 = TriangleTableEn[nIndex][0][TriIndex + 2] ; - int i1 = TriangleTableEn[nIndex][0][TriIndex + 1] ; - int i2 = TriangleTableEn[nIndex][0][TriIndex] ; - - // Il triangolo è pronto - Triangle3d CurrentTriangle ; - CurrentTriangle.Set( ptIntPoint[i0], ptIntPoint[i1], ptIntPoint[i2]) ; - CurrentTriangle.Validate() ; - - // Aggiungo triangolo - lstTria.emplace_back( CurrentTriangle) ; - } - } - } - } - - return true ; -} - //---------------------------------------------------------------------------- bool VolZmap::MarchingCubes( int nBlock, TRIA3DLIST& lstTria) const { - if ( nBlock < 0 || nBlock >= int( m_BlockToUpdate.size())) + if ( nBlock < 0 || nBlock >= int( m_nNumBlock)) return false ; - Point3d ptMapOrig = m_MapFrame[0].Orig() ; - // Calcolo posizione del blocco nel reticolo int nIBlock = ( nBlock % int( m_nFracLin[0] * m_nFracLin[1])) % int( m_nFracLin[0]) ; int nJBlock = ( nBlock % int( m_nFracLin[0] * m_nFracLin[1])) / int( m_nFracLin[0]) ; @@ -667,18 +581,18 @@ VolZmap::MarchingCubes( int nBlock, TRIA3DLIST& lstTria) const // Calcolo limiti per l'indice i int nStartI = nIBlock * int( m_nDexNumPBlock) - 1 ; - int nEndI = ( nIBlock + 1 == int( m_nFracLin[0]) ? - int( m_nNx[0]) : ( nIBlock + 1) * int( m_nDexNumPBlock)) ; + int nEndI = ( nIBlock + 1 == int( m_nFracLin[0]) ? + int( m_nNx[0]) : ( nIBlock + 1) * int( m_nDexNumPBlock)) ; // Calcolo limiti per l'indice j int nStartJ = nJBlock * int( m_nDexNumPBlock) - 1 ; - int nEndJ = ( nJBlock + 1 == int( m_nFracLin[1]) ? - int( m_nNy[0]) : ( nJBlock + 1) * int( m_nDexNumPBlock)) ; + int nEndJ = ( nJBlock + 1 == int( m_nFracLin[1]) ? + int( m_nNy[0]) : ( nJBlock + 1) * int( m_nDexNumPBlock)) ; // Calcolo limiti per l'indice k int nStartK = nKBlock * int( m_nDexNumPBlock) - 1 ; - int nEndK = ( nKBlock + 1 == int( m_nFracLin[2]) ? - int( m_nNy[1]) : ( nKBlock + 1) * int( m_nDexNumPBlock)) ; + int nEndK = ( nKBlock + 1 == int( m_nFracLin[2]) ? + int( m_nNy[1]) : ( nKBlock + 1) * int( m_nDexNumPBlock)) ; // Ciclo su tutti i voxel dello Zmap for ( int i = nStartI ; i < nEndI ; ++ i) { @@ -773,24 +687,27 @@ VolZmap::MarchingCubes( int nBlock, TRIA3DLIST& lstTria) const //---------------------------------------------------------------------------- bool -VolZmap::ExtMarchingCubes( const int nLimits[], size_t nBlockNumber, TRIA3DLIST& lstTria, TriHolder& triHold) const +VolZmap::ExtMarchingCubes( int nBlock, TRIA3DLIST& lstTria, TriHolder& triHold) const { - Point3d ptMapOrig = m_MapFrame[0].Orig() ; + if ( nBlock < 0 || nBlock >= int( m_nNumBlock)) + return false ; + + // Calcolo i limiti sugli indici dei voxel del blocco + int nIJK[3] ; + // Vettore indici i,j,k del blocco + GetBlockIJKFromN( nBlock, nIJK) ; + // Vettore limiti sugli indici dei voxel del blocco + int nLimits[6] ; + GetBlockLimitsIJK( nIJK, nLimits) ; // 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) { - /*if ( !( i == 2 && j == 14 && k == 24)) - continue ;*/ // Riconoscimento dei voxel di frontiera - bool bBoundary = false ; int nVoxIndexes[3] = { i, j, k} ; - - if ( IsAVoxelOnBoundary( nLimits, nVoxIndexes, true)) - bBoundary = true ; - + bool bBoundary = IsAVoxelOnBoundary( nLimits, nVoxIndexes, true) ; // Indici i,j,k dei vertici int IndexCorner[8][3] = { @@ -804,9 +721,8 @@ VolZmap::ExtMarchingCubes( const int nLimits[], size_t nBlockNumber, TRIA3DLIST& { i, j + 1, k + 1} } ; - int nIndex = 0 ; - // Classificazione dei vertici: interni o esterni al materiale + int nIndex = 0 ; if ( IsThereMat( i, j, k)) nIndex |= ( 1 << 0) ; if ( IsThereMat( i + 1, j, k)) @@ -1289,7 +1205,6 @@ VolZmap::ExtMarchingCubes( const int nLimits[], size_t nBlockNumber, TRIA3DLIST& Point3d ptSol( dUnknownVector( 0) + vtO.x, dUnknownVector( 1) + vtO.y, dUnknownVector( 2) + vtO.z) ; - Triangle3d CurrentTriangle ; TRIA3DVECTOR triContainer ; @@ -1312,10 +1227,9 @@ VolZmap::ExtMarchingCubes( const int nLimits[], size_t nBlockNumber, TRIA3DLIST& triContainer.emplace_back( CurrentTriangle) ; // Controllo sull'inversione delle normali - if ( ( CurrentTriangle.GetN() * CompoVert[nj].vtNorm < - 0.01 && - CurrentTriangle.GetN() * CompoVert[ni].vtNorm < - 0.01) && - ( ! bInvNormB)) { - + if ( ( CurrentTriangle.GetN() * CompoVert[nj].vtNorm < - 0.01 && // - 0.01 passando a - 0.002 si tolgono i 4 restanti + CurrentTriangle.GetN() * CompoVert[ni].vtNorm < - 0.01) && // - 0.01 passando a - 0.002 si tolgono i 4 restanti + ! bInvNormB) { ptSol = ptGravityCenter ; bInvNormB = true ; triContainer.resize( 0) ; @@ -1328,7 +1242,7 @@ VolZmap::ExtMarchingCubes( const int nLimits[], size_t nBlockNumber, TRIA3DLIST& // interno o esterno al voxel a cui appartiene. bool bInsideVoxel = IsPointInsideVoxelApprox( i, j, k, ptSol) ; - // Proprietà gemoetriche locali + // Proprietà geometriche locali // del campo vettoriale: se vero // procediamo con un ulteriore // analisi per eliminare triangoli @@ -1361,23 +1275,22 @@ VolZmap::ExtMarchingCubes( const int nLimits[], size_t nBlockNumber, TRIA3DLIST& // esito negativo. if ( nFeatureType == 2 && bLocProp && ! ( bOutside || bInsideVoxel || bInvNormB)) { - + if ( abs( nMin1 - nMin2) == 1 || - abs( nMin1 - nMin2) == 3) { + abs( nMin1 - nMin2) == 3) { int nSum = nMin1 + nMin2 ; int nTriNum = ( nSum == 3 && ( nMin1 == 3 || nMin2 == 3) ? - max( nMin1, nMin2) : min( nMin1, nMin2)) ; + max( nMin1, nMin2) : min( nMin1, nMin2)) ; double dDot1 = triContainer[nTriNum].GetN() * CompoVert[nMin1].vtNorm ; double dDot2 = triContainer[nTriNum].GetN() * CompoVert[nMin2].vtNorm ; if ( ( dDot1 < - 0.2 && dDot2 > - EPS_ZERO) || - ( dDot2 < - 0.2 && dDot1 > - EPS_ZERO)) { + ( dDot2 < - 0.2 && dDot1 > - EPS_ZERO)) { int nNm = dDot1 < - 0.2 ? nMin1 : nMin2 ; int nNp = dDot1 < - 0.2 ? nMin2 : nMin1 ; - Vector3d vtVV = CompoVert[nNp].ptInt - CompoVert[nNm].ptInt ; Point3d ptSolZmFrame = ptSol ; @@ -1409,14 +1322,14 @@ VolZmap::ExtMarchingCubes( const int nLimits[], size_t nBlockNumber, TRIA3DLIST& for ( int nc = 0 ; nc < nVertComp ; ++ nc) { int nd = ( nc + 1 < nVertComp) ? nc + 1 : 0 ; - // Il triangolo è pronto + // Il triangolo è pronto CurrentTriangle.Set( ptSol, CompoVert[nd].ptInt, CompoVert[nc].ptInt) ; CurrentTriangle.Validate( true) ; - // Aggiungo triangolo al vettore temporaneo + // Aggiungo triangolo al vettore temporaneo triContainer.emplace_back( CurrentTriangle) ; } - } + } } } else { @@ -1431,17 +1344,17 @@ VolZmap::ExtMarchingCubes( const int nLimits[], size_t nBlockNumber, TRIA3DLIST& if ( CompoVert[nPrev1].vtNorm * CompoVert[nMin1].vtNorm > 0) { - bool bTestOnVert = abs( nPrev1 - nMin2) == 0 || abs( nPrev1 - nMin2) == 3 ; + bool bTestOnVert = abs( nPrev1 - nMin2) == 1 || abs( nPrev1 - nMin2) == 3 ; nNeighbourIndex = bTestOnVert ? nPrev1 : nNext1 ; nStartIndex = nMin2 ; } else { - bool bTestOnVert = abs( nPrev2 - nMin1) == 0 || abs( nPrev2 - nMin1) == 3 ; + bool bTestOnVert = abs( nPrev2 - nMin1) == 1 || abs( nPrev2 - nMin1) == 3 ; nNeighbourIndex = bTestOnVert ? nPrev2 : nNext2 ; nStartIndex = nMin1 ; } - + Vector3d vtVV = CompoVert[nNeighbourIndex].ptInt - CompoVert[nStartIndex].ptInt ; double dVVLen = vtVV.Len() ; vtVV.Normalize() ; @@ -1449,9 +1362,9 @@ VolZmap::ExtMarchingCubes( const int nLimits[], size_t nBlockNumber, TRIA3DLIST& Vector3d vtVF = ptSol - CompoVert[nStartIndex].ptInt ; Vector3d vtNewVF = ( vtVF * vtVV) * vtVV ; - double dNewVFLen = vtNewVF.Len() ; + double dNewVFLen = vtNewVF.Len() ; - if ( dNewVFLen > dVVLen) { + if ( dNewVFLen > dVVLen || vtVV * vtNewVF < 0) { ptSol = CompoVert[nNeighbourIndex].ptInt ; } @@ -1459,23 +1372,21 @@ VolZmap::ExtMarchingCubes( const int nLimits[], size_t nBlockNumber, TRIA3DLIST& ptSol = CompoVert[nStartIndex].ptInt + vtNewVF ; } - triContainer.resize( 0) ; for ( int nc = 0 ; nc < nVertComp ; ++ nc) { int nd = ( nc + 1 < nVertComp) ? nc + 1 : 0 ; - // Il triangolo è pronto + // Il triangolo è pronto CurrentTriangle.Set( ptSol, CompoVert[nd].ptInt, CompoVert[nc].ptInt) ; CurrentTriangle.Validate( true) ; - // Aggiungo triangolo al vettore temporaneo + // Aggiungo triangolo al vettore temporaneo triContainer.emplace_back( CurrentTriangle) ; - } - } + } + } } - // Valuto normali: questo è ancora un controllo // sulle normali, se risultano in tutti i punti @@ -1512,7 +1423,7 @@ VolZmap::ExtMarchingCubes( const int nLimits[], size_t nBlockNumber, TRIA3DLIST& // Se l'area è non nulla aggiungo il triangolo // a uno dei due vettori. - if ( triContainer[ni].GetArea() > EPS_SMALL) { + if ( triContainer[ni].GetArea() > EPS_SMALL * EPS_SMALL) { int nVoxIJK[3] = { i, j, k} ; @@ -1558,7 +1469,6 @@ VolZmap::ExtMarchingCubes( const int nLimits[], size_t nBlockNumber, TRIA3DLIST& // nel vector. int nCurrent = int( triHold.size()) - 1 ; - // Aggiungo vertice della componente // connessa alla lista dei vertici. triHold[nCurrent].ptCompoVert.emplace_back( ptSol) ; @@ -1586,36 +1496,36 @@ VolZmap::ExtMarchingCubes( const int nLimits[], size_t nBlockNumber, TRIA3DLIST& // la struttura dati del voxel. if ( nBorderFeatureInVoxel == 1) { - m_InterBlockTria[nBlockNumber].resize( m_InterBlockTria[nBlockNumber].size() + 1) ; + m_InterBlockTria[nBlock].resize( m_InterBlockTria[nBlock].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 ; + int nCurrent = int( m_InterBlockTria[nBlock].size()) - 1 ; - m_InterBlockTria[nBlockNumber][nCurrent].i = i ; - m_InterBlockTria[nBlockNumber][nCurrent].j = j ; - m_InterBlockTria[nBlockNumber][nCurrent].k = k ; + m_InterBlockTria[nBlock][nCurrent].i = i ; + m_InterBlockTria[nBlock][nCurrent].j = j ; + m_InterBlockTria[nBlock][nCurrent].k = k ; } // Indice che identifica la posizione del voxel // nel vector. - int nCurrent = int( m_InterBlockTria[nBlockNumber].size()) - 1 ; + int nCurrent = int( m_InterBlockTria[nBlock].size()) - 1 ; // Aggiungo vertice della componente // connessa alla lista dei vertici. - m_InterBlockTria[nBlockNumber][nCurrent].ptCompoVert.emplace_back( ptSol) ; + m_InterBlockTria[nBlock][nCurrent].ptCompoVert.emplace_back( ptSol) ; - int nOldFeatureNum = int( m_InterBlockTria[nBlockNumber][nCurrent].vCompoTria.size()) ; + int nOldFeatureNum = int( m_InterBlockTria[nBlock][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) ; + m_InterBlockTria[nBlock][nCurrent].vCompoTria.resize( nNewFeatureNum) ; for ( int ni = 0 ; ni < int( vBorderTriaTemp.size()) ; ++ ni) - m_InterBlockTria[nBlockNumber][nCurrent].vCompoTria[nOldFeatureNum].emplace_back( vBorderTriaTemp[ni]) ; + m_InterBlockTria[nBlock][nCurrent].vCompoTria[nOldFeatureNum].emplace_back( vBorderTriaTemp[ni]) ; } } // ExtMC non confermato, si passa a MC @@ -1634,9 +1544,8 @@ VolZmap::ExtMarchingCubes( const int nLimits[], size_t nBlockNumber, TRIA3DLIST& CurrentTriangle.Set( VecField[i0].ptInt, VecField[i1].ptInt, VecField[i2].ptInt) ; CurrentTriangle.Validate( true) ; - // Se il triangolo non è degenere lo aggiungo alla lista - if ( CurrentTriangle.GetArea() > EPS_SMALL) - lstTria.emplace_back( CurrentTriangle) ; + // Aggiungo alla lista + lstTria.emplace_back( CurrentTriangle) ; } } } @@ -1657,9 +1566,8 @@ VolZmap::ExtMarchingCubes( const int nLimits[], size_t nBlockNumber, TRIA3DLIST& CurrentTriangle.Set( VecField[i0].ptInt, VecField[i1].ptInt, VecField[i2].ptInt) ; CurrentTriangle.Validate( true) ; - // Se il triangolo non è degenere lo aggiungo alla lista - if ( CurrentTriangle.GetArea() > EPS_SMALL) - lstTria.emplace_back( CurrentTriangle) ; + // Aggiungo alla lista + lstTria.emplace_back( CurrentTriangle) ; } } nTableOffset += nVertComp ; @@ -1675,28 +1583,18 @@ VolZmap::ExtMarchingCubes( const int nLimits[], size_t nBlockNumber, TRIA3DLIST& 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) { + for ( int n1 = 0 ; n1 < nVoxelNum ; ++ n1) { - // Determino estremi del ciclo sui voxel interno - int nSt2 = nSt1 + 1 ; - int nEn2 = nVoxelNum ; - - for ( int n2 = nSt2 ; n2 < nEn2 ; ++ n2) { + for ( int n2 = n1 ; n2 < nVoxelNum ; ++ 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) { + if ( ( TriHold[n2].i >= TriHold[n1].i && TriHold[n2].i <= TriHold[n1].i + 1) || + ( TriHold[n2].j >= TriHold[n1].j && TriHold[n2].j <= TriHold[n1].j + 1) || + ( TriHold[n2].k >= TriHold[n1].k && TriHold[n2].k <= TriHold[n1].k + 1)) { // Numero delle componenti connesse nei due voxel int nNumCompo1 = int( TriHold[n1].ptCompoVert.size()) ; @@ -1707,7 +1605,7 @@ VolZmap::FlipEdgesII( TriHolder& TriHold) const // Ciclo sulle componenti for ( ; nCompo1 < nNumCompo1 ; ++ nCompo1) { - int nCompo2 = 0 ; + int nCompo2 = ( n1 == n2 ? nCompo1 + 1 : 0) ; for ( ; nCompo2 < nNumCompo2 ; ++ nCompo2) { @@ -1730,13 +1628,13 @@ VolZmap::FlipEdgesII( TriHolder& TriHold) const Point3d ptP1 = TriHold[n1].vCompoTria[nCompo1][nTri1].GetP( nVert1) ; Point3d ptP2 = TriHold[n2].vCompoTria[nCompo2][nTri2].GetP( nVert2) ; - if ( SqDist( ptP1, ptP2) < dSqEps) { + if ( AreSamePointEpsilon( ptP1, ptP2, EPS_SMALL)) { Point3d ptVert1 = TriHold[n1].ptCompoVert[nCompo1] ; Point3d ptVert2 = TriHold[n2].ptCompoVert[nCompo2] ; - if ( ! ( AreSamePointApprox( ptP1, ptVert1) || - AreSamePointApprox( ptP2, ptVert2))) { + if ( ! ( AreSamePointEpsilon( ptP1, ptVert1, EPS_SMALL) || + AreSamePointEpsilon( ptP2, ptVert2, EPS_SMALL))) { SharedIndex.emplace_back( nVert1) ; SharedIndex.emplace_back( nVert2) ; @@ -1753,20 +1651,10 @@ VolZmap::FlipEdgesII( TriHolder& TriHold) const // Si deve operare la modifica dei triangoli if ( SharedIndex.size() > 2) { - + + // verifico che i due lai adiacenti siano controversi 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) ; TriHold[n1].vCompoTria[nCompo1][nTri1].SetP( SharedIndex[0], TriHold[n2].ptCompoVert[nCompo2]) ; @@ -1796,28 +1684,26 @@ VolZmap::FlipEdgesII( TriHolder& TriHold) const //---------------------------------------------------------------------------- bool -VolZmap::FlipEdgesBB( TriaMatrix& InterTria) const { - - double dSqEps = EPS_SMALL * EPS_SMALL ; - +VolZmap::FlipEdgesBB( TriaMatrix& InterTria) const +{ // Numero di blocchi size_t nBlocksNum = InterTria.size() ; // ciclo sui blocchi for ( size_t tFB = 0 ; tFB < nBlocksNum - 1 ; ++ tFB) { + + int nFBijk[3] ; + GetBlockIJKFromN( int( tFB), nFBijk) ; + 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 @@ -1861,13 +1747,13 @@ VolZmap::FlipEdgesBB( TriaMatrix& InterTria) const { 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) { + if ( AreSamePointEpsilon( ptPF, ptPL, EPS_SMALL)) { Point3d ptVertF = InterTria[tFB][tVFB].ptCompoVert[tCmpF] ; Point3d ptVertL = InterTria[tLB][tVLB].ptCompoVert[tCmpL] ; - if ( ! ( AreSamePointApprox( ptPF, ptVertF) || - AreSamePointApprox( ptPL, ptVertL))) { + if ( ! ( AreSamePointEpsilon( ptPF, ptVertF, EPS_SMALL) || + AreSamePointEpsilon( ptPL, ptVertL, EPS_SMALL))) { SharedIndex.emplace_back( nVertF) ; SharedIndex.emplace_back( nVertL) ; @@ -1885,21 +1771,10 @@ VolZmap::FlipEdgesBB( TriaMatrix& InterTria) const { // Si deve operare la modifica dei triangoli if ( SharedIndex.size() > 2) { + // verifico che i due lai adiacenti siano controversi 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], @@ -1910,7 +1785,6 @@ VolZmap::FlipEdgesBB( TriaMatrix& InterTria) const { bModified = true ; break ; - } } } diff --git a/VolTriZmapVolume.cpp b/VolTriZmapVolume.cpp index 6ee9574..6f9dd7a 100644 --- a/VolTriZmapVolume.cpp +++ b/VolTriZmapVolume.cpp @@ -185,15 +185,15 @@ VolZmap::SubtractIntervals( unsigned int nGrid, unsigned int nI, unsigned int nJ int nMinK = max( int( floor( ( ( dMin - 0.5 * m_dStep) / m_dStep - EPS_SMALL))), 0) ; int nMaxK = max( int( floor( ( ( dMax + 0.5 * m_dStep) / m_dStep + EPS_SMALL))), 0) ; - int nMinZBlock = max( 0, nMinK / int( m_nDexNumPBlock)) ; + int nMinZBlock = ( m_nMapNum == 1 ? 0 : max( 0, nMinK / int( m_nDexNumPBlock))) ; int nMaxZBlock = min( int( m_nFracLin[2] - 1), nMaxK / int( m_nDexNumPBlock)) ; - + for ( int k = nMinZBlock ; k <= nMaxZBlock ; ++ k) { - int nBlockNum = k * nLayerBlock + nYBlock * m_nFracLin[0] + nXBlock ; m_BlockToUpdate[nBlockNum] = true ; } } + else if ( nGrid == 1) { int nYBlock = min( nI / m_nDexNumPBlock, m_nFracLin[1] - 1) ; @@ -210,6 +210,7 @@ VolZmap::SubtractIntervals( unsigned int nGrid, unsigned int nI, unsigned int nJ m_BlockToUpdate[nBlockNum] = true ; } } + else if ( nGrid == 2) { int nXBlock = min( nJ / m_nDexNumPBlock, m_nFracLin[0] - 1) ; @@ -222,7 +223,6 @@ 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) { - int nBlockNum = nZBlock * nLayerBlock + k * m_nFracLin[0] + nXBlock ; m_BlockToUpdate[nBlockNum] = true ; } diff --git a/VolZmap.cpp b/VolZmap.cpp index 0cc6aaf..27164f0 100644 --- a/VolZmap.cpp +++ b/VolZmap.cpp @@ -578,12 +578,3 @@ VolZmap::IsMapConnected( void) } return true ; } - - - - - - - - - diff --git a/VolZmap.h b/VolZmap.h index a458ca4..a823539 100644 --- a/VolZmap.h +++ b/VolZmap.h @@ -122,9 +122,8 @@ class VolZmap : public IVolZmap, public IGeoObjRW bool CalcDexelPrisms( int nPos1, int nPos2, TRIA3DLIST& lstTria) const ; bool AddDexelSideFace( int nPos, int nPosAdj, const Point3d& ptP, const Point3d& ptQ, 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[], size_t nBlockNumber, TRIA3DLIST& lstTria, TriHolder& triHold) const ; + bool ExtMarchingCubes( int nBlock, TRIA3DLIST& lstTria, TriHolder& triHold) const ; bool FlipEdgesII( TriHolder& TriHold) const ; bool FlipEdgesBB( TriaMatrix& InterTria) const ; bool IsThereMat( int nI, int nJ, int nK) const ;