From 38590bf38112d1e85491f6eee977bdffbebdfca9 Mon Sep 17 00:00:00 2001 From: Dario Sassi Date: Mon, 18 Sep 2017 06:14:33 +0000 Subject: [PATCH] EgtGeomKernel 1.8i2 : - migliorie a Zmap (grandi rettangoli con EMC). --- EgtGeomKernel.rc | Bin 11710 -> 11710 bytes VolTriZmapGraphics.cpp | 1670 ++++++++++++++++++++++++---------------- VolZmap.h | 54 +- 3 files changed, 1014 insertions(+), 710 deletions(-) diff --git a/EgtGeomKernel.rc b/EgtGeomKernel.rc index 6037d2834ac2592ecb904eb9dea3eeb451914846..79d7dd073a096a1c032355ad4d39e1769b31311e 100644 GIT binary patch delta 94 zcmdlNy)SyhFE&P_&A-_cnHh~HD{|{@_Trkr0u;H;XNwSVW8B;$>;>dw2zN+>g;Df- LFmBFL4&ed-Vf-5U delta 94 zcmdlNy)SyhFE&QQ&A-_cnHdcyD{|{@_Trkr0u;H;XNwSVW8B;$>;>dw2zN+>g;Df- LFmBFL4&ed-U|t&X diff --git a/VolTriZmapGraphics.cpp b/VolTriZmapGraphics.cpp index a2e95a4..865099c 100644 --- a/VolTriZmapGraphics.cpp +++ b/VolTriZmapGraphics.cpp @@ -34,6 +34,8 @@ struct VectorField { Vector3d vtNorm ; } ; +// ------------------------- UNORDERED MAP PER LA GESTIONE DEI TRIANGOLI GRANDI --------------------------------------------------- + // ------------------------- TABELLA BLOCCHI ADIACENTI ---------------------------------------------------------------------------- static int NeighbourTable[8][4] = { {0, -1, -1, -1}, @@ -48,6 +50,7 @@ static int NeighbourTable[8][4] = { // ------------------------- FUNZIONE TEST SULLE NORMALI -------------------------------------------------------------------------- enum FatureType { NO_FEATURE = 0, CORNER = 1, EDGE = 2} ; +enum CanonicDir { X_PLUS = 1, X_MINUS = -1, Y_PLUS = 2, Y_MINUS = -2, Z_PLUS = 3, Z_MINUS = -3} ; //---------------------------------------------------------------------------- int @@ -68,7 +71,7 @@ TestOnNormal( const VectorField CompoVert[], int nCompoElem) } // Se la massima deviazione non supera il limite non è feature - const double SHARP_COS = 0.9 ; // 0.8 ; + const double SHARP_COS = 0.9 ; if ( dMinCosTheta >= SHARP_COS) return NO_FEATURE ; @@ -84,13 +87,50 @@ TestOnNormal( const VectorField CompoVert[], int nCompoElem) dMaxAbsCos = dAbsCurrentCos ; } // se esiste normale diretta quasi come potenziale edge, allora corner - const double CORNER_COS = 0.7 ; // 0.7 ; // 0.5 ; + const double CORNER_COS = 0.7 ; if ( dMaxAbsCos > CORNER_COS) return CORNER ; else return EDGE ; } +//---------------------------------------------------------------------------- +bool +CanonicPlaneTest( const VectorField CompoVert[], int nDir, double& dPos) +{ + // Verifico posizione dei punti + int nI ; + switch ( nDir) { + case X_PLUS : case X_MINUS : nI = 0 ; break ; + case Y_PLUS : case Y_MINUS : nI = 1 ; break ; + case Z_PLUS : case Z_MINUS : nI = 2 ; break ; + } + double dPos0 = CompoVert[0].ptInt.v[nI] ; + double dPos1 = CompoVert[1].ptInt.v[nI] ; + double dPos2 = CompoVert[2].ptInt.v[nI] ; + double dPos3 = CompoVert[3].ptInt.v[nI] ; + dPos = ( dPos0 + dPos1 + dPos2 + dPos3) / 4 ; + if ( abs( dPos0 - dPos) > EPS_SMALL || abs( dPos1 - dPos) > EPS_SMALL || + abs( dPos2 - dPos) > EPS_SMALL || abs( dPos3 - dPos) > EPS_SMALL) + return false ; + // Verifico direzione normali + Vector3d vtN ; + switch ( nDir) { + case X_PLUS : vtN = X_AX ; break ; + case X_MINUS : vtN = -X_AX ; break ; + case Y_PLUS : vtN = Y_AX ; break ; + case Y_MINUS : vtN = - Y_AX ; break ; + case Z_PLUS : vtN = Z_AX ; break ; + case Z_MINUS : vtN = - Z_AX ; break ; + } + for ( int i = 0 ; i < 4 ; ++ i) { + if ( CompoVert[i].vtNorm * vtN < 0.99) + return false ; + } + // Superati tutti i test + return true ; +} + //---------------------------------------------------------------------------- bool DotTest( const VectorField CompoVert[], int nCompoElem, Vector3d& vtAvg, double dThreshold = 0) @@ -118,10 +158,8 @@ DotTest( const VectorField CompoVert[], int nCompoElem, Vector3d& vtAvg, double return true ; } - - + // ------------------------- VISUALIZZAZIONE -------------------------------------------------------------------------------------- - //---------------------------------------------------------------------------- bool VolZmap::GetDexelLines( int nDir, int nPos1, int nPos2, POLYLINELIST& lstPL) const @@ -360,7 +398,7 @@ VolZmap::GetTriangles( bool bAllBlocks, INTVECTOR& nModifiedBlocks, TRIA3DLISTVE 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 * EPS_SMALL) { + if ( m_InterBlockTria[t][t1].vCompoTria[t2][t3].GetArea() > SQ_EPS_SMALL) { vLstTria[nPos].emplace_back( InterBlockTria[t][t1].vCompoTria[t2][t3]) ; } } @@ -607,7 +645,15 @@ VolZmap::MarchingCubes( int nBlock, TRIA3DLIST& lstTria) const for ( int i = nStartI ; i < nEndI ; ++ i) { for ( int j = nStartJ ; j < nEndJ ; ++ j) { for ( int k = nStartK ; k < nEndK ; ++ k) { - + + // Classificazione dei vertici: interni o esterni al materiale + int nIndex = CalcIndex( i, j , k) ; + + // Se vi è qualche intersezione fra segmenti e superficie + // continuo altrimenti passo al prossimo voxel + if ( EdgeTable[nIndex] == 0) + continue ; + // Indici i,j,k dei vertici int IndexCorner[8][3] = { { i, j, k}, @@ -618,31 +664,7 @@ VolZmap::MarchingCubes( int nBlock, TRIA3DLIST& lstTria) const { 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 }, @@ -657,11 +679,11 @@ VolZmap::MarchingCubes( int nBlock, TRIA3DLIST& lstTria) const continue ; int n1 = intersections[EdgeIndex][0] ; - int n2 = intersections[EdgeIndex][1] ; + 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]) ; } @@ -674,13 +696,13 @@ VolZmap::MarchingCubes( int nBlock, TRIA3DLIST& lstTria) const 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) ; @@ -695,9 +717,10 @@ VolZmap::MarchingCubes( int nBlock, TRIA3DLIST& lstTria) const } //---------------------------------------------------------------------------- -bool +bool VolZmap::ExtMarchingCubes( int nBlock, TRIA3DLIST& lstTria, TriHolder& triHold) const -{ +{ + // Controllo sulla validità del blocco if ( nBlock < 0 || nBlock >= int( m_nNumBlock)) return false ; @@ -709,15 +732,32 @@ VolZmap::ExtMarchingCubes( int nBlock, TRIA3DLIST& lstTria, TriHolder& triHold) int nLimits[6] ; GetBlockLimitsIJK( nIJK, nLimits) ; + // Unordered Map per la riduzione del numero di triangoli + int nDim = m_nDexNumPBlock * m_nDexNumPBlock ; + VoxelContainer VoxContXZInf( nDim) ; + VoxelContainer VoxContXZSup( nDim) ; + VoxelContainer VoxContXYInf( nDim) ; + VoxelContainer VoxContXYSup( nDim) ; + VoxelContainer VoxContYZInf( nDim) ; + VoxelContainer VoxContYZSup( nDim) ; + // 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) { - + for ( int k = nLimits[4] ; k < nLimits[5] ; ++ k) { + + // Classificazione dei vertici: interni o esterni al materiale + int nIndex = CalcIndex( i, j , k) ; + + // Se vi è qualche intersezione fra segmenti e superficie + // continuo altrimenti passo al prossimo voxel. + if ( EdgeTable[nIndex] == 0) + continue ; + // Riconoscimento dei voxel di frontiera int nVoxIndexes[3] = { i, j, k} ; bool bBoundary = IsAVoxelOnBoundary( nLimits, nVoxIndexes, true) ; - + // Indici i,j,k dei vertici int IndexCorner[8][3] = { { i, j, k}, @@ -728,31 +768,7 @@ VolZmap::ExtMarchingCubes( int nBlock, TRIA3DLIST& lstTria, TriHolder& triHold) { 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 }, { 3, 2 }, { 0, 3 }, { 4, 5 }, { 5, 6 }, @@ -761,46 +777,37 @@ VolZmap::ExtMarchingCubes( int nBlock, TRIA3DLIST& lstTria, TriHolder& triHold) // Array 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. + // vero se i campi sono regolari 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))) + if ( ( EdgeTable[nIndex] & ( 1 << EdgeIndex)) == 0) continue ; int n1 = intersections[EdgeIndex][0] ; - int n2 = intersections[EdgeIndex][1] ; + 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, + 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 - // in caso di configurazione standard. + + // Determino il numero di componenti connesse nel voxel in caso di configurazione standard int nComponents = TriangleTableEn[nIndex][1][0] ; // Matrici di campi vettoriali: - // CompoVert[i] ha i vertici della base del triangle fan - // della (i+1)-esima componente connessa; - // CompoTriVert[i] ha i vertici di tutti i triangoli, nel - // nel caso di assenza di sharp feature, della (i+1)-esima - // componente connessa. + // CompoVert[i] ha i vertici della base del triangle fan della (i+1)-esima componente connessa; + // CompoTriVert[i] ha i vertici di tutti i triangoli, nel caso di assenza di sharp feature, + // della (i+1)-esima componente connessa. VectorField CompoVert[6][12] ; VectorField CompoTriVert[6][17] ; @@ -809,8 +816,7 @@ VolZmap::ExtMarchingCubes( int nBlock, TRIA3DLIST& lstTria, TriHolder& triHold) // della base del fan della (i+1)-esima componente connessa. int nVertComp[6] ; - // Matrice di indici dei punti: serve per - // la gestione del caso + // Matrice di indici dei punti: serve per la gestione del caso int nIndArrey[6][4] ; int nExtTabOff = nComponents ; @@ -818,20 +824,20 @@ VolZmap::ExtMarchingCubes( int nBlock, TRIA3DLIST& lstTria, TriHolder& triHold) // Carico le matrici CompoVert e CompoTriVert for ( int nCompCount = 1 ; nCompCount <= nComponents ; ++ nCompCount) { - + // Numero vertici per componenti nVertComp[nCompCount - 1] = TriangleTableEn[nIndex][1][nCompCount] ; - + // Riempio il nCompCount-esimo vettore di vertici della base del fan for ( int nVertCount = 0 ; nVertCount < nVertComp[nCompCount - 1] ; ++ nVertCount) CompoVert[nCompCount - 1][nVertCount] = VecField[TriangleTableEn[nIndex][1][nVertCount + nExtTabOff + 1]] ; - + // Serve per la gestione del caso ... if ( nVertComp[nCompCount - 1] == 4) { for ( int nVertCount = 0 ; nVertCount < nVertComp[nCompCount - 1] ; ++ nVertCount) nIndArrey[nCompCount - 1][nVertCount] = TriangleTableEn[nIndex][1][nVertCount + nExtTabOff + 1] ; } - + // Riempio il nCompCount-esimo vettore di vertici dei triangoli in assenza di // sharp feature: in una mesh di triangoli con n vertici vi sono n - 2 triangoli. for ( int nVertCount = 0 ; nVertCount < 3 * ( nVertComp[nCompCount - 1] - 2) ; nVertCount += 3) { @@ -839,19 +845,81 @@ VolZmap::ExtMarchingCubes( int nBlock, TRIA3DLIST& lstTria, TriHolder& triHold) CompoTriVert[nCompCount - 1][nVertCount+1] = VecField[TriangleTableEn[nIndex][0][nStdTabOff + nVertCount+1]] ; CompoTriVert[nCompCount - 1][nVertCount+2] = VecField[TriangleTableEn[nIndex][0][nStdTabOff + nVertCount]] ; } - - // Aggiorno gli offsets per raggiungere i + + // Aggiorno gli offsets per raggiungere i // vertici della componente successiva. nExtTabOff += nVertComp[nCompCount - 1] ; - nStdTabOff += 3 * ( nVertComp[nCompCount - 1] - 2) ; + nStdTabOff += 3 * ( nVertComp[nCompCount - 1] - 2) ; } - + + // Controllo se il voxel ha una sola faccia che giace in un piano canonico e quindi ha gestione speciale + // Faccia XY normale Z+ + if ( nIndex == 15) { + double dPos ; + if ( CanonicPlaneTest( CompoVert[0], Z_PLUS, dPos)) { + int nN ; + GetVoxNFromIJK( i, j, k, nN) ; + VoxContXYSup.emplace( nN, dPos) ; + continue ; + } + } + // Faccia YZ normale X+ + else if ( nIndex == 153) { + double dPos ; + if ( CanonicPlaneTest( CompoVert[0], X_PLUS, dPos)) { + int nN ; + GetVoxNFromIJK( i, j, k, nN) ; + VoxContYZSup.emplace( nN, dPos) ; + continue ; + } + } + // Faccia ZX normale Y+ + else if ( nIndex == 51) { + double dPos ; + if ( CanonicPlaneTest( CompoVert[0], Y_PLUS, dPos)) { + int nN ; + GetVoxNFromIJK( i, j, k, nN) ; + VoxContXZSup.emplace( nN, dPos) ; + continue ; + } + } + // Faccia YX normale Z- + else if ( nIndex == 240) { + double dPos ; + if ( CanonicPlaneTest( CompoVert[0], Z_MINUS, dPos)) { + int nN ; + GetVoxNFromIJK( i, j, k, nN) ; + VoxContXYInf.emplace( nN, dPos) ; + continue ; + } + } + // Faccia ZY normale X- + else if ( nIndex == 102) { + double dPos ; + if ( CanonicPlaneTest( CompoVert[0], X_MINUS, dPos)) { + int nN ; + GetVoxNFromIJK( i, j, k, nN) ; + VoxContYZInf.emplace( nN, dPos) ; + continue ; + } + } + // Faccia XZ normale Y- + else if ( nIndex == 204) { + double dPos ; + if ( CanonicPlaneTest( CompoVert[0], Y_MINUS, dPos)) { + int nN ; + GetVoxNFromIJK( i, j, k, nN) ; + VoxContXZInf.emplace( nN, dPos) ; + continue ; + } + } + // Test sulla topologia: dal momento che il nostro test // si fonda sugli angoli compresi fra le normali, esso ha // senso solo se il campo è regolare. if ( bReg) { if ( nAllConfig[nIndex] == 3) { - + Vector3d vtCmpAvg0, vtCmpAvg1 ; // Verifico se i versori delle componenti sono tutti @@ -876,12 +944,12 @@ VolZmap::ExtMarchingCubes( int nBlock, TRIA3DLIST& lstTria, TriHolder& triHold) double dThreshold = 0.7 ; if ( ! ( bTest0 && bTest1) || ( bTest0 && bTest1 && dScProd > dThreshold)) { - + int nt = 0 ; while ( nIndexVsIndex3[nt][0] != nIndex) ++ nt ; - + int nRotCase = nIndexVsIndex3[nt][1] ; nComponents = Cases3Plus[nRotCase][1][0] ; @@ -895,7 +963,7 @@ VolZmap::ExtMarchingCubes( int nBlock, TRIA3DLIST& lstTria, TriHolder& triHold) // Numero vertici per componenti nVertComp[nC - 1] = Cases3Plus[nRotCase][1][nC] ; - + // Matrice dei vertici della base del fan for ( int nFanVert = 0 ; nFanVert < nVertComp[nC - 1] ; ++ nFanVert) @@ -903,29 +971,30 @@ VolZmap::ExtMarchingCubes( int nBlock, TRIA3DLIST& lstTria, TriHolder& triHold) // Matrici dei vertici dei triangoli in assenza di sharp feature for ( int nTriVert = 0 ; nTriVert < 3 * ( nVertComp[nC - 1] - 2) ; nTriVert += 3) { - + CompoTriVert[nC - 1][nTriVert] = VecField[Cases3Plus[nRotCase][0][nStdTabOff + nTriVert+2]] ; CompoTriVert[nC - 1][nTriVert+1] = VecField[Cases3Plus[nRotCase][0][nStdTabOff + nTriVert+1]] ; CompoTriVert[nC - 1][nTriVert+2] = VecField[Cases3Plus[nRotCase][0][nStdTabOff + nTriVert]] ; } - + // Aggiorno gli offsets per raggiungere i // vertici della componente successiva. nExtTabOff += nVertComp[nC - 1] ; - nStdTabOff += 3 * ( nVertComp[nC - 1] - 2) ; - } + nStdTabOff += 3 * ( nVertComp[nC - 1] - 2) ; + } - } + } } - else if ( nAllConfig[nIndex] == 6) { - + + else if ( nAllConfig[nIndex] == 6) { + // Procedura analoga a quella della configurazione 3 Vector3d vtCmpAvg0, vtCmpAvg1 ; bool bTest0 = DotTest( CompoVert[0], 4, vtCmpAvg0, 0.7) ; bool bTest1 = DotTest( CompoVert[1], 3, vtCmpAvg1, 0.7) ; - + double dScProd = - 2 ; if ( bTest0 && bTest1) @@ -942,7 +1011,6 @@ VolZmap::ExtMarchingCubes( int nBlock, TRIA3DLIST& lstTria, TriHolder& triHold) int nRotCase = nIndexVsIndex6[nt][1] ; - // Costruzione dei triangoli for ( int TriIndex = 0 ; TriIndex < 15 ; TriIndex += 3) { @@ -952,17 +1020,18 @@ VolZmap::ExtMarchingCubes( int nBlock, TRIA3DLIST& lstTria, TriHolder& triHold) int i2 = Cases6Plus[nRotCase][TriIndex] ; Triangle3d CurrentTriangle ; - + // Il triangolo è pronto CurrentTriangle.Set( VecField[i0].ptInt, VecField[i1].ptInt, VecField[i2].ptInt) ; bool bV = CurrentTriangle.Validate( true) ; - // Aggiungo alla lista + // Aggiungo alla lista lstTria.emplace_back( CurrentTriangle) ; } continue ; } - } + } + else if ( nAllConfig[nIndex] == 10) { Vector3d vtCmpAvg0, vtCmpAvg1 ; @@ -977,10 +1046,9 @@ VolZmap::ExtMarchingCubes( int nBlock, TRIA3DLIST& lstTria, TriHolder& triHold) if ( ! ( bTest0 && bTest1)) { int nt = 0 ; - while ( nIndexVsIndex10[nt][0] != nIndex) ++ nt ; - + int nRotCase = nIndexVsIndex10[nt][1] ; // Riaggiorno gli offsets @@ -989,33 +1057,28 @@ VolZmap::ExtMarchingCubes( int nBlock, TRIA3DLIST& lstTria, TriHolder& triHold) // Modifico le matrici for ( int nC = 1 ; nC <= 2 ; ++ nC) { - // Numero vertici per componenti nVertComp[nC - 1] = Cases10Plus[nRotCase][1][nC] ; - + // Matrice dei vertici della base del fan for ( int nFanVert = 0 ; nFanVert < 4 ; ++ nFanVert) - CompoVert[nC - 1][nFanVert] = VecField[Cases10Plus[nRotCase][1][nFanVert + nExtTabOff + 1]] ; // Matrici dei vertici dei triangoli in assenza di sharp feature for ( int nTriVert = 0 ; nTriVert < 6 ; nTriVert += 3) { - CompoTriVert[nC - 1][nTriVert] = VecField[Cases10Plus[nRotCase][0][nStdTabOff + nTriVert+2]] ; CompoTriVert[nC - 1][nTriVert+1] = VecField[Cases10Plus[nRotCase][0][nStdTabOff + nTriVert+1]] ; CompoTriVert[nC - 1][nTriVert+2] = VecField[Cases10Plus[nRotCase][0][nStdTabOff + nTriVert]] ; } - - // Aggiorno gli offsets per raggiungere i - // vertici della componente successiva. + + // Aggiorno gli offsets per raggiungere i vertici della componente successiva. nExtTabOff += nVertComp[nC - 1] ; - nStdTabOff += 3 * ( nVertComp[nC - 1] - 2) ; - } - } + nStdTabOff += 3 * ( nVertComp[nC - 1] - 2) ; + } + } } } - - + // Numero di feature nel voxel: al più vi è una feature per componente connessa. int nInnerFeatureInVoxel = 0 ; int nBorderFeatureInVoxel = 0 ; @@ -1027,7 +1090,7 @@ VolZmap::ExtMarchingCubes( int nBlock, TRIA3DLIST& lstTria, TriHolder& triHold) // Se i componenti sono regolari valuto le normali per stabilire se eseguire ExtMC o MC if ( bReg) nFeatureType = TestOnNormal( CompoVert[nCompCount - 1], nVertComp[nCompCount - 1]) ; - + // Controllo per il caso piano su una griglia // con versori normali a due a due paralleli. bool bGridControl = true ; @@ -1035,10 +1098,10 @@ VolZmap::ExtMarchingCubes( int nBlock, TRIA3DLIST& lstTria, TriHolder& triHold) if ( nFeatureType != NO_FEATURE) { if ( nVertComp[nCompCount - 1] == 4) { - + // Ordino i 4 indici in senso crescente for ( int nSrtInd1 = 0 ; nSrtInd1 < nVertComp[nCompCount - 1] - 1 ; ++ nSrtInd1) { - for ( int nSrtInd2 = nSrtInd1 + 1 ; nSrtInd2 < nVertComp[nCompCount - 1] ; ++ nSrtInd2) { + for ( int nSrtInd2 = nSrtInd1 + 1 ; nSrtInd2 < nVertComp[nCompCount - 1] ; ++ nSrtInd2) { if ( nIndArrey[nCompCount - 1][nSrtInd1] > nIndArrey[nCompCount - 1][nSrtInd2]) swap( nIndArrey[nCompCount - 1][nSrtInd1], nIndArrey[nCompCount - 1][nSrtInd2]) ; } @@ -1060,18 +1123,12 @@ VolZmap::ExtMarchingCubes( int nBlock, TRIA3DLIST& lstTria, TriHolder& triHold) nIndArrey[nCompCount - 1][2] == 10 && nIndArrey[nCompCount - 1][3] == 11) || ( nIndArrey[nCompCount - 1][0] == 5 && nIndArrey[nCompCount - 1][1] == 7 && nIndArrey[nCompCount - 1][2] == 8 && nIndArrey[nCompCount - 1][3] == 9 )) { - + VectorField LocVecF[12], LocCompV[12] ; for ( int LocInd = 0 ; LocInd < 12 ; ++ LocInd) { - LocVecF[LocInd] = VecField[LocInd] ; LocCompV[LocInd] = CompoVert[nCompCount - 1][LocInd] ; - - LocVecF[LocInd].ptInt.ToLoc( m_MapFrame[0]) ; - LocVecF[LocInd].vtNorm.ToLoc( m_MapFrame[0]) ; - LocCompV[LocInd].ptInt.ToLoc( m_MapFrame[0]) ; - LocCompV[LocInd].vtNorm.ToLoc( m_MapFrame[0]) ; } if ( ( AreSameVectorApprox( LocVecF[nIndArrey[nCompCount - 1][0]].vtNorm, LocVecF[nIndArrey[nCompCount - 1][1]].vtNorm) && @@ -1088,17 +1145,17 @@ VolZmap::ExtMarchingCubes( int nBlock, TRIA3DLIST& lstTria, TriHolder& triHold) abs( LocCompV[1].ptInt.x - ptBarycenter.x) < EPS_SMALL && abs( LocCompV[2].ptInt.x - ptBarycenter.x) < EPS_SMALL && abs( LocCompV[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 ; } @@ -1109,17 +1166,17 @@ VolZmap::ExtMarchingCubes( int nBlock, TRIA3DLIST& lstTria, TriHolder& triHold) abs( LocCompV[1].ptInt.y - ptBarycenter.y) < EPS_SMALL && abs( LocCompV[2].ptInt.y - ptBarycenter.y) < EPS_SMALL && abs( LocCompV[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 ; } @@ -1130,17 +1187,17 @@ VolZmap::ExtMarchingCubes( int nBlock, TRIA3DLIST& lstTria, TriHolder& triHold) abs( LocCompV[1].ptInt.z - ptBarycenter.z) < EPS_SMALL && abs( LocCompV[2].ptInt.z - ptBarycenter.z) < EPS_SMALL && abs( LocCompV[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 ; } @@ -1161,14 +1218,8 @@ VolZmap::ExtMarchingCubes( int nBlock, TRIA3DLIST& lstTria, TriHolder& triHold) VectorField LocVecF[12], LocCompV[12] ; for ( int LocInd = 0 ; LocInd < 12 ; ++ LocInd) { - LocVecF[LocInd] = VecField[LocInd] ; LocCompV[LocInd] = CompoVert[nCompCount - 1][LocInd] ; - - LocVecF[LocInd].ptInt.ToLoc( m_MapFrame[0]) ; - LocVecF[LocInd].vtNorm.ToLoc( m_MapFrame[0]) ; - LocCompV[LocInd].ptInt.ToLoc( m_MapFrame[0]) ; - LocCompV[LocInd].vtNorm.ToLoc( m_MapFrame[0]) ; } if ( ( AreSameVectorApprox( LocVecF[nIndArrey[nCompCount - 1][0]].vtNorm, LocVecF[nIndArrey[nCompCount - 1][2]].vtNorm) && @@ -1185,17 +1236,17 @@ VolZmap::ExtMarchingCubes( int nBlock, TRIA3DLIST& lstTria, TriHolder& triHold) abs( LocCompV[1].ptInt.x - ptBarycenter.x) < EPS_SMALL && abs( LocCompV[2].ptInt.x - ptBarycenter.x) < EPS_SMALL && abs( LocCompV[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 ; } @@ -1206,17 +1257,17 @@ VolZmap::ExtMarchingCubes( int nBlock, TRIA3DLIST& lstTria, TriHolder& triHold) abs( LocCompV[1].ptInt.y - ptBarycenter.y) < EPS_SMALL && abs( LocCompV[2].ptInt.y - ptBarycenter.y) < EPS_SMALL && abs( LocCompV[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 ; } @@ -1227,34 +1278,34 @@ VolZmap::ExtMarchingCubes( int nBlock, TRIA3DLIST& lstTria, TriHolder& triHold) abs( LocCompV[1].ptInt.z - ptBarycenter.z) < EPS_SMALL && abs( LocCompV[2].ptInt.z - ptBarycenter.z) < EPS_SMALL && abs( LocCompV[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 ; } } - } + } } } } - + // Flag ExtMC bool bExtMC = ( nFeatureType != NO_FEATURE) && bGridControl ; // Extended MC - if ( bExtMC) { - + if ( bExtMC) { + // Passo al sistema di riferimento del baricentro Point3d ptGravityCenter( 0, 0, 0) ; for ( int ni = 0 ; ni < nVertComp[nCompCount - 1] ; ++ ni) @@ -1264,7 +1315,7 @@ VolZmap::ExtMarchingCubes( int nBlock, TRIA3DLIST& lstTria, TriHolder& triHold) Vector3d vtTrasf[12] ; for ( int ni = 0 ; ni < nVertComp[nCompCount - 1] ; ++ ni) vtTrasf[ni] = CompoVert[nCompCount - 1][ni].ptInt - ptGravityCenter ; - + // Preparo le matrici per il sistema typedef Eigen::Matrix dSystemMatrix ; typedef Eigen::Matrix dSystemVector ; @@ -1273,7 +1324,7 @@ VolZmap::ExtMarchingCubes( int nBlock, TRIA3DLIST& lstTria, TriHolder& triHold) dSystemMatrix dMatrixN( nVertComp[nCompCount - 1], 3) ; dSystemVector dKnownVector( nVertComp[nCompCount - 1], 1) ; -#if 0 +#if 0 // Studio del caso 4 punti su un piano int nEqual = 0 ; int nPosD ; @@ -1281,23 +1332,23 @@ VolZmap::ExtMarchingCubes( int nBlock, TRIA3DLIST& lstTria, TriHolder& triHold) if ( nVertComp[nCompCount - 1] == 4 && nFeatureType == EDGE) { int nPosEq ; - for ( int ni = 0 ; ni < 2 ; ++ ni) { + for ( int ni = 0 ; ni < 2 ; ++ ni) { for ( int nj = ni + 1 ; nj < nVertComp[nCompCount - 1] ; ++ nj) { if ( AreSameVectorApprox( CompoVert[nCompCount - 1][ni].vtNorm, CompoVert[nCompCount - 1][nj].vtNorm)) { nEqual ++ ; nPosEq = ni ; - } + } } - if ( nEqual == 2) - break ; + if ( nEqual == 2) + break ; else nEqual = 0 ; - } + } if ( nEqual == 2) { - - for ( int ni = 0 ; ni < nVertComp[nCompCount - 1] ; ++ ni) + + for ( int ni = 0 ; ni < nVertComp[nCompCount - 1] ; ++ ni) if ( ! AreSameVectorApprox( CompoVert[nCompCount - 1][ni].vtNorm, CompoVert[nCompCount - 1][nPosEq].vtNorm)) { @@ -1307,14 +1358,14 @@ VolZmap::ExtMarchingCubes( int nBlock, TRIA3DLIST& lstTria, TriHolder& triHold) } } } - + double dDot = abs( ( CompoVert[nCompCount - 1][1].ptInt - CompoVert[nCompCount - 1][0].ptInt) * ( ( CompoVert[nCompCount - 1][2].ptInt - CompoVert[nCompCount - 1][1].ptInt) ^ ( CompoVert[nCompCount - 1][3].ptInt - CompoVert[nCompCount - 1][2].ptInt))) ; - + // Caso superficie piana if ( false && nVertComp[nCompCount - 1] == 4 && nEqual == 2 && dDot < EPS_SMALL) { - + for ( int ni = 0 ; ni < nVertComp[nCompCount - 1] ; ++ ni) { if ( ni != nPosD) { dMatrixN( ni, 0) = CompoVert[nCompCount - 1][ni].vtNorm.x ; @@ -1326,7 +1377,7 @@ VolZmap::ExtMarchingCubes( int nBlock, TRIA3DLIST& lstTria, TriHolder& triHold) dMatrixN( ni, 0) = vtE.x ; dMatrixN( ni, 1) = vtE.y ; dMatrixN( ni, 2) = vtE.z ; - dKnownVector( ni) = vtE *vtTrasf[ni] ; + dKnownVector( ni) = vtE *vtTrasf[ni] ; } } } @@ -1337,7 +1388,7 @@ VolZmap::ExtMarchingCubes( int nBlock, TRIA3DLIST& lstTria, TriHolder& triHold) dMatrixN( ni, 1) = CompoVert[nCompCount - 1][ni].vtNorm.y ; dMatrixN( ni, 2) = CompoVert[nCompCount - 1][ni].vtNorm.z ; dKnownVector( ni) = CompoVert[nCompCount - 1][ni].vtNorm * vtTrasf[ni] ; - } + } } #else // Caso generale @@ -1346,9 +1397,9 @@ VolZmap::ExtMarchingCubes( int nBlock, TRIA3DLIST& lstTria, TriHolder& triHold) dMatrixN( ni, 1) = CompoVert[nCompCount - 1][ni].vtNorm.y ; dMatrixN( ni, 2) = CompoVert[nCompCount - 1][ni].vtNorm.z ; dKnownVector( ni) = CompoVert[nCompCount - 1][ni].vtNorm * vtTrasf[ni] ; - } + } #endif - + // calcolo SVD DecomposerSVD svd( dMatrixN, Eigen::ComputeThinU | Eigen::ComputeThinV) ; dSystemMatrix dMatrixV = svd.matrixV() ; @@ -1358,42 +1409,42 @@ VolZmap::ExtMarchingCubes( int nBlock, TRIA3DLIST& lstTria, TriHolder& triHold) if ( nFeatureType == EDGE) { double dThres = 0.5 * ( dSingularValue( 1) + dSingularValue( 2)) / dSingularValue( 0) ; svd.setThreshold( dThres) ; - } + } // risolvo il sistema con SVD, quindi ai minimi quadrati dSystemVector dUnknownVector( 3, 1) ; dUnknownVector = svd.solve( dKnownVector) ; - // Vettore Baricentro-Feature + // Vettore Baricentro-Feature Vector3d vtFeature( dUnknownVector( 0), dUnknownVector( 1), dUnknownVector( 2)) ; - - // Esprimo la soluzione nel sistema di riferimento in cui è immerso quello dello z-map + + // Esprimo la soluzione nel sistema di riferimento z-map Point3d ptSol = ptGravityCenter + vtFeature ; // Limito la feature a una distanza dal baricentro pari alla diagonale del voxel double dDistFeature = vtFeature.Len() ; - const double MAX_DIST = sqrt( 3) * m_dStep ; - bool bOutside = ( dDistFeature > MAX_DIST) ; - + const double MAX_DIST = sqrt( 3) * m_dStep ; + bool bOutside = ( dDistFeature > MAX_DIST) ; + TRIA3DVECTOR triContainer ; - - // Costruisco triangoli di prova + + // Costruisco triangoli di prova for ( int ni = 0 ; ni < nVertComp[nCompCount - 1] ; ++ ni) { - int nj = ( ni + 1 < nVertComp[nCompCount - 1]) ? ni + 1 : 0 ; + int nj = ( ni + 1 < nVertComp[nCompCount - 1]) ? ni + 1 : 0 ; // Il triangolo è pronto Triangle3d CurrentTriangle ; CurrentTriangle.Set( ptSol, CompoVert[nCompCount - 1][nj].ptInt, CompoVert[nCompCount - 1][ni].ptInt) ; CurrentTriangle.Validate( true) ; // Aggiungo triangolo al vettore temporaneo - triContainer.emplace_back( CurrentTriangle) ; + triContainer.emplace_back( CurrentTriangle) ; } - + // Controllo delle inversioni dei triangoli bool bDangerInversion = false ; - + // Caso ventaglio con tre vertici di base if ( nVertComp[nCompCount - 1] == 3) { - + // Controllo se esiste almeno un triangolo con normale avente prodotto scalare // negativo con la normale di almeno uno dei vertici di base del ventaglio. bool bInversione = false ; @@ -1403,8 +1454,8 @@ VolZmap::ExtMarchingCubes( int nBlock, TRIA3DLIST& lstTria, TriHolder& triHold) double dDI = triContainer[ni].GetN() * CompoVert[nCompCount - 1][ni].vtNorm ; double dDJ = triContainer[ni].GetN() * CompoVert[nCompCount - 1][nj].vtNorm ; - - if ( dDI < - EPS_SMALL || dDJ < - EPS_SMALL) { + + if ( dDI < - EPS_SMALL || dDJ < - EPS_SMALL) { bInversione = true ; break ; } @@ -1412,48 +1463,42 @@ VolZmap::ExtMarchingCubes( int nBlock, TRIA3DLIST& lstTria, TriHolder& triHold) // Se tale triangolo esiste procedo if ( bInversione) { - + // Conto le coppie di normali con angolo compreso maggiore di 90° int nNegDot = 0 ; for ( int ni = 0 ; ni < nVertComp[nCompCount - 1] - 1 ; ++ ni) { - for ( int nj = ni + 1 ; nj < nVertComp[nCompCount - 1] ; ++ nj) { + for ( int nj = ni + 1 ; nj < nVertComp[nCompCount - 1] ; ++ nj) { if ( CompoVert[nCompCount - 1][ni].vtNorm * CompoVert[nCompCount - 1][nj].vtNorm < - EPS_SMALL) nNegDot ++ ; } } - + if ( nNegDot == nVertComp[nCompCount - 1] - 1) { - // Cerco se esiste un punto in cui la normale ha prodotto scalre negativo // con le normali di entrambi i triangoli che lo contengono bool bInversione2 = false ; for ( int ni = 0 ; ni < nVertComp[nCompCount - 1] ; ++ ni) { - int nj = ( ni == 0 ? nVertComp[nCompCount - 1] - 1 : ni - 1) ; - double dDLast = triContainer[ni].GetN() * CompoVert[nCompCount - 1][ni].vtNorm ; double dDPrev = triContainer[nj].GetN() * CompoVert[nCompCount - 1][ni].vtNorm ; - if ( ( dDLast < EPS_SMALL && dDPrev < EPS_SMALL) || ( dDLast < - 0.75 || dDPrev < - 0.75)) { bInversione2 = true ; break ; } } - + // Se tale normale esiste if ( bInversione2) { - // Soluzione del sistema nel sistema Zmap + // Soluzione del sistema nel sistema Zmap Point3d ptSolZMapFrame = ptSol ; - ptSolZMapFrame.ToLoc( m_MapFrame[0]) ; - + // Se la soluzione non cade nel voxel di appartenenza vedo se può // essere riportata dentro muovendosi lungo la linea di feature. if ( ! IsPointInsideVoxelApprox( i, j, k, ptSolZMapFrame)) { - + Vector3d vtNullSpace( dMatrixV( 0, 2), dMatrixV( 1, 2), dMatrixV( 2, 2)) ; - vtNullSpace.ToLoc( m_MapFrame[0]) ; double dParInt1, dParInt2 ; Point3d ptVoxMin( ( i + 0.5) * m_dStep, ( j + 0.5) * m_dStep, ( k + 0.5) * m_dStep) ; Point3d ptVoxMax( ( i + 1.5) * m_dStep, ( j + 1.5) * m_dStep, ( k + 1.5) * m_dStep) ; @@ -1467,20 +1512,18 @@ VolZmap::ExtMarchingCubes( int nBlock, TRIA3DLIST& lstTria, TriHolder& triHold) dParInt2 + ( dParInt1 - dParInt2) / 100 ; Point3d ptNewSol = ptSolZMapFrame + dPar * vtNullSpace ; - - ptNewSol.ToGlob( m_MapFrame[0]) ; ptSol = ptNewSol ; - // Costruisco triangoli di prova + // Costruisco triangoli di prova for ( int ni = 0 ; ni < nVertComp[nCompCount - 1] ; ++ ni) { - int nj = ( ni + 1 < nVertComp[nCompCount - 1]) ? ni + 1 : 0 ; + int nj = ( ni + 1 < nVertComp[nCompCount - 1]) ? ni + 1 : 0 ; // Il triangolo è pronto Triangle3d CurrentTriangle ; CurrentTriangle.Set( ptSol, CompoVert[nCompCount - 1][nj].ptInt, CompoVert[nCompCount - 1][ni].ptInt) ; CurrentTriangle.Validate( true) ; // Aggiungo triangolo al vettore temporaneo - triContainer.emplace_back( CurrentTriangle) ; + triContainer.emplace_back( CurrentTriangle) ; } } @@ -1488,53 +1531,37 @@ VolZmap::ExtMarchingCubes( int nBlock, TRIA3DLIST& lstTria, TriHolder& triHold) // si passa alla routine standard se il voxel // in cui cade è pieno. else { - + int nAdjVoxI, nAdjVoxJ, nAdjVoxK ; if ( GetPointVoxel( ptSolZMapFrame, nAdjVoxI, nAdjVoxJ, nAdjVoxK)) { // Classificazione del voxel adiacente - int nAdjIndex = 0 ; - if ( IsThereMat( nAdjVoxI, nAdjVoxJ, nAdjVoxK)) - nAdjIndex |= ( 1 << 0) ; - if ( IsThereMat( nAdjVoxI + 1, nAdjVoxJ, nAdjVoxK)) - nAdjIndex |= ( 1 << 1) ; - if ( IsThereMat( nAdjVoxI + 1, nAdjVoxJ + 1, nAdjVoxK)) - nAdjIndex |= ( 1 << 2) ; - if ( IsThereMat( nAdjVoxI, nAdjVoxJ + 1, nAdjVoxK)) - nAdjIndex |= ( 1 << 3) ; - if ( IsThereMat( nAdjVoxI, nAdjVoxJ, nAdjVoxK + 1)) - nAdjIndex |= ( 1 << 4) ; - if ( IsThereMat( nAdjVoxI + 1, nAdjVoxJ, nAdjVoxK + 1)) - nAdjIndex |= ( 1 << 5) ; - if ( IsThereMat( nAdjVoxI + 1, nAdjVoxJ + 1, nAdjVoxK + 1)) - nAdjIndex |= ( 1 << 6) ; - if ( IsThereMat( nAdjVoxI, nAdjVoxJ + 1, nAdjVoxK + 1)) - nAdjIndex |= ( 1 << 7) ; - + int nAdjIndex = CalcIndex( nAdjVoxI, nAdjVoxJ, nAdjVoxK) ; + // Se il voxel è pieno if ( EdgeTable[nAdjIndex] != 0) - bDangerInversion = true ; + bDangerInversion = true ; } - } - } + } + } } } } } - + // Ventaglio con base a quattro vertici else if ( nVertComp[nCompCount - 1] == 4) { - + // Controllo se esiste almeno un triangolo con normale avente prodotto scalare // negativo con la normale di almeno uno dei vertici di base del ventaglio. bool bInversione = false ; for ( int ni = 0 ; ni < nVertComp[nCompCount - 1] ; ++ ni) { - int nj = ( ni + 1 < nVertComp[nCompCount - 1]) ? ni + 1 : 0 ; + int nj = ( ni + 1 < nVertComp[nCompCount - 1]) ? ni + 1 : 0 ; double dDI = triContainer[ni].GetN() * CompoVert[nCompCount - 1][ni].vtNorm ; double dDJ = triContainer[ni].GetN() * CompoVert[nCompCount - 1][nj].vtNorm ; - + if ( dDI < - EPS_SMALL || dDJ < - EPS_SMALL) bInversione = true ; } @@ -1553,15 +1580,14 @@ VolZmap::ExtMarchingCubes( int nBlock, TRIA3DLIST& lstTria, TriHolder& triHold) // Caso in cui tale numero è 3 if ( nNegDot == nVertComp[nCompCount - 1] - 1) { - + Point3d ptSolZMapFrame = ptSol ; - ptSolZMapFrame.ToLoc( m_MapFrame[0]) ; - + // Se la feature non cade nel suo voxel if ( ! IsPointInsideVoxelApprox( i, j, k, ptSolZMapFrame)) { Vector3d vtNullSpace( dMatrixV( 0, 2), dMatrixV( 1, 2), dMatrixV( 2, 2)) ; - vtNullSpace.ToLoc( m_MapFrame[0]) ; + double dParInt1, dParInt2 ; Point3d ptVoxMin( ( i + 0.5) * m_dStep, ( j + 0.5) * m_dStep, ( k + 0.5) * m_dStep) ; Point3d ptVoxMax( ( i + 1.5) * m_dStep, ( j + 1.5) * m_dStep, ( k + 1.5) * m_dStep) ; @@ -1573,22 +1599,20 @@ VolZmap::ExtMarchingCubes( int nBlock, TRIA3DLIST& lstTria, TriHolder& triHold) double dPar = abs( dParInt1) < abs( dParInt2) ? dParInt1 + ( dParInt2 - dParInt1) / 100 : dParInt2 + ( dParInt1 - dParInt2) / 100 ; - + Point3d ptNewSol = ptSolZMapFrame + dPar * vtNullSpace ; - - ptNewSol.ToGlob( m_MapFrame[0]) ; ptSol = ptNewSol ; - // Costruisco triangoli di prova + // Costruisco triangoli di prova for ( int ni = 0 ; ni < nVertComp[nCompCount - 1] ; ++ ni) { - int nj = ( ni + 1 < nVertComp[nCompCount - 1]) ? ni + 1 : 0 ; + int nj = ( ni + 1 < nVertComp[nCompCount - 1]) ? ni + 1 : 0 ; // Il triangolo è pronto Triangle3d CurrentTriangle ; CurrentTriangle.Set( ptSol, CompoVert[nCompCount - 1][nj].ptInt, CompoVert[nCompCount - 1][ni].ptInt) ; CurrentTriangle.Validate( true) ; // Aggiungo triangolo al vettore temporaneo - triContainer.emplace_back( CurrentTriangle) ; + triContainer.emplace_back( CurrentTriangle) ; } } @@ -1598,67 +1622,48 @@ VolZmap::ExtMarchingCubes( int nBlock, TRIA3DLIST& lstTria, TriHolder& triHold) int nAdjVoxI, nAdjVoxJ, nAdjVoxK ; if ( GetPointVoxel( ptSolZMapFrame, nAdjVoxI, nAdjVoxJ, nAdjVoxK)) { // Classificazione del voxel adiacente - int nAdjIndex = 0 ; - if ( IsThereMat( nAdjVoxI, nAdjVoxJ, nAdjVoxK)) - nAdjIndex |= ( 1 << 0) ; - if ( IsThereMat( nAdjVoxI + 1, nAdjVoxJ, nAdjVoxK)) - nAdjIndex |= ( 1 << 1) ; - if ( IsThereMat( nAdjVoxI + 1, nAdjVoxJ + 1, nAdjVoxK)) - nAdjIndex |= ( 1 << 2) ; - if ( IsThereMat( nAdjVoxI, nAdjVoxJ + 1, nAdjVoxK)) - nAdjIndex |= ( 1 << 3) ; - if ( IsThereMat( nAdjVoxI, nAdjVoxJ, nAdjVoxK + 1)) - nAdjIndex |= ( 1 << 4) ; - if ( IsThereMat( nAdjVoxI + 1, nAdjVoxJ, nAdjVoxK + 1)) - nAdjIndex |= ( 1 << 5) ; - if ( IsThereMat( nAdjVoxI + 1, nAdjVoxJ + 1, nAdjVoxK + 1)) - nAdjIndex |= ( 1 << 6) ; - if ( IsThereMat( nAdjVoxI, nAdjVoxJ + 1, nAdjVoxK + 1)) - nAdjIndex |= ( 1 << 7) ; + int nAdjIndex = CalcIndex( nAdjVoxI, nAdjVoxJ, nAdjVoxK) ; // Se il voxel è pieno if ( EdgeTable[nAdjIndex] != 0) - bDangerInversion = true ; + bDangerInversion = true ; } - } - } + } + } } // Caso in cui il numero di coppie di normali con prodotto // scalare negativo non è 3 else { Point3d ptSolZMapFrame = ptSol ; - ptSolZMapFrame.ToLoc( m_MapFrame[0]) ; if ( ! IsPointInsideVoxelApprox( i, j, k, ptSolZMapFrame)) { Vector3d vtNullSpace( dMatrixV( 0, 2), dMatrixV( 1, 2), dMatrixV( 2, 2)) ; - vtNullSpace.ToLoc( m_MapFrame[0]) ; + double dParInt1, dParInt2 ; Point3d ptVoxMin( ( i + 0.5) * m_dStep, ( j + 0.5) * m_dStep, ( k + 0.5) * m_dStep) ; Point3d ptVoxMax( ( i + 1.5) * m_dStep, ( j + 1.5) * m_dStep, ( k + 1.5) * m_dStep) ; if ( IntersLineBox( ptSolZMapFrame, vtNullSpace, ptVoxMin, ptVoxMax, dParInt1, dParInt2)) { - + triContainer.resize( 0) ; double dPar = abs( dParInt1) < abs( dParInt2) ? dParInt1 + ( dParInt2 - dParInt1) / 100: dParInt2 + ( dParInt1 - dParInt2) / 100 ; Point3d ptNewSol = ptSolZMapFrame + dPar * vtNullSpace ; - - ptNewSol.ToGlob( m_MapFrame[0]) ; ptSol = ptNewSol ; - // Costruisco triangoli di prova + // Costruisco triangoli di prova for ( int ni = 0 ; ni < nVertComp[nCompCount - 1] ; ++ ni) { - int nj = ( ni + 1 < nVertComp[nCompCount - 1]) ? ni + 1 : 0 ; + int nj = ( ni + 1 < nVertComp[nCompCount - 1]) ? ni + 1 : 0 ; // Il triangolo è pronto Triangle3d CurrentTriangle ; CurrentTriangle.Set( ptSol, CompoVert[nCompCount - 1][nj].ptInt, CompoVert[nCompCount - 1][ni].ptInt) ; CurrentTriangle.Validate( true) ; // Aggiungo triangolo al vettore temporaneo - triContainer.emplace_back( CurrentTriangle) ; + triContainer.emplace_back( CurrentTriangle) ; } } @@ -1672,15 +1677,15 @@ VolZmap::ExtMarchingCubes( int nBlock, TRIA3DLIST& lstTria, TriHolder& triHold) for ( int nNI = 0 ; nNI < 3 ; ++ nNI) { for ( int nNJ = nNI + 1 ; nNJ < 4 ; ++ nNJ) { if ( AreSameVectorApprox( CompoVert[nCompCount - 1][nNI].vtNorm, - CompoVert[nCompCount - 1][nNJ].vtNorm)) { + CompoVert[nCompCount - 1][nNJ].vtNorm)) { ++ nCouple ; - if ( nCouple == 1) { + if ( nCouple == 1) { nCoupleIndex[0] = nNI ; nCoupleIndex[1] = nNJ ; } else if ( nCouple == 2) { nCoupleIndex[2] = nNI ; - nCoupleIndex[3] = nNJ ; + nCoupleIndex[3] = nNJ ; } } } @@ -1754,10 +1759,10 @@ VolZmap::ExtMarchingCubes( int nBlock, TRIA3DLIST& lstTria, TriHolder& triHold) } } } - } - } + } + } } - + // Valuto normali: questo è ancora un controllo // sulle normali, se risultano in tutti i punti // approssimativamente uguali passiamo alla @@ -1776,12 +1781,12 @@ VolZmap::ExtMarchingCubes( int nBlock, TRIA3DLIST& lstTria, TriHolder& triHold) bPlane = false ; break ; } - } + } if ( ! bPlane) break ; } - + // Se la feature non è fuori dai limiti // e i triangoli formano una superficie // non piana confermo ExtMC @@ -1792,35 +1797,37 @@ VolZmap::ExtMarchingCubes( int nBlock, TRIA3DLIST& lstTria, TriHolder& triHold) 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 * EPS_SMALL) { + if ( triContainer[ni].GetArea() > SQ_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]) ; - } - } + } + } + + // Porto la soluzione nel sistema in cui è immerso lo Zmap + ptSol.ToGlob( m_MapFrame[0]) ; // 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, @@ -1831,28 +1838,31 @@ VolZmap::ExtMarchingCubes( int nBlock, TRIA3DLIST& lstTria, TriHolder& triHold) 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 ; + 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 + // 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) + for ( int ni = 0 ; ni < int( vInnerTriaTemp.size()) ; ++ ni) { + // Riporto le coordinate nel sistema in cui è immerso lo Zmap + vInnerTriaTemp[ni].ToGlob( m_MapFrame[0]) ; triHold[nCurrent].vCompoTria[nOldFeatureNum].emplace_back( vInnerTriaTemp[ni]) ; + } } // Triangoli di frontiera if ( vBorderTriaTemp.size() > 0) { - + // Aggiorno il numero di feature. ++ nBorderFeatureInVoxel ; @@ -1860,7 +1870,7 @@ VolZmap::ExtMarchingCubes( int nBlock, TRIA3DLIST& lstTria, TriHolder& triHold) // voxel, ridimensiono il vettore che contiene // la struttura dati del voxel. if ( nBorderFeatureInVoxel == 1) { - + m_InterBlockTria[nBlock].resize( m_InterBlockTria[nBlock].size() + 1) ; // Questi dati dipendono solo dal voxel, @@ -1871,21 +1881,24 @@ VolZmap::ExtMarchingCubes( int nBlock, TRIA3DLIST& lstTria, TriHolder& triHold) m_InterBlockTria[nBlock][nCurrent].i = i ; m_InterBlockTria[nBlock][nCurrent].j = j ; m_InterBlockTria[nBlock][nCurrent].k = k ; - } + } // Indice che identifica la posizione del voxel nel vector - int nCurrent = int( m_InterBlockTria[nBlock].size()) - 1 ; + int nCurrent = int( m_InterBlockTria[nBlock].size()) - 1 ; // Aggiungo vertice della componente connessa alla lista dei vertici - m_InterBlockTria[nBlock][nCurrent].ptCompoVert.emplace_back( ptSol) ; + m_InterBlockTria[nBlock][nCurrent].ptCompoVert.emplace_back( ptSol) ; int nOldFeatureNum = int( m_InterBlockTria[nBlock][nCurrent].vCompoTria.size()) ; int nNewFeatureNum = nOldFeatureNum + 1 ; // Aggiungo una componente per il vettore dei triangoli della componente connessa m_InterBlockTria[nBlock][nCurrent].vCompoTria.resize( nNewFeatureNum) ; - for ( int ni = 0 ; ni < int( vBorderTriaTemp.size()) ; ++ ni) - m_InterBlockTria[nBlock][nCurrent].vCompoTria[nOldFeatureNum].emplace_back( vBorderTriaTemp[ni]) ; + for ( int ni = 0 ; ni < int( vBorderTriaTemp.size()) ; ++ ni) { + // Riporto le coordinate nel sistema in cui è immerso lo Zmap + vBorderTriaTemp[ni].ToGlob( m_MapFrame[0]) ; + m_InterBlockTria[nBlock][nCurrent].vCompoTria[nOldFeatureNum].emplace_back( vBorderTriaTemp[ni]) ; + } } } @@ -1893,43 +1906,55 @@ VolZmap::ExtMarchingCubes( int nBlock, TRIA3DLIST& lstTria, TriHolder& triHold) else { // Costruzione dei triangoli for ( int TriIndex = 0; TriIndex < ( nVertComp[nCompCount - 1] - 2) * 3 ; TriIndex += 3) { - // Il triangolo è pronto + // Il triangolo è pronto Triangle3d CurrentTriangle ; CurrentTriangle.Set( CompoTriVert[nCompCount - 1][TriIndex].ptInt, CompoTriVert[nCompCount - 1][TriIndex+1].ptInt, CompoTriVert[nCompCount - 1][TriIndex+2].ptInt) ; bool bV = CurrentTriangle.Validate( true) ; - // Aggiungo alla lista - lstTria.emplace_back( CurrentTriangle) ; - } - } + // Riporto lew coordinate nel sistema in cui è immerso lo Zmap + CurrentTriangle.ToGlob( m_MapFrame[0]) ; + // Aggiungo alla lista + lstTria.emplace_back( CurrentTriangle) ; + } + } } // Standard MC - else { + else { // Costruzione dei triangoli for ( int TriIndex = 0; TriIndex < ( nVertComp[nCompCount - 1] - 2) * 3 ; TriIndex += 3) { - // Il triangolo è pronto + // Il triangolo è pronto Triangle3d CurrentTriangle ; CurrentTriangle.Set( CompoTriVert[nCompCount - 1][TriIndex].ptInt, CompoTriVert[nCompCount - 1][TriIndex+1].ptInt, CompoTriVert[nCompCount - 1][TriIndex+2].ptInt) ; bool bV = CurrentTriangle.Validate( true) ; - // Aggiungo alla lista - lstTria.emplace_back( CurrentTriangle) ; - } - } - } + // Riporto lew coordinate nel sistema in cui è immerso lo Zmap + CurrentTriangle.ToGlob( m_MapFrame[0]) ; + // Aggiungo alla lista + lstTria.emplace_back( CurrentTriangle) ; + } + } + } } } } + // Processo i Voxel con possibile superficie piana + ProcessVoxContXY( VoxContXYInf, false, lstTria) ; + ProcessVoxContXY( VoxContXYSup, true, lstTria) ; + ProcessVoxContYZ( VoxContYZInf, false, lstTria) ; + ProcessVoxContYZ( VoxContYZSup, true, lstTria) ; + ProcessVoxContXZ( VoxContXZInf, false, lstTria) ; + ProcessVoxContXZ( VoxContXZSup, true, lstTria) ; + return true ; } //---------------------------------------------------------------------------- -bool -VolZmap::FlipEdgesII( TriHolder& TriHold) const +bool +VolZmap::FlipEdgesII( TriHolder& TriHold) const { // Numero di voxel in cui si presentano sharp feature int nVoxelNum = int( TriHold.size()) ; @@ -1960,21 +1985,21 @@ VolZmap::FlipEdgesII( TriHolder& TriHold) const // 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) ; + Point3d ptP2 = TriHold[n2].vCompoTria[nCompo2][nTri2].GetP( nVert2) ; if ( AreSamePointEpsilon( ptP1, ptP2, EPS_ZERO)) { @@ -1983,7 +2008,7 @@ VolZmap::FlipEdgesII( TriHolder& TriHold) const if ( ! ( AreSamePointEpsilon( ptP1, ptVert1, EPS_ZERO) || AreSamePointEpsilon( ptP2, ptVert2, EPS_ZERO))) { - + SharedIndex.emplace_back( nVert1) ; SharedIndex.emplace_back( nVert2) ; } @@ -2000,33 +2025,33 @@ VolZmap::FlipEdgesII( TriHolder& TriHold) const // Si deve operare la modifica dei triangoli if ( SharedIndex.size() > 2) { - // verifico che i due lati adiacenti siano controversi + // verifico che i due lati adiacenti siano controversi int nProd = ( SharedIndex[2] - SharedIndex[0]) * ( SharedIndex[3] - SharedIndex[1]) ; if ( nProd != 1 && nProd != - 2 && nProd != 4 && nProd != 0) { - + 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 ; - } - } + break ; + } + } } - if ( bModified) + if ( bModified) break ; } } } } - } + } } - + return true ; } @@ -2040,11 +2065,11 @@ VolZmap::FlipEdgesBB( TriaMatrix& InterTria) const // ciclo sui blocchi for ( size_t tFB = 0 ; tFB < nBlocksNum - 1 ; ++ tFB) { - int nFBijk[3] ; + int nFBijk[3] ; GetBlockIJKFromN( int( tFB), nFBijk) ; for ( size_t tLB = tFB + 1 ; tLB < nBlocksNum ; ++ tLB) { - + int nLBijk[3] ; GetBlockIJKFromN( int( tLB), nLBijk) ; @@ -2070,30 +2095,30 @@ VolZmap::FlipEdgesBB( TriaMatrix& InterTria) const continue ; // Numero di componenti connesse dei voxel - size_t nCompoVFBNum = InterTria[tFB][tVFB].ptCompoVert.size() ; - size_t nCompoVLBNum = InterTria[tLB][tVLB].ptCompoVert.size() ; + 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() ; + 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) ; + Point3d ptPL = InterTria[tLB][tVLB].vCompoTria[tCmpL][tTriLB].GetP( nVertL) ; if ( AreSamePointEpsilon( ptPF, ptPL, EPS_ZERO)) { @@ -2102,24 +2127,24 @@ VolZmap::FlipEdgesBB( TriaMatrix& InterTria) const if ( ! ( AreSamePointEpsilon( ptPF, ptVertF, EPS_ZERO) || AreSamePointEpsilon( ptPL, ptVertL, EPS_ZERO))) { - + SharedIndex.emplace_back( nVertF) ; SharedIndex.emplace_back( nVertL) ; } } - if ( SharedIndex.size() > 2) + if ( SharedIndex.size() > 2) break ; } - if ( SharedIndex.size() > 2) + if ( SharedIndex.size() > 2) break ; } // Si deve operare la modifica dei triangoli if ( SharedIndex.size() > 2) { - - // verifico che i due lai adiacenti siano controversi + + // verifico che i due lai adiacenti siano controversi int nProd = ( SharedIndex[2] - SharedIndex[0]) * ( SharedIndex[3] - SharedIndex[1]) ; if ( nProd != 1 && nProd != - 2 && nProd != 4 && nProd != 0) { @@ -2127,16 +2152,16 @@ VolZmap::FlipEdgesBB( TriaMatrix& InterTria) const 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) + if ( bModified) break ; } } @@ -2157,42 +2182,36 @@ VolZmap::IsThereMat( int nI, int nJ, int nK) const 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 ; - + // ciclo sulle griglie int nCount = 0 ; - for ( int nGrid = 0 ; nGrid < int ( m_nMapNum) ; ++ nGrid) { - + // assegnazione dati vertice dipendenti dalla griglia unsigned int nGrI, nGrJ ; - - if ( nGrid == 0) { + double dZ ; + switch ( nGrid) { + case 0 : nGrI = nI ; nGrJ = nJ ; - } - else if ( nGrid == 1) { + dZ = ( nK + 0.5) * m_dStep ; + break ; + case 1 : nGrI = nJ ; nGrJ = nK ; - } - else { + dZ = ( nI + 0.5) * m_dStep ; + break ; + case 2 : nGrI = nK ; nGrJ = nI ; + dZ = ( nJ + 0.5) * m_dStep ; + break ; } - - unsigned int nPos = nGrJ * m_nNx[nGrid] + nGrI ; - size_t nDexSize = m_Values[nGrid][nPos].size() ; + // verifica spillone su vertice size_t nIndex = 0 ; - + unsigned int nPos = nGrJ * m_nNx[nGrid] + nGrI ; + size_t nDexSize = m_Values[nGrid][nPos].size() ; while ( nIndex < nDexSize) { - - if ( dZ[nGrid] > m_Values[nGrid][nPos][nIndex].dMin - dEps && - dZ[nGrid] < m_Values[nGrid][nPos][nIndex].dMax + dEps) { - + if ( dZ > m_Values[nGrid][nPos][nIndex].dMin - 2 * EPS_SMALL && + dZ < m_Values[nGrid][nPos][nIndex].dMax + 2 * EPS_SMALL) { ++ nCount ; break ; } @@ -2202,6 +2221,30 @@ VolZmap::IsThereMat( int nI, int nJ, int nK) const return ( nCount == 3) ; } +//---------------------------------------------------------------------------- +int +VolZmap::CalcIndex( int nI, int nJ, int nK) const +{ + int nIndex = 0 ; + if ( IsThereMat( nI, nJ, nK)) + nIndex |= ( 1 << 0) ; + if ( IsThereMat( nI + 1, nJ, nK)) + nIndex |= ( 1 << 1) ; + if ( IsThereMat( nI + 1, nJ + 1, nK)) + nIndex |= ( 1 << 2) ; + if ( IsThereMat( nI, nJ + 1, nK)) + nIndex |= ( 1 << 3) ; + if ( IsThereMat( nI, nJ, nK + 1)) + nIndex |= ( 1 << 4) ; + if ( IsThereMat( nI + 1, nJ, nK + 1)) + nIndex |= ( 1 << 5) ; + if ( IsThereMat( nI + 1, nJ + 1, nK + 1)) + nIndex |= ( 1 << 6) ; + if ( IsThereMat( nI, nJ + 1, nK + 1)) + nIndex |= ( 1 << 7) ; + return nIndex ; +} + //---------------------------------------------------------------------------- bool VolZmap::IntersPos( int nVec1[], int nVec2[], Point3d& ptInt) const @@ -2231,14 +2274,14 @@ VolZmap::IntersPos( int nVec1[], int nVec2[], Point3d& ptInt) const bFound = true ; break ; } - else if ( dx1 > dMinX - EPS_SMALL && dx1 < dMaxX + EPS_SMALL && dx2 > dMaxX + EPS_SMALL) { - ptInt.x = dx1 ; + 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 ; + ptInt.x = ( dMinX + dMaxX) / 2 ; } else if ( nVec1[1] != nVec2[1]) { @@ -2262,19 +2305,19 @@ VolZmap::IntersPos( int nVec1[], int nVec2[], Point3d& ptInt) const double dy2 = m_Values[2][nDexel][j].dMax ; if ( dy1 < dMinY - EPS_SMALL && dy2 > dMinY - EPS_SMALL && dy2 < dMaxY + EPS_SMALL) { - ptInt.y = dy2 ; + ptInt.y = dy2 ; bFound = true ; break ; } else if ( dy1 > dMinY - EPS_SMALL && dy1 < dMaxY + EPS_SMALL && dy2 > dMaxY + EPS_SMALL) { - ptInt.y = dy1 ; + ptInt.y = dy1 ; bFound = true ; break ; } - } + } if ( ! bFound) ptInt.y = ( dMinY + dMaxY) / 2 ; - } + } else if ( nVec1[2] != nVec2[2]) { @@ -2292,21 +2335,21 @@ VolZmap::IntersPos( int nVec1[], int nVec2[], Point3d& ptInt) const bool bFound = false ; for ( size_t k = 0 ; k < nSize ; k += 1) { - + double dz1 = m_Values[0][nDexel][k].dMin ; double dz2 = m_Values[0][nDexel][k].dMax ; if ( dz1 < dMinZ - EPS_SMALL && dz2 > dMinZ - EPS_SMALL && dz2 < dMaxZ + EPS_SMALL) { - ptInt.z = dz2 ; + ptInt.z = dz2 ; bFound = true ; break ; } else if ( dz1 > dMinZ - EPS_SMALL && dz1 < dMaxZ + EPS_SMALL && dz2 > dMaxZ + EPS_SMALL) { - ptInt.z = dz1 ; + ptInt.z = dz1 ; bFound = true ; break ; } - } + } if ( ! bFound) ptInt.z = ( dMinZ + dMaxZ) / 2 ; } @@ -2315,18 +2358,18 @@ VolZmap::IntersPos( int nVec1[], int nVec2[], Point3d& ptInt) const } //---------------------------------------------------------------------------- -bool +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 ; @@ -2344,17 +2387,17 @@ VolZmap::IntersPos( int nVec1[], int nVec2[], bool bFirstCorner, Point3d& ptInt, while ( n >= 0 && dX > dMinX - dEps) { if ( dX < dMaxX + dEps) { - + ptInt.x = dX ; vtNormal = m_Values[1][nDexel][n].vtMaxN ; bFound = true ; break ; } - + n -= 1 ; if ( n >= 0) - dX = m_Values[1][nDexel][n].dMax ; + dX = m_Values[1][nDexel][n].dMax ; } } else { @@ -2365,19 +2408,19 @@ VolZmap::IntersPos( int nVec1[], int nVec2[], bool bFirstCorner, Point3d& ptInt, while ( n < nSize && dX < dMaxX + dEps) { if ( dX > dMinX - dEps) { - + ptInt.x = dX ; vtNormal = m_Values[1][nDexel][n].vtMinN ; bFound = true ; break ; } - + n += 1 ; if ( n < nSize) - dX = m_Values[1][nDexel][n].dMin ; + dX = m_Values[1][nDexel][n].dMin ; } - } + } if ( ! bFound) @@ -2399,58 +2442,46 @@ VolZmap::IntersPos( int nVec1[], int nVec2[], bool bFirstCorner, Point3d& ptInt, size_t nSize = m_Values[2][nDexel].size() ; if ( bFirstCorner) { - int n = int( nSize) - 1 ; double dY = m_Values[2][nDexel][n].dMax ; - while ( n >= 0 && dY > dMinY - dEps) { - if ( dY < dMaxY + dEps) { - ptInt.y = dY ; vtNormal = m_Values[2][nDexel][n].vtMaxN ; bFound = true ; break ; } - n -= 1 ; - - if ( n >= 0) - dY = m_Values[2][nDexel][n].dMax ; + if ( n >= 0) + dY = m_Values[2][nDexel][n].dMax ; } } - else { + else { size_t n = 0 ; double dY = m_Values[2][nDexel][n].dMin ; - while ( n < nSize && dY < dMaxY + dEps) { - if ( dY > dMinY - dEps) { - ptInt.y = dY ; vtNormal = m_Values[2][nDexel][n].vtMinN ; bFound = true ; break ; } - n += 1 ; - if ( n < nSize) - dY = m_Values[2][nDexel][n].dMin ; + dY = m_Values[2][nDexel][n].dMin ; } - } + } 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 ; @@ -2461,56 +2492,81 @@ VolZmap::IntersPos( int nVec1[], int nVec2[], bool bFirstCorner, Point3d& ptInt, size_t nSize = m_Values[0][nDexel].size() ; if ( bFirstCorner) { - int n = int( nSize) - 1 ; double dZ = m_Values[0][nDexel][n].dMax ; - - while ( n >= 0 && dZ > dMinZ - dEps) { - + while ( n >= 0 && dZ > dMinZ - dEps) { if ( dZ < dMaxZ + dEps) { - ptInt.z = dZ ; vtNormal = m_Values[0][nDexel][n].vtMaxN ; bFound = true ; break ; } - n -= 1 ; - if ( n >= 0) - dZ = m_Values[0][nDexel][n].dMax ; + dZ = m_Values[0][nDexel][n].dMax ; } } - else { + else { size_t n = 0 ; double dZ = m_Values[0][nDexel][n].dMin ; - while ( n < nSize && dZ < dMaxZ + dEps) { - if ( dZ > dMinZ - dEps) { - ptInt.z = dZ ; vtNormal = m_Values[0][nDexel][n].vtMinN ; bFound = true ; break ; } - n += 1 ; - if ( n < nSize) - dZ = m_Values[0][nDexel][n].dMin ; + dZ = m_Values[0][nDexel][n].dMin ; } } if ( ! bFound) - ptInt.z = 0.5 * ( dMinZ + dMaxZ) ; } return bFound ; } +//---------------------------------------------------------------------------- +bool +VolZmap::GetVoxIJKFromN( int nN, int& nI, int& nJ, int& nK) const +{ + // Controllo sulla validità del voxel + if ( nN < 0 || nN >= int( ( m_nNx[0] + 1) * ( m_nNy[0] + 1) * ( m_nNy[1] + 1))) + return false ; + + // Calcolo gli indici + nI = ( nN % ( ( m_nNx[0] + 1) * ( m_nNy[0] + 1))) % ( m_nNx[0] + 1) - 1 ; + nJ = ( nN % ( ( m_nNx[0] + 1) * ( m_nNy[0] + 1))) / ( m_nNx[0] + 1) - 1 ; + nK = ( nN / ( ( m_nNx[0] + 1) * ( m_nNy[0] + 1))) - 1 ; + + // Controllo la sensatezza del risultato ottenuto + return ( nI >= - 1 && nI < int( m_nNx[0]) && + nJ >= - 1 && nJ < int( m_nNy[0]) && + nK >= - 1 && nK < int( m_nNy[1]) ) ; +} + +//---------------------------------------------------------------------------- +bool +VolZmap::GetVoxNFromIJK( int nI, int nJ, int nK, int& nN) const +{ + // Controllo sull'esistenza del voxel + if ( nI <= - 2 || nI >= int( m_nNx[0]) || + nJ <= - 2 || nJ >= int( m_nNy[0]) || + nK <= - 2 || nK >= int( m_nNy[1]) ) + return false ; + + // Calcolo il numero di Voxel + nN = ( m_nNx[0] + 1) * ( m_nNy[0] + 1) * ( nK + 1) + + ( m_nNx[0] + 1) * ( nJ + 1) + nI + 1 ; + + // controllo sulla sensatezza del numero ottenuto + return ( nN >= 0 && nN < int( ( m_nNx[0] + 1) * ( m_nNy[0] + 1) * ( m_nNy[1] + 1))) ; +} + //---------------------------------------------------------------------------- bool VolZmap::GetBlockIJKFromN( int nBlock, int nIJK[]) const @@ -2644,15 +2700,10 @@ VolZmap::GetVoxelBlockIJK( const int nVoxIJK[], int nBlockIJK[]) const 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))) ; + 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 ; } @@ -2660,7 +2711,7 @@ VolZmap::GetVoxelBlockIJK( const int nVoxIJK[], int nBlockIJK[]) const //---------------------------------------------------------------------------- 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 || @@ -2679,214 +2730,147 @@ VolZmap::GetAdjBlockToBlock( int nBlockN, int nDeltaI, int nDeltaJ, int nDeltaK, //---------------------------------------------------------------------------- bool -VolZmap::IsAVoxelOnBoundary( const int nLimits[], const int nIJK[], bool bType) const +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])) + 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])) + 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, + // Cerchiamo i voxel che confinano con quelli di altri 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) + 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. + int nAdjBlockIJK[3] ; 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])) { - + 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. + // 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) + 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 +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()) + if ( trTria.GetArea() < EPS_SMALL) return false ; - // Se il voxel sta sulla frontiera contino, - if ( IsAVoxelOnBoundary( nBlockLimits, nVoxIJK, true)) { - - double dSqEps = EPS_SMALL * EPS_SMALL ; + // Se il voxel non è sulla frontiera esco + if ( ! IsAVoxelOnBoundary( nBlockLimits, nVoxIJK, true)) + return false ; - Point3d ptFirstGrPt, ptSecondGrPt ; - - int nCount = 0 ; - - if ( SqDist( trTria.GetP( 0), ptVert) > dSqEps) { - - ptFirstGrPt = trTria.GetP( 0) ; + Point3d ptFirstGrPt, ptSecondGrPt ; + int nCount = 0 ; + if ( SqDist( trTria.GetP( 0), ptVert) > SQ_EPS_SMALL) { + ptFirstGrPt = trTria.GetP( 0) ; + nCount ++ ; + } + if ( SqDist( trTria.GetP( 1), ptVert) > SQ_EPS_SMALL) { + if ( nCount == 0) { + ptFirstGrPt = trTria.GetP( 1) ; 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 ++ ; - } + 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 + if ( SqDist( trTria.GetP( 2), ptVert) > SQ_EPS_SMALL) { + if ( nCount == 1) { + ptSecondGrPt = trTria.GetP( 2) ; + nCount ++ ; + } + } + + // Se non ho trovato due punti, errore + if ( nCount != 2) return false ; + + // Verifico se stanno sulla griglia + 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 false ; } //---------------------------------------------------------------------------- @@ -2894,70 +2878,418 @@ 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 ; - } - } + // Se il voxel non è di frontiera, esco + if ( ! bBorderBox) return false ; + + // Se il triangolo ha area nulla, esco + if ( trTria.GetArea() < EPS_SMALL) + return false ; + + // Determino i punti del triangolo sulla griglia + int nNotVert[2] ; + int nCount = 0 ; + for ( int nV = 0 ; nV < 3 ; ++ nV) { + if ( SqDist( trTria.GetP( nV), ptVert) > SQ_EPS_SMALL) { + if ( nCount == 0) + nNotVert[nCount] = nV ; + else if ( nCount == 1) + nNotVert[nCount] = nV ; + nCount ++ ; + } } - // altrimenti non ha senso continuare - else + + // Se non ho trovato due punti, 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 sulla 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 ; +} + +//---------------------------------------------------------------------------- +bool +VolZmap::ProcessVoxContXY( VoxelContainer& VoxContXY, bool bPlus, TRIA3DLIST& lstTria) const +{ + VoxelContainer::iterator it = VoxContXY.begin() ; + while ( it != VoxContXY.end()) { + + int nN = ( *it).first ; + int nI, nJ, nK ; + GetVoxIJKFromN( nN, nI, nJ, nK) ; + + // Costruzione del primo rettangolo: un singolo voxel + int nMinI = nI ; + int nMinJ = nJ ; + int nMaxI = nI ; + int nMaxJ = nJ ; + double dCordZ = ( *it).second ; + + // Flag sul ritrovamento di un rettangolo più grande. + bool bOkI = true ; + bool bOkJ = true ; + while ( bOkI || bOkJ) { + // Se precedente espansione ok e lato I più lungo o non ok J + if ( bOkI && ( nMaxI - nMinI >= nMaxJ - nMinJ || ! bOkJ)) { + // Analizzo linea superiore + bool bSupJ = true ; + for ( int i = nMinI ; bSupJ && i <= nMaxI ; ++ i) { + if ( ! Find( VoxContXY, i, nMaxJ + 1, nK, dCordZ)) + bSupJ = false ; + } + // Se linea superiore accettata, elimino i voxel + if ( bSupJ) { + for ( int i = nMinI ; i <= nMaxI ; ++ i) + Remove( VoxContXY, i, nMaxJ + 1, nK) ; + ++ nMaxJ ; + } + // Analizzo linea inferiore + bool bInfJ = true ; + for ( int i = nMinI ; bInfJ && i <= nMaxI ; ++ i) { + if ( ! Find( VoxContXY, i, nMinJ - 1, nK, dCordZ)) + bInfJ = false ; + } + // Se linea inferiore accettata, elimino i voxel + if ( bInfJ) { + for ( int i = nMinI ; i <= nMaxI ; ++ i) + Remove( VoxContXY, i, nMinJ - 1, nK) ; + -- nMinJ ; + } + bOkI = bSupJ || bInfJ ; + } + // altrimenti se precedente espansione J ok + else if ( bOkJ) { + // Analizzo linea destra + bool bSupI = true ; + for ( int j = nMinJ ; bSupI && j <= nMaxJ ; ++ j) { + if ( ! Find( VoxContXY, nMaxI + 1, j, nK, dCordZ)) + bSupI = false ; + } + // Se linea destra accettata, elimino i voxel + if ( bSupI) { + for ( int j = nMinJ ; j <= nMaxJ ; ++ j) + Remove( VoxContXY, nMaxI + 1, j, nK) ; + ++ nMaxI ; + } + // Analizzo linea sinistra + bool bInfI = true ; + for ( int j = nMinJ ; bInfI && j <= nMaxJ ; ++ j) { + if ( ! Find( VoxContXY, nMinI - 1, j, nK, dCordZ)) + bInfI = false ; + } + // Se linea sinistra accettata, elimino i voxel + if ( bInfI) { + for ( int j = nMinJ ; j <= nMaxJ ; ++ j) + Remove( VoxContXY, nMinI - 1, j, nK) ; + -- nMinI ; + } + bOkJ = bSupI || bInfI ; + } + } + + Point3d ptT0( ( nMinI + 0.5) * m_dStep, ( ( nMinJ + 0.5) * m_dStep), dCordZ) ; + Point3d ptT1( ( nMaxI + 1.5) * m_dStep, ( ( nMinJ + 0.5) * m_dStep), dCordZ) ; + Point3d ptT2( ( nMaxI + 1.5) * m_dStep, ( ( nMaxJ + 1.5) * m_dStep), dCordZ) ; + Point3d ptT3( ( nMinI + 0.5) * m_dStep, ( ( nMaxJ + 1.5) * m_dStep), dCordZ) ; + + ptT0.ToGlob( m_MapFrame[0]) ; + ptT1.ToGlob( m_MapFrame[0]) ; + ptT2.ToGlob( m_MapFrame[0]) ; + ptT3.ToGlob( m_MapFrame[0]) ; + + Triangle3d Tria0, Tria1 ; + if ( bPlus) { + Tria0.Set( ptT0, ptT1, ptT3) ; + Tria1.Set( ptT1, ptT2, ptT3) ; + } + else { + Tria0.Set( ptT0, ptT3, ptT1) ; + Tria1.Set( ptT1, ptT3, ptT2) ; + } + bool bV0 = Tria0.Validate( true) ; + bool bV1 = Tria1.Validate( true) ; + // Aggiungo alla lista + lstTria.emplace_back( Tria0) ; + lstTria.emplace_back( Tria1) ; + + // Elimino il voxel da cui sono partito a ingrandire. + VoxContXY.erase( nN) ; + + // Passo al primo voxel rimasto + it = VoxContXY.begin() ; + } + + return true ; +} + +//---------------------------------------------------------------------------- +bool +VolZmap::ProcessVoxContYZ( VoxelContainer& VoxContYZ, bool bPlus, TRIA3DLIST& lstTria) const +{ + VoxelContainer::iterator it = VoxContYZ.begin() ; + while ( it != VoxContYZ.end()) { + + int nN = ( *it).first ; + int nI, nJ, nK ; + GetVoxIJKFromN( nN, nI, nJ, nK) ; + + // Costruzione del primo rettangolo: un singolo voxel + int nMinJ = nJ ; + int nMinK = nK ; + int nMaxJ = nJ ; + int nMaxK = nK ; + double dCordX = ( *it).second ; + + // Flag sul ritrovamento di un rettangolo più grande. + bool bOkJ = true ; + bool bOkK = true ; + while ( bOkJ || bOkK) { + // Se precedente espansione ok e lato J più lungo o non ok K + if ( bOkJ && ( nMaxJ - nMinJ >= nMaxK - nMinK || ! bOkK)) { + // Analizzo linea superiore + bool bSupK = true ; + for ( int j = nMinJ ; bSupK && j <= nMaxJ ; ++ j) { + if ( ! Find( VoxContYZ, nI, j, nMaxK + 1, dCordX)) + bSupK = false ; + } + // Se linea superiore accettata, elimino i voxel + if ( bSupK) { + for ( int j = nMinJ ; j <= nMaxJ ; ++ j) + Remove( VoxContYZ, nI, j, nMaxK + 1) ; + ++ nMaxK ; + } + // Analizzo linea inferiore + bool bInfK = true ; + for ( int j = nMinJ ; bInfK && j <= nMaxJ ; ++ j) { + if ( ! Find( VoxContYZ, nI, j, nMinK - 1, dCordX)) + bInfK = false ; + } + // Se linea inferiore accettata, elimino i voxel + if ( bInfK) { + for ( int j = nMinJ ; j <= nMaxJ ; ++ j) + Remove( VoxContYZ, nI, j, nMinK - 1) ; + -- nMinK ; + } + bOkJ = bSupK || bInfK ; + } + // altrimenti se precedente espansione K ok + else if ( bOkK) { + // Analizzo linea destra + bool bSupJ = true ; + for ( int k = nMinK ; bSupJ && k <= nMaxK ; ++ k) { + if ( ! Find( VoxContYZ, nI, nMaxJ + 1, k, dCordX)) + bSupJ = false ; + } + // Se linea destra accettata, elimino i voxel + if ( bSupJ) { + for ( int k = nMinK ; k <= nMaxK ; ++ k) + Remove( VoxContYZ, nI, nMaxJ + 1, k) ; + ++ nMaxJ ; + } + // Analizzo linea sinistra + bool bInfJ = true ; + for ( int k = nMinK ; bInfJ && k <= nMaxK ; ++ k) { + if ( ! Find( VoxContYZ, nI, nMinJ - 1, k, dCordX)) + bInfJ = false ; + } + // Se linea sinistra accettata, elimino i voxel + if ( bInfJ) { + for ( int k = nMinK ; k <= nMaxK ; ++ k) + Remove( VoxContYZ, nI, nMinJ - 1, k) ; + -- nMinJ ; + } + bOkK = bSupJ || bInfJ ; + } + } + + Point3d ptT0( dCordX, ( nMinJ + 0.5) * m_dStep, ( ( nMinK + 0.5) * m_dStep)) ; + Point3d ptT1( dCordX, ( nMaxJ + 1.5) * m_dStep, ( ( nMinK + 0.5) * m_dStep)) ; + Point3d ptT2( dCordX, ( nMaxJ + 1.5) * m_dStep, ( ( nMaxK + 1.5) * m_dStep)) ; + Point3d ptT3( dCordX, ( nMinJ + 0.5) * m_dStep, ( ( nMaxK + 1.5) * m_dStep)) ; + + ptT0.ToGlob( m_MapFrame[0]) ; + ptT1.ToGlob( m_MapFrame[0]) ; + ptT2.ToGlob( m_MapFrame[0]) ; + ptT3.ToGlob( m_MapFrame[0]) ; + + Triangle3d Tria0, Tria1 ; + if ( bPlus) { + Tria0.Set( ptT0, ptT1, ptT3) ; + Tria1.Set( ptT1, ptT2, ptT3) ; + } + else { + Tria0.Set( ptT0, ptT3, ptT1) ; + Tria1.Set( ptT1, ptT3, ptT2) ; + } + bool bV0 = Tria0.Validate( true) ; + bool bV1 = Tria1.Validate( true) ; + // Aggiungo alla lista + lstTria.emplace_back( Tria0) ; + lstTria.emplace_back( Tria1) ; + + // Elimino il voxel da cui sono partito a ingrandire. + VoxContYZ.erase( nN) ; + + // Passo al primo voxel rimasto + it = VoxContYZ.begin() ; + } + + return true ; +} + +//---------------------------------------------------------------------------- +bool +VolZmap::ProcessVoxContXZ( VoxelContainer& VoxContXZ, bool bPlus, TRIA3DLIST& lstTria) const +{ + VoxelContainer::iterator it = VoxContXZ.begin() ; + while ( it != VoxContXZ.end()) { + + int nN = ( *it).first ; + int nI, nJ, nK ; + GetVoxIJKFromN( nN, nI, nJ, nK) ; + + // Costruzione del primo rettangolo: un singolo voxel + int nMinI = nI ; + int nMinK = nK ; + int nMaxI = nI ; + int nMaxK = nK ; + double dCordY = ( *it).second ; + + // Flag sul ritrovamento di un rettangolo più grande. + bool bOkI = true ; + bool bOkK = true ; + while ( bOkI || bOkK) { + // Se lato I più lungo e precedente espansione ok + if ( bOkI && ( nMaxI - nMinI >= nMaxK - nMinK || ! bOkK)) { + // Analizzo linea superiore + bool bSupK = true ; + for ( int i = nMinI ; bSupK && i <= nMaxI ; ++ i) { + if ( ! Find( VoxContXZ, i, nJ, nMaxK + 1, dCordY)) + bSupK = false ; + } + // Se linea superiore accettata, elimino i voxel + if ( bSupK) { + for ( int i = nMinI ; i <= nMaxI ; ++ i) + Remove( VoxContXZ, i, nJ, nMaxK + 1) ; + ++ nMaxK ; + } + // Analizzo linea inferiore + bool bInfK = true ; + for ( int i = nMinI ; bInfK && i <= nMaxI ; ++ i) { + if ( ! Find( VoxContXZ, i, nJ, nMinK - 1, dCordY)) + bInfK = false ; + } + // Se linea inferiore accettata, elimino i voxel + if ( bInfK) { + for ( int i = nMinI ; i <= nMaxI ; ++ i) + Remove( VoxContXZ, i, nJ, nMinK - 1) ; + -- nMinK ; + } + bOkI = bSupK || bInfK ; + } + // altrimenti se precedente espansione K ok + else if ( bOkK) { + // Analizzo linea destra + bool bSupI = true ; + for ( int k = nMinK ; bSupI && k <= nMaxK ; ++ k) { + if ( ! Find( VoxContXZ, nMaxI + 1, nJ, k, dCordY)) + bSupI = false ; + } + // Se linea destra accettata, elimino i voxel + if ( bSupI) { + for ( int k = nMinK ; k <= nMaxK ; ++ k) + Remove( VoxContXZ, nMaxI + 1, nJ, k) ; + ++ nMaxI ; + } + // Analizzo linea sinistra + bool bInfI = true ; + for ( int k = nMinK ; bInfI && k <= nMaxK ; ++ k) { + if ( ! Find( VoxContXZ, nMinI - 1, nJ, k, dCordY)) + bInfI = false ; + } + // Se linea sinistra accettata, elimino i voxel + if ( bInfI) { + for ( int k = nMinK ; k <= nMaxK ; ++ k) + Remove( VoxContXZ, nMinI - 1, nJ, k) ; + -- nMinI ; + } + bOkK = bSupI || bInfI ; + } + } + + Point3d ptT0( ( nMinI + 0.5) * m_dStep, dCordY, ( ( nMinK + 0.5) * m_dStep)) ; + Point3d ptT1( ( nMaxI + 1.5) * m_dStep, dCordY, ( ( nMinK + 0.5) * m_dStep)) ; + Point3d ptT2( ( nMaxI + 1.5) * m_dStep, dCordY, ( ( nMaxK + 1.5) * m_dStep)) ; + Point3d ptT3( ( nMinI + 0.5) * m_dStep, dCordY, ( ( nMaxK + 1.5) * m_dStep)) ; + + ptT0.ToGlob( m_MapFrame[0]) ; + ptT1.ToGlob( m_MapFrame[0]) ; + ptT2.ToGlob( m_MapFrame[0]) ; + ptT3.ToGlob( m_MapFrame[0]) ; + + Triangle3d Tria0, Tria1 ; + if ( bPlus) { + Tria0.Set( ptT0, ptT3, ptT1) ; + Tria1.Set( ptT1, ptT3, ptT2) ; + } + else { + Tria0.Set( ptT0, ptT1, ptT3) ; + Tria1.Set( ptT1, ptT2, ptT3) ; + } + bool bV0 = Tria0.Validate( true) ; + bool bV1 = Tria1.Validate( true) ; + // Aggiungo alla lista + lstTria.emplace_back( Tria0) ; + lstTria.emplace_back( Tria1) ; + + // Elimino il voxel da cui sono partito a ingrandire. + VoxContXZ.erase( nN) ; + + // Passo al primo voxel rimasto + it = VoxContXZ.begin() ; + } + + return true ; +} + +//---------------------------------------------------------------------------- +bool +VolZmap::Find( const VoxelContainer& VoxCont, int nI, int nJ, int nK, double dPos) const +{ + // indice globale del voxel + int nN ; + if ( ! GetVoxNFromIJK( nI, nJ, nK, nN)) + return false ; + // cerco il voxel nel contenitore + auto iter = VoxCont.find( nN) ; + return ( iter != VoxCont.end() && abs( dPos - ( *iter).second) < EPS_SMALL) ; +} + +//---------------------------------------------------------------------------- +bool +VolZmap::Remove( VoxelContainer& VoxCont, int nI, int nJ, int nK) const +{ + int nN ; + if ( ! GetVoxNFromIJK( nI, nJ, nK, nN)) + return false ; + return ( VoxCont.erase( nN) > 0) ; } diff --git a/VolZmap.h b/VolZmap.h index 23d229e..cc54e0c 100644 --- a/VolZmap.h +++ b/VolZmap.h @@ -22,6 +22,7 @@ #include "/EgtDev/Include/EGkVector3d.h" #include "/EgtDev/Include/ENkPolynomialRoots.h" #include "/EgtDev/Include/EgkIntersLinesurfTm.h" +#include //---------------------------------------------------------------------------- struct TriaStruct { @@ -35,8 +36,7 @@ typedef std::vector TriHolder ; // Vettori di TriHolder con sharp feature di frontiera: // il primo indice individua il blocco, il secondo il voxel -typedef std::vector TriaMatrix ; - +typedef std::vector TriaMatrix ; //---------------------------------------------------------------------------- class VolZmap : public IVolZmap, public IGeoObjRW @@ -118,6 +118,7 @@ class VolZmap : public IVolZmap, public IGeoObjRW CONUSMILL = 4, // con parte terminale conica MORTISER = 5, // mortasatrice CHISEL = 6} ; // scalpello + typedef std::unordered_map VoxelContainer ; private : bool CopyFrom( const VolZmap& clSrc) ; @@ -132,7 +133,7 @@ class VolZmap : public IVolZmap, public IGeoObjRW bool FlipEdgesII( TriHolder& TriHold) const ; bool FlipEdgesBB( TriaMatrix& InterTria) const ; bool IsThereMat( int nI, int nJ, int nK) const ; - bool IsThereMat( const int nMatr[][3], int nNum, double & dHx, double & dHy, double & dHz) const ; + int CalcIndex( int nI, int nJ, int nK) const ; bool IntersPos( int nVec1[], int nVec2[], Point3d & ptInt) const ; bool IntersPos( int nVec1[], int nVec2[], bool bFirstCorner, Point3d& ptInt, Vector3d& vtNormal) const ; bool IsPointInsideVoxel( int nI, int nJ, int nK, const Point3d& ptP) const ; @@ -152,63 +153,39 @@ class VolZmap : public IVolZmap, public IGeoObjRW // SOTTRAZIONI // UTENSILI // Asse di simmetria parallelo a Z - - // Cilindro Sfera bool CylBall_ZDrilling( unsigned int nGrid, const Point3d& ptS, const Point3d& ptE, const Vector3d& vtToolDir) ; bool CylBall_ZPerp( unsigned int nGrid, const Point3d& ptS, const Point3d& ptE, const Vector3d& vtToolDir) ; bool CylBall_ZMilling( unsigned int nGrid, const Point3d& ptS, const Point3d& ptE, const Vector3d& vtToolDir) ; - - // Coni bool Conus_ZDrilling( unsigned int nGrid, const Point3d & ptS, const Point3d& ptE, const Vector3d& vtToolDir) ; bool Conus_ZPerp( unsigned int nGrid, const Point3d& ptS, const Point3d& ptE, const Vector3d& vtToolDir) ; bool Conus_ZMilling( unsigned int nGrid, const Point3d& ptS, const Point3d& ptE, const Vector3d& vtToolDir) ; - - // Mortasatrice bool Mrt_ZDrilling( unsigned int nGrid, const Point3d& ptS, const Point3d& ptE, const Vector3d& vtToolDir, const Vector3d& vtAux) ; bool Mrt_ZMilling( unsigned int nGrid, const Point3d& ptS, const Point3d& ptE, const Vector3d& vtToolDir, const Vector3d& vtAux) ; - - // Chisel bool Chs_ZDrilling( unsigned int nGrid, const Point3d& ptS, const Point3d& ptE, const Vector3d& vtToolDir, const Vector3d& vtAux) ; bool Chs_ZMilling( unsigned int nGrid, const Point3d& ptS, const Point3d& ptE, const Vector3d& vtToolDir, const Vector3d& vtAux) ; - - // Utensile generico bool GenTool_ZDrilling( unsigned int nGrid, const Point3d& ptS, const Point3d& ptE, const Vector3d& vtToolDir) ; bool GenTool_ZMilling( unsigned int nGrid, const Point3d& ptS, const Point3d& ptE, const Vector3d& vtToolDir) ; - // Asse di simmetria nel piano bool CylBall_XYDrilling( unsigned int nGrid, const Point3d& ptS, const Point3d& ptE, const Vector3d& vtToolDir) ; bool CylBall_XYPerp( unsigned int nGrid, const Point3d& ptS, const Point3d& ptE, const Vector3d& vtToolDir) ; bool CylBall_XYMilling( unsigned int nGrid, const Point3d& ptS, const Point3d& ptE, const Vector3d& vtToolDir) ; - bool Conus_XYDrilling( unsigned int nGrid, const Point3d& ptS, const Point3d& ptE, const Vector3d& vtToolDir) ; bool Conus_XYPerp( unsigned int nGrid, const Point3d& ptS, const Point3d& ptE, const Vector3d& vtToolDir) ; bool Conus_XYMilling( unsigned int nGrid, const Point3d& ptS, const Point3d& ptE, const Vector3d& vtToolDir) ; - // Asse di simmetria con orientazione generica - // Cilindro e sfera bool CylBall_Drilling( unsigned int nGrid, const Point3d& ptS, const Point3d& ptE, const Vector3d& vtToolDir) ; bool CylBall_Milling( unsigned int nGrid, const Point3d& ptS, const Point3d& ptE, const Vector3d& vtToolDir) ; - - // Coni bool Conus_Drilling( unsigned int nGrid, const Point3d& ptS, const Point3d& ptE, const Vector3d& vtToolDir) ; bool Conus_Milling( unsigned int nGrid, const Point3d& ptS, const Point3d& ptE, const Vector3d& vtToolDir) ; - - // Mortasatrice bool Mrt_Drilling( unsigned int nGrid, const Point3d& ptS, const Point3d& ptE, const Vector3d& vtToolDir, const Vector3d& vtAux) ; bool Mrt_Milling( unsigned int nGrid, const Point3d& ptS, const Point3d& ptE, const Vector3d& vtToolDir, const Vector3d& vtAux) ; // E' in realtà un Perp - - // Chisel bool Chs_Drilling( unsigned int nGrid, const Point3d& ptS, const Point3d& ptE, const Vector3d& vtToolDir, const Vector3d& vtAux) ; bool Chs_Milling( unsigned int nGrid, const Point3d& ptS, const Point3d& ptE, const Vector3d& vtToolDir, const Vector3d& vtAux) ; // E' in realtà un Perp - - // Utensile generico bool GenTool_Drilling( unsigned int nGrid, const Point3d& ptS, const Point3d& ptE, const Vector3d& vtToolDir) ; bool GenTool_Milling( unsigned int nGrid, const Point3d& ptS, const Point3d& ptE, const Vector3d& vtToolDir) ; // COMPONENTI // Asse di simmetria diretto come l'asse Z - - // Drilling bool CompCyl_ZDrilling( unsigned int nGrid, const Point3d& ptLs, const Point3d& ptLe, const Vector3d& vtToolDir, double dHei, double dRad) ; bool CompConus_ZDrilling( unsigned int nGrid, const Point3d& ptLs, const Point3d& ptLe, const Vector3d& vtToolDir, @@ -216,8 +193,6 @@ class VolZmap : public IVolZmap, public IGeoObjRW bool CompPar_ZDrilling( unsigned int nGrid, double dLenX, double dLenY, double dLenZ, const Point3d& ptS, const Point3d& ptE, const Vector3d& vtToolDir, const Vector3d& vtAux) ; - - // Milling bool CompCyl_ZMilling( unsigned int nGrid, const Point3d& ptS, const Point3d& ptE, const Vector3d& vtToolDir, double dHei, double dRad) ; bool CompConus_ZMilling( unsigned int nGrid, const Point3d& ptS, const Point3d& ptE, const Vector3d& vtToolDir, @@ -225,11 +200,7 @@ class VolZmap : public IVolZmap, public IGeoObjRW bool CompPar_ZMilling( unsigned int nGrid, double dLenX, double dLenY, double dLenZ, const Point3d& ptS, const Point3d& ptE, const Vector3d& vtToolDir, const Vector3d& vtAux) ; // E' in realtà MillingPerp - - // Asse di simmetria con orientazione generica - - // Drilling bool CompCyl_Drilling( unsigned int nGrid, const Point3d& ptS, const Point3d& ptE, const Vector3d& vtToolDir, double dHei, double dRad, bool bTapB, bool bTapT) ; bool CompConus_Drilling( unsigned int nGrid, const Point3d& ptS, const Point3d& ptE, const Vector3d& vtToolDir, @@ -237,8 +208,6 @@ class VolZmap : public IVolZmap, public IGeoObjRW bool CompPar_Drilling( unsigned int nGrid, double dLenX, double dLenY, double dLenZ, const Point3d& ptS, const Point3d& ptE, const Vector3d& vtToolDir, const Vector3d& vtAux) ; - - // Milling bool CompCyl_Milling( unsigned int nGrid, const Point3d& ptS, const Point3d& ptE, const Vector3d& vtToolDir, double dHei, double dRad, bool bTapB, bool bTapT) ; bool CompConus_Milling( unsigned int nGrid, const Point3d& ptS, const Point3d& ptE, const Vector3d& vtToolDir, @@ -246,10 +215,8 @@ class VolZmap : public IVolZmap, public IGeoObjRW bool CompPar_Milling( unsigned int nGrid, double dLenX, double dLenY, double dLenZ, const Point3d& ptS, const Point3d& ptE, const Vector3d& vtToolDir, const Vector3d& vtAux) ; // E' in realtà MillingPerp - // Generica traslazione sfera bool CompBall_Milling( unsigned int nGrid, const Point3d& ptS, const Point3d& ptE, double dRad) ; - // BBox per utensili e solidi semplici in movimenti di traslazione pura; in futuro si dovrà implementare le routine che prevedono rotazioni inline bool BoundingBox( unsigned int nGrid, const Point3d& ptP1, const Point3d& ptP2, const Vector3d& vtV, unsigned int& nStI, unsigned int& nStJ, unsigned int& nEnI, unsigned int& nEnJ) ; @@ -259,8 +226,6 @@ class VolZmap : public IVolZmap, public IGeoObjRW inline bool BBoxParaComp( unsigned int nGrid, double dLenX, double dLenY, double dLenZ, const Point3d& ptS, const Point3d& ptE, const Vector3d& vtD, const Vector3d& vtA, unsigned int& nStI, unsigned int& nStJ, unsigned int& nEnI, unsigned int& nEnJ) ; - - // Intersezioni bool IntersLineBox( const Point3d& ptP, const Vector3d& vtV, const Point3d& ptMin, const Point3d& ptMax, double& dU1, double& dU2) const ; @@ -290,7 +255,9 @@ class VolZmap : public IVolZmap, public IGeoObjRW bool IntersLineMyPolyhedron( const Point3d& ptLineSt, const Vector3d& vtLineDir, const Frame3d& PolyFrame, double dLenX, double dLenY, double dLenZ, double dDeltaX, Point3d& ptInt1, Point3d& ptInt2, Vector3d& vtN1, Vector3d& vtN2) ; - + // Passaggio da N a ijk per i voxel + bool GetVoxIJKFromN( int nN, int& nI, int& nJ, int& nK) const ; + bool GetVoxNFromIJK( int nI, int nJ, int nK, int& nN) const ; // Funzioni di gestione dei blocchi bool GetBlockIJKFromN( int nBlock, int nIJK[]) const ; bool GetBlockNFromIJK( int nIJK[], int& nBlock) const ; @@ -302,6 +269,12 @@ class VolZmap : public IVolZmap, public IGeoObjRW const int nBlockLimits[], const int nVoxIJK[]) const ; bool IsATriangleOnBorder( const Triangle3d& trTria, const Point3d& ptVert, const int nBlockLimits[], const int nVoxIJK[], bool bBorderBox) const ; + // Funzioni per facce canoniche con grandi triangoli + bool ProcessVoxContXY( VoxelContainer& VoxContXY, bool bPlus, TRIA3DLIST& lstTria) const ; + bool ProcessVoxContYZ( VoxelContainer& VoxContYZ, bool bPlus, TRIA3DLIST& lstTria) const ; + bool ProcessVoxContXZ( VoxelContainer& VoxContXZ, bool bPlus, TRIA3DLIST& lstTria) const ; + bool Find( const VoxelContainer& VoxCont, int nI, int nJ, int nK, double dPos) const ; + bool Remove( VoxelContainer& VoxCont, int nI, int nJ, int nK) const ; // Connessione Zmap bool IsMapConnected( void) ; @@ -360,7 +333,6 @@ class VolZmap : public IVolZmap, public IGeoObjRW double m_dMrtChsWidth ; double m_dMrtChsThickness ; - mutable TriaMatrix m_InterBlockTria ; } ;