diff --git a/EgtGeomKernel.rc b/EgtGeomKernel.rc index b1ef3ea..6bc37a9 100644 Binary files a/EgtGeomKernel.rc and b/EgtGeomKernel.rc differ diff --git a/VolTriZmapCalculus.cpp b/VolTriZmapCalculus.cpp index 7ef952f..bac53c0 100644 --- a/VolTriZmapCalculus.cpp +++ b/VolTriZmapCalculus.cpp @@ -71,6 +71,352 @@ VolZmap::IntersLineBox( const Point3d& ptP, const Vector3d& vtV, return ( dU2 >= dU1) ; } +/* +//---------------------------------------------------------------------------- +bool +VolZmap::IntersLineVoxel( const Point3d& ptP, const Vector3d& vtV, int nIndI, int nIndJ, int nIndK, + int& nFace1, int& nFace2, double& dU1, double& dU2) const +{ + // Controllo sull'ammissibilità del voxel + if ( nIndI < - 1 || nIndI >= int( m_nNx[0]) || + nIndJ < - 1 || nIndJ >= int( m_nNy[0]) || + nIndK < - 1 || nIndK >= int( m_nNy[1])) + return false ; + + Point3d ptInt ; + int nIntNum = 0 ; + double dU ; + + // Intersezione con le facce 1 e 3 + if ( abs( vtV.y) > EPS_ZERO) { + + // Intersezione con la prima faccia + dU1 = ( ( nIndJ + 0.5) * m_dStep - ptP.y) / vtV.y ; + ptInt = ptP + dU1 * vtV ; + + if ( ptInt.x > ( nIndI + 0.5) * m_dStep - EPS_SMALL && + ptInt.x < ( nIndI + 1.5) * m_dStep + EPS_SMALL && + ptInt.z > ( nIndK + 0.5) * m_dStep - EPS_SMALL && + ptInt.z < ( nIndK + 1.5) * m_dStep + EPS_SMALL) { + + nFace1 = 1 ; + ++ nIntNum ; + } + + // Intersezione con la terza faccia + dU = ( ( nIndJ + 1.5) * m_dStep - ptP.y) / vtV.y ; + ptInt = ptP + dU * vtV ; + + if ( ptInt.x > ( nIndI + 0.5) * m_dStep - EPS_SMALL && + ptInt.x < ( nIndI + 1.5) * m_dStep + EPS_SMALL && + ptInt.z > ( nIndK + 0.5) * m_dStep - EPS_SMALL && + ptInt.z < ( nIndK + 1.5) * m_dStep + EPS_SMALL) { + + if ( nIntNum == 0) { + + dU1 = dU ; + nFace1 = 3 ; + ++ nIntNum ; + } + else if ( nIntNum == 1) { + + dU2 = dU ; + nFace2 = 3 ; + ++ nIntNum ; + } + } + } + + // Intersezione con le facce 2 e 4 + if ( abs( vtV.x) > EPS_ZERO) { + + // Intersezione con la seconda faccia + dU = ( ( nIndI + 0.5) * m_dStep - ptP.x) / vtV.x ; + ptInt = ptP + dU * vtV ; + + if ( ptInt.y > ( nIndJ + 0.5) * m_dStep - EPS_SMALL && + ptInt.y < ( nIndJ + 1.5) * m_dStep + EPS_SMALL && + ptInt.z > ( nIndK + 0.5) * m_dStep - EPS_SMALL && + ptInt.z < ( nIndK + 1.5) * m_dStep + EPS_SMALL) { + + if ( nIntNum == 0) { + + dU1 = dU ; + nFace1 = 2 ; + ++ nIntNum ; + } + else if ( nIntNum == 1) { + + dU2 = dU ; + nFace2 = 2 ; + ++ nIntNum ; + } + } + + // Intersezione con la quarta faccia + dU = ( ( nIndI + 1.5) * m_dStep - ptP.x) / vtV.x ; + ptInt = ptP + dU * vtV ; + + if ( ptInt.y > ( nIndJ + 0.5) * m_dStep - EPS_SMALL && + ptInt.y < ( nIndJ + 1.5) * m_dStep + EPS_SMALL && + ptInt.z > ( nIndK + 0.5) * m_dStep - EPS_SMALL && + ptInt.z < ( nIndK + 1.5) * m_dStep + EPS_SMALL) { + + if ( nIntNum == 0) { + + dU1 = dU ; + nFace1 = 4 ; + ++ nIntNum ; + } + else if ( nIntNum == 1) { + + dU2 = dU ; + nFace2 = 4 ; + ++ nIntNum ; + } + } + } + + // Intersezione con le facce 5 e 6 + if ( abs( vtV.z) > EPS_ZERO) { + + // Intersezione con la quinta faccia + dU = ( ( nIndK + 0.5) * m_dStep - ptP.z) / vtV.z ; + ptInt = ptP + dU * vtV ; + + if ( ptInt.x > ( nIndI + 0.5) * m_dStep - EPS_SMALL && + ptInt.x < ( nIndI + 1.5) * m_dStep + EPS_SMALL && + ptInt.y > ( nIndJ + 0.5) * m_dStep - EPS_SMALL && + ptInt.y < ( nIndJ + 1.5) * m_dStep + EPS_SMALL) { + + if ( nIntNum == 0) { + + dU1 = dU ; + nFace1 = 5 ; + ++ nIntNum ; + } + else if ( nIntNum == 1) { + + dU2 = dU ; + nFace2 = 5 ; + ++ nIntNum ; + } + } + + // Intersezione con la sesta faccia + dU = ( ( nIndK + 1.5) * m_dStep - ptP.z) / vtV.z ; + ptInt = ptP + dU * vtV ; + + if ( ptInt.x > ( nIndI + 0.5) * m_dStep - EPS_SMALL && + ptInt.x < ( nIndI + 1.5) * m_dStep + EPS_SMALL && + ptInt.y > ( nIndJ + 0.5) * m_dStep - EPS_SMALL && + ptInt.y < ( nIndJ + 1.5) * m_dStep + EPS_SMALL) { + + + if ( nIntNum == 0) { + + dU1 = dU ; + nFace1 = 6 ; + ++ nIntNum ; + } + else if ( nIntNum == 1) { + + dU2 = dU ; + nFace2 = 6 ; + ++ nIntNum ; + } + } + } + + if ( dU1 > dU2) { + + swap( dU1, dU2) ; + swap( nFace1, nFace2) ; + } + return ( nIntNum == 2) ; +} */ + +//---------------------------------------------------------------------------- +bool +VolZmap::IntersLineVoxel( const Point3d& ptP, const Vector3d& vtV, int nIndI, int nIndJ, int nIndK, + int& nFaceF, int& nFaceL, double& dUF, double& dUL) const +{ + // Controllo sull'ammissibilità del voxel + if ( nIndI < - 1 || nIndI >= int( m_nNx[0]) || + nIndJ < - 1 || nIndJ >= int( m_nNy[0]) || + nIndK < - 1 || nIndK >= int( m_nNy[1])) + return false ; + + Point3d ptInt, ptIntF, ptIntL ; + double dSqEps = EPS_SMALL * EPS_SMALL ; + int nIntNum = 0 ; + + double dU1 = ( ( nIndJ + 0.5) * m_dStep - ptP.y) / vtV.y ; + double dU2 = ( ( nIndI + 0.5) * m_dStep - ptP.x) / vtV.x ; + double dU3 = ( ( nIndJ + 1.5) * m_dStep - ptP.y) / vtV.y ; + double dU4 = ( ( nIndI + 1.5) * m_dStep - ptP.x) / vtV.x ; + double dU5 = ( ( nIndK + 0.5) * m_dStep - ptP.z) / vtV.z ; + double dU6 = ( ( nIndK + 1.5) * m_dStep - ptP.z) / vtV.z ; + + // Intersezione con le facce 1 e 3 + if ( abs( vtV.y) > EPS_ZERO) { + + // Intersezione con la prima faccia + ptInt = ptP + dU1 * vtV ; + + if ( ptInt.x >= ( nIndI + 0.5) * m_dStep && + ptInt.x <= ( nIndI + 1.5) * m_dStep && + ptInt.z >= ( nIndK + 0.5) * m_dStep && + ptInt.z <= ( nIndK + 1.5) * m_dStep) { + + dUF = dU1 ; + nFaceF = 1 ; + ptIntF = ptInt ; + ++ nIntNum ; + } + + // Intersezione con la terza faccia + ptInt = ptP + dU3 * vtV ; + + if ( ptInt.x >= ( nIndI + 0.5) * m_dStep && + ptInt.x <= ( nIndI + 1.5) * m_dStep && + ptInt.z >= ( nIndK + 0.5) * m_dStep && + ptInt.z <= ( nIndK + 1.5) * m_dStep) { + + if ( nIntNum == 0) { + + dUF = dU3 ; + nFaceF = 3 ; + ptIntF = ptInt ; + ++ nIntNum ; + } + else if ( nIntNum == 1 && + ( ptIntF - ptInt).SqLen() > dSqEps) { + + dUL = dU3 ; + nFaceL = 3 ; + ptIntL = ptInt ; + ++ nIntNum ; + } + } + } + + // Intersezione con le facce 2 e 4 + if ( abs( vtV.x) > EPS_ZERO) { + + // Intersezione con la seconda faccia + ptInt = ptP + dU2 * vtV ; + + if ( ptInt.y >= ( nIndJ + 0.5) * m_dStep && + ptInt.y <= ( nIndJ + 1.5) * m_dStep && + ptInt.z >= ( nIndK + 0.5) * m_dStep && + ptInt.z <= ( nIndK + 1.5) * m_dStep) { + + if ( nIntNum == 0) { + + dUF = dU2 ; + nFaceF = 2 ; + ptIntF = ptInt ; + ++ nIntNum ; + } + else if ( nIntNum == 1 && + ( ptIntF - ptInt).SqLen() > dSqEps) { + + dUL = dU2 ; + nFaceL = 2 ; + ptIntL = ptInt ; + ++ nIntNum ; + } + } + + // Intersezione con la quarta faccia + ptInt = ptP + dU4 * vtV ; + + if ( ptInt.y >= ( nIndJ + 0.5) * m_dStep && + ptInt.y <= ( nIndJ + 1.5) * m_dStep && + ptInt.z >= ( nIndK + 0.5) * m_dStep && + ptInt.z <= ( nIndK + 1.5) * m_dStep) { + + if ( nIntNum == 0) { + + dUF = dU4 ; + nFaceF = 4 ; + ptIntF = ptInt ; + ++ nIntNum ; + } + else if ( nIntNum == 1 && + ( ptIntF - ptInt).SqLen() > dSqEps) { + + dUL = dU4 ; + nFaceL = 4 ; + ptIntL = ptInt ; + ++ nIntNum ; + } + } + } + + // Intersezione con le facce 5 e 6 + if ( abs( vtV.z) > EPS_ZERO) { + + // Intersezione con la quinta faccia + ptInt = ptP + dU5 * vtV ; + + if ( ptInt.x >= ( nIndI + 0.5) * m_dStep && + ptInt.x <= ( nIndI + 1.5) * m_dStep && + ptInt.y >= ( nIndJ + 0.5) * m_dStep && + ptInt.y <= ( nIndJ + 1.5) * m_dStep) { + + if ( nIntNum == 0) { + + dUF = dU5 ; + nFaceF = 5 ; + ptIntF = ptInt ; + ++ nIntNum ; + } + else if ( nIntNum == 1 && + ( ptIntF - ptInt).SqLen() > dSqEps) { + + dUL = dU5 ; + nFaceL = 5 ; + ptIntL = ptInt ; + ++ nIntNum ; + } + } + + // Intersezione con la sesta faccia + ptInt = ptP + dU6 * vtV ; + + if ( ptInt.x >= ( nIndI + 0.5) * m_dStep && + ptInt.x <= ( nIndI + 1.5) * m_dStep && + ptInt.y >= ( nIndJ + 0.5) * m_dStep && + ptInt.y <= ( nIndJ + 1.5) * m_dStep) { + + + if ( nIntNum == 0) { + + dUF = dU6 ; + nFaceF = 6 ; + ptIntF = ptInt ; + ++ nIntNum ; + } + else if ( nIntNum == 1 && + ( ptIntF - ptInt).SqLen() > dSqEps) { + + dUL = dU6 ; + nFaceL = 6 ; + ptIntL = ptInt ; + ++ nIntNum ; + } + } + } + + if ( dUF > dUL) { + + swap( dUF, dUL) ; + swap( nFaceF, nFaceL) ; + } + return ( nIntNum == 2) ; +} //---------------------------------------------------------------------------- bool @@ -1038,7 +1384,7 @@ VolZmap::IntersLineEllipticalCylinder( const Vector3d& vtLineDir, const Point3d& Vector3d vtTan1( 0, - vtCirc1.z, vtCirc1.y) ; Vector3d vtCross1 = vtTan1 ^ vtMv ; - vtN1 = ( vtCross1 * vtCirc1 > 0 ? vtCross1 : - vtCross1) ; + vtN1 = ( vtCross1 * vtCirc1 > - EPS_ZERO ? vtCross1 : - vtCross1) ; if ( vtTest2.y > 0) { @@ -1053,7 +1399,7 @@ VolZmap::IntersLineEllipticalCylinder( const Vector3d& vtLineDir, const Point3d& Vector3d vtTan2( 0, - vtCirc2.z, vtCirc2.y) ; Vector3d vtCross2 = vtTan2 ^ vtMv ; - vtN2 = ( vtCross2 * vtCirc2 > 0 ? vtCross2 : - vtCross2) ; + vtN2 = ( vtCross2 * vtCirc2 > - EPS_ZERO ? vtCross2 : - vtCross2) ; diff --git a/VolTriZmapGraphics.cpp b/VolTriZmapGraphics.cpp index 59182bf..bd4a5eb 100644 --- a/VolTriZmapGraphics.cpp +++ b/VolTriZmapGraphics.cpp @@ -52,15 +52,66 @@ static int NeighbourTable[8][4] = { // ------------------------- FUNZIONE TEST SULLE NORMALI -------------------------------------------------------------------------- -enum FatureType { NoFeature = 0, Corner = 1, Edge = 2} ; +enum FatureType { NoFeature = 0, Corner = 1, Edge = 2} ; //, Edge2 = 3} ; + +//---------------------------------------------------------------------------- +//bool +//TestOnNormal( const VectorField CompoVert[], int nCompoElem, int& FeatureType, Vector3d& vtFeature) +//{ +// int nI, nJ ; +// double dMinCosTheta = 1.001 ; +// const double dCosThetaSharp = sqrt( 3) / 2 ; 0.9 ; +// +// for ( int i = 0 ; i < nCompoElem ; ++ i) { +// +// for ( int j = i + 1 ; j < nCompoElem ; ++ j) { +// +// double dCurrentCos = CompoVert[i].vtNorm * CompoVert[j].vtNorm ; +// +// if ( dCurrentCos < dMinCosTheta) { +// +// nI = i ; +// nJ = j ; +// dMinCosTheta = dCurrentCos ; +// } +// } +// } +// +// if ( dMinCosTheta >= dCosThetaSharp) { +// +// FeatureType = NoFeature ; +// return false ; +// } +// +// Vector3d vtK = CompoVert[nI].vtNorm ^ CompoVert[nJ].vtNorm ; +// +// const double dCosPhiCorner = 0.5;//0.7;//0.5 ; +// +// for ( int i = 0 ; i < nCompoElem ; ++ i) { +// +// double dAbsCurrentCos = abs( CompoVert[i].vtNorm * vtK) ; +// +// if ( dAbsCurrentCos > dCosPhiCorner) { +// +// FeatureType = Corner ; +// vtFeature = vtK ; +// return true ; +// } +// } +// +// FeatureType = Edge ; +// vtFeature = vtK ; +// +// return true ; +//} //---------------------------------------------------------------------------- bool -TestOnNormal( const VectorField CompoVert[], int nCompoElem, int& FeatureType) +TestOnNormal( const VectorField CompoVert[], int nCompoElem, int& FeatureType, Vector3d& vtFeature, int& nMin1, int& nMin2) { int nI, nJ ; double dMinCosTheta = 1.001 ; - const double dCosThetaSharp = 0.9 ; + const double dCosThetaSharp = 0.8; for ( int i = 0 ; i < nCompoElem ; ++ i) { @@ -84,121 +135,37 @@ TestOnNormal( const VectorField CompoVert[], int nCompoElem, int& FeatureType) } Vector3d vtK = CompoVert[nI].vtNorm ^ CompoVert[nJ].vtNorm ; + double dDotTest1 = CompoVert[nI].vtNorm * vtK ; + double dDotTest2 = CompoVert[nJ].vtNorm * vtK ; - const double dCosPhiCorner = 0.5 ; + nMin1 = nI ; + nMin2 = nJ ; + + const double dCosPhiCorner = 0.5 ; + double dMaxAbsCos = 0 ; for ( int i = 0 ; i < nCompoElem ; ++ i) { double dAbsCurrentCos = abs( CompoVert[i].vtNorm * vtK) ; - if ( dAbsCurrentCos > dCosPhiCorner) { - - FeatureType = Corner ; - return true ; - } + if ( dAbsCurrentCos > dMaxAbsCos) + + dMaxAbsCos = dAbsCurrentCos ; } - FeatureType = Edge ; + if ( dMaxAbsCos > dCosPhiCorner) + + FeatureType = Corner ; + + else + + FeatureType = Edge ; + vtFeature = vtK ; + return true ; } -//---------------------------------------------------------------------------- -bool -TestOnSlice( const VectorField CompoVert[], int nLimArrey, - int i, - int j, - int k, double dStep) -{ - // Coordinate dei piani - double dX0 = ( i + 0.5) * dStep ; - double dX1 = ( i + 1.5) * dStep ; - double dY0 = ( j + 0.5) * dStep ; - double dY1 = ( j + 1.5) * dStep ; - double dZ0 = ( k + 0.5) * dStep ; - double dZ1 = ( k + 1.5) * dStep ; - - - // Analisi faccia 1 - int nCount = 0 ; - double dDelta = abs( CompoVert[0].ptInt.x - dX0) ; - - while ( dDelta < EPS_SMALL && nCount < nLimArrey) { - - ++ nCount ; - dDelta = abs( CompoVert[nCount].ptInt.x - dX0) ; - } - - if ( nCount == nLimArrey) - return false ; - - // Analisi faccia 2 - nCount = 0 ; - dDelta = abs( CompoVert[0].ptInt.x - dX1) ; - - while ( dDelta < EPS_SMALL && nCount < nLimArrey) { - - ++ nCount ; - dDelta = abs( CompoVert[nCount].ptInt.x - dX1) ; - } - - if ( nCount == nLimArrey) - return false ; - - // Analisi faccia 3 - nCount = 0 ; - dDelta = abs( CompoVert[0].ptInt.y - dY0) ; - - while ( dDelta < EPS_SMALL && nCount < nLimArrey) { - - ++ nCount ; - dDelta = abs( CompoVert[nCount].ptInt.y - dY0) ; - } - - if ( nCount == nLimArrey) - return false ; - - // Analisi faccia 4 - nCount = 0 ; - dDelta = abs( CompoVert[0].ptInt.y - dY1) ; - - while ( dDelta < EPS_SMALL && nCount < nLimArrey) { - - ++ nCount ; - dDelta = abs( CompoVert[nCount].ptInt.y - dY1) ; - } - - if ( nCount == nLimArrey) - return false ; - - // Analisi faccia 5 - nCount = 0 ; - dDelta = abs( CompoVert[0].ptInt.z - dZ0) ; - - while ( dDelta < EPS_SMALL && nCount < nLimArrey) { - - ++ nCount ; - dDelta = abs( CompoVert[nCount].ptInt.z - dZ0) ; - } - - if ( nCount == nLimArrey) - return false ; - - // Analisi faccia 6 - nCount = 0 ; - dDelta = abs( CompoVert[0].ptInt.z - dZ1) ; - - while ( dDelta < EPS_SMALL && nCount < nLimArrey) { - - ++ nCount ; - dDelta = abs( CompoVert[nCount].ptInt.z - dZ1) ; - } - - if ( nCount == nLimArrey) - return false ; - - return true ; -} // ------------------------- VISUALIZZAZIONE -------------------------------------------------------------------------------------- @@ -304,26 +271,6 @@ VolZmap::GetAllTriangles( TRIA3DLIST& lstTria) const } } } - //else { - // - // //std::vector vecTria ; - // //vecTria.resize( int( m_BlockToUpdate.size())) ; - - // //for ( int i = 0 ; i < int( m_BlockToUpdate.size()) ; ++ i) { - // - // //if ( m_BlockToUpdate[i]) - // // ExtMarchingCubes( i, lstTria, vecTria[i]) ; } - - // TriHolder triHold ; - - // ExtMarchingCubes( 0, lstTria, triHold) ; - // FlipEdges( triHold) ; - // - // for ( int i = 0 ; i < int( triHold.size()) ; ++ i) - // for ( int j = 0 ; j < int( triHold[i].vecTria.size()) ; ++ j) - - // lstTria.emplace_back( triHold[i].vecTria[j]) ; - //} else MarchingCubes( lstTria) ; @@ -332,125 +279,105 @@ VolZmap::GetAllTriangles( TRIA3DLIST& lstTria) const //---------------------------------------------------------------------------- bool -VolZmap::GetBlockTriangles( int nBlock, TRIA3DLIST& lstTria) const +VolZmap::GetTriangles( bool bAllBlocks, INTVECTOR& nModifiedBlocks, TRIA3DLISTVECTOR& vLstTria) const { - // Controllo sulla validità del blocco - if ( nBlock < 0 || nBlock >= int( m_BlockToUpdate.size())) - return false ; - // Caso di singola mappa if ( m_nMapNum == 1) { const int MAX_DIM_CHUNK = 128 ; - // Calcolo posizione del blocco nella griglia - int nIBlock = nBlock % int( m_nFracLin[0]) ; - int nJBlock = nBlock / int( m_nFracLin[0]) ; + nModifiedBlocks.resize( m_nNumBlock) ; + vLstTria.reserve( m_nNumBlock) ; + + // Ciclo sui blocchi + for ( size_t t = 0 ; t < m_nNumBlock ; ++ t) { - // Calcolo limiti per l'indice i - int nStartI = nIBlock * int( m_nDexNumPBlock) ; - int nEndI = ( nIBlock + 1 == int( m_nFracLin[0]) ? - int( m_nNx[0]) : ( nIBlock + 1) * int( m_nDexNumPBlock)) ; + // Se il blocco deve essere aggiornato, eseguo + if ( bAllBlocks || m_BlockToUpdate[t]) { - // Calcolo limiti per l'indice j - int nStartJ = nJBlock * int( m_nDexNumPBlock) ; - int nEndJ = ( nJBlock + 1 == int( m_nFracLin[1]) ? - int( m_nNy[0]) : ( nJBlock + 1) * int( m_nDexNumPBlock)) ; + // preparo lista + vLstTria.emplace_back() ; + nModifiedBlocks[t] = int( vLstTria.size()) - 1 ; + + // Calcolo posizione del blocco nella griglia + int nIBlock = int( t) % int( m_nFracLin[0]) ; + int nJBlock = int( t) / int( m_nFracLin[0]) ; - // Ciclo su i e j - for ( int i = nStartI ; i < nEndI ; i += MAX_DIM_CHUNK) { - int nDimChunkX = min( MAX_DIM_CHUNK, nEndI - i) ; - for ( int j = nStartJ ; j < nEndJ ; j += MAX_DIM_CHUNK) { - int nDimChunkY = min( MAX_DIM_CHUNK, nEndJ - j) ; - GetChunkPrisms( i, j, nDimChunkX, nDimChunkY, MAX_DIM_CHUNK, lstTria) ; - } + // Calcolo limiti per l'indice i + int nStartI = nIBlock * int( m_nDexNumPBlock) ; + int nEndI = ( nIBlock + 1 == int( m_nFracLin[0]) ? + int( m_nNx[0]) : ( nIBlock + 1) * int( m_nDexNumPBlock)) ; + + // Calcolo limiti per l'indice j + int nStartJ = nJBlock * int( m_nDexNumPBlock) ; + int nEndJ = ( nJBlock + 1 == int( m_nFracLin[1]) ? + int( m_nNy[0]) : ( nJBlock + 1) * int( m_nDexNumPBlock)) ; + + // Ciclo su i e j + for ( int i = nStartI ; i < nEndI ; i += MAX_DIM_CHUNK) { + int nDimChunkX = min( MAX_DIM_CHUNK, nEndI - i) ; + for ( int j = nStartJ ; j < nEndJ ; j += MAX_DIM_CHUNK) { + int nDimChunkY = min( MAX_DIM_CHUNK, nEndJ - j) ; + GetChunkPrisms( i, j, nDimChunkX, nDimChunkY, MAX_DIM_CHUNK, vLstTria.back()) ; + } + } + + m_BlockToUpdate[t] = false ; + } + else + nModifiedBlocks[t] = -1 ; } } + // Caso con tre mappe else { - // Calcolo i limiti sugli indici dei voxel del blocco - // Vettore indici i,j,k del blocco - int nIJK[3] ; - GetBlockIJK( nIJK, nBlock) ; - // Vettore limiti sugli indici dei voxel del blocco - int LimitIndexes[6] ; - GetBlockLimitsIJK( nIJK, LimitIndexes) ; + nModifiedBlocks.resize( m_nNumBlock) ; + vLstTria.reserve( m_nNumBlock) ; - TriHolder triHold ; + std::vector VecTriHold ; + VecTriHold.resize( m_nNumBlock) ; - ExtMarchingCubes( LimitIndexes, lstTria, triHold) ; - FlipEdges( triHold) ; - - // Valuto se esistono voxel contenenti sharp feature - if ( ! triHold.size()) + // Ciclo sui blocchi + for ( size_t t = 0 ; t < m_nNumBlock ; ++ t) { + // Calcolo i limiti sugli indici dei voxel del blocco + int nIJK[3] ; + // Vettore indici i,j,k del blocco + GetBlockIJK( int( t), nIJK) ; + // Vettore limiti sugli indici dei voxel del blocco + int LimitIndexes[6] ; + GetBlockLimitsIJK( nIJK, LimitIndexes) ; + // Se il blocco deve essere processato, eseguo + if ( bAllBlocks || m_BlockToUpdate[t]) { + vLstTria.emplace_back() ; + nModifiedBlocks[t] = int( vLstTria.size()) - 1 ; + ExtMarchingCubes( LimitIndexes, vLstTria.back(), VecTriHold[t]) ; + m_BlockToUpdate[t] = false ; + } + else + nModifiedBlocks[t] = -1 ; + } - ; - - else { + // Eseguo il flipping + FlipEdges( VecTriHold) ; - // Cerco fra i voxel con sharp feature quelli - // di frontiera e quelli ad essi adiacenti. - - // Ciclo sui voxel con sharp feature - for ( size_t t = 0 ; t < triHold.size() ; ++ t) { - - // Voxel di frontiera - if ( triHold[t].i == LimitIndexes[0] || triHold[t].i == LimitIndexes[1] - 1 || - triHold[t].j == LimitIndexes[2] || triHold[t].j == LimitIndexes[3] - 1 || - triHold[t].k == LimitIndexes[4] || triHold[t].k == LimitIndexes[5] - 1) { - - // Ciclo sui voxel adiacenti - for ( int nI = triHold[t].i - 1 ; nI <= triHold[t].i + 1 ; ++ nI) { - for ( int nJ = triHold[t].j - 1 ; nJ <= triHold[t].j + 1 ; ++ nJ) { - for ( int nK = triHold[t].k - 1 ; nK <= triHold[t].k + 1 ; ++ nK) { - - // Voxel adiacente appartenente allo - // Zmap e non appartenente al blocco - if ( ( nI >= -1 && nI <= int( m_nNx[0]) - 1 && - nJ >= -1 && nJ <= int( m_nNy[0]) - 1 && - nK >= -1 && nK <= int( m_nNy[1]) - 1) && - ! ( nI >= LimitIndexes[0] && nI < LimitIndexes[1] && - nJ >= LimitIndexes[2] && nJ < LimitIndexes[3] && - nK >= LimitIndexes[4] && nK < LimitIndexes[5])) { - - std::vector AdjVoxel ; - - AdjVoxel.emplace_back( nI) ; - AdjVoxel.emplace_back( nJ) ; - AdjVoxel.emplace_back( nK) ; - - TriHolder NewTriHold ; - - ExtMarchingCubes( AdjVoxel, NewTriHold) ; - - if ( NewTriHold.size()) - FlipEdgesLocalFlipEdges( triHold[t], NewTriHold[0]) ; - } - } - } + // ciclo sui blocchi + for ( size_t t = 0 ; t < m_nNumBlock ; ++ t) { + // ciclo sui voxel del blocco + for ( size_t t1 = 0 ; t1 < VecTriHold[t].size() ; ++ t1) { + // ciclo sulle componenti connesse del voxel + for ( size_t t2 = 0 ; t2 < VecTriHold[t][t1].vCompoTria.size() ; ++ t2) { + // ciclo sui triangoli delle componenti connesse + for ( size_t t3 = 0 ; t3 < VecTriHold[t][t1].vCompoTria[t2].size() ; ++ t3) { + // aggiungo un singolo triangolo alla lista + if ( nModifiedBlocks[t] >= 0) + vLstTria[nModifiedBlocks[t]].emplace_back( VecTriHold[t][t1].vCompoTria[t2][t3]) ; } } - } - } - - for ( size_t t1 = 0 ; t1 < triHold.size() ; ++ t1) - for ( size_t t2 = 0 ; t2 < triHold[t1].vCompoTria.size() ; ++ t2) - for ( size_t t3 = 0 ; t3 < triHold[t1].vCompoTria[t2].size() ; ++ t3) - lstTria.emplace_back( triHold[t1].vCompoTria[t2][t3]) ; - + } + } } - m_BlockToUpdate[nBlock] = false ; - - return true ; -} - -//---------------------------------------------------------------------------- -bool -VolZmap::GetBlockInfo( std::vector & bModified) const -{ - bModified = m_BlockToUpdate ; return true ; } @@ -874,6 +801,7 @@ VolZmap::ExtMarchingCubes( const int nLimits[], TRIA3DLIST& lstTria, TriHolder& for ( int j = nLimits[2] ; j < nLimits[3] ; ++ j) { for ( int k = nLimits[4] ; k < nLimits[5] ; ++ k) { + // Indici i,j,k dei vertici int IndexCorner[8][3] = { { i, j, k}, @@ -919,7 +847,12 @@ VolZmap::ExtMarchingCubes( const int nLimits[], TRIA3DLIST& lstTria, TriHolder& // Arrey di strutture punto di intersezione // e normale alla superficie in esso. VectorField VecField[12] ; - + + // Flag sulla regolatrità dei campi scalare e vettoriale: + // se i campi sono regolari esso resta vero, altrimenti + // assume il valore falso. + bool bReg = true ; + // Ciclo sui segmenti for ( int EdgeIndex = 0 ; EdgeIndex < 12 ; ++ EdgeIndex) { // Se il segmento non attraversa la superficie passo al successivo @@ -931,16 +864,21 @@ VolZmap::ExtMarchingCubes( const int nLimits[], TRIA3DLIST& lstTria, TriHolder& bool bN1 = ( ( nIndex & ( 1 << n1)) != 0) ; - // Determino con precisione il punto di intersezione sullo spigolo. - NewIntersPos( IndexCorner[n1], IndexCorner[n2], bN1, - VecField[EdgeIndex].ptInt, - VecField[EdgeIndex].vtNorm) ; - + // Determino con precisione il punto di intersezione sullo spigolo, + // se i campi scalare e vettoriale non sono regolari bReg diviene falso. + if ( ! IntersPos( IndexCorner[n1], IndexCorner[n2], bN1, + VecField[EdgeIndex].ptInt, + VecField[EdgeIndex].vtNorm)) + + bReg = false ; + + // Riporto punti e normali nel sistema locale in cui + // è immerso lo Zmap col suo sistema di riferimento. VecField[EdgeIndex].ptInt.ToGlob( m_MapFrame[0]) ; VecField[EdgeIndex].vtNorm.ToGlob( m_MapFrame[0]) ; } - - // Determino il numero di componenti connesse + + // Determino il numero di componenti connesse nel voxel int nComponents = TriangleTableEn[nIndex][1][0] ; // Serve nel ciclo che salva i punti e vettori di @@ -951,7 +889,8 @@ VolZmap::ExtMarchingCubes( const int nLimits[], TRIA3DLIST& lstTria, TriHolder& // della seconda e così via. int nTableOffset = nComponents ; - // Numero di feature nel voxel: al più vi è una feature per componente connessa + // Numero di feature nel voxel: al più vi è + // una feature per componente connessa. int nFeatureInVoxel = 0 ; // Ciclo sulle componenti @@ -972,19 +911,33 @@ VolZmap::ExtMarchingCubes( const int nLimits[], TRIA3DLIST& lstTria, TriHolder& // Valuto le relazioni reciproche fra le normali e // se i punti sono su un piano di griglia. - bool bNormal = TestOnNormal( CompoVert, nVertComp, nFeatureType) ; - bool bExt = bNormal ; + Vector3d vtFeatureAxis ; + bool bNormal = false ; + + // Indici dei punti ove le normali + // formano il massimo angolo + int nMin1, nMin2 ; + + // Se i campi sono regolari valuto le normali + // per stabilire se eseguire ExtMC o MC, + // altrimenti eseguo MC. + if ( bReg) + + bNormal = TestOnNormal( CompoVert, nVertComp, nFeatureType, vtFeatureAxis, nMin1, nMin2) ; + + // Flag ExtMC + bool bExtMC = bNormal ; // Extended MC - if ( bExt) { + if ( bExtMC) { // Passo al sistema di riferimento del baricentro Point3d ptGravityCenter( 0, 0, 0) ; for ( int ni = 0 ; ni < nVertComp ; ++ ni) - ptGravityCenter += CompoVert[ni].ptInt ; + ptGravityCenter = ptGravityCenter + CompoVert[ni].ptInt ; - ptGravityCenter /= nVertComp ; + ptGravityCenter = ptGravityCenter / nVertComp ; Vector3d vtO = ptGravityCenter - ORIG ; @@ -1004,7 +957,7 @@ VolZmap::ExtMarchingCubes( const int nLimits[], TRIA3DLIST& lstTria, TriHolder& dMatrixN.resize( nVertComp, 3) ; dKnownVector.resize( nVertComp, 1) ; dUnknownVector.resize( 3, 1) ; - + // Studio del caso 4 punti su un piano int nEqual = 0 ; int nPosD ; @@ -1043,7 +996,7 @@ VolZmap::ExtMarchingCubes( const int nLimits[], TRIA3DLIST& lstTria, TriHolder& ( ( CompoVert[2].ptInt - CompoVert[1].ptInt) ^ ( CompoVert[3].ptInt - CompoVert[2].ptInt))) ; - // Caso quattro punti su un piano + // Caso superficie piana if ( nVertComp == 4 && nEqual == 2 && dDot < EPS_SMALL) { for ( int ni = 0 ; ni < nVertComp ; ++ ni) { @@ -1061,7 +1014,7 @@ VolZmap::ExtMarchingCubes( const int nLimits[], TRIA3DLIST& lstTria, TriHolder& } } } - // caso generale + // Caso generale else { for ( int ni = 0 ; ni < nVertComp ; ++ ni) { dMatrixN( ni, 0) = CompoVert[ni].vtNorm.x ; @@ -1077,13 +1030,12 @@ VolZmap::ExtMarchingCubes( const int nLimits[], TRIA3DLIST& lstTria, TriHolder& dMatrixV = svd.matrixV() ; dSingularValue = svd.singularValues() ; - double s0 = dSingularValue( 0) ; - double s1 = dSingularValue( 1) ; - double s2 = dSingularValue( 2) ; - - if ( nFeatureType == 2 ) - dSingularValue( 2) = 0 ; + // Se la feature è un edge annullo + // il valore singolare minore. + if ( nFeatureType == Edge) + dSingularValue( 2) = 0 ; + // Back substitution: risolvo il sistema USV*x = b // Calcolo U^T b double vTemp[3] ; @@ -1092,7 +1044,7 @@ VolZmap::ExtMarchingCubes( const int nLimits[], TRIA3DLIST& lstTria, TriHolder& double s = 0 ; - if ( abs( dSingularValue( ni)) > EPS_SMALL) { + if ( dSingularValue( ni) > 0) { for ( int nj = 0 ; nj < nVertComp ; ++ nj) s += dMatrixU( nj, ni) * dKnownVector( nj) ; @@ -1120,73 +1072,212 @@ VolZmap::ExtMarchingCubes( const int nLimits[], TRIA3DLIST& lstTria, TriHolder& dUnknownVector( 2)) ; double dDistFeature = vtFeature.Len() ; - const double dMaxDist = sqrt( 3) * m_dStep ; - if ( dDistFeature > dMaxDist) { - - // Costruzione dei triangoli - for ( int TriIndex = 0 ; TriangleTableEn[nIndex][0][TriIndex] != - 1 ; TriIndex += 3) { + // Flag sulla distanza del vertice dal + // baricentro del sistema di punti + bool bOutside = false ; - // Costruzione triangolo - int i0 = TriangleTableEn[nIndex][0][TriIndex + 2] ; - int i1 = TriangleTableEn[nIndex][0][TriIndex + 1] ; - int i2 = TriangleTableEn[nIndex][0][TriIndex] ; - - Triangle3d CurrentTriangle ; - - // Il triangolo è pronto - CurrentTriangle.Set( VecField[i0].ptInt, VecField[i1].ptInt, VecField[i2].ptInt) ; - CurrentTriangle.Validate( true) ; - - // Se il triangolo non è degenere lo aggiungo alla lista - if ( ! ( AreSamePointApprox( CurrentTriangle.GetP( 0), CurrentTriangle.GetP( 1)) || - AreSamePointApprox( CurrentTriangle.GetP( 1), CurrentTriangle.GetP( 2)) || - AreSamePointApprox( CurrentTriangle.GetP( 2), CurrentTriangle.GetP( 0)))) - lstTria.emplace_back( CurrentTriangle) ; - } - - continue ; - } - + if ( dDistFeature > dMaxDist) + bOutside = true ; + // Esprimo la soluzione nel sistema di riferimento dello z-map. Point3d ptSol( dUnknownVector( 0) + vtO.x, dUnknownVector( 1) + vtO.y, dUnknownVector( 2) + vtO.z) ; + Triangle3d CurrentTriangle ; TRIA3DVECTOR triContainer ; - bool bInvNormal = false ; + // Flag sull'inversione delle normali + bool bInvNormB = false ; - for ( int ni = 0 ; ni < nVertComp ; ++ ni) { + // Questo controllo viene eseguito solo se + // il vertice è distante dal baricentro + // entro la soglia stabilita. + if ( ! bOutside) { + for ( int ni = 0 ; ni < nVertComp ; ++ ni) { + + int nj = ( ni + 1 < nVertComp) ? ni + 1 : 0 ; + // Il triangolo è pronto + CurrentTriangle.Set( ptSol, CompoVert[nj].ptInt, CompoVert[ni].ptInt) ; + CurrentTriangle.Validate( true) ; + + // Aggiungo triangolo al vettore temporaneo + triContainer.emplace_back( CurrentTriangle) ; - // Il triangolo è pronto - int nj = ( ni + 1 < nVertComp) ? ni + 1 : 0 ; - CurrentTriangle.Set( ptSol, CompoVert[nj].ptInt, CompoVert[ni].ptInt) ; - CurrentTriangle.Validate( true) ; - - // Test sulla degenerazione - if ( CurrentTriangle.GetArea() < EPS_SMALL) - continue ; - - // Controllo sull'inversione delle normali - if ( CurrentTriangle.GetN() * CompoVert[nj].vtNorm < - 0.7 || - CurrentTriangle.GetN() * CompoVert[ni].vtNorm < - 0.7) { - //string sInfo ; - //sInfo = "TriaN=" + ToString( CurrentTriangle.GetN()) + - // "VertJ=" + ToString( CompoVert[nj].vtNorm) + - // "VertI=" + ToString( CompoVert[ni].vtNorm) ; - //LOG_INFO( GetEGkLogger(), sInfo.c_str()) - //bInvNormal = true ; + // Controllo sull'inversione delle normali + if ( ( CurrentTriangle.GetN() * CompoVert[nj].vtNorm < - 0.01 && + CurrentTriangle.GetN() * CompoVert[ni].vtNorm < - 0.01) && + ( ! bInvNormB)) { + + ptSol = ptGravityCenter ; + bInvNormB = true ; + triContainer.resize( 0) ; + ni = -1 ; + } } - - // Aggiungo triangolo al vettore temporaneo - triContainer.emplace_back( CurrentTriangle) ; } - // Valuto normali - int nContSize = triContainer.size() ; + // Questo flag esprime se il vertice è, entro la tolleranza, + // interno o esterno al voxel a cui appartiene. + bool bInsideVoxel = IsPointInsideVoxelApprox( i, j, k, ptSol) ; + + // Proprietà gemoetriche locali + // del campo vettoriale: se vero + // procediamo con un ulteriore + // analisi per eliminare triangoli + // invertiti. + bool bLocProp = false ; + + if ( nVertComp == 4) { + + int nNegDotNum = 0 ; + + for ( int nLocIndI = 0 ; nLocIndI < 3 ; ++ nLocIndI) { + for ( int nLocIndJ = nLocIndI + 1 ; nLocIndJ < 4 ; ++ nLocIndJ) { + + if ( CompoVert[nLocIndI].vtNorm * CompoVert[nLocIndJ].vtNorm < 0) { + + nNegDotNum ++ ; + } + } + } + + if ( nNegDotNum == 3) + + bLocProp = true ; + } + + // Se è necessario si effettua l'ulteriore analisi + // Questo controllo si effettua se la feature non + // esce dai limiti, non esce dal voxel e il primo + // controllo sull'inversione delle normali ha dato + // esito negativo. + if ( nFeatureType == 2 && bLocProp && + ! ( bOutside || bInsideVoxel || bInvNormB)) { + + if ( abs( nMin1 - nMin2) == 1 || + abs( nMin1 - nMin2) == 3) { + + int nSum = nMin1 + nMin2 ; + int nTriNum = ( nSum == 3 && ( nMin1 == 3 || nMin2 == 3) ? + max( nMin1, nMin2) : min( nMin1, nMin2)) ; + + double dDot1 = triContainer[nTriNum].GetN() * CompoVert[nMin1].vtNorm ; + double dDot2 = triContainer[nTriNum].GetN() * CompoVert[nMin2].vtNorm ; + + if ( ( dDot1 < - 0.2 && dDot2 > - EPS_ZERO) || + ( dDot2 < - 0.2 && dDot1 > - EPS_ZERO)) { + + int nNm = dDot1 < - 0.2 ? nMin1 : nMin2 ; + int nNp = dDot1 < - 0.2 ? nMin2 : nMin1 ; + + + Vector3d vtVV = CompoVert[nNp].ptInt - CompoVert[nNm].ptInt ; + Point3d ptSolZmFrame = ptSol ; + + double dLenVV = vtVV.Len() ; + + ptSolZmFrame.ToLoc( m_MapFrame[0]) ; + + if ( ! IsPointInsideVoxel( i, j, k, ptSol)) { + + Vector3d vtVS = ptSol - CompoVert[nNm].ptInt ; + + vtVV.Normalize() ; + + Vector3d vtVSNew = ( vtVS * vtVV) * vtVV ; + + ptSol = CompoVert[nNm].ptInt + vtVSNew ; + + double dLNm = ( ptSol - CompoVert[nNm].ptInt).Len() ; + double dLNp = ( ptSol - CompoVert[nNp].ptInt).Len() ; + + if ( dLNm > dLenVV || dLNp > dLenVV) { + + ptSol = dLNm < dLNp ? CompoVert[nNm].ptInt : CompoVert[nNp].ptInt ; + } + + triContainer.resize( 0) ; + + for ( int nc = 0 ; nc < nVertComp ; ++ nc) { + + int nd = ( nc + 1 < nVertComp) ? nc + 1 : 0 ; + // Il triangolo è pronto + CurrentTriangle.Set( ptSol, CompoVert[nd].ptInt, CompoVert[nc].ptInt) ; + CurrentTriangle.Validate( true) ; + + // Aggiungo triangolo al vettore temporaneo + triContainer.emplace_back( CurrentTriangle) ; + } + } + } + } + else { + + int nPrev1 = nMin1 == 0 ? 3 : nMin1 - 1 ; + int nPrev2 = nMin2 == 0 ? 3 : nMin2 - 1 ; + int nNext1 = nMin1 == 3 ? 0 : nMin1 + 1 ; + int nNext2 = nMin2 == 3 ? 0 : nMin2 + 1 ; + + int nNeighbourIndex ; + int nStartIndex ; + + if ( CompoVert[nPrev1].vtNorm * CompoVert[nMin1].vtNorm > 0) { + + bool bTestOnVert = abs( nPrev1 - nMin2) == 0 || abs( nPrev1 - nMin2) == 3 ; + nNeighbourIndex = bTestOnVert ? nPrev1 : nNext1 ; + nStartIndex = nMin2 ; + } + else { + + bool bTestOnVert = abs( nPrev2 - nMin1) == 0 || abs( nPrev2 - nMin1) == 3 ; + nNeighbourIndex = bTestOnVert ? nPrev2 : nNext2 ; + nStartIndex = nMin1 ; + } + + Vector3d vtVV = CompoVert[nNeighbourIndex].ptInt - CompoVert[nStartIndex].ptInt ; + double dVVLen = vtVV.Len() ; + vtVV.Normalize() ; + + Vector3d vtVF = ptSol - CompoVert[nStartIndex].ptInt ; + Vector3d vtNewVF = ( vtVF * vtVV) * vtVV ; + + double dNewVFLen = vtNewVF.Len() ; + + if ( dNewVFLen > dVVLen) { + + ptSol = CompoVert[nNeighbourIndex].ptInt ; + } + else { + + ptSol = CompoVert[nStartIndex].ptInt + vtNewVF ; + } + + + triContainer.resize( 0) ; + + for ( int nc = 0 ; nc < nVertComp ; ++ nc) { + + int nd = ( nc + 1 < nVertComp) ? nc + 1 : 0 ; + // Il triangolo è pronto + CurrentTriangle.Set( ptSol, CompoVert[nd].ptInt, CompoVert[nc].ptInt) ; + CurrentTriangle.Validate( true) ; + + // Aggiungo triangolo al vettore temporaneo + triContainer.emplace_back( CurrentTriangle) ; + } + } + } + + + // Valuto normali: questo è ancora un controllo + // sulle normali, se risultano in tutti i punti + // approssimativamente uguali passiamo alla + // routine standard + int nContSize = int( triContainer.size()) ; bool bPlane = true ; @@ -1201,9 +1292,15 @@ VolZmap::ExtMarchingCubes( const int nLimits[], TRIA3DLIST& lstTria, TriHolder& break ; } } + + if ( ! bPlane) + break ; } - if ( ! ( bPlane || bInvNormal)) { + // Se la feaure non è fuori dai limiti + // e i triangoli formano una superficie + // non piana confermo ExtMC + if ( ! ( bOutside || bPlane)) { // Aggiorno il numero di feature. ++ nFeatureInVoxel ; @@ -1216,7 +1313,7 @@ VolZmap::ExtMarchingCubes( const int nLimits[], TRIA3DLIST& lstTria, TriHolder& triHold.resize( triHold.size() + 1) ; // Questi dati dipendono solo dal voxel, - // quindi sono validi validi per tutte le + // quindi sono validi per tutte le // componenti che vi appartengono. int nCurrent = int( triHold.size()) - 1 ; @@ -1234,19 +1331,24 @@ VolZmap::ExtMarchingCubes( const int nLimits[], TRIA3DLIST& lstTria, TriHolder& // connessa alla lista dei vertici. triHold[nCurrent].ptCompoVert.emplace_back( ptSol) ; - int nOldFeatureNum = triHold[nCurrent].vCompoTria.size() ; + int nOldFeatureNum = int( triHold[nCurrent].vCompoTria.size()) ; int nNewFeatureNum = nOldFeatureNum + 1 ; // Aggiungo una componente per il vettore // dei triangoli della componente connessa. triHold[nCurrent].vCompoTria.resize( nNewFeatureNum) ; - for ( int ni = 0 ; ni < nContSize ; ++ ni) + for ( int ni = 0 ; ni < nContSize ; ++ ni) { + // Se l'area è non nulla aggiungo il triangolo + if ( triContainer[ni].GetArea() > EPS_SMALL) { + triHold[nCurrent].vCompoTria[nOldFeatureNum].emplace_back( triContainer[ni]) ; + } + } } - - else { + // ExtMC non confermato, si passa a MC + else { // Costruzione dei triangoli for ( int TriIndex = 0 ; TriangleTableEn[nIndex][0][TriIndex] != - 1 ; TriIndex += 3) { @@ -1262,15 +1364,13 @@ VolZmap::ExtMarchingCubes( const int nLimits[], TRIA3DLIST& lstTria, TriHolder& CurrentTriangle.Validate( true) ; // Se il triangolo non è degenere lo aggiungo alla lista - if ( ! ( AreSamePointApprox( CurrentTriangle.GetP( 0), CurrentTriangle.GetP( 1)) || - AreSamePointApprox( CurrentTriangle.GetP( 1), CurrentTriangle.GetP( 2)) || - AreSamePointApprox( CurrentTriangle.GetP( 2), CurrentTriangle.GetP( 0)))) + if ( CurrentTriangle.GetArea() > EPS_SMALL) lstTria.emplace_back( CurrentTriangle) ; } } } // Standard MC - else { + else { // Costruzione dei triangoli for ( int TriIndex = 0 ; TriangleTableEn[nIndex][0][TriIndex] != - 1 ; TriIndex += 3) { @@ -1287,9 +1387,7 @@ VolZmap::ExtMarchingCubes( const int nLimits[], TRIA3DLIST& lstTria, TriHolder& CurrentTriangle.Validate( true) ; // Se il triangolo non è degenere lo aggiungo alla lista - if ( ! ( AreSamePointApprox( CurrentTriangle.GetP( 0), CurrentTriangle.GetP( 1)) || - AreSamePointApprox( CurrentTriangle.GetP( 1), CurrentTriangle.GetP( 2)) || - AreSamePointApprox( CurrentTriangle.GetP( 2), CurrentTriangle.GetP( 0)))) + if ( CurrentTriangle.GetArea() > EPS_SMALL) lstTria.emplace_back( CurrentTriangle) ; } } @@ -1303,608 +1401,145 @@ VolZmap::ExtMarchingCubes( const int nLimits[], TRIA3DLIST& lstTria, TriHolder& return true ; } -//---------------------------------------------------------------------------- -bool -VolZmap::ExtMarchingCubes( const std::vector VoxelsIndexes, TriHolder& triHold) const -{ - // Controllo sulla validità della dimensione del - // vector: essa deve essere un multiplo di 3 - if ( ( VoxelsIndexes.size() % 3)) - return false ; - - Point3d ptMapOrig = m_MapFrame[0].Orig() ; - - // Ciclo sui voxel - for ( size_t t = 0 ; t < VoxelsIndexes.size() - 2 ; t += 3) { - - // Indici del voxel - int i = VoxelsIndexes[t] ; - int j = VoxelsIndexes[t + 1] ; - int k = VoxelsIndexes[t + 2] ; - - // Indici i,j,k dei vertici - int IndexCorner[8][3] = { - { i, j, k}, - { i + 1, j, k}, - { i + 1, j + 1, k}, - { i, j + 1, k}, - { i, j, k + 1}, - { i + 1, j, k + 1}, - { i + 1, j + 1, k + 1}, - { i, j + 1, k + 1} - } ; - - int nIndex = 0 ; - - // Classificazione dei vertici: interni o esterni al materiale - if ( IsThereMat( i, j, k)) - nIndex |= ( 1 << 0) ; - if ( IsThereMat( i + 1, j, k)) - nIndex |= ( 1 << 1) ; - if ( IsThereMat( i + 1, j + 1, k)) - nIndex |= ( 1 << 2) ; - if ( IsThereMat( i, j + 1, k)) - nIndex |= ( 1 << 3) ; - if ( IsThereMat( i, j, k + 1)) - nIndex |= ( 1 << 4) ; - if ( IsThereMat( i + 1, j, k + 1)) - nIndex |= ( 1 << 5) ; - if ( IsThereMat( i + 1, j + 1, k + 1)) - nIndex |= ( 1 << 6) ; - if ( IsThereMat( i, j + 1, k + 1)) - nIndex |= ( 1 << 7) ; - - // Se vi è qualche intersezione fra segmenti e superficie - // continuo altrimenti passo al prossimo voxel. - if ( EdgeTable[nIndex] == 0) - continue ; - - static int intersections[12][2] = { - { 0, 1 }, { 1, 2 }, { 2, 3 }, { 0, 3 }, { 4, 5 }, { 5, 6 }, - { 6, 7 }, { 4, 7 }, { 0, 4 }, { 1, 5 }, { 2, 6 }, { 3, 7 } - } ; - - // Arrey di strutture punto di intersezione - // e normale alla superficie in esso. - VectorField VecField[12] ; - - // Ciclo sui segmenti - for ( int EdgeIndex = 0 ; EdgeIndex < 12 ; ++ EdgeIndex) { - // Se il segmento non attraversa la superficie passo al successivo - if ( ! ( EdgeTable[nIndex] & ( 1 << EdgeIndex))) - continue ; - - int n1 = intersections[EdgeIndex][0] ; - int n2 = intersections[EdgeIndex][1] ; - - // Determino con precisione il punto di intersezione sullo spigolo. - IntersPos( IndexCorner[n1], IndexCorner[n2], - VecField[EdgeIndex].ptInt, - VecField[EdgeIndex].vtNorm) ; - - VecField[EdgeIndex].ptInt.ToGlob( m_MapFrame[0]) ; - VecField[EdgeIndex].vtNorm.ToGlob( m_MapFrame[0]) ; - } - - // Determino il numero di componenti connesse - int nComponents = TriangleTableEn[nIndex][1][0] ; - - // Serve nel ciclo che salva i punti e vettori di - // una componente nell'arrey di compentenza: La tabella - // fornisce numero di componenti, numero di vertici per - // componenti per OGNUNA delle componenti e in fine - // elenca i vertici della prima componente, seguiti da quelli - // della seconda e così via. - int nTableOffset = nComponents ; - - // Numero di feature nel voxel: al più vi - // è una feature per componente connessa - int nFeatureInVoxel = 0 ; - - // Ciclo sulle componenti - for ( int nCompCount = 1 ; nCompCount <= nComponents ; ++ nCompCount) { - - // Numero vertici per componenti - int nVertComp = TriangleTableEn[nIndex][1][nCompCount] ; - - // Vettore di Vector3d - VectorField CompoVert[12] ; - - // Riempio il vettore - for ( int nVertCount = 0 ; nVertCount < nVertComp ; ++ nVertCount) - // Nota che il primo elemento dell'array (0-esimo) non viene inizializzato - CompoVert[nVertCount] = VecField[TriangleTableEn[nIndex][1][nVertCount + nTableOffset + 1]] ; - - int nFeatureType ; - - // Valuto le relazioni reciproche fra le normali - bool bExt = TestOnNormal( CompoVert, nVertComp, nFeatureType) ; - - // Extended MC - if ( bExt) { - - // Aggiorno il numero di feature. - ++ nFeatureInVoxel ; - - // Se siamo in presenza della prima feature del - // voxel, ridimensiono il vettore che contiene - // la struttura dati del voxel. - if ( nFeatureInVoxel == 1) { - - triHold.resize( triHold.size() + 1) ; - - // Questi dati dipendono solo dal voxel, - // quindi sono validi validi per tutte le - // componenti che vi appartengono. - int nCurrent = int( triHold.size()) - 1 ; - - triHold[nCurrent].i = i ; - triHold[nCurrent].j = j ; - triHold[nCurrent].k = k ; - } - - // Indice che identifica la posizione del voxel - // nel vector. - int nCurrent = int( triHold.size()) - 1 ; - - // Passo al sistema di riferimento del baricentro - Point3d ptGravityCenter( 0, 0, 0) ; - - for ( int ni = 0 ; ni < nVertComp ; ++ ni) - - ptGravityCenter += CompoVert[ni].ptInt ; - - ptGravityCenter /= nVertComp ; - - Vector3d vtO = ptGravityCenter - ORIG ; - - Point3d ptTrasf[12] ; - - for ( int ni = 0 ; ni < nVertComp ; ++ ni) - - ptTrasf[ni] = CompoVert[ni].ptInt - vtO ; - - - // Preparo le matrici per il sistema - typedef Eigen::Matrix dSystemMatrix ; - typedef Eigen::Matrix dSystemVector ; - typedef Eigen::JacobiSVD DecomposerSVD ; - - dSystemMatrix dMatrixN, dMatrixU, dMatrixV ; - dSystemVector dKnownVector, dUnknownVector, dSingularValue ; - - dMatrixN.resize( nVertComp, 3) ; - dKnownVector.resize( nVertComp, 1) ; - dUnknownVector.resize( 3, 1) ; - - // Studio del caso 4 punti su un piano - int nEqual = 0 ; - int nPosD ; - Vector3d vtD, vtE ; - - if ( nVertComp == 4 && nFeatureType == 2) { - - int nPosEq ; - - for ( int ni = 0 ; ni < 2 ; ++ ni) { - - for ( int nj = ni + 1 ; nj < nVertComp ; ++ nj) { - - if ( AreSameVectorApprox( CompoVert[ni].vtNorm, - CompoVert[nj].vtNorm)) { - - nEqual ++ ; - nPosEq = ni ; - } - } - - if ( nEqual == 2) - break ; - else - nEqual = 0 ; - } - - if ( nEqual == 2) { - - for ( int ni = 0 ; ni < nVertComp ; ++ ni) - - if ( ! AreSameVectorApprox( CompoVert[ni].vtNorm, - CompoVert[nPosEq].vtNorm)) { - - nPosD = ni ; - vtD = CompoVert[ni].vtNorm ; - vtE = CompoVert[nPosEq].vtNorm ; - } - } - } - - double dDot = abs( ( CompoVert[1].ptInt - CompoVert[0].ptInt) * - ( ( CompoVert[2].ptInt - CompoVert[1].ptInt) ^ - ( CompoVert[3].ptInt - CompoVert[2].ptInt))) ; - - - - // Caso quattro punti su un piano - if ( nVertComp == 4 && nEqual == 2 && dDot < EPS_SMALL) { - - for ( int ni = 0 ; ni < nVertComp ; ++ ni) { - - if ( ni != nPosD) { - - dMatrixN( ni, 0) = CompoVert[ni].vtNorm.x ; - dMatrixN( ni, 1) = CompoVert[ni].vtNorm.y ; - dMatrixN( ni, 2) = CompoVert[ni].vtNorm.z ; - - dKnownVector( ni) = CompoVert[ni].vtNorm * ( ptTrasf[ni] - ORIG) ; - } - else { - - dMatrixN( ni, 0) = vtE.x ; - dMatrixN( ni, 1) = vtE.y ; - dMatrixN( ni, 2) = vtE.z ; - - dKnownVector( ni) = vtE * ( ptTrasf[ni] - ORIG) ; - } - } - } - // caso generale - else { - - for ( int ni = 0 ; ni < nVertComp ; ++ ni) { - - dMatrixN( ni, 0) = CompoVert[ni].vtNorm.x ; - dMatrixN( ni, 1) = CompoVert[ni].vtNorm.y ; - dMatrixN( ni, 2) = CompoVert[ni].vtNorm.z ; - - dKnownVector( ni) = CompoVert[ni].vtNorm * ( ptTrasf[ni] - ORIG) ; - } - } - - DecomposerSVD svd( dMatrixN, Eigen::ComputeThinU | Eigen::ComputeThinV) ; - - dMatrixU = svd.matrixU() ; - dMatrixV = svd.matrixV() ; - dSingularValue = svd.singularValues() ; - - double s0 = dSingularValue( 0) ; - double s1 = dSingularValue( 1) ; - double s2 = dSingularValue( 2) ; - - if ( nFeatureType == 2 ) - dSingularValue( 2) = 0 ; - - // Back substitution: risolvo il sistema USV*x = b - // Calcolo U^T b - double vTemp[3] ; - - for ( int ni = 0 ; ni < 3 ; ++ ni) { - - double s = 0 ; - - if ( abs( dSingularValue( ni)) > EPS_SMALL) { - - for ( int nj = 0 ; nj < nVertComp ; ++ nj) - - s += dMatrixU( nj, ni) * dKnownVector( nj) ; - - s /= dSingularValue( ni) ; - } - vTemp[ni] = s ; - } - - // Moltiplico per V - for ( int ni = 0 ; ni < 3 ; ++ ni) { - - double s = 0 ; - - for ( int nj = 0 ; nj < 3 ; ++ nj) - s += dMatrixV( ni, nj) * vTemp[nj] ; - - dUnknownVector( ni) = s ; - } - - // Limito la feature entro una distanza di 3 - // volte la diagonale del voxel dal baricentro. - Vector3d vtFeature( dUnknownVector( 0), - dUnknownVector( 1), - dUnknownVector( 2)) ; - - double dDistFeature = vtFeature.Len() ; - - const double dMaxDist = 2 * sqrt( 3) * m_dStep / 2 ; - - if ( dDistFeature > dMaxDist) { - - vtFeature = ( dMaxDist / dDistFeature) * vtFeature ; - - dUnknownVector( 0) = vtFeature.x ; - dUnknownVector( 1) = vtFeature.y ; - dUnknownVector( 2) = vtFeature.z ; - } - - // Esprimo la soluzione nel sistema di riferimento dello z-map - Point3d ptSol( dUnknownVector( 0) + vtO.x, - dUnknownVector( 1) + vtO.y, - dUnknownVector( 2) + vtO.z) ; - - // Aggiungo vertice della componente - // connessa alla lista dei vertici. - triHold[nCurrent].ptCompoVert.emplace_back( ptSol) ; - - int nOldFeatureNum = triHold[nCurrent].vCompoTria.size() ; - int nNewFeatureNum = nOldFeatureNum + 1 ; - - // Aggiungo una componente per il vettore - // dei triangoli della componente connessa. - triHold[nCurrent].vCompoTria.resize( nNewFeatureNum) ; - - Triangle3d CurrentTriangle ; - TRIA3DVECTOR triContainer ; - - for ( int ni = 0 ; ni < nVertComp - 1 ; ++ ni) { - - // Il triangolo è pronto - CurrentTriangle.Set( ptSol, CompoVert[ni+1].ptInt, CompoVert[ni].ptInt) ; - CurrentTriangle.Validate( true) ; - - // Test sulla degenerazione - if ( CurrentTriangle.GetArea() < EPS_SMALL /*AreSamePointApprox( CurrentTriangle.GetP( 0), CurrentTriangle.GetP( 1)) || - AreSamePointApprox( CurrentTriangle.GetP( 1), CurrentTriangle.GetP( 2)) || - AreSamePointApprox( CurrentTriangle.GetP( 2), CurrentTriangle.GetP( 0))*/) - - continue ; - - // Aggiungo triangolo - triHold[nCurrent].vCompoTria[nOldFeatureNum].emplace_back( CurrentTriangle) ; - } - - // Ultimo triangolo - CurrentTriangle.Set( ptSol, CompoVert[0].ptInt, CompoVert[nVertComp - 1].ptInt) ; - CurrentTriangle.Validate( true) ; - - // Test sulla degenerazione - if ( CurrentTriangle.GetArea() < EPS_SMALL /*AreSamePointApprox( CurrentTriangle.GetP( 0), CurrentTriangle.GetP( 1)) || - AreSamePointApprox( CurrentTriangle.GetP( 1), CurrentTriangle.GetP( 2)) || - AreSamePointApprox( CurrentTriangle.GetP( 2), CurrentTriangle.GetP( 0))*/) - - continue ; - - // Aggiungo ultimo triangolo - triHold[nCurrent].vCompoTria[nOldFeatureNum].emplace_back( CurrentTriangle) ; - } - - nTableOffset += nVertComp ; - } - } - return true ; -} - //---------------------------------------------------------------------------- bool -VolZmap::FlipEdges( TriHolder& triHold) const +VolZmap::FlipEdges( std::vector& VecTriHold) const { - // Numero di voxel in cui si presentano sharp feature - int nVoxelNum = int( triHold.size()) ; + double dSqEps = EPS_SMALL * EPS_SMALL ; - // Ciclo su tali voxel - for ( int n = 0 ; n < nVoxelNum ; ++ n) { - for ( int m = n ; m < nVoxelNum ; ++ m) { + // Ciclo sui blocchi + for ( size_t t1 = 0 ; t1 < m_nNumBlock ; ++ t1) { + + // Determino i,j,k del primo blocco + int nIJK1[3] ; + GetBlockIJK( int( t1), nIJK1) ; + + for ( size_t t2 = t1 ; t2 < m_nNumBlock ; ++ t2) { + + // Determino i,j,k del secondo blocco + int nIJK2[3] ; + GetBlockIJK( int( t2), nIJK2) ; + + // Se i blocchi sono adiacenti o coincidenti proseguo + if ( ( nIJK2[0] >= nIJK1[0] - 1 && nIJK2[0] <= nIJK1[0] + 1) || + ( nIJK2[1] >= nIJK1[1] - 1 && nIJK2[1] <= nIJK1[1] + 1) || + ( nIJK2[2] >= nIJK1[2] - 1 && nIJK2[2] <= nIJK1[2] + 1)) { + + // Numero di voxel in cui si presentano sharp feature + int nVoxelNum1 = int( VecTriHold[t1].size()) ; + int nVoxelNum2 = int( VecTriHold[t2].size()) ; + + // Determino estremi del ciclo sui voxel esterno + int nSt1 = 0 ; + int nEn1 = nVoxelNum1 - ( t1 == t2 ? 1 : 0) ; + + // Ciclo su tali voxel + for ( int n1 = nSt1 ; n1 < nEn1 ; ++ n1) { + + // Determino estremi del ciclo sui voxel interno + int nSt2 = ( t1 == t2 ? nSt1 : 0) ; + int nEn2 = nVoxelNum2 ; + + for ( int n2 = n1 ; n2 < nVoxelNum2 ; ++ n2) { - // Voxel adiacenti o coincidenti - if ( m == n || - ( ( triHold[m].i < int( m_nNx[0]) && - triHold[m].j < int( m_nNy[0]) && - triHold[m].k < int( m_nNy[1])) && - ( triHold[m].i == triHold[n].i + 1 || - triHold[m].j == triHold[n].j + 1 || - triHold[m].k == triHold[n].k + 1))) { + // Se i voxel sono adiacenti proseguo + if ( VecTriHold[t2][n2].i == VecTriHold[t1][n1].i - 1 || + VecTriHold[t2][n2].i == VecTriHold[t1][n1].i + 1 || + VecTriHold[t2][n2].j == VecTriHold[t1][n1].j - 1 || + VecTriHold[t2][n2].j == VecTriHold[t1][n1].j + 1 || + VecTriHold[t2][n2].k == VecTriHold[t1][n1].k - 1 || + VecTriHold[t2][n2].k == VecTriHold[t1][n1].k + 1) { - // Numero delle componenti connesse nei due voxel - int nNumCompoN = triHold[n].ptCompoVert.size() ; - int nNumCompoM = triHold[m].ptCompoVert.size() ; + // Numero delle componenti connesse nei due voxel + int nNumCompo1 = int( VecTriHold[t1][n1].ptCompoVert.size()) ; + int nNumCompo2 = int( VecTriHold[t2][n2].ptCompoVert.size()) ; - int nCompoN = 0 ; + int nCompo1 = 0 ; - // Ciclo sulle componenti - for ( ; nCompoN < nNumCompoN ; ++ nCompoN) { + // Ciclo sulle componenti + for ( ; nCompo1 < nNumCompo1 ; ++ nCompo1) { - int nCompoM = ( m == n ? nCompoN + 1 : 0) ; + int nCompo2 = ( t1 == t2 && n1 == n2 ? nCompo1 + 1 : 0) ; - for ( ; nCompoM < nNumCompoM ; ++ nCompoM) { + for ( ; nCompo2 < nNumCompo2 ; ++ nCompo2) { - // Numero di triangoli per le componenti connesse - int nNumN = int( triHold[n].vCompoTria[nCompoN].size()) ; - int nNumM = int( triHold[m].vCompoTria[nCompoM].size()) ; + // Numero di triangoli per le componenti connesse + int nTriNum1 = int( VecTriHold[t1][n1].vCompoTria[nCompo1].size()) ; + int nTriNum2 = int( VecTriHold[t2][n2].vCompoTria[nCompo2].size()) ; - for ( int triN = 0 ; triN < nNumN ; ++ triN) { + for ( int nTri1 = 0 ; nTri1 < nTriNum1 ; ++ nTri1) { - bool bModified = false ; + bool bModified = false ; - for ( int triM = 0 ; triM < nNumM ; ++ triM) { + for ( int nTri2 = 0 ; nTri2 < nTriNum2 ; ++ nTri2) { - std::vector SharedIndex ; + INTVECTOR SharedIndex ; - for ( int vertN = 0 ; vertN < 3 ; ++ vertN) { + for ( int nVert1 = 0 ; nVert1 < 3 ; ++ nVert1) { - for ( int vertM = 0 ; vertM < 3 ; ++ vertM) { + for ( int nVert2 = 0 ; nVert2 < 3 ; ++ nVert2) { - Point3d ptN = triHold[n].vCompoTria[nCompoN][triN].GetP( vertN) ; - Point3d ptM = triHold[m].vCompoTria[nCompoM][triM].GetP( vertM) ; + Point3d ptP1 = VecTriHold[t1][n1].vCompoTria[nCompo1][nTri1].GetP( nVert1) ; + Point3d ptP2 = VecTriHold[t2][n2].vCompoTria[nCompo2][nTri2].GetP( nVert2) ; - if ( SqDist( ptN, ptM) < EPS_SMALL * EPS_SMALL) { + if ( SqDist( ptP1, ptP2) < dSqEps) { - Point3d ptVertN = triHold[n].ptCompoVert[nCompoN] ; - Point3d ptVertM = triHold[m].ptCompoVert[nCompoM] ; + Point3d ptVert1 = VecTriHold[t1][n1].ptCompoVert[nCompo1] ; + Point3d ptVert2 = VecTriHold[t2][n2].ptCompoVert[nCompo2] ; - if ( ! ( AreSamePointApprox( ptN, ptVertN) || - AreSamePointApprox( ptM, ptVertM))) { + if ( ! ( AreSamePointApprox( ptP1, ptVert1) || + AreSamePointApprox( ptP2, ptVert2))) { - SharedIndex.emplace_back( vertN) ; - SharedIndex.emplace_back( vertM) ; - } - } - + SharedIndex.emplace_back( nVert1) ; + SharedIndex.emplace_back( nVert2) ; + } + } - if ( SharedIndex.size() > 2) + if ( SharedIndex.size() > 2) + break ; + } + + if ( SharedIndex.size() > 2) + break ; + } + + // Si deve operare la modifica dei triangoli + if ( SharedIndex.size() > 2) { + + int nProd = ( SharedIndex[2] - SharedIndex[0]) * ( SharedIndex[3] - SharedIndex[1]) ; + + // --- + if ( nProd != 1 && nProd != - 2 && nProd != 4) { + + VecTriHold[t1][n1].vCompoTria[nCompo1][nTri1].SetP( SharedIndex[0], + VecTriHold[t2][n2].ptCompoVert[nCompo2]) ; + VecTriHold[t2][n2].vCompoTria[nCompo2][nTri2].SetP( SharedIndex[3], + VecTriHold[t1][n1].ptCompoVert[nCompo1]) ; + + VecTriHold[t1][n1].vCompoTria[nCompo1][nTri1].Validate( true) ; + VecTriHold[t2][n2].vCompoTria[nCompo2][nTri2].Validate( true) ; + + bModified = true ; + break ; + } + } + } + + if ( bModified) break ; } - - if ( SharedIndex.size() > 2) - break ; } - // Si deve operare la modifica dei triangoli - if ( SharedIndex.size() > 2) { - - // Controllo sulle normali - Point3d ptMFeature = triHold[m].ptCompoVert[nCompoM] ; - Point3d ptNFeature = triHold[n].ptCompoVert[nCompoN] ; - - Vector3d vtFeature = ptMFeature - ptNFeature ; - - vtFeature.Normalize() ; - - double dDot = vtFeature * triHold[m].vCompoTria[nCompoM][triM].GetN() ; - - // Ulteriore controllo sulle normali - Point3d ptP0 = ptNFeature ; - Point3d ptP1 = triHold[n].vCompoTria[nCompoN][triN].GetP( SharedIndex[0]) ; - Point3d ptP2 = triHold[n].vCompoTria[nCompoN][triN].GetP( SharedIndex[2]) ; - Point3d ptP3 = ptMFeature ; - - double dPlane = ( ( ptP1 - ptP0) ^ ( ptP2 - ptP0)) * ( ptP3 - ptP0) ; - - Vector3d vtN = triHold[n].vCompoTria[nCompoN][triN].GetN() ; - Vector3d vtM = triHold[m].vCompoTria[nCompoM][triM].GetN() ; - - double dDotNM = vtN * vtM ; - - int nProd = ( SharedIndex[2] - SharedIndex[0]) * ( SharedIndex[3] - SharedIndex[1]) ; - - // --- - if ( nProd != 1 && nProd != - 2 && nProd != 4) { - - triHold[n].vCompoTria[nCompoN][triN].SetP( SharedIndex[0], - triHold[m].ptCompoVert[nCompoM]) ; - triHold[m].vCompoTria[nCompoM][triM].SetP( SharedIndex[3], - triHold[n].ptCompoVert[nCompoN]) ; - - triHold[n].vCompoTria[nCompoN][triN].Validate( true) ; - triHold[m].vCompoTria[nCompoM][triM].Validate( true) ; - - bModified = true ; - break ; - } - } } - - if ( bModified) - break ; } - } + } } - } - } + } + } } - + return true ; } -//---------------------------------------------------------------------------- -bool -VolZmap::FlipEdgesLocalFlipEdges( TriaStruct& triStrCurr, TriaStruct& triStrAdj) const -{ - // Numero delle componenti connesse nei due voxel - size_t nCurVoxCompNum = triStrCurr.ptCompoVert.size() ; - size_t nAdjVoxCompNum = triStrAdj.ptCompoVert.size() ; - - // Ciclo sulle componenti connesse - for ( size_t nCVCurComp = 0 ; nCVCurComp < nCurVoxCompNum ; ++ nCVCurComp) { - for ( size_t nAVCurComp = 0 ; nAVCurComp < nAdjVoxCompNum ; ++ nAVCurComp) { - - // Numero di triangoli per le componenti connesse - size_t nCVTriaNum = triStrCurr.vCompoTria[nCVCurComp].size() ; - size_t nAVTriaNum = triStrAdj.vCompoTria[nAVCurComp].size() ; - - for ( size_t nCVT = 0 ; nCVT < nCVTriaNum ; ++ nCVT) { - - bool bModified = false ; - - for ( size_t nAVT = 0 ; nAVT < nAVTriaNum ; ++ nAVT) { - - std::vector SharedIndex ; - - for ( size_t nCurVoxVert = 0 ; nCurVoxVert < 3 ; ++ nCurVoxVert) { - for ( size_t nAdjVoxVert = 0 ; nAdjVoxVert < 3 ; ++ nAdjVoxVert) { - - Point3d ptC = triStrCurr.vCompoTria[nCVCurComp][nCVT].GetP( nCurVoxVert) ; - Point3d ptA = triStrAdj.vCompoTria[nAVCurComp][nAVT].GetP( nAdjVoxVert) ; - - if ( SqDist( ptC, ptA) < EPS_SMALL * EPS_SMALL) { - - SharedIndex.emplace_back( int( nCurVoxVert)) ; - SharedIndex.emplace_back( int( nAdjVoxVert)) ; - } - - if ( SharedIndex.size() > 2) - break ; - } - - if ( SharedIndex.size() > 2) - break ; - } - // Si deve operare la modifica dei triangoli - if ( SharedIndex.size() > 2) { - - // Controllo sulle normali - Point3d ptCurFeature = triStrCurr.ptCompoVert[nCVCurComp] ; - Point3d ptAdjFeature = triStrAdj.ptCompoVert[nAVCurComp] ; - - Vector3d vtFeature = ptAdjFeature - ptCurFeature ; - - vtFeature.Normalize() ; - - double dDot = vtFeature * triStrAdj.vCompoTria[nAVCurComp][nAVT].GetN() ; - - // Ulteriore controllo sulle normali - Point3d ptP0 = ptCurFeature ; - Point3d ptP1 = triStrCurr.vCompoTria[nCVCurComp][nCVT].GetP( SharedIndex[0]) ; - Point3d ptP2 = triStrCurr.vCompoTria[nCVCurComp][nCVT].GetP( SharedIndex[2]) ; - Point3d ptP3 = ptAdjFeature ; - - double dPlane = ( ( ptP1 - ptP0) ^ ( ptP2 - ptP0)) * ( ptP3 - ptP0) ; - - Vector3d vtN = triStrCurr.vCompoTria[nCVCurComp][nCVT].GetN() ; - Vector3d vtM = triStrAdj.vCompoTria[nAVCurComp][nAVT].GetN() ; - - double dDotNM = vtN * vtM ; - - int nProd = ( SharedIndex[2] - SharedIndex[0]) * ( SharedIndex[3] - SharedIndex[1]) ; - - if ( nProd != 1 && nProd != - 2 && nProd != 4 && - ( dDot < - EPS_SMALL || ( dPlane < EPS_SMALL && dDotNM < - 0.95))) { - - triStrCurr.vCompoTria[nCVCurComp][nCVT].SetP( SharedIndex[0], - triStrAdj.ptCompoVert[nAVCurComp]) ; - triStrAdj.vCompoTria[nAVCurComp][nAVT].SetP( SharedIndex[3], - triStrCurr.ptCompoVert[nCVCurComp]) ; - - triStrCurr.vCompoTria[nCVCurComp][nCVT].Validate( true) ; - triStrAdj.vCompoTria[nAVCurComp][nAVT].Validate( true) ; - - bModified = true ; - break ; - } - } - } - - if ( bModified) - break ; - } - } - } - return true ; -} - //---------------------------------------------------------------------------- bool VolZmap::IsThereMat( int nI, int nJ, int nK) const @@ -2071,221 +1706,13 @@ VolZmap::IntersPos( int nVec1[], int nVec2[], Point3d& ptInt) const return true ; } -//---------------------------------------------------------------------------- -bool -VolZmap::IntersPos( int nVec1[], int nVec2[], Point3d& ptInt, Vector3d& vtNormal) const -{ - double Eps = EPS_SMALL ; - - if ( nVec1[0] != nVec2[0]) { - - int nMinI = min( nVec1[0], nVec2[0]) ; - int nMaxI = max( nVec1[0], nVec2[0]) ; - - double dMinX = ( nMinI + 0.5) * m_dStep ; - double dMaxX = ( nMaxI + 0.5) * m_dStep ; - - ptInt.y = ( nVec1[1] + 0.5) * m_dStep ; - ptInt.z = ( nVec1[2] + 0.5) * m_dStep ; - - unsigned int nDexel = nVec1[2] * m_nNx[1] + nVec1[1] ; - size_t nSize = m_Values[1][nDexel].size() ; - - bool bFound = false ; - for ( size_t i = 0 ; i < nSize ; i += 2) { - - double dx1 = m_Values[1][nDexel][i].dZVal ; - double dx2 = m_Values[1][nDexel][i+1].dZVal ; - - if ( dx1 < dMinX - Eps && dx2 > dMinX - Eps && dx2 < dMaxX + Eps) { - ptInt.x = dx2 ; - vtNormal = m_Values[1][nDexel][i+1].vtN ; - bFound = true ; - break ; - } - else if ( dx1 > dMinX - Eps && dx1 < dMaxX + Eps && dx2 > dMaxX + Eps) { - ptInt.x = dx1 ; - vtNormal = m_Values[1][nDexel][i].vtN ; - bFound = true ; - break ; - }/* - if ( dx1 < dMinX + Eps && dx2 > dMinX && dx2 < dMaxX - Eps) { - ptInt.x = dx2 ; - vtNormal = m_Values[1][nDexel][i+1].vtN ; - bFound = true ; - break ; - } - else if ( dx1 > dMinX + Eps && dx1 < dMaxX + Eps && dx2 > dMaxX + Eps) { - ptInt.x = dx1 ; - vtNormal = m_Values[1][nDexel][i].vtN ; - bFound = true ; - break ; - }*/ - } - if ( ! bFound) { - - ptInt.x = ( dMinX + dMaxX) / 2 ; - /* - if ( IsThereMat( nVec1[0], nVec1[1], nVec1[2])) { - - double dY = ( nVec2[1] + 0.5) * m_dStep ; - double dZ = ( nVec2[2] + 0.5) * m_dStep ; - - unsigned int YDirMinIndex = 0 ; - unsigned int ZDirMinIndex = 0 ; - - double dYMinDist = DBL_MAX ; - double dZMinDist = DBL_MAX ; - - unsigned int nYDirGrid = nVec2[0] * m_nNx[2] + nVec2[2] ; - unsigned int nZDirGrid = nVec2[1] * m_nNx[0] + nVec2[0] ; - - for ( unsigned int t = 0 ; t < m_Values[0][nZDirGrid].size() ; ++ t) { - - double dMZ1 = m_Values[0][nZDirGrid][t].dZVal ; - - if ( abs( dMZ1 - dZ) < dZMinDist) { - - ZDirMinIndex = t ; - dZMinDist = abs( dMZ1 - dZ) ; - } - } - - for ( unsigned int t = 0 ; t < m_Values[2][nYDirGrid].size() ; ++ t) { - - double dMY1 = m_Values[2][nYDirGrid][t].dZVal ; - - if ( abs( dMY1 - dY) < dYMinDist) { - - YDirMinIndex = t ; - dYMinDist = abs( dMY1 - dY) ; - } - } - - unsigned int AbsoluteMinIndex ; - - if ( dZMinDist < dYMinDist) - - vtNormal = m_Values[0][nZDirGrid][ZDirMinIndex].vtN ; - else - vtNormal = m_Values[0][nYDirGrid][YDirMinIndex].vtN ; - } - else { - - }*/ - } - } - - else if ( nVec1[1] != nVec2[1]) { - - int nMinJ = min( nVec1[1], nVec2[1]) ; - int nMaxJ = max( nVec1[1], nVec2[1]) ; - - double dMinY = ( nMinJ + 0.5) * m_dStep ; - double dMaxY = ( nMaxJ + 0.5) * m_dStep ; - - ptInt.x = ( nVec1[0] + 0.5) * m_dStep ; - ptInt.z = ( nVec1[2] + 0.5) * m_dStep ; - - unsigned int nDexel = nVec1[0] * m_nNx[2] + nVec1[2] ; - size_t nSize = m_Values[2][nDexel].size() ; - - bool bFound = false ; - for ( size_t j = 0 ; j < nSize ; j += 2) { - - double dy1 = m_Values[2][nDexel][j].dZVal ; - double dy2 = m_Values[2][nDexel][j+1].dZVal ; - - if ( dy1 < dMinY - Eps && dy2 > dMinY - Eps && dy2 < dMaxY + Eps) { - ptInt.y = dy2 ; - vtNormal = m_Values[2][nDexel][j+1].vtN ; - bFound = true ; - break ; - } - else if ( dy1 > dMinY - Eps && dy1 < dMaxY + Eps && dy2 > dMaxY + Eps) { - ptInt.y = dy1 ; - vtNormal = m_Values[2][nDexel][j].vtN ; - bFound = true ; - break ; - }/* - if ( dy1 < dMinY + Eps && dy2 > dMinY && dy2 < dMaxY - Eps) { - ptInt.y = dy2 ; - vtNormal = m_Values[2][nDexel][j+1].vtN ; - bFound = true ; - break ; - } - else if ( dy1 > dMinY + Eps && dy1 < dMaxY + Eps && dy2 > dMaxY + Eps) { - ptInt.y = dy1 ; - vtNormal = m_Values[2][nDexel][j].vtN ; - bFound = true ; - break ; - }*/ - } - if ( ! bFound) { - ptInt.y = ( dMinY + dMaxY) / 2 ; - // Versore Normale ??? - } - } - - else if ( nVec1[2] != nVec2[2]) { - - int nMinK = min( nVec1[2], nVec2[2]) ; - int nMaxK = max( nVec1[2], nVec2[2]) ; - - double dMinZ = ( nMinK + 0.5) * m_dStep ; - double dMaxZ = ( nMaxK + 0.5) * m_dStep ; - - ptInt.x = ( nVec1[0] + 0.5) * m_dStep ; - ptInt.y = ( nVec1[1] + 0.5) * m_dStep ; - - unsigned int nDexel = nVec1[1] * m_nNx[0] + nVec1[0] ; - size_t nSize = m_Values[0][nDexel].size() ; - - bool bFound = false ; - for ( size_t k = 0 ; k < nSize ; k += 2) { - - double dz1 = m_Values[0][nDexel][k].dZVal ; - double dz2 = m_Values[0][nDexel][k+1].dZVal ; - - if ( dz1 < dMinZ - Eps && dz2 > dMinZ - Eps && dz2 < dMaxZ + Eps) { - ptInt.z = dz2 ; - vtNormal = m_Values[0][nDexel][k+1].vtN ; - bFound = true ; - break ; - } - else if ( dz1 > dMinZ - Eps && dz1 < dMaxZ + Eps && dz2 > dMaxZ + Eps) { - ptInt.z = dz1 ; - vtNormal = m_Values[0][nDexel][k].vtN ; - bFound = true ; - break ; - }/* - if ( dz1 < dMinZ + Eps && dz2 > dMinZ && dz2 < dMaxZ - Eps) { - ptInt.z = dz2 ; - vtNormal = m_Values[0][nDexel][k+1].vtN ; - bFound = true ; - break ; - } - else if ( dz1 > dMinZ + Eps && dz1 < dMaxZ + Eps && dz2 > dMaxZ + Eps) { - ptInt.z = dz1 ; - vtNormal = m_Values[0][nDexel][k].vtN ; - bFound = true ; - break ; - }*/ - } - if ( ! bFound) { - ptInt.z = ( dMinZ + dMaxZ) / 2 ; - // Versore Normale ??? - } - } - - return true ; -} - //---------------------------------------------------------------------------- bool -VolZmap::NewIntersPos( int nVec1[], int nVec2[], bool bFirstCorner, Point3d& ptInt, Vector3d& vtNormal) const +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]) { @@ -2312,6 +1739,7 @@ VolZmap::NewIntersPos( int nVec1[], int nVec2[], bool bFirstCorner, Point3d& ptI ptInt.x = dX ; vtNormal = m_Values[1][nDexel][n].vtN ; + bFound = true ; break ; } @@ -2333,6 +1761,7 @@ VolZmap::NewIntersPos( int nVec1[], int nVec2[], bool bFirstCorner, Point3d& ptI ptInt.x = dX ; vtNormal = m_Values[1][nDexel][n].vtN ; + bFound = true ; break ; } @@ -2342,7 +1771,11 @@ VolZmap::NewIntersPos( int nVec1[], int nVec2[], bool bFirstCorner, Point3d& ptI n += 2 ; dX = m_Values[1][nDexel][n].dZVal ; } - } + } + + if ( ! bFound) + + ptInt.x = 0.5 * ( dMinX + dMaxX) ; } else if ( nVec1[1] != nVec2[1]) { @@ -2370,6 +1803,7 @@ VolZmap::NewIntersPos( int nVec1[], int nVec2[], bool bFirstCorner, Point3d& ptI ptInt.y = dY ; vtNormal = m_Values[2][nDexel][n].vtN ; + bFound = true ; break ; } @@ -2391,6 +1825,7 @@ VolZmap::NewIntersPos( int nVec1[], int nVec2[], bool bFirstCorner, Point3d& ptI ptInt.y = dY ; vtNormal = m_Values[2][nDexel][n].vtN ; + bFound = true ; break ; } @@ -2400,7 +1835,11 @@ VolZmap::NewIntersPos( int nVec1[], int nVec2[], bool bFirstCorner, Point3d& ptI n += 2 ; dY = m_Values[2][nDexel][n].dZVal ; } - } + } + + if ( ! bFound) + + ptInt.y = 0.5 * ( dMinY + dMaxY) ; } else if ( nVec1[2] != nVec2[2]) { @@ -2428,6 +1867,7 @@ VolZmap::NewIntersPos( int nVec1[], int nVec2[], bool bFirstCorner, Point3d& ptI ptInt.z = dZ ; vtNormal = m_Values[0][nDexel][n].vtN ; + bFound = true ; break ; } @@ -2449,6 +1889,7 @@ VolZmap::NewIntersPos( int nVec1[], int nVec2[], bool bFirstCorner, Point3d& ptI ptInt.z = dZ ; vtNormal = m_Values[0][nDexel][n].vtN ; + bFound = true ; break ; } @@ -2458,18 +1899,22 @@ VolZmap::NewIntersPos( int nVec1[], int nVec2[], bool bFirstCorner, Point3d& ptI n += 2 ; dZ = m_Values[0][nDexel][n].dZVal ; } - } + } + + if ( ! bFound) + + ptInt.z = 0.5 * ( dMinZ + dMaxZ) ; } - return true ; + return bFound ; } //---------------------------------------------------------------------------- bool -VolZmap::GetBlockIJK( int nIJK[], int nBlock) const +VolZmap::GetBlockIJK( int nBlock, int nIJK[]) const { // Controllo sulla validità del blocco - if ( nBlock < 0 || nBlock >= int( m_nNumBlock)) + if ( nBlock < 0 || nBlock >= int( m_nNumBlock)) return false ; // Calcolo posizione del blocco nel reticolo @@ -2488,7 +1933,6 @@ VolZmap::GetBlockLimitsIJK( const int nIJK[], int nLimits[]) const if ( nIJK[0] < 0 || nIJK[0] >= int( m_nFracLin[0]) || nIJK[1] < 0 || nIJK[1] >= int( m_nFracLin[1]) || nIJK[2] < 0 || nIJK[2] >= int( m_nFracLin[2])) - return false ; // Calcolo limiti per l'indice i @@ -2508,3 +1952,62 @@ VolZmap::GetBlockLimitsIJK( const int nIJK[], int nLimits[]) const return true ; } + +//---------------------------------------------------------------------------- +bool +VolZmap::IsPointInsideVoxel( int nI, int nJ, int nK, const Point3d& ptP) const +{ + // Controllo sull'ammissibilità del voxel + if ( nI < - 1 || nI >= int( m_nNx[0]) || + nJ < - 1 || nJ >= int( m_nNy[0]) || + nK < - 1 || nK >= int( m_nNy[1])) + return false ; + + int nPointI = int( floor( ( ptP.x - 0.5 * m_dStep) / m_dStep)) ; + int nPointJ = int( floor( ( ptP.y - 0.5 * m_dStep) / m_dStep)) ; + int nPointK = int( floor( ( ptP.z - 0.5 * m_dStep) / m_dStep)) ; + + return ( nPointI == nI && + nPointJ == nJ && + nPointK == nK) ; +} + +//---------------------------------------------------------------------------- +bool +VolZmap::IsPointInsideVoxelApprox( int nI, int nJ, int nK, const Point3d& ptP) const +{ + // Controllo sull'ammissibilità del voxel + if ( nI < - 1 || nI >= int( m_nNx[0]) || + nJ < - 1 || nJ >= int( m_nNy[0]) || + nK < - 1 || nK >= int( m_nNy[1])) + return false ; + + Point3d ptPZmap = ptP ; + ptPZmap.ToLoc( m_MapFrame[0]) ; + + bool bI = ptPZmap.x > ( nI + 0.5) * m_dStep - EPS_SMALL && + ptPZmap.x < ( nI + 1.5) * m_dStep + EPS_SMALL ; + bool bJ = ptPZmap.y > ( nJ + 0.5) * m_dStep - EPS_SMALL && + ptPZmap.y < ( nJ + 1.5) * m_dStep + EPS_SMALL ; + bool bK = ptPZmap.z > ( nK + 0.5) * m_dStep - EPS_SMALL && + ptPZmap.z < ( nK + 1.5) * m_dStep + EPS_SMALL ; + + return ( bI && bJ && bK) ; +} + +//---------------------------------------------------------------------------- +bool +VolZmap::GetPointVoxel( Point3d& ptP, int& nVoxI, int& nVoxJ, int& nVoxK) const +{ + nVoxI = int( floor( ( ptP.x - 0.5 * m_dStep) / m_dStep)) ; + nVoxJ = int( floor( ( ptP.y - 0.5 * m_dStep) / m_dStep)) ; + nVoxK = int( floor( ( ptP.z - 0.5 * m_dStep) / m_dStep)) ; + + return ( nVoxI >= -1 && nVoxI < int( m_nNx[0])) && + ( nVoxJ >= -1 && nVoxJ < int( m_nNy[0])) && + ( nVoxK >= -1 && nVoxK < int( m_nNy[1])) ; +} + + + + diff --git a/VolTriZmapVolume.cpp b/VolTriZmapVolume.cpp index fdd4eb1..5a8358f 100644 --- a/VolTriZmapVolume.cpp +++ b/VolTriZmapVolume.cpp @@ -169,50 +169,6 @@ VolZmap::SubtractIntervals( unsigned int nGrid, unsigned int nI, unsigned int nJ } } } - /* - // Blocco adiacente in X - if ( nXBlock > 0 && - ( nI % m_nDexNumPBlock == 0) && - ( nI / m_nDexNumPBlock <= m_nFracLin[0] - 1)) { - - int nAdjXBlock = nXBlock - 1 ; - - for ( int k = nMinZBlock ; k <= nMaxZBlock ; ++ k) { - - int nBlockNum = k * nLayerBlock + nYBlock * m_nFracLin[0] + nAdjXBlock ; - m_BlockToUpdate[nBlockNum] = true ; - } - } - // Blocco adiacente in Y - if ( nYBlock > 0 && - ( nJ % m_nDexNumPBlock == 0) && - ( nJ / m_nDexNumPBlock <= m_nFracLin[1] - 1)) { - - int nAdjYBlock = nYBlock - 1 ; - - for ( int k = nMinZBlock ; k <= nMaxZBlock ; ++ k) { - - int nBlockNum = k * nLayerBlock + nAdjYBlock * m_nFracLin[0] + nXBlock ; - m_BlockToUpdate[nBlockNum] = true ; - } - } - // Blocco adiacente in XY - if ( ( nXBlock > 0 && - ( nI % m_nDexNumPBlock == 0) && - ( nI / m_nDexNumPBlock <= m_nFracLin[0] - 1)) && - ( nYBlock > 0 && - ( nJ % m_nDexNumPBlock == 0) && - ( nJ / m_nDexNumPBlock <= m_nFracLin[1] - 1))) { - - int nAdjXBlock = nXBlock - 1 ; - int nAdjYBlock = nYBlock - 1 ; - - for ( int k = nMinZBlock ; k <= nMaxZBlock ; ++ k) { - - int nBlockNum = k * nLayerBlock + nAdjYBlock * m_nFracLin[0] + nAdjXBlock ; - m_BlockToUpdate[nBlockNum] = true ; - } - }*/ } else if ( nGrid == 1) { @@ -238,50 +194,6 @@ VolZmap::SubtractIntervals( unsigned int nGrid, unsigned int nI, unsigned int nJ } } } - /* - // Blocco adiacente in Y - if ( nYBlock > 0 && - ( nI % m_nDexNumPBlock == 0) && - ( nI / m_nDexNumPBlock <= m_nFracLin[1] - 1)) { - - int nAdjYBlock = nYBlock - 1 ; - - for ( int k = nMinXBlock ; k <= nMaxXBlock ; ++ k) { - - int nBlockNum = nZBlock * nLayerBlock + nAdjYBlock * m_nFracLin[0] + k ; - m_BlockToUpdate[nBlockNum] = true ; - } - } - // Blocco adiacente in Z - if ( nZBlock > 0 && - ( nJ % m_nDexNumPBlock == 0) && - ( nJ / m_nDexNumPBlock <= m_nFracLin[2] - 1)) { - - int nAdjZBlock = nZBlock - 1 ; - - for ( int k = nMinXBlock ; k <= nMaxXBlock ; ++ k) { - - int nBlockNum = nAdjZBlock * nLayerBlock + nYBlock * m_nFracLin[0] + k ; - m_BlockToUpdate[nBlockNum] = true ; - } - } - // Blocco adiacente in YZ - if ( ( nYBlock > 0 && - ( nI % m_nDexNumPBlock == 0) && - ( nI / m_nDexNumPBlock <= m_nFracLin[1] - 1)) && - ( nZBlock > 0 && - ( nJ % m_nDexNumPBlock == 0) && - ( nJ / m_nDexNumPBlock <= m_nFracLin[2] - 1))) { - - int nAdjYBlock = nYBlock - 1 ; - int nAdjZBlock = nZBlock - 1 ; - - for ( int k = nMinXBlock ; k <= nMaxXBlock ; ++ k) { - - int nBlockNum = nAdjZBlock * nLayerBlock + nAdjYBlock * m_nFracLin[0] + k ; - m_BlockToUpdate[nBlockNum] = true ; - } - } */ } else if ( nGrid == 2) { @@ -307,50 +219,6 @@ VolZmap::SubtractIntervals( unsigned int nGrid, unsigned int nI, unsigned int nJ } } } - /* - // Blocchi adiacenti in X - if ( nXBlock > 0 && - ( nJ % m_nDexNumPBlock == 0) && - ( nJ / m_nDexNumPBlock <= m_nFracLin[0] - 1)) { - - int nAdjXBlock = nXBlock - 1 ; - - for ( int k = nMinYBlock ; k <= nMaxYBlock ; ++ k) { - - int nBlockNum = nZBlock * nLayerBlock + k * m_nFracLin[0] + nAdjXBlock ; - m_BlockToUpdate[nBlockNum] = true ; - } - } - // Blocchi adiacenti in Z - if ( nZBlock > 0 && - ( nI % m_nDexNumPBlock == 0) && - ( nI / m_nDexNumPBlock <= m_nFracLin[2] - 1)) { - - int nAdjZBlock = nZBlock - 1 ; - - for ( int k = nMinYBlock ; k <= nMaxYBlock ; ++ k) { - - int nBlockNum = nAdjZBlock * nLayerBlock + k * m_nFracLin[0] + nXBlock ; - m_BlockToUpdate[nBlockNum] = true ; - } - } - // Blocchi adiacenti in XZ - if ( ( nXBlock > 0 && - ( nJ % m_nDexNumPBlock == 0) && - ( nJ / m_nDexNumPBlock <= m_nFracLin[0] - 1)) && - ( nZBlock > 0 && - ( nI % m_nDexNumPBlock == 0) && - ( nI / m_nDexNumPBlock <= m_nFracLin[2] - 1))) { - - int nAdjXBlock = nXBlock - 1 ; - int nAdjZBlock = nZBlock - 1 ; - - for ( int k = nMinYBlock ; k <= nMaxYBlock ; ++ k) { - - int nBlockNum = nAdjZBlock * nLayerBlock + k * m_nFracLin[0] + nAdjXBlock ; - m_BlockToUpdate[nBlockNum] = true ; - } - }*/ } m_OGrMgr.Reset() ; @@ -1162,7 +1030,7 @@ VolZmap::CylBall_ZMilling( unsigned int nGrid, const Point3d & ptS, const Point3 Vector3d vtTan( - vtCirc.y, vtCirc.x, 0) ; // Tangente alla circonferenza Vector3d vtCross = vtTan ^ vtMove ; - vtNmax = ( vtCross * vtCirc > 0 ? vtCross : - vtCross) ; + vtNmax = ( vtCross * vtCirc > - EPS_ZERO ? vtCross : - vtCross) ; vtNmax.Normalize() ; } // Minimo @@ -1180,7 +1048,7 @@ VolZmap::CylBall_ZMilling( unsigned int nGrid, const Point3d & ptS, const Point3 Vector3d vtTan( - vtCirc.y, vtCirc.x, 0) ; // Tangente alla circonferenza Vector3d vtCross = vtTan ^ vtMove ; - vtNmin = ( vtCross * vtCirc > 0 ? vtCross : - vtCross) ; + vtNmin = ( vtCross * vtCirc > - EPS_ZERO ? vtCross : - vtCross) ; vtNmin.Normalize() ; } @@ -2636,6 +2504,12 @@ VolZmap::Conus_XYPerp( unsigned int nGrid, const Point3d& ptS, const Point3d& pt Vector3d vtUpCross = vtMove ^ vtUpTan ; Vector3d vtDwCross = - vtMove ^ vtDwTan ; + if ( vtUpCross.z > 0) + vtUpCross = - vtUpCross ; + + if ( vtDwCross.z < 0) + vtDwCross = - vtDwCross ; + // Descrizione piani tangenti al cono Vector3d vtR0Up = ptUp - ORIG ; Vector3d vtR0Dw = ptDw - ORIG ; @@ -3882,7 +3756,7 @@ VolZmap::CompCyl_ZMilling( unsigned int nGrid, const Point3d & ptS, const Point3 Vector3d vtTan( - vtCirc.y, vtCirc.x, 0) ; Vector3d vtCross = vtTan ^ vtMove ; - vtMax = ( vtCross * vtCirc > 0 ? vtCross : - vtCross) ; + vtMax = ( vtCross * vtCirc > - EPS_ZERO ? vtCross : - vtCross) ; vtMax.Normalize() ; } // Minimo @@ -3900,7 +3774,7 @@ VolZmap::CompCyl_ZMilling( unsigned int nGrid, const Point3d & ptS, const Point3 Vector3d vtTan( - vtCirc.y, vtCirc.x, 0) ; Vector3d vtCross = vtTan ^ vtMove ; - vtMin = ( vtCross * vtCirc > 0 ? vtCross : - vtCross) ; + vtMin = ( vtCross * vtCirc > - EPS_ZERO ? vtCross : - vtCross) ; vtMin.Normalize() ; } SubtractIntervals( nGrid, i, j, dMin, dMax, vtMin, vtMax) ; @@ -4171,9 +4045,9 @@ VolZmap::CompConus_ZMilling( unsigned int nGrid, const Point3d & ptS, const Poin Vector3d vtCirc = - dIDL_0 * vtV2 - dIDO * vtV3 ; Vector3d vtTan( - vtCirc.y, vtCirc.x, 0) ; Vector3d vtCross = vtTan ^ vtUmv ; - + vtP = ( vtCross * vtCirc > 0 ? vtCross : - vtCross) ; - vtP.Normalize() ; + vtP.Normalize() ; } // Limiti nella direzione negativa di vtV1 @@ -4229,9 +4103,10 @@ VolZmap::CompConus_ZMilling( unsigned int nGrid, const Point3d & ptS, const Poin Vector3d vtCirc = - dIDL_0 * vtV2 - dIDO * vtV3 ; Vector3d vtTan( - vtCirc.y, vtCirc.x, 0) ; Vector3d vtCross = vtTan ^ vtUmv ; - - vtP = ( vtCross * vtCirc > 0 ? vtCross : - vtCross) ; - vtP.Normalize() ; + //////////////////////////////////////////////////////////////////////////////////////// + vtP = ( vtCross * vtCirc > - EPS_ZERO ? vtCross : - vtCross) ; // vtCross * vtCirc o vtCross * vtMove? + //////////////////////////////////////////////////////////////////////////////////////// + vtP.Normalize() ; } // Limiti nella direzione negativa di vtV1 @@ -4272,9 +4147,9 @@ VolZmap::CompConus_ZMilling( unsigned int nGrid, const Point3d & ptS, const Poin Vector3d vtCirc = - dIDL_0 * vtV2 - dIDO * vtV3 ; Vector3d vtTan( - vtCirc.y, vtCirc.x, 0) ; Vector3d vtCross = vtTan ^ vtUmv ; - - vtM = ( vtCross * vtMove > 0 ? vtCross : - vtCross) ; - vtM.Normalize() ; + + vtM = ( vtCross * vtMove > - EPS_ZERO ? vtCross : - vtCross) ; + vtM.Normalize() ; } else if ( dIDO >= dMinRad * dSin && dIDO < dMaxRad * dSin) { @@ -4294,9 +4169,10 @@ VolZmap::CompConus_ZMilling( unsigned int nGrid, const Point3d & ptS, const Poin Vector3d vtCirc = - dIDL_0 * vtV2 - dIDO * vtV3 ; Vector3d vtTan( - vtCirc.y, vtCirc.x, 0) ; Vector3d vtCross = vtTan ^ vtUmv ; - vtM = ( vtCross * vtMove > 0 ? vtCross : - vtCross) ; - vtM.Normalize() ; + vtM = ( vtCross * vtMove > - EPS_ZERO ? vtCross : - vtCross) ; + vtM.Normalize() ; + } else { @@ -4323,8 +4199,9 @@ VolZmap::CompConus_ZMilling( unsigned int nGrid, const Point3d & ptS, const Poin Vector3d vtTan( - vtCirc.y, vtCirc.x, 0) ; Vector3d vtCross = vtTan ^ vtMove ; - vtM = ( vtCross * vtMove > 0 ? vtCross : - vtCross) ; - vtM.Normalize() ; + + vtM = ( vtCross * vtMove > - EPS_ZERO ? vtCross : - vtCross) ; + vtM.Normalize() ; } } diff --git a/VolZmap.h b/VolZmap.h index 6a48f53..45d103c 100644 --- a/VolZmap.h +++ b/VolZmap.h @@ -74,8 +74,7 @@ class VolZmap : public IVolZmap, public IGeoObjRW bool CreateFromTriMesh( const ISurfTriMesh& Surf, double dPrec, bool bTriDex) override ; bool GetAllTriangles( TRIA3DLIST& lstTria) const override ; int GetBlockCount( void) const override ; - bool GetBlockInfo( std::vector & bModified) const override ; - bool GetBlockTriangles( int nBlock, TRIA3DLIST& lstTria) const override ; + bool GetTriangles( bool bAllBlocks, INTVECTOR& nModifiedBlocks, TRIA3DLISTVECTOR& vLstTria) const override ; bool GetDexelLines( int nDir, int nPos1, int nPos2, POLYLINELIST& lstPL) const override ; bool SetTolerances( double dLinTol, double dAngTolDeg = 90) override ; bool SetStdTool( const std::string& pToolName, double dH, double dR, double dCornR) override ; @@ -117,15 +116,15 @@ class VolZmap : public IVolZmap, public IGeoObjRW const Vector3d& vtZ, const Vector3d& vtNorm, TRIA3DLIST& lstTria) const ; bool MarchingCubes( TRIA3DLIST& lstTria) const ; bool MarchingCubes( int nBlock, TRIA3DLIST& lstTria) const ; - bool ExtMarchingCubes( const int nLimits[], TRIA3DLIST& lstTria, TriHolder& triHold) const ; - bool ExtMarchingCubes( const std::vector VoxelsIndexes, TriHolder& triHold) const ; - bool FlipEdges( TriHolder& triHold) const ; - bool FlipEdgesLocalFlipEdges( TriaStruct& triStrCurr, TriaStruct& triStrAdj) const ; + bool ExtMarchingCubes( const int nLimits[], TRIA3DLIST& lstTria, TriHolder& triHold) const ; + bool FlipEdges( std::vector& VecTriHold) 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 ; bool IntersPos( int nVec1[], int nVec2[], Point3d & ptInt) const ; - bool IntersPos( int nVec1[], int nVec2[], Point3d & ptInt, Vector3d & vtNormal) const ; - bool NewIntersPos( int nVec1[], int nVec2[], bool bFirstCorner, Point3d& ptInt, Vector3d& vtNormal) 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 ; + bool IsPointInsideVoxelApprox( int nI, int nJ, int nK, const Point3d& ptP) const ; + bool GetPointVoxel( Point3d& ptP, int& nVoxI, int& nVoxJ, int& nVoxK) const ; // OPERAZIONI SU INTERVALLI bool SubtractIntervals( unsigned int nGrid, unsigned int nI, unsigned int nJ, @@ -215,6 +214,8 @@ class VolZmap : public IVolZmap, public IGeoObjRW // Intersezioni bool IntersLineBox( const Point3d& ptP, const Vector3d& vtV, const Point3d& ptMin, const Point3d& ptMax, double& dU1, double& dU2) const ; + bool IntersLineVoxel( const Point3d& ptP, const Vector3d& vtV, int nIndI, int nIndJ, int nIndK, + int& nFace1, int& nFace2, double& dU1, double& dU2) const ; bool IntersLineZMapBBox( unsigned int nGrid, const Point3d& ptP, const Vector3d& vtV, double& dU1, double& dU2) ; bool IntersLineDexel( unsigned int nGrid, const Point3d& ptP, const Vector3d& vtV, unsigned int nI, unsigned int nJ, double& dU1, double& dU2) ; @@ -244,13 +245,13 @@ class VolZmap : public IVolZmap, public IGeoObjRW Point3d& ptInt1, Point3d& ptInt2, Vector3d& vtN1, Vector3d& vtN2) ; // Funzioni di gestione dei blocchi - bool GetBlockIJK( int nIJK[], int nBlock) const ; + bool GetBlockIJK( int nBlock, int nIJK[]) const ; bool GetBlockLimitsIJK( const int nIJK[], int nLimits[]) const ; private : enum Status { ERR = 0, OK = 1, TO_VERIFY = 2} ; static const int N_MAPS = 3 ; - static const int N_DEXBLOCK = 32;//20;//32 ; + static const int N_DEXBLOCK = 32 ; // 10000 ;//20 ;//32 ; private : ObjGraphicsMgr m_OGrMgr ; // gestore grafica dell'oggetto