//---------------------------------------------------------------------------- // 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 "\EgtDev\Include\EGkIntervals.h" #include "\EgtDev\Include\EgtNumUtils.h" #include "MC_Tables.h" using namespace std ; // ------------------------- VISUALIZZAZIONE -------------------------------------------------------------------------------------- //---------------------------------------------------------------------------- bool VolZmap::GetDexelLines( int nDir, int nPos1, int nPos2, POLYLINELIST& lstPL) const { // Controllo l'ammissibilità della griglia if ( nDir < 0 || nDir > 2) return false ; // Verifiche sugli indici if ( nPos1 < 0 || nPos1 >= int( m_nVNx[nDir]) || nPos2 < 0 || nPos2 >= int( m_nVNy[nDir])) return false ; int nPos = nPos1 + nPos2 * m_nVNx[nDir] ; if ( nPos < 0 || nPos >= int( m_TriZValues[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 partensa 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_TriZValues[nDir][nPos].size()) ; j += 2) { // aggiungo polilinea a lista lstPL.emplace_back() ; // inserisco punti estremi lstPL.back().AddUPoint( 0, ptP + m_TriZValues[nDir][nPos][j-1] * m_MapFrame[nDir].VersZ()) ; lstPL.back().AddUPoint( 1, ptP + m_TriZValues[nDir][nPos][j] * m_MapFrame[nDir].VersZ()) ; } 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_nVNx[0]) ; i += MAX_DIM_CHUNK) { int nDimChunkX = min( MAX_DIM_CHUNK, int( m_nVNx[0]) - i) ; for ( int j = 0 ; j < int( m_nVNy[0]) ; j += MAX_DIM_CHUNK) { int nDimChunkY = min( MAX_DIM_CHUNK, int( m_nVNy[0]) - j) ; GetChunkPrisms( i, j, nDimChunkX, nDimChunkY, MAX_DIM_CHUNK, lstTria) ; } } } else MarchingCubes( lstTria) ; return true ; } //---------------------------------------------------------------------------- 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_nVNx[0] ; if ( nPos > int( m_nVDim[0]) || int( m_TriZValues[0][nPos].size()) != 2) bIsSimple = false ; else if ( i == 0 && j == 0) { dBotZ = m_TriZValues[0][nPos][0] ; dTopZ = m_TriZValues[0][nPos][1] ; } else if ( abs( m_TriZValues[0][nPos][0] - dBotZ) > EPS_SMALL || abs( m_TriZValues[0][nPos][1] - 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_nVNx[0]) || nPos2 < 0 || nPos2 + nDim2 > int( m_nVNy[0])) return false ; int nPos = nPos1 + nPos2 * m_nVNx[0] ; if ( nPos < 0 || nPos >= int( m_nVDim[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_TriZValues[0][nPos][1] * m_MapFrame[0].VersZ() ; Vector3d vtDZb = m_TriZValues[0][nPos][0] * 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_nVNx[0] ; int nPosEst = ( nPos1 + nDim1 - 1 < int( m_nVNx[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_nVNx[0] + i ; int nPosNord = ( nPos2 + nDim2 - 1 < int( m_nVNy[0] - 1) ? nPosD + m_nVNx[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_nVNx[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_nVNx[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_nVNx[0]) || nPos2 < 0 || nPos2 >= int( m_nVNy[0])) return false ; int nPos = nPos1 + nPos2 * m_nVNx[0] ; if ( nPos < 0 || nPos >= int( m_nVDim[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_TriZValues[0][nPos].size()) ; i += 2) { Vector3d vtDZt = m_TriZValues[0][nPos][i] * m_MapFrame[0].VersZ() ; Vector3d vtDZb = m_TriZValues[0][nPos][i-1] * 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_nVNx[0] - 1) ? nPos + 1 : - 1) ; AddDexelSideFace( nPos, nPosEst, ptP2, ptP3, m_MapFrame[0].VersZ(), m_MapFrame[0].VersX(), lstTria) ; int nPosNord = ( nPos2 < int( m_nVNy[0] - 1) ? nPos + m_nVNx[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_nVNx[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_TriZValues[0][nPos].size()) ; i += 2) intFace.Add( m_TriZValues[0][nPos][i-1], m_TriZValues[0][nPos][i]) ; if ( nPosAdj > 0) { for ( int i = 1 ; i < int( m_TriZValues[0][nPosAdj].size()) ; i += 2) intFace.Subtract( m_TriZValues[0][nPosAdj][i-1], m_TriZValues[0][nPosAdj][i]) ; } 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) { // Limiti superiori sui tre indici unsigned int nLimI = m_nVNx[0] ; unsigned int nLimJ = m_nVNy[0] ; unsigned int nLimK = m_nVNy[1] ; struct Corner { int nCornerNumber ; unsigned int nI, nJ, nK ; } ; // Ciclo su tutti i voxel dello Zmap for ( unsigned int k = 0 ; k < nLimK ; ++ k) { for ( unsigned int i = 0 ; i < nLimI ; ++ i) { for ( unsigned int j = 0 ; j < nLimJ ; ++ j) { } } } return true ; } */ /* //---------------------------------------------------------------------------- bool VolZmap::MarchingCubes( TRIA3DLIST& lstTria) { unsigned int nLimI = m_nVNx[0] ; unsigned int nLimJ = m_nVNy[0] ; unsigned int nLimK = m_nVNy[1] ; // Ciclo su tutti i voxel dello Zmap for ( unsigned int k = 0 ; k < nLimK ; ++ k) { for ( unsigned int i = 0 ; i < nLimI ; ++ i) { for ( unsigned int j = 0 ; j < nLimJ ; ++ j) { // 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 ; bool CornerTF[8] = { false, false, false, false, false, false, false, false,} ; // Classificazione dei vertici: interni o esterni al materiale if ( IsThereMat( i, j, k)) { nIndex |= ( 1 << 0) ; CornerTF[0] = true ; } if ( IsThereMat( i + 1, j, k)) { nIndex |= ( 1 << 1) ; CornerTF[1] = true ; } if ( IsThereMat( i + 1, j + 1, k)) { nIndex |= ( 1 << 2) ; CornerTF[2] = true ; } if ( IsThereMat( i, j + 1, k)) { nIndex |= ( 1 << 3) ; CornerTF[3] = true ; } if ( IsThereMat( i, j, k + 1)) { nIndex |= ( 1 << 4) ; CornerTF[4] = true ; } if ( IsThereMat( i + 1, j, k + 1)) { nIndex |= ( 1 << 5) ; CornerTF[5] = true ; } if ( IsThereMat( i + 1, j + 1, k + 1)) { nIndex |= ( 1 << 6) ; CornerTF[6] = true ; } if ( IsThereMat( i, j + 1, k + 1)) { nIndex |= ( 1 << 7) ; CornerTF[7] = true ; } // Se vi è qualche intersezione fra segmenti e superficie // continuo altrimenti passo al prossimo voxel if ( EdgeTable[nIndex] != 0) { 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 i = 0 ; i < 12 ; ++ i) { // Se il segmento non attraversa la superficie // passo al successivo if ( ! ( EdgeTable[nIndex] & ( 1 << i))) continue ; int n1 = intersections[i][0]; int n2 = intersections[i][1]; // Determino con precisione il punto di intersezione sullo spigolo IntersPos( IndexCorner[n1], IndexCorner[n2], ptIntPoint[i]) ; } // Costruzione dei triangoli for ( int i = 0 ; TriangleTable[nIndex][i] != - 1 ; ++ i) { // Costruzione triangolo int i0 = TriangleTable[nIndex][i] ; int i1 = TriangleTable[nIndex][i+1] ; int i2 = TriangleTable[nIndex][i+2] ; Triangle3d CurrentTriangle ; Vector3d vtV1 = ptIntPoint[i1] - ptIntPoint[i0] ; Vector3d vtV2 = ptIntPoint[i2] - ptIntPoint[i0] ; Vector3d vtN = vtV1 ^ vtV2 ; vtN.Normalize() ; int nCorner = intersections[i0][0] ; Point3d ptCorner( IndexCorner[nCorner][0] * m_dStep, IndexCorner[nCorner][1] * m_dStep, IndexCorner[nCorner][2] * m_dStep) ; Vector3d vtT = ptCorner - ptIntPoint[i0] ; vtT.Normalize() ; if ( CornerTF[nCorner]) { if ( vtN * vtT < 0) vtN = - vtN ; } else { if ( vtN * vtT > 0) vtN = - vtN ; }*/ /* oppure: if( nIndex & ( 1 << nCorner)) { if( vtN * vtT < 0) } else { if ( vtN * vtT > 0) vtN = - vtN ; } */ /* CurrentTriangle.Set( ptIntPoint[i0], ptIntPoint[i1], ptIntPoint[i2], vtN) ; // Aggiungo triangolo lstTria.emplace_back( CurrentTriangle) ; } } } } } return true ; } //---------------------------------------------------------------------------- bool VolZmap::IsThereMat( unsigned int nI, unsigned int nJ, unsigned int nK) { double dZ[3] ; dZ[0] = ( nI + 0.5) * m_dStep ; dZ[1] = ( nK + 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_nVNx[nGrid] + nGrI ; unsigned int nDexSize = m_TriZValues[nGrid][nPos].size() ; unsigned int nIndex = 0 ; while ( nIndex < nDexSize) { if ( dZ[nGrid] > m_TriZValues[nGrid][nPos][nIndex] + EPS_SMALL || dZ[nGrid] < m_TriZValues[nGrid][nPos][nIndex + 1] - EPS_SMALL) { ++ nCount ; break ; } ++ nIndex ; } } if ( nCount > 1) return true ; return false ; } //---------------------------------------------------------------------------- bool VolZmap::IntersPos( int nVec1[], int nVec2[], Point3d & ptInt) { if ( nVec1[0] != nVec2[0]) { int nMinI = min( nVec1[0], nVec2[0]) ; int nMaxI = max( nVec1[0], nVec2[0]) ; double dMinX = nMinI * m_dStep ; double dMaxX = nMaxI * m_dStep ; unsigned int nDexel = nVec1[2] * m_nVNx[1] + nVec1[1] ; unsigned int nSize = m_TriZValues[1][nDexel].size() ; ptInt.y = nVec1[1] * m_dStep ; ptInt.z = nVec1[2] * m_dStep ; for ( unsigned int i = 0 ; i < nSize ; i += 2) { double dx1 = m_TriZValues[1][nDexel][i] ; double dx2 = m_TriZValues[1][nDexel][i+1] ; if ( dx1 < dMinX && dx2 > dMinX && dx2 < dMaxX) { ptInt.x = dx2 ; break ; } else if ( dx1 > dMinX && dx1 < dMaxX && dx2 > dMaxX) { ptInt.x = dx1 ; break ; } } } else if ( nVec1[1] != nVec2[1]) { int nMinJ = min( nVec1[1], nVec2[1]) ; int nMaxJ = max( nVec1[1], nVec2[1]) ; double dMinY = nMinJ * m_dStep ; double dMaxY = nMaxJ * m_dStep ; unsigned int nDexel = nVec1[2] * m_nVNx[2] + nVec1[0] ; unsigned int nSize = m_TriZValues[2][nDexel].size() ; ptInt.x = nVec1[0] * m_dStep ; ptInt.z = nVec1[2] * m_dStep ; for ( unsigned int j = 0 ; j < nSize ; j += 2) { double dy1 = m_TriZValues[2][nDexel][j] ; double dy2 = m_TriZValues[2][nDexel][j+1] ; if ( dy1 < dMinY && dy2 > dMinY && dy2 < dMaxY) { ptInt.y = dy2 ; break ; } else if ( dy1 > dMinY && dy1 < dMaxY && dy2 > dMaxY) { ptInt.y = dy1 ; break ; } } } else if ( nVec1[2] != nVec2[2]) { int nMinK = min( nVec1[2], nVec2[2]) ; int nMaxK = max( nVec1[2], nVec2[2]) ; double dMinZ = nMinK * m_dStep ; double dMaxZ = nMaxK * m_dStep ; unsigned int nDexel = nVec1[1] * m_nVNx[0] + nVec1[0] ; unsigned int nSize = m_TriZValues[0][nDexel].size() ; ptInt.x = nVec1[0] * m_dStep ; ptInt.y = nVec1[1] * m_dStep ; for ( unsigned int k = 0 ; k < nSize ; k += 2) { double dz1 = m_TriZValues[0][nDexel][k] ; double dz2 = m_TriZValues[0][nDexel][k+1] ; if ( dz1 < dMinZ && dz2 > dMinZ && dz2 < dMaxZ) { ptInt.z = dz2 ; break ; } else if ( dz1 > dMinZ && dz1 < dMaxZ && dz2 > dMaxZ) { ptInt.z = dz1 ; break ; } } } return true ; } */ // Prova bool function( const Vector3d & vtV) { if ( 100 - vtV * vtV < 0) return true ; return false ; } bool midpoint( int nVec1[], int nVec2[], Point3d & ptInt, double dStep) { Point3d pt1( ( nVec1[0] + 0.5) * dStep, ( nVec1[1] + 0.5) * dStep, ( nVec1[2] + 0.5) * dStep) ; Point3d pt2( ( nVec2[0] + 0.5) * dStep, ( nVec2[1] + 0.5) * dStep, ( nVec2[2] + 0.5) * dStep) ; ptInt = pt1 + 0.5 * ( pt2 - pt1) ; return true ; } //---------------------------------------------------------------------------- bool VolZmap::MarchingCubes( TRIA3DLIST& lstTria) const { Point3d ptMapOrig = m_MapFrame[0].Orig() ; int nLimI = int( m_nVNx[0]) ; int nLimJ = int( m_nVNy[0]) ; int nLimK = int( m_nVNy[1]) ; // Ciclo su tutti i voxel dello Zmap for ( int i = - 1 ; i < nLimI ; ++ i) { for ( int j = - 1 ; j < nLimJ ; ++ j) { for ( int k = - 1 ; k < nLimK ; ++ k) { if ( i == 0 && j == 6 && k == 8) { double bau = 1 ; } // Indici i,j,k dei vertici int IndexCorner[8][3] = { { i, j, k}, { i + 1, j, k}, { i + 1, j + 1, k}, { i, j + 1, k}, { i, j, k + 1}, { i + 1, j, k + 1}, { i + 1, j + 1, k + 1}, { i, j + 1, k + 1} } ; int nIndex = 0 ; // Classificazione dei vertici: interni o esterni al materiale if ( IsThereMat( i, j, k)) nIndex |= ( 1 << 0) ; if ( IsThereMat( i + 1, j, k)) nIndex |= ( 1 << 1) ; if ( IsThereMat( i + 1, j + 1, k)) nIndex |= ( 1 << 2) ; if ( IsThereMat( i, j + 1, k)) nIndex |= ( 1 << 3) ; if ( IsThereMat( i, j, k + 1)) nIndex |= ( 1 << 4) ; if ( IsThereMat( i + 1, j, k + 1)) nIndex |= ( 1 << 5) ; if ( IsThereMat( i + 1, j + 1, k + 1)) nIndex |= ( 1 << 6) ; if ( IsThereMat( i, j + 1, k + 1)) nIndex |= ( 1 << 7) ; // Se vi è qualche intersezione fra segmenti e superficie // continuo altrimenti passo al prossimo voxel if ( EdgeTable[nIndex] == 0) continue ; static int intersections[12][2] = { { 0, 1 }, { 1, 2 }, { 2, 3 }, { 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]) ; // midpoint( IndexCorner[n1], IndexCorner[n2], ptIntPoint[EdgeIndex], m_dStep) ; ptIntPoint[EdgeIndex] = ptIntPoint[EdgeIndex] ; ptIntPoint[EdgeIndex].ToGlob( m_MapFrame[0]) ; } // Costruzione dei triangoli for ( int TriIndex = 0 ; TriangleTable[nIndex][TriIndex] != - 1 ; TriIndex += 3) { // Costruzione triangolo int i0 = TriangleTable[nIndex][TriIndex + 2] ; int i1 = TriangleTable[nIndex][TriIndex + 1] ; int i2 = TriangleTable[nIndex][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::IsThereMat( int nI, int nJ, int nK) const { if ( nI == - 1 || nI == m_nVNx[0] || nJ == - 1 || nJ == m_nVNy[0] || nK == - 1 || nK == m_nVNy[1]) return false ; 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_nVNx[nGrid] + nGrI ; size_t nDexSize = m_TriZValues[nGrid][nPos].size() ; unsigned int nIndex = 0 ; while ( nIndex < nDexSize) { if ( dZ[nGrid] > m_TriZValues[nGrid][nPos][nIndex] && dZ[nGrid] < m_TriZValues[nGrid][nPos][nIndex + 1]) { ++ nCount ; break ; } nIndex += 2 ; } } //if ( nCount > 0) if ( nCount == 3) return true ; return false ; } //---------------------------------------------------------------------------- bool VolZmap::IsThereMat( const int nMatr[][3], int nNum, double & dHx, double & dHy, double & dHz) const { return true ; } //---------------------------------------------------------------------------- bool VolZmap::IntersPos( int nVec1[], int nVec2[], Point3d & ptInt) const { 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 ; unsigned int nDexel = nVec1[2] * m_nVNx[1] + nVec1[1] ; size_t nSize = m_TriZValues[1][nDexel].size() ; ptInt.y = ( nVec1[1] + 0.5) * m_dStep ; ptInt.z = ( nVec1[2] + 0.5) * m_dStep ; unsigned int i ; for ( i = 0 ; i < nSize ; i += 2) { double dx1 = m_TriZValues[1][nDexel][i] ; double dx2 = m_TriZValues[1][nDexel][i+1] ; if ( dx1 <= dMinX && dx2 >= dMinX && dx2 <= dMaxX) { ptInt.x = dx2 ; break ; } else if ( dx1 >= dMinX && dx1 <= dMaxX && dx2 >= dMaxX) { ptInt.x = dx1 ; break ; } } if ( i == nSize) 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 ; // 2 0 unsigned int nDexel = nVec1[0] * m_nVNx[2] + nVec1[2] ; size_t nSize = m_TriZValues[2][nDexel].size() ; ptInt.x = ( nVec1[0] + 0.5) * m_dStep ; ptInt.z = ( nVec1[2] + 0.5) * m_dStep ; unsigned int j ; for ( j = 0 ; j < nSize ; j += 2) { double dy1 = m_TriZValues[2][nDexel][j] ; double dy2 = m_TriZValues[2][nDexel][j+1] ; if ( dy1 <= dMinY && dy2 >= dMinY && dy2 <= dMaxY) { ptInt.y = dy2 ; break ; } else if ( dy1 >= dMinY && dy1 <= dMaxY && dy2 >= dMaxY) { ptInt.y = dy1 ; break ; } } if ( j == nSize) 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 ; unsigned int nDexel = nVec1[1] * m_nVNx[0] + nVec1[0] ; size_t nSize = m_TriZValues[0][nDexel].size() ; ptInt.x = ( nVec1[0] + 0.5) * m_dStep ; ptInt.y = ( nVec1[1] + 0.5) * m_dStep ; unsigned int k ; for ( k = 0 ; k < nSize ; k += 2) { double dz1 = m_TriZValues[0][nDexel][k] ; double dz2 = m_TriZValues[0][nDexel][k+1] ; if ( dz1 <= dMinZ && dz2 >= dMinZ && dz2 <= dMaxZ) { ptInt.z = dz2 ; break ; } else if ( dz1 >= dMinZ && dz1 <= dMaxZ && dz2 >= dMaxZ) { ptInt.z = dz1 ; break ; } } if ( k == nSize) ptInt.z = 0.5 * ( dMinZ + dMaxZ) ; } return true ; } /* Point3d ptProva[8] ; Vector3d vtProva[8] ; ptProva[0].x = ( i + 0.5) * m_dStep ; ptProva[0].y = ( j + 0.5) * m_dStep ; ptProva[0].z = ( k + 0.5) * m_dStep ; ptProva[1].x = ( i + 1.5) * m_dStep ; ptProva[1].y = ( j + 0.5) * m_dStep ; ptProva[1].z = ( k + 0.5) * m_dStep ; ptProva[2].x = ( i + 1.5) * m_dStep ; ptProva[2].y = ( j + 1.5) * m_dStep ; ptProva[2].z = ( k + 0.5) * m_dStep ; ptProva[3].x = ( i + 0.5) * m_dStep ; ptProva[3].y = ( j + 1.5) * m_dStep ; ptProva[3].z = ( k + 0.5) * m_dStep ; ptProva[4].x = ( i + 0.5) * m_dStep ; ptProva[4].y = ( j + 0.5) * m_dStep ; ptProva[4].z = ( k + 1.5) * m_dStep ; ptProva[5].x = ( i + 1.5) * m_dStep ; ptProva[5].y = ( j + 0.5) * m_dStep ; ptProva[5].z = ( k + 1.5) * m_dStep ; ptProva[6].x = ( i + 1.5) * m_dStep ; ptProva[6].y = ( j + 1.5) * m_dStep ; ptProva[6].z = ( k + 1.5) * m_dStep ; ptProva[7].x = ( i + 0.5) * m_dStep ; ptProva[7].y = ( j + 1.5) * m_dStep ; ptProva[7].z = ( k + 1.5) * m_dStep ; Point3d ptC( 10, 10, 10) ; for ( int boh = 0 ; boh < 8 ; ++ boh) vtProva[boh] = ptProva[boh] - ptC ; */ /* if ( function( vtProva[0])) nIndex |= ( 1 << 0) ; if ( function( vtProva[1])) nIndex |= ( 1 << 1) ; if ( function( vtProva[2])) nIndex |= ( 1 << 2) ; if ( function( vtProva[3])) nIndex |= ( 1 << 3) ; if ( function( vtProva[4])) nIndex |= ( 1 << 4) ; if ( function( vtProva[5])) nIndex |= ( 1 << 5) ; if ( function( vtProva[6])) nIndex |= ( 1 << 6) ; if ( function( vtProva[7])) nIndex |= ( 1 << 7) ; */ /* Vector3d vtV1 = ptIntPoint[i1] - ptIntPoint[i0] ; Vector3d vtV2 = ptIntPoint[i2] - ptIntPoint[i0] ; Vector3d vtN = vtV1 ^ vtV2 ; vtN.Normalize() ; int nCorner = intersections[i0][0] ; Point3d ptCorner( IndexCorner[nCorner][0] * m_dStep, IndexCorner[nCorner][1] * m_dStep, IndexCorner[nCorner][2] * m_dStep) ; Vector3d vtT = ptCorner - ptIntPoint[i0] ; vtT.Normalize() ; if( nIndex & ( 1 << nCorner)) { if( vtN * vtT < 0) vtN = - vtN ; } else { if ( vtN * vtT > 0) vtN = - vtN ; }*/