From e00d048d5d0875c66aaee7c78ba084fa585e064a Mon Sep 17 00:00:00 2001 From: Dario Sassi Date: Fri, 23 Mar 2018 11:16:56 +0000 Subject: [PATCH] EgtGeomkernel : - modifiche a Zmap per visualizzazione con catene di feature. --- VolZmap.cpp | 2 + VolZmap.h | 87 +-- VolZmapCalculus.cpp | 4 +- VolZmapCreation.cpp | 11 +- VolZmapGraphics.cpp | 1496 +++++++++++++++++++++++-------------------- 5 files changed, 855 insertions(+), 745 deletions(-) diff --git a/VolZmap.cpp b/VolZmap.cpp index 658a68a..d88ae3b 100644 --- a/VolZmap.cpp +++ b/VolZmap.cpp @@ -1082,6 +1082,8 @@ VolZmap::ClonePart( int nPart) const for ( unsigned int nCount = 0 ; nCount < pVolume->m_nNumBlock ; ++ nCount) pVolume->m_BlockToUpdate[nCount] = true ; + // Dimensiono raccolta di voxel di confine + pVolume->m_InterBlockVox.resize( pVolume->m_nNumBlock) ; // Dimensiono raccolta triangoli di feature tra blocchi pVolume->m_InterBlockTria.resize( pVolume->m_nNumBlock) ; diff --git a/VolZmap.h b/VolZmap.h index 9c7af85..5a87906 100644 --- a/VolZmap.h +++ b/VolZmap.h @@ -27,26 +27,11 @@ #include #include -//---------------------------------------------------------------------------- -struct TriaStruct { - int i, j, k ; - PNTVECTOR ptCompoVert ; - std::vector vCompoTria ; - std::vector> vbFlipped ; -} ; - -// Vettore di TriaStruct con sharp fature interni a un blocco -typedef std::vector TriHolder ; - -// Vettore di TriHolder con sharp feature di frontiera: -// il primo indice individua il blocco, il secondo il voxel -typedef std::vector TriaMatrix ; - -// ------------------------- STRUTTURA VERTICE TRIANGOLO - NORMALE ALLA SUPERFICIE ------------------------------------------------ -struct VectorField { - Point3d ptInt ; - Vector3d vtNorm ; - int nToolFlag ; +// ------------------------- STRUTTURE ----------------------------------------------------------- +struct AppliedVector { + Point3d ptPApp ; + Vector3d vtVec ; + int nPropIndex ; } ; //---------------------------------------------------------------------------- @@ -143,11 +128,38 @@ class VolZmap : public IVolZmap, public IGeoObjRW HeigthAndColor( void) : nTool( 0), dHeigth( 0) {} HeigthAndColor( int nT, double dH) : nTool( nT), dHeigth( dH) {} } ; - typedef std::unordered_map VoxelContainer ; - // Struttura voxel - struct Voxel { + typedef std::unordered_map FlatVoxelContainer ; + // Struttura indici di voxel + struct VoxelIndexes { int nI, nJ, nK ; } ; + // Struttura per componenti connesse + struct ConComp { + Point3d ptVert ; + AppliedVector CompVecField[7] ; + int nVertNum ; + bool bInside ; + bool bCorner ; + } ; + // Struttura per Voxel + struct Voxel { + int i, j, k ; + int nNumComp ; + ConComp Compo[4] ; + } ; + // Triangoli + struct TriaStruct { + int i, j, k ; + PNTVECTOR ptCompoVert ; + std::vector vCompoTria ; + std::vector vbFlipped ; + } ; + // Vettore di TriaStruct con sharp feature interni a un blocco + typedef std::vector TriHolder ; + // Vettore di TriHolder con sharp feature di frontiera: il primo indice individua il blocco, il secondo il voxel + typedef std::vector TriaMatrix ; + // Tavola hash di Voxel + typedef std::unordered_map VoxelContainer ; private : bool CopyFrom( const VolZmap& clSrc) ; @@ -158,27 +170,30 @@ class VolZmap : public IVolZmap, public IGeoObjRW bool AddDexelSideFace( int nPos, int nPosAdj, const Point3d& ptP, const Point3d& ptQ, const Vector3d& vtZ, const Vector3d& vtNorm, TRIA3DEXLIST& lstTria) const ; bool ProcessCube( int nVoxI, int nVoxJ, int nVoxK, int& nCubeIndex, TRIA3DEXLIST& lstTria) const ; - bool ProcessCube( int nVoxI, int nVoxJ, int nVoxK, TRIA3DEXLIST& lstTria, TriHolder& triHold, bool bExt) const ; ///// ora + bool ProcessCube( int nVoxI, int nVoxJ, int nVoxK, TRIA3DEXLIST& lstTria, TriHolder& triHold, bool bExt) const ; bool ProcessCell( int nGrid, int nCellI, int nCellJ, const Plane3d& plPlane, std::vector& vLine) const ; - bool ExtMarchingCubes( int nBlock, TRIA3DEXLIST& lstTria, TriHolder& triHold) const ; - bool ExtMarchingCubes( std::vector& vVox, TRIA3DEXLIST& lstTria, bool bEnh) const ; + bool ExtMarchingCubes( int nBlock, TRIA3DEXLIST& lstTria, VoxelContainer& vVox) const ; + bool ExtMarchingCubes( std::vector& vVox, TRIA3DEXLIST& lstTria, bool bEnh) const ; + bool RegulateFeaturesChain( std::vector& vVecVox) const ; + bool CreateSharpFeatureTriangle( int nBlock, const VoxelContainer& vVox, TriHolder& triHold) const ; bool FlipEdgesII( TriHolder& TriHold) const ; bool FlipEdgesBB( TriaMatrix& InterTria) const ; bool IsThereMat( int nI, int nJ, int nK) const ; bool InOut( const Plane3d& plPlane, int nGrid, int nI, int nJ) const ; int CalcIndex( int nI, int nJ, int nK) const ; int CalcIndexForPlaneCells( const Plane3d& plPlane, int nGrid, int nCellI, int nCellJ) const ; - bool IntersPos( int nVec1[], int nVec2[], bool bFirstCorner, VectorField& vfField) const ; + bool IntersPos( int nVec1[], int nVec2[], bool bFirstCorner, AppliedVector& vfField) const ; bool IsPointInsideVoxelApprox( int nI, int nJ, int nK, const Point3d& ptP, double dPrec = EPS_SMALL) const ; bool GetPointVoxel( const Point3d& ptP, int& nVoxI, int& nVoxJ, int& nVoxK) const ; bool IsZInsideInterval( int nGrid, int nDex, double dZ) const ; - + // Per regolarizzazione della catena di feature + bool FindAdjComp( const std::vector& vVecVox, int nCurBlock, int nCurVox, int nCurComp, + INTVECTOR& vAdjBlockVoxComp, INTVECTOR& vAdjBordBlockVoxComp) const ; // OPERAZIONI SU INTERVALLI bool SubtractIntervals( unsigned int nGrid, unsigned int nI, unsigned int nJ, double dMin, double dMax, const Vector3d& vtNMin, const Vector3d& vtNMax) ; bool AddIntervals( unsigned int nGrid, unsigned int nI, unsigned int nJ, double dMin, double dMax, const Vector3d& vtNMin, const Vector3d& vtNMax) ; - // SOTTRAZIONI // UTENSILI // Asse di simmetria parallelo a Z @@ -295,11 +310,11 @@ class VolZmap : public IVolZmap, public IGeoObjRW bool IsATriangleOnBorder( const Triangle3dEx& trTria, const Point3d& ptVert, const int nBlockLimits[], const int nVoxIJK[]) const ; // Funzioni per facce canoniche con grandi triangoli - bool ProcessVoxContXY( VoxelContainer& VoxContXY, bool bPlus, TRIA3DEXLIST& lstTria) const ; - bool ProcessVoxContYZ( VoxelContainer& VoxContYZ, bool bPlus, TRIA3DEXLIST& lstTria) const ; - bool ProcessVoxContXZ( VoxelContainer& VoxContXZ, bool bPlus, TRIA3DEXLIST& lstTria) const ; - bool Find( const VoxelContainer& VoxCont, int nI, int nJ, int nK, double dPos, int nTool) const ; - bool Remove( VoxelContainer& VoxCont, int nI, int nJ, int nK) const ; + bool ProcessVoxContXY( FlatVoxelContainer& VoxContXY, bool bPlus, TRIA3DEXLIST& lstTria) const ; + bool ProcessVoxContYZ( FlatVoxelContainer& VoxContYZ, bool bPlus, TRIA3DEXLIST& lstTria) const ; + bool ProcessVoxContXZ( FlatVoxelContainer& VoxContXZ, bool bPlus, TRIA3DEXLIST& lstTria) const ; + bool Find( const FlatVoxelContainer& VoxCont, int nI, int nJ, int nK, double dPos, int nTool) const ; + bool Remove( FlatVoxelContainer& VoxCont, int nI, int nJ, int nK) const ; // Connessione Zmap struct IntervalIndexes { size_t tMap ; @@ -340,18 +355,18 @@ class VolZmap : public IVolZmap, public IGeoObjRW Vector3d vtMaxN ; int nCompo ; } ; - std::vector> m_Values[N_MAPS] ; // dexel delle 3 griglie unsigned int m_nVoxNumPerBlock ; // Numero di voxel per blocco unsigned int m_nFracLin[3] ; // Numero di blocchi per ogni asse unsigned int m_nNumBlock ; // Numero totale di blocchi - mutable std::vector m_BlockToUpdate ; + mutable BOOLVECTOR m_BlockToUpdate ; int m_nConnectedCompoCount ; // Se == - 1 il numero di componenti non è noto // Se >= 0 è il numero di componenti connesse + mutable std::vector m_InterBlockVox ; mutable TriaMatrix m_InterBlockTria ; Tool m_Tool ; diff --git a/VolZmapCalculus.cpp b/VolZmapCalculus.cpp index 6c5798b..87c1328 100644 --- a/VolZmapCalculus.cpp +++ b/VolZmapCalculus.cpp @@ -420,7 +420,7 @@ VolZmap::GetDepthWithVoxel( const Point3d& ptP, const Vector3d& vtD, double& dIn int nPlaneK = ( vtD.z >= 0 ? 1 : 0) ; // Ciclo sui voxel - vector vVox ; + vector vVox ; while ( IsValidVoxel( nVoxI, nVoxJ, nVoxK)) { // Estremi della diagonale del voxel corrente Point3d ptMin( ( nVoxI * N_DEXVOXRATIO + 0.5) * m_dStep, @@ -431,7 +431,7 @@ VolZmap::GetDepthWithVoxel( const Point3d& ptP, const Vector3d& vtD, double& dIn ( ( nVoxK + 1) * N_DEXVOXRATIO + 0.5) * m_dStep) ; // Studio il voxel corrente if ( IntersLineBox( ptP, vtD, ptMin, ptMax)) { - Voxel NewVox ; + VoxelIndexes NewVox ; NewVox.nI = nVoxI ; NewVox.nJ = nVoxJ ; NewVox.nK = nVoxK ; diff --git a/VolZmapCreation.cpp b/VolZmapCreation.cpp index eab5a9a..30c97f9 100644 --- a/VolZmapCreation.cpp +++ b/VolZmapCreation.cpp @@ -151,9 +151,10 @@ VolZmap::Create( const Point3d& ptO, double dLengthX, double dLengthY, double dL for ( unsigned int nCount = 0 ; nCount < m_nNumBlock ; ++ nCount) m_BlockToUpdate[nCount] = true ; - // Dimensiono raccolta triangoli di feature tra blocchi + // Dimensiono raccolta di voxel e triangoli di feature tra blocchi + m_InterBlockVox.resize( m_nNumBlock) ; m_InterBlockTria.resize( m_nNumBlock) ; - + // Aggiornamento dello stato m_nStatus = OK ; @@ -469,7 +470,8 @@ VolZmap::CreateFromFlatRegion( const ISurfFlatRegion& Surf, double dDimZ, double for ( unsigned int nCount = 0 ; nCount < m_nNumBlock ; ++ nCount) m_BlockToUpdate[nCount] = true ; - // Dimensiono raccolta triangoli di feature tra blocchi + // Dimensiono raccolta di voxel e triangoli di feature tra blocchi + m_InterBlockVox.resize( m_nNumBlock) ; m_InterBlockTria.resize( m_nNumBlock) ; // Aggiornamento dello stato @@ -676,7 +678,8 @@ VolZmap::CreateFromTriMesh( const ISurfTriMesh& Surf, double dPrec, bool bTriDex for ( unsigned int nCount = 0 ; nCount < m_nNumBlock ; ++ nCount) m_BlockToUpdate[nCount] = true ; - // Dimensiono raccolta triangoli di feature tra blocchi + // Dimensiono raccolta di voxel e triangoli di feature tra blocchi + m_InterBlockVox.resize( m_nNumBlock) ; m_InterBlockTria.resize( m_nNumBlock) ; // Aggiornamento dello stato diff --git a/VolZmapGraphics.cpp b/VolZmapGraphics.cpp index 3842efa..e119cb0 100644 --- a/VolZmapGraphics.cpp +++ b/VolZmapGraphics.cpp @@ -47,14 +47,14 @@ enum CanonicDir { X_PLUS = 1, X_MINUS = -1, Y_PLUS = 2, Y_MINUS = -2, Z_PLUS = 3 //---------------------------------------------------------------------------- int -TestOnNormal( const VectorField CompoVert[], int nCompoElem) +TestOnNormal( const AppliedVector CompoField[], int nCompoElem) { // Cerco la massima deviazione tra le normali nei punti della parte connessa int nI, nJ ; double dMinCosTheta = 2 ; for ( int i = 0 ; i < nCompoElem ; ++ i) { for ( int j = i + 1 ; j < nCompoElem ; ++ j) { - double dCurrCos = CompoVert[i].vtNorm * CompoVert[j].vtNorm ; + double dCurrCos = CompoField[i].vtVec * CompoField[j].vtVec ; if ( dCurrCos < dMinCosTheta) { dMinCosTheta = dCurrCos ; nI = i ; @@ -70,12 +70,12 @@ TestOnNormal( const VectorField CompoVert[], int nCompoElem) // Verifico se Edge o Corner // direzione perpendicolare alle normali con massima differenza (potenziale edge) - Vector3d vtK = CompoVert[nI].vtNorm ^ CompoVert[nJ].vtNorm ; + Vector3d vtK = CompoField[nI].vtVec ^ CompoField[nJ].vtVec ; vtK.Normalize() ; // cerco normale con massima vicinanza al potenziale edge double dMaxAbsCos = 0 ; for ( int i = 0 ; i < nCompoElem ; ++ i) { - double dAbsCurrentCos = abs( CompoVert[i].vtNorm * vtK) ; + double dAbsCurrentCos = abs( CompoField[i].vtVec * vtK) ; if ( dAbsCurrentCos > dMaxAbsCos) dMaxAbsCos = dAbsCurrentCos ; } @@ -89,7 +89,7 @@ TestOnNormal( const VectorField CompoVert[], int nCompoElem) //---------------------------------------------------------------------------- bool -CanonicPlaneTest( const VectorField CompoVert[], int nDir, double& dPos, int& nTool) +CanonicPlaneTest( const AppliedVector CompoField[], int nDir, double& dPos, int& nTool) { // Verifico posizione dei punti int nI ; @@ -98,10 +98,10 @@ CanonicPlaneTest( const VectorField CompoVert[], int nDir, double& dPos, int& nT 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] ; + double dPos0 = CompoField[0].ptPApp.v[nI] ; + double dPos1 = CompoField[1].ptPApp.v[nI] ; + double dPos2 = CompoField[2].ptPApp.v[nI] ; + double dPos3 = CompoField[3].ptPApp.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) @@ -110,7 +110,7 @@ CanonicPlaneTest( const VectorField CompoVert[], int nDir, double& dPos, int& nT Vector3d vtN ; switch ( nDir) { case X_PLUS : vtN = X_AX ; break ; - case X_MINUS : 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 ; @@ -118,13 +118,13 @@ CanonicPlaneTest( const VectorField CompoVert[], int nDir, double& dPos, int& nT } int nDifferent = 0 ; for ( int i = 0 ; i < 4 ; ++ i) { - if ( CompoVert[i].vtNorm * vtN < 0.999) + if ( CompoField[i].vtVec * vtN < 0.999) return false ; for ( int j = i + 1 ; j < 4 ; ++ j) { - if ( CompoVert[i].nToolFlag != CompoVert[j].nToolFlag) + if ( CompoField[i].nPropIndex != CompoField[j].nPropIndex) ++ nDifferent ; else - nTool = CompoVert[i].nToolFlag ; + nTool = CompoField[i].nPropIndex ; } } if ( nDifferent > 3) @@ -135,13 +135,13 @@ CanonicPlaneTest( const VectorField CompoVert[], int nDir, double& dPos, int& nT //---------------------------------------------------------------------------- bool -DotTest( const VectorField CompoVert[], int nCompoElem, Vector3d& vtAvg, double dThreshold = 0) +DotTest( const AppliedVector CompoField[], int nCompoElem, Vector3d& vtAvg, double dThreshold = 0) { // Cerco la massima deviazione tra le normali nei punti della parte connessa double dMinCosTheta = 2 ; for ( int i = 0 ; i < nCompoElem ; ++ i) { for ( int j = i + 1 ; j < nCompoElem ; ++ j) { - double dCurrCos = CompoVert[i].vtNorm * CompoVert[j].vtNorm ; + double dCurrCos = CompoField[i].vtVec * CompoField[j].vtVec ; if ( dCurrCos < dMinCosTheta) { dMinCosTheta = dCurrCos ; } @@ -155,12 +155,123 @@ DotTest( const VectorField CompoVert[], int nCompoElem, Vector3d& vtAvg, double // determino media delle normali vtAvg = V_NULL ; for ( int i = 0 ; i < nCompoElem ; ++ i) - vtAvg += CompoVert[i].vtNorm ; + vtAvg += CompoField[i].vtVec ; vtAvg /= nCompoElem ; return true ; } +//---------------------------------------------------------------------------- +bool +VolZmap::FindAdjComp( const std::vector& vVecVox, int nCurBlock, int nCurVox, int nCurComp, + INTVECTOR& vAdjBlockVoxComp, INTVECTOR& vAdjBordBlockVoxComp) const +{ + // Controllo sulla validità del blocco + if ( nCurBlock < 0 || nCurBlock >= int( m_nNumBlock)) + return false ; + // Calcolo gli indici ijk minimi e massimi dei voxel del blocco + // Vettore indici i,j,k del blocco + int nIJK[3] ; + GetBlockIJKFromN( nCurBlock, nIJK) ; + // Vettore limiti sugli indici dei voxel del blocco + int nLimits[6] ; + GetBlockLimitsIJK( nIJK, nLimits) ; + // Voxel corrente + // Determino gli indici ijk del voxel corrente + int nCurIJK[3] ; + GetVoxIJKFromN( nCurVox, nCurIJK[0], nCurIJK[1], nCurIJK[2]) ; + bool bInnerVox = ! IsAVoxelOnBoundary( nLimits, nCurIJK, true) ; + Voxel CurVox = bInnerVox ? ( * vVecVox[nCurBlock].find( nCurVox)).second : + ( * m_InterBlockVox[nCurBlock].find( nCurVox)).second ; + // Ciclo su tutti i voxel adiacenti + for ( int nI = - 1 ; nI <= 1 ; ++ nI) { + for ( int nJ = - 1 ; nJ <= 1 ; ++ nJ) { + for ( int nK = - 1 ; nK <= 1 ; ++ nK) { + // Salto il voxel stesso + if ( nI == 0 && nJ == 0 && nK == 0) + continue ; + // Indici del voxel adiacente + int nAdjVoxIJK[3] = { nCurIJK[0] + nI, nCurIJK[1] + nJ, nCurIJK[2] + nK} ; + // Risalgo all'indice del voxel adiacente (chiave) a partire dagli indici ijk + int nAdjVoxInd ; + if ( GetVoxNFromIJK( nAdjVoxIJK[0], nAdjVoxIJK[1], nAdjVoxIJK[2], nAdjVoxInd)) { + // Risalgo al blocco del voxel adiacente + int nAdjBlockIJK[3] ; + // Se tale blocco esiste + if ( GetVoxelBlockIJK( nAdjVoxIJK, nAdjBlockIJK)) { + // Risalgo all'indice del blocco adiacente dai suoi indici ijk + int nAdjBlockInd ; + GetBlockNFromIJK( nAdjBlockIJK, nAdjBlockInd) ; + // Se tale blocco è stato aggiornato, il voxel adiacente contiene componenti feature ed è interno + if ( m_BlockToUpdate[nAdjBlockInd] && + vVecVox[nAdjBlockInd].find( nAdjVoxInd) != vVecVox[nAdjBlockInd].end()) { + // Accedo al voxel adiacente + Voxel AdjVox = ( * vVecVox[nAdjBlockInd].find( nAdjVoxInd)).second ; + // Valuto la connessione fra le componenti dei due voxels + bool bConnected = false ; + for ( int nAdjComp = 0 ; nAdjComp < AdjVox.nNumComp ; ++ nAdjComp) { + for ( int nCurV = 0 ; nCurV < CurVox.Compo[nCurComp].nVertNum ; ++ nCurV) { + int nNextCurV = ( nCurV + 1) % CurVox.Compo[nCurComp].nVertNum ; + for ( int nAdjV = 0 ; nAdjV < AdjVox.Compo[nAdjComp].nVertNum ; ++ nAdjV) { + int nVextAdjV = ( nAdjV + 1) % AdjVox.Compo[nAdjComp].nVertNum ; + Point3d ptV = CurVox.Compo[nCurComp].CompVecField[nCurV].ptPApp ; + Point3d ptVN = CurVox.Compo[nCurComp].CompVecField[nNextCurV].ptPApp ; + Point3d ptAdV = AdjVox.Compo[nAdjComp].CompVecField[nAdjV].ptPApp ; + Point3d ptAdVN = AdjVox.Compo[nAdjComp].CompVecField[nVextAdjV].ptPApp ; + if ( AreSamePointExact( ptV, ptAdVN) &&/*||*/ AreSamePointExact( ptVN, ptAdV)) { + bConnected = true ; + break ; + } + } + if ( bConnected) + break ; + } + if ( bConnected) { + vAdjBlockVoxComp.emplace_back( nAdjBlockInd) ; + vAdjBlockVoxComp.emplace_back( nAdjVoxInd) ; + vAdjBlockVoxComp.emplace_back( nAdjComp) ; + } + } + } + // A prescindere dal fatto che il blocco sia stato aggiornato, + // se il voxel è di frontiera lo analizzo. + if ( m_InterBlockVox[nAdjBlockInd].find( nAdjVoxInd) != m_InterBlockVox[nAdjBlockInd].end()) { + // Accedo al voxel adiacente + Voxel AdjVox = ( * m_InterBlockVox[nAdjBlockInd].find( nAdjVoxInd)).second ; + // Valuto la connessione fra le componenti dei due voxels + bool bConnected = false ; + for ( int nAdjComp = 0 ; nAdjComp < AdjVox.nNumComp ; ++ nAdjComp) { + for ( int nCurV = 0 ; nCurV < CurVox.Compo[nCurComp].nVertNum ; ++ nCurV) { + int nNextCurV = ( nCurV + 1) % CurVox.Compo[nCurComp].nVertNum ; + for ( int nAdjV = 0 ; nAdjV < AdjVox.Compo[nAdjComp].nVertNum ; ++ nAdjV) { + int nVextAdjV = ( nAdjV + 1) % AdjVox.Compo[nAdjComp].nVertNum ; + Point3d ptV = CurVox.Compo[nCurComp].CompVecField[nCurV].ptPApp ; + Point3d ptVN = CurVox.Compo[nCurComp].CompVecField[nNextCurV].ptPApp ; + Point3d ptAdV = AdjVox.Compo[nAdjComp].CompVecField[nAdjV].ptPApp ; + Point3d ptAdVN = AdjVox.Compo[nAdjComp].CompVecField[nVextAdjV].ptPApp ; + if ( AreSamePointExact( ptV, ptAdVN) &&/*||*/ AreSamePointExact( ptVN, ptAdV)) { + bConnected = true ; + break ; + } + } + if ( bConnected) + break ; + } + if ( bConnected) { + vAdjBordBlockVoxComp.emplace_back( nAdjBlockInd) ; + vAdjBordBlockVoxComp.emplace_back( nAdjVoxInd) ; + vAdjBordBlockVoxComp.emplace_back( nAdjComp) ; + } + } + } + } + } + } + } + } + return true ; +} + // ------------------------- VISUALIZZAZIONE -------------------------------------------------------------------------------------- //---------------------------------------------------------------------------- bool @@ -337,6 +448,9 @@ VolZmap::GetTriangles( bool bAllBlocks, INTVECTOR& nModifiedBlocks, TRIA3DEXLIST nModifiedBlocks.resize( m_nNumBlock + 1) ; vLstTria.reserve( m_nNumBlock + 1) ; + //std::vector> VoxMatrix ; + std::vector vVoxContainerVec ; + vVoxContainerVec.resize( m_nNumBlock) ; TriaMatrix VecTriHold ; VecTriHold.resize( m_nNumBlock) ; @@ -344,28 +458,36 @@ VolZmap::GetTriangles( bool bAllBlocks, INTVECTOR& nModifiedBlocks, TRIA3DEXLIST // Calcolo i triangoli sui blocchi for ( size_t t = 0 ; t < m_nNumBlock ; ++ t) { - // Se il blocco deve essere processato if ( bAllBlocks || m_BlockToUpdate[t]) { - // processo ... vLstTria.emplace_back() ; - nModifiedBlocks[t] = int( vLstTria.size()) - 1 ; - m_InterBlockTria[t].clear() ; - ExtMarchingCubes( int( t), vLstTria.back(), VecTriHold[t]) ; - // Flipping fra voxel interni - FlipEdgesII( VecTriHold[t]) ; - bCalcInterBlock = true ; - m_BlockToUpdate[t] = false ; + nModifiedBlocks[t] = int( vLstTria.size()) - 1 ; + ExtMarchingCubes( int( t), vLstTria.back(), vVoxContainerVec[t]) ; } else nModifiedBlocks[t] = -1 ; } + // Regolarizzo la catena + RegulateFeaturesChain( vVoxContainerVec) ; + + // Costruisco i triangoli di feature + for ( size_t t = 0 ; t < m_nNumBlock ; ++ t) { + // Se il blocco è stato processato + if ( bAllBlocks || m_BlockToUpdate[t]) { + CreateSharpFeatureTriangle( int( t), vVoxContainerVec[t], VecTriHold[t]) ; + // Flipping fra voxel interni + FlipEdgesII( VecTriHold[t]) ; + bCalcInterBlock = true ; + m_BlockToUpdate[t] = false ; + } + } + // Calcolo i triangoli di frontiera tra feature di blocchi diversi - // copio i triangoli di frontiera in una matrice gemella - // di m_InterBlockTria per avere sempre a disposizione - // i triangoli non flippati. + // copio i triangoli di frontiera in una matrice gemella + // di m_InterBlockTria per avere sempre a disposizione + // i triangoli non flippati. TriaMatrix InterBlockTria ; if ( bCalcInterBlock) { InterBlockTria = m_InterBlockTria ; @@ -402,7 +524,6 @@ VolZmap::GetTriangles( bool bAllBlocks, INTVECTOR& nModifiedBlocks, TRIA3DEXLIST if ( bCalcInterBlock) { vLstTria.resize( vLstTria.size() + 1) ; size_t nPos = size_t( vLstTria.size() - 1) ; - for ( size_t t = 0 ; t < InterBlockTria.size() ; ++ t) { for ( size_t t1 = 0 ; t1 < InterBlockTria[t].size() ; ++ t1) { for ( size_t t2 = 0 ; t2 < InterBlockTria[t][t1].vCompoTria.size() ; ++ t2) { @@ -424,10 +545,8 @@ VolZmap::GetTriangles( bool bAllBlocks, INTVECTOR& nModifiedBlocks, TRIA3DEXLIST } } } - nModifiedBlocks.back() = int( nPos) ; } - else nModifiedBlocks.back() = - 1 ; } @@ -676,7 +795,7 @@ VolZmap::ProcessCube( int nVoxI, int nVoxJ, int nVoxK, TRIA3DEXLIST& lstTria, Tr } ; // Array di strutture punto di intersezione e normale alla superficie in esso. - VectorField VecField[12] ; + AppliedVector VecField[12] ; // Flag di regolarità dei campi scalare e vettoriale bool bReg = true ; @@ -704,16 +823,13 @@ VolZmap::ProcessCube( int nVoxI, int nVoxJ, int nVoxK, TRIA3DEXLIST& lstTria, Tr // 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] ; + AppliedVector CompoVert[4][7] ; + AppliedVector CompoTriVert[4][17] ; // Arrey numero di vertici della base del fan per componente // connessa: nVertComp[i] contiene il numero di vertici // della base del fan della (i+1)-esima componente connessa. - int nVertComp[6] ; - - // Matrice di indici dei punti: serve per la gestione del caso - int nIndArrey[6][4] ; + int nVertComp[4] ; int nExtTabOff = nComponents ; int nStdTabOff = 0 ; @@ -727,13 +843,6 @@ VolZmap::ProcessCube( int nVoxI, int nVoxJ, int nVoxK, TRIA3DEXLIST& lstTria, Tr // Riempio il nCompCount-esimo vettore di vertici della base del fan for ( int nVertCount = 0 ; nVertCount < nVertComp[nComp] ; ++ nVertCount) CompoVert[nComp][nVertCount] = VecField[TriangleTableEn[nIndex][1][nVertCount + nExtTabOff + 1]] ; - - // Serve per la gestione del caso ... - if ( nVertComp[nComp] == 4) { - for ( int nVert = 0 ; nVert < nVertComp[nComp] ; ++ nVert) - nIndArrey[nComp][nVert] = TriangleTableEn[nIndex][1][nVert + 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 nVert = 0 ; nVert < 3 * ( nVertComp[nComp] - 2) ; nVert += 3) { @@ -751,19 +860,14 @@ VolZmap::ProcessCube( int nVoxI, int nVoxJ, int nVoxK, TRIA3DEXLIST& lstTria, Tr // 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. - // Il flag bSecondConfig6 serve perché nel caso di configurazione 6 - // con una sola componente connessa non si deve eseguire l'Extended MC. - bool bSecondConfig6 = false ; if ( bReg) { if ( nAllConfig[nIndex] == 3) { - - Vector3d vtCmpAvg0, vtCmpAvg1 ; - // Verifico se i versori delle componenti sono tutti // più o meno concordi (per esserlo non devono esserci // due vettori di una medesima componente con prodotto // scalare inferiore a 0.7). + Vector3d vtCmpAvg0, vtCmpAvg1 ; bool bTest0 = DotTest( CompoVert[0], 3, vtCmpAvg0, 0.7) ; bool bTest1 = DotTest( CompoVert[1], 3, vtCmpAvg1, 0.7) ; @@ -775,27 +879,24 @@ VolZmap::ProcessCube( int nVoxI, int nVoxJ, int nVoxK, TRIA3DEXLIST& lstTria, Tr // (il prodotto scalare fra versori ha modulo non superiore // a uno). double dScProd = - 2 ; - if ( bTest0 && bTest1) dScProd = vtCmpAvg0 * vtCmpAvg1 ; double dThreshold = 0.7 ; - + // Si passa alla seconda topologia if ( ! ( bTest0 && bTest1) || ( bTest0 && bTest1 && dScProd > dThreshold)) { - + // Ricerca del caso corrispondente della nuova topologia int nt = 0 ; while ( nIndexVsIndex3[nt][0] != nIndex) ++ nt ; int nRotCase = nIndexVsIndex3[nt][1] ; - + // Riaggiorno numero di componenti nComponents = Cases3Plus[nRotCase][1][0] ; - // Riaggiorno gli offsets nExtTabOff = nComponents ; nStdTabOff = 0 ; - // Modifico le matrici for ( int nC = 1 ; nC <= nComponents ; ++ nC) { @@ -805,15 +906,12 @@ VolZmap::ProcessCube( int nVoxI, int nVoxJ, int nVoxK, TRIA3DEXLIST& lstTria, Tr // Matrice dei vertici della base del fan for ( int nFanVert = 0 ; nFanVert < nVertComp[nC - 1] ; ++ nFanVert) CompoVert[nC - 1][nFanVert] = VecField[Cases3Plus[nRotCase][1][nFanVert + nExtTabOff + 1]] ; - // 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) ; @@ -825,93 +923,74 @@ VolZmap::ProcessCube( int nVoxI, int nVoxJ, int nVoxK, TRIA3DEXLIST& lstTria, Tr // 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) dScProd = vtCmpAvg0 * vtCmpAvg1 ; double dThreshold = 0.7 ; - + // Si passa alla seconda topologia if ( ! ( bTest0 && bTest1) || ( bTest0 && bTest1 && dScProd > dThreshold)) { - + // Ricerca del caso corrispondente della nuova topologia int nt = 0 ; - while ( nIndexVsIndex6[nt][0] != nIndex) ++ nt ; - int nRotCase = nIndexVsIndex6[nt][1] ; - + // Riaggiorno numero di componenti nComponents = 1 ; - // Riaggiorno gli offsets nExtTabOff = nComponents ; nStdTabOff = 0 ; - // Modifico le matrici for ( int nC = 1 ; nC <= nComponents ; ++ nC) { - // Numero vertici per componenti nVertComp[nC - 1] = 7 ; - // 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) ; } - bSecondConfig6 = true ; + return true ; } } else if ( nAllConfig[nIndex] == 10) { - - Vector3d vtCmpAvg0, vtCmpAvg1 ; - // Verifico se i versori delle componenti sono tutti // più o meno concordi (per esserlo non devono esserci // due vettori di una medesima componente con prodotto // scalare inferiore a 0). decidere se 0.0 o 0.7 + Vector3d vtCmpAvg0, vtCmpAvg1 ; bool bTest0 = DotTest( CompoVert[0], 4, vtCmpAvg0) ; bool bTest1 = DotTest( CompoVert[1], 4, vtCmpAvg1) ; - + // Si passa alla seconda topologia if ( ! ( bTest0 && bTest1)) { - + // Ricerca del caso corrispondente della nuova topologia int nt = 0 ; while ( nIndexVsIndex10[nt][0] != nIndex) ++ nt ; - int nRotCase = nIndexVsIndex10[nt][1] ; - // Riaggiorno gli offsets nExtTabOff = 2 ; nStdTabOff = 0 ; - // 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. nExtTabOff += nVertComp[nC - 1] ; nStdTabOff += 3 * ( nVertComp[nC - 1] - 2) ; @@ -933,7 +1012,7 @@ VolZmap::ProcessCube( int nVoxI, int nVoxJ, int nVoxK, TRIA3DEXLIST& lstTria, Tr nFeatureType = TestOnNormal( CompoVert[nComp], nVertComp[nComp]) ; // Flag ExtMC - bool bExtMC = ( nFeatureType != NO_FEATURE) && ( ! bSecondConfig6) && bEnh ; + bool bExtMC = ( nFeatureType != NO_FEATURE) && bEnh ; // Extended MC if ( bExtMC) { @@ -941,12 +1020,12 @@ VolZmap::ProcessCube( int nVoxI, int nVoxJ, int nVoxK, TRIA3DEXLIST& lstTria, Tr // Passo al sistema di riferimento del baricentro Point3d ptGravityCenter( 0, 0, 0) ; for ( int ni = 0 ; ni < nVertComp[nComp] ; ++ ni) - ptGravityCenter = ptGravityCenter + CompoVert[nComp][ni].ptInt ; + ptGravityCenter = ptGravityCenter + CompoVert[nComp][ni].ptPApp ; ptGravityCenter = ptGravityCenter / nVertComp[nComp] ; Vector3d vtTrasf[12] ; for ( int ni = 0 ; ni < nVertComp[nComp] ; ++ ni) - vtTrasf[ni] = CompoVert[nComp][ni].ptInt - ptGravityCenter ; + vtTrasf[ni] = CompoVert[nComp][ni].ptPApp - ptGravityCenter ; // Preparo le matrici per il sistema typedef Eigen::Matrix dSystemMatrix ; @@ -958,10 +1037,10 @@ VolZmap::ProcessCube( int nVoxI, int nVoxJ, int nVoxK, TRIA3DEXLIST& lstTria, Tr // Caso generale for ( int ni = 0 ; ni < nVertComp[nComp] ; ++ ni) { - dMatrixN( ni, 0) = CompoVert[nComp][ni].vtNorm.x ; - dMatrixN( ni, 1) = CompoVert[nComp][ni].vtNorm.y ; - dMatrixN( ni, 2) = CompoVert[nComp][ni].vtNorm.z ; - dKnownVector( ni) = CompoVert[nComp][ni].vtNorm * vtTrasf[ni] ; + dMatrixN( ni, 0) = CompoVert[nComp][ni].vtVec.x ; + dMatrixN( ni, 1) = CompoVert[nComp][ni].vtVec.y ; + dMatrixN( ni, 2) = CompoVert[nComp][ni].vtVec.z ; + dKnownVector( ni) = CompoVert[nComp][ni].vtVec * vtTrasf[ni] ; } // calcolo SVD @@ -989,8 +1068,10 @@ VolZmap::ProcessCube( int nVoxI, int nVoxJ, int nVoxK, TRIA3DEXLIST& lstTria, Tr // nel caso essa esca dal reticolo la limito entro una // distanza dal baricentro pari alla diagonale del voxel. bool bOutside = false ; + int nFtVxI, nFtVxJ, nFtVxK ; if ( GetPointVoxel( ptSol, nFtVxI, nFtVxJ, nFtVxK)) { + if ( abs( nFtVxI - nVoxI) > 2 || abs( nFtVxJ - nVoxJ) > 2 || abs( nFtVxK - nVoxK) > 2) @@ -1008,7 +1089,7 @@ VolZmap::ProcessCube( int nVoxI, int nVoxJ, int nVoxK, TRIA3DEXLIST& lstTria, Tr int nj = ( ni + 1 < nVertComp[nComp]) ? ni + 1 : 0 ; // Il triangolo è pronto Triangle3d CurrentTriangle ; - CurrentTriangle.Set( ptSol, CompoVert[nComp][nj].ptInt, CompoVert[nComp][ni].ptInt) ; + CurrentTriangle.Set( ptSol, CompoVert[nComp][nj].ptPApp, CompoVert[nComp][ni].ptPApp) ; CurrentTriangle.Validate( true) ; // Aggiungo triangolo al vettore temporaneo triContainer.emplace_back( CurrentTriangle) ; @@ -1027,8 +1108,8 @@ VolZmap::ProcessCube( int nVoxI, int nVoxJ, int nVoxK, TRIA3DEXLIST& lstTria, Tr int nj = ( ni + 1 < nVertComp[nComp]) ? ni + 1 : 0 ; - double dDI = triContainer[ni].GetN() * CompoVert[nComp][ni].vtNorm ; - double dDJ = triContainer[ni].GetN() * CompoVert[nComp][nj].vtNorm ; + double dDI = triContainer[ni].GetN() * CompoVert[nComp][ni].vtVec ; + double dDJ = triContainer[ni].GetN() * CompoVert[nComp][nj].vtVec ; if ( dDI < - EPS_SMALL || dDJ < - EPS_SMALL) { bInversione = true ; @@ -1043,7 +1124,7 @@ VolZmap::ProcessCube( int nVoxI, int nVoxJ, int nVoxK, TRIA3DEXLIST& lstTria, Tr int nNegDot = 0 ; for ( int ni = 0 ; ni < nVertComp[nComp] - 1 ; ++ ni) { for ( int nj = ni + 1 ; nj < nVertComp[nComp] ; ++ nj) { - if ( CompoVert[nComp][ni].vtNorm * CompoVert[nComp][nj].vtNorm < - EPS_SMALL) + if ( CompoVert[nComp][ni].vtVec * CompoVert[nComp][nj].vtVec < - EPS_SMALL) nNegDot ++ ; } } @@ -1054,8 +1135,8 @@ VolZmap::ProcessCube( int nVoxI, int nVoxJ, int nVoxK, TRIA3DEXLIST& lstTria, Tr bool bInversione2 = false ; for ( int ni = 0 ; ni < nVertComp[nComp] ; ++ ni) { int nj = ( ni == 0 ? nVertComp[nComp] - 1 : ni - 1) ; - double dDLast = triContainer[ni].GetN() * CompoVert[nComp][ni].vtNorm ; - double dDPrev = triContainer[nj].GetN() * CompoVert[nComp][ni].vtNorm ; + double dDLast = triContainer[ni].GetN() * CompoVert[nComp][ni].vtVec ; + double dDPrev = triContainer[nj].GetN() * CompoVert[nComp][ni].vtVec ; if ( ( dDLast < EPS_SMALL && dDPrev < EPS_SMALL) || ( dDLast < - 0.75 || dDPrev < - 0.75)) { bInversione2 = true ; @@ -1097,7 +1178,7 @@ VolZmap::ProcessCube( int nVoxI, int nVoxJ, int nVoxK, TRIA3DEXLIST& lstTria, Tr int nj = ( ni + 1 < nVertComp[nComp]) ? ni + 1 : 0 ; // Il triangolo è pronto Triangle3d CurrentTriangle ; - CurrentTriangle.Set( ptSol, CompoVert[nComp][nj].ptInt, CompoVert[nComp][ni].ptInt) ; + CurrentTriangle.Set( ptSol, CompoVert[nComp][nj].ptPApp, CompoVert[nComp][ni].ptPApp) ; CurrentTriangle.Validate( true) ; // Aggiungo triangolo al vettore temporaneo triContainer.emplace_back( CurrentTriangle) ; @@ -1108,8 +1189,10 @@ VolZmap::ProcessCube( int nVoxI, int nVoxJ, int nVoxK, TRIA3DEXLIST& lstTria, Tr // 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 = CalcIndex( nAdjVoxI, nAdjVoxJ, nAdjVoxK) ; // Se il voxel è pieno @@ -1132,8 +1215,8 @@ VolZmap::ProcessCube( int nVoxI, int nVoxJ, int nVoxK, TRIA3DEXLIST& lstTria, Tr int nj = ( ni + 1 < nVertComp[nComp]) ? ni + 1 : 0 ; - double dDI = triContainer[ni].GetN() * CompoVert[nComp][ni].vtNorm ; - double dDJ = triContainer[ni].GetN() * CompoVert[nComp][nj].vtNorm ; + double dDI = triContainer[ni].GetN() * CompoVert[nComp][ni].vtVec ; + double dDJ = triContainer[ni].GetN() * CompoVert[nComp][nj].vtVec ; if ( dDI < - EPS_SMALL || dDJ < - EPS_SMALL) bInversione = true ; @@ -1143,6 +1226,7 @@ VolZmap::ProcessCube( int nVoxI, int nVoxJ, int nVoxK, TRIA3DEXLIST& lstTria, Tr if ( bInversione) { // Se la feature non cade nel suo voxel if ( ! IsPointInsideVoxelApprox( nVoxI, nVoxJ, nVoxK, ptSol)) { + Vector3d vtNullSpace( dMatrixV( 0, 2), dMatrixV( 1, 2), dMatrixV( 2, 2)) ; double dParInt1, dParInt2 ; Point3d ptVoxMin( ( nVoxI * N_DEXVOXRATIO + 0.5) * m_dStep, @@ -1152,17 +1236,21 @@ VolZmap::ProcessCube( int nVoxI, int nVoxJ, int nVoxK, TRIA3DEXLIST& lstTria, Tr // Se è possibile riportarla dentro e il voxel in cui cade è pieno, la riporto nel suo voxel // lungo la sua linea if ( IntersLineBox( ptSol, vtNullSpace, ptVoxMin, ptVoxMax, dParInt1, dParInt2)) { + triContainer.resize( 0) ; + double dPar = abs( dParInt1) < abs( dParInt2) ? dParInt1 + ( dParInt2 - dParInt1) / 100 : dParInt2 + ( dParInt1 - dParInt2) / 100 ; + Point3d ptNewSol = ptSol + dPar * vtNullSpace ; ptSol = ptNewSol ; + // Costruisco triangoli di prova for ( int ni = 0 ; ni < nVertComp[nComp] ; ++ ni) { int nj = ( ni + 1 < nVertComp[nComp]) ? ni + 1 : 0 ; // Il triangolo è pronto Triangle3d CurrentTriangle ; - CurrentTriangle.Set( ptSol, CompoVert[nComp][nj].ptInt, CompoVert[nComp][ni].ptInt) ; + CurrentTriangle.Set( ptSol, CompoVert[nComp][nj].ptPApp, CompoVert[nComp][ni].ptPApp) ; CurrentTriangle.Validate( true) ; // Aggiungo triangolo al vettore temporaneo triContainer.emplace_back( CurrentTriangle) ; @@ -1255,9 +1343,9 @@ VolZmap::ProcessCube( int nVoxI, int nVoxJ, int nVoxK, TRIA3DEXLIST& lstTria, Tr for ( int TriIndex = 0; TriIndex < ( nVertComp[nComp] - 2) * 3 ; TriIndex += 3) { // Il triangolo è pronto Triangle3d CurrentTriangle ; - CurrentTriangle.Set( CompoTriVert[nComp][TriIndex].ptInt, - CompoTriVert[nComp][TriIndex+1].ptInt, - CompoTriVert[nComp][TriIndex+2].ptInt) ; + CurrentTriangle.Set( CompoTriVert[nComp][TriIndex].ptPApp, + CompoTriVert[nComp][TriIndex+1].ptPApp, + CompoTriVert[nComp][TriIndex+2].ptPApp) ; bool bV = CurrentTriangle.Validate( true) ; // Aggiungo alla lista lstTria.emplace_back( CurrentTriangle) ; @@ -1271,9 +1359,9 @@ VolZmap::ProcessCube( int nVoxI, int nVoxJ, int nVoxK, TRIA3DEXLIST& lstTria, Tr for ( int TriIndex = 0; TriIndex < ( nVertComp[nComp] - 2) * 3 ; TriIndex += 3) { // Il triangolo è pronto Triangle3d CurrentTriangle ; - CurrentTriangle.Set( CompoTriVert[nComp][TriIndex].ptInt, - CompoTriVert[nComp][TriIndex+1].ptInt, - CompoTriVert[nComp][TriIndex+2].ptInt) ; + CurrentTriangle.Set( CompoTriVert[nComp][TriIndex].ptPApp, + CompoTriVert[nComp][TriIndex+1].ptPApp, + CompoTriVert[nComp][TriIndex+2].ptPApp) ; bool bV = CurrentTriangle.Validate( true) ; // Aggiungo alla lista lstTria.emplace_back( CurrentTriangle) ; @@ -1286,33 +1374,39 @@ VolZmap::ProcessCube( int nVoxI, int nVoxJ, int nVoxK, TRIA3DEXLIST& lstTria, Tr //---------------------------------------------------------------------------- bool -VolZmap::ExtMarchingCubes( int nBlock, TRIA3DEXLIST& lstTria, TriHolder& triHold) const +VolZmap::ExtMarchingCubes( int nBlock, TRIA3DEXLIST& lstTria, VoxelContainer& vVox) const { // Controllo sulla validità del blocco if ( nBlock < 0 || nBlock >= int( m_nNumBlock)) return false ; - // Calcolo i limiti sugli indici dei voxel del blocco - int nIJK[3] ; + // Calcolo i limiti sugli indici dei voxel del blocco // Vettore indici i,j,k del blocco + int nIJK[3] ; GetBlockIJKFromN( nBlock, nIJK) ; // Vettore limiti sugli indici dei voxel del blocco int nLimits[6] ; GetBlockLimitsIJK( nIJK, nLimits) ; + // Pulisco il contenitore dei voxel di frontiera + m_InterBlockVox[nBlock].clear() ; + // Unordered Map per la riduzione del numero di triangoli int nDim = m_nVoxNumPerBlock * m_nVoxNumPerBlock ; - VoxelContainer VoxContXZInf( nDim) ; - VoxelContainer VoxContXZSup( nDim) ; - VoxelContainer VoxContXYInf( nDim) ; - VoxelContainer VoxContXYSup( nDim) ; - VoxelContainer VoxContYZInf( nDim) ; - VoxelContainer VoxContYZSup( nDim) ; - + FlatVoxelContainer VoxContXZInf( nDim) ; + FlatVoxelContainer VoxContXZSup( nDim) ; + FlatVoxelContainer VoxContXYInf( nDim) ; + FlatVoxelContainer VoxContXYSup( nDim) ; + FlatVoxelContainer VoxContYZInf( nDim) ; + FlatVoxelContainer VoxContYZSup( nDim) ; + // Ciclo su tutti i voxel del blocco 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) { + + /* if ( ! ( j >= 2 && j <= 8) ) + continue ;*/ // Classificazione dei vertici: interni o esterni al materiale int nIndex = CalcIndex( i, j, k) ; @@ -1340,7 +1434,7 @@ VolZmap::ExtMarchingCubes( int nBlock, TRIA3DEXLIST& lstTria, TriHolder& triHold } ; // Array di strutture punto di intersezione e normale alla superficie in esso. - VectorField VecField[12] ; + AppliedVector VecField[12] ; // Flag di regolarità dei campi scalare e vettoriale bool bReg = true ; @@ -1368,36 +1462,25 @@ VolZmap::ExtMarchingCubes( int nBlock, TRIA3DEXLIST& lstTria, TriHolder& triHold // 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] ; + AppliedVector CompoVert[4][7] ; + AppliedVector CompoTriVert[4][17] ; // Arrey numero di vertici della base del fan per componente // connessa: nVertComp[i] contiene il numero di vertici // della base del fan della (i+1)-esima componente connessa. - int nVertComp[6] ; - - // Matrice di indici dei punti: serve per la gestione del caso - int nIndArrey[6][4] ; + int nVertComp[4] ; int nExtTabOff = nComponents ; int nStdTabOff = 0 ; // Carico le matrici CompoVert e CompoTriVert for ( int nComp = 0 ; nComp < nComponents ; ++ nComp) { - // Numero vertici per componenti nVertComp[nComp] = TriangleTableEn[nIndex][1][nComp+1] ; // Riempio il nCompCount-esimo vettore di vertici della base del fan for ( int nVert = 0 ; nVert < nVertComp[nComp] ; ++ nVert) CompoVert[nComp][nVert] = VecField[TriangleTableEn[nIndex][1][nVert + nExtTabOff + 1]] ; - - // Serve per la gestione del caso ... - if ( nVertComp[nComp] == 4) { - for ( int nVert = 0 ; nVert < nVertComp[nComp] ; ++ nVert) - nIndArrey[nComp][nVert] = TriangleTableEn[nIndex][1][nVert + 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 nVert = 0 ; nVert < 3 * ( nVertComp[nComp] - 2) ; nVert += 3) { @@ -1469,35 +1552,30 @@ VolZmap::ExtMarchingCubes( int nBlock, TRIA3DEXLIST& lstTria, TriHolder& triHold // Test sulla topologia: richiede campo regolare perchè si fonda su angoli tra normali if ( bReg) { - // Configurazione 3 if ( nAllConfig[nIndex] == 3) { - // Verifico concordanza tra i versori di una stessa componente // (ogni coppia di vettori di una medesima componente deve avere prodotto scalare non inferiore a 0.7) Vector3d vtCmpAvg0, vtCmpAvg1 ; bool bTest0 = DotTest( CompoVert[0], 3, vtCmpAvg0, 0.7) ; bool bTest1 = DotTest( CompoVert[1], 3, vtCmpAvg1, 0.7) ; - // Se i versori di entrambe le componenti sono concordi ha senso parlare di vettori medi. // Altrimenti assegno al loro prodotto scalare un valore assurdo come -2 double dScProd = - 2 ; if ( bTest0 && bTest1) dScProd = vtCmpAvg0 * vtCmpAvg1 ; - + // Si passa alla seconda topologia if ( ! ( bTest0 && bTest1) || ( bTest0 && bTest1 && dScProd > 0.7)) { - + // Ricerca del caso corrispondente della nuova topologia int nt = 0 ; while ( nIndexVsIndex3[nt][0] != nIndex) ++ nt ; int nRotCase = nIndexVsIndex3[nt][1] ; - + // Aggiorno numero di componenti nComponents = Cases3Plus[nRotCase][1][0] ; - // Riaggiorno gli offsets nExtTabOff = nComponents ; nStdTabOff = 0 ; - // Modifico le matrici for ( int nC = 1 ; nC <= nComponents ; ++ nC) { // Numero vertici per componenti @@ -1520,36 +1598,27 @@ VolZmap::ExtMarchingCubes( int nBlock, TRIA3DEXLIST& lstTria, TriHolder& triHold // Configurazione 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) - dScProd = vtCmpAvg0 * vtCmpAvg1 ; - bool bFirstTest = ! ( bTest0 && bTest1) || ( bTest0 && bTest1 && dScProd > 0.7) ;*/ - // Altro metodo + // Test sulla topologia double dDotSum = 0 ; for ( int nFV = 0 ; nFV < 4 ; ++ nFV) { for ( int nTV = 0 ; nTV < 3 ; ++ nTV) { - dDotSum += CompoVert[0][nFV].vtNorm * CompoVert[1][nTV].vtNorm ; + dDotSum += CompoVert[0][nFV].vtVec * CompoVert[1][nTV].vtVec ; } } for ( int nFVI = 0 ; nFVI < 3 ; ++ nFVI) { for ( int nFVJ = nFVI + 1 ; nFVJ < 4 ; ++ nFVJ) { - dDotSum -= CompoVert[0][nFVI].vtNorm * CompoVert[0][nFVJ].vtNorm ; + dDotSum -= CompoVert[0][nFVI].vtVec * CompoVert[0][nFVJ].vtVec ; } } for ( int nTVI = 0 ; nTVI < 2 ; ++ nTVI) { for ( int nTVJ = nTVI + 1 ; nTVJ < 3 ; ++ nTVJ) { - dDotSum -= CompoVert[1][nTVI].vtNorm * CompoVert[1][nTVJ].vtNorm ; + dDotSum -= CompoVert[1][nTVI].vtVec * CompoVert[1][nTVJ].vtVec ; } } bool bTestOnSum = dDotSum > - 4 ; - + // Si deve passare alla seconda topologia if ( bTestOnSum) { - + // Ricerca del caso corrispondente della nuova topologia int nt = 0 ; while ( nIndexVsIndex6[nt][0] != nIndex) ++ nt ; @@ -1563,32 +1632,32 @@ VolZmap::ExtMarchingCubes( int nBlock, TRIA3DEXLIST& lstTria, TriHolder& triHold int i2 = Cases6Plus[nRotCase][TriIndex] ; // Costruzione triangolo Triangle3dEx CurrentTriangle ; - CurrentTriangle.Set( VecField[i0].ptInt, VecField[i1].ptInt, VecField[i2].ptInt) ; + CurrentTriangle.Set( VecField[i0].ptPApp, VecField[i1].ptPApp, VecField[i2].ptPApp) ; CurrentTriangle.Validate( true) ; - CurrentTriangle.SetVertexNorm( 0, VecField[i0].vtNorm) ; - CurrentTriangle.SetVertexNorm( 1, VecField[i1].vtNorm) ; - CurrentTriangle.SetVertexNorm( 2, VecField[i2].vtNorm) ; + CurrentTriangle.SetVertexNorm( 0, VecField[i0].vtVec) ; + CurrentTriangle.SetVertexNorm( 1, VecField[i1].vtVec) ; + CurrentTriangle.SetVertexNorm( 2, VecField[i2].vtVec) ; // Setto il numero di utensile (conta solo positivo, nullo o negativo) - int nTool0 = Clamp( VecField[i0].nToolFlag, -1, 1) ; - int nTool1 = Clamp( VecField[i1].nToolFlag, -1, 1) ; - int nTool2 = Clamp( VecField[i2].nToolFlag, -1, 1) ; + int nTool0 = Clamp( VecField[i0].nPropIndex, -1, 1) ; + int nTool1 = Clamp( VecField[i1].nPropIndex, -1, 1) ; + int nTool2 = Clamp( VecField[i2].nPropIndex, -1, 1) ; if ( nTool0 == nTool1 || nTool0 == nTool2) CurrentTriangle.SetGrade( nTool0) ; else if ( nTool1 == nTool2) CurrentTriangle.SetGrade( nTool1) ; // Valido il triangolo e setto le normali del campo vettoriale ai corrispondenti vertici if ( CurrentTriangle.Validate( true)) { - Vector3d vtVertNorm = VecField[i0].vtNorm ; + Vector3d vtVertNorm = VecField[i0].vtVec ; if ( CurrentTriangle.GetN() * vtVertNorm > 0.6) CurrentTriangle.SetVertexNorm( 0, vtVertNorm) ; else CurrentTriangle.SetVertexNorm( 0, CurrentTriangle.GetN()) ; - vtVertNorm = VecField[i1].vtNorm ; + vtVertNorm = VecField[i1].vtVec ; if ( CurrentTriangle.GetN() * vtVertNorm > 0.6) CurrentTriangle.SetVertexNorm( 1, vtVertNorm) ; else CurrentTriangle.SetVertexNorm( 1, CurrentTriangle.GetN()) ; - vtVertNorm = VecField[i2].vtNorm ; + vtVertNorm = VecField[i2].vtVec ; if ( CurrentTriangle.GetN() * vtVertNorm > 0.6) CurrentTriangle.SetVertexNorm( 2, vtVertNorm) ; else @@ -1602,42 +1671,36 @@ VolZmap::ExtMarchingCubes( int nBlock, TRIA3DEXLIST& lstTria, TriHolder& triHold continue ; } } - // Configurazione 10 else if ( nAllConfig[nIndex] == 10) { - // Verifico concordanza tra i versori di una stessa componente // (ogni coppia di vettori di una medesima componente deve avere prodotto scalare non inferiore a 0.0) Vector3d vtCmpAvg0, vtCmpAvg1 ; bool bTest0 = DotTest( CompoVert[0], 4, vtCmpAvg0, 0.0) ; bool bTest1 = DotTest( CompoVert[1], 4, vtCmpAvg1, 0.0) ; - + // Si passa alla seconda topologia if ( ! bTest0 || ! bTest1) { + // Ricerca del caso corrispondente della nuova topologia int nt = 0 ; while ( nIndexVsIndex10[nt][0] != nIndex) ++ nt ; - // Riaggiorno gli offsets nExtTabOff = 2 ; nStdTabOff = 0 ; - // Modifico le matrici int nRotCase = nIndexVsIndex10[nt][1] ; 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. nExtTabOff += nVertComp[nC - 1] ; nStdTabOff += 3 * ( nVertComp[nC - 1] - 2) ; @@ -1646,9 +1709,11 @@ VolZmap::ExtMarchingCubes( int nBlock, TRIA3DEXLIST& lstTria, TriHolder& triHold } } + Voxel VoxConf ; + VoxConf.nNumComp = 0 ; + // Numero di feature nel voxel: al più vi è una feature per componente connessa. - int nInnerFeatureInVoxel = 0 ; - int nBorderFeatureInVoxel = 0 ; + int nFeatureInVoxel = 0 ; // Ciclo sulle componenti for ( int nComp = 0 ; nComp < nComponents ; ++ nComp) { @@ -1658,21 +1723,18 @@ VolZmap::ExtMarchingCubes( int nBlock, TRIA3DEXLIST& lstTria, TriHolder& triHold if ( bReg) nFeatureType = TestOnNormal( CompoVert[nComp], nVertComp[nComp]) ; - // Flag ExtMC - bool bExtMC = ( nFeatureType != NO_FEATURE) ; - // Extended MC - if ( bExtMC) { + if ( nFeatureType != NO_FEATURE) { // Passo al sistema di riferimento del baricentro Point3d ptGravityCenter( 0, 0, 0) ; for ( int ni = 0 ; ni < nVertComp[nComp] ; ++ ni) - ptGravityCenter += CompoVert[nComp][ni].ptInt ; + ptGravityCenter += CompoVert[nComp][ni].ptPApp ; ptGravityCenter = ptGravityCenter / nVertComp[nComp] ; Vector3d vtTrasf[12] ; for ( int ni = 0 ; ni < nVertComp[nComp] ; ++ ni) - vtTrasf[ni] = CompoVert[nComp][ni].ptInt - ptGravityCenter ; + vtTrasf[ni] = CompoVert[nComp][ni].ptPApp - ptGravityCenter ; // Preparo le matrici per il sistema typedef Eigen::Matrix dSystemMatrix ; @@ -1685,7 +1747,7 @@ VolZmap::ExtMarchingCubes( int nBlock, TRIA3DEXLIST& lstTria, TriHolder& triHold // medio le normali adiacenti molto vicine (delta angolare inferiore a 22.5 deg) Vector3d vtNorm[12] ; for ( int ni = 0 ; ni < nVertComp[nComp] ; ++ ni) - vtNorm[ni] = CompoVert[nComp][ni].vtNorm ; + vtNorm[ni] = CompoVert[nComp][ni].vtVec ; for ( int ni = 0 ; ni < nVertComp[nComp] ; ++ ni) { int nj = ( ni + 1) % nVertComp[nComp] ; if ( vtNorm[ni] * vtNorm[nj] > 0.92) { @@ -1699,7 +1761,7 @@ VolZmap::ExtMarchingCubes( int nBlock, TRIA3DEXLIST& lstTria, TriHolder& triHold } } - // Caso generale + // Definisco la matrice del sistema for ( int ni = 0 ; ni < nVertComp[nComp] ; ++ ni) { dMatrixN( ni, 0) = vtNorm[ni].x ; dMatrixN( ni, 1) = vtNorm[ni].y ; @@ -1721,474 +1783,40 @@ VolZmap::ExtMarchingCubes( int nBlock, TRIA3DEXLIST& lstTria, TriHolder& triHold // risolvo il sistema con SVD, quindi ai minimi quadrati dSystemVector dUnknownVector( 3, 1) ; dUnknownVector = svd.solve( dKnownVector) ; - // Vettore Baricentro-Feature Vector3d vtFeature( dUnknownVector( 0), dUnknownVector( 1), dUnknownVector( 2)) ; - // Esprimo la soluzione nel sistema di riferimento z-map Point3d ptSol = ptGravityCenter + vtFeature ; - - // Limito la feature entro i secondi voxel vicini, - // nel caso essa esca dal reticolo la limito entro una - // distanza dal baricentro pari alla diagonale del voxel. - bool bOutside = false ; - int nFtVxI, nFtVxJ, nFtVxK ; - if ( GetPointVoxel( ptSol, nFtVxI, nFtVxJ, nFtVxK)) { - if ( abs( nFtVxI - i) > 2 || - abs( nFtVxJ - j) > 2 || - abs( nFtVxK - k) > 2) - bOutside = true ; - } - else { - double dDistFeature = vtFeature.Len() ; - const double MAX_DIST = sqrt( 3) * N_DEXVOXRATIO * m_dStep ; - bOutside = ( dDistFeature > MAX_DIST) ; - } - - // Costruisco triangoli del fan - TRIA3DEXVECTOR triContainer ; - for ( int ni = 0 ; ni < nVertComp[nComp] ; ++ ni) { - int nj = ( ni + 1 < nVertComp[nComp]) ? ni + 1 : 0 ; - // Il triangolo è pronto - Triangle3dEx CurrentTriangle ; - // Setto i punti - CurrentTriangle.Set( ptSol, CompoVert[nComp][nj].ptInt, CompoVert[nComp][ni].ptInt) ; - // Setto il numero di utensile ai vertici di base del fan - CurrentTriangle.SetAttrib( 1, CompoVert[nComp][nj].nToolFlag) ; - CurrentTriangle.SetAttrib( 2, CompoVert[nComp][ni].nToolFlag) ; - // Setto il numero di utensile al triangolo nel complesso - if ( CurrentTriangle.GetAttrib( 1) < 0 || - CurrentTriangle.GetAttrib( 2) < 0) - CurrentTriangle.SetGrade( - 1) ; - else if ( CurrentTriangle.GetAttrib( 1) > 0 || - CurrentTriangle.GetAttrib( 2) > 0) - CurrentTriangle.SetGrade( 1) ; - else - CurrentTriangle.SetGrade( 0) ; - // Setto le normali a ogni vertice - CurrentTriangle.SetVertexNorm( 1, CompoVert[nComp][nj].vtNorm) ; - CurrentTriangle.SetVertexNorm( 2, CompoVert[nComp][ni].vtNorm) ; - if ( CurrentTriangle.GetVertexNorm( 1) * CurrentTriangle.GetVertexNorm( 2) > 0.5) - CurrentTriangle.SetVertexNorm( 0, 0.5 * ( CurrentTriangle.GetVertexNorm( 1) + - CurrentTriangle.GetVertexNorm( 2))) ; - // Valido il triangolo - CurrentTriangle.Validate( true) ; - // Aggiungo triangolo al vettore temporaneo - triContainer.emplace_back( CurrentTriangle) ; - } - - // Controllo delle inversioni dei triangoli - bool bDangerInversion = false ; - - // Caso ventaglio con tre vertici di base - if ( nVertComp[nComp] == 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 ; - for ( int ni = 0 ; ni < nVertComp[nComp] ; ++ ni) { - int nj = ( ni + 1 < nVertComp[nComp]) ? ni + 1 : 0 ; - double dDI = triContainer[ni].GetN() * CompoVert[nComp][ni].vtNorm ; - double dDJ = triContainer[ni].GetN() * CompoVert[nComp][nj].vtNorm ; - if ( dDI < - EPS_SMALL || dDJ < - EPS_SMALL) { - bInversione = true ; - break ; - } - } - // 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[nComp] - 1 ; ++ ni) { - for ( int nj = ni + 1 ; nj < nVertComp[nComp] ; ++ nj) { - if ( CompoVert[nComp][ni].vtNorm * CompoVert[nComp][nj].vtNorm < - EPS_SMALL) - nNegDot ++ ; - } - } - if ( nNegDot == nVertComp[nComp] - 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[nComp] ; ++ ni) { - int nj = ( ni == 0 ? nVertComp[nComp] - 1 : ni - 1) ; - double dDLast = triContainer[ni].GetN() * CompoVert[nComp][ni].vtNorm ; - double dDPrev = triContainer[nj].GetN() * CompoVert[nComp][ni].vtNorm ; - if ( ( dDLast < EPS_SMALL && dDPrev < EPS_SMALL) || - ( dDLast < - 0.75 || dDPrev < - 0.75)) { - bInversione2 = true ; - break ; - } - } - // Se tale normale esiste - if ( bInversione2) { - // 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, ptSol)) { - Vector3d vtNullSpace( dMatrixV( 0, 2), dMatrixV( 1, 2), dMatrixV( 2, 2)) ; - double dParInt1, dParInt2 ; - Point3d ptVoxMin( ( i * N_DEXVOXRATIO + 0.5) * m_dStep, - ( j * N_DEXVOXRATIO + 0.5) * m_dStep, ( k * N_DEXVOXRATIO + 0.5) * m_dStep) ; - Point3d ptVoxMax( ( ( i + 1) * N_DEXVOXRATIO + 0.5) * m_dStep, - ( ( j + 1) * N_DEXVOXRATIO + 0.5) * m_dStep, ( ( k + 1) * N_DEXVOXRATIO + 0.5) * m_dStep) ; - // Caso in cui può essere riportata dentro: se il voxel in cui cade la feature è pieno - // la riporto muovendola lungo la sua linea - if ( IntersLineBox( ptSol, vtNullSpace, ptVoxMin, ptVoxMax, dParInt1, dParInt2)) { - triContainer.resize( 0) ; - double dPar = abs( dParInt1) < abs( dParInt2) ? dParInt1 + ( dParInt2 - dParInt1) / 100 : - dParInt2 + ( dParInt1 - dParInt2) / 100 ; - ptSol += dPar * vtNullSpace ; - // Costruisco triangoli di prova - for ( int ni = 0 ; ni < nVertComp[nComp] ; ++ ni) { - int nj = ( ni + 1 < nVertComp[nComp]) ? ni + 1 : 0 ; - // Il triangolo è pronto - Triangle3dEx CurrentTriangle ; - CurrentTriangle.Set( ptSol, CompoVert[nComp][nj].ptInt, CompoVert[nComp][ni].ptInt) ; - CurrentTriangle.SetAttrib( 1, CompoVert[nComp][nj].nToolFlag) ; - CurrentTriangle.SetAttrib( 2, CompoVert[nComp][ni].nToolFlag) ; - CurrentTriangle.SetVertexNorm( 1, CompoVert[nComp][nj].vtNorm) ; - CurrentTriangle.SetVertexNorm( 2, CompoVert[nComp][ni].vtNorm) ; - if ( CurrentTriangle.GetVertexNorm( 1) * CurrentTriangle.GetVertexNorm( 2) > 0.5) - CurrentTriangle.SetVertexNorm( 0, 0.5 * ( CurrentTriangle.GetVertexNorm( 1) + - CurrentTriangle.GetVertexNorm( 2))) ; - // Valido il triangolo - CurrentTriangle.Validate( true) ; - // Assegno i colori - if ( CurrentTriangle.GetAttrib( 1) < 0 || - CurrentTriangle.GetAttrib( 2) < 0) - CurrentTriangle.SetGrade( - 1) ; - else if ( CurrentTriangle.GetAttrib( 1) > 0 || - CurrentTriangle.GetAttrib( 2) > 0) - CurrentTriangle.SetGrade( 1) ; - else - CurrentTriangle.SetGrade( 0) ; - // Aggiungo triangolo al vettore temporaneo - triContainer.emplace_back( CurrentTriangle) ; - } - } - // Caso in cui non può essere riportata dentro: - // si passa alla routine standard se il voxel - // in cui cade è pieno. - else { - int nAdjVoxI, nAdjVoxJ, nAdjVoxK ; - if ( GetPointVoxel( ptSol, nAdjVoxI, nAdjVoxJ, nAdjVoxK)) { - // Classificazione del voxel adiacente - int nAdjIndex = CalcIndex( nAdjVoxI, nAdjVoxJ, nAdjVoxK) ; - // Se il voxel è pieno - if ( EdgeTable[nAdjIndex] != 0) - bDangerInversion = true ; - } - } - } - } - } - } - } - // Ventaglio con base a quattro vertici - else if ( nVertComp[nComp] == 4) { - // Controllo preliminare sulle normali - bool bInversione = false ; - for ( int ni = 0 ; ni < nVertComp[nComp] ; ++ ni) { - int nj = ( ni + 1 < nVertComp[nComp]) ? ni + 1 : 0 ; - double dDI = triContainer[ni].GetN() * CompoVert[nComp][ni].vtNorm ; - double dDJ = triContainer[ni].GetN() * CompoVert[nComp][nj].vtNorm ; - if ( dDI < - EPS_SMALL || dDJ < - EPS_SMALL) { - bInversione = true ; - break ; - } - } - // Si necessita di un ulteriore controllo - if ( bInversione) { - PolygonPlane Polygon ; - Point3d ptP0 = CompoVert[nComp][0].ptInt ; - Point3d ptP1 = CompoVert[nComp][1].ptInt ; - Point3d ptP2 = CompoVert[nComp][2].ptInt ; - Point3d ptP3 = CompoVert[nComp][3].ptInt ; - Polygon.AddPoint( ptP0) ; - Polygon.AddPoint( ptP1) ; - Polygon.AddPoint( ptP2) ; - Polygon.AddPoint( ptP3) ; - Point3d ptCentr ; - Vector3d vtPolyN ; - bool bCentFound = Polygon.GetCentroid( ptCentr) ; - bool bNormFound = Polygon.GetNormal( vtPolyN) ; - if ( bCentFound && bNormFound) { - Point3d ptPoly0 = ptP0 + ( ptP0 - ptCentr) - ( ptP0 - ptCentr) * vtPolyN * vtPolyN ; - Point3d ptPoly1 = ptP1 + ( ptP1 - ptCentr) - ( ptP1 - ptCentr) * vtPolyN * vtPolyN ; - Point3d ptPoly2 = ptP2 + ( ptP2 - ptCentr) - ( ptP2 - ptCentr) * vtPolyN * vtPolyN ; - Point3d ptPoly3 = ptP3 + ( ptP3 - ptCentr) - ( ptP3 - ptCentr) * vtPolyN * vtPolyN ; - Point3d ptSolP = ptSol + ( ptSol - ptCentr) - ( ptSol - ptCentr) * vtPolyN * vtPolyN ; - Vector3d vtSeg0 = ptPoly1 - ptPoly0 ; - Vector3d vtSeg1 = ptPoly2 - ptPoly1 ; - Vector3d vtSeg2 = ptPoly3 - ptPoly2 ; - Vector3d vtSeg3 = ptPoly0 - ptPoly3 ; - Vector3d vtSol0 = ptSolP - ptPoly0 ; - Vector3d vtSol1 = ptSolP - ptPoly1 ; - Vector3d vtSol2 = ptSolP - ptPoly2 ; - Vector3d vtSol3 = ptSolP - ptPoly3 ; - vtSeg0.Normalize() ; - vtSeg1.Normalize() ; - vtSeg2.Normalize() ; - vtSeg3.Normalize() ; - /*vtSol0.Normalize() ; - vtSol1.Normalize() ; - vtSol2.Normalize() ; - vtSol3.Normalize() ;*/ - double dIn[4] ; - dIn[0] = ( vtSeg0 ^ vtSol0) * vtPolyN ; - dIn[1] = ( vtSeg1 ^ vtSol1) * vtPolyN ; - dIn[2] = ( vtSeg2 ^ vtSol2) * vtPolyN ; - dIn[3] = ( vtSeg3 ^ vtSol3) * vtPolyN ; - double dThr = 0.1 * N_DEXVOXRATIO * m_dStep ; - bool bSpecial = false ; - for ( int nE = 0 ; nE < 4 ; ++ nE) { - if ( dIn[nE] < dThr) { - int nF = nE + 1 < 4 ? nE + 1 : 0 ; - double dDE = triContainer[nE].GetN() * CompoVert[nComp][nE].vtNorm ; - double dDF = triContainer[nE].GetN() * CompoVert[nComp][nF].vtNorm ; - double dLim = - 0.5/*EPS_SMALL*/ ; - if ( dDE < dLim || dDF < dLim) { - bSpecial = true ; - break ; - } - } - } - if ( bSpecial) { - // Se la feature non cade nel suo voxel - if ( ! IsPointInsideVoxelApprox( i, j, k, ptSol)) { - // Vettore di direzione feature - Vector3d vtNullSpace( dMatrixV( 0, 2), dMatrixV( 1, 2), dMatrixV( 2, 2)) ; - double dParInt1, dParInt2 ; - Point3d ptVoxMin( ( i * N_DEXVOXRATIO + 0.5) * m_dStep, - ( j * N_DEXVOXRATIO + 0.5) * m_dStep, ( k * N_DEXVOXRATIO + 0.5) * m_dStep) ; - Point3d ptVoxMax( ( ( i + 1) * N_DEXVOXRATIO + 0.5) * m_dStep, - ( ( j + 1) * N_DEXVOXRATIO + 0.5) * m_dStep, ( ( k + 1) * N_DEXVOXRATIO + 0.5) * m_dStep) ; - // Se è possibile riportarla dentro e il voxel in cui cade è pieno, la riporto nel suo voxel - // lungo la sua linea - if ( IntersLineBox( ptSol, vtNullSpace, ptVoxMin, ptVoxMax, dParInt1, dParInt2)) { - triContainer.resize( 0) ; - double dPar = abs( dParInt1) < abs( dParInt2) ? dParInt1 + ( dParInt2 - dParInt1) / 100 : - dParInt2 + ( dParInt1 - dParInt2) / 100 ; - ptSol += dPar * vtNullSpace ; - // Costruisco triangoli di prova - for ( int ni = 0 ; ni < nVertComp[nComp] ; ++ ni) { - int nj = ( ni + 1 < nVertComp[nComp]) ? ni + 1 : 0 ; - // Il triangolo è pronto - Triangle3dEx CurrentTriangle ; - CurrentTriangle.Set( ptSol, CompoVert[nComp][nj].ptInt, CompoVert[nComp][ni].ptInt) ; - CurrentTriangle.SetAttrib( 1, CompoVert[nComp][nj].nToolFlag) ; - CurrentTriangle.SetAttrib( 2, CompoVert[nComp][ni].nToolFlag) ; - CurrentTriangle.SetVertexNorm( 1, CompoVert[nComp][nj].vtNorm) ; - CurrentTriangle.SetVertexNorm( 2, CompoVert[nComp][ni].vtNorm) ; - if ( CurrentTriangle.GetVertexNorm( 1) * CurrentTriangle.GetVertexNorm( 2) > 0.5) - CurrentTriangle.SetVertexNorm( 0, 0.5 * ( CurrentTriangle.GetVertexNorm( 1) + - CurrentTriangle.GetVertexNorm( 2))) ; - CurrentTriangle.Validate( true) ; - // Assegno i colori - if ( CurrentTriangle.GetAttrib( 1) < 0 || - CurrentTriangle.GetAttrib( 2) < 0) - CurrentTriangle.SetGrade( - 1) ; - else if ( CurrentTriangle.GetAttrib( 1) > 0 || - CurrentTriangle.GetAttrib( 2) > 0) - CurrentTriangle.SetGrade( 1) ; - else - CurrentTriangle.SetGrade( 0) ; - // Aggiungo triangolo al vettore temporaneo - triContainer.emplace_back( CurrentTriangle) ; - } - } - - // Se non è possibile riportarla dentro e il voxel in - // cui cade è pieno passo alla routine standard - else { - int nAdjVoxI, nAdjVoxJ, nAdjVoxK ; - if ( GetPointVoxel( ptSol, nAdjVoxI, nAdjVoxJ, nAdjVoxK)) { - // Classificazione del voxel adiacente - int nAdjIndex = CalcIndex( nAdjVoxI, nAdjVoxJ, nAdjVoxK) ; - // Se il voxel è pieno - if ( EdgeTable[nAdjIndex] != 0) - bDangerInversion = true ; - } - } - } - } - } - } - } - - // Verifico se il fan di triangoli è piano - bool bPlane = true ; - int nContSize = int( triContainer.size()) ; - if ( nContSize > 0) { - Vector3d vtI = triContainer[0].GetN() ; - for ( int nj = 1 ; nj < nContSize ; ++ nj) { - Vector3d vtJ = triContainer[nj].GetN() ; - if ( ! AreSameVectorApprox( vtI, vtJ)) { - bPlane = false ; - break ; - } - } - } - - // Se feature nei limiti e triangoli non in un piano, confermo ExtMC - if ( ! bOutside && ! bPlane && ! bDangerInversion) { - - // Riconoscimento se voxel di frontiera - int nVoxIndexes[3] = { i, j, k} ; - bool bBoundary = IsAVoxelOnBoundary( nLimits, nVoxIndexes, true) ; - - TRIA3DEXVECTOR vInnerTriaTemp, vBorderTriaTemp ; - - for ( int ni = 0 ; ni < nContSize ; ++ ni) { - - // Se l'area è non nulla aggiungo il triangolo a uno dei due vettori. - if ( triContainer[ni].GetArea() > SQ_EPS_SMALL) { - int nVoxIJK[3] = { i, j, k} ; - Point3d ptTemp = ptSol ; - Triangle3dEx trTemp = triContainer[ni] ; - bool bTriangleOnBorder = IsATriangleOnBorder( trTemp, ptTemp, nLimits, nVoxIJK) ; - if ( ! bBoundary || ! bTriangleOnBorder) - 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) ; - - // 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, - int nCurrent = int( triHold.size()) - 1 ; - triHold[nCurrent].i = i ; - triHold[nCurrent].j = j ; - triHold[nCurrent].k = k ; - } - - // Indice che identifica la posizione del voxel nel vector. - int nCurrent = int( triHold.size()) - 1 ; - - // Aggiungo vertice della componente - // connessa alla lista dei vertici. - triHold[nCurrent].ptCompoVert.emplace_back( ptSol) ; - - int nOldFeatureNum = int( triHold[nCurrent].vCompoTria.size()) ; - int nNewFeatureNum = nOldFeatureNum + 1 ; - - // Aggiungo una componente per il vettore - // dei triangoli della componente connessa. - triHold[nCurrent].vCompoTria.resize( nNewFeatureNum) ; - triHold[nCurrent].vbFlipped.resize( nNewFeatureNum) ; - for ( int ni = 0 ; ni < int( vInnerTriaTemp.size()) ; ++ ni) { - // Riporto le coordinate nel sistema in cui è immerso lo Zmap - vInnerTriaTemp[ni].ToGlob( m_MapFrame) ; - triHold[nCurrent].vCompoTria[nOldFeatureNum].emplace_back( vInnerTriaTemp[ni]) ; - } - triHold[nCurrent].vbFlipped[nOldFeatureNum].resize( vInnerTriaTemp.size()) ; - for ( int ni = 0 ; ni < int( vInnerTriaTemp.size()) ; ++ ni) { - triHold[nCurrent].vbFlipped[nOldFeatureNum][ni] = false ; - } - } - - // Triangoli di frontiera - if ( vBorderTriaTemp.size() > 0) { - - // Aggiorno il numero di feature. - ++ nBorderFeatureInVoxel ; - - // Se siamo in presenza della prima feature del - // voxel, ridimensiono il vettore che contiene - // la struttura dati del voxel. - if ( nBorderFeatureInVoxel == 1) { - // Aggiungo un voxel al blocco per i triangoli - m_InterBlockTria[nBlock].resize( m_InterBlockTria[nBlock].size() + 1) ; - // Questi dati dipendono solo dal voxel - int nCurrent = int( m_InterBlockTria[nBlock].size()) - 1 ; - 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 ; - - // Aggiungo vertice della componente connessa alla lista dei vertici - 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) { - // Riporto le coordinate nel sistema in cui è immerso lo Zmap - vBorderTriaTemp[ni].ToGlob( m_MapFrame) ; - m_InterBlockTria[nBlock][nCurrent].vCompoTria[nOldFeatureNum].emplace_back( vBorderTriaTemp[ni]) ; - } - } - } - - // ExtMC non confermato, si passa a MC - else { - vector vTria ; - // Costruzione dei triangoli - for ( int TriIndex = 0; TriIndex < ( nVertComp[nComp] - 2) * 3 ; TriIndex += 3) { - // Il triangolo è pronto - Triangle3dEx CurrentTriangle ; - CurrentTriangle.Set( CompoTriVert[nComp][TriIndex].ptInt, - CompoTriVert[nComp][TriIndex+1].ptInt, - CompoTriVert[nComp][TriIndex+2].ptInt) ; - // Setto il numero di utensile (conta solo positivo, nullo o negativo) - int nTool0 = Clamp( CompoTriVert[nComp][TriIndex].nToolFlag, -1, 1) ; - int nTool1 = Clamp( CompoTriVert[nComp][TriIndex+1].nToolFlag, -1, 1) ; - int nTool2 = Clamp( CompoTriVert[nComp][TriIndex+2].nToolFlag, -1, 1) ; - // Setto il numero dei colori - if ( nTool0 == nTool1 || nTool0 == nTool2) - CurrentTriangle.SetGrade( nTool0) ; - else if ( nTool1 == nTool2) - CurrentTriangle.SetGrade( nTool1) ; - // Valido il triangolo e setto le normali del campo vettoriale ai corrispondenti vertici - if ( CurrentTriangle.Validate( true)) { - for ( int nV = 0 ; nV < 3 ; ++ nV) { - const Vector3d& vtVertNorm = CompoTriVert[nComp][TriIndex+nV].vtNorm ; - if ( CurrentTriangle.GetN() * vtVertNorm > 0.6) - CurrentTriangle.SetVertexNorm( nV, vtVertNorm) ; - else - CurrentTriangle.SetVertexNorm( nV, CurrentTriangle.GetN()) ; - } - } - // Riporto le coordinate nel sistema in cui è immerso lo Zmap - CurrentTriangle.ToGlob( m_MapFrame) ; - // Aggiungo alla lista - vTria.emplace_back( CurrentTriangle) ; - } - // Controllo i colori di configurazioni 2 e 8 - if ( ( nAllConfig[nIndex] == 2 || nAllConfig[nIndex] == 8) && - vTria[0].GetN() * vTria[1].GetN() > 0.8) { - if ( vTria[0].GetGrade() < 0 || vTria[1].GetGrade() < 0) { - vTria[0].SetGrade( -1) ; - vTria[1].SetGrade( -1) ; - } - else if ( vTria[0].GetGrade() > 0 || vTria[1].GetGrade() > 0) { - vTria[0].SetGrade( 1) ; - vTria[1].SetGrade( 1) ; - } - } - for ( int nT = 0 ; nT < int( vTria.size()) ; ++ nT) - lstTria.emplace_back( vTria[nT]) ; + + size_t tOldCompo = VoxConf.nNumComp ; + ++ VoxConf.nNumComp ; + VoxConf.Compo[tOldCompo].nVertNum = nVertComp[nComp] ; + for ( int nV = 0 ; nV < nVertComp[nComp] ; ++ nV) { + VoxConf.Compo[tOldCompo].CompVecField[nV] = CompoVert[nComp][nV] ; } + // Esperimento + //for ( int nVtI = 0 ; nVtI < VoxConf[tOldSize].Compo[tOldCompo].nVertNum ; ++ nVtI) { + // int nVtJ = ( nVtI + 1) % nVertComp[nComp] ; + // Point3d ptPI = VoxConf[tOldSize].Compo[tOldCompo].CompVecField[nVtI].ptPApp ; + // Point3d ptPJ = VoxConf[tOldSize].Compo[tOldCompo].CompVecField[nVtJ].ptPApp ; + // Vector3d vtNI = VoxConf[tOldSize].Compo[tOldCompo].CompVecField[nVtI].vtVec ; + // Vector3d vtNJ = VoxConf[tOldSize].Compo[tOldCompo].CompVecField[nVtJ].vtVec ; + // Vector3d vtN = ( ptPJ - ptSol) ^ ( ptPI - ptPJ) ; + // double dFactor = 50 * EPS_SMALL * N_DEXVOXRATIO * m_dStep * vtN.Len() ; + // vtN.Normalize() ; + // if ( AreOppositeVectorApprox( vtN, vtNI) && AreOppositeVectorApprox( vtN, vtNJ)) { + // // Vers 1 + // //ptSol += dFactor * vtN ; + // // Vers 2 + // ptSol -= dFactor * 0.5 * ( vtNI + vtNJ) ; + // break ; + // } + //} + // Fine esperimento + VoxConf.Compo[tOldCompo].ptVert = ptSol ; + VoxConf.Compo[tOldCompo].bInside = IsPointInsideVoxelApprox( i, j, k, ptSol, 0) ; + VoxConf.Compo[tOldCompo].bCorner = ( nFeatureType == CORNER) ; } - // Standard MC else { vector vTria ; @@ -2196,13 +1824,13 @@ VolZmap::ExtMarchingCubes( int nBlock, TRIA3DEXLIST& lstTria, TriHolder& triHold for ( int TriIndex = 0; TriIndex < ( nVertComp[nComp] - 2) * 3 ; TriIndex += 3) { // Il triangolo è pronto Triangle3dEx CurrentTriangle ; - CurrentTriangle.Set( CompoTriVert[nComp][TriIndex].ptInt, - CompoTriVert[nComp][TriIndex+1].ptInt, - CompoTriVert[nComp][TriIndex+2].ptInt) ; + CurrentTriangle.Set( CompoTriVert[nComp][TriIndex].ptPApp, + CompoTriVert[nComp][TriIndex+1].ptPApp, + CompoTriVert[nComp][TriIndex+2].ptPApp) ; // Setto il numero di utensile (conta solo positivo, nullo o negativo) - int nTool0 = Clamp( CompoTriVert[nComp][TriIndex].nToolFlag, -1, 1) ; - int nTool1 = Clamp( CompoTriVert[nComp][TriIndex+1].nToolFlag, -1, 1) ; - int nTool2 = Clamp( CompoTriVert[nComp][TriIndex+2].nToolFlag, -1, 1) ; + int nTool0 = Clamp( CompoTriVert[nComp][TriIndex].nPropIndex, -1, 1) ; + int nTool1 = Clamp( CompoTriVert[nComp][TriIndex+1].nPropIndex, -1, 1) ; + int nTool2 = Clamp( CompoTriVert[nComp][TriIndex+2].nPropIndex, -1, 1) ; if ( nTool0 == nTool1 || nTool0 == nTool2) CurrentTriangle.SetGrade( nTool0) ; else if ( nTool1 == nTool2) @@ -2210,14 +1838,13 @@ VolZmap::ExtMarchingCubes( int nBlock, TRIA3DEXLIST& lstTria, TriHolder& triHold // Valido il triangolo e setto le normali del campo vettoriale ai corrispondenti vertici if ( CurrentTriangle.Validate( true)) { for ( int nV = 0 ; nV < 3 ; ++ nV) { - const Vector3d& vtVertNorm = CompoTriVert[nComp][TriIndex+nV].vtNorm ; - if ( CurrentTriangle.GetN() * vtVertNorm > 0.6) + const Vector3d& vtVertNorm = CompoTriVert[nComp][TriIndex+nV].vtVec ; + if ( CurrentTriangle.GetN() * vtVertNorm > 0.5) CurrentTriangle.SetVertexNorm( nV, vtVertNorm) ; else CurrentTriangle.SetVertexNorm( nV, CurrentTriangle.GetN()) ; } - } - // si può inserire controllo sui colori per i casi 2 e 8 come per ext non confermato + } // Riporto le coordinate nel sistema in cui è immerso lo Zmap CurrentTriangle.ToGlob( m_MapFrame) ; // Aggiungo alla lista @@ -2237,12 +1864,26 @@ VolZmap::ExtMarchingCubes( int nBlock, TRIA3DEXLIST& lstTria, TriHolder& triHold } for ( int nT = 0 ; nT < int( vTria.size()) ; ++ nT) lstTria.emplace_back( vTria[nT]) ; - } + } + } + // Se nel voxel abbiamo trovato feature + // aggiorniamo i contenitori. + if ( VoxConf.nNumComp > 0) { + VoxConf.i = i ; + VoxConf.j = j ; + VoxConf.k = k ; + int nIJK[3] = { i, j, k} ; + int nKey ; + GetVoxNFromIJK( i, j, k, nKey) ; + if ( IsAVoxelOnBoundary( nLimits, nIJK, true)) + m_InterBlockVox[nBlock].emplace( nKey, VoxConf) ; + else + vVox.emplace( nKey, VoxConf) ; } } } } - + // Processo i Voxel con possibile superficie piana ProcessVoxContXY( VoxContXYInf, false, lstTria) ; ProcessVoxContXY( VoxContXYSup, true, lstTria) ; @@ -2254,12 +1895,461 @@ VolZmap::ExtMarchingCubes( int nBlock, TRIA3DEXLIST& lstTria, TriHolder& triHold return true ; } +//---------------------------------------------------------------------------- +bool +VolZmap::RegulateFeaturesChain( std::vector& vVecVox) const +{ + // Ciclo sui blocchi + for ( int nBlock = 0 ; nBlock < int( m_nNumBlock) ; ++ nBlock) { + // Se il blocco è da aggiornare + if ( m_BlockToUpdate[nBlock]) { + // Ciclo sui voxel interni + for ( auto itVox = vVecVox[nBlock].begin() ; itVox != vVecVox[nBlock].end() ; ++ itVox) { + int nVox ; + Voxel& CurVox = itVox->second ; + GetVoxNFromIJK( CurVox.i, CurVox.j, CurVox.k, nVox) ; + // Ciclo sulle componenti + for ( int nComp = 0 ; nComp < CurVox.nNumComp ; ++ nComp) { + if ( nBlock == 1 && nVox == 41024 && nComp == 0) + int a = 1 ; + // Vertice fuori dal suo voxel + if ( ! CurVox.Compo[nComp].bInside) { + // Caso corner + if ( CurVox.Compo[nComp].bCorner) { + // Cerco i primi vicini + std::vector vNearFirst, vBordNearFirst ; + FindAdjComp( vVecVox, nBlock, nVox, nComp, vNearFirst, vBordNearFirst) ; + // Ciclo sui primi vicini + int nInnSizeF = int( vNearFirst.size()) ; + int nBorSizeF = int( vBordNearFirst.size()) ; + for ( int nF = 0 ; nF < nInnSizeF + nBorSizeF ; nF += 3) { + // Indice e vettore corrente dei primi vicini + int nFCur = nF < nInnSizeF ? nF : nF - nInnSizeF ; + std::vector vVecNFCur = nF < nInnSizeF ? vNearFirst : vBordNearFirst ; + // Cerco i secondi vicini + std::vector vNearSecond, vBordNearSecond ; + FindAdjComp( vVecVox, vVecNFCur[nFCur], vVecNFCur[nFCur+1], vVecNFCur[nFCur+2], + vNearSecond, vBordNearSecond) ; + // Ciclo sui secondi vicini + int nInnSizeS = int( vNearSecond.size()) ; + int nBorSizeS = int( vBordNearSecond.size()) ; + for ( int nS = 0 ; nS < nInnSizeS + nBorSizeS ; nS += 3) { + // Indice e vettore corrente dei secondi vicini + int nNSCur = nS < nInnSizeS ? nS : nS - nInnSizeS ; + std::vector vVecNSCur = nS < nInnSizeS ? vNearSecond : vBordNearSecond ; + // Escludo dai secondi i primi vicini + bool bFirst = false ; + if ( vVecNSCur[nNSCur] == nBlock && + vVecNSCur[nNSCur+1] == nVox && + vVecNSCur[nNSCur+2] == nComp) + bFirst = true ; + for ( int nOldF = 0 ; nOldF < nInnSizeF + nBorSizeF ; nOldF += 3) { + if ( nOldF < nInnSizeF) { + if ( vVecNSCur[nNSCur] == vNearFirst[nOldF] && + vVecNSCur[nNSCur+1] == vNearFirst[nOldF+1] && + vVecNSCur[nNSCur+2] == vNearFirst[nOldF+2]) + bFirst = true ; + } + else { + if ( vVecNSCur[nNSCur] == vBordNearFirst[nOldF-nInnSizeF] && + vVecNSCur[nNSCur+1] == vBordNearFirst[nOldF-nInnSizeF+1] && + vVecNSCur[nNSCur+2] == vBordNearFirst[nOldF-nInnSizeF+2]) + bFirst = true ; + } + } + // Se trovo un secondo fra i primi salto l'iterazione + if ( bFirst) + continue ; + // Se necessario regolarizzo la catena + Voxel& VoxNearFirst = nF < nInnSizeF ? vVecVox[vVecNFCur[nFCur]].find( vVecNFCur[nFCur+1])->second : + m_InterBlockVox[vVecNFCur[nFCur]].find( vVecNFCur[nFCur+1])->second ; + Voxel& VoxNearSecond = nS < nInnSizeS ? vVecVox[vVecNSCur[nNSCur]].find( vVecNSCur[nNSCur+1])->second : + m_InterBlockVox[vVecNSCur[nNSCur]].find( vVecNSCur[nNSCur+1])->second ; + Point3d ptCurV = CurVox.Compo[nComp].ptVert ; + Point3d ptFirst = VoxNearFirst.Compo[vVecNFCur[nFCur+2]].ptVert ; + Point3d ptSecond = VoxNearSecond.Compo[vVecNSCur[nNSCur+2]].ptVert ; + Vector3d vtF = ptFirst - ptCurV ; + Vector3d vtS = ptSecond - ptCurV ; + vtF.Normalize() ; + vtS.Normalize() ; + // Calcolo i baricentri dei vertici dei fan + Point3d ptCurBar ; + for ( int n = 0 ; n < CurVox.Compo[nComp].nVertNum ; ++ n) + ptCurBar += CurVox.Compo[nComp].CompVecField[n].ptPApp ; + ptCurBar /= CurVox.Compo[nComp].nVertNum ; + Point3d ptFirstBar ; + for ( int n = 0 ; n < VoxNearFirst.Compo[vVecNFCur[nFCur+2]].nVertNum ; ++ n) + ptFirstBar += VoxNearFirst.Compo[vVecNFCur[nFCur+2]].CompVecField[n].ptPApp ; + ptFirstBar /= VoxNearFirst.Compo[vVecNFCur[nFCur+2]].nVertNum ; + Point3d ptSecondBar ; + for ( int n = 0 ; n < VoxNearSecond.Compo[vVecNSCur[nNSCur+2]].nVertNum ; ++ n) + ptSecondBar += VoxNearSecond.Compo[vVecNSCur[nNSCur+2]].CompVecField[n].ptPApp ; + ptSecondBar /= VoxNearSecond.Compo[vVecNSCur[nNSCur+2]].nVertNum ; + Vector3d vtBarCF = ptFirstBar - ptCurBar ; + Vector3d vtBarCS = ptSecondBar - ptCurBar ; + Vector3d vtBarSF = ptSecondBar - ptFirstBar ; + vtBarCF.Normalize() ; + vtBarCS.Normalize() ; + vtBarSF.Normalize() ; + if ( vtBarCF * vtBarCS > 0.5 && vtBarCF * vtBarSF > 0.5 && vtBarCS * vtBarSF > 0.5 && + vtF * vtS < - 0.9 && ! VoxNearFirst.Compo[vVecNFCur[nFCur+2]].bCorner) { + VoxNearFirst.Compo[vVecNFCur[nFCur+2]].ptVert = 0.5 * ( ptCurV + ptSecond) ; + } + } + } + } + // Caso feature + else { + std::vector vNearInn, vNearBord ; + FindAdjComp( vVecVox, nBlock, nVox, nComp, vNearInn, vNearBord) ; + int nSizeInn = int( vNearInn.size()) ; + int nSizeBord = int( vNearBord.size() ); + if ( nSizeInn + nSizeBord == 6) { + Voxel* pVoxSt = new Voxel ; + Voxel* pVoxEn = new Voxel ; + if ( nSizeInn == 6) { + *pVoxSt = vVecVox[vNearInn[0]].find( vNearInn[1])->second ; + *pVoxEn = vVecVox[vNearInn[3]].find( vNearInn[4])->second ; + } + else if ( nSizeBord == 6) { + *pVoxSt = m_InterBlockVox[vNearBord[0]].find( vNearBord[1])->second ; + *pVoxEn = m_InterBlockVox[vNearBord[3]].find( vNearBord[4])->second ; + } + else { + *pVoxSt = vVecVox[vNearInn[0]].find( vNearInn[1])->second ; + *pVoxEn = m_InterBlockVox[vNearBord[0]].find( vNearBord[1])->second ; + } + Voxel& AdjVoxSt = * pVoxSt ; + Voxel& AdjVoxEn = * pVoxEn ; + Point3d ptPCur = CurVox.Compo[nComp].ptVert ; + Point3d ptPSt ; + Point3d ptPEn ; + if ( nSizeInn == 6) { + ptPSt = AdjVoxSt.Compo[vNearInn[2]].ptVert ; + ptPEn = AdjVoxEn.Compo[vNearInn[5]].ptVert ; + } + else if ( nSizeBord == 6) { + ptPSt = AdjVoxSt.Compo[vNearBord[2]].ptVert ; + ptPEn = AdjVoxEn.Compo[vNearBord[5]].ptVert ; + } + else { + ptPSt = AdjVoxSt.Compo[vNearInn[2]].ptVert ; + ptPEn = AdjVoxEn.Compo[vNearBord[2]].ptVert ; + } + Vector3d vtStCurr = ptPCur - ptPSt ; + Vector3d vtStEn = ptPEn - ptPSt ; + Vector3d vtCurrEn = ptPEn - ptPCur ; + vtStCurr.Normalize() ; + vtStEn.Normalize() ; + vtCurrEn.Normalize() ; + Point3d ptMid = 0.5 * ( ptPSt + ptPEn) ; + double dMidU = ( ptMid - ptPSt) * vtStEn ; + double dCurU = ( ptPCur - ptPSt) * vtStEn ; + Point3d ptNew ; + Point3d ptPLine ; + Vector3d vtDLine ; + if ( dMidU < dCurU) { + ptPLine = ptPEn ; + vtDLine = - vtCurrEn ; + ptNew = ptPEn + ( ptMid - ptPEn) * vtCurrEn * vtCurrEn ; + } + else { + ptPLine = ptPSt ; + vtDLine = vtStCurr ; + ptNew = ptPSt + ( ptMid - ptPSt) * vtStCurr * vtStCurr ; + } + Point3d ptCubeInf( ( CurVox.i * N_DEXVOXRATIO + 0.5) * m_dStep, + ( CurVox.j * N_DEXVOXRATIO + 0.5) * m_dStep, + ( CurVox.k * N_DEXVOXRATIO + 0.5) * m_dStep) ; + Point3d ptCubeSup( ( ( CurVox.i + 1) * N_DEXVOXRATIO + 0.5) * m_dStep, + ( ( CurVox.j + 1) * N_DEXVOXRATIO + 0.5) * m_dStep, + ( ( CurVox.k + 1) * N_DEXVOXRATIO + 0.5) * m_dStep) ; + double dU1, dU2 ; + if ( 1 - abs( vtStCurr * vtCurrEn ) < EPS_ZERO && + IntersLineBox( ptPLine, vtDLine, ptCubeInf, ptCubeSup, dU1, dU2)) { + double dU = abs( dU1) < abs( dU2) ? dU1 + ( dU2 - dU1) / 2 : dU2 + ( dU1 - dU2) / 2 ; + ptNew = ptPLine + dU * vtDLine ; + } + bool bNewInside = IsPointInsideVoxelApprox( CurVox.i, CurVox.j, CurVox.k, ptNew, 0) ; + if ( abs( vtStCurr * vtStEn) > 0.95 && abs( vtStCurr * vtCurrEn) > 0.95 && + abs( vtStEn * vtCurrEn) > 0.95 /*&& bNewInside*/) { + CurVox.Compo[nComp].ptVert = ptNew ; + } + } + } + } + } + } + // Ciclo sui voxel di frontiera + for ( auto itVox = m_InterBlockVox[nBlock].begin() ; itVox != m_InterBlockVox[nBlock].end() ; ++ itVox) { + int nVox ; + Voxel& CurVox = itVox->second ; + GetVoxNFromIJK( CurVox.i, CurVox.j, CurVox.k, nVox) ; + // Ciclo sulle componenti + for ( int nComp = 0 ; nComp < CurVox.nNumComp ; ++ nComp) { + // Vertice fuori dal suo voxel + if ( ! CurVox.Compo[nComp].bInside) { + // Caso feature + if ( ! CurVox.Compo[nComp].bCorner) { + std::vector vNearInn, vNearBord ; + FindAdjComp( vVecVox, nBlock, nVox, nComp, vNearInn, vNearBord) ; + int nSizeInn = int( vNearInn.size()) ; + int nSizeBord = int( vNearBord.size() ); + if ( nSizeInn + nSizeBord == 6) { + Voxel* pVoxSt = new Voxel ; + Voxel* pVoxEn = new Voxel ; + if ( nSizeInn == 6) { + *pVoxSt = vVecVox[vNearInn[0]].find( vNearInn[1])->second ; + *pVoxEn = vVecVox[vNearInn[3]].find( vNearInn[4])->second ; + } + else if ( nSizeBord == 6) { + *pVoxSt = m_InterBlockVox[vNearBord[0]].find( vNearBord[1])->second ; + *pVoxEn = m_InterBlockVox[vNearBord[3]].find( vNearBord[4])->second ; + } + else { + *pVoxSt = vVecVox[vNearInn[0]].find( vNearInn[1])->second ; + *pVoxEn = m_InterBlockVox[vNearBord[0]].find( vNearBord[1])->second ; + } + Voxel& AdjVoxSt = * pVoxSt ; + Voxel& AdjVoxEn = * pVoxEn ; + Point3d ptPCur = CurVox.Compo[nComp].ptVert ; + Point3d ptPSt ; + Point3d ptPEn ; + if ( nSizeInn == 6) { + ptPSt = AdjVoxSt.Compo[vNearInn[2]].ptVert ; + ptPEn = AdjVoxEn.Compo[vNearInn[5]].ptVert ; + } + else if ( nSizeBord == 6) { + ptPSt = AdjVoxSt.Compo[vNearBord[2]].ptVert ; + ptPEn = AdjVoxEn.Compo[vNearBord[5]].ptVert ; + } + else { + ptPSt = AdjVoxSt.Compo[vNearInn[2]].ptVert ; + ptPEn = AdjVoxEn.Compo[vNearBord[2]].ptVert ; + } + Vector3d vtStCurr = ptPCur - ptPSt ; + Vector3d vtStEn = ptPEn - ptPSt ; + Vector3d vtCurrEn = ptPEn - ptPCur ; + vtStCurr.Normalize() ; + vtStEn.Normalize() ; + vtCurrEn.Normalize() ; + Point3d ptMid = 0.5 * ( ptPSt + ptPEn) ; + double dMidU = ( ptMid - ptPSt) * vtStEn ; + double dCurU = ( ptPCur - ptPSt) * vtStEn ; + Point3d ptNew ; + Point3d ptPLine ; + Vector3d vtDLine ; + if ( dMidU < dCurU) { + ptPLine = ptPEn ; + vtDLine = - vtCurrEn ; + ptNew = ptPEn + ( ptMid - ptPEn) * vtCurrEn * vtCurrEn ; + } + else { + ptPLine = ptPSt ; + vtDLine = vtStCurr ; + ptNew = ptPSt + ( ptMid - ptPSt) * vtStCurr * vtStCurr ; + } + Point3d ptCubeInf( ( CurVox.i * N_DEXVOXRATIO + 0.5) * m_dStep, + ( CurVox.j * N_DEXVOXRATIO + 0.5) * m_dStep, + ( CurVox.k * N_DEXVOXRATIO + 0.5) * m_dStep) ; + Point3d ptCubeSup( ( ( CurVox.i + 1) * N_DEXVOXRATIO + 0.5) * m_dStep, + ( ( CurVox.j + 1) * N_DEXVOXRATIO + 0.5) * m_dStep, + ( ( CurVox.k + 1) * N_DEXVOXRATIO + 0.5) * m_dStep) ; + double dU1, dU2 ; + if ( 1 - abs( vtStCurr * vtCurrEn ) < EPS_ZERO && + IntersLineBox( ptPLine, vtDLine, ptCubeInf, ptCubeSup, dU1, dU2)) { + double dU = abs( dU1) < abs( dU2) ? dU1 + ( dU2 - dU1) / 2 : dU2 + ( dU1 - dU2) / 2 ; + ptNew = ptPLine + dU * vtDLine ; + } + bool bNewInside = IsPointInsideVoxelApprox( CurVox.i, CurVox.j, CurVox.k, ptNew, 0) ; + if ( abs( vtStCurr * vtStEn) > 0.95 && abs( vtStCurr * vtCurrEn) > 0.95 && + abs( vtStEn * vtCurrEn) > 0.95 /*&& bNewInside*/) { + CurVox.Compo[nComp].ptVert = ptNew ; + } + } + } + } + } + } + } + } + return true ; +} + +//---------------------------------------------------------------------------- +// A partire dalla struttura dati di tipo voxel con sharp features crea i triangoli corrispondenti. +// Accetta come parametri il numero del blocco, l'unordered map di voxel e l'apposito contenitore +// ove salva i triangoli. +bool +VolZmap::CreateSharpFeatureTriangle( int nBlock, const VoxelContainer& vVoxel, TriHolder& triHold) const +{ + // Calcolo i limiti sugli indici dei voxel del blocco + // Vettore indici i,j,k del blocco + int nBlockIJK[3] ; + GetBlockIJKFromN( nBlock, nBlockIJK) ; + // Vettore limiti sugli indici dei voxel del blocco + int nLimits[6] ; + GetBlockLimitsIJK( nBlockIJK, nLimits) ; + + // Ciclo sui voxel interni + for ( auto it = vVoxel.begin() ; it != vVoxel.end() ; ++ it) { + size_t tOldSize = triHold.size() ; + triHold.resize( tOldSize + 1) ; + triHold[tOldSize].i = it->second.i ; + triHold[tOldSize].j = it->second.j ; + triHold[tOldSize].k = it->second.k ; + // Ciclo sulle componenti connesse del voxel + for ( int nComp = 0 ; nComp < it->second.nNumComp ; ++ nComp) { + triHold[tOldSize].ptCompoVert.emplace_back( it->second.Compo[nComp].ptVert) ; + triHold[tOldSize].ptCompoVert.back().ToGlob( m_MapFrame) ; + size_t tOldCompNum = triHold[tOldSize].vCompoTria.size() ; + triHold[tOldSize].vCompoTria.resize( tOldCompNum + 1) ; + triHold[tOldSize].vbFlipped.resize( tOldCompNum + 1) ; + // ciclo sui vertici della componente connessa + int nNumVert = it->second.Compo[nComp].nVertNum ; + for ( int nVert = 0 ; nVert < nNumVert ; ++ nVert) { + int nNextVert = ( nVert + 1 < nNumVert ? nVert + 1 : 0) ; + // Definisco il triangolo + Triangle3dEx CurrTri ; + CurrTri.Set( it->second.Compo[nComp].ptVert, + it->second.Compo[nComp].CompVecField[nNextVert].ptPApp, + it->second.Compo[nComp].CompVecField[nVert].ptPApp) ; + // Setto il numero di utensile ai vertici di base del fan + CurrTri.SetAttrib( 1, it->second.Compo[nComp].CompVecField[nNextVert].nPropIndex) ; + CurrTri.SetAttrib( 2, it->second.Compo[nComp].CompVecField[nVert].nPropIndex) ; + // Setto il numero di utensile al triangolo nel complesso + if ( CurrTri.GetAttrib( 1) < 0 || + CurrTri.GetAttrib( 2) < 0) + CurrTri.SetGrade( - 1) ; + else if ( CurrTri.GetAttrib( 1) > 0 || + CurrTri.GetAttrib( 2) > 0) + CurrTri.SetGrade( 1) ; + else + CurrTri.SetGrade( 0) ; + // Setto le normali a ogni vertice + CurrTri.SetVertexNorm( 1, it->second.Compo[nComp].CompVecField[nNextVert].vtVec) ; + CurrTri.SetVertexNorm( 2, it->second.Compo[nComp].CompVecField[nVert].vtVec) ; + if ( CurrTri.GetVertexNorm( 1) * CurrTri.GetVertexNorm( 2) > 0.5) + CurrTri.SetVertexNorm( 0, 0.5 * ( CurrTri.GetVertexNorm( 1) + + CurrTri.GetVertexNorm( 2))) ; + // Valido il triangolo + CurrTri.Validate( true) ; + triHold[tOldSize].vCompoTria[tOldCompNum].emplace_back( CurrTri) ; + triHold[tOldSize].vCompoTria[tOldCompNum].back().ToGlob( m_MapFrame) ; + size_t tTri = triHold[tOldSize].vbFlipped[tOldCompNum].size() ; + triHold[tOldSize].vbFlipped[tOldCompNum].resize( tTri + 1) ; + triHold[tOldSize].vbFlipped[tOldCompNum][tTri] = false ; + } + } + } + + // Pulisco il contenitore dei voxel di frontiera + m_InterBlockTria[nBlock].clear() ; + + // Ciclo sui voxel di frontiera + for ( auto itVox = m_InterBlockVox[nBlock].begin() ; itVox != m_InterBlockVox[nBlock].end() ; ++ itVox) { + // Indici del voxel + int nVoxIJK[3] = { itVox->second.i, + itVox->second.j, + itVox->second.k} ; + // Ridimensiono il contenitore dei triangoli interni + int nOldSizeInn = int( triHold.size()) ; + triHold.resize( nOldSizeInn + 1) ; + triHold[nOldSizeInn].i = nVoxIJK[0] ; + triHold[nOldSizeInn].j = nVoxIJK[1] ; + triHold[nOldSizeInn].k = nVoxIJK[2] ; + // Ridimensiono il contenitore dei triangoli di frontiera + int nOldSizeBor = int( m_InterBlockTria[nBlock].size()) ; + m_InterBlockTria[nBlock].resize( nOldSizeBor + 1) ; + m_InterBlockTria[nBlock][nOldSizeBor].i = nVoxIJK[0] ; + m_InterBlockTria[nBlock][nOldSizeBor].j = nVoxIJK[1] ; + m_InterBlockTria[nBlock][nOldSizeBor].k = nVoxIJK[2] ; + // Ciclo sulle componenti connesse del voxel + for ( int nComp = 0 ; nComp < itVox->second.nNumComp ; ++ nComp) { + bool bNewCompInn = true ; + bool bNewCompBor = true ; + // ciclo sui vertici della componente connessa + int nNumVert = itVox->second.Compo[nComp].nVertNum ; + for ( int nVert = 0 ; nVert < nNumVert ; ++ nVert) { + int nNextVert = ( nVert + 1 < nNumVert ? nVert + 1 : 0) ; + // Definisco il triangolo + Triangle3dEx CurrTri ; + CurrTri.Set( itVox->second.Compo[nComp].ptVert, + itVox->second.Compo[nComp].CompVecField[nNextVert].ptPApp, + itVox->second.Compo[nComp].CompVecField[nVert].ptPApp) ; + // Setto il numero di utensile ai vertici di base del fan + CurrTri.SetAttrib( 1, itVox->second.Compo[nComp].CompVecField[nNextVert].nPropIndex) ; + CurrTri.SetAttrib( 2, itVox->second.Compo[nComp].CompVecField[nVert].nPropIndex) ; + // Setto il numero di utensile al triangolo nel complesso + if ( CurrTri.GetAttrib( 1) < 0 || + CurrTri.GetAttrib( 2) < 0) + CurrTri.SetGrade( - 1) ; + else if ( CurrTri.GetAttrib( 1) > 0 || + CurrTri.GetAttrib( 2) > 0) + CurrTri.SetGrade( 1) ; + else + CurrTri.SetGrade( 0) ; + // Setto le normali a ogni vertice + CurrTri.SetVertexNorm( 1, itVox->second.Compo[nComp].CompVecField[nNextVert].vtVec) ; + CurrTri.SetVertexNorm( 2, itVox->second.Compo[nComp].CompVecField[nVert].vtVec) ; + if ( CurrTri.GetVertexNorm( 1) * CurrTri.GetVertexNorm( 2) > 0.5) + CurrTri.SetVertexNorm( 0, 0.5 * ( CurrTri.GetVertexNorm( 1) + + CurrTri.GetVertexNorm( 2))) ; + // Valido il triangolo + CurrTri.Validate( true) ; + // Smisto i triangoli fra di frontiera e interni + bool bTriOnBorder = IsATriangleOnBorder( CurrTri, itVox->second.Compo[nComp].ptVert, + nLimits, nVoxIJK) ; + // Triangolo di frontiera + if ( bTriOnBorder) { + if ( bNewCompBor) { + Point3d ptVert = itVox->second.Compo[nComp].ptVert ; + ptVert.ToGlob( m_MapFrame) ; + m_InterBlockTria[nBlock][nOldSizeBor].ptCompoVert.emplace_back( ptVert) ; + size_t tOldComp = m_InterBlockTria[nBlock][nOldSizeBor].vCompoTria.size() ; + m_InterBlockTria[nBlock][nOldSizeBor].vCompoTria.resize( tOldComp + 1) ; + m_InterBlockTria[nBlock][nOldSizeBor].vbFlipped.resize( tOldComp + 1) ; + bNewCompBor = false ; + } + size_t tCurrSz = m_InterBlockTria[nBlock][nOldSizeBor].vCompoTria.size() ; + m_InterBlockTria[nBlock][nOldSizeBor].vCompoTria[tCurrSz - 1].emplace_back( CurrTri) ; + m_InterBlockTria[nBlock][nOldSizeBor].vCompoTria[tCurrSz - 1].back().ToGlob( m_MapFrame) ; + size_t tTri = m_InterBlockTria[nBlock][nOldSizeBor].vbFlipped[tCurrSz - 1].size() ; + m_InterBlockTria[nBlock][nOldSizeBor].vbFlipped[tCurrSz - 1].resize( tTri + 1) ; + m_InterBlockTria[nBlock][nOldSizeBor].vbFlipped[tCurrSz - 1][tTri] = false ; + } + // Triangolo interno + else { + if ( bNewCompInn) { + Point3d ptVert = itVox->second.Compo[nComp].ptVert ; + ptVert.ToGlob( m_MapFrame) ; + triHold[nOldSizeInn].ptCompoVert.emplace_back( ptVert) ; + size_t tOldComp = triHold[nOldSizeInn].vCompoTria.size() ; + triHold[nOldSizeInn].vCompoTria.resize( tOldComp + 1) ; + triHold[nOldSizeInn].vbFlipped.resize( tOldComp + 1) ; + bNewCompInn = false ; + } + size_t tCurrSz = triHold[nOldSizeInn].vCompoTria.size() ; + triHold[nOldSizeInn].vCompoTria[tCurrSz - 1].emplace_back( CurrTri) ; + triHold[nOldSizeInn].vCompoTria[tCurrSz - 1].back().ToGlob( m_MapFrame) ; + size_t tTri = triHold[nOldSizeInn].vbFlipped[tCurrSz - 1].size() ; + triHold[nOldSizeInn].vbFlipped[tCurrSz - 1].resize( tTri + 1) ; + triHold[nOldSizeInn].vbFlipped[tCurrSz - 1][tTri] = false ; + } + } + } + } + + return true ; +} + + //---------------------------------------------------------------------------- // Calcola i triangoli dei voxel contenuti nel vettore, se si vuole il riconoscimento // delle sharp-feature passare true in bEhn. In questo caso vengono analizzai anche i // voxel adiacenti a quelli contenenti sharp-feature, e viene eseguito il flipping. bool -VolZmap::ExtMarchingCubes( vector& vVox, TRIA3DEXLIST& lstTria, bool bEnh) const +VolZmap::ExtMarchingCubes( vector& vVox, TRIA3DEXLIST& lstTria, bool bEnh) const { // Contenitore dei triangoli costituenti sharp-feaure TriHolder triHold ; @@ -2294,7 +2384,7 @@ VolZmap::ExtMarchingCubes( vector& vVox, TRIA3DEXLIST& lstTria, bool bEnh continue ; } if ( nOldV == nOriginalSize) { - Voxel NewVox ; + VoxelIndexes NewVox ; NewVox.nI = vVox[nV].nI + i ; NewVox.nJ = vVox[nV].nJ + j ; NewVox.nK = vVox[nV].nK + k ; @@ -2640,7 +2730,7 @@ VolZmap::FlipEdgesBB( TriaMatrix& InterTria) const bool VolZmap::IsThereMat( int nI, int nJ, int nK) const { - // Trasformo gl indici della griglia voxel in quelli della griglia dexel + // Trasformo gli indici della griglia voxel in quelli della griglia dexel nI *= N_DEXVOXRATIO ; nJ *= N_DEXVOXRATIO ; nK *= N_DEXVOXRATIO ; @@ -2716,7 +2806,7 @@ VolZmap::CalcIndex( int nI, int nJ, int nK) const //---------------------------------------------------------------------------- bool -VolZmap::IntersPos( int nVec1[], int nVec2[], bool bFirstCorner, VectorField& vfField) const +VolZmap::IntersPos( int nVec1[], int nVec2[], bool bFirstCorner, AppliedVector& vfField) const { const double dEps = 2 * EPS_SMALL ; @@ -2730,8 +2820,8 @@ VolZmap::IntersPos( int nVec1[], int nVec2[], bool bFirstCorner, VectorField& vf double dMinX = ( nMinI + 0.5) * m_dStep ; double dMaxX = ( nMaxI + 0.5) * m_dStep ; - vfField.ptInt.y = ( nVec1[1] * N_DEXVOXRATIO + 0.5) * m_dStep ; - vfField.ptInt.z = ( nVec1[2] * N_DEXVOXRATIO + 0.5) * m_dStep ; + vfField.ptPApp.y = ( nVec1[1] * N_DEXVOXRATIO + 0.5) * m_dStep ; + vfField.ptPApp.z = ( nVec1[2] * N_DEXVOXRATIO + 0.5) * m_dStep ; size_t nDexel = ( nVec1[2] * m_nNx[1] + nVec1[1]) * N_DEXVOXRATIO ; int nSize = int( m_Values[1][nDexel].size()) ; @@ -2741,9 +2831,9 @@ VolZmap::IntersPos( int nVec1[], int nVec2[], bool bFirstCorner, VectorField& vf double dX = m_Values[1][nDexel][n].dMax ; while ( n >= 0 && dX > dMinX - dEps) { if ( dX < dMaxX + dEps) { - vfField.ptInt.x = Clamp( dX, dMinX + EPS_SMALL, dMaxX - EPS_SMALL) ; - vfField.vtNorm = m_Values[1][nDexel][n].vtMaxN ; - vfField.nToolFlag = m_Values[1][nDexel][n].nToolMax ; + vfField.ptPApp.x = Clamp( dX, dMinX + EPS_SMALL, dMaxX - EPS_SMALL) ; + vfField.vtVec = m_Values[1][nDexel][n].vtMaxN ; + vfField.nPropIndex = m_Values[1][nDexel][n].nToolMax ; bFound = true ; break ; } @@ -2758,9 +2848,9 @@ VolZmap::IntersPos( int nVec1[], int nVec2[], bool bFirstCorner, VectorField& vf double dX = m_Values[1][nDexel][n].dMin ; while ( n < nSize && dX < dMaxX + dEps) { if ( dX > dMinX - dEps) { - vfField.ptInt.x = Clamp( dX, dMinX + EPS_SMALL, dMaxX - EPS_SMALL) ; - vfField.vtNorm = m_Values[1][nDexel][n].vtMinN ; - vfField.nToolFlag = m_Values[1][nDexel][n].nToolMin ; + vfField.ptPApp.x = Clamp( dX, dMinX + EPS_SMALL, dMaxX - EPS_SMALL) ; + vfField.vtVec = m_Values[1][nDexel][n].vtMinN ; + vfField.nPropIndex = m_Values[1][nDexel][n].nToolMin ; bFound = true ; break ; } @@ -2771,7 +2861,7 @@ VolZmap::IntersPos( int nVec1[], int nVec2[], bool bFirstCorner, VectorField& vf } if ( ! bFound) - vfField.ptInt.x = 0.5 * ( dMinX + dMaxX) ; + vfField.ptPApp.x = 0.5 * ( dMinX + dMaxX) ; } else if ( nVec1[1] != nVec2[1]) { @@ -2782,8 +2872,8 @@ VolZmap::IntersPos( int nVec1[], int nVec2[], bool bFirstCorner, VectorField& vf double dMinY = ( nMinJ + 0.5) * m_dStep ; double dMaxY = ( nMaxJ + 0.5) * m_dStep ; - vfField.ptInt.x = ( nVec1[0] * N_DEXVOXRATIO + 0.5) * m_dStep ; - vfField.ptInt.z = ( nVec1[2] * N_DEXVOXRATIO + 0.5) * m_dStep ; + vfField.ptPApp.x = ( nVec1[0] * N_DEXVOXRATIO + 0.5) * m_dStep ; + vfField.ptPApp.z = ( nVec1[2] * N_DEXVOXRATIO + 0.5) * m_dStep ; size_t nDexel = ( nVec1[0] * m_nNx[2] + nVec1[2]) * N_DEXVOXRATIO ; int nSize = int( m_Values[2][nDexel].size()) ; @@ -2793,9 +2883,9 @@ VolZmap::IntersPos( int nVec1[], int nVec2[], bool bFirstCorner, VectorField& vf double dY = m_Values[2][nDexel][n].dMax ; while ( n >= 0 && dY > dMinY - dEps) { if ( dY < dMaxY + dEps) { - vfField.ptInt.y = Clamp( dY, dMinY + EPS_SMALL, dMaxY - EPS_SMALL) ; - vfField.vtNorm = m_Values[2][nDexel][n].vtMaxN ; - vfField.nToolFlag = m_Values[2][nDexel][n].nToolMax; + vfField.ptPApp.y = Clamp( dY, dMinY + EPS_SMALL, dMaxY - EPS_SMALL) ; + vfField.vtVec = m_Values[2][nDexel][n].vtMaxN ; + vfField.nPropIndex = m_Values[2][nDexel][n].nToolMax; bFound = true ; break ; } @@ -2810,9 +2900,9 @@ VolZmap::IntersPos( int nVec1[], int nVec2[], bool bFirstCorner, VectorField& vf double dY = m_Values[2][nDexel][n].dMin ; while ( n < nSize && dY < dMaxY + dEps) { if ( dY > dMinY - dEps) { - vfField.ptInt.y = Clamp( dY, dMinY + EPS_SMALL, dMaxY - EPS_SMALL) ; - vfField.vtNorm = m_Values[2][nDexel][n].vtMinN ; - vfField.nToolFlag = m_Values[2][nDexel][n].nToolMin ; + vfField.ptPApp.y = Clamp( dY, dMinY + EPS_SMALL, dMaxY - EPS_SMALL) ; + vfField.vtVec = m_Values[2][nDexel][n].vtMinN ; + vfField.nPropIndex = m_Values[2][nDexel][n].nToolMin ; bFound = true ; break ; } @@ -2823,7 +2913,7 @@ VolZmap::IntersPos( int nVec1[], int nVec2[], bool bFirstCorner, VectorField& vf } if ( ! bFound) - vfField.ptInt.y = 0.5 * ( dMinY + dMaxY) ; + vfField.ptPApp.y = 0.5 * ( dMinY + dMaxY) ; } else if ( nVec1[2] != nVec2[2]) { @@ -2834,8 +2924,8 @@ VolZmap::IntersPos( int nVec1[], int nVec2[], bool bFirstCorner, VectorField& vf double dMinZ = ( nMinK + 0.5) * m_dStep ; double dMaxZ = ( nMaxK + 0.5) * m_dStep ; - vfField.ptInt.x = ( nVec1[0] * N_DEXVOXRATIO + 0.5) * m_dStep ; - vfField.ptInt.y = ( nVec1[1] * N_DEXVOXRATIO + 0.5) * m_dStep ; + vfField.ptPApp.x = ( nVec1[0] * N_DEXVOXRATIO + 0.5) * m_dStep ; + vfField.ptPApp.y = ( nVec1[1] * N_DEXVOXRATIO + 0.5) * m_dStep ; size_t nDexel = ( nVec1[1] * m_nNx[0] + nVec1[0]) * N_DEXVOXRATIO ; int nSize = int( m_Values[0][nDexel].size()) ; @@ -2845,9 +2935,9 @@ VolZmap::IntersPos( int nVec1[], int nVec2[], bool bFirstCorner, VectorField& vf double dZ = m_Values[0][nDexel][n].dMax ; while ( n >= 0 && dZ > dMinZ - dEps) { if ( dZ < dMaxZ + dEps) { - vfField.ptInt.z = Clamp( dZ, dMinZ + EPS_SMALL, dMaxZ - EPS_SMALL) ; - vfField.vtNorm = m_Values[0][nDexel][n].vtMaxN ; - vfField.nToolFlag = m_Values[0][nDexel][n].nToolMax ; + vfField.ptPApp.z = Clamp( dZ, dMinZ + EPS_SMALL, dMaxZ - EPS_SMALL) ; + vfField.vtVec = m_Values[0][nDexel][n].vtMaxN ; + vfField.nPropIndex = m_Values[0][nDexel][n].nToolMax ; bFound = true ; break ; } @@ -2862,9 +2952,9 @@ VolZmap::IntersPos( int nVec1[], int nVec2[], bool bFirstCorner, VectorField& vf double dZ = m_Values[0][nDexel][n].dMin ; while ( n < nSize && dZ < dMaxZ + dEps) { if ( dZ > dMinZ - dEps) { - vfField.ptInt.z = Clamp( dZ, dMinZ + EPS_SMALL, dMaxZ - EPS_SMALL) ; - vfField.vtNorm = m_Values[0][nDexel][n].vtMinN ; - vfField.nToolFlag = m_Values[0][nDexel][n].nToolMin ; + vfField.ptPApp.z = Clamp( dZ, dMinZ + EPS_SMALL, dMaxZ - EPS_SMALL) ; + vfField.vtVec = m_Values[0][nDexel][n].vtMinN ; + vfField.nPropIndex = m_Values[0][nDexel][n].nToolMin ; bFound = true ; break ; } @@ -2875,7 +2965,7 @@ VolZmap::IntersPos( int nVec1[], int nVec2[], bool bFirstCorner, VectorField& vf } if ( ! bFound) - vfField.ptInt.z = 0.5 * ( dMinZ + dMaxZ) ; + vfField.ptPApp.z = 0.5 * ( dMinZ + dMaxZ) ; } return bFound ; @@ -2902,8 +2992,8 @@ VolZmap::IsValidVoxel( int nI, int nJ, int nK) const int nVoxNumY = int( m_nNy[0] / N_DEXVOXRATIO + ( m_nNy[0] % N_DEXVOXRATIO == 0 ? 1 : 2)) ; int nVoxNumZ = int( m_nNy[1] / N_DEXVOXRATIO + ( m_nNy[1] % N_DEXVOXRATIO == 0 ? 1 : 2)) ; // Verifico la validità del voxel - return ( nI >= - 1 && nI < nVoxNumX - 1 && - nJ >= - 1 && nJ < nVoxNumY - 1 && + return ( nI >= - 1 && nI < nVoxNumX - 1 && + nJ >= - 1 && nJ < nVoxNumY - 1 && nK >= - 1 && nK < nVoxNumZ - 1 ) ; } @@ -2926,8 +3016,8 @@ VolZmap::GetVoxIJKFromN( int nN, int& nI, int& nJ, int& nK) const nK = ( nN / ( nVoxNumX * nVoxNumY)) - 1 ; // Controllo la sensatezza del risultato ottenuto - return ( nI >= - 1 && nI < nVoxNumX - 1 && - nJ >= - 1 && nJ < nVoxNumY - 1 && + return ( nI >= - 1 && nI < nVoxNumX - 1 && + nJ >= - 1 && nJ < nVoxNumY - 1 && nK >= - 1 && nK < nVoxNumZ - 1 ) ; } @@ -3227,9 +3317,9 @@ VolZmap::IsATriangleOnBorder( const Triangle3dEx& trTria, const Point3d& ptVert, //---------------------------------------------------------------------------- bool -VolZmap::ProcessVoxContXY( VoxelContainer& VoxContXY, bool bPlus, TRIA3DEXLIST& lstTria) const +VolZmap::ProcessVoxContXY( FlatVoxelContainer& VoxContXY, bool bPlus, TRIA3DEXLIST& lstTria) const { - VoxelContainer::iterator it = VoxContXY.begin() ; + FlatVoxelContainer::iterator it = VoxContXY.begin() ; while ( it != VoxContXY.end()) { int nN = ( *it).first ; @@ -3349,9 +3439,9 @@ VolZmap::ProcessVoxContXY( VoxelContainer& VoxContXY, bool bPlus, TRIA3DEXLIST& //---------------------------------------------------------------------------- bool -VolZmap::ProcessVoxContYZ( VoxelContainer& VoxContYZ, bool bPlus, TRIA3DEXLIST& lstTria) const +VolZmap::ProcessVoxContYZ( FlatVoxelContainer& VoxContYZ, bool bPlus, TRIA3DEXLIST& lstTria) const { - VoxelContainer::iterator it = VoxContYZ.begin() ; + FlatVoxelContainer::iterator it = VoxContYZ.begin() ; while ( it != VoxContYZ.end()) { int nN = ( *it).first ; @@ -3471,9 +3561,9 @@ VolZmap::ProcessVoxContYZ( VoxelContainer& VoxContYZ, bool bPlus, TRIA3DEXLIST& //---------------------------------------------------------------------------- bool -VolZmap::ProcessVoxContXZ( VoxelContainer& VoxContXZ, bool bPlus, TRIA3DEXLIST& lstTria) const +VolZmap::ProcessVoxContXZ( FlatVoxelContainer& VoxContXZ, bool bPlus, TRIA3DEXLIST& lstTria) const { - VoxelContainer::iterator it = VoxContXZ.begin() ; + FlatVoxelContainer::iterator it = VoxContXZ.begin() ; while ( it != VoxContXZ.end()) { int nN = ( *it).first ; @@ -3593,7 +3683,7 @@ VolZmap::ProcessVoxContXZ( VoxelContainer& VoxContXZ, bool bPlus, TRIA3DEXLIST& //---------------------------------------------------------------------------- bool -VolZmap::Find( const VoxelContainer& VoxCont, int nI, int nJ, int nK, double dPos, int nTool) const +VolZmap::Find( const FlatVoxelContainer& VoxCont, int nI, int nJ, int nK, double dPos, int nTool) const { // indice globale del voxel int nN ; @@ -3608,7 +3698,7 @@ VolZmap::Find( const VoxelContainer& VoxCont, int nI, int nJ, int nK, double dPo //---------------------------------------------------------------------------- bool -VolZmap::Remove( VoxelContainer& VoxCont, int nI, int nJ, int nK) const +VolZmap::Remove( FlatVoxelContainer& VoxCont, int nI, int nJ, int nK) const { int nN ; if ( ! GetVoxNFromIJK( nI, nJ, nK, nN))