//---------------------------------------------------------------------------- // EgalTech 2015-2016 //---------------------------------------------------------------------------- // File : VolZmap.cpp Data : 22.01.15 Versione : 1.6a4 // Contenuto : Implementazione della classe Volume Zmap (tre griglie) // // // // Modifiche : 22.01.15 DS Creazione modulo. // // //---------------------------------------------------------------------------- //--------------------------- Include ---------------------------------------- #include "stdafx.h" #include "CurveLine.h" #include "VolZmap.h" #include "GeoConst.h" #include "IntersLineSurfTm.h" #include "MC_Tables.h" #include "\EgtDev\Include\EGkIntervals.h" #include "\EgtDev\Include\EgtNumUtils.h" #include "\EgtDev\Include\EGkStringUtils3d.h" #include "\EgtDev\Extern\Eigen\Core" #include "\EgtDev\Extern\Eigen\SVD" using namespace std ; // ------------------------- STRUTTURA VERTICE TRIANGOLO - NORMALE ALLA SUPERFICIE ------------------------------------------------ struct VectorField { Point3d ptInt ; Vector3d vtNorm ; } ; // ------------------------- TABELLA BLOCCHI ADIACENTI ---------------------------------------------------------------------------- static int NeighbourTable[8][4] = { {0, -1, -1, -1}, {1, 1, -1, -1}, {1, 1, 2, -1}, {2, 1, 2, -1}, {1, 3, -1, -1}, {2, 1, 3, -1}, {2, 2, 3, -1}, {3, 1, 2, 3} } ; // ------------------------- FUNZIONE TEST SULLE NORMALI -------------------------------------------------------------------------- enum FatureType { NoFeature = 0, Corner = 1, Edge = 2} ; //---------------------------------------------------------------------------- bool TestOnNormal( const VectorField CompoVert[], int nCompoElem, int& FeatureType) { int nI, nJ ; double dMinCosTheta = 1.001 ; const double dCosThetaSharp = 0.9 ; for ( int i = 0 ; i < nCompoElem ; ++ i) { for ( int j = i + 1 ; j < nCompoElem ; ++ j) { double dCurrentCos = CompoVert[i].vtNorm * CompoVert[j].vtNorm ; if ( dCurrentCos < dMinCosTheta) { nI = i ; nJ = j ; dMinCosTheta = dCurrentCos ; } } } if ( dMinCosTheta >= dCosThetaSharp) { FeatureType = NoFeature ; return false ; } Vector3d vtK = CompoVert[nI].vtNorm ^ CompoVert[nJ].vtNorm ; const double dCosPhiCorner = 0.5 ; for ( int i = 0 ; i < nCompoElem ; ++ i) { double dAbsCurrentCos = abs( CompoVert[i].vtNorm * vtK) ; if ( dAbsCurrentCos > dCosPhiCorner) { FeatureType = Corner ; return true ; } } FeatureType = Edge ; return true ; } //---------------------------------------------------------------------------- bool TestOnSlice( const VectorField CompoVert[], int nLimArrey, int i, int j, int k, double dStep) { // Coordinate dei piani double dX0 = ( i + 0.5) * dStep ; double dX1 = ( i + 1.5) * dStep ; double dY0 = ( j + 0.5) * dStep ; double dY1 = ( j + 1.5) * dStep ; double dZ0 = ( k + 0.5) * dStep ; double dZ1 = ( k + 1.5) * dStep ; // Analisi faccia 1 int nCount = 0 ; double dDelta = abs( CompoVert[0].ptInt.x - dX0) ; while ( dDelta < EPS_SMALL && nCount < nLimArrey) { ++ nCount ; dDelta = abs( CompoVert[nCount].ptInt.x - dX0) ; } if ( nCount == nLimArrey) return false ; // Analisi faccia 2 nCount = 0 ; dDelta = abs( CompoVert[0].ptInt.x - dX1) ; while ( dDelta < EPS_SMALL && nCount < nLimArrey) { ++ nCount ; dDelta = abs( CompoVert[nCount].ptInt.x - dX1) ; } if ( nCount == nLimArrey) return false ; // Analisi faccia 3 nCount = 0 ; dDelta = abs( CompoVert[0].ptInt.y - dY0) ; while ( dDelta < EPS_SMALL && nCount < nLimArrey) { ++ nCount ; dDelta = abs( CompoVert[nCount].ptInt.y - dY0) ; } if ( nCount == nLimArrey) return false ; // Analisi faccia 4 nCount = 0 ; dDelta = abs( CompoVert[0].ptInt.y - dY1) ; while ( dDelta < EPS_SMALL && nCount < nLimArrey) { ++ nCount ; dDelta = abs( CompoVert[nCount].ptInt.y - dY1) ; } if ( nCount == nLimArrey) return false ; // Analisi faccia 5 nCount = 0 ; dDelta = abs( CompoVert[0].ptInt.z - dZ0) ; while ( dDelta < EPS_SMALL && nCount < nLimArrey) { ++ nCount ; dDelta = abs( CompoVert[nCount].ptInt.z - dZ0) ; } if ( nCount == nLimArrey) return false ; // Analisi faccia 6 nCount = 0 ; dDelta = abs( CompoVert[0].ptInt.z - dZ1) ; while ( dDelta < EPS_SMALL && nCount < nLimArrey) { ++ nCount ; dDelta = abs( CompoVert[nCount].ptInt.z - dZ1) ; } if ( nCount == nLimArrey) return false ; return true ; } // ------------------------- VISUALIZZAZIONE -------------------------------------------------------------------------------------- //---------------------------------------------------------------------------- bool VolZmap::GetDexelLines( int nDir, int nPos1, int nPos2, POLYLINELIST& lstPL) const { // Se richiesti spilloni ( 0 <= nDir < 3) if ( nDir < 3) { // Controllo l'ammissibilità della griglia if ( nDir < 0 || nDir > 2) return false ; // Verifiche sugli indici if ( nPos1 < 0 || nPos1 >= int( m_nNx[nDir]) || nPos2 < 0 || nPos2 >= int( m_nNy[nDir])) return false ; int nPos = nPos1 + nPos2 * m_nNx[nDir] ; if ( nPos < 0 || nPos >= int( m_Values[nDir].size())) return false ; // Calcolo coordinate punto double dX = m_dStep * ( 0.5 + nPos1) ; double dY = m_dStep * ( 0.5 + nPos2) ; // Determino il punto di partenza sulla griglia Point3d ptP = m_MapFrame[nDir].Orig() + dX * m_MapFrame[nDir].VersX() + dY * m_MapFrame[nDir].VersY() ; // Creo le polilinee for ( int j = 1 ; j < int( m_Values[nDir][nPos].size()) ; j += 2) { // aggiungo polilinea a lista lstPL.emplace_back() ; // inserisco punti estremi lstPL.back().AddUPoint( 0, ptP + m_Values[nDir][nPos][j-1].dZVal * m_MapFrame[nDir].VersZ()) ; lstPL.back().AddUPoint( 1, ptP + m_Values[nDir][nPos][j].dZVal * m_MapFrame[nDir].VersZ()) ; } return true ; } // altrimenti richieste normali ( 3 <= nDir < 6) else { // riporto a indice griglia nDir -= 3 ; // Controllo l'ammissibilità della griglia if ( nDir < 0 || nDir > 2) return false ; // Verifiche sugli indici if ( nPos1 < 0 || nPos1 >= int( m_nNx[nDir]) || nPos2 < 0 || nPos2 >= int( m_nNy[nDir])) return false ; int nPos = nPos1 + nPos2 * m_nNx[nDir] ; if ( nPos < 0 || nPos >= int( m_Values[nDir].size())) return false ; // Calcolo coordinate punto double dX = m_dStep * ( 0.5 + nPos1) ; double dY = m_dStep * ( 0.5 + nPos2) ; // Determino il punto di partenza sulla griglia Point3d ptP = m_MapFrame[nDir].Orig() + dX * m_MapFrame[nDir].VersX() + dY * m_MapFrame[nDir].VersY() ; // Creo le polilinee for ( int j = 1 ; j < int( m_Values[nDir][nPos].size()) ; j += 2) { // aggiungo polilinea a lista lstPL.emplace_back() ; // calcolo e inserisco punto inizio spillone Point3d ptQ = ptP + m_Values[nDir][nPos][j-1].dZVal * m_MapFrame[nDir].VersZ() ; lstPL.back().AddUPoint( 0, ptQ) ; // calcolo e inserisco punto su termine sua normale Vector3d vtV = m_Values[nDir][nPos][j-1].vtN ; vtV.ToGlob( m_MapFrame[0]) ; lstPL.back().AddUPoint( 1, ptQ + vtV * m_dStep / 4) ; // aggiungo polilinea a lista lstPL.emplace_back() ; // calcolo e inserisco punto fine spillone Point3d ptR = ptP + m_Values[nDir][nPos][j].dZVal * m_MapFrame[nDir].VersZ() ; lstPL.back().AddUPoint( 0, ptR) ; // calcolo e inserisco punto su termine sua normale Vector3d vtW = m_Values[nDir][nPos][j].vtN ; vtW.ToGlob( m_MapFrame[0]) ; lstPL.back().AddUPoint( 1, ptR + vtW * m_dStep / 4) ; } return true ; } } //---------------------------------------------------------------------------- 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) ; } } } //else { // // //std::vector vecTria ; // //vecTria.resize( int( m_BlockToUpdate.size())) ; // //for ( int i = 0 ; i < int( m_BlockToUpdate.size()) ; ++ i) { // // //if ( m_BlockToUpdate[i]) // // ExtMarchingCubes( i, lstTria, vecTria[i]) ; } // TriHolder triHold ; // ExtMarchingCubes( 0, lstTria, triHold) ; // FlipEdges( triHold) ; // // for ( int i = 0 ; i < int( triHold.size()) ; ++ i) // for ( int j = 0 ; j < int( triHold[i].vecTria.size()) ; ++ j) // lstTria.emplace_back( triHold[i].vecTria[j]) ; //} else MarchingCubes( lstTria) ; return true ; } //---------------------------------------------------------------------------- bool VolZmap::GetBlockTriangles( int nBlock, TRIA3DLIST& lstTria) const { // Controllo sulla validità del blocco if ( nBlock < 0 || nBlock >= int( m_BlockToUpdate.size())) return false ; // Caso di singola mappa if ( m_nMapNum == 1) { const int MAX_DIM_CHUNK = 128 ; // Calcolo posizione del blocco nella griglia int nIBlock = nBlock % int( m_nFracLin[0]) ; int nJBlock = nBlock / int( m_nFracLin[0]) ; // Calcolo limiti per l'indice i int nStartI = nIBlock * 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) ; int nEndJ = ( nJBlock + 1 == int( m_nFracLin[1]) ? int( m_nNy[0]) : ( nJBlock + 1) * int( m_nDexNumPBlock)) ; // Ciclo su i e j for ( int i = nStartI ; i < nEndI ; i += MAX_DIM_CHUNK) { int nDimChunkX = min( MAX_DIM_CHUNK, nEndI - i) ; for ( int j = nStartJ ; j < nEndJ ; j += MAX_DIM_CHUNK) { int nDimChunkY = min( MAX_DIM_CHUNK, nEndJ - j) ; GetChunkPrisms( i, j, nDimChunkX, nDimChunkY, MAX_DIM_CHUNK, lstTria) ; } } } // Caso con tre mappe else { // Calcolo i limiti sugli indici dei voxel del blocco // Vettore indici i,j,k del blocco int nIJK[3] ; GetBlockIJK( nIJK, nBlock) ; // Vettore limiti sugli indici dei voxel del blocco int LimitIndexes[6] ; GetBlockLimitsIJK( nIJK, LimitIndexes) ; TriHolder triHold ; ExtMarchingCubes( LimitIndexes, lstTria, triHold) ; FlipEdges( triHold) ; // Valuto se esistono voxel contenenti sharp feature if ( ! triHold.size()) ; else { // Cerco fra i voxel con sharp feature quelli // di frontiera e quelli ad essi adiacenti. // Ciclo sui voxel con sharp feature for ( size_t t = 0 ; t < triHold.size() ; ++ t) { // Voxel di frontiera if ( triHold[t].i == LimitIndexes[0] || triHold[t].i == LimitIndexes[1] - 1 || triHold[t].j == LimitIndexes[2] || triHold[t].j == LimitIndexes[3] - 1 || triHold[t].k == LimitIndexes[4] || triHold[t].k == LimitIndexes[5] - 1) { // Ciclo sui voxel adiacenti for ( int nI = triHold[t].i - 1 ; nI <= triHold[t].i + 1 ; ++ nI) { for ( int nJ = triHold[t].j - 1 ; nJ <= triHold[t].j + 1 ; ++ nJ) { for ( int nK = triHold[t].k - 1 ; nK <= triHold[t].k + 1 ; ++ nK) { // Voxel adiacente appartenente allo // Zmap e non appartenente al blocco if ( ( nI >= -1 && nI <= int( m_nNx[0]) - 1 && nJ >= -1 && nJ <= int( m_nNy[0]) - 1 && nK >= -1 && nK <= int( m_nNy[1]) - 1) && ! ( nI >= LimitIndexes[0] && nI < LimitIndexes[1] && nJ >= LimitIndexes[2] && nJ < LimitIndexes[3] && nK >= LimitIndexes[4] && nK < LimitIndexes[5])) { std::vector AdjVoxel ; AdjVoxel.emplace_back( nI) ; AdjVoxel.emplace_back( nJ) ; AdjVoxel.emplace_back( nK) ; TriHolder NewTriHold ; ExtMarchingCubes( AdjVoxel, NewTriHold) ; if ( NewTriHold.size()) FlipEdgesLocalFlipEdges( triHold[t], NewTriHold[0]) ; } } } } } } } for ( size_t t1 = 0 ; t1 < triHold.size() ; ++ t1) for ( size_t t2 = 0 ; t2 < triHold[t1].vCompoTria.size() ; ++ t2) for ( size_t t3 = 0 ; t3 < triHold[t1].vCompoTria[t2].size() ; ++ t3) lstTria.emplace_back( triHold[t1].vCompoTria[t2][t3]) ; } m_BlockToUpdate[nBlock] = false ; return true ; } //---------------------------------------------------------------------------- bool VolZmap::GetBlockInfo( std::vector & bModified) const { bModified = m_BlockToUpdate ; return true ; } //---------------------------------------------------------------------------- int VolZmap::GetBlockCount( void) const { return m_nNumBlock ; } //---------------------------------------------------------------------------- bool VolZmap::GetChunkPrisms( int nPos1, int nPos2, int nDim1, int nDim2, int nDimChk, TRIA3DLIST& lstTria) const { // determino se è un semplice parallelepipedo bool bIsSimple = true ; double dBotZ ; double dTopZ ; for ( int i = 0 ; i < nDim1 && bIsSimple ; ++ i) { for ( int j = 0 ; j < nDim2 && bIsSimple ; ++ j) { int nPos = ( nPos1 + i) + ( nPos2 + j) * m_nNx[0] ; if ( nPos > int( m_nDim[0]) || int( m_Values[0][nPos].size()) != 2) bIsSimple = false ; else if ( i == 0 && j == 0) { dBotZ = m_Values[0][nPos][0].dZVal ; dTopZ = m_Values[0][nPos][1].dZVal ; } else if ( abs( m_Values[0][nPos][0].dZVal - dBotZ) > EPS_SMALL || abs( m_Values[0][nPos][1].dZVal - dTopZ) > EPS_SMALL) bIsSimple = false ; } } // se semplice parallelepipedo if ( bIsSimple) { CalcChunkPrisms( nPos1, nPos2, nDim1, nDim2, lstTria) ; } // se chunk di dimensioni accettabili else if ( nDimChk >= 4) { int nNewDimChk = nDimChk / 2 ; for ( int i = nPos1 ; i < int( nPos1 + nDim1) ; i += nNewDimChk) { int nDimChunkX = min( nNewDimChk, int( nPos1 + nDim1) - i) ; for ( int j = nPos2 ; j < int( nPos2 + nDim2) ; j += nNewDimChk) { int nDimChunkY = min( nNewDimChk, int( nPos2 + nDim2) - j) ; GetChunkPrisms( i, j, nDimChunkX, nDimChunkY, nNewDimChk, lstTria) ; } } } // altrimenti else { // elaboro ogni singolo dexel for ( int i = 0 ; i < nDim1 ; ++ i) { for ( int j = 0 ; j < nDim2 ; ++ j) { CalcDexelPrisms( nPos1 + i, nPos2 + j, lstTria) ; } } } return true ; } //---------------------------------------------------------------------------- bool VolZmap::CalcChunkPrisms( int nPos1, int nPos2, int nDim1, int nDim2, TRIA3DLIST& lstTria) const { // verifiche sugli indici if ( nPos1 < 0 || nPos1 + nDim1 > int( m_nNx[0]) || nPos2 < 0 || nPos2 + nDim2 > int( m_nNy[0])) return false ; int nPos = nPos1 + nPos2 * m_nNx[0] ; if ( nPos < 0 || nPos >= int( m_nDim[0])) return false ; // calcolo coordinate punti double dX = m_dStep * nPos1 ; double dY = m_dStep * nPos2 ; Point3d ptP1 = m_MapFrame[0].Orig() + dX * m_MapFrame[0].VersX() + dY * m_MapFrame[0].VersY() ; Point3d ptP2 = ptP1 + nDim1 * m_dStep * m_MapFrame[0].VersX() ; Point3d ptP3 = ptP2 + nDim2 * m_dStep * m_MapFrame[0].VersY() ; Point3d ptP4 = ptP1 + nDim2 * m_dStep * m_MapFrame[0].VersY() ; // creo le facce sopra e sotto Vector3d vtDZt = m_Values[0][nPos][1].dZVal * m_MapFrame[0].VersZ() ; Vector3d vtDZb = m_Values[0][nPos][0].dZVal * m_MapFrame[0].VersZ() ; // faccia superiore P1t->P2t->P3t->P4t : sempre visibile lstTria.emplace_back() ; lstTria.back().Set( ptP1 + vtDZt, ptP2 + vtDZt, ptP3 + vtDZt, m_MapFrame[0].VersZ()) ; lstTria.emplace_back() ; lstTria.back().Set( ptP3 + vtDZt, ptP4 + vtDZt, ptP1 + vtDZt, m_MapFrame[0].VersZ()) ; // faccia inferiore P1b->P4b->P3b->P2b : sempre visibile lstTria.emplace_back() ; lstTria.back().Set( ptP1 + vtDZb, ptP4 + vtDZb, ptP3 + vtDZb, - m_MapFrame[0].VersZ()) ; lstTria.emplace_back() ; lstTria.back().Set( ptP3 + vtDZb, ptP2 + vtDZb, ptP1 + vtDZb, - m_MapFrame[0].VersZ()) ; // creo le facce laterali for ( int j = 0 ; j < nDim2 ; ++ j) { int nPosD = nPos + nDim1 - 1 + j * m_nNx[0] ; int nPosEst = ( nPos1 + nDim1 - 1 < int( m_nNx[0] - 1) ? nPosD + 1 : - 1) ; Point3d ptP2D = ptP2 + j * m_dStep * m_MapFrame[0].VersY() ; Point3d ptP3D = ptP2D + m_dStep * m_MapFrame[0].VersY() ; AddDexelSideFace( nPosD, nPosEst, ptP2D, ptP3D, m_MapFrame[0].VersZ(), m_MapFrame[0].VersX(), lstTria) ; } for ( int i = 0 ; i < nDim1 ; ++ i) { int nPosD = nPos + ( nDim2 - 1) * m_nNx[0] + i ; int nPosNord = ( nPos2 + nDim2 - 1 < int( m_nNy[0] - 1) ? nPosD + m_nNx[0] : - 1) ; Point3d ptP4D = ptP4 + i * m_dStep * m_MapFrame[0].VersX() ; Point3d ptP3D = ptP4D + m_dStep * m_MapFrame[0].VersX() ; AddDexelSideFace( nPosD, nPosNord, ptP3D, ptP4D, m_MapFrame[0].VersZ(), m_MapFrame[0].VersY(), lstTria) ; } for ( int j = 0 ; j < nDim2 ; ++ j) { int nPosD = nPos + j * m_nNx[0] ; int nPosWest = ( nPos1 > 0 ? nPosD - 1 : - 1) ; Point3d ptP1D = ptP1 + j * m_dStep * m_MapFrame[0].VersY() ; Point3d ptP4D = ptP1D + m_dStep * m_MapFrame[0].VersY() ; AddDexelSideFace( nPosD, nPosWest, ptP4D, ptP1D, m_MapFrame[0].VersZ(), - m_MapFrame[0].VersX(), lstTria) ; } for ( int i = 0 ; i < nDim1 ; ++ i) { int nPosD = nPos + i ; int nPosSud = ( nPos2 > 0 ? nPosD - m_nNx[0] : - 1) ; Point3d ptP1D = ptP1 + i * m_dStep * m_MapFrame[0].VersX() ; Point3d ptP2D = ptP1D + m_dStep * m_MapFrame[0].VersX() ; AddDexelSideFace( nPosD, nPosSud, ptP1D, ptP2D, m_MapFrame[0].VersZ(), - m_MapFrame[0].VersY(), lstTria) ; } return true ; } //---------------------------------------------------------------------------- bool VolZmap::CalcDexelPrisms( int nPos1, int nPos2, TRIA3DLIST& lstTria) const { // verifiche sugli indici if ( nPos1 < 0 || nPos1 >= int( m_nNx[0]) || nPos2 < 0 || nPos2 >= int( m_nNy[0])) return false ; int nPos = nPos1 + nPos2 * m_nNx[0] ; if ( nPos < 0 || nPos >= int( m_nDim[0])) return false ; // calcolo coordinate punto double dX = m_dStep * nPos1 ; double dY = m_dStep * nPos2 ; Point3d ptP1 = m_MapFrame[0].Orig() + dX * m_MapFrame[0].VersX() + dY * m_MapFrame[0].VersY() ; Point3d ptP2 = ptP1 + m_dStep * m_MapFrame[0].VersX() ; Point3d ptP3 = ptP2 + m_dStep * m_MapFrame[0].VersY() ; Point3d ptP4 = ptP1 + m_dStep * m_MapFrame[0].VersY() ; // creo le facce sopra e sotto di ogni intervallo (sempre visibili) for ( int i = 1 ; i < int( m_Values[0][nPos].size()) ; i += 2) { Vector3d vtDZt = m_Values[0][nPos][i].dZVal * m_MapFrame[0].VersZ() ; Vector3d vtDZb = m_Values[0][nPos][i-1].dZVal * m_MapFrame[0].VersZ() ; // faccia superiore P1t->P2t->P3t->P4t : sempre visibile lstTria.emplace_back() ; lstTria.back().Set( ptP1 + vtDZt, ptP2 + vtDZt, ptP3 + vtDZt, m_MapFrame[0].VersZ()) ; lstTria.emplace_back() ; lstTria.back().Set( ptP3 + vtDZt, ptP4 + vtDZt, ptP1 + vtDZt, m_MapFrame[0].VersZ()) ; // faccia inferiore P1b->P4b->P3b->P2b : sempre visibile lstTria.emplace_back() ; lstTria.back().Set( ptP1 + vtDZb, ptP4 + vtDZb, ptP3 + vtDZb, - m_MapFrame[0].VersZ()) ; lstTria.emplace_back() ; lstTria.back().Set( ptP3 + vtDZb, ptP2 + vtDZb, ptP1 + vtDZb, - m_MapFrame[0].VersZ()) ; } // creo le facce laterali int nPosEst = ( nPos1 < int( m_nNx[0] - 1) ? nPos + 1 : - 1) ; AddDexelSideFace( nPos, nPosEst, ptP2, ptP3, m_MapFrame[0].VersZ(), m_MapFrame[0].VersX(), lstTria) ; int nPosNord = ( nPos2 < int( m_nNy[0] - 1) ? nPos + m_nNx[0] : - 1) ; AddDexelSideFace( nPos, nPosNord, ptP3, ptP4, m_MapFrame[0].VersZ(), m_MapFrame[0].VersY(), lstTria) ; int nPosWest = ( nPos1 > 0 ? nPos - 1 : - 1) ; AddDexelSideFace( nPos, nPosWest, ptP4, ptP1, m_MapFrame[0].VersZ(), - m_MapFrame[0].VersX(), lstTria) ; int nPosSud = ( nPos2 > 0 ? nPos - m_nNx[0] : - 1) ; AddDexelSideFace( nPos, nPosSud, ptP1, ptP2, m_MapFrame[0].VersZ(), - m_MapFrame[0].VersY(), lstTria) ; return true ; } //---------------------------------------------------------------------------- bool VolZmap::AddDexelSideFace( int nPos, int nPosAdj, const Point3d& ptP, const Point3d& ptQ, const Vector3d& vtZ, const Vector3d& vtNorm, TRIA3DLIST& lstTria) const { Intervals intFace ; for ( int i = 1 ; i < int( m_Values[0][nPos].size()) ; i += 2) intFace.Add( m_Values[0][nPos][i-1].dZVal, m_Values[0][nPos][i].dZVal) ; if ( nPosAdj > 0) { for ( int i = 1 ; i < int( m_Values[0][nPosAdj].size()) ; i += 2) intFace.Subtract( m_Values[0][nPosAdj][i-1].dZVal, m_Values[0][nPosAdj][i].dZVal) ; } double dMin, dMax ; bool bFound = intFace.GetFirst( dMin, dMax) ; while ( bFound) { Vector3d vtDZt = dMax * vtZ ; Vector3d vtDZb = dMin * vtZ ; lstTria.emplace_back() ; lstTria.back().Set( ptP + vtDZb, ptQ + vtDZb, ptQ + vtDZt, vtNorm) ; lstTria.emplace_back() ; lstTria.back().Set( ptQ + vtDZt, ptP + vtDZt, ptP + vtDZb, vtNorm) ; bFound = intFace.GetNext( dMin, dMax) ; } 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())) 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]) ; int nKBlock = ( nBlock / int( m_nFracLin[0] * m_nFracLin[1])) ; // 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)) ; // 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)) ; // 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)) ; // Ciclo su tutti i voxel dello Zmap for ( int i = nStartI ; i < nEndI ; ++ i) { for ( int j = nStartJ ; j < nEndJ ; ++ j) { for ( int k = nStartK ; k < nEndK ; ++ 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 } } ; // Ciclo sui segmenti Point3d ptIntPoint[12] ; 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] ; Triangle3d CurrentTriangle ; Vector3d vtN = ( ptIntPoint[i1] - ptIntPoint[i0]) ^ ( ptIntPoint[i2] - ptIntPoint[i1]) ; vtN.Normalize() ; vtN.ToGlob( m_MapFrame[0]) ; // Il triangolo è pronto CurrentTriangle.Set( ptIntPoint[i0], ptIntPoint[i1], ptIntPoint[i2], vtN) ; // Aggiungo triangolo lstTria.emplace_back( CurrentTriangle) ; } } } } return true ; } //---------------------------------------------------------------------------- bool VolZmap::ExtMarchingCubes( const int nLimits[], TRIA3DLIST& lstTria, TriHolder& triHold) const { Point3d ptMapOrig = m_MapFrame[0].Orig() ; // 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) { // 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} } ; int nIndex = 0 ; // Classificazione dei vertici: interni o esterni al materiale 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 }, { 3, 2 }, { 0, 3 }, { 4, 5 }, { 5, 6 }, { 7, 6 }, { 4, 7 }, { 0, 4 }, { 1, 5 }, { 2, 6 }, { 3, 7 } } ; // Arrey di strutture punto di intersezione // e normale alla superficie in esso. VectorField VecField[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] ; bool bN1 = ( ( nIndex & ( 1 << n1)) != 0) ; // Determino con precisione il punto di intersezione sullo spigolo. NewIntersPos( IndexCorner[n1], IndexCorner[n2], bN1, VecField[EdgeIndex].ptInt, VecField[EdgeIndex].vtNorm) ; VecField[EdgeIndex].ptInt.ToGlob( m_MapFrame[0]) ; VecField[EdgeIndex].vtNorm.ToGlob( m_MapFrame[0]) ; } // Determino il numero di componenti connesse int nComponents = TriangleTableEn[nIndex][1][0] ; // Serve nel ciclo che salva i punti e vettori di // una componente nell'arrey di compentenza: La tabella // fornisce numero di componenti, numero di vertici per // componenti per OGNUNA delle componenti e in fine // elenca i vertici della prima componente, seguiti da quelli // della seconda e così via. int nTableOffset = nComponents ; // Numero di feature nel voxel: al più vi è una feature per componente connessa int nFeatureInVoxel = 0 ; // Ciclo sulle componenti for ( int nCompCount = 1 ; nCompCount <= nComponents ; ++ nCompCount) { // Numero vertici per componenti int nVertComp = TriangleTableEn[nIndex][1][nCompCount] ; // Vettore di Vector3d VectorField CompoVert[12] ; // Riempio il vettore for ( int nVertCount = 0 ; nVertCount < nVertComp ; ++ nVertCount) // Nota che il primo elemento dell'array (0-esimo) non viene inizializzato CompoVert[nVertCount] = VecField[TriangleTableEn[nIndex][1][nVertCount + nTableOffset + 1]] ; int nFeatureType ; // Valuto le relazioni reciproche fra le normali e // se i punti sono su un piano di griglia. bool bNormal = TestOnNormal( CompoVert, nVertComp, nFeatureType) ; bool bExt = bNormal ; // Extended MC if ( bExt) { // Passo al sistema di riferimento del baricentro Point3d ptGravityCenter( 0, 0, 0) ; for ( int ni = 0 ; ni < nVertComp ; ++ ni) ptGravityCenter += CompoVert[ni].ptInt ; ptGravityCenter /= nVertComp ; Vector3d vtO = ptGravityCenter - ORIG ; Point3d ptTrasf[12] ; for ( int ni = 0 ; ni < nVertComp ; ++ ni) ptTrasf[ni] = CompoVert[ni].ptInt - vtO ; // Preparo le matrici per il sistema typedef Eigen::Matrix dSystemMatrix ; typedef Eigen::Matrix dSystemVector ; typedef Eigen::JacobiSVD DecomposerSVD ; dSystemMatrix dMatrixN, dMatrixU, dMatrixV ; dSystemVector dKnownVector, dUnknownVector, dSingularValue ; dMatrixN.resize( nVertComp, 3) ; dKnownVector.resize( nVertComp, 1) ; dUnknownVector.resize( 3, 1) ; // Studio del caso 4 punti su un piano int nEqual = 0 ; int nPosD ; Vector3d vtD, vtE ; if ( nVertComp == 4 && nFeatureType == 2) { int nPosEq ; for ( int ni = 0 ; ni < 2 ; ++ ni) { for ( int nj = ni + 1 ; nj < nVertComp ; ++ nj) { if ( AreSameVectorApprox( CompoVert[ni].vtNorm, CompoVert[nj].vtNorm)) { nEqual ++ ; nPosEq = ni ; } } if ( nEqual == 2) break ; else nEqual = 0 ; } if ( nEqual == 2) { for ( int ni = 0 ; ni < nVertComp ; ++ ni) if ( ! AreSameVectorApprox( CompoVert[ni].vtNorm, CompoVert[nPosEq].vtNorm)) { nPosD = ni ; vtD = CompoVert[ni].vtNorm ; vtE = CompoVert[nPosEq].vtNorm ; } } } double dDot = abs( ( CompoVert[1].ptInt - CompoVert[0].ptInt) * ( ( CompoVert[2].ptInt - CompoVert[1].ptInt) ^ ( CompoVert[3].ptInt - CompoVert[2].ptInt))) ; // Caso quattro punti su un piano if ( nVertComp == 4 && nEqual == 2 && dDot < EPS_SMALL) { for ( int ni = 0 ; ni < nVertComp ; ++ ni) { if ( ni != nPosD) { dMatrixN( ni, 0) = CompoVert[ni].vtNorm.x ; dMatrixN( ni, 1) = CompoVert[ni].vtNorm.y ; dMatrixN( ni, 2) = CompoVert[ni].vtNorm.z ; dKnownVector( ni) = CompoVert[ni].vtNorm * ( ptTrasf[ni] - ORIG) ; } else { dMatrixN( ni, 0) = vtE.x ; dMatrixN( ni, 1) = vtE.y ; dMatrixN( ni, 2) = vtE.z ; dKnownVector( ni) = vtE * ( ptTrasf[ni] - ORIG) ; } } } // caso generale else { for ( int ni = 0 ; ni < nVertComp ; ++ ni) { dMatrixN( ni, 0) = CompoVert[ni].vtNorm.x ; dMatrixN( ni, 1) = CompoVert[ni].vtNorm.y ; dMatrixN( ni, 2) = CompoVert[ni].vtNorm.z ; dKnownVector( ni) = CompoVert[ni].vtNorm * ( ptTrasf[ni] - ORIG) ; } } DecomposerSVD svd( dMatrixN, Eigen::ComputeThinU | Eigen::ComputeThinV) ; dMatrixU = svd.matrixU() ; dMatrixV = svd.matrixV() ; dSingularValue = svd.singularValues() ; double s0 = dSingularValue( 0) ; double s1 = dSingularValue( 1) ; double s2 = dSingularValue( 2) ; if ( nFeatureType == 2 ) dSingularValue( 2) = 0 ; // Back substitution: risolvo il sistema USV*x = b // Calcolo U^T b double vTemp[3] ; for ( int ni = 0 ; ni < 3 ; ++ ni) { double s = 0 ; if ( abs( dSingularValue( ni)) > EPS_SMALL) { for ( int nj = 0 ; nj < nVertComp ; ++ nj) s += dMatrixU( nj, ni) * dKnownVector( nj) ; s /= dSingularValue( ni) ; } vTemp[ni] = s ; } // Moltiplico per V for ( int ni = 0 ; ni < 3 ; ++ ni) { double s = 0 ; for ( int nj = 0 ; nj < 3 ; ++ nj) s += dMatrixV( ni, nj) * vTemp[nj] ; dUnknownVector( ni) = s ; } // Limito la feature entro una distanza di 3 // volte la diagonale del voxel dal baricentro. Vector3d vtFeature( dUnknownVector( 0), dUnknownVector( 1), dUnknownVector( 2)) ; double dDistFeature = vtFeature.Len() ; const double dMaxDist = sqrt( 3) * m_dStep ; if ( dDistFeature > dMaxDist) { // 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] ; Triangle3d CurrentTriangle ; // Il triangolo è pronto CurrentTriangle.Set( VecField[i0].ptInt, VecField[i1].ptInt, VecField[i2].ptInt) ; CurrentTriangle.Validate( true) ; // Se il triangolo non è degenere lo aggiungo alla lista if ( ! ( AreSamePointApprox( CurrentTriangle.GetP( 0), CurrentTriangle.GetP( 1)) || AreSamePointApprox( CurrentTriangle.GetP( 1), CurrentTriangle.GetP( 2)) || AreSamePointApprox( CurrentTriangle.GetP( 2), CurrentTriangle.GetP( 0)))) lstTria.emplace_back( CurrentTriangle) ; } continue ; } // Esprimo la soluzione nel sistema di riferimento dello z-map. Point3d ptSol( dUnknownVector( 0) + vtO.x, dUnknownVector( 1) + vtO.y, dUnknownVector( 2) + vtO.z) ; Triangle3d CurrentTriangle ; TRIA3DVECTOR triContainer ; bool bInvNormal = false ; for ( int ni = 0 ; ni < nVertComp ; ++ ni) { // Il triangolo è pronto int nj = ( ni + 1 < nVertComp) ? ni + 1 : 0 ; CurrentTriangle.Set( ptSol, CompoVert[nj].ptInt, CompoVert[ni].ptInt) ; CurrentTriangle.Validate( true) ; // Test sulla degenerazione if ( CurrentTriangle.GetArea() < EPS_SMALL) continue ; // Controllo sull'inversione delle normali if ( CurrentTriangle.GetN() * CompoVert[nj].vtNorm < - 0.7 || CurrentTriangle.GetN() * CompoVert[ni].vtNorm < - 0.7) { //string sInfo ; //sInfo = "TriaN=" + ToString( CurrentTriangle.GetN()) + // "VertJ=" + ToString( CompoVert[nj].vtNorm) + // "VertI=" + ToString( CompoVert[ni].vtNorm) ; //LOG_INFO( GetEGkLogger(), sInfo.c_str()) //bInvNormal = true ; } // Aggiungo triangolo al vettore temporaneo triContainer.emplace_back( CurrentTriangle) ; } // Valuto normali int nContSize = triContainer.size() ; bool bPlane = true ; for ( int ni = 0 ; ni < nContSize - 1 ; ++ ni) { for ( int nj = ni + 1 ; nj < nContSize ; ++ nj) { Vector3d vtI = triContainer[ni].GetN() ; Vector3d vtJ = triContainer[nj].GetN() ; if ( ! AreSameVectorApprox( vtI, vtJ)) { bPlane = false ; break ; } } } if ( ! ( bPlane || bInvNormal)) { // Aggiorno il numero di feature. ++ nFeatureInVoxel ; // Se siamo in presenza della prima feature del // voxel, ridimensiono il vettore che contiene // la struttura dati del voxel. if ( nFeatureInVoxel == 1) { triHold.resize( triHold.size() + 1) ; // Questi dati dipendono solo dal voxel, // quindi sono validi 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 = 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 < nContSize ; ++ ni) triHold[nCurrent].vCompoTria[nOldFeatureNum].emplace_back( triContainer[ni]) ; } else { // 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] ; Triangle3d CurrentTriangle ; // Il triangolo è pronto CurrentTriangle.Set( VecField[i0].ptInt, VecField[i1].ptInt, VecField[i2].ptInt) ; CurrentTriangle.Validate( true) ; // Se il triangolo non è degenere lo aggiungo alla lista if ( ! ( AreSamePointApprox( CurrentTriangle.GetP( 0), CurrentTriangle.GetP( 1)) || AreSamePointApprox( CurrentTriangle.GetP( 1), CurrentTriangle.GetP( 2)) || AreSamePointApprox( CurrentTriangle.GetP( 2), CurrentTriangle.GetP( 0)))) lstTria.emplace_back( CurrentTriangle) ; } } } // Standard MC else { // 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] ; Triangle3d CurrentTriangle ; // Il triangolo è pronto CurrentTriangle.Set( VecField[i0].ptInt, VecField[i1].ptInt, VecField[i2].ptInt) ; CurrentTriangle.Validate( true) ; // Se il triangolo non è degenere lo aggiungo alla lista if ( ! ( AreSamePointApprox( CurrentTriangle.GetP( 0), CurrentTriangle.GetP( 1)) || AreSamePointApprox( CurrentTriangle.GetP( 1), CurrentTriangle.GetP( 2)) || AreSamePointApprox( CurrentTriangle.GetP( 2), CurrentTriangle.GetP( 0)))) lstTria.emplace_back( CurrentTriangle) ; } } nTableOffset += nVertComp ; } } } } return true ; } //---------------------------------------------------------------------------- bool VolZmap::ExtMarchingCubes( const std::vector VoxelsIndexes, TriHolder& triHold) const { // Controllo sulla validità della dimensione del // vector: essa deve essere un multiplo di 3 if ( ( VoxelsIndexes.size() % 3)) return false ; Point3d ptMapOrig = m_MapFrame[0].Orig() ; // Ciclo sui voxel for ( size_t t = 0 ; t < VoxelsIndexes.size() - 2 ; t += 3) { // Indici del voxel int i = VoxelsIndexes[t] ; int j = VoxelsIndexes[t + 1] ; int k = VoxelsIndexes[t + 2] ; // 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} } ; int nIndex = 0 ; // Classificazione dei vertici: interni o esterni al materiale 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 }, { 0, 3 }, { 4, 5 }, { 5, 6 }, { 6, 7 }, { 4, 7 }, { 0, 4 }, { 1, 5 }, { 2, 6 }, { 3, 7 } } ; // Arrey di strutture punto di intersezione // e normale alla superficie in esso. VectorField VecField[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], VecField[EdgeIndex].ptInt, VecField[EdgeIndex].vtNorm) ; VecField[EdgeIndex].ptInt.ToGlob( m_MapFrame[0]) ; VecField[EdgeIndex].vtNorm.ToGlob( m_MapFrame[0]) ; } // Determino il numero di componenti connesse int nComponents = TriangleTableEn[nIndex][1][0] ; // Serve nel ciclo che salva i punti e vettori di // una componente nell'arrey di compentenza: La tabella // fornisce numero di componenti, numero di vertici per // componenti per OGNUNA delle componenti e in fine // elenca i vertici della prima componente, seguiti da quelli // della seconda e così via. int nTableOffset = nComponents ; // Numero di feature nel voxel: al più vi // è una feature per componente connessa int nFeatureInVoxel = 0 ; // Ciclo sulle componenti for ( int nCompCount = 1 ; nCompCount <= nComponents ; ++ nCompCount) { // Numero vertici per componenti int nVertComp = TriangleTableEn[nIndex][1][nCompCount] ; // Vettore di Vector3d VectorField CompoVert[12] ; // Riempio il vettore for ( int nVertCount = 0 ; nVertCount < nVertComp ; ++ nVertCount) // Nota che il primo elemento dell'array (0-esimo) non viene inizializzato CompoVert[nVertCount] = VecField[TriangleTableEn[nIndex][1][nVertCount + nTableOffset + 1]] ; int nFeatureType ; // Valuto le relazioni reciproche fra le normali bool bExt = TestOnNormal( CompoVert, nVertComp, nFeatureType) ; // Extended MC if ( bExt) { // Aggiorno il numero di feature. ++ nFeatureInVoxel ; // Se siamo in presenza della prima feature del // voxel, ridimensiono il vettore che contiene // la struttura dati del voxel. if ( nFeatureInVoxel == 1) { triHold.resize( triHold.size() + 1) ; // Questi dati dipendono solo dal voxel, // quindi sono validi 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 ; // Passo al sistema di riferimento del baricentro Point3d ptGravityCenter( 0, 0, 0) ; for ( int ni = 0 ; ni < nVertComp ; ++ ni) ptGravityCenter += CompoVert[ni].ptInt ; ptGravityCenter /= nVertComp ; Vector3d vtO = ptGravityCenter - ORIG ; Point3d ptTrasf[12] ; for ( int ni = 0 ; ni < nVertComp ; ++ ni) ptTrasf[ni] = CompoVert[ni].ptInt - vtO ; // Preparo le matrici per il sistema typedef Eigen::Matrix dSystemMatrix ; typedef Eigen::Matrix dSystemVector ; typedef Eigen::JacobiSVD DecomposerSVD ; dSystemMatrix dMatrixN, dMatrixU, dMatrixV ; dSystemVector dKnownVector, dUnknownVector, dSingularValue ; dMatrixN.resize( nVertComp, 3) ; dKnownVector.resize( nVertComp, 1) ; dUnknownVector.resize( 3, 1) ; // Studio del caso 4 punti su un piano int nEqual = 0 ; int nPosD ; Vector3d vtD, vtE ; if ( nVertComp == 4 && nFeatureType == 2) { int nPosEq ; for ( int ni = 0 ; ni < 2 ; ++ ni) { for ( int nj = ni + 1 ; nj < nVertComp ; ++ nj) { if ( AreSameVectorApprox( CompoVert[ni].vtNorm, CompoVert[nj].vtNorm)) { nEqual ++ ; nPosEq = ni ; } } if ( nEqual == 2) break ; else nEqual = 0 ; } if ( nEqual == 2) { for ( int ni = 0 ; ni < nVertComp ; ++ ni) if ( ! AreSameVectorApprox( CompoVert[ni].vtNorm, CompoVert[nPosEq].vtNorm)) { nPosD = ni ; vtD = CompoVert[ni].vtNorm ; vtE = CompoVert[nPosEq].vtNorm ; } } } double dDot = abs( ( CompoVert[1].ptInt - CompoVert[0].ptInt) * ( ( CompoVert[2].ptInt - CompoVert[1].ptInt) ^ ( CompoVert[3].ptInt - CompoVert[2].ptInt))) ; // Caso quattro punti su un piano if ( nVertComp == 4 && nEqual == 2 && dDot < EPS_SMALL) { for ( int ni = 0 ; ni < nVertComp ; ++ ni) { if ( ni != nPosD) { dMatrixN( ni, 0) = CompoVert[ni].vtNorm.x ; dMatrixN( ni, 1) = CompoVert[ni].vtNorm.y ; dMatrixN( ni, 2) = CompoVert[ni].vtNorm.z ; dKnownVector( ni) = CompoVert[ni].vtNorm * ( ptTrasf[ni] - ORIG) ; } else { dMatrixN( ni, 0) = vtE.x ; dMatrixN( ni, 1) = vtE.y ; dMatrixN( ni, 2) = vtE.z ; dKnownVector( ni) = vtE * ( ptTrasf[ni] - ORIG) ; } } } // caso generale else { for ( int ni = 0 ; ni < nVertComp ; ++ ni) { dMatrixN( ni, 0) = CompoVert[ni].vtNorm.x ; dMatrixN( ni, 1) = CompoVert[ni].vtNorm.y ; dMatrixN( ni, 2) = CompoVert[ni].vtNorm.z ; dKnownVector( ni) = CompoVert[ni].vtNorm * ( ptTrasf[ni] - ORIG) ; } } DecomposerSVD svd( dMatrixN, Eigen::ComputeThinU | Eigen::ComputeThinV) ; dMatrixU = svd.matrixU() ; dMatrixV = svd.matrixV() ; dSingularValue = svd.singularValues() ; double s0 = dSingularValue( 0) ; double s1 = dSingularValue( 1) ; double s2 = dSingularValue( 2) ; if ( nFeatureType == 2 ) dSingularValue( 2) = 0 ; // Back substitution: risolvo il sistema USV*x = b // Calcolo U^T b double vTemp[3] ; for ( int ni = 0 ; ni < 3 ; ++ ni) { double s = 0 ; if ( abs( dSingularValue( ni)) > EPS_SMALL) { for ( int nj = 0 ; nj < nVertComp ; ++ nj) s += dMatrixU( nj, ni) * dKnownVector( nj) ; s /= dSingularValue( ni) ; } vTemp[ni] = s ; } // Moltiplico per V for ( int ni = 0 ; ni < 3 ; ++ ni) { double s = 0 ; for ( int nj = 0 ; nj < 3 ; ++ nj) s += dMatrixV( ni, nj) * vTemp[nj] ; dUnknownVector( ni) = s ; } // Limito la feature entro una distanza di 3 // volte la diagonale del voxel dal baricentro. Vector3d vtFeature( dUnknownVector( 0), dUnknownVector( 1), dUnknownVector( 2)) ; double dDistFeature = vtFeature.Len() ; const double dMaxDist = 2 * sqrt( 3) * m_dStep / 2 ; if ( dDistFeature > dMaxDist) { vtFeature = ( dMaxDist / dDistFeature) * vtFeature ; dUnknownVector( 0) = vtFeature.x ; dUnknownVector( 1) = vtFeature.y ; dUnknownVector( 2) = vtFeature.z ; } // Esprimo la soluzione nel sistema di riferimento dello z-map Point3d ptSol( dUnknownVector( 0) + vtO.x, dUnknownVector( 1) + vtO.y, dUnknownVector( 2) + vtO.z) ; // Aggiungo vertice della componente // connessa alla lista dei vertici. triHold[nCurrent].ptCompoVert.emplace_back( ptSol) ; int nOldFeatureNum = triHold[nCurrent].vCompoTria.size() ; int nNewFeatureNum = nOldFeatureNum + 1 ; // Aggiungo una componente per il vettore // dei triangoli della componente connessa. triHold[nCurrent].vCompoTria.resize( nNewFeatureNum) ; Triangle3d CurrentTriangle ; TRIA3DVECTOR triContainer ; for ( int ni = 0 ; ni < nVertComp - 1 ; ++ ni) { // Il triangolo è pronto CurrentTriangle.Set( ptSol, CompoVert[ni+1].ptInt, CompoVert[ni].ptInt) ; CurrentTriangle.Validate( true) ; // Test sulla degenerazione if ( CurrentTriangle.GetArea() < EPS_SMALL /*AreSamePointApprox( CurrentTriangle.GetP( 0), CurrentTriangle.GetP( 1)) || AreSamePointApprox( CurrentTriangle.GetP( 1), CurrentTriangle.GetP( 2)) || AreSamePointApprox( CurrentTriangle.GetP( 2), CurrentTriangle.GetP( 0))*/) continue ; // Aggiungo triangolo triHold[nCurrent].vCompoTria[nOldFeatureNum].emplace_back( CurrentTriangle) ; } // Ultimo triangolo CurrentTriangle.Set( ptSol, CompoVert[0].ptInt, CompoVert[nVertComp - 1].ptInt) ; CurrentTriangle.Validate( true) ; // Test sulla degenerazione if ( CurrentTriangle.GetArea() < EPS_SMALL /*AreSamePointApprox( CurrentTriangle.GetP( 0), CurrentTriangle.GetP( 1)) || AreSamePointApprox( CurrentTriangle.GetP( 1), CurrentTriangle.GetP( 2)) || AreSamePointApprox( CurrentTriangle.GetP( 2), CurrentTriangle.GetP( 0))*/) continue ; // Aggiungo ultimo triangolo triHold[nCurrent].vCompoTria[nOldFeatureNum].emplace_back( CurrentTriangle) ; } nTableOffset += nVertComp ; } } return true ; } //---------------------------------------------------------------------------- bool VolZmap::FlipEdges( TriHolder& triHold) const { // Numero di voxel in cui si presentano sharp feature int nVoxelNum = int( triHold.size()) ; // Ciclo su tali voxel for ( int n = 0 ; n < nVoxelNum ; ++ n) { for ( int m = n ; m < nVoxelNum ; ++ m) { // Voxel adiacenti o coincidenti if ( m == n || ( ( triHold[m].i < int( m_nNx[0]) && triHold[m].j < int( m_nNy[0]) && triHold[m].k < int( m_nNy[1])) && ( triHold[m].i == triHold[n].i + 1 || triHold[m].j == triHold[n].j + 1 || triHold[m].k == triHold[n].k + 1))) { // Numero delle componenti connesse nei due voxel int nNumCompoN = triHold[n].ptCompoVert.size() ; int nNumCompoM = triHold[m].ptCompoVert.size() ; int nCompoN = 0 ; // Ciclo sulle componenti for ( ; nCompoN < nNumCompoN ; ++ nCompoN) { int nCompoM = ( m == n ? nCompoN + 1 : 0) ; for ( ; nCompoM < nNumCompoM ; ++ nCompoM) { // Numero di triangoli per le componenti connesse int nNumN = int( triHold[n].vCompoTria[nCompoN].size()) ; int nNumM = int( triHold[m].vCompoTria[nCompoM].size()) ; for ( int triN = 0 ; triN < nNumN ; ++ triN) { bool bModified = false ; for ( int triM = 0 ; triM < nNumM ; ++ triM) { std::vector SharedIndex ; for ( int vertN = 0 ; vertN < 3 ; ++ vertN) { for ( int vertM = 0 ; vertM < 3 ; ++ vertM) { Point3d ptN = triHold[n].vCompoTria[nCompoN][triN].GetP( vertN) ; Point3d ptM = triHold[m].vCompoTria[nCompoM][triM].GetP( vertM) ; if ( SqDist( ptN, ptM) < EPS_SMALL * EPS_SMALL) { Point3d ptVertN = triHold[n].ptCompoVert[nCompoN] ; Point3d ptVertM = triHold[m].ptCompoVert[nCompoM] ; if ( ! ( AreSamePointApprox( ptN, ptVertN) || AreSamePointApprox( ptM, ptVertM))) { SharedIndex.emplace_back( vertN) ; SharedIndex.emplace_back( vertM) ; } } if ( SharedIndex.size() > 2) break ; } if ( SharedIndex.size() > 2) break ; } // Si deve operare la modifica dei triangoli if ( SharedIndex.size() > 2) { // Controllo sulle normali Point3d ptMFeature = triHold[m].ptCompoVert[nCompoM] ; Point3d ptNFeature = triHold[n].ptCompoVert[nCompoN] ; Vector3d vtFeature = ptMFeature - ptNFeature ; vtFeature.Normalize() ; double dDot = vtFeature * triHold[m].vCompoTria[nCompoM][triM].GetN() ; // Ulteriore controllo sulle normali Point3d ptP0 = ptNFeature ; Point3d ptP1 = triHold[n].vCompoTria[nCompoN][triN].GetP( SharedIndex[0]) ; Point3d ptP2 = triHold[n].vCompoTria[nCompoN][triN].GetP( SharedIndex[2]) ; Point3d ptP3 = ptMFeature ; double dPlane = ( ( ptP1 - ptP0) ^ ( ptP2 - ptP0)) * ( ptP3 - ptP0) ; Vector3d vtN = triHold[n].vCompoTria[nCompoN][triN].GetN() ; Vector3d vtM = triHold[m].vCompoTria[nCompoM][triM].GetN() ; double dDotNM = vtN * vtM ; int nProd = ( SharedIndex[2] - SharedIndex[0]) * ( SharedIndex[3] - SharedIndex[1]) ; // --- if ( nProd != 1 && nProd != - 2 && nProd != 4) { triHold[n].vCompoTria[nCompoN][triN].SetP( SharedIndex[0], triHold[m].ptCompoVert[nCompoM]) ; triHold[m].vCompoTria[nCompoM][triM].SetP( SharedIndex[3], triHold[n].ptCompoVert[nCompoN]) ; triHold[n].vCompoTria[nCompoN][triN].Validate( true) ; triHold[m].vCompoTria[nCompoM][triM].Validate( true) ; bModified = true ; break ; } } } if ( bModified) break ; } } } } } } return true ; } //---------------------------------------------------------------------------- bool VolZmap::FlipEdgesLocalFlipEdges( TriaStruct& triStrCurr, TriaStruct& triStrAdj) const { // Numero delle componenti connesse nei due voxel size_t nCurVoxCompNum = triStrCurr.ptCompoVert.size() ; size_t nAdjVoxCompNum = triStrAdj.ptCompoVert.size() ; // Ciclo sulle componenti connesse for ( size_t nCVCurComp = 0 ; nCVCurComp < nCurVoxCompNum ; ++ nCVCurComp) { for ( size_t nAVCurComp = 0 ; nAVCurComp < nAdjVoxCompNum ; ++ nAVCurComp) { // Numero di triangoli per le componenti connesse size_t nCVTriaNum = triStrCurr.vCompoTria[nCVCurComp].size() ; size_t nAVTriaNum = triStrAdj.vCompoTria[nAVCurComp].size() ; for ( size_t nCVT = 0 ; nCVT < nCVTriaNum ; ++ nCVT) { bool bModified = false ; for ( size_t nAVT = 0 ; nAVT < nAVTriaNum ; ++ nAVT) { std::vector SharedIndex ; for ( size_t nCurVoxVert = 0 ; nCurVoxVert < 3 ; ++ nCurVoxVert) { for ( size_t nAdjVoxVert = 0 ; nAdjVoxVert < 3 ; ++ nAdjVoxVert) { Point3d ptC = triStrCurr.vCompoTria[nCVCurComp][nCVT].GetP( nCurVoxVert) ; Point3d ptA = triStrAdj.vCompoTria[nAVCurComp][nAVT].GetP( nAdjVoxVert) ; if ( SqDist( ptC, ptA) < EPS_SMALL * EPS_SMALL) { SharedIndex.emplace_back( int( nCurVoxVert)) ; SharedIndex.emplace_back( int( nAdjVoxVert)) ; } if ( SharedIndex.size() > 2) break ; } if ( SharedIndex.size() > 2) break ; } // Si deve operare la modifica dei triangoli if ( SharedIndex.size() > 2) { // Controllo sulle normali Point3d ptCurFeature = triStrCurr.ptCompoVert[nCVCurComp] ; Point3d ptAdjFeature = triStrAdj.ptCompoVert[nAVCurComp] ; Vector3d vtFeature = ptAdjFeature - ptCurFeature ; vtFeature.Normalize() ; double dDot = vtFeature * triStrAdj.vCompoTria[nAVCurComp][nAVT].GetN() ; // Ulteriore controllo sulle normali Point3d ptP0 = ptCurFeature ; Point3d ptP1 = triStrCurr.vCompoTria[nCVCurComp][nCVT].GetP( SharedIndex[0]) ; Point3d ptP2 = triStrCurr.vCompoTria[nCVCurComp][nCVT].GetP( SharedIndex[2]) ; Point3d ptP3 = ptAdjFeature ; double dPlane = ( ( ptP1 - ptP0) ^ ( ptP2 - ptP0)) * ( ptP3 - ptP0) ; Vector3d vtN = triStrCurr.vCompoTria[nCVCurComp][nCVT].GetN() ; Vector3d vtM = triStrAdj.vCompoTria[nAVCurComp][nAVT].GetN() ; double dDotNM = vtN * vtM ; int nProd = ( SharedIndex[2] - SharedIndex[0]) * ( SharedIndex[3] - SharedIndex[1]) ; if ( nProd != 1 && nProd != - 2 && nProd != 4 && ( dDot < - EPS_SMALL || ( dPlane < EPS_SMALL && dDotNM < - 0.95))) { triStrCurr.vCompoTria[nCVCurComp][nCVT].SetP( SharedIndex[0], triStrAdj.ptCompoVert[nAVCurComp]) ; triStrAdj.vCompoTria[nAVCurComp][nAVT].SetP( SharedIndex[3], triStrCurr.ptCompoVert[nCVCurComp]) ; triStrCurr.vCompoTria[nCVCurComp][nCVT].Validate( true) ; triStrAdj.vCompoTria[nAVCurComp][nAVT].Validate( true) ; bModified = true ; break ; } } } if ( bModified) break ; } } } return true ; } //---------------------------------------------------------------------------- bool VolZmap::IsThereMat( int nI, int nJ, int nK) const { if ( nI == - 1 || nI == m_nNx[0] || nJ == - 1 || nJ == m_nNy[0] || nK == - 1 || nK == m_nNy[1]) return false ; double dEps = 2 * EPS_SMALL ; double dZ[3] ; dZ[0] = ( nK + 0.5) * m_dStep ; dZ[1] = ( nI + 0.5) * m_dStep ; dZ[2] = ( nJ + 0.5) * m_dStep ; int nCount = 0 ; for ( int nGrid = 0 ; nGrid < int ( m_nMapNum) ; ++ nGrid) { unsigned int nGrI, nGrJ ; if ( nGrid == 0) { nGrI = nI ; nGrJ = nJ ; } else if ( nGrid == 1) { nGrI = nJ ; nGrJ = nK ; } else { nGrI = nK ; nGrJ = nI ; } unsigned int nPos = nGrJ * m_nNx[nGrid] + nGrI ; size_t nDexSize = m_Values[nGrid][nPos].size() ; size_t nIndex = 0 ; while ( nIndex < nDexSize) { if ( dZ[nGrid] > m_Values[nGrid][nPos][nIndex].dZVal - dEps && dZ[nGrid] < m_Values[nGrid][nPos][nIndex + 1].dZVal + dEps) { ++ nCount ; break ; } nIndex += 2 ; } } return ( nCount == 3) ; } //---------------------------------------------------------------------------- bool VolZmap::IntersPos( int nVec1[], int nVec2[], Point3d& ptInt) const { if ( nVec1[0] != nVec2[0]) { ptInt.y = ( nVec1[1] + 0.5) * m_dStep ; ptInt.z = ( nVec1[2] + 0.5) * m_dStep ; int nMinI = min( nVec1[0], nVec2[0]) ; int nMaxI = max( nVec1[0], nVec2[0]) ; double dMinX = ( nMinI + 0.5) * m_dStep ; double dMaxX = ( nMaxI + 0.5) * m_dStep ; unsigned int nDexel = nVec1[2] * m_nNx[1] + nVec1[1] ; size_t nSize = m_Values[1][nDexel].size() ; bool bFound = false ; for ( size_t i = 0 ; i < nSize ; i += 2) { double dx1 = m_Values[1][nDexel][i].dZVal ; double dx2 = m_Values[1][nDexel][i+1].dZVal ; if ( dx1 < dMinX - EPS_SMALL && dx2 > dMinX - EPS_SMALL && dx2 < dMaxX + EPS_SMALL) { ptInt.x = dx2 ; bFound = true ; break ; } else if ( dx1 > dMinX - EPS_SMALL && dx1 < dMaxX + EPS_SMALL && dx2 > dMaxX + EPS_SMALL) { ptInt.x = dx1 ; bFound = true ; break ; } } if ( ! bFound) ptInt.x = ( dMinX + dMaxX) / 2 ; } else if ( nVec1[1] != nVec2[1]) { ptInt.x = ( nVec1[0] + 0.5) * m_dStep ; ptInt.z = ( nVec1[2] + 0.5) * m_dStep ; int nMinJ = min( nVec1[1], nVec2[1]) ; int nMaxJ = max( nVec1[1], nVec2[1]) ; double dMinY = ( nMinJ + 0.5) * m_dStep ; double dMaxY = ( nMaxJ + 0.5) * m_dStep ; unsigned int nDexel = nVec1[0] * m_nNx[2] + nVec1[2] ; size_t nSize = m_Values[2][nDexel].size() ; bool bFound = false ; for ( size_t j = 0 ; j < nSize ; j += 2) { double dy1 = m_Values[2][nDexel][j].dZVal ; double dy2 = m_Values[2][nDexel][j+1].dZVal ; if ( dy1 < dMinY - EPS_SMALL && dy2 > dMinY - EPS_SMALL && dy2 < dMaxY + EPS_SMALL) { ptInt.y = dy2 ; bFound = true ; break ; } else if ( dy1 > dMinY - EPS_SMALL && dy1 < dMaxY + EPS_SMALL && dy2 > dMaxY + EPS_SMALL) { ptInt.y = dy1 ; bFound = true ; break ; } } if ( ! bFound) ptInt.y = ( dMinY + dMaxY) / 2 ; } else if ( nVec1[2] != nVec2[2]) { ptInt.x = ( nVec1[0] + 0.5) * m_dStep ; ptInt.y = ( nVec1[1] + 0.5) * m_dStep ; int nMinK = min( nVec1[2], nVec2[2]) ; int nMaxK = max( nVec1[2], nVec2[2]) ; double dMinZ = ( nMinK + 0.5) * m_dStep ; double dMaxZ = ( nMaxK + 0.5) * m_dStep ; unsigned int nDexel = nVec1[1] * m_nNx[0] + nVec1[0] ; size_t nSize = m_Values[0][nDexel].size() ; bool bFound = false ; for ( size_t k = 0 ; k < nSize ; k += 2) { double dz1 = m_Values[0][nDexel][k].dZVal ; double dz2 = m_Values[0][nDexel][k+1].dZVal ; if ( dz1 < dMinZ - EPS_SMALL && dz2 > dMinZ - EPS_SMALL && dz2 < dMaxZ + EPS_SMALL) { ptInt.z = dz2 ; bFound = true ; break ; } else if ( dz1 > dMinZ - EPS_SMALL && dz1 < dMaxZ + EPS_SMALL && dz2 > dMaxZ + EPS_SMALL) { ptInt.z = dz1 ; bFound = true ; break ; } } if ( ! bFound) ptInt.z = ( dMinZ + dMaxZ) / 2 ; } return true ; } //---------------------------------------------------------------------------- bool VolZmap::IntersPos( int nVec1[], int nVec2[], Point3d& ptInt, Vector3d& vtNormal) const { double Eps = EPS_SMALL ; if ( nVec1[0] != nVec2[0]) { int nMinI = min( nVec1[0], nVec2[0]) ; int nMaxI = max( nVec1[0], nVec2[0]) ; double dMinX = ( nMinI + 0.5) * m_dStep ; double dMaxX = ( nMaxI + 0.5) * m_dStep ; ptInt.y = ( nVec1[1] + 0.5) * m_dStep ; ptInt.z = ( nVec1[2] + 0.5) * m_dStep ; unsigned int nDexel = nVec1[2] * m_nNx[1] + nVec1[1] ; size_t nSize = m_Values[1][nDexel].size() ; bool bFound = false ; for ( size_t i = 0 ; i < nSize ; i += 2) { double dx1 = m_Values[1][nDexel][i].dZVal ; double dx2 = m_Values[1][nDexel][i+1].dZVal ; if ( dx1 < dMinX - Eps && dx2 > dMinX - Eps && dx2 < dMaxX + Eps) { ptInt.x = dx2 ; vtNormal = m_Values[1][nDexel][i+1].vtN ; bFound = true ; break ; } else if ( dx1 > dMinX - Eps && dx1 < dMaxX + Eps && dx2 > dMaxX + Eps) { ptInt.x = dx1 ; vtNormal = m_Values[1][nDexel][i].vtN ; bFound = true ; break ; }/* if ( dx1 < dMinX + Eps && dx2 > dMinX && dx2 < dMaxX - Eps) { ptInt.x = dx2 ; vtNormal = m_Values[1][nDexel][i+1].vtN ; bFound = true ; break ; } else if ( dx1 > dMinX + Eps && dx1 < dMaxX + Eps && dx2 > dMaxX + Eps) { ptInt.x = dx1 ; vtNormal = m_Values[1][nDexel][i].vtN ; bFound = true ; break ; }*/ } if ( ! bFound) { ptInt.x = ( dMinX + dMaxX) / 2 ; /* if ( IsThereMat( nVec1[0], nVec1[1], nVec1[2])) { double dY = ( nVec2[1] + 0.5) * m_dStep ; double dZ = ( nVec2[2] + 0.5) * m_dStep ; unsigned int YDirMinIndex = 0 ; unsigned int ZDirMinIndex = 0 ; double dYMinDist = DBL_MAX ; double dZMinDist = DBL_MAX ; unsigned int nYDirGrid = nVec2[0] * m_nNx[2] + nVec2[2] ; unsigned int nZDirGrid = nVec2[1] * m_nNx[0] + nVec2[0] ; for ( unsigned int t = 0 ; t < m_Values[0][nZDirGrid].size() ; ++ t) { double dMZ1 = m_Values[0][nZDirGrid][t].dZVal ; if ( abs( dMZ1 - dZ) < dZMinDist) { ZDirMinIndex = t ; dZMinDist = abs( dMZ1 - dZ) ; } } for ( unsigned int t = 0 ; t < m_Values[2][nYDirGrid].size() ; ++ t) { double dMY1 = m_Values[2][nYDirGrid][t].dZVal ; if ( abs( dMY1 - dY) < dYMinDist) { YDirMinIndex = t ; dYMinDist = abs( dMY1 - dY) ; } } unsigned int AbsoluteMinIndex ; if ( dZMinDist < dYMinDist) vtNormal = m_Values[0][nZDirGrid][ZDirMinIndex].vtN ; else vtNormal = m_Values[0][nYDirGrid][YDirMinIndex].vtN ; } else { }*/ } } else if ( nVec1[1] != nVec2[1]) { int nMinJ = min( nVec1[1], nVec2[1]) ; int nMaxJ = max( nVec1[1], nVec2[1]) ; double dMinY = ( nMinJ + 0.5) * m_dStep ; double dMaxY = ( nMaxJ + 0.5) * m_dStep ; ptInt.x = ( nVec1[0] + 0.5) * m_dStep ; ptInt.z = ( nVec1[2] + 0.5) * m_dStep ; unsigned int nDexel = nVec1[0] * m_nNx[2] + nVec1[2] ; size_t nSize = m_Values[2][nDexel].size() ; bool bFound = false ; for ( size_t j = 0 ; j < nSize ; j += 2) { double dy1 = m_Values[2][nDexel][j].dZVal ; double dy2 = m_Values[2][nDexel][j+1].dZVal ; if ( dy1 < dMinY - Eps && dy2 > dMinY - Eps && dy2 < dMaxY + Eps) { ptInt.y = dy2 ; vtNormal = m_Values[2][nDexel][j+1].vtN ; bFound = true ; break ; } else if ( dy1 > dMinY - Eps && dy1 < dMaxY + Eps && dy2 > dMaxY + Eps) { ptInt.y = dy1 ; vtNormal = m_Values[2][nDexel][j].vtN ; bFound = true ; break ; }/* if ( dy1 < dMinY + Eps && dy2 > dMinY && dy2 < dMaxY - Eps) { ptInt.y = dy2 ; vtNormal = m_Values[2][nDexel][j+1].vtN ; bFound = true ; break ; } else if ( dy1 > dMinY + Eps && dy1 < dMaxY + Eps && dy2 > dMaxY + Eps) { ptInt.y = dy1 ; vtNormal = m_Values[2][nDexel][j].vtN ; bFound = true ; break ; }*/ } if ( ! bFound) { ptInt.y = ( dMinY + dMaxY) / 2 ; // Versore Normale ??? } } else if ( nVec1[2] != nVec2[2]) { int nMinK = min( nVec1[2], nVec2[2]) ; int nMaxK = max( nVec1[2], nVec2[2]) ; double dMinZ = ( nMinK + 0.5) * m_dStep ; double dMaxZ = ( nMaxK + 0.5) * m_dStep ; ptInt.x = ( nVec1[0] + 0.5) * m_dStep ; ptInt.y = ( nVec1[1] + 0.5) * m_dStep ; unsigned int nDexel = nVec1[1] * m_nNx[0] + nVec1[0] ; size_t nSize = m_Values[0][nDexel].size() ; bool bFound = false ; for ( size_t k = 0 ; k < nSize ; k += 2) { double dz1 = m_Values[0][nDexel][k].dZVal ; double dz2 = m_Values[0][nDexel][k+1].dZVal ; if ( dz1 < dMinZ - Eps && dz2 > dMinZ - Eps && dz2 < dMaxZ + Eps) { ptInt.z = dz2 ; vtNormal = m_Values[0][nDexel][k+1].vtN ; bFound = true ; break ; } else if ( dz1 > dMinZ - Eps && dz1 < dMaxZ + Eps && dz2 > dMaxZ + Eps) { ptInt.z = dz1 ; vtNormal = m_Values[0][nDexel][k].vtN ; bFound = true ; break ; }/* if ( dz1 < dMinZ + Eps && dz2 > dMinZ && dz2 < dMaxZ - Eps) { ptInt.z = dz2 ; vtNormal = m_Values[0][nDexel][k+1].vtN ; bFound = true ; break ; } else if ( dz1 > dMinZ + Eps && dz1 < dMaxZ + Eps && dz2 > dMaxZ + Eps) { ptInt.z = dz1 ; vtNormal = m_Values[0][nDexel][k].vtN ; bFound = true ; break ; }*/ } if ( ! bFound) { ptInt.z = ( dMinZ + dMaxZ) / 2 ; // Versore Normale ??? } } return true ; } //---------------------------------------------------------------------------- bool VolZmap::NewIntersPos( int nVec1[], int nVec2[], bool bFirstCorner, Point3d& ptInt, Vector3d& vtNormal) const { double dEps = 2 * EPS_SMALL ; if ( nVec1[0] != nVec2[0]) { int nMinI = min( nVec1[0], nVec2[0]) ; int nMaxI = max( nVec1[0], nVec2[0]) ; double dMinX = ( nMinI + 0.5) * m_dStep ; double dMaxX = ( nMaxI + 0.5) * m_dStep ; ptInt.y = ( nVec1[1] + 0.5) * m_dStep ; ptInt.z = ( nVec1[2] + 0.5) * m_dStep ; size_t nDexel = nVec1[2] * m_nNx[1] + nVec1[1] ; size_t nSize = m_Values[1][nDexel].size() ; if ( bFirstCorner) { size_t n = nSize - 1 ; double dX = m_Values[1][nDexel][n].dZVal ; while ( n > 0 && dX > dMinX - dEps) { if ( dX < dMaxX + dEps) { ptInt.x = dX ; vtNormal = m_Values[1][nDexel][n].vtN ; break ; } if ( n == 1) break ; n -= 2 ; dX = m_Values[1][nDexel][n].dZVal ; } } else { size_t n = 0 ; double dX = m_Values[1][nDexel][0].dZVal ; while ( n <= nSize - 2 && dX < dMaxX + dEps) { if ( dX > dMinX - dEps) { ptInt.x = dX ; vtNormal = m_Values[1][nDexel][n].vtN ; break ; } if ( n == nSize - 2) break ; n += 2 ; dX = m_Values[1][nDexel][n].dZVal ; } } } else if ( nVec1[1] != nVec2[1]) { int nMinJ = min( nVec1[1], nVec2[1]) ; int nMaxJ = max( nVec1[1], nVec2[1]) ; double dMinY = ( nMinJ + 0.5) * m_dStep ; double dMaxY = ( nMaxJ + 0.5) * m_dStep ; ptInt.x = ( nVec1[0] + 0.5) * m_dStep ; ptInt.z = ( nVec1[2] + 0.5) * m_dStep ; size_t nDexel = nVec1[0] * m_nNx[2] + nVec1[2] ; size_t nSize = m_Values[2][nDexel].size() ; if ( bFirstCorner) { size_t n = nSize - 1 ; double dY = m_Values[2][nDexel][n].dZVal ; while ( n > 0 && dY > dMinY - dEps) { if ( dY < dMaxY + dEps) { ptInt.y = dY ; vtNormal = m_Values[2][nDexel][n].vtN ; break ; } if ( n == 1) break ; n -= 2 ; dY = m_Values[2][nDexel][n].dZVal ; } } else { size_t n = 0 ; double dY = m_Values[2][nDexel][0].dZVal ; while ( n <= nSize - 2 && dY < dMaxY + dEps) { if ( dY > dMinY - dEps) { ptInt.y = dY ; vtNormal = m_Values[2][nDexel][n].vtN ; break ; } if ( n == nSize - 2) break ; n += 2 ; dY = m_Values[2][nDexel][n].dZVal ; } } } else if ( nVec1[2] != nVec2[2]) { int nMinK = min( nVec1[2], nVec2[2]) ; int nMaxK = max( nVec1[2], nVec2[2]) ; double dMinZ = ( nMinK + 0.5) * m_dStep ; double dMaxZ = ( nMaxK + 0.5) * m_dStep ; ptInt.x = ( nVec1[0] + 0.5) * m_dStep ; ptInt.y = ( nVec1[1] + 0.5) * m_dStep ; size_t nDexel = nVec1[1] * m_nNx[0] + nVec1[0] ; size_t nSize = m_Values[0][nDexel].size() ; if ( bFirstCorner) { size_t n = nSize - 1 ; double dZ = m_Values[0][nDexel][n].dZVal ; while ( n > 0 && dZ > dMinZ - dEps) { if ( dZ < dMaxZ + dEps) { ptInt.z = dZ ; vtNormal = m_Values[0][nDexel][n].vtN ; break ; } if ( n == 1) break ; n -= 2 ; dZ = m_Values[0][nDexel][n].dZVal ; } } else { size_t n = 0 ; double dZ = m_Values[0][nDexel][0].dZVal ; while ( n <= nSize - 2 && dZ < dMaxZ + dEps) { if ( dZ > dMinZ - dEps) { ptInt.z = dZ ; vtNormal = m_Values[0][nDexel][n].vtN ; break ; } if ( n == nSize - 2) break ; n += 2 ; dZ = m_Values[0][nDexel][n].dZVal ; } } } return true ; } //---------------------------------------------------------------------------- bool VolZmap::GetBlockIJK( int nIJK[], int nBlock) const { // Controllo sulla validità del blocco if ( nBlock < 0 || nBlock >= int( m_nNumBlock)) return false ; // Calcolo posizione del blocco nel reticolo nIJK[0] = ( nBlock % int( m_nFracLin[0] * m_nFracLin[1])) % int( m_nFracLin[0]) ; nIJK[1] = ( nBlock % int( m_nFracLin[0] * m_nFracLin[1])) / int( m_nFracLin[0]) ; nIJK[2] = ( nBlock / int( m_nFracLin[0] * m_nFracLin[1])) ; return true ; } //---------------------------------------------------------------------------- bool VolZmap::GetBlockLimitsIJK( const int nIJK[], int nLimits[]) 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 ; // Calcolo limiti per l'indice i nLimits[0] = ( nIJK[0] == 0 ? - 1 : nIJK[0] * int( m_nDexNumPBlock)) ; nLimits[1] = ( nIJK[0] + 1 == int( m_nFracLin[0]) ? int( m_nNx[0]) : ( nIJK[0] + 1) * int( m_nDexNumPBlock)) ; // Calcolo limiti per l'indice j nLimits[2] =( nIJK[1] == 0 ? - 1 : nIJK[1] * int( m_nDexNumPBlock)) ; nLimits[3] = ( nIJK[1] + 1 == int( m_nFracLin[1]) ? int( m_nNy[0]) : ( nIJK[1] + 1) * int( m_nDexNumPBlock)) ; // Calcolo limiti per l'indice k nLimits[4] = ( nIJK[2] == 0 ? - 1 : nIJK[2] * int( m_nDexNumPBlock)) ; nLimits[5] = ( nIJK[2] + 1 == int( m_nFracLin[2]) ? int( m_nNy[1]) : ( nIJK[2] + 1) * int( m_nDexNumPBlock)) ; return true ; }