//---------------------------------------------------------------------------- // 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} ; //, Edge2 = 3} ; //---------------------------------------------------------------------------- //bool //TestOnNormal( const VectorField CompoVert[], int nCompoElem, int& FeatureType, Vector3d& vtFeature) //{ // int nI, nJ ; // double dMinCosTheta = 1.001 ; // const double dCosThetaSharp = sqrt( 3) / 2 ; 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;//0.7;//0.5 ; // // for ( int i = 0 ; i < nCompoElem ; ++ i) { // // double dAbsCurrentCos = abs( CompoVert[i].vtNorm * vtK) ; // // if ( dAbsCurrentCos > dCosPhiCorner) { // // FeatureType = Corner ; // vtFeature = vtK ; // return true ; // } // } // // FeatureType = Edge ; // vtFeature = vtK ; // // return true ; //} //---------------------------------------------------------------------------- bool TestOnNormal( const VectorField CompoVert[], int nCompoElem, int& FeatureType, Vector3d& vtFeature, int& nMin1, int& nMin2) { int nI, nJ ; double dMinCosTheta = 1.001 ; const double dCosThetaSharp = 0.8; 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 ; double dDotTest1 = CompoVert[nI].vtNorm * vtK ; double dDotTest2 = CompoVert[nJ].vtNorm * vtK ; nMin1 = nI ; nMin2 = nJ ; const double dCosPhiCorner = 0.5 ; double dMaxAbsCos = 0 ; for ( int i = 0 ; i < nCompoElem ; ++ i) { double dAbsCurrentCos = abs( CompoVert[i].vtNorm * vtK) ; if ( dAbsCurrentCos > dMaxAbsCos) dMaxAbsCos = dAbsCurrentCos ; } if ( dMaxAbsCos > dCosPhiCorner) FeatureType = Corner ; else FeatureType = Edge ; vtFeature = vtK ; 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 MarchingCubes( lstTria) ; return true ; } //---------------------------------------------------------------------------- bool VolZmap::GetTriangles( bool bAllBlocks, INTVECTOR& nModifiedBlocks, TRIA3DLISTVECTOR& vLstTria) const { // Caso di singola mappa if ( m_nMapNum == 1) { const int MAX_DIM_CHUNK = 128 ; nModifiedBlocks.resize( m_nNumBlock) ; vLstTria.reserve( m_nNumBlock) ; // Ciclo sui blocchi for ( size_t t = 0 ; t < m_nNumBlock ; ++ t) { // Se il blocco deve essere aggiornato, eseguo if ( bAllBlocks || m_BlockToUpdate[t]) { // preparo lista vLstTria.emplace_back() ; nModifiedBlocks[t] = int( vLstTria.size()) - 1 ; // Calcolo posizione del blocco nella griglia int nIBlock = int( t) % int( m_nFracLin[0]) ; int nJBlock = int( t) / 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, vLstTria.back()) ; } } m_BlockToUpdate[t] = false ; } else nModifiedBlocks[t] = -1 ; } } // Caso con tre mappe else { nModifiedBlocks.resize( m_nNumBlock + 1) ; vLstTria.reserve( m_nNumBlock + 1) ; TriaMatrix VecTriHold ; VecTriHold.resize( m_nNumBlock) ; // Ciclo sui blocchi for ( size_t t = 0 ; t < m_nNumBlock ; ++ t) { // Calcolo i limiti sugli indici dei voxel del blocco int nIJK[3] ; // Vettore indici i,j,k del blocco GetBlockIJKFromN( int( t), nIJK) ; // Vettore limiti sugli indici dei voxel del blocco int LimitIndexes[6] ; GetBlockLimitsIJK( nIJK, LimitIndexes) ; // Se il blocco deve essere processato, eseguo if ( bAllBlocks || m_BlockToUpdate[t]) { vLstTria.emplace_back() ; nModifiedBlocks[t] = int( vLstTria.size()) - 1 ; m_InterBlockTria[t].clear() ; ExtMarchingCubes( LimitIndexes, t, vLstTria.back(), VecTriHold[t]) ; // Flipping fra voxel interni FlipEdgesII( VecTriHold[t]) ; //for ( int nDeltaI = -1 ; nDeltaI <= 1 ; ++ nDeltaI) { // for ( int nDeltaJ = -1 ; nDeltaJ <= 1 ; ++ nDeltaJ) { // for ( int nDeltaK = -1 ; nDeltaK <= 1 ; ++ nDeltaK) { // // // Flipping fra voxel dello stesso blocco // if ( nDeltaI == 0 && nDeltaJ == 0 && nDeltaK == 0) { // // // Flipping fra voxel interni // FlipEdgesII( VecTriHold[t]) ; // // Flipping fra voxel interni e di frontiera // //FlipEdgesIB( VecTriHold[t], t) ; // } // // Flipping fra voxel adiacenti // else { // int nAdjN ; // // if ( GetAdjBlockToBlock( int( t), nDeltaI, nDeltaJ, nDeltaK, nAdjN)) { // // // Qui bisogna fare il flipping fra i blocchi: t indica il blocco // // corrente appena aggiornato. Guardiamo il suo vettore di voxel bislacchi // // e cicliamo su tali voxel. Per ogni voxel bislacco guardiamo i voxel adiacenti. // // Se un voxel adiacente è nel blocco adiacente corrente ricalcoliamo i triangoli e flippiamo. // size_t nBoundaryVoxelNum = m_InterBlockTria[t].size() ; // for ( size_t a = 0 ; a < nBoundaryVoxelNum ; ++ a) { // // // vedere se il voxel appartiene al blocco adiacente // int nBVoxI = m_InterBlockTria[t][a].i ; // int nBVoxJ = m_InterBlockTria[t][a].j ; // int nBVoxK = m_InterBlockTria[t][a].k ; // // Indici blocco adiacente // int nAdjBlockIJK[3] ; // GetBlockIJKFromN( nAdjN, nAdjBlockIJK) ; // // Limiti indici di voxel per il blocco adiacente // int nLimitAdjBlock[6] ; // GetBlockLimitsIJK( nAdjBlockIJK, nLimitAdjBlock) ; // // // Se il blocco adiacente esiste eseguo il flipping // if ( nBVoxI >= nLimitAdjBlock[0] && nBVoxI < nLimitAdjBlock[1] && // nBVoxJ >= nLimitAdjBlock[2] && nBVoxJ < nLimitAdjBlock[3] && // nBVoxK >= nLimitAdjBlock[4] && nBVoxK < nLimitAdjBlock[5]) { // // m_InterBlockTria[nAdjN].clear() ; // } // } // } // } // } // } //} m_BlockToUpdate[t] = false ; } else nModifiedBlocks[t] = -1 ; } // Copio i triangoli di frontiera in una matrice gemella // di m_InterBlockTria per avere sempre a disposizione // i triangoli non flippati. TriaMatrix InterBlockTria ; InterBlockTria = m_InterBlockTria ; FlipEdgesBB( InterBlockTria) ; // ciclo sui blocchi for ( size_t t = 0 ; t < m_nNumBlock ; ++ t) { // ciclo sui voxel del blocco 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) 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]) ; } } } } } nModifiedBlocks[nModifiedBlocks.size()-1] = int( vLstTria.size() - 1) ; } return true ; } //---------------------------------------------------------------------------- int VolZmap::GetBlockCount( void) const { return m_nNumBlock + 1 ; } //---------------------------------------------------------------------------- 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[], size_t nBlockNumber, 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) { if ( !(i == 31 && j == -1 && k == 49)) ;// continue ; // Riconoscimento dei voxel di frontiera bool bBoundary = false ; int nVoxIndexes[3] = { i, j, k} ; if ( IsAVoxelOnBoundary( nLimits, nVoxIndexes, true)) bBoundary = true ; // Indici i,j,k dei vertici int IndexCorner[8][3] = { { i, j, k}, { 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] ; // Flag sulla regolatrità dei campi scalare e vettoriale: // se i campi sono regolari esso resta vero, altrimenti // assume il valore falso. bool bReg = true ; // 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, // se i campi scalare e vettoriale non sono regolari bReg diviene falso. if ( ! IntersPos( IndexCorner[n1], IndexCorner[n2], bN1, VecField[EdgeIndex].ptInt, VecField[EdgeIndex].vtNorm)) bReg = false ; // Riporto punti e normali nel sistema locale in cui // è immerso lo Zmap col suo sistema di riferimento. VecField[EdgeIndex].ptInt.ToGlob( m_MapFrame[0]) ; VecField[EdgeIndex].vtNorm.ToGlob( m_MapFrame[0]) ; } // Determino il numero di componenti connesse nel voxel 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 nInnerFeatureInVoxel = 0 ; int nBorderFeatureInVoxel = 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]] ; // Controllo per il caso piano su una griglia // con versori normali a due a due paralleli. bool bGridControl = true ; if ( nVertComp == 4) { int nIndArrey[4] ; // Riempio un vettore con gli indici dei punti for ( int nCpInd = 0 ; nCpInd < nVertComp ; ++ nCpInd) nIndArrey[nCpInd] = TriangleTableEn[nIndex][1][nCpInd + nTableOffset + 1] ; // Ordino i 4 indici in senso crescente for ( int nSrtInd1 = 0 ; nSrtInd1 < nVertComp - 1 ; ++ nSrtInd1) { for ( int nSrtInd2 = nSrtInd1 + 1 ; nSrtInd2 < nVertComp ; ++ nSrtInd2) { if ( nIndArrey[nSrtInd1] > nIndArrey[nSrtInd2]) swap( nIndArrey[nSrtInd1], nIndArrey[nSrtInd2]) ; } } if ( ( nIndArrey[0] == 0 && nIndArrey[1] == 2 && nIndArrey[2] == 9 && nIndArrey[3] == 10) || ( nIndArrey[0] == 4 && nIndArrey[1] == 6 && nIndArrey[2] == 9 && nIndArrey[3] == 10) || ( nIndArrey[0] == 4 && nIndArrey[1] == 6 && nIndArrey[2] == 8 && nIndArrey[3] == 11) || ( nIndArrey[0] == 0 && nIndArrey[1] == 2 && nIndArrey[2] == 8 && nIndArrey[3] == 11) || ( nIndArrey[0] == 1 && nIndArrey[1] == 3 && nIndArrey[2] == 8 && nIndArrey[3] == 9 ) || ( nIndArrey[0] == 1 && nIndArrey[1] == 3 && nIndArrey[2] == 10 && nIndArrey[3] == 11) || ( nIndArrey[0] == 5 && nIndArrey[1] == 7 && nIndArrey[2] == 10 && nIndArrey[3] == 11) || ( nIndArrey[0] == 5 && nIndArrey[1] == 7 && nIndArrey[2] == 8 && nIndArrey[3] == 9 )) { if ( AreSameVectorApprox( VecField[nIndArrey[0]].vtNorm, VecField[nIndArrey[1]].vtNorm) && AreSameVectorApprox( VecField[nIndArrey[2]].vtNorm, VecField[nIndArrey[3]].vtNorm) && abs( VecField[nIndArrey[0]].vtNorm * VecField[nIndArrey[2]].vtNorm) < EPS_SMALL) { Point3d ptBarycenter = ( CompoVert[0].ptInt + CompoVert[1].ptInt + CompoVert[2].ptInt + CompoVert[3].ptInt) / 4 ; if ( abs( CompoVert[0].ptInt.x - ptBarycenter.x) < EPS_SMALL && abs( CompoVert[1].ptInt.x - ptBarycenter.x) < EPS_SMALL && abs( CompoVert[2].ptInt.x - ptBarycenter.x) < EPS_SMALL && abs( CompoVert[3].ptInt.x - ptBarycenter.x) < EPS_SMALL) { double dXBar = ptBarycenter.x / m_dStep - 0.5 ; int nBarLimSup = int( m_nNx[0]) ; int nBarInd = 0 ; while ( nBarInd < nBarLimSup) { double dXInd = double( nBarInd) ; if ( abs( dXInd - dXBar) < EPS_SMALL) { bGridControl = false ; break ; } ++ nBarInd ; } } else if ( abs( CompoVert[0].ptInt.y - ptBarycenter.y) < EPS_SMALL && abs( CompoVert[1].ptInt.y - ptBarycenter.y) < EPS_SMALL && abs( CompoVert[2].ptInt.y - ptBarycenter.y) < EPS_SMALL && abs( CompoVert[3].ptInt.y - ptBarycenter.y) < EPS_SMALL) { double dYBar = ptBarycenter.y / m_dStep - 0.5 ; int nBarLimSup = int( m_nNy[0]) ; int nBarInd = 0 ; while ( nBarInd < nBarLimSup) { double dYInd = double( nBarInd) ; if ( abs( dYInd - dYBar) < EPS_SMALL) { bGridControl = false ; break ; } ++ nBarInd ; } } else if ( abs( CompoVert[0].ptInt.z - ptBarycenter.z) < EPS_SMALL && abs( CompoVert[1].ptInt.z - ptBarycenter.z) < EPS_SMALL && abs( CompoVert[2].ptInt.z - ptBarycenter.z) < EPS_SMALL && abs( CompoVert[3].ptInt.z - ptBarycenter.z) < EPS_SMALL) { double dZBar = ptBarycenter.z / m_dStep - 0.5 ; int nBarLimSup = int( m_nNy[1]) ; int nBarInd = 0 ; while ( nBarInd < nBarLimSup) { double dZInd = double( nBarInd) ; if ( abs( dZInd - dZBar) < EPS_SMALL) { bGridControl = false ; break ; } ++ nBarInd ; } } } } else if ( ( nIndArrey[0] == 0 && nIndArrey[1] == 1 && nIndArrey[2] == 4 && nIndArrey[3] == 5) || ( nIndArrey[0] == 1 && nIndArrey[1] == 2 && nIndArrey[2] == 5 && nIndArrey[3] == 6) || ( nIndArrey[0] == 2 && nIndArrey[1] == 3 && nIndArrey[2] == 6 && nIndArrey[3] == 7) || ( nIndArrey[0] == 0 && nIndArrey[1] == 3 && nIndArrey[2] == 4 && nIndArrey[3] == 7)) { if ( AreSameVectorApprox( VecField[nIndArrey[0]].vtNorm, VecField[nIndArrey[2]].vtNorm) && AreSameVectorApprox( VecField[nIndArrey[1]].vtNorm, VecField[nIndArrey[3]].vtNorm) && abs( VecField[nIndArrey[0]].vtNorm * VecField[nIndArrey[1]].vtNorm) < EPS_SMALL) { Point3d ptBarycenter = ( CompoVert[0].ptInt + CompoVert[1].ptInt + CompoVert[2].ptInt + CompoVert[3].ptInt) / 4 ; if ( abs( CompoVert[0].ptInt.x - ptBarycenter.x) < EPS_SMALL && abs( CompoVert[1].ptInt.x - ptBarycenter.x) < EPS_SMALL && abs( CompoVert[2].ptInt.x - ptBarycenter.x) < EPS_SMALL && abs( CompoVert[3].ptInt.x - ptBarycenter.x) < EPS_SMALL) { double dXBar = ptBarycenter.x / m_dStep - 0.5 ; int nBarLimSup = int( m_nNx[0]) ; int nBarInd = 0 ; while ( nBarInd < nBarLimSup) { double dXInd = double( nBarInd) ; if ( abs( dXInd - dXBar) < EPS_SMALL) { bGridControl = false ; break ; } ++ nBarInd ; } } else if ( abs( CompoVert[0].ptInt.y - ptBarycenter.y) < EPS_SMALL && abs( CompoVert[1].ptInt.y - ptBarycenter.y) < EPS_SMALL && abs( CompoVert[2].ptInt.y - ptBarycenter.y) < EPS_SMALL && abs( CompoVert[3].ptInt.y - ptBarycenter.y) < EPS_SMALL) { double dYBar = ptBarycenter.y / m_dStep - 0.5 ; int nBarLimSup = int( m_nNy[0]) ; int nBarInd = 0 ; while ( nBarInd < nBarLimSup) { double dYInd = double( nBarInd) ; if ( abs( dYInd - dYBar) < EPS_SMALL) { bGridControl = false ; break ; } ++ nBarInd ; } } } } } int nFeatureType ; // Valuto le relazioni reciproche fra le normali e // se i punti sono su un piano di griglia. Vector3d vtFeatureAxis ; bool bNormal = false ; // Indici dei punti ove le normali // formano il massimo angolo int nMin1, nMin2 ; // Se i campi sono regolari valuto le normali // per stabilire se eseguire ExtMC o MC, // altrimenti eseguo MC. if ( bReg) bNormal = TestOnNormal( CompoVert, nVertComp, nFeatureType, vtFeatureAxis, nMin1, nMin2) ; // Flag ExtMC bool bExtMC = bNormal && bGridControl ; // Extended MC if ( bExtMC) { // Passo al sistema di riferimento del baricentro Point3d ptGravityCenter( 0, 0, 0) ; for ( int ni = 0 ; ni < nVertComp ; ++ ni) ptGravityCenter = ptGravityCenter + CompoVert[ni].ptInt ; ptGravityCenter = 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 superficie piana 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() ; // Se la feature è un edge annullo // il valore singolare minore. if ( nFeatureType == Edge) dSingularValue( 2) = 0 ; // Back substitution: risolvo il sistema USV*x = b // Calcolo U^T b double vTemp[3] ; for ( int ni = 0 ; ni < 3 ; ++ ni) { double s = 0 ; if ( dSingularValue( ni) > 0) { 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 ; // Flag sulla distanza del vertice dal // baricentro del sistema di punti bool bOutside = false ; if ( dDistFeature > dMaxDist) bOutside = true ; // 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 ; // Flag sull'inversione delle normali bool bInvNormB = false ; // Questo controllo viene eseguito solo se // il vertice è distante dal baricentro // entro la soglia stabilita. if ( ! bOutside) { for ( int ni = 0 ; ni < nVertComp ; ++ ni) { int nj = ( ni + 1 < nVertComp) ? ni + 1 : 0 ; // Il triangolo è pronto CurrentTriangle.Set( ptSol, CompoVert[nj].ptInt, CompoVert[ni].ptInt) ; CurrentTriangle.Validate( true) ; // Aggiungo triangolo al vettore temporaneo 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)) { ptSol = ptGravityCenter ; bInvNormB = true ; triContainer.resize( 0) ; ni = -1 ; } } } // Questo flag esprime se il vertice è, entro la tolleranza, // interno o esterno al voxel a cui appartiene. bool bInsideVoxel = IsPointInsideVoxelApprox( i, j, k, ptSol) ; // Proprietà gemoetriche locali // del campo vettoriale: se vero // procediamo con un ulteriore // analisi per eliminare triangoli // invertiti. bool bLocProp = false ; if ( nVertComp == 4) { int nNegDotNum = 0 ; for ( int nLocIndI = 0 ; nLocIndI < 3 ; ++ nLocIndI) { for ( int nLocIndJ = nLocIndI + 1 ; nLocIndJ < 4 ; ++ nLocIndJ) { if ( CompoVert[nLocIndI].vtNorm * CompoVert[nLocIndJ].vtNorm < 0) { nNegDotNum ++ ; } } } if ( nNegDotNum == 3) bLocProp = true ; } // Se è necessario si effettua l'ulteriore analisi // Questo controllo si effettua se la feature non // esce dai limiti, non esce dal voxel e il primo // controllo sull'inversione delle normali ha dato // esito negativo. if ( nFeatureType == 2 && bLocProp && ! ( bOutside || bInsideVoxel || bInvNormB)) { if ( abs( nMin1 - nMin2) == 1 || abs( nMin1 - nMin2) == 3) { int nSum = nMin1 + nMin2 ; int nTriNum = ( nSum == 3 && ( nMin1 == 3 || nMin2 == 3) ? 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)) { 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 ; double dLenVV = vtVV.Len() ; ptSolZmFrame.ToLoc( m_MapFrame[0]) ; if ( ! IsPointInsideVoxel( i, j, k, ptSol)) { Vector3d vtVS = ptSol - CompoVert[nNm].ptInt ; vtVV.Normalize() ; Vector3d vtVSNew = ( vtVS * vtVV) * vtVV ; ptSol = CompoVert[nNm].ptInt + vtVSNew ; double dLNm = ( ptSol - CompoVert[nNm].ptInt).Len() ; double dLNp = ( ptSol - CompoVert[nNp].ptInt).Len() ; if ( dLNm > dLenVV || dLNp > dLenVV) { ptSol = dLNm < dLNp ? CompoVert[nNm].ptInt : CompoVert[nNp].ptInt ; } triContainer.resize( 0) ; for ( int nc = 0 ; nc < nVertComp ; ++ nc) { int nd = ( nc + 1 < nVertComp) ? nc + 1 : 0 ; // Il triangolo è pronto CurrentTriangle.Set( ptSol, CompoVert[nd].ptInt, CompoVert[nc].ptInt) ; CurrentTriangle.Validate( true) ; // Aggiungo triangolo al vettore temporaneo triContainer.emplace_back( CurrentTriangle) ; } } } } else { int nPrev1 = nMin1 == 0 ? 3 : nMin1 - 1 ; int nPrev2 = nMin2 == 0 ? 3 : nMin2 - 1 ; int nNext1 = nMin1 == 3 ? 0 : nMin1 + 1 ; int nNext2 = nMin2 == 3 ? 0 : nMin2 + 1 ; int nNeighbourIndex ; int nStartIndex ; if ( CompoVert[nPrev1].vtNorm * CompoVert[nMin1].vtNorm > 0) { bool bTestOnVert = abs( nPrev1 - nMin2) == 0 || abs( nPrev1 - nMin2) == 3 ; nNeighbourIndex = bTestOnVert ? nPrev1 : nNext1 ; nStartIndex = nMin2 ; } else { bool bTestOnVert = abs( nPrev2 - nMin1) == 0 || abs( nPrev2 - nMin1) == 3 ; nNeighbourIndex = bTestOnVert ? nPrev2 : nNext2 ; nStartIndex = nMin1 ; } Vector3d vtVV = CompoVert[nNeighbourIndex].ptInt - CompoVert[nStartIndex].ptInt ; double dVVLen = vtVV.Len() ; vtVV.Normalize() ; Vector3d vtVF = ptSol - CompoVert[nStartIndex].ptInt ; Vector3d vtNewVF = ( vtVF * vtVV) * vtVV ; double dNewVFLen = vtNewVF.Len() ; if ( dNewVFLen > dVVLen) { ptSol = CompoVert[nNeighbourIndex].ptInt ; } else { 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 CurrentTriangle.Set( ptSol, CompoVert[nd].ptInt, CompoVert[nc].ptInt) ; CurrentTriangle.Validate( true) ; // Aggiungo triangolo al vettore temporaneo triContainer.emplace_back( CurrentTriangle) ; } } } // Valuto normali: questo è ancora un controllo // sulle normali, se risultano in tutti i punti // approssimativamente uguali passiamo alla // routine standard int nContSize = int( 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) break ; } // Se la feature non è fuori dai limiti // e i triangoli formano una superficie // non piana confermo ExtMC if ( ! ( bOutside || bPlane)) { // Aggiorno il numero di feature. // ++ nInnerFeatureInVoxel ; //// Se siamo in presenza della prima feature del //// voxel, ridimensiono il vettore che contiene //// la struttura dati del voxel. // if ( nInnerFeatureInVoxel == 1) { // triHold.resize( triHold.size() + 1) ; // // Questi dati dipendono solo dal voxel, // // quindi sono validi per tutte le // // componenti che vi appartengono. // int nCurrent = int( triHold.size()) - 1 ; // triHold[nCurrent].i = i ; // triHold[nCurrent].j = j ; // triHold[nCurrent].k = k ; // } //// Indice che identifica la posizione del voxel //// nel vector. // int nCurrent = int( triHold.size()) - 1 ; //// Aggiungo vertice della componente //// connessa alla lista dei vertici. // triHold[nCurrent].ptCompoVert.emplace_back( ptSol) ; // int nOldFeatureNum = int( triHold[nCurrent].vCompoTria.size()) ; // int nNewFeatureNum = nOldFeatureNum + 1 ; //// Aggiungo una componente per il vettore //// dei triangoli della componente connessa. // triHold[nCurrent].vCompoTria.resize( nNewFeatureNum) ; TRIA3DVECTOR vInnerTriaTemp, vBorderTriaTemp ; for ( int ni = 0 ; ni < nContSize ; ++ ni) { // Se l'area è non nulla aggiungo il triangolo // a uno dei due vettori. if ( triContainer[ni].GetArea() > EPS_SMALL) { int nVoxIJK[3] = { i, j, k} ; Point3d ptTemp = ptSol ; ptTemp.ToLoc( m_MapFrame[0]) ; Triangle3d trTemp = triContainer[ni] ; trTemp.ToLoc( m_MapFrame[0]) ; if ( ! IsATriangleOnBorder( trTemp, ptTemp, nLimits, nVoxIJK, bBoundary)) vInnerTriaTemp.emplace_back( triContainer[ni]) ; else vBorderTriaTemp.emplace_back( triContainer[ni]) ; } } // Metto i triangoli negli appositi contenitori: // Triangoli interni if ( vInnerTriaTemp.size() > 0) { // Aggiorno il numero di feature. ++ nInnerFeatureInVoxel ; // Se siamo in presenza della prima feature del // voxel, ridimensiono il vettore che contiene // la struttura dati del voxel. if ( nInnerFeatureInVoxel == 1) { triHold.resize( triHold.size() + 1) ; // Questi dati dipendono solo dal voxel, // quindi sono validi per tutte le // componenti che vi appartengono. int nCurrent = int( triHold.size()) - 1 ; triHold[nCurrent].i = i ; triHold[nCurrent].j = j ; triHold[nCurrent].k = k ; } // Indice che identifica la posizione del voxel // nel vector. int nCurrent = int( triHold.size()) - 1 ; // Aggiungo vertice della componente // connessa alla lista dei vertici. triHold[nCurrent].ptCompoVert.emplace_back( ptSol) ; int nOldFeatureNum = int( triHold[nCurrent].vCompoTria.size()) ; int nNewFeatureNum = nOldFeatureNum + 1 ; // Aggiungo una componente per il vettore // dei triangoli della componente connessa. triHold[nCurrent].vCompoTria.resize( nNewFeatureNum) ; for ( int ni = 0 ; ni < int( vInnerTriaTemp.size()) ; ++ ni) triHold[nCurrent].vCompoTria[nOldFeatureNum].emplace_back( vInnerTriaTemp[ni]) ; } // Triangoli di frontiera if ( vBorderTriaTemp.size() > 0) { // Aggiorno il numero di feature. ++ nBorderFeatureInVoxel ; // Se siamo in presenza della prima feature del // voxel, ridimensiono il vettore che contiene // la struttura dati del voxel. if ( nBorderFeatureInVoxel == 1) { m_InterBlockTria[nBlockNumber].resize( m_InterBlockTria[nBlockNumber].size() + 1) ; // Questi dati dipendono solo dal voxel, // quindi sono validi per tutte le // componenti che vi appartengono. int nCurrent = int( m_InterBlockTria[nBlockNumber].size()) - 1 ; m_InterBlockTria[nBlockNumber][nCurrent].i = i ; m_InterBlockTria[nBlockNumber][nCurrent].j = j ; m_InterBlockTria[nBlockNumber][nCurrent].k = k ; } // Indice che identifica la posizione del voxel // nel vector. int nCurrent = int( m_InterBlockTria[nBlockNumber].size()) - 1 ; // Aggiungo vertice della componente // connessa alla lista dei vertici. m_InterBlockTria[nBlockNumber][nCurrent].ptCompoVert.emplace_back( ptSol) ; int nOldFeatureNum = int( m_InterBlockTria[nBlockNumber][nCurrent].vCompoTria.size()) ; int nNewFeatureNum = nOldFeatureNum + 1 ; // Aggiungo una componente per il vettore // dei triangoli della componente connessa. m_InterBlockTria[nBlockNumber][nCurrent].vCompoTria.resize( nNewFeatureNum) ; for ( int ni = 0 ; ni < int( vBorderTriaTemp.size()) ; ++ ni) m_InterBlockTria[nBlockNumber][nCurrent].vCompoTria[nOldFeatureNum].emplace_back( vBorderTriaTemp[ni]) ; } } // ExtMC non confermato, si passa a MC else { // 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 ( CurrentTriangle.GetArea() > EPS_SMALL) 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 ( CurrentTriangle.GetArea() > EPS_SMALL) lstTria.emplace_back( CurrentTriangle) ; } } nTableOffset += nVertComp ; } } } } return true ; } //---------------------------------------------------------------------------- //bool //VolZmap::FlipEdges( std::vector& VecTriHold) const //{ // double dSqEps = EPS_SMALL * EPS_SMALL ; // // // Ciclo sui blocchi // for ( size_t t1 = 0 ; t1 < m_nNumBlock ; ++ t1) { // // // Determino i,j,k del primo blocco // int nIJK1[3] ; // GetBlockIJK( int( t1), nIJK1) ; // // for ( size_t t2 = t1 ; t2 < m_nNumBlock ; ++ t2) { // // // Determino i,j,k del secondo blocco // int nIJK2[3] ; // GetBlockIJK( int( t2), nIJK2) ; // // // Se i blocchi sono adiacenti o coincidenti proseguo // if ( ( nIJK2[0] >= nIJK1[0] - 1 && nIJK2[0] <= nIJK1[0] + 1) || // ( nIJK2[1] >= nIJK1[1] - 1 && nIJK2[1] <= nIJK1[1] + 1) || // ( nIJK2[2] >= nIJK1[2] - 1 && nIJK2[2] <= nIJK1[2] + 1)) { // // // Numero di voxel in cui si presentano sharp feature // int nVoxelNum1 = int( VecTriHold[t1].size()) ; // int nVoxelNum2 = int( VecTriHold[t2].size()) ; // // // Determino estremi del ciclo sui voxel esterno // int nSt1 = 0 ; // int nEn1 = nVoxelNum1 - ( t1 == t2 ? 1 : 0) ; // // // Ciclo su tali voxel // for ( int n1 = nSt1 ; n1 < nEn1 ; ++ n1) { // // // Determino estremi del ciclo sui voxel interno // int nSt2 = ( t1 == t2 ? nSt1 : 0) ; // int nEn2 = nVoxelNum2 ; // // for ( int n2 = n1 ; n2 < nVoxelNum2 ; ++ n2) { // // // Se i voxel sono adiacenti proseguo // if ( VecTriHold[t2][n2].i == VecTriHold[t1][n1].i - 1 || // VecTriHold[t2][n2].i == VecTriHold[t1][n1].i + 1 || // VecTriHold[t2][n2].j == VecTriHold[t1][n1].j - 1 || // VecTriHold[t2][n2].j == VecTriHold[t1][n1].j + 1 || // VecTriHold[t2][n2].k == VecTriHold[t1][n1].k - 1 || // VecTriHold[t2][n2].k == VecTriHold[t1][n1].k + 1) { // // // Numero delle componenti connesse nei due voxel // int nNumCompo1 = int( VecTriHold[t1][n1].ptCompoVert.size()) ; // int nNumCompo2 = int( VecTriHold[t2][n2].ptCompoVert.size()) ; // // int nCompo1 = 0 ; // // // Ciclo sulle componenti // for ( ; nCompo1 < nNumCompo1 ; ++ nCompo1) { // // int nCompo2 = ( t1 == t2 && n1 == n2 ? nCompo1 + 1 : 0) ; // // for ( ; nCompo2 < nNumCompo2 ; ++ nCompo2) { // // // Numero di triangoli per le componenti connesse // int nTriNum1 = int( VecTriHold[t1][n1].vCompoTria[nCompo1].size()) ; // int nTriNum2 = int( VecTriHold[t2][n2].vCompoTria[nCompo2].size()) ; // // for ( int nTri1 = 0 ; nTri1 < nTriNum1 ; ++ nTri1) { // // bool bModified = false ; // // for ( int nTri2 = 0 ; nTri2 < nTriNum2 ; ++ nTri2) { // // INTVECTOR SharedIndex ; // // for ( int nVert1 = 0 ; nVert1 < 3 ; ++ nVert1) { // // for ( int nVert2 = 0 ; nVert2 < 3 ; ++ nVert2) { // // Point3d ptP1 = VecTriHold[t1][n1].vCompoTria[nCompo1][nTri1].GetP( nVert1) ; // Point3d ptP2 = VecTriHold[t2][n2].vCompoTria[nCompo2][nTri2].GetP( nVert2) ; // // if ( SqDist( ptP1, ptP2) < dSqEps) { // // Point3d ptVert1 = VecTriHold[t1][n1].ptCompoVert[nCompo1] ; // Point3d ptVert2 = VecTriHold[t2][n2].ptCompoVert[nCompo2] ; // // if ( ! ( AreSamePointApprox( ptP1, ptVert1) || // AreSamePointApprox( ptP2, ptVert2))) { // // SharedIndex.emplace_back( nVert1) ; // SharedIndex.emplace_back( nVert2) ; // } // } // // if ( SharedIndex.size() > 2) // break ; // } // // if ( SharedIndex.size() > 2) // break ; // } // // // Si deve operare la modifica dei triangoli // if ( SharedIndex.size() > 2) { // // int nProd = ( SharedIndex[2] - SharedIndex[0]) * ( SharedIndex[3] - SharedIndex[1]) ; // // // --- // if ( nProd != 1 && nProd != - 2 && nProd != 4) { // // VecTriHold[t1][n1].vCompoTria[nCompo1][nTri1].SetP( SharedIndex[0], // VecTriHold[t2][n2].ptCompoVert[nCompo2]) ; // VecTriHold[t2][n2].vCompoTria[nCompo2][nTri2].SetP( SharedIndex[3], // VecTriHold[t1][n1].ptCompoVert[nCompo1]) ; // // VecTriHold[t1][n1].vCompoTria[nCompo1][nTri1].Validate( true) ; // VecTriHold[t2][n2].vCompoTria[nCompo2][nTri2].Validate( true) ; // // bModified = true ; // break ; // } // } // } // // if ( bModified) // break ; // } // } // } // } // } // } // } // } // } // // return true ; //} //---------------------------------------------------------------------------- bool VolZmap::FlipEdges( TriaMatrix& VecTriHold) const { double dSqEps = EPS_SMALL * EPS_SMALL ; // Ciclo sui blocchi for ( size_t t1 = 0 ; t1 < m_nNumBlock ; ++ t1) { // Determino i,j,k del primo blocco int nIJK1[3] ; GetBlockIJKFromN( int( t1), nIJK1) ; for ( size_t t2 = t1 ; t2 < m_nNumBlock ; ++ t2) { // Determino i,j,k del secondo blocco int nIJK2[3] ; GetBlockIJKFromN( int( t2), nIJK2) ; // Se i blocchi sono adiacenti o coincidenti proseguo if ( ( nIJK2[0] >= nIJK1[0] - 1 && nIJK2[0] <= nIJK1[0] + 1) || ( nIJK2[1] >= nIJK1[1] - 1 && nIJK2[1] <= nIJK1[1] + 1) || ( nIJK2[2] >= nIJK1[2] - 1 && nIJK2[2] <= nIJK1[2] + 1)) { // Numero di voxel in cui si presentano sharp feature int nVoxelNum1 = int( VecTriHold[t1].size()) ; int nVoxelNum2 = int( VecTriHold[t2].size()) ; // Determino estremi del ciclo sui voxel esterno int nSt1 = 0 ; int nEn1 = nVoxelNum1 - ( t1 == t2 ? 1 : 0) ; // Ciclo su tali voxel for ( int n1 = nSt1 ; n1 < nEn1 ; ++ n1) { // Determino estremi del ciclo sui voxel interno int nSt2 = ( t1 == t2 ? n1 + 1 : 0) ; int nEn2 = nVoxelNum2 ; for ( int n2 = nSt2 ; n2 < nVoxelNum2 ; ++ n2) { // Se i voxel sono adiacenti proseguo if ( VecTriHold[t2][n2].i == VecTriHold[t1][n1].i - 1 || VecTriHold[t2][n2].i == VecTriHold[t1][n1].i + 1 || VecTriHold[t2][n2].j == VecTriHold[t1][n1].j - 1 || VecTriHold[t2][n2].j == VecTriHold[t1][n1].j + 1 || VecTriHold[t2][n2].k == VecTriHold[t1][n1].k - 1 || VecTriHold[t2][n2].k == VecTriHold[t1][n1].k + 1) { // Numero delle componenti connesse nei due voxel int nNumCompo1 = int( VecTriHold[t1][n1].ptCompoVert.size()) ; int nNumCompo2 = int( VecTriHold[t2][n2].ptCompoVert.size()) ; int nCompo1 = 0 ; // Ciclo sulle componenti for ( ; nCompo1 < nNumCompo1 ; ++ nCompo1) { int nCompo2 = ( t1 == t2 && n1 == n2 ? nCompo1 + 1 : 0) ; for ( ; nCompo2 < nNumCompo2 ; ++ nCompo2) { // Numero di triangoli per le componenti connesse int nTriNum1 = int( VecTriHold[t1][n1].vCompoTria[nCompo1].size()) ; int nTriNum2 = int( VecTriHold[t2][n2].vCompoTria[nCompo2].size()) ; for ( int nTri1 = 0 ; nTri1 < nTriNum1 ; ++ nTri1) { bool bModified = false ; for ( int nTri2 = 0 ; nTri2 < nTriNum2 ; ++ nTri2) { INTVECTOR SharedIndex ; for ( int nVert1 = 0 ; nVert1 < 3 ; ++ nVert1) { for ( int nVert2 = 0 ; nVert2 < 3 ; ++ nVert2) { Point3d ptP1 = VecTriHold[t1][n1].vCompoTria[nCompo1][nTri1].GetP( nVert1) ; Point3d ptP2 = VecTriHold[t2][n2].vCompoTria[nCompo2][nTri2].GetP( nVert2) ; if ( SqDist( ptP1, ptP2) < dSqEps) { Point3d ptVert1 = VecTriHold[t1][n1].ptCompoVert[nCompo1] ; Point3d ptVert2 = VecTriHold[t2][n2].ptCompoVert[nCompo2] ; if ( ! ( AreSamePointApprox( ptP1, ptVert1) || AreSamePointApprox( ptP2, ptVert2))) { SharedIndex.emplace_back( nVert1) ; SharedIndex.emplace_back( nVert2) ; } } if ( SharedIndex.size() > 2) break ; } if ( SharedIndex.size() > 2) break ; } // Si deve operare la modifica dei triangoli if ( SharedIndex.size() > 2) { int nProd = ( SharedIndex[2] - SharedIndex[0]) * ( SharedIndex[3] - SharedIndex[1]) ; // --- if ( nProd != 1 && nProd != - 2 && nProd != 4) { VecTriHold[t1][n1].vCompoTria[nCompo1][nTri1].SetP( SharedIndex[0], VecTriHold[t2][n2].ptCompoVert[nCompo2]) ; VecTriHold[t2][n2].vCompoTria[nCompo2][nTri2].SetP( SharedIndex[3], VecTriHold[t1][n1].ptCompoVert[nCompo1]) ; VecTriHold[t1][n1].vCompoTria[nCompo1][nTri1].Validate( true) ; VecTriHold[t2][n2].vCompoTria[nCompo2][nTri2].Validate( true) ; if ( t1 != t2) { int bau = 1 ; } bModified = true ; break ; } } } if ( bModified) break ; } } } } } } } } } return true ; } //---------------------------------------------------------------------------- //bool //VolZmap::FlipEdges( std::vector& VecTriHold) const //{ // double dSqEps = EPS_SMALL * EPS_SMALL ; // size_t nEnd = VecTriHold.size() ; //m_nNumBlock ; // // Ciclo sui blocchi // for ( size_t t = 0 ; t < nEnd ; ++ t) { // // // Numero di voxel in cui si presentano sharp feature // int nVoxelNum = int( VecTriHold[t].size()) ; // // // Determino estremi del ciclo sui voxel esterno // int nSt1 = 0 ; // int nEn1 = nVoxelNum - 1 ; // // // Ciclo su tali voxel // for ( int n1 = nSt1 ; n1 < nEn1 ; ++ n1) { // // // Determino estremi del ciclo sui voxel interno // int nSt2 = nSt1 + 1 ; // int nEn2 = nVoxelNum ; // // for ( int n2 = nSt2 ; n2 < nEn2 ; ++ n2) { // // // Se i voxel sono adiacenti proseguo // if ( VecTriHold[t][n2].i == VecTriHold[t][n1].i + 1 || // VecTriHold[t][n2].j == VecTriHold[t][n1].j + 1 || // VecTriHold[t][n2].k == VecTriHold[t][n1].k + 1) { // // // Numero delle componenti connesse nei due voxel // int nNumCompo1 = int( VecTriHold[t][n1].ptCompoVert.size()) ; // int nNumCompo2 = int( VecTriHold[t][n2].ptCompoVert.size()) ; // // int nCompo1 = 0 ; // // // Ciclo sulle componenti // for ( ; nCompo1 < nNumCompo1 ; ++ nCompo1) { // // int nCompo2 = 0 ; // // for ( ; nCompo2 < nNumCompo2 ; ++ nCompo2) { // // // Numero di triangoli per le componenti connesse // int nTriNum1 = int( VecTriHold[t][n1].vCompoTria[nCompo1].size()) ; // int nTriNum2 = int( VecTriHold[t][n2].vCompoTria[nCompo2].size()) ; // // for ( int nTri1 = 0 ; nTri1 < nTriNum1 ; ++ nTri1) { // // bool bModified = false ; // // for ( int nTri2 = 0 ; nTri2 < nTriNum2 ; ++ nTri2) { // // INTVECTOR SharedIndex ; // // for ( int nVert1 = 0 ; nVert1 < 3 ; ++ nVert1) { // // for ( int nVert2 = 0 ; nVert2 < 3 ; ++ nVert2) { // // Point3d ptP1 = VecTriHold[t][n1].vCompoTria[nCompo1][nTri1].GetP( nVert1) ; // Point3d ptP2 = VecTriHold[t][n2].vCompoTria[nCompo2][nTri2].GetP( nVert2) ; // // if ( SqDist( ptP1, ptP2) < dSqEps) { // // Point3d ptVert1 = VecTriHold[t][n1].ptCompoVert[nCompo1] ; // Point3d ptVert2 = VecTriHold[t][n2].ptCompoVert[nCompo2] ; // // if ( ! ( AreSamePointApprox( ptP1, ptVert1) || // AreSamePointApprox( ptP2, ptVert2))) { // // SharedIndex.emplace_back( nVert1) ; // SharedIndex.emplace_back( nVert2) ; // } // } // // if ( SharedIndex.size() > 2) // break ; // } // // if ( SharedIndex.size() > 2) // break ; // } // // // Si deve operare la modifica dei triangoli // if ( SharedIndex.size() > 2) { // // int nProd = ( SharedIndex[2] - SharedIndex[0]) * ( SharedIndex[3] - SharedIndex[1]) ; // // // --- // if ( nProd != 1 && nProd != - 2 && nProd != 4) { // // Triangle3d TriVox1 = VecTriHold[t][n1].vCompoTria[nCompo1][nTri1] ; // Triangle3d TriVox2 = VecTriHold[t][n2].vCompoTria[nCompo2][nTri2] ; // // TriVox1.SetP( SharedIndex[0], VecTriHold[t][n2].ptCompoVert[nCompo2]) ; // TriVox2.SetP( SharedIndex[3], VecTriHold[t][n1].ptCompoVert[nCompo1]) ; // // TriVox1.Validate( true) ; // TriVox2.Validate( true) ; // // questo dentro il commento è un'assurdità // /*bool bContinue = true ; // // for ( int nControl = 0 ; nControl < nTriNum1 ; ++ nControl) { // // Vector3d vtControl = VecTriHold[t][n1].vCompoTria[nCompo1][nControl].GetN() ; // Vector3d vtVox = TriVox1.GetN() ; // // if ( vtControl * vtVox < - 0.9) { // bContinue = false ; // break ; // } // } // // if ( bContinue) { // // for ( int nControl = 0 ; nControl < nTriNum2 ; ++ nControl) { // // Vector3d vtControl = VecTriHold[t][n2].vCompoTria[nCompo2][nControl].GetN() ; // Vector3d vtVox = TriVox2.GetN() ; // // if ( vtControl * vtVox < - 0.9) { // bContinue = false ; // break ; // } // } // } // // if ( bContinue) {*/ // // VecTriHold[t][n1].vCompoTria[nCompo1][nTri1].SetP( SharedIndex[0], // VecTriHold[t][n2].ptCompoVert[nCompo2]) ; // VecTriHold[t][n2].vCompoTria[nCompo2][nTri2].SetP( SharedIndex[3], // VecTriHold[t][n1].ptCompoVert[nCompo1]) ; // // VecTriHold[t][n1].vCompoTria[nCompo1][nTri1].Validate( true) ; // VecTriHold[t][n2].vCompoTria[nCompo2][nTri2].Validate( true) ; // // bModified = true ; // break ; // /*}*/ // } // } // } // // if ( bModified) // break ; // } // } // } // } // } // } // } // return true ; //} //---------------------------------------------------------------------------- bool VolZmap::FlipEdgesII( TriHolder& TriHold) const { double dSqEps = EPS_SMALL * EPS_SMALL ; // Numero di voxel in cui si presentano sharp feature int nVoxelNum = int( TriHold.size()) ; // Determino estremi del ciclo sui voxel esterno int nSt1 = 0 ; int nEn1 = nVoxelNum - 1 ; // Ciclo su tali voxel for ( int n1 = nSt1 ; n1 < nEn1 ; ++ n1) { // Determino estremi del ciclo sui voxel interno int nSt2 = nSt1 + 1 ; int nEn2 = nVoxelNum ; for ( int n2 = nSt2 ; n2 < nEn2 ; ++ n2) { // Se i voxel sono adiacenti proseguo if ( TriHold[n2].i == TriHold[n1].i + 1 || TriHold[n2].j == TriHold[n1].j + 1 || TriHold[n2].k == TriHold[n1].k + 1) { // Numero delle componenti connesse nei due voxel int nNumCompo1 = int( TriHold[n1].ptCompoVert.size()) ; int nNumCompo2 = int( TriHold[n2].ptCompoVert.size()) ; int nCompo1 = 0 ; // Ciclo sulle componenti for ( ; nCompo1 < nNumCompo1 ; ++ nCompo1) { int nCompo2 = 0 ; for ( ; nCompo2 < nNumCompo2 ; ++ nCompo2) { // Numero di triangoli per le componenti connesse int nTriNum1 = int( TriHold[n1].vCompoTria[nCompo1].size()) ; int nTriNum2 = int( TriHold[n2].vCompoTria[nCompo2].size()) ; for ( int nTri1 = 0 ; nTri1 < nTriNum1 ; ++ nTri1) { bool bModified = false ; for ( int nTri2 = 0 ; nTri2 < nTriNum2 ; ++ nTri2) { INTVECTOR SharedIndex ; for ( int nVert1 = 0 ; nVert1 < 3 ; ++ nVert1) { for ( int nVert2 = 0 ; nVert2 < 3 ; ++ nVert2) { Point3d ptP1 = TriHold[n1].vCompoTria[nCompo1][nTri1].GetP( nVert1) ; Point3d ptP2 = TriHold[n2].vCompoTria[nCompo2][nTri2].GetP( nVert2) ; if ( SqDist( ptP1, ptP2) < dSqEps) { Point3d ptVert1 = TriHold[n1].ptCompoVert[nCompo1] ; Point3d ptVert2 = TriHold[n2].ptCompoVert[nCompo2] ; if ( ! ( AreSamePointApprox( ptP1, ptVert1) || AreSamePointApprox( ptP2, ptVert2))) { SharedIndex.emplace_back( nVert1) ; SharedIndex.emplace_back( nVert2) ; } } if ( SharedIndex.size() > 2) break ; } if ( SharedIndex.size() > 2) break ; } // Si deve operare la modifica dei triangoli if ( SharedIndex.size() > 2) { int nProd = ( SharedIndex[2] - SharedIndex[0]) * ( SharedIndex[3] - SharedIndex[1]) ; // --- if ( nProd != 1 && nProd != - 2 && nProd != 4) { Triangle3d TriVox1 = TriHold[n1].vCompoTria[nCompo1][nTri1] ; Triangle3d TriVox2 = TriHold[n2].vCompoTria[nCompo2][nTri2] ; TriVox1.SetP( SharedIndex[0], TriHold[n2].ptCompoVert[nCompo2]) ; TriVox2.SetP( SharedIndex[3], TriHold[n1].ptCompoVert[nCompo1]) ; TriVox1.Validate( true) ; TriVox2.Validate( true) ; // questo dentro il commento è un'assurdità /*bool bContinue = true ; for ( int nControl = 0 ; nControl < nTriNum1 ; ++ nControl) { Vector3d vtControl = VecTriHold[t][n1].vCompoTria[nCompo1][nControl].GetN() ; Vector3d vtVox = TriVox1.GetN() ; if ( vtControl * vtVox < - 0.9) { bContinue = false ; break ; } } if ( bContinue) { for ( int nControl = 0 ; nControl < nTriNum2 ; ++ nControl) { Vector3d vtControl = VecTriHold[t][n2].vCompoTria[nCompo2][nControl].GetN() ; Vector3d vtVox = TriVox2.GetN() ; if ( vtControl * vtVox < - 0.9) { bContinue = false ; break ; } } } if ( bContinue) {*/ TriHold[n1].vCompoTria[nCompo1][nTri1].SetP( SharedIndex[0], TriHold[n2].ptCompoVert[nCompo2]) ; TriHold[n2].vCompoTria[nCompo2][nTri2].SetP( SharedIndex[3], TriHold[n1].ptCompoVert[nCompo1]) ; TriHold[n1].vCompoTria[nCompo1][nTri1].Validate( true) ; TriHold[n2].vCompoTria[nCompo2][nTri2].Validate( true) ; bModified = true ; break ; /*}*/ } } } if ( bModified) break ; } } } } } } return true ; } //---------------------------------------------------------------------------- bool VolZmap::FlipEdgesBB( TriaMatrix& InterTria) const { double dSqEps = EPS_SMALL * EPS_SMALL ; // Numero di blocchi size_t nBlocksNum = InterTria.size() ; // ciclo sui blocchi for ( size_t tFB = 0 ; tFB < nBlocksNum - 1 ; ++ tFB) { for ( size_t tLB = tFB + 1 ; tLB < nBlocksNum ; ++ tLB) { int nFBijk[3] ; int nLBijk[3] ; GetBlockIJKFromN( int( tFB), nFBijk) ; GetBlockIJKFromN( int( tLB), nLBijk) ; // Se i blocchi non sono adiacenti salto l'iterazione if ( ! ( abs( nFBijk[0] - nLBijk[0]) <= 1 && abs( nFBijk[1] - nLBijk[1]) <= 1 && abs( nFBijk[2] - nLBijk[2]) <= 1)) continue ; // Numero di voxel nei blocchi correnti size_t nVoxelNumFB = InterTria[tFB].size() ; size_t nVoxelNumLB = InterTria[tLB].size() ; // Ciclo sui voxel dei due blocchi for ( size_t tVFB = 0 ; tVFB < nVoxelNumFB ; ++ tVFB) { for ( size_t tVLB = 0 ; tVLB < nVoxelNumLB ; ++ tVLB) { // Se i voxel non sono adiacenti salto l'iterazione if ( ! ( abs( InterTria[tFB][tVFB].i - InterTria[tLB][tVLB].i) <= 1 && abs( InterTria[tFB][tVFB].j - InterTria[tLB][tVLB].j) <= 1 && abs( InterTria[tFB][tVFB].k - InterTria[tLB][tVLB].k) <= 1)) continue ; // Numero di componenti connesse dei voxel size_t nCompoVFBNum = InterTria[tFB][tVFB].ptCompoVert.size() ; size_t nCompoVLBNum = InterTria[tLB][tVLB].ptCompoVert.size() ; // Ciclo sulle componenti connesse for ( size_t tCmpF = 0 ; tCmpF < nCompoVFBNum ; ++ tCmpF) { for ( size_t tCmpL = 0 ; tCmpL < nCompoVLBNum ; ++ tCmpL) { // Numero di triangoli delle componenti connesse size_t nTriFBNum = InterTria[tFB][tVFB].vCompoTria[tCmpF].size() ; size_t nTriLBNum = InterTria[tLB][tVLB].vCompoTria[tCmpL].size() ; for ( size_t tTriFB = 0 ; tTriFB < nTriFBNum ; ++ tTriFB) { bool bModified = false ; for ( size_t tTriLB = 0 ; tTriLB < nTriLBNum ; ++ tTriLB) { INTVECTOR SharedIndex ; for ( int nVertF = 0 ; nVertF < 3 ; ++ nVertF) { for ( int nVertL = 0 ; nVertL < 3 ; ++ nVertL) { Point3d ptPF = InterTria[tFB][tVFB].vCompoTria[tCmpF][tTriFB].GetP( nVertF) ; Point3d ptPL = InterTria[tLB][tVLB].vCompoTria[tCmpL][tTriLB].GetP( nVertL) ; if ( SqDist( ptPF, ptPL) < dSqEps) { Point3d ptVertF = InterTria[tFB][tVFB].ptCompoVert[tCmpF] ; Point3d ptVertL = InterTria[tLB][tVLB].ptCompoVert[tCmpL] ; if ( ! ( AreSamePointApprox( ptPF, ptVertF) || AreSamePointApprox( ptPL, ptVertL))) { SharedIndex.emplace_back( nVertF) ; SharedIndex.emplace_back( nVertL) ; } } if ( SharedIndex.size() > 2) break ; } if ( SharedIndex.size() > 2) break ; } // Si deve operare la modifica dei triangoli if ( SharedIndex.size() > 2) { int nProd = ( SharedIndex[2] - SharedIndex[0]) * ( SharedIndex[3] - SharedIndex[1]) ; if ( nProd != 1 && nProd != - 2 && nProd != 4) { Triangle3d TriVoxF = InterTria[tFB][tVFB].vCompoTria[tCmpF][tTriFB] ; Triangle3d TriVoxL = InterTria[tLB][tVLB].vCompoTria[tCmpL][tTriLB] ; TriVoxF.SetP( SharedIndex[0], InterTria[tLB][tVLB].ptCompoVert[tCmpL]) ; TriVoxL.SetP( SharedIndex[3], InterTria[tFB][tVFB].ptCompoVert[tCmpF]) ; TriVoxF.Validate( true) ; TriVoxL.Validate( true) ; InterTria[tFB][tVFB].vCompoTria[tCmpF][tTriFB].SetP( SharedIndex[0], InterTria[tLB][tVLB].ptCompoVert[tCmpL]) ; InterTria[tLB][tVLB].vCompoTria[tCmpL][tTriLB].SetP( SharedIndex[3], InterTria[tFB][tVFB].ptCompoVert[tCmpF]) ; InterTria[tFB][tVFB].vCompoTria[tCmpF][tTriFB].Validate( true) ; InterTria[tLB][tVLB].vCompoTria[tCmpL][tTriLB].Validate( true) ; bModified = true ; break ; } } } if ( bModified) break ; } } } } } } } return true ; } //---------------------------------------------------------------------------- bool VolZmap::FlipEdgesIB( TriHolder& InnerVox, size_t nBlock) const { double dSqEps = EPS_SMALL * EPS_SMALL ; size_t nNumI = InnerVox.size() ; size_t nNumB = m_InterBlockTria[nBlock].size() ; for ( size_t tB = 0 ; tB < nNumB ; ++ tB) { for ( size_t tI = 0 ; tI < nNumI ; ++ tI) { if ( abs( InnerVox[tI].i - m_InterBlockTria[nBlock][tB].i) <= 1 && abs( InnerVox[tI].j - m_InterBlockTria[nBlock][tB].j) <= 1 && abs( InnerVox[tI].k - m_InterBlockTria[nBlock][tB].k) <= 1) { // Numero delle componenti connessVecTre nei due voxel int nNumCompoB = int( m_InterBlockTria[nBlock][tB].ptCompoVert.size()) ; int nNumCompoI = int( InnerVox[tI].ptCompoVert.size()) ; int nCompoB = 0 ; // Ciclo sulle componenti for ( ; nCompoB < nNumCompoB ; ++ nCompoB) { int nCompoI = 0 ; for ( ; nCompoI < nNumCompoI ; ++ nCompoI) { // Numero di triangoli per le componenti connesse int nTriNumB = int( m_InterBlockTria[nBlock][tB].vCompoTria[nCompoB].size()) ; int nTriNumI = int( InnerVox[tI].vCompoTria[nCompoI].size()) ; for ( int nTriB = 0 ; nTriB < nTriNumB ; ++ nTriB) { bool bModified = false ; for ( int nTriI = 0 ; nTriI < nTriNumI ; ++ nTriI) { INTVECTOR SharedIndex ; for ( int nVertB = 0 ; nVertB < 3 ; ++ nVertB) { for ( int nVertI = 0 ; nVertI < 3 ; ++ nVertI) { Point3d ptPB = m_InterBlockTria[nBlock][tB].vCompoTria[nCompoB][nTriB].GetP( nVertB) ; Point3d ptPI = InnerVox[tI].vCompoTria[nCompoI][nTriI].GetP( nVertI) ; if ( SqDist( ptPB, ptPI) < dSqEps) { Point3d ptVertB = m_InterBlockTria[nBlock][tB].ptCompoVert[nCompoB] ; Point3d ptVertI = InnerVox[tI].ptCompoVert[nCompoI] ; if ( ! ( AreSamePointApprox( ptPB, ptVertB) || AreSamePointApprox( ptPI, ptVertI))) { SharedIndex.emplace_back( nVertB) ; SharedIndex.emplace_back( nVertI) ; } } if ( SharedIndex.size() > 2) break ; } if ( SharedIndex.size() > 2) break ; } // Si deve operare la modifica dei triangoli if ( SharedIndex.size() > 2) { int nProd = ( SharedIndex[2] - SharedIndex[0]) * ( SharedIndex[3] - SharedIndex[1]) ; // --- if ( nProd != 1 && nProd != - 2 && nProd != 4) { Triangle3d TriVoxB = m_InterBlockTria[nBlock][tB].vCompoTria[nCompoB][nTriB] ; Triangle3d TriVoxI = InnerVox[tI].vCompoTria[nCompoI][nTriI] ; TriVoxB.SetP( SharedIndex[0], InnerVox[tI].ptCompoVert[nCompoI]) ; TriVoxI.SetP( SharedIndex[3], m_InterBlockTria[nBlock][tB].ptCompoVert[nCompoB]) ; TriVoxB.Validate( true) ; TriVoxI.Validate( true) ; m_InterBlockTria[nBlock][tB].vCompoTria[nCompoB][nTriB].SetP( SharedIndex[0], InnerVox[tI].ptCompoVert[nCompoI]) ; InnerVox[tI].vCompoTria[nCompoI][nTriI].SetP( SharedIndex[3], m_InterBlockTria[nBlock][tB].ptCompoVert[nCompoB]) ; m_InterBlockTria[nBlock][tB].vCompoTria[nCompoB][nTriB].Validate( true) ; InnerVox[tI].vCompoTria[nCompoI][nTriI].Validate( true) ; bModified = true ; break ; } } } if ( bModified) break ; } } } } } } return true ; } //---------------------------------------------------------------------------- bool VolZmap::FlipEdgesBB( size_t nCurBlock, size_t nAdjBlock) const { double dSqEps = EPS_SMALL * EPS_SMALL ; return true ; } //---------------------------------------------------------------------------- bool VolZmap::FlipEdgesIB( TriHolder& InnerVox, TriHolder& BoundaryVox) const { double dSqEps = EPS_SMALL * EPS_SMALL ; size_t nNumI = InnerVox.size() ; size_t nNumB = BoundaryVox.size() ; for ( size_t tB = 0 ; tB < nNumB ; ++ tB) { for ( size_t tI = 0 ; tI < nNumI ; ++ tI) { if ( abs( InnerVox[tI].i - BoundaryVox[tB].i) <= 1 && abs( InnerVox[tI].j - BoundaryVox[tB].j) <= 1 && abs( InnerVox[tI].k - BoundaryVox[tB].k) <= 1) { // Numero delle componenti connessVecTre nei due voxel int nNumCompoB = int( BoundaryVox[tB].ptCompoVert.size()) ; int nNumCompoI = int( InnerVox[tI].ptCompoVert.size()) ; int nCompoB = 0 ; // Ciclo sulle componenti for ( ; nCompoB < nNumCompoB ; ++ nCompoB) { int nCompoI = 0 ; for ( ; nCompoI < nNumCompoI ; ++ nCompoI) { // Numero di triangoli per le componenti connesse int nTriNumB = int( BoundaryVox[tB].vCompoTria[nCompoB].size()) ; int nTriNumI = int( InnerVox[tI].vCompoTria[nCompoI].size()) ; for ( int nTriB = 0 ; nTriB < nTriNumB ; ++ nTriB) { bool bModified = false ; for ( int nTriI = 0 ; nTriI < nTriNumI ; ++ nTriI) { INTVECTOR SharedIndex ; for ( int nVertB = 0 ; nVertB < 3 ; ++ nVertB) { for ( int nVertI = 0 ; nVertI < 3 ; ++ nVertI) { Point3d ptPB = BoundaryVox[tB].vCompoTria[nCompoB][nTriB].GetP( nVertB) ; Point3d ptPI = InnerVox[tI].vCompoTria[nCompoI][nTriI].GetP( nVertI) ; if ( SqDist( ptPB, ptPI) < dSqEps) { Point3d ptVertB = BoundaryVox[tB].ptCompoVert[nCompoB] ; Point3d ptVertI = InnerVox[tI].ptCompoVert[nCompoI] ; if ( ! ( AreSamePointApprox( ptPB, ptVertB) || AreSamePointApprox( ptPI, ptVertI))) { SharedIndex.emplace_back( nVertB) ; SharedIndex.emplace_back( nVertI) ; } } if ( SharedIndex.size() > 2) break ; } if ( SharedIndex.size() > 2) break ; } // Si deve operare la modifica dei triangoli if ( SharedIndex.size() > 2) { int nProd = ( SharedIndex[2] - SharedIndex[0]) * ( SharedIndex[3] - SharedIndex[1]) ; // --- if ( nProd != 1 && nProd != - 2 && nProd != 4) { Triangle3d TriVoxB = BoundaryVox[tB].vCompoTria[nCompoB][nTriB] ; Triangle3d TriVoxI = InnerVox[tI].vCompoTria[nCompoI][nTriI] ; TriVoxB.SetP( SharedIndex[0], InnerVox[tI].ptCompoVert[nCompoI]) ; TriVoxI.SetP( SharedIndex[3], BoundaryVox[tB].ptCompoVert[nCompoB]) ; TriVoxB.Validate( true) ; TriVoxI.Validate( true) ; BoundaryVox[tB].vCompoTria[nCompoB][nTriB].SetP( SharedIndex[0], InnerVox[tI].ptCompoVert[nCompoI]) ; InnerVox[tI].vCompoTria[nCompoI][nTriI].SetP( SharedIndex[3], BoundaryVox[tB].ptCompoVert[nCompoB]) ; BoundaryVox[tB].vCompoTria[nCompoB][nTriB].Validate( true) ; InnerVox[tI].vCompoTria[nCompoI][nTriI].Validate( true) ; bModified = true ; break ; } } } if ( bModified) break ; } } } } } } return true ; } //---------------------------------------------------------------------------- bool VolZmap::IsThereMat( int nI, int nJ, int nK) const { 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[], bool bFirstCorner, Point3d& ptInt, Vector3d& vtNormal) const { double dEps = 2 * EPS_SMALL ; bool bFound = false ; 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 ; bFound = true ; 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 ; bFound = true ; break ; } if ( n == nSize - 2) break ; n += 2 ; dX = m_Values[1][nDexel][n].dZVal ; } } if ( ! bFound) ptInt.x = 0.5 * ( dMinX + dMaxX) ; } 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 ; bFound = true ; 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 ; bFound = true ; break ; } if ( n == nSize - 2) break ; n += 2 ; dY = m_Values[2][nDexel][n].dZVal ; } } if ( ! bFound) ptInt.y = 0.5 * ( dMinY + dMaxY) ; } 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 ; bFound = true ; 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 ; bFound = true ; break ; } if ( n == nSize - 2) break ; n += 2 ; dZ = m_Values[0][nDexel][n].dZVal ; } } if ( ! bFound) ptInt.z = 0.5 * ( dMinZ + dMaxZ) ; } return bFound ; } //---------------------------------------------------------------------------- bool VolZmap::GetBlockIJKFromN( int nBlock, int nIJK[]) 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::GetBlockNFromIJK( int nIJK[], int& nBlock) const { // Controllo sulla validità degli indici i, j, k del blocco if ( nIJK[0] < 0 || nIJK[0] >= int( m_nFracLin[0]) || nIJK[1] < 0 || nIJK[1] >= int( m_nFracLin[1]) || nIJK[2] < 0 || nIJK[2] >= int( m_nFracLin[2])) return false ; // Determino il blocco nBlock = m_nFracLin[0] * m_nFracLin[1] * nIJK[2] + m_nFracLin[0] * nIJK[1] + nIJK[0] ; return true ; } //---------------------------------------------------------------------------- bool VolZmap::GetBlockLimitsIJK( const int nIJK[], int nLimits[]) const { // 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 ; } //---------------------------------------------------------------------------- bool VolZmap::IsPointInsideVoxel( int nI, int nJ, int nK, const Point3d& ptP) const { // Controllo sull'ammissibilità del voxel if ( nI < - 1 || nI >= int( m_nNx[0]) || nJ < - 1 || nJ >= int( m_nNy[0]) || nK < - 1 || nK >= int( m_nNy[1])) return false ; int nPointI = int( floor( ( ptP.x - 0.5 * m_dStep) / m_dStep)) ; int nPointJ = int( floor( ( ptP.y - 0.5 * m_dStep) / m_dStep)) ; int nPointK = int( floor( ( ptP.z - 0.5 * m_dStep) / m_dStep)) ; return ( nPointI == nI && nPointJ == nJ && nPointK == nK) ; } //---------------------------------------------------------------------------- bool VolZmap::IsPointInsideVoxelApprox( int nI, int nJ, int nK, const Point3d& ptP) const { // Controllo sull'ammissibilità del voxel if ( nI < - 1 || nI >= int( m_nNx[0]) || nJ < - 1 || nJ >= int( m_nNy[0]) || nK < - 1 || nK >= int( m_nNy[1])) return false ; Point3d ptPZmap = ptP ; ptPZmap.ToLoc( m_MapFrame[0]) ; bool bI = ptPZmap.x > ( nI + 0.5) * m_dStep - EPS_SMALL && ptPZmap.x < ( nI + 1.5) * m_dStep + EPS_SMALL ; bool bJ = ptPZmap.y > ( nJ + 0.5) * m_dStep - EPS_SMALL && ptPZmap.y < ( nJ + 1.5) * m_dStep + EPS_SMALL ; bool bK = ptPZmap.z > ( nK + 0.5) * m_dStep - EPS_SMALL && ptPZmap.z < ( nK + 1.5) * m_dStep + EPS_SMALL ; return ( bI && bJ && bK) ; } //---------------------------------------------------------------------------- bool VolZmap::GetPointVoxel( Point3d& ptP, int& nVoxI, int& nVoxJ, int& nVoxK) const { nVoxI = int( floor( ( ptP.x - 0.5 * m_dStep) / m_dStep)) ; nVoxJ = int( floor( ( ptP.y - 0.5 * m_dStep) / m_dStep)) ; nVoxK = int( floor( ( ptP.z - 0.5 * m_dStep) / m_dStep)) ; return ( nVoxI >= -1 && nVoxI < int( m_nNx[0])) && ( nVoxJ >= -1 && nVoxJ < int( m_nNy[0])) && ( nVoxK >= -1 && nVoxK < int( m_nNy[1])) ; } //---------------------------------------------------------------------------- bool VolZmap::GetVoxelBlockIJK( const int nVoxIJK[], int nBlockIJK[]) const { // Controllo sull'ammissibilità del voxel if ( nVoxIJK[0] <= -2 || nVoxIJK[0] >= int( m_nNx[0]) || nVoxIJK[1] <= -2 || nVoxIJK[1] >= int( m_nNy[0]) || nVoxIJK[2] <= -2 || nVoxIJK[2] >= int( m_nNy[1])) return false ; // Divisioni intere int nIntRatio0 = nVoxIJK[0] / m_nDexNumPBlock ; int nIntRatio1 = nVoxIJK[1] / m_nDexNumPBlock ; int nIntRatio2 = nVoxIJK[2] / m_nDexNumPBlock ; // Resti delle divisioni intere /*int nMod0 = nVoxIJK[0] % m_nDexNumPBlock ; int nMod1 = nVoxIJK[1] % m_nDexNumPBlock ; int nMod2 = nVoxIJK[2] % m_nDexNumPBlock ;*/ // Calcolo indici del blocco nBlockIJK[0] = ( nVoxIJK[0] == -1 ? 0 : max( 0, nIntRatio0 - ( nIntRatio0 == m_nFracLin[0] ? 1 : 0))) ; nBlockIJK[1] = ( nVoxIJK[1] == -1 ? 0 : max( 0, nIntRatio1 - ( nIntRatio1 == m_nFracLin[1] ? 1 : 0))) ; nBlockIJK[2] = ( nVoxIJK[2] == -1 ? 0 : max( 0, nIntRatio2 - ( nIntRatio2 == m_nFracLin[2] ? 1 : 0))) ; return true ; } //---------------------------------------------------------------------------- bool VolZmap::GetAdjBlockToBlock( int nBlockN, int nDeltaI, int nDeltaJ, int nDeltaK, int& nAdjBlockN) const { // Test sulla validità degli incrementi su i,j,k if ( nDeltaI < - 1 || nDeltaI > 1 || nDeltaJ < - 1 || nDeltaJ > 1 || nDeltaK < - 1 || nDeltaK > 1) return false ; // Determino blocco adiacente nAdjBlockN = nBlockN ; nAdjBlockN += nDeltaI ; nAdjBlockN += nDeltaJ * m_nFracLin[0] ; nAdjBlockN += nDeltaK * m_nFracLin[0] * m_nFracLin[1] ; // Se il blocco adiacente esiste restituisco vero, altrimenti falso. return ( nAdjBlockN > -1 && nAdjBlockN < int( m_nNumBlock)) ; } //---------------------------------------------------------------------------- bool VolZmap::IsAVoxelOnBoundary( const int nLimits[], const int nIJK[], bool bType) const { // Test sulla validità dei limiti if ( nLimits[0] < -1 || nLimits[0] > int( m_nNx[0]) || nLimits[1] < -1 || nLimits[1] > int( m_nNx[0]) || nLimits[2] < -1 || nLimits[2] > int( m_nNy[0]) || nLimits[3] < -1 || nLimits[3] > int( m_nNy[0]) || nLimits[4] < -1 || nLimits[4] > int( m_nNy[1]) || nLimits[5] < -1 || nLimits[5] > int( m_nNy[1])) return false ; // Controllo sull'ammissibilità del voxel if ( nIJK[0] <= -2 || nIJK[0] >= int( m_nNx[0]) || nIJK[1] <= -2 || nIJK[1] >= int( m_nNx[1]) || nIJK[2] <= -2 || nIJK[2] >= int( m_nNy[1])) return false ; // Se bType è vero cerchiamo i voxel che // confinano con quelli di atri blocchi, if ( bType) { int nCurrentBlockIJK[3] ; GetVoxelBlockIJK( nIJK, nCurrentBlockIJK) ; // Condizione necessaria è che il voxel sia sulla frontiera if ( IsAVoxelOnBoundary( nLimits, nIJK, false)) { // Ciclo sulle posizioni possibili dei voxel adiacenti for ( int nDeltaI = -1 ; nDeltaI <= 1 ; ++ nDeltaI) { for ( int nDeltaJ = -1 ; nDeltaJ <= 1 ; ++ nDeltaJ) { for ( int nDeltaK = -1 ; nDeltaK <= 1 ; ++ nDeltaK) { if ( nDeltaI == 0 && nDeltaJ == 0 && nDeltaK == 0) continue ; // Indici dei voxel adiacenti int nAdjIJK[3] = { nIJK[0] + nDeltaI, nIJK[1] + nDeltaJ, nIJK[2] + nDeltaK} ; int nAdjBlockIJK[3] ; // Determino (se esiste) il blocco in cui cade il voxel adiacente. if ( GetVoxelBlockIJK( nAdjIJK, nAdjBlockIJK)) { if ( nAdjBlockIJK[0] > -1 && nAdjBlockIJK[0] < int( m_nFracLin[0]) && nAdjBlockIJK[1] > -1 && nAdjBlockIJK[1] < int( m_nFracLin[1]) && nAdjBlockIJK[2] > -1 && nAdjBlockIJK[2] < int( m_nFracLin[2]) && ( nCurrentBlockIJK[0] != nAdjBlockIJK[0] || nCurrentBlockIJK[1] != nAdjBlockIJK[1] || nCurrentBlockIJK[2] != nAdjBlockIJK[2])) { return true ; } } } } } } } // altrimenti cerchiamo i voxel che // sono sulla frontiera del blocco. else { if ( nIJK[0] == nLimits[0] || nIJK[0] == nLimits[1] -1 || nIJK[1] == nLimits[2] || nIJK[1] == nLimits[3] -1 || nIJK[2] == nLimits[4] || nIJK[2] == nLimits[5] -1) return true ; } return false ; } //---------------------------------------------------------------------------- bool VolZmap::IsATriangleOnBorder( const Triangle3d& trTria, const Point3d& ptVert, const int nBlockLimits[], const int nVoxIJK[]) const { // Se il triangolo ha area nulla non proseguiamo if ( trTria.GetArea()) return false ; // Se il voxel sta sulla frontiera contino, if ( IsAVoxelOnBoundary( nBlockLimits, nVoxIJK, true)) { double dSqEps = EPS_SMALL * EPS_SMALL ; Point3d ptFirstGrPt, ptSecondGrPt ; int nCount = 0 ; if ( SqDist( trTria.GetP( 0), ptVert) > dSqEps) { ptFirstGrPt = trTria.GetP( 0) ; nCount ++ ; } if ( SqDist( trTria.GetP( 1), ptVert) > dSqEps) { if ( nCount == 0) { ptFirstGrPt = trTria.GetP( 1) ; nCount ++ ; } else if ( nCount == 1) { ptSecondGrPt = trTria.GetP( 1) ; nCount ++ ; } } if ( SqDist( trTria.GetP( 2), ptVert) > dSqEps) { if ( nCount == 1) { ptSecondGrPt = trTria.GetP( 2) ; nCount ++ ; } } /* Altro modo: std::vector vPt ; for ( int nC = 0 ; nC < 3 ; ++ nC) if ( SqDist( trTria.GetP( 0), ptVert) > dSqEps) vPt.emplace_back( trTria.GetP( 0)) ; if ( vPt.size() == 2) ......*/ if ( nVoxIJK[0] == nBlockLimits[0]) { double dXGrid = ( nBlockLimits[0] + 0.5) * m_dStep ; if ( abs( ptFirstGrPt.x - dXGrid) < EPS_SMALL && abs( ptSecondGrPt.x - dXGrid) < EPS_SMALL) return true ; } if ( nVoxIJK[1] == nBlockLimits[2]) { double dYGrid = ( nBlockLimits[2] + 0.5) * m_dStep ; if ( abs( ptFirstGrPt.y - dYGrid) < EPS_SMALL && abs( ptSecondGrPt.y - dYGrid) < EPS_SMALL) return true ; } if ( nVoxIJK[2] == nBlockLimits[4]) { double dZGrid = ( nBlockLimits[4] + 0.5) * m_dStep ; if ( abs( ptFirstGrPt.z - dZGrid) < EPS_SMALL && abs( ptSecondGrPt.z - dZGrid) < EPS_SMALL) return true ; } if ( nVoxIJK[0] + 1 == nBlockLimits[1]) { double dXGrid = ( nBlockLimits[1] + 0.5) * m_dStep ; if ( abs( ptFirstGrPt.x - dXGrid) < EPS_SMALL && abs( ptSecondGrPt.x - dXGrid) < EPS_SMALL) return true ; } if ( nVoxIJK[1] + 1 == nBlockLimits[3]) { double dYGrid = ( nBlockLimits[3] + 0.5) * m_dStep ; if ( abs( ptFirstGrPt.y - dYGrid) < EPS_SMALL && abs( ptSecondGrPt.y - dYGrid) < EPS_SMALL) return true ; } if ( nVoxIJK[2] + 1 == nBlockLimits[5]) { double dZGrid = ( nBlockLimits[5] + 0.5) * m_dStep ; if ( abs( ptFirstGrPt.z - dZGrid) < EPS_SMALL && abs( ptSecondGrPt.z - dZGrid) < EPS_SMALL) return true ; } return true ; /*Altro modo for ( int nC = 0 ; nC < 3 ; ++ nC) { if ( nVoxIJK[nC] == nBlockLimits[2*nC]) { double dGrid = ( nBlockLimits[2*nC] + 0.5) * m_dStep ; if ( abs( ptFirstGrPt.v[nC] - dGrid) < EPS_SMALL && abs( ptSecondGrPt.v[nC] - dGrid) < EPS_SMALL) return true ; } if ( nVoxIJK[nC] + 1 == nBlockLimits[2*nC+1]) { double dGrid = ( nBlockLimits[2*nC+1] + 0.5) * m_dStep ; if ( abs( ptFirstGrPt.v[nC] - dGrid) < EPS_SMALL && abs( ptSecondGrPt.v[nC] - dGrid) < EPS_SMALL) return true ; } }*/ } // altrimenti non può starci il triangolo else return false ; } //---------------------------------------------------------------------------- bool VolZmap::IsATriangleOnBorder( const Triangle3d& trTria, const Point3d& ptVert, const int nBlockLimits[], const int nVoxIJK[], bool bBorderBox) const { // Se il voxel è di frontiera continuo, if ( bBorderBox) { // Se il triangolo ha area nulla non proseguiamo if ( trTria.GetArea() < EPS_SMALL) return false ; // Determino i punti del triangolo sulla griglia double dSqEps = EPS_SMALL * EPS_SMALL ; int nNotVert[2] ; int nCount = 0 ; for ( int nV = 0 ; nV < 3 ; ++ nV) { if ( SqDist( trTria.GetP( nV), ptVert) > dSqEps) { if ( nCount == 0) { nNotVert[nCount] = nV ; } else if ( nCount == 1) { nNotVert[nCount] = nV ; } nCount ++ ; } } // Se nCount != 2 deve esserci un errore if ( nCount != 2) return false ; // Punti del triangolo sulla griglia Point3d ptFirstGrPt = trTria.GetP( nNotVert[0]) ; Point3d ptSecondGrPt = trTria.GetP( nNotVert[1]) ; // Verifico se tali punti sono sula griglia for ( int nC = 0 ; nC < 3 ; ++ nC) { if ( nVoxIJK[nC] == nBlockLimits[2*nC]) { double dGrid = ( nBlockLimits[2*nC] + 0.5) * m_dStep ; if ( abs( ptFirstGrPt.v[nC] - dGrid) < EPS_SMALL && abs( ptSecondGrPt.v[nC] - dGrid) < EPS_SMALL) return true ; } if ( nVoxIJK[nC] + 1 == nBlockLimits[2*nC+1]) { double dGrid = ( nBlockLimits[2*nC+1] + 0.5) * m_dStep ; if ( abs( ptFirstGrPt.v[nC] - dGrid) < EPS_SMALL && abs( ptSecondGrPt.v[nC] - dGrid) < EPS_SMALL) return true ; } } return false ; } // altrimenti non ha senso continuare else return false ; }