From 2ecd82c61f97919aaa20af93e0908d50605aa2f3 Mon Sep 17 00:00:00 2001 From: Dario Sassi Date: Mon, 23 Oct 2017 06:47:22 +0000 Subject: [PATCH] EgtGeomKernel : - migliorata GetDepth di VolZmap, con flag per scelta algoritmo. --- EGkDllMain.cpp | 10 +- GeomDB.cpp | 6 +- HashGrids2d.cpp | 4 +- HashGrids3d.cpp | 4 +- NgeReader.cpp | 6 +- SurfFlatRegionBooleans.cpp | 4 +- SurfTriMeshFaceting.cpp | 2 +- Triangulate.cpp | 4 +- UserObjDefault.cpp | 2 +- VolTriZmapCalculus.cpp | 262 ++++--- VolTriZmapGraphics.cpp | 1425 ++++++++++++++++++++++++++++-------- VolZmap.h | 15 +- VolZmapTool.cpp | 2 +- 13 files changed, 1306 insertions(+), 440 deletions(-) diff --git a/EGkDllMain.cpp b/EGkDllMain.cpp index c05440c..b1cf275 100644 --- a/EGkDllMain.cpp +++ b/EGkDllMain.cpp @@ -20,6 +20,8 @@ #include "\EgtDev\Include\EgtTrace.h" #include "\EgtDev\Include\SELkLockId.h" +using namespace std ; + //--------------------------- Costanti ---------------------------------------- #if defined( _WIN64) #if defined( _DEBUG) @@ -112,7 +114,7 @@ SetEGkKeyType( int nType) } //----------------------------------------------------------------------------- -const std::string& +const string& GetEGkKey( void) { SetLockType( s_nKeyType) ; @@ -121,7 +123,7 @@ GetEGkKey( void) //----------------------------------------------------------------------------- void -InitFontManager( const std::string& sNfeFontDir, const std::string& sDefaultFont) +InitFontManager( const string& sNfeFontDir, const string& sDefaultFont) { // recupero il font manager FontManager& fntMgr = FontManager::GetFontManager() ; @@ -131,7 +133,7 @@ InitFontManager( const std::string& sNfeFontDir, const std::string& sDefaultFont } //----------------------------------------------------------------------------- -const std::string& +const string& GetNfeFontDir( void) { // recupero il font manager @@ -141,7 +143,7 @@ GetNfeFontDir( void) } //----------------------------------------------------------------------------- -const std::string& +const string& GetDefaultFont( void) { // recupero il font manager diff --git a/GeomDB.cpp b/GeomDB.cpp index 3a6bd40..e6ddaf3 100644 --- a/GeomDB.cpp +++ b/GeomDB.cpp @@ -2917,7 +2917,7 @@ GeomDB::AddMaterial( const string& sName, const Material& matM) //---------------------------------------------------------------------------- int -GeomDB::FindMaterial( const std::string& sName) const +GeomDB::FindMaterial( const string& sName) const { return m_MatManager.FindMaterial( sName) ; } @@ -2955,7 +2955,7 @@ GeomDB::GetMaterialData( int nMat, Material& matM) const //---------------------------------------------------------------------------- bool -GeomDB::GetMaterialName( int nMat, std::string& sName) const +GeomDB::GetMaterialName( int nMat, string& sName) const { return m_MatManager.GetMaterialName( nMat, sName) ; } @@ -2994,7 +2994,7 @@ GeomDB::ModifyMaterialData( int nMat, const Material& matM) //---------------------------------------------------------------------------- bool -GeomDB::ModifyMaterialName( int nMat, const std::string& sName) +GeomDB::ModifyMaterialName( int nMat, const string& sName) { return m_MatManager.ModifyMaterialName( nMat, sName) ; } diff --git a/HashGrids2d.cpp b/HashGrids2d.cpp index 6c83eb0..cbbbb50 100644 --- a/HashGrids2d.cpp +++ b/HashGrids2d.cpp @@ -44,7 +44,7 @@ class HashGrid2d : m_Objs( nullptr), m_neighborOffset( nullptr), m_occupiedCellsId( 0) {} } ; - typedef std::vector CellVector ; + typedef vector CellVector ; public : explicit HashGrid2d( double dCellSpan) ; @@ -695,7 +695,7 @@ HashGrids2d::addGrid( ObjData& obj) // If no hash grid yet exists in the hierarchy, an initial hash grid is created // based on the body's size. - pGrid = new HashGrid2d( size * std::sqrt( hierarchyFactor)) ; + pGrid = new HashGrid2d( size * sqrt( hierarchyFactor)) ; } else { // Check the hierarchy for a hash grid with suitably sized cells - if such a grid does not diff --git a/HashGrids3d.cpp b/HashGrids3d.cpp index a6f0532..a5009db 100644 --- a/HashGrids3d.cpp +++ b/HashGrids3d.cpp @@ -45,7 +45,7 @@ class HashGrid3d : m_Objs( nullptr), m_neighborOffset( nullptr), m_occupiedCellsId( 0) {} } ; - typedef std::vector CellVector ; + typedef vector CellVector ; public : explicit HashGrid3d( double dCellSpan) ; @@ -738,7 +738,7 @@ HashGrids3d::addGrid( ObjData& obj) // If no hash grid yet exists in the hierarchy, an initial hash grid is created // based on the body's size. - pGrid = new HashGrid3d( size * std::sqrt( hierarchyFactor)) ; + pGrid = new HashGrid3d( size * sqrt( hierarchyFactor)) ; } else { // Check the hierarchy for a hash grid with suitably sized cells - if such a grid does not diff --git a/NgeReader.cpp b/NgeReader.cpp index 7b58801..5877176 100644 --- a/NgeReader.cpp +++ b/NgeReader.cpp @@ -22,7 +22,7 @@ using namespace std ; //---------------------------------------------------------------------------- bool -NgeReader::Init( const std::string& sFileIn) +NgeReader::Init( const string& sFileIn) { switch ( NgeType( sFileIn)) { case NGE_ASCII : @@ -57,7 +57,7 @@ NgeReader::Close( void) //---------------------------------------------------------------------------- int -NgeReader::NgeType( const std::string& sFile) +NgeReader::NgeType( const string& sFile) { // apertura del file di ingresso ifstream InFile ; @@ -98,7 +98,7 @@ NgeReader::GetCurrPos( void) //---------------------------------------------------------------------------- bool -NgeReader::GetToken( std::string& sToken, const char* szSep, bool bEndL) +NgeReader::GetToken( string& sToken, const char* szSep, bool bEndL) { // se necessario, lettura nuova linea if ( m_iPosStart == string::npos) { diff --git a/SurfFlatRegionBooleans.cpp b/SurfFlatRegionBooleans.cpp index af13c4c..ae844d8 100644 --- a/SurfFlatRegionBooleans.cpp +++ b/SurfFlatRegionBooleans.cpp @@ -378,8 +378,8 @@ SurfFlatRegion* SurfFlatRegion::MyNewSurfFromLoops( PCRV_DEQUE& vpLoop) { // verifico che le curve siano chiuse e ne determino l'area - typedef std::pair INDAREA ; // coppia indice, area - typedef std::vector INDAREAVECTOR ; // vettore di coppie indice, area + typedef pair INDAREA ; // coppia indice, area + typedef vector INDAREAVECTOR ; // vettore di coppie indice, area INDAREAVECTOR vArea ; vArea.reserve( vpLoop.size()) ; for ( int i = 0 ; i < int( vpLoop.size()) ; ++ i) { diff --git a/SurfTriMeshFaceting.cpp b/SurfTriMeshFaceting.cpp index e9ed8fc..36f7e8b 100644 --- a/SurfTriMeshFaceting.cpp +++ b/SurfTriMeshFaceting.cpp @@ -653,7 +653,7 @@ SurfTriMesh::CloneFacet( int nF) const if ( ! pSTM->Init( nTria + 2, nTria, 1)) return nullptr ; // tabella hash per indici vertici (vecchi e nuovi indici) - typedef std::unordered_map VVMAP ; + typedef unordered_map VVMAP ; VVMAP PntMap ; PntMap.rehash( nTria + 2) ; // inserisco i triangoli nella nuova superficie diff --git a/Triangulate.cpp b/Triangulate.cpp index c7ed7a6..98de39d 100644 --- a/Triangulate.cpp +++ b/Triangulate.cpp @@ -702,8 +702,8 @@ bool Triangulate::SortInternalLoops( const POLYLINEVECTOR& vPL, INTVECTOR& vOrd) { // riempio vettore con indice e massimo specifico per ogni loop interno - typedef std::pair INDMAX ; // coppia indice, massimo - typedef std::vector INDMAXVECTOR ; // vettore di coppie indice, massimo + typedef pair INDMAX ; // coppia indice, massimo + typedef vector INDMAXVECTOR ; // vettore di coppie indice, massimo INDMAXVECTOR vMax ; vMax.reserve( vPL.size()) ; vMax.emplace_back( 0, 0.0) ; diff --git a/UserObjDefault.cpp b/UserObjDefault.cpp index 1820fae..039f9b0 100644 --- a/UserObjDefault.cpp +++ b/UserObjDefault.cpp @@ -54,7 +54,7 @@ UserObjDefault::Clone( void) const } //---------------------------------------------------------------------------- -const std::string& +const string& UserObjDefault::GetClassName( void) const { return m_sName ; diff --git a/VolTriZmapCalculus.cpp b/VolTriZmapCalculus.cpp index 4aa8dc6..9a8430e 100644 --- a/VolTriZmapCalculus.cpp +++ b/VolTriZmapCalculus.cpp @@ -79,6 +79,19 @@ VolZmap::IntersLineBox( const Point3d& ptP, const Vector3d& vtV, return ( dU2 >= dU1) ; } +//---------------------------------------------------------------------------- +// La retta è rappresentata da un punto e dal versore espressi nel riferimento dello Zmap. +// Se non vi è intersezione, viene restituito falso +bool +VolZmap::IntersLineZMapLattice( const Point3d& ptP, const Vector3d& vtV, double& dU1, double& dU2) const +{ + // Punti estremi della diagonale del box del reticolo + Point3d ptMin( - 0.5 * m_dStep, - 0.5 * m_dStep, - 0.5 * m_dStep) ; + Point3d ptMax( ( m_nNx[0] + 0.5) * m_dStep, ( m_nNy[0] + 0.5) * m_dStep, ( m_nNy[1] + 0.5) * m_dStep) ; + + return IntersLineBox( ptP, vtV, ptMin, ptMax, dU1, dU2) ; +} + //---------------------------------------------------------------------------- bool VolZmap::IntersLineZMapBBox( const Point3d& ptP, const Vector3d& vtV, unsigned int nGrid, double& dU1, double& dU2) const @@ -130,7 +143,7 @@ VolZmap::IntersLineDexel( const Point3d& ptP, const Vector3d& vtV, unsigned int //---------------------------------------------------------------------------- bool -VolZmap::IntersRayDexel( const Point3d& ptP, const Vector3d& vtV, unsigned int nGrid, unsigned int nI, unsigned int nJ, +VolZmap::IntersRayDexel( const Point3d& ptP, const Vector3d& vtV, unsigned int nGrid, unsigned int nI, unsigned int nJ, double& dU1, double& dU2) const { // Determino l'indice del dexel e il numero di suoi intervalli @@ -172,25 +185,34 @@ VolZmap::IntersRayDexel( const Point3d& ptP, const Vector3d& vtV, unsigned int n } //---------------------------------------------------------------------------- +// Calcola la profondità del materiale lungo un raggio, identificato da un punto e un versoreore. +// Punto e versore devono essere espressi nel sistema di riferimento locale, in cui è immerso quello dello Zmap. // InLength = distanza di ingresso (se -1 il punto è interno, se -2 il punto è esterno e il raggio non interseca lo Zmap) // OutLength = distanza di uscita bool -VolZmap::GetDepth( const Point3d& ptPGlob, const Vector3d& vtDir, double& dInLength, double& dOutLength) const +VolZmap::GetDepth( const Point3d& ptPLoc, const Vector3d& vtDLoc, double& dInLength, double& dOutLength, bool bExact) const { -#if 1 - GetDepthWithVoxel( ptPGlob, vtDir, dInLength, dOutLength) ; - return true ; -#else + Point3d ptP = ptPLoc ; + Vector3d vtD = vtDLoc ; + ptP.ToLoc( m_MapFrame[0]) ; + vtD.ToLoc( m_MapFrame[0]) ; + if ( bExact) + return GetDepthWithVoxel( ptP, vtD, dInLength, dOutLength, true) ; + else + return GetDepthWithDexel( ptP, vtD, dInLength, dOutLength) ; +} + +//---------------------------------------------------------------------------- +// Calcola la profondità del materiale usando i parallelepipedi corrispondenti ai segmenti di dexel di una griglia. +// Punto e versore devono essere espressi nel sistema di riferimento Zmap. Il vettore deve essere normalizzato. +// InLength = distanza di ingresso (se -1 il punto è interno, se -2 il punto è esterno e il raggio non interseca lo Zmap) +// OutLength = distanza di uscita +bool +VolZmap::GetDepthWithDexel( const Point3d& ptP, const Vector3d& vtV, double& dInLength, double& dOutLength) const +{ int nGrid = 0 ; - // Porto il raggio nel riferimento intrinseco - Point3d ptP = ptPGlob ; - ptP.ToLoc( m_MapFrame[nGrid]) ; - Vector3d vtV = vtDir ; - vtV.ToLoc( m_MapFrame[nGrid]) ; - vtV.Normalize() ; - // Intersezione fra semiretta e BBox dello Zmap double dU1, dU2 ; bool bTest = IntersLineZMapBBox( ptP, vtV, nGrid, dU1, dU2) && ( dU1 > 0 || dU2 > 0) ; @@ -286,23 +308,17 @@ VolZmap::GetDepth( const Point3d& ptPGlob, const Vector3d& vtDir, double& dInLen dInLength = - 1 ; return true ; -#endif } //---------------------------------------------------------------------------- -// Punti e vettori devono essere espressi nel sistema locale (quello in cui è immerso lo Zmap) +// Calcola la profondità del materiale usando i triangoli della rappresentazione grafica. +// Punto e versore devono essere espressi nel sistema di riferimento Zmap. +// Se si vuole che vengano utilizzati i triangoli con sharp-featue bEnh deve valere true, false altrimenti. // InLength = distanza di ingresso (se -1 il punto è interno, se -2 il punto è esterno e il raggio non interseca lo Zmap) -// OutLength = distanza di uscita +// OutLength = distanza di uscita. bool -VolZmap::GetDepthWithVoxel( const Point3d& ptPLoc, const Vector3d& vtDir, double& dInLength, double& dOutLength) const +VolZmap::GetDepthWithVoxel( const Point3d& ptP, const Vector3d& vtD, double& dInLength, double& dOutLength, bool bEnh) const { - // Porto punto e vettore della retta nel sistema Zmap - Point3d ptP = ptPLoc ; - ptP.ToLoc( m_MapFrame[0]) ; - Vector3d vtD = vtDir ; - vtD.ToLoc( m_MapFrame[0]) ; - vtD.Normalize() ; - // Parametri di intersezione retta BBox dello Zmap double dU1, dU2 ; @@ -322,106 +338,120 @@ VolZmap::GetDepthWithVoxel( const Point3d& ptPLoc, const Vector3d& vtDir, double if ( ! GetPointVoxel( ptP + dU1 * vtD, nVoxI, nVoxJ, nVoxK)) return false ; } - - // Indici dell'ultimo voxel - int nVxEndI, nVxEndJ, nVxEndK ; - // Determino il voxel finale - if ( ! GetPointVoxel( ptP + dU2 * vtD, nVxEndI, nVxEndJ, nVxEndK)) - return false ; + + // Vettore di voxel + vector vVox ; - // Struttura studio dell'intersezione + // Incrementi degli indici per sguire la retta + int nDeltaI = ( vtD.x > 0 ? 1 : ( vtD.x < 0 ? - 1 : 0)) ; + int nDeltaJ = ( vtD.y > 0 ? 1 : ( vtD.y < 0 ? - 1 : 0)) ; + int nDeltaK = ( vtD.z > 0 ? 1 : ( vtD.z < 0 ? - 1 : 0)) ; + // Scelgo i piani del voxel con cui intersecare la retta, + // per determinare il voxel successivo + int nPlaneI = ( vtD.x >= 0 ? 1 : 0) ; + int nPlaneJ = ( vtD.y >= 0 ? 1 : 0) ; + int nPlaneK = ( vtD.z >= 0 ? 1 : 0) ; + + // Ciclo sui voxel + while ( IsValidVoxel( nVoxI, nVoxJ, nVoxK)) { + // Estremi della diagonale del voxel corrente + Point3d ptMin( ( nVoxI + 0.5) * m_dStep, + ( nVoxJ + 0.5) * m_dStep, + ( nVoxK + 0.5) * m_dStep) ; + Point3d ptMax( ( nVoxI + 1.5) * m_dStep, + ( nVoxJ + 1.5) * m_dStep, + ( nVoxK + 1.5) * m_dStep) ; + // Studio il voxel corrente + if ( IntersLineBox( ptP, vtD, ptMin, ptMax)) { + Voxel NewVox ; + NewVox.nI = nVoxI ; + NewVox.nJ = nVoxJ ; + NewVox.nK = nVoxK ; + vVox.emplace_back( NewVox) ; + } + // Interseco la retta con i piani frontiera del voxel + double dMaxTX = ( abs( vtD.x) > EPS_ZERO ? abs( ( ( nVoxI + nPlaneI + 0.5) * m_dStep - ptP.x) / vtD.x) : INFINITO) ; + double dMaxTY = ( abs( vtD.y) > EPS_ZERO ? abs( ( ( nVoxJ + nPlaneJ + 0.5) * m_dStep - ptP.y) / vtD.y) : INFINITO) ; + double dMaxTZ = ( abs( vtD.z) > EPS_ZERO ? abs( ( ( nVoxK + nPlaneK + 0.5) * m_dStep - ptP.z) / vtD.z) : INFINITO) ; + + if ( dMaxTX < dMaxTY) { + if ( dMaxTX < dMaxTZ) + nVoxI += nDeltaI ; + else + nVoxK += nDeltaK ; + } + else { + if ( dMaxTY < dMaxTZ) + nVoxJ += nDeltaJ ; + else + nVoxK += nDeltaK ; + } + } + + TRIA3DLIST lstTria ; + ExtMarchingCubes( vVox, lstTria, bEnh) ; + + // Struttura studio dell'intersezione struct LineTriaInt { double dPar ; double dDot ; } ; - std::vector> IntMatrix ; + vector> IntMatrix ; - int nStI = min( nVoxI, nVxEndI) ; - int nStJ = min( nVoxJ, nVxEndJ) ; - int nStK = min( nVoxK, nVxEndK) ; - int nEnI = max( nVoxI, nVxEndI) ; - int nEnJ = max( nVoxJ, nVxEndJ) ; - int nEnK = max( nVoxK, nVxEndK) ; - - for ( int nI = nStI ; nI <= nEnI ; ++ nI) { - for ( int nJ = nStJ ; nJ <= nEnJ ; ++ nJ) { - for ( int nK = nStK ; nK <= nEnK ; ++ nK) { - - Point3d ptMin( ( nI + 0.5) * m_dStep, - ( nJ + 0.5) * m_dStep, - ( nK + 0.5) * m_dStep) ; - Point3d ptMax( ( nI + 1.5) * m_dStep, - ( nJ + 1.5) * m_dStep, - ( nK + 1.5) * m_dStep) ; - - if ( IsValidVoxel( nI, nJ, nK) && - IntersLineBox( ptP, vtD, ptMin, ptMax)) { - - // Analisi del voxel - int nCbType ; - TRIA3DLIST lstTria ; - ProcessCube( nI, nJ, nK, nCbType, lstTria) ; - - // Se il voxel contiene triangoli - if ( nCbType == VOX_ON_BOUNDARY) { - - // Ciclo sui triangoli del voxel - for ( auto it = lstTria.begin() ; it != lstTria.end() ; ++it ) { - - Triangle3d CurrTria = *it ; - Point3d ptLineTria1, ptLineTria2 ; - // Studio dell'intersezione della retta con il triangolo corrente - int nIntType = IntersLineTria( ptP, vtD, 1.5 * dU2, CurrTria, ptLineTria1, ptLineTria2) ; - // Se non ci sono intersezioni passo al prossimo triangolo - if ( nIntType == ILTT_NO) - continue ; - // se altrimenti c'è una sola intersezione - else if ( nIntType == ILTT_VERT || - nIntType == ILTT_EDGE || - nIntType == ILTT_IN) { - LineTriaInt NewInt ; - NewInt.dPar = ( ptLineTria1 - ptP) * vtD ; - NewInt.dDot = vtD * CurrTria.GetN() ; - std::vector vSing ; - vSing.emplace_back( NewInt) ; + // Ciclo sui triangoli del voxel + for ( auto it = lstTria.begin() ; it != lstTria.end() ; ++it ) { + // Triangolo corrente e suoi punti di intersezione con la retta + Triangle3d CurrTria = *it ; + Point3d ptLineTria1, ptLineTria2 ; + // Studio dell'intersezione della retta con il triangolo corrente + int nIntType = IntersLineTria( ptP, vtD, 1.5 * dU2, CurrTria, ptLineTria1, ptLineTria2) ; + // Se non ci sono intersezioni passo al prossimo triangolo + if ( nIntType == ILTT_NO) + continue ; + // se altrimenti c'è una sola intersezione + else if ( nIntType == ILTT_VERT || + nIntType == ILTT_EDGE || + nIntType == ILTT_IN) { + LineTriaInt NewInt ; + NewInt.dPar = ( ptLineTria1 - ptP) * vtD ; + NewInt.dDot = vtD * CurrTria.GetN() ; + vector vSing ; + vSing.emplace_back( NewInt) ; - IntMatrix.emplace_back( vSing) ; - } - // altrimenti ci sono due intersezioni - else { - LineTriaInt NewInt1, NewInt2 ; - std::vector vCouple ; - double dP1 = ( ptLineTria1 - ptP) * vtD ; - double dP2 = ( ptLineTria2 - ptP) * vtD ; - NewInt1.dPar = ( dP1 < dP2 ? dP1 : dP2) ; - NewInt2.dPar = ( dP1 < dP2 ? dP2 : dP1) ; - NewInt1.dDot = vtD * CurrTria.GetN() ; - NewInt2.dDot = NewInt1.dDot ; + IntMatrix.emplace_back( vSing) ; + } + // altrimenti ci sono due intersezioni + else { + LineTriaInt NewInt1, NewInt2 ; + vector vCouple ; + double dP1 = ( ptLineTria1 - ptP) * vtD ; + double dP2 = ( ptLineTria2 - ptP) * vtD ; + NewInt1.dPar = ( dP1 < dP2 ? dP1 : dP2) ; + NewInt2.dPar = ( dP1 < dP2 ? dP2 : dP1) ; + NewInt1.dDot = vtD * CurrTria.GetN() ; + NewInt2.dDot = NewInt1.dDot ; - vCouple.emplace_back( NewInt1) ; - vCouple.emplace_back( NewInt2) ; - IntMatrix.emplace_back( vCouple) ; - } - } - } - } - } + vCouple.emplace_back( NewInt1) ; + vCouple.emplace_back( NewInt2) ; + IntMatrix.emplace_back( vCouple) ; } } - + // Ordino le intersezioni in base al parametro distanza con segno da ptP for ( int nN = 0 ; nN < int( IntMatrix.size()) - 1 ; ++ nN) { + for ( int nM = nN + 1 ; nM < int( IntMatrix.size()) ; ++ nM) { - double dFP = ( IntMatrix[nN].size() == 2 ? 0.5 * ( IntMatrix[nN][0].dPar + IntMatrix[nN][1].dPar) : - IntMatrix[nN][0].dPar) ; - double dLP = ( IntMatrix[nN+1].size() == 2 ? 0.5 * ( IntMatrix[nN+1][0].dPar + IntMatrix[nN+1][1].dPar) : - IntMatrix[nN+1][0].dPar) ; + double dFP = ( IntMatrix[nN].size() == 2 ? 0.5 * ( IntMatrix[nN][0].dPar + IntMatrix[nN][1].dPar) : + IntMatrix[nN][0].dPar) ; + double dLP = ( IntMatrix[nM].size() == 2 ? 0.5 * ( IntMatrix[nM][0].dPar + IntMatrix[nM][1].dPar) : + IntMatrix[nM][0].dPar) ; - if ( dFP > dLP) - swap( IntMatrix[nN], IntMatrix[nN+1]) ; + if ( dFP > dLP) + swap( IntMatrix[nN], IntMatrix[nM]) ; + } } - std::vector vInt ; + vector vInt ; for ( int nN = 0 ; nN < int( IntMatrix.size()) ; ++ nN) { vInt.emplace_back( IntMatrix[nN][0]) ; @@ -437,7 +467,7 @@ VolZmap::GetDepthWithVoxel( const Point3d& ptPLoc, const Vector3d& vtDir, double int nFirstPosN ; int nN = 0 ; for ( ; nN < int( vInt.size()) ; ++ nN) { - if ( vInt[nN].dPar >= 0) { + if ( /*vInt[nN].dPar >= 0*/vInt[nN].dPar > - EPS_SMALL) { // SCEGLIERE FRA LE DUE OPZIONI nFirstPosN = nN ; break ; } @@ -454,8 +484,14 @@ VolZmap::GetDepthWithVoxel( const Point3d& ptPLoc, const Vector3d& vtDir, double dInLength = -1 ; } else if ( nFirstPosN == 0) { - if ( vInt[nFirstPosN].dDot > - EPS_ZERO) - dInLength = -1 ; + + if ( vInt[nFirstPosN].dDot > EPS_ZERO) // + dInLength = -1 ; + else if ( GetPointVoxel( ptP, nVoxI, nVoxJ, nVoxK)) { + int nCubeType = CalcIndex( nVoxI, nVoxJ, nVoxK) ; + if ( nCubeType == 255) + dInLength = -1 ; + } } for ( int nN = nFirstPosN ; nN < int( vInt.size()) ; ++ nN) { @@ -1001,8 +1037,8 @@ VolZmap::IntersLineEllipticalCylinder( const Vector3d& vtLineDir, const Point3d& ptP.ToLoc( CircFrame) ; vtV.ToLoc( CircFrame) ; - std::vector vdCoef(3) ; - std::vector vdRoots ; + vector vdCoef(3) ; + vector vdRoots ; vdCoef[0] = dSqCoef * ptP.x * ptP.x + ptP.y * ptP.y + ptP.z * ptP.z - 2 * dObCoef * ptP.x * ptP.y - dSqRad ; vdCoef[1] = 2 * ( dSqCoef * vtV.x * ptP.x + vtV.y * ptP.y + vtV.z * ptP.z - dObCoef * ( vtV.x * ptP.y + vtV.y * ptP.x)) ; diff --git a/VolTriZmapGraphics.cpp b/VolTriZmapGraphics.cpp index 5c6b03b..ab37e44 100644 --- a/VolTriZmapGraphics.cpp +++ b/VolTriZmapGraphics.cpp @@ -351,7 +351,7 @@ VolZmap::GetTriangles( bool bAllBlocks, INTVECTOR& nModifiedBlocks, TRIA3DLISTVE #if 1 ExtMarchingCubes( int( t), vLstTria.back(), VecTriHold[t]) ; // Flipping fra voxel interni - FlipEdgesII( VecTriHold[t]) ; + FlipEdgesII( VecTriHold[t], true) ; bCalcInterBlock = true ; #else MarchingCubes( int( t), vLstTria.back()) ; @@ -615,9 +615,11 @@ VolZmap::AddDexelSideFace( int nPos, int nPosAdj, const Point3d& ptP, const Poin } //---------------------------------------------------------------------------- -// I triangoli sono espressi nel sistema Zmap +// Calcola i triangoli del voxel con indici nVoxI, nVoxJ, nVoxK; se si vuole +// il riconoscimento delle sharp-feature bEnh deve valere true. +// I triangoli formanti sharp-feature vengono messi nel TriHolder, gli altri nella lista. bool -VolZmap::ProcessCube( int nVoxI, int nVoxJ, int nVoxK, int& nCubeIndex, TRIA3DLIST& lstTria) const +VolZmap::ProcessCube( int nVoxI, int nVoxJ, int nVoxK, TRIA3DLIST& lstTria, TriHolder& triHold, bool bEnh) const { // Se il voxel non esiste, vi è un errore. if ( nVoxI + 1 < 0 || nVoxI + 1 > int( m_nNx[0]) || @@ -625,17 +627,13 @@ VolZmap::ProcessCube( int nVoxI, int nVoxJ, int nVoxK, int& nCubeIndex, TRIA3DLI nVoxK + 1 < 0 || nVoxK + 1 > int( m_nNy[1])) return false ; - // Classificazione dei vertici: interni o esterni al materiale o di frontiera - int nIndex = CalcIndex( nVoxI, nVoxJ, nVoxK) ; - if ( nIndex == 0) - nCubeIndex = VOX_EXTERN ; - else if ( nIndex == 255) - nCubeIndex = VOX_INNER ; - else - nCubeIndex = VOX_ON_BOUNDARY ; - - if ( nCubeIndex == VOX_ON_BOUNDARY) { + // Classificazione dei vertici: interni o esterni al materiale + int nIndex = CalcIndex( nVoxI, nVoxJ , nVoxK) ; + // Se vi è qualche intersezione fra segmenti e superficie + // continuo altrimenti passo al prossimo voxel. + if ( EdgeTable[nIndex] != 0) { + // Indici i,j,k dei vertici int IndexCorner[8][3] = { { nVoxI, nVoxJ, nVoxK}, @@ -676,220 +674,954 @@ VolZmap::ProcessCube( int nVoxI, int nVoxJ, int nVoxK, int& nCubeIndex, TRIA3DLI bReg = false ; } - // Determino il numero di componenti connesse nel voxel in caso di configurazione standard + // Determino il numero di componenti connesse nel voxel in caso di configurazione standard int nComponents = TriangleTableEn[nIndex][1][0] ; - // Matrici di campi vettoriali: - // CompoVert[i] ha i vertici della base del triangle fan della (i+1)-esima componente connessa; - // CompoTriVert[i] ha i vertici di tutti i triangoli, nel caso di assenza di sharp feature, - // della (i+1)-esima componente connessa. + // Matrici di campi vettoriali: + // CompoVert[i] ha i vertici della base del triangle fan della (i+1)-esima componente connessa; + // CompoTriVert[i] ha i vertici di tutti i triangoli, nel caso di assenza di sharp feature, + // della (i+1)-esima componente connessa. VectorField CompoVert[6][12] ; VectorField CompoTriVert[6][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. + // 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 nExtTabOff = nComponents ; int nStdTabOff = 0 ; - // Carico le matrici CompoVert e CompoTriVert + // Carico le matrici CompoVert e CompoTriVert for ( int nCompCount = 1 ; nCompCount <= nComponents ; ++ nCompCount) { - // Numero vertici per componenti + // Numero vertici per componenti nVertComp[nCompCount - 1] = TriangleTableEn[nIndex][1][nCompCount] ; - - // 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. + + // Riempio il nCompCount-esimo vettore di vertici della base del fan + for ( int nVertCount = 0 ; nVertCount < nVertComp[nCompCount - 1] ; ++ nVertCount) + CompoVert[nCompCount - 1][nVertCount] = VecField[TriangleTableEn[nIndex][1][nVertCount + nExtTabOff + 1]] ; + + // Serve per la gestione del caso ... + if ( nVertComp[nCompCount - 1] == 4) { + for ( int nVertCount = 0 ; nVertCount < nVertComp[nCompCount - 1] ; ++ nVertCount) + nIndArrey[nCompCount - 1][nVertCount] = TriangleTableEn[nIndex][1][nVertCount + nExtTabOff + 1] ; + } + + // Riempio il nCompCount-esimo vettore di vertici dei triangoli in assenza di + // sharp feature: in una mesh di triangoli con n vertici vi sono n - 2 triangoli. for ( int nVertCount = 0 ; nVertCount < 3 * ( nVertComp[nCompCount - 1] - 2) ; nVertCount += 3) { CompoTriVert[nCompCount - 1][nVertCount] = VecField[TriangleTableEn[nIndex][0][nStdTabOff + nVertCount+2]] ; CompoTriVert[nCompCount - 1][nVertCount+1] = VecField[TriangleTableEn[nIndex][0][nStdTabOff + nVertCount+1]] ; CompoTriVert[nCompCount - 1][nVertCount+2] = VecField[TriangleTableEn[nIndex][0][nStdTabOff + nVertCount]] ; } - // Aggiorno gli offsets per raggiungere i - // vertici della componente successiva. + // Aggiorno gli offsets per raggiungere i + // vertici della componente successiva. nExtTabOff += nVertComp[nCompCount - 1] ; nStdTabOff += 3 * ( nVertComp[nCompCount - 1] - 2) ; } - // Test sulla topologia: dal momento che il nostro test - // si fonda sugli angoli compresi fra le normali, esso ha - // senso solo se il campo è regolare. - if ( bReg && nAllConfig[nIndex] == 3) { + // 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 ; - Vector3d vtCmpAvg0, vtCmpAvg1 ; + if ( bReg) { + if ( nAllConfig[nIndex] == 3) { - // 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). - bool bTest0 = DotTest( CompoVert[0], 3, vtCmpAvg0, 0.7) ; - bool bTest1 = DotTest( CompoVert[1], 3, vtCmpAvg1, 0.7) ; + Vector3d vtCmpAvg0, vtCmpAvg1 ; - // Se i versori di entrambe le componenti sono concordi - // ha senso parlare di vettori medi, altrimenti non ha - // senso. Se non ha senso parlare di vettori medi non - // ha senso parlare di prodotti scalari fra loro, - // quindi pongo il loro prodotto a un valore assurdo -2 - // (il prodotto scalare fra versori ha modulo non superiore - // a uno). - double dScProd = - 2 ; + // 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). + bool bTest0 = DotTest( CompoVert[0], 3, vtCmpAvg0, 0.7) ; + bool bTest1 = DotTest( CompoVert[1], 3, vtCmpAvg1, 0.7) ; - if ( bTest0 && bTest1) - dScProd = vtCmpAvg0 * vtCmpAvg1 ; + // Se i versori di entrambe le componenti sono concordi + // ha senso parlare di vettori medi, altrimenti non ha + // senso. Se non ha senso parlare di vettori medi non + // ha senso parlare di prodotti scalari fra loro, + // quindi pongo il loro prodotto a un valore assurdo -2 + // (il prodotto scalare fra versori ha modulo non superiore + // a uno). + double dScProd = - 2 ; - double dThreshold = 0.7 ; + if ( bTest0 && bTest1) + dScProd = vtCmpAvg0 * vtCmpAvg1 ; - if ( ! ( bTest0 && bTest1) || ( bTest0 && bTest1 && dScProd > dThreshold)) { + double dThreshold = 0.7 ; - int nt = 0 ; + if ( ! ( bTest0 && bTest1) || ( bTest0 && bTest1 && dScProd > dThreshold)) { - while ( nIndexVsIndex3[nt][0] != nIndex) - ++ nt ; + int nt = 0 ; - int nRotCase = nIndexVsIndex3[nt][1] ; + while ( nIndexVsIndex3[nt][0] != nIndex) + ++ nt ; - nComponents = Cases3Plus[nRotCase][1][0] ; + int nRotCase = nIndexVsIndex3[nt][1] ; - // Riaggiorno gli offsets - nExtTabOff = nComponents ; - nStdTabOff = 0 ; + nComponents = Cases3Plus[nRotCase][1][0] ; - // Modifico le matrici - for ( int nC = 1 ; nC <= nComponents ; ++ nC) { + // Riaggiorno gli offsets + nExtTabOff = nComponents ; + nStdTabOff = 0 ; - // Numero vertici per componenti - nVertComp[nC - 1] = Cases3Plus[nRotCase][1][nC] ; + // Modifico le matrici + for ( int nC = 1 ; nC <= nComponents ; ++ nC) { - // Matrici dei vertici dei triangoli in assenza di sharp feature - for ( int nTriVert = 0 ; nTriVert < 3 * ( nVertComp[nC - 1] - 2) ; nTriVert += 3) { + // Numero vertici per componenti + nVertComp[nC - 1] = Cases3Plus[nRotCase][1][nC] ; - 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]] ; + // 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) ; } + } + } - // Aggiorno gli offsets per raggiungere i - // vertici della componente successiva. - nExtTabOff += nVertComp[nC - 1] ; - nStdTabOff += 3 * ( nVertComp[nC - 1] - 2) ; + 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 ; + + double dThreshold = 0.7 ; + + if ( ! ( bTest0 && bTest1) || ( bTest0 && bTest1 && dScProd > dThreshold)) { + + int nt = 0 ; + + while ( nIndexVsIndex6[nt][0] != nIndex) + ++ nt ; + + int nRotCase = nIndexVsIndex6[nt][1] ; + + 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 ; + } + } + + 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 + bool bTest0 = DotTest( CompoVert[0], 4, vtCmpAvg0) ; + bool bTest1 = DotTest( CompoVert[1], 4, vtCmpAvg1) ; + + if ( ! ( bTest0 && bTest1)) { + + 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) ; + } } } } - else if ( bReg && nAllConfig[nIndex] == 6) { + // Numero di feature nel voxel: al più vi è una feature per componente connessa. + int nInnerFeatureInVoxel = 0 ; + int nBorderFeatureInVoxel = 0 ; - // 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 ; - - if ( ! ( bTest0 && bTest1) || ( bTest0 && bTest1 && dScProd > dThreshold)) { - - int nt = 0 ; - - while ( nIndexVsIndex6[nt][0] != nIndex) - ++ nt ; - - int nRotCase = nIndexVsIndex6[nt][1] ; - - 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) ; - } - } - } - - else if ( bReg && 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 - bool bTest0 = DotTest( CompoVert[0], 4, vtCmpAvg0) ; - bool bTest1 = DotTest( CompoVert[1], 4, vtCmpAvg1) ; - - if ( ! ( bTest0 && bTest1)) { - - 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] ; - - - // 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) ; - } - } - } - // Ciclo sulle componenti + // Ciclo sulle componenti for ( int nCompCount = 1 ; nCompCount <= nComponents ; ++ nCompCount) { - // Costruzione dei triangoli - for ( int TriIndex = 0; TriIndex < ( nVertComp[nCompCount - 1] - 2) * 3 ; TriIndex += 3) { - // Il triangolo è pronto - Triangle3d CurrentTriangle ; - CurrentTriangle.Set( CompoTriVert[nCompCount - 1][TriIndex].ptInt, - CompoTriVert[nCompCount - 1][TriIndex+1].ptInt, - CompoTriVert[nCompCount - 1][TriIndex+2].ptInt) ; - bool bV = CurrentTriangle.Validate( true) ; - // Aggiungo alla lista - lstTria.emplace_back( CurrentTriangle) ; - } + + int nFeatureType = NO_FEATURE ; + // Se i componenti sono regolari valuto le normali per stabilire se eseguire ExtMC o MC + if ( bReg) + nFeatureType = TestOnNormal( CompoVert[nCompCount - 1], nVertComp[nCompCount - 1]) ; + + // Controllo per il caso piano su una griglia + // con versori normali a due a due paralleli. + bool bGridControl = true ; + + if ( nFeatureType != NO_FEATURE) { + + if ( nVertComp[nCompCount - 1] == 4) { + + // Ordino i 4 indici in senso crescente + for ( int nSrtInd1 = 0 ; nSrtInd1 < nVertComp[nCompCount - 1] - 1 ; ++ nSrtInd1) { + for ( int nSrtInd2 = nSrtInd1 + 1 ; nSrtInd2 < nVertComp[nCompCount - 1] ; ++ nSrtInd2) { + if ( nIndArrey[nCompCount - 1][nSrtInd1] > nIndArrey[nCompCount - 1][nSrtInd2]) + swap( nIndArrey[nCompCount - 1][nSrtInd1], nIndArrey[nCompCount - 1][nSrtInd2]) ; + } + } + + if ( ( nIndArrey[nCompCount - 1][0] == 0 && nIndArrey[nCompCount - 1][1] == 2 && + nIndArrey[nCompCount - 1][2] == 9 && nIndArrey[nCompCount - 1][3] == 10) || + ( nIndArrey[nCompCount - 1][0] == 4 && nIndArrey[nCompCount - 1][1] == 6 && + nIndArrey[nCompCount - 1][2] == 9 && nIndArrey[nCompCount - 1][3] == 10) || + ( nIndArrey[nCompCount - 1][0] == 4 && nIndArrey[nCompCount - 1][1] == 6 && + nIndArrey[nCompCount - 1][2] == 8 && nIndArrey[nCompCount - 1][3] == 11) || + ( nIndArrey[nCompCount - 1][0] == 0 && nIndArrey[nCompCount - 1][1] == 2 && + nIndArrey[nCompCount - 1][2] == 8 && nIndArrey[nCompCount - 1][3] == 11) || + ( nIndArrey[nCompCount - 1][0] == 1 && nIndArrey[nCompCount - 1][1] == 3 && + nIndArrey[nCompCount - 1][2] == 8 && nIndArrey[nCompCount - 1][3] == 9 ) || + ( nIndArrey[nCompCount - 1][0] == 1 && nIndArrey[nCompCount - 1][1] == 3 && + nIndArrey[nCompCount - 1][2] == 10 && nIndArrey[nCompCount - 1][3] == 11) || + ( nIndArrey[nCompCount - 1][0] == 5 && nIndArrey[nCompCount - 1][1] == 7 && + nIndArrey[nCompCount - 1][2] == 10 && nIndArrey[nCompCount - 1][3] == 11) || + ( nIndArrey[nCompCount - 1][0] == 5 && nIndArrey[nCompCount - 1][1] == 7 && + nIndArrey[nCompCount - 1][2] == 8 && nIndArrey[nCompCount - 1][3] == 9 )) { + + VectorField LocVecF[12], LocCompV[12] ; + + for ( int LocInd = 0 ; LocInd < 12 ; ++ LocInd) { + LocVecF[LocInd] = VecField[LocInd] ; + LocCompV[LocInd] = CompoVert[nCompCount - 1][LocInd] ; + } + + if ( ( AreSameVectorApprox( LocVecF[nIndArrey[nCompCount - 1][0]].vtNorm, LocVecF[nIndArrey[nCompCount - 1][1]].vtNorm) && + abs( LocVecF[nIndArrey[nCompCount - 1][0]].vtNorm *LocVecF[nIndArrey[nCompCount - 1][2]].vtNorm) < EPS_SMALL && + abs( LocVecF[nIndArrey[nCompCount - 1][0]].vtNorm * LocVecF[nIndArrey[nCompCount - 1][3]].vtNorm) < EPS_SMALL) || + ( AreSameVectorApprox( LocVecF[nIndArrey[nCompCount - 1][2]].vtNorm, LocVecF[nIndArrey[nCompCount - 1][3]].vtNorm) && + abs( LocVecF[nIndArrey[nCompCount - 1][2]].vtNorm * LocVecF[nIndArrey[nCompCount - 1][0]].vtNorm) < EPS_SMALL && + abs( LocVecF[nIndArrey[nCompCount - 1][2]].vtNorm * LocVecF[nIndArrey[nCompCount - 1][1]].vtNorm) < EPS_SMALL)) { + + Point3d ptBarycenter = ( LocCompV[0].ptInt + LocCompV[1].ptInt + + LocCompV[2].ptInt + LocCompV[3].ptInt) / 4 ; + + if ( abs( LocCompV[0].ptInt.x - ptBarycenter.x) < EPS_SMALL && + abs( LocCompV[1].ptInt.x - ptBarycenter.x) < EPS_SMALL && + abs( LocCompV[2].ptInt.x - ptBarycenter.x) < EPS_SMALL && + abs( LocCompV[3].ptInt.x - ptBarycenter.x) < EPS_SMALL) { + + double dXBar = ptBarycenter.x / m_dStep - 0.5 ; + int nBarLimSup = int( m_nNx[0]) ; + int nBarInd = 0 ; + + while ( nBarInd < nBarLimSup) { + + double dXInd = double( nBarInd) ; + + if ( abs( dXInd - dXBar) < EPS_SMALL) { + + bGridControl = false ; + break ; + } + ++ nBarInd ; + } + } + else if ( abs( LocCompV[0].ptInt.y - ptBarycenter.y) < EPS_SMALL && + abs( LocCompV[1].ptInt.y - ptBarycenter.y) < EPS_SMALL && + abs( LocCompV[2].ptInt.y - ptBarycenter.y) < EPS_SMALL && + abs( LocCompV[3].ptInt.y - ptBarycenter.y) < EPS_SMALL) { + + double dYBar = ptBarycenter.y / m_dStep - 0.5 ; + int nBarLimSup = int( m_nNy[0]) ; + int nBarInd = 0 ; + + while ( nBarInd < nBarLimSup) { + + double dYInd = double( nBarInd) ; + + if ( abs( dYInd - dYBar) < EPS_SMALL) { + + bGridControl = false ; + break ; + } + ++ nBarInd ; + } + } + else if ( abs( LocCompV[0].ptInt.z - ptBarycenter.z) < EPS_SMALL && + abs( LocCompV[1].ptInt.z - ptBarycenter.z) < EPS_SMALL && + abs( LocCompV[2].ptInt.z - ptBarycenter.z) < EPS_SMALL && + abs( LocCompV[3].ptInt.z - ptBarycenter.z) < EPS_SMALL) { + + double dZBar = ptBarycenter.z / m_dStep - 0.5 ; + int nBarLimSup = int( m_nNy[1]) ; + int nBarInd = 0 ; + + while ( nBarInd < nBarLimSup) { + + double dZInd = double( nBarInd) ; + + if ( abs( dZInd - dZBar) < EPS_SMALL) { + + bGridControl = false ; + break ; + } + ++ nBarInd ; + } + } + } + } + else if ( ( nIndArrey[nCompCount - 1][0] == 0 && nIndArrey[nCompCount - 1][1] == 1 && + nIndArrey[nCompCount - 1][2] == 4 && nIndArrey[nCompCount - 1][3] == 5) || + ( nIndArrey[nCompCount - 1][0] == 1 && nIndArrey[nCompCount - 1][1] == 2 && + nIndArrey[nCompCount - 1][2] == 5 && nIndArrey[nCompCount - 1][3] == 6) || + ( nIndArrey[nCompCount - 1][0] == 2 && nIndArrey[nCompCount - 1][1] == 3 && + nIndArrey[nCompCount - 1][2] == 6 && nIndArrey[nCompCount - 1][3] == 7) || + ( nIndArrey[nCompCount - 1][0] == 0 && nIndArrey[nCompCount - 1][1] == 3 && + nIndArrey[nCompCount - 1][2] == 4 && nIndArrey[nCompCount - 1][3] == 7)) { + + VectorField LocVecF[12], LocCompV[12] ; + + for ( int LocInd = 0 ; LocInd < 12 ; ++ LocInd) { + LocVecF[LocInd] = VecField[LocInd] ; + LocCompV[LocInd] = CompoVert[nCompCount - 1][LocInd] ; + } + + if ( ( AreSameVectorApprox( LocVecF[nIndArrey[nCompCount - 1][0]].vtNorm, LocVecF[nIndArrey[nCompCount - 1][2]].vtNorm) && + abs( LocVecF[nIndArrey[nCompCount - 1][0]].vtNorm * LocVecF[nIndArrey[nCompCount - 1][1]].vtNorm) < EPS_SMALL && + abs( LocVecF[nIndArrey[nCompCount - 1][0]].vtNorm * LocVecF[nIndArrey[nCompCount - 1][3]].vtNorm) < EPS_SMALL) || + ( AreSameVectorApprox( LocVecF[nIndArrey[nCompCount - 1][1]].vtNorm, LocVecF[nIndArrey[nCompCount - 1][3]].vtNorm) && + abs( LocVecF[nIndArrey[nCompCount - 1][1]].vtNorm * LocVecF[nIndArrey[nCompCount - 1][0]].vtNorm) < EPS_SMALL && + abs( LocVecF[nIndArrey[nCompCount - 1][1]].vtNorm * LocVecF[nIndArrey[nCompCount - 1][2]].vtNorm) < EPS_SMALL)) { + + Point3d ptBarycenter = ( LocCompV[0].ptInt + LocCompV[1].ptInt + + LocCompV[2].ptInt + LocCompV[3].ptInt) / 4 ; + + if ( abs( LocCompV[0].ptInt.x - ptBarycenter.x) < EPS_SMALL && + abs( LocCompV[1].ptInt.x - ptBarycenter.x) < EPS_SMALL && + abs( LocCompV[2].ptInt.x - ptBarycenter.x) < EPS_SMALL && + abs( LocCompV[3].ptInt.x - ptBarycenter.x) < EPS_SMALL) { + + double dXBar = ptBarycenter.x / m_dStep - 0.5 ; + int nBarLimSup = int( m_nNx[0]) ; + int nBarInd = 0 ; + + while ( nBarInd < nBarLimSup) { + + double dXInd = double( nBarInd) ; + + if ( abs( dXInd - dXBar) < EPS_SMALL) { + + bGridControl = false ; + break ; + } + ++ nBarInd ; + } + } + else if ( abs( LocCompV[0].ptInt.y - ptBarycenter.y) < EPS_SMALL && + abs( LocCompV[1].ptInt.y - ptBarycenter.y) < EPS_SMALL && + abs( LocCompV[2].ptInt.y - ptBarycenter.y) < EPS_SMALL && + abs( LocCompV[3].ptInt.y - ptBarycenter.y) < EPS_SMALL) { + + double dYBar = ptBarycenter.y / m_dStep - 0.5 ; + int nBarLimSup = int( m_nNy[0]) ; + int nBarInd = 0 ; + + while ( nBarInd < nBarLimSup) { + + double dYInd = double( nBarInd) ; + + if ( abs( dYInd - dYBar) < EPS_SMALL) { + + bGridControl = false ; + break ; + } + ++ nBarInd ; + } + } + else if ( abs( LocCompV[0].ptInt.z - ptBarycenter.z) < EPS_SMALL && + abs( LocCompV[1].ptInt.z - ptBarycenter.z) < EPS_SMALL && + abs( LocCompV[2].ptInt.z - ptBarycenter.z) < EPS_SMALL && + abs( LocCompV[3].ptInt.z - ptBarycenter.z) < EPS_SMALL) { + + double dZBar = ptBarycenter.z / m_dStep - 0.5 ; + int nBarLimSup = int( m_nNy[1]) ; + int nBarInd = 0 ; + + while ( nBarInd < nBarLimSup) { + + double dZInd = double( nBarInd) ; + + if ( abs( dZInd - dZBar) < EPS_SMALL) { + + bGridControl = false ; + break ; + } + ++ nBarInd ; + } + } + } + } + } + } + + // Flag ExtMC + bool bExtMC = ( nFeatureType != NO_FEATURE && bGridControl) && ( ! bSecondConfig6) && bEnh ; + + // Extended MC + if ( bExtMC) { + + // Passo al sistema di riferimento del baricentro + Point3d ptGravityCenter( 0, 0, 0) ; + for ( int ni = 0 ; ni < nVertComp[nCompCount - 1] ; ++ ni) + ptGravityCenter = ptGravityCenter + CompoVert[nCompCount - 1][ni].ptInt ; + ptGravityCenter = ptGravityCenter / nVertComp[nCompCount - 1] ; + + Vector3d vtTrasf[12] ; + for ( int ni = 0 ; ni < nVertComp[nCompCount - 1] ; ++ ni) + vtTrasf[ni] = CompoVert[nCompCount - 1][ni].ptInt - ptGravityCenter ; + + // Preparo le matrici per il sistema + typedef Eigen::Matrix dSystemMatrix ; + typedef Eigen::Matrix dSystemVector ; + typedef Eigen::JacobiSVD DecomposerSVD ; + + dSystemMatrix dMatrixN( nVertComp[nCompCount - 1], 3) ; + dSystemVector dKnownVector( nVertComp[nCompCount - 1], 1) ; + + // Caso generale + for ( int ni = 0 ; ni < nVertComp[nCompCount - 1] ; ++ ni) { + dMatrixN( ni, 0) = CompoVert[nCompCount - 1][ni].vtNorm.x ; + dMatrixN( ni, 1) = CompoVert[nCompCount - 1][ni].vtNorm.y ; + dMatrixN( ni, 2) = CompoVert[nCompCount - 1][ni].vtNorm.z ; + dKnownVector( ni) = CompoVert[nCompCount - 1][ni].vtNorm * vtTrasf[ni] ; + } + + // calcolo SVD + DecomposerSVD svd( dMatrixN, Eigen::ComputeThinU | Eigen::ComputeThinV) ; + dSystemMatrix dMatrixV = svd.matrixV() ; + dSystemVector dSingularValue = svd.singularValues() ; + + // Se la feature è un edge annullo il valore singolare minore. + if ( nFeatureType == EDGE) { + double dThres = 0.5 * ( dSingularValue( 1) + dSingularValue( 2)) / dSingularValue( 0) ; + svd.setThreshold( dThres) ; + } + + // risolvo il sistema con SVD, quindi ai minimi quadrati + dSystemVector dUnknownVector( 3, 1) ; + dUnknownVector = svd.solve( dKnownVector) ; + + // Vettore Baricentro-Feature + Vector3d vtFeature( dUnknownVector( 0), dUnknownVector( 1), dUnknownVector( 2)) ; + + // Esprimo la soluzione nel sistema di riferimento z-map + Point3d ptSol = ptGravityCenter + vtFeature ; + + // Limito la feature a una distanza dal baricentro pari alla diagonale del voxel + double dDistFeature = vtFeature.Len() ; + const double MAX_DIST = sqrt( 3) * m_dStep ; + bool bOutside = ( dDistFeature > MAX_DIST) ; + + TRIA3DVECTOR triContainer ; + + // Costruisco triangoli di prova + for ( int ni = 0 ; ni < nVertComp[nCompCount - 1] ; ++ ni) { + int nj = ( ni + 1 < nVertComp[nCompCount - 1]) ? ni + 1 : 0 ; + // Il triangolo è pronto + Triangle3d CurrentTriangle ; + CurrentTriangle.Set( ptSol, CompoVert[nCompCount - 1][nj].ptInt, CompoVert[nCompCount - 1][ni].ptInt) ; + CurrentTriangle.Validate( true) ; + // Aggiungo triangolo al vettore temporaneo + triContainer.emplace_back( CurrentTriangle) ; + } + + // Controllo delle inversioni dei triangoli + bool bDangerInversion = false ; + + // Caso ventaglio con tre vertici di base + if ( nVertComp[nCompCount - 1] == 3) { + + // Controllo se esiste almeno un triangolo con normale avente prodotto scalare + // negativo con la normale di almeno uno dei vertici di base del ventaglio. + bool bInversione = false ; + for ( int ni = 0 ; ni < nVertComp[nCompCount - 1] ; ++ ni) { + + int nj = ( ni + 1 < nVertComp[nCompCount - 1]) ? ni + 1 : 0 ; + + double dDI = triContainer[ni].GetN() * CompoVert[nCompCount - 1][ni].vtNorm ; + double dDJ = triContainer[ni].GetN() * CompoVert[nCompCount - 1][nj].vtNorm ; + + if ( dDI < - EPS_SMALL || dDJ < - EPS_SMALL) { + bInversione = true ; + 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[nCompCount - 1] - 1 ; ++ ni) { + for ( int nj = ni + 1 ; nj < nVertComp[nCompCount - 1] ; ++ nj) { + if ( CompoVert[nCompCount - 1][ni].vtNorm * CompoVert[nCompCount - 1][nj].vtNorm < - EPS_SMALL) + nNegDot ++ ; + } + } + + if ( nNegDot == nVertComp[nCompCount - 1] - 1) { + // Cerco se esiste un punto in cui la normale ha prodotto scalre negativo + // con le normali di entrambi i triangoli che lo contengono + bool bInversione2 = false ; + for ( int ni = 0 ; ni < nVertComp[nCompCount - 1] ; ++ ni) { + int nj = ( ni == 0 ? nVertComp[nCompCount - 1] - 1 : ni - 1) ; + double dDLast = triContainer[ni].GetN() * CompoVert[nCompCount - 1][ni].vtNorm ; + double dDPrev = triContainer[nj].GetN() * CompoVert[nCompCount - 1][ni].vtNorm ; + if ( ( dDLast < EPS_SMALL && dDPrev < EPS_SMALL) || + ( dDLast < - 0.75 || dDPrev < - 0.75)) { + bInversione2 = true ; + break ; + } + } + + // Se tale normale esiste + if ( bInversione2) { + + // Soluzione del sistema nel sistema Zmap + Point3d ptSolZMapFrame = ptSol ; + + // Se la soluzione non cade nel voxel di appartenenza vedo se può + // essere riportata dentro muovendosi lungo la linea di feature. + if ( ! IsPointInsideVoxelApprox( nVoxI, nVoxJ, nVoxK, ptSolZMapFrame)) { + + Vector3d vtNullSpace( dMatrixV( 0, 2), dMatrixV( 1, 2), dMatrixV( 2, 2)) ; + double dParInt1, dParInt2 ; + Point3d ptVoxMin( ( nVoxI + 0.5) * m_dStep, ( nVoxJ + 0.5) * m_dStep, ( nVoxK + 0.5) * m_dStep) ; + Point3d ptVoxMax( ( nVoxI + 1.5) * m_dStep, ( nVoxJ + 1.5) * m_dStep, ( nVoxK + 1.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( ptSolZMapFrame, vtNullSpace, ptVoxMin, ptVoxMax, dParInt1, dParInt2)) { + + triContainer.resize( 0) ; + + double dPar = abs( dParInt1) < abs( dParInt2) ? dParInt1 + ( dParInt2 - dParInt1) / 100 : + dParInt2 + ( dParInt1 - dParInt2) / 100 ; + + Point3d ptNewSol = ptSolZMapFrame + dPar * vtNullSpace ; + + ptSol = ptNewSol ; + + // Costruisco triangoli di prova + for ( int ni = 0 ; ni < nVertComp[nCompCount - 1] ; ++ ni) { + int nj = ( ni + 1 < nVertComp[nCompCount - 1]) ? ni + 1 : 0 ; + // Il triangolo è pronto + Triangle3d CurrentTriangle ; + CurrentTriangle.Set( ptSol, CompoVert[nCompCount - 1][nj].ptInt, CompoVert[nCompCount - 1][ni].ptInt) ; + CurrentTriangle.Validate( true) ; + // Aggiungo triangolo al vettore temporaneo + triContainer.emplace_back( CurrentTriangle) ; + } + } + + // 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( ptSolZMapFrame, 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[nCompCount - 1] == 4) { + + // Controllo se esiste almeno un triangolo con normale avente prodotto scalare + // negativo con la normale di almeno uno dei vertici di base del ventaglio. + bool bInversione = false ; + for ( int ni = 0 ; ni < nVertComp[nCompCount - 1] ; ++ ni) { + + int nj = ( ni + 1 < nVertComp[nCompCount - 1]) ? ni + 1 : 0 ; + + double dDI = triContainer[ni].GetN() * CompoVert[nCompCount - 1][ni].vtNorm ; + double dDJ = triContainer[ni].GetN() * CompoVert[nCompCount - 1][nj].vtNorm ; + + if ( dDI < - EPS_SMALL || dDJ < - EPS_SMALL) + bInversione = true ; + } + // Se tale triangolo esiste continuo i test + if ( bInversione) { + + // Conto il numero di coppie di normali con prodotto scalare negativo + int nNegDot = 0 ; + for ( int ni = 0 ; ni < nVertComp[nCompCount - 1] - 1 ; ++ ni) { + for ( int nj = ni + 1 ; nj < nVertComp[nCompCount - 1] ; ++ nj) { + double dDot = CompoVert[nCompCount - 1][ni].vtNorm * CompoVert[nCompCount - 1][nj].vtNorm ; + if ( dDot < - EPS_SMALL) + nNegDot ++ ; + } + } + + // Caso in cui tale numero è 3 + if ( nNegDot == nVertComp[nCompCount - 1] - 1) { + + Point3d ptSolZMapFrame = ptSol ; + + // Se la feature non cade nel suo voxel + if ( ! IsPointInsideVoxelApprox( nVoxI, nVoxJ, nVoxK, ptSolZMapFrame)) { + + Vector3d vtNullSpace( dMatrixV( 0, 2), dMatrixV( 1, 2), dMatrixV( 2, 2)) ; + + double dParInt1, dParInt2 ; + Point3d ptVoxMin( ( nVoxI + 0.5) * m_dStep, ( nVoxJ + 0.5) * m_dStep, ( nVoxK + 0.5) * m_dStep) ; + Point3d ptVoxMax( ( nVoxI + 1.5) * m_dStep, ( nVoxJ + 1.5) * m_dStep, ( nVoxK + 1.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( ptSolZMapFrame, vtNullSpace, ptVoxMin, ptVoxMax, dParInt1, dParInt2)) { + + triContainer.resize( 0) ; + + double dPar = abs( dParInt1) < abs( dParInt2) ? dParInt1 + ( dParInt2 - dParInt1) / 100 : + dParInt2 + ( dParInt1 - dParInt2) / 100 ; + + Point3d ptNewSol = ptSolZMapFrame + dPar * vtNullSpace ; + + ptSol = ptNewSol ; + + // Costruisco triangoli di prova + for ( int ni = 0 ; ni < nVertComp[nCompCount - 1] ; ++ ni) { + int nj = ( ni + 1 < nVertComp[nCompCount - 1]) ? ni + 1 : 0 ; + // Il triangolo è pronto + Triangle3d CurrentTriangle ; + CurrentTriangle.Set( ptSol, CompoVert[nCompCount - 1][nj].ptInt, CompoVert[nCompCount - 1][ni].ptInt) ; + CurrentTriangle.Validate( true) ; + // Aggiungo triangolo al vettore temporaneo + triContainer.emplace_back( CurrentTriangle) ; + } + } + + // Se non è possibile riportarla dentro e il voxel in + // cui cade è pieno passo alla routine standard + 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 + if ( EdgeTable[nAdjIndex] != 0) + bDangerInversion = true ; + } + } + } + } + // Caso in cui il numero di coppie di normali con prodotto + // scalare negativo non è 3 + else { + + Point3d ptSolZMapFrame = ptSol ; + + if ( ! IsPointInsideVoxelApprox( nVoxI, nVoxJ, nVoxK, ptSolZMapFrame)) { + + Vector3d vtNullSpace( dMatrixV( 0, 2), dMatrixV( 1, 2), dMatrixV( 2, 2)) ; + + double dParInt1, dParInt2 ; + Point3d ptVoxMin( ( nVoxI + 0.5) * m_dStep, ( nVoxJ + 0.5) * m_dStep, ( nVoxK + 0.5) * m_dStep) ; + Point3d ptVoxMax( ( nVoxI + 1.5) * m_dStep, ( nVoxJ + 1.5) * m_dStep, ( nVoxK + 1.5) * m_dStep) ; + + if ( IntersLineBox( ptSolZMapFrame, vtNullSpace, ptVoxMin, ptVoxMax, dParInt1, dParInt2)) { + + triContainer.resize( 0) ; + + double dPar = abs( dParInt1) < abs( dParInt2) ? dParInt1 + ( dParInt2 - dParInt1) / 100: + dParInt2 + ( dParInt1 - dParInt2) / 100 ; + + Point3d ptNewSol = ptSolZMapFrame + dPar * vtNullSpace ; + + ptSol = ptNewSol ; + + // Costruisco triangoli di prova + for ( int ni = 0 ; ni < nVertComp[nCompCount - 1] ; ++ ni) { + int nj = ( ni + 1 < nVertComp[nCompCount - 1]) ? ni + 1 : 0 ; + // Il triangolo è pronto + Triangle3d CurrentTriangle ; + CurrentTriangle.Set( ptSol, CompoVert[nCompCount - 1][nj].ptInt, CompoVert[nCompCount - 1][ni].ptInt) ; + CurrentTriangle.Validate( true) ; + // Aggiungo triangolo al vettore temporaneo + triContainer.emplace_back( CurrentTriangle) ; + } + } + + else { + int nCouple = 0 ; + int nCoupleIndex[4] ; + + // Valuto il numero di coppie di vettori + // quasi coincidenti e per ogni coppia salvo gli + // indici dei vettori + for ( int nNI = 0 ; nNI < 3 ; ++ nNI) { + for ( int nNJ = nNI + 1 ; nNJ < 4 ; ++ nNJ) { + if ( AreSameVectorApprox( CompoVert[nCompCount - 1][nNI].vtNorm, + CompoVert[nCompCount - 1][nNJ].vtNorm)) { + ++ nCouple ; + if ( nCouple == 1) { + nCoupleIndex[0] = nNI ; + nCoupleIndex[1] = nNJ ; + } + else if ( nCouple == 2) { + nCoupleIndex[2] = nNI ; + nCoupleIndex[3] = nNJ ; + } + } + } + } + + // caso due coppie + if ( nCouple == 2) { + // vedo se c'è un triangolo con normale invertita rispetto a quelle + // si entrambi i vertici, se esiste si passerà a std MC + for ( int ni = 0 ; ni < 4 ; ++ ni) { + int nj = ( ni == 3 ? 0 : ni + 1) ; + if ( triContainer[ni].GetN() * CompoVert[nCompCount - 1][ni].vtNorm < - 0.95 && + triContainer[ni].GetN() * CompoVert[nCompCount - 1][nj].vtNorm < - 0.95) { + bDangerInversion = true ; + break ; + } + } + } + // caso una coppia + else if ( nCouple == 1) { + // cerco gli indici dei vettori non appartenenti alla coppia + for ( int ni = 0 ; ni < 4 ; ++ ni) { + if ( ni != nCoupleIndex[0] && ni != nCoupleIndex[1]) { + nCoupleIndex[2] = ni ; + break ; + } + } + for ( int ni = 0 ; ni < 4 ; ++ ni) { + if ( ni != nCoupleIndex[0] && ni != nCoupleIndex[1] && ni != nCoupleIndex[2]) { + nCoupleIndex[3] = ni ; + break ; + } + } + // Media dei vettori coppia + Vector3d vtAv01 = 0.5 * ( CompoVert[nCompCount - 1][nCoupleIndex[0]].vtNorm + + CompoVert[nCompCount - 1][nCoupleIndex[1]].vtNorm) ; + // vettore nello spazio genenrato dai due non appartenenti alla coppia + Vector3d vtAv23 = 0.5 * ( CompoVert[nCompCount - 1][nCoupleIndex[2]].vtNorm + + CompoVert[nCompCount - 1][nCoupleIndex[3]].vtNorm) ; + + vtAv01.Normalize() ; + vtAv23.Normalize() ; + double dDAvAV = vtAv01 * vtAv23 ; + // se angolo grande si esegue std MC + if ( abs( vtAv01 * vtAv23) < EPS_SMALL) { + for ( int ni = 0 ; ni < 4 ; ++ ni) { + int nj = ni == 3 ? 0 : ni + 1 ; + if ( triContainer[ni].GetN() * CompoVert[nCompCount - 1][ni].vtNorm < - 0.95 && + triContainer[ni].GetN() * CompoVert[nCompCount - 1][nj].vtNorm < - 0.95) { + bDangerInversion = true ; + break ; + } + } + } + else { + + double dD23 = CompoVert[nCompCount - 1][nCoupleIndex[2]].vtNorm * + CompoVert[nCompCount - 1][nCoupleIndex[3]].vtNorm ; + if ( dD23 > 0.7 && dDAvAV < 0.7) { + + for ( int ni = 0 ; ni < 4 ; ++ ni) { + int nj = ni == 3 ? 0 : ni + 1 ; + if ( triContainer[ni].GetN() * CompoVert[nCompCount - 1][ni].vtNorm < - 0.9 && + triContainer[ni].GetN() * CompoVert[nCompCount - 1][nj].vtNorm < - 0.9) { + bDangerInversion = true ; + break ; + } + } + } + } + } + } + } + } + } + } + + // Valuto normali: questo è ancora un controllo + // sulle normali, se risultano in tutti i punti + // approssimativamente uguali passiamo alla + // routine standard + int nContSize = int( triContainer.size()) ; + + bool bPlane = true ; + + for ( int ni = 0 ; ni < nContSize - 1 ; ++ ni) { + for ( int nj = ni + 1 ; nj < nContSize ; ++ nj) { + + Vector3d vtI = triContainer[ni].GetN() ; + Vector3d vtJ = triContainer[nj].GetN() ; + + if ( ! AreSameVectorApprox( vtI, vtJ)) { + bPlane = false ; + break ; + } + } + + if ( ! bPlane) + break ; + } + + // Se feature nei limiti e triangoli non in un piano, confermo ExtMC + if ( ! bOutside && ! bPlane && ! bDangerInversion) { + + TRIA3DVECTOR vInnerTriaTemp ; + + 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] = { nVoxI, nVoxJ, nVoxK} ; + Point3d ptTemp = ptSol ; + Triangle3d trTemp = triContainer[ni] ; + vInnerTriaTemp.emplace_back( triContainer[ni]) ; + } + } + + 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 = nVoxI ; + triHold[nCurrent].j = nVoxJ ; + triHold[nCurrent].k = nVoxK ; + } + + // Indice che identifica la posizione del voxel nel vector. + int nCurrent = int( triHold.size()) - 1 ; + + // Aggiungo vertice della componente connessa alla lista dei vertici. + triHold[nCurrent].ptCompoVert.emplace_back( ptSol) ; + + int nOldFeatureNum = int( triHold[nCurrent].vCompoTria.size()) ; + int nNewFeatureNum = nOldFeatureNum + 1 ; + + // Aggiungo una componente per il vettore dei triangoli della componente connessa. + triHold[nCurrent].vCompoTria.resize( nNewFeatureNum) ; + for ( int ni = 0 ; ni < int( vInnerTriaTemp.size()) ; ++ ni) { + triHold[nCurrent].vCompoTria[nOldFeatureNum].emplace_back( vInnerTriaTemp[ni]) ; + } + } + } + + // ExtMC non confermato, si passa a MC + else { + // Costruzione dei triangoli + for ( int TriIndex = 0; TriIndex < ( nVertComp[nCompCount - 1] - 2) * 3 ; TriIndex += 3) { + // Il triangolo è pronto + Triangle3d CurrentTriangle ; + CurrentTriangle.Set( CompoTriVert[nCompCount - 1][TriIndex].ptInt, + CompoTriVert[nCompCount - 1][TriIndex+1].ptInt, + CompoTriVert[nCompCount - 1][TriIndex+2].ptInt) ; + bool bV = CurrentTriangle.Validate( true) ; + // Aggiungo alla lista + lstTria.emplace_back( CurrentTriangle) ; + } + } + } + + // Standard MC + else { + // Costruzione dei triangoli + for ( int TriIndex = 0; TriIndex < ( nVertComp[nCompCount - 1] - 2) * 3 ; TriIndex += 3) { + // Il triangolo è pronto + Triangle3d CurrentTriangle ; + CurrentTriangle.Set( CompoTriVert[nCompCount - 1][TriIndex].ptInt, + CompoTriVert[nCompCount - 1][TriIndex+1].ptInt, + CompoTriVert[nCompCount - 1][TriIndex+2].ptInt) ; + bool bV = CurrentTriangle.Validate( true) ; + // Aggiungo alla lista + lstTria.emplace_back( CurrentTriangle) ; + } + } } } - return true ; } @@ -1301,28 +2033,23 @@ VolZmap::ExtMarchingCubes( int nBlock, TRIA3DLIST& lstTria, TriHolder& triHold) 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 + // più o meno concordi (soglia su prodotto scalare 0, forse 0.7?) + Vector3d vtCmpAvg0, vtCmpAvg1 ; bool bTest0 = DotTest( CompoVert[0], 4, vtCmpAvg0) ; bool bTest1 = DotTest( CompoVert[1], 4, vtCmpAvg1) ; - if ( ! ( bTest0 && bTest1)) { - + if ( ! bTest0 || ! bTest1) { int nt = 0 ; while ( nIndexVsIndex10[nt][0] != nIndex) ++ nt ; - int nRotCase = nIndexVsIndex10[nt][1] ; - // 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] ; @@ -1418,17 +2145,15 @@ VolZmap::ExtMarchingCubes( int nBlock, TRIA3DLIST& lstTria, TriHolder& triHold) int nBarInd = 0 ; while ( nBarInd < nBarLimSup) { - double dXInd = double( nBarInd) ; - if ( abs( dXInd - dXBar) < EPS_SMALL) { - bGridControl = false ; break ; } ++ nBarInd ; } } + else if ( abs( LocCompV[0].ptInt.y - ptBarycenter.y) < EPS_SMALL && abs( LocCompV[1].ptInt.y - ptBarycenter.y) < EPS_SMALL && abs( LocCompV[2].ptInt.y - ptBarycenter.y) < EPS_SMALL && @@ -1439,17 +2164,15 @@ VolZmap::ExtMarchingCubes( int nBlock, TRIA3DLIST& lstTria, TriHolder& triHold) int nBarInd = 0 ; while ( nBarInd < nBarLimSup) { - double dYInd = double( nBarInd) ; - if ( abs( dYInd - dYBar) < EPS_SMALL) { - bGridControl = false ; break ; } ++ nBarInd ; } } + else if ( abs( LocCompV[0].ptInt.z - ptBarycenter.z) < EPS_SMALL && abs( LocCompV[1].ptInt.z - ptBarycenter.z) < EPS_SMALL && abs( LocCompV[2].ptInt.z - ptBarycenter.z) < EPS_SMALL && @@ -1460,11 +2183,8 @@ VolZmap::ExtMarchingCubes( int nBlock, TRIA3DLIST& lstTria, TriHolder& triHold) int nBarInd = 0 ; while ( nBarInd < nBarLimSup) { - double dZInd = double( nBarInd) ; - if ( abs( dZInd - dZBar) < EPS_SMALL) { - bGridControl = false ; break ; } @@ -1473,6 +2193,7 @@ VolZmap::ExtMarchingCubes( int nBlock, TRIA3DLIST& lstTria, TriHolder& triHold) } } } + else if ( ( nIndArrey[nCompCount - 1][0] == 0 && nIndArrey[nCompCount - 1][1] == 1 && nIndArrey[nCompCount - 1][2] == 4 && nIndArrey[nCompCount - 1][3] == 5) || ( nIndArrey[nCompCount - 1][0] == 1 && nIndArrey[nCompCount - 1][1] == 2 && @@ -1509,17 +2230,15 @@ VolZmap::ExtMarchingCubes( int nBlock, TRIA3DLIST& lstTria, TriHolder& triHold) int nBarInd = 0 ; while ( nBarInd < nBarLimSup) { - double dXInd = double( nBarInd) ; - if ( abs( dXInd - dXBar) < EPS_SMALL) { - bGridControl = false ; break ; } ++ nBarInd ; } } + else if ( abs( LocCompV[0].ptInt.y - ptBarycenter.y) < EPS_SMALL && abs( LocCompV[1].ptInt.y - ptBarycenter.y) < EPS_SMALL && abs( LocCompV[2].ptInt.y - ptBarycenter.y) < EPS_SMALL && @@ -1530,17 +2249,15 @@ VolZmap::ExtMarchingCubes( int nBlock, TRIA3DLIST& lstTria, TriHolder& triHold) int nBarInd = 0 ; while ( nBarInd < nBarLimSup) { - double dYInd = double( nBarInd) ; - if ( abs( dYInd - dYBar) < EPS_SMALL) { - bGridControl = false ; break ; } ++ nBarInd ; } } + else if ( abs( LocCompV[0].ptInt.z - ptBarycenter.z) < EPS_SMALL && abs( LocCompV[1].ptInt.z - ptBarycenter.z) < EPS_SMALL && abs( LocCompV[2].ptInt.z - ptBarycenter.z) < EPS_SMALL && @@ -1551,11 +2268,8 @@ VolZmap::ExtMarchingCubes( int nBlock, TRIA3DLIST& lstTria, TriHolder& triHold) int nBarInd = 0 ; while ( nBarInd < nBarLimSup) { - double dZInd = double( nBarInd) ; - if ( abs( dZInd - dZBar) < EPS_SMALL) { - bGridControl = false ; break ; } @@ -1591,73 +2305,6 @@ VolZmap::ExtMarchingCubes( int nBlock, TRIA3DLIST& lstTria, TriHolder& triHold) dSystemMatrix dMatrixN( nVertComp[nCompCount - 1], 3) ; dSystemVector dKnownVector( nVertComp[nCompCount - 1], 1) ; -#if 0 - // Studio del caso 4 punti su un piano - int nEqual = 0 ; - int nPosD ; - Vector3d vtD, vtE ; - - if ( nVertComp[nCompCount - 1] == 4 && nFeatureType == EDGE) { - int nPosEq ; - for ( int ni = 0 ; ni < 2 ; ++ ni) { - for ( int nj = ni + 1 ; nj < nVertComp[nCompCount - 1] ; ++ nj) { - - if ( AreSameVectorApprox( CompoVert[nCompCount - 1][ni].vtNorm, - CompoVert[nCompCount - 1][nj].vtNorm)) { - nEqual ++ ; - nPosEq = ni ; - } - } - if ( nEqual == 2) - break ; - else - nEqual = 0 ; - } - if ( nEqual == 2) { - - for ( int ni = 0 ; ni < nVertComp[nCompCount - 1] ; ++ ni) - - if ( ! AreSameVectorApprox( CompoVert[nCompCount - 1][ni].vtNorm, - CompoVert[nCompCount - 1][nPosEq].vtNorm)) { - nPosD = ni ; - vtD = CompoVert[nCompCount - 1][ni].vtNorm ; - vtE = CompoVert[nCompCount - 1][nPosEq].vtNorm ; - } - } - } - - double dDot = abs( ( CompoVert[nCompCount - 1][1].ptInt - CompoVert[nCompCount - 1][0].ptInt) * - ( ( CompoVert[nCompCount - 1][2].ptInt - CompoVert[nCompCount - 1][1].ptInt) ^ - ( CompoVert[nCompCount - 1][3].ptInt - CompoVert[nCompCount - 1][2].ptInt))) ; - - // Caso superficie piana - if ( false && nVertComp[nCompCount - 1] == 4 && nEqual == 2 && dDot < EPS_SMALL) { - - for ( int ni = 0 ; ni < nVertComp[nCompCount - 1] ; ++ ni) { - if ( ni != nPosD) { - dMatrixN( ni, 0) = CompoVert[nCompCount - 1][ni].vtNorm.x ; - dMatrixN( ni, 1) = CompoVert[nCompCount - 1][ni].vtNorm.y ; - dMatrixN( ni, 2) = CompoVert[nCompCount - 1][ni].vtNorm.z ; - dKnownVector( ni) = CompoVert[nCompCount - 1][ni].vtNorm * vtTrasf[ni] ; - } - else { - dMatrixN( ni, 0) = vtE.x ; - dMatrixN( ni, 1) = vtE.y ; - dMatrixN( ni, 2) = vtE.z ; - dKnownVector( ni) = vtE *vtTrasf[ni] ; - } - } - } - // Caso generale - else { - for ( int ni = 0 ; ni < nVertComp[nCompCount - 1] ; ++ ni) { - dMatrixN( ni, 0) = CompoVert[nCompCount - 1][ni].vtNorm.x ; - dMatrixN( ni, 1) = CompoVert[nCompCount - 1][ni].vtNorm.y ; - dMatrixN( ni, 2) = CompoVert[nCompCount - 1][ni].vtNorm.z ; - dKnownVector( ni) = CompoVert[nCompCount - 1][ni].vtNorm * vtTrasf[ni] ; - } - } -#else // Caso generale for ( int ni = 0 ; ni < nVertComp[nCompCount - 1] ; ++ ni) { dMatrixN( ni, 0) = CompoVert[nCompCount - 1][ni].vtNorm.x ; @@ -1665,7 +2312,6 @@ VolZmap::ExtMarchingCubes( int nBlock, TRIA3DLIST& lstTria, TriHolder& triHold) dMatrixN( ni, 2) = CompoVert[nCompCount - 1][ni].vtNorm.z ; dKnownVector( ni) = CompoVert[nCompCount - 1][ni].vtNorm * vtTrasf[ni] ; } -#endif // calcolo SVD DecomposerSVD svd( dMatrixN, Eigen::ComputeThinU | Eigen::ComputeThinV) ; @@ -2030,34 +2676,25 @@ VolZmap::ExtMarchingCubes( int nBlock, TRIA3DLIST& lstTria, TriHolder& triHold) } } - // Valuto normali: questo è ancora un controllo - // sulle normali, se risultano in tutti i punti - // approssimativamente uguali passiamo alla - // routine standard + // Controllo sulle normali, se risultano in tutti i punti + // approssimativamente uguali passiamo alla routine standard int nContSize = int( triContainer.size()) ; - bool bPlane = true ; - for ( int ni = 0 ; ni < nContSize - 1 ; ++ ni) { for ( int nj = ni + 1 ; nj < nContSize ; ++ nj) { - Vector3d vtI = triContainer[ni].GetN() ; Vector3d vtJ = triContainer[nj].GetN() ; - if ( ! AreSameVectorApprox( vtI, vtJ)) { bPlane = false ; break ; } } - if ( ! bPlane) break ; } - // Se la feature non è fuori dai limiti - // e i triangoli formano una superficie - // non piana confermo ExtMC - if ( ! ( bOutside || bPlane) && ! bDangerInversion) { + // 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} ; @@ -2069,13 +2706,9 @@ VolZmap::ExtMarchingCubes( int nBlock, TRIA3DLIST& lstTria, TriHolder& triHold) // 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 ; - Triangle3d trTemp = triContainer[ni] ; - if ( ! bBoundary || ! IsATriangleOnBorder( trTemp, ptTemp, nLimits, nVoxIJK)) vInnerTriaTemp.emplace_back( triContainer[ni]) ; else @@ -2098,14 +2731,9 @@ VolZmap::ExtMarchingCubes( int nBlock, TRIA3DLIST& lstTria, TriHolder& triHold) // voxel, ridimensiono il vettore che contiene // la struttura dati del voxel. if ( nInnerFeatureInVoxel == 1) { - triHold.resize( triHold.size() + 1) ; - // Questi dati dipendono solo dal voxel, - // quindi sono validi per tutte le - // componenti che vi appartengono. int nCurrent = int( triHold.size()) - 1 ; - triHold[nCurrent].i = i ; triHold[nCurrent].j = j ; triHold[nCurrent].k = k ; @@ -2141,14 +2769,9 @@ VolZmap::ExtMarchingCubes( int nBlock, TRIA3DLIST& lstTria, TriHolder& triHold) // voxel, ridimensiono il vettore che contiene // la struttura dati del voxel. if ( nBorderFeatureInVoxel == 1) { - m_InterBlockTria[nBlock].resize( m_InterBlockTria[nBlock].size() + 1) ; - - // Questi dati dipendono solo dal voxel, - // quindi sono validi per tutte le - // componenti che vi appartengono. + // 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 ; @@ -2224,8 +2847,84 @@ VolZmap::ExtMarchingCubes( int nBlock, TRIA3DLIST& lstTria, TriHolder& triHold) } //---------------------------------------------------------------------------- +// 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::FlipEdgesII( TriHolder& TriHold) const +VolZmap::ExtMarchingCubes( vector& vVox, TRIA3DLIST& lstTria, bool bEnh) const +{ + // Contenitore dei triangoli costituenti sharp-feaure + TriHolder triHold ; + int nOriginalSize = int( vVox.size()) ; + for ( int nV = 0 ; nV < int( vVox.size()) ; ++ nV) { + // I voxel intersecati dalla retta hanno indice nel vettore + // nV < OriginalSize; quelli adiacenti avranno indice + // nV >= OriginalSize + if ( nV < nOriginalSize) { + // Se a questa iterazione troviamo una sharp-feature + // la dimensione del TriHolder aumenta + int nCurrentSize = int( triHold.size()) ; + ProcessCube( vVox[nV].nI, vVox[nV].nJ, vVox[nV].nK, lstTria, triHold, bEnh) ; + int nNewSize = int( triHold.size()) ; + + // Se abbiamo trovato una sharp-feature + if ( nNewSize > nCurrentSize) { + // Ciclo fra tutti i possibili voxel adiacenti + for ( int i = - 1 ; i <= 1 ; ++ i) { + for ( int j = - 1 ; j <= 1 ; ++ j) { + for ( int k = - 1 ; k <= 1 ; ++ k) { + if ( i == 0 && j == 0 && k == 0) + continue ; + if ( IsValidVoxel( vVox[nV].nI + i, vVox[nV].nJ + j, vVox[nV].nK + k)) { + // Se il voxel adiacente è già nella lista di quelli attraversati dalla retta + // lo ignoro, altrimenti lo metto in coda nel TriHolder. + int nOldV = 0 ; + for ( ; nOldV < nOriginalSize ; ++ nOldV) { + + if ( vVox[nOldV].nI == vVox[nV].nI + i || + vVox[nOldV].nJ == vVox[nV].nJ + j || + vVox[nOldV].nK == vVox[nV].nK + k) + continue ; + } + if ( nOldV == nOriginalSize) { + Voxel NewVox ; + NewVox.nI = vVox[nV].nI + i ; + NewVox.nJ = vVox[nV].nJ + j ; + NewVox.nK = vVox[nV].nK + k ; + vVox.emplace_back( NewVox) ; + } + } + } + } + } + } + } + // Se il voxel non è attraversato dalla retta, ma fa parte di + // quelli adiacenti lo processiamo normalmente + else + ProcessCube( vVox[nV].nI, vVox[nV].nJ, vVox[nV].nK, lstTria, triHold, bEnh) ; + } + // Eseguo il flipping + FlipEdgesII( triHold, false) ; + // Aggiungo alla lista di triangoli, tutti quelli che formano una sharp-feature. + for ( int nVox = 0 ; nVox < int( triHold.size()) ; ++ nVox) { + for ( int nCompo = 0 ; nCompo < int( triHold[nVox].vCompoTria.size()) ; ++ nCompo) { + for ( int nTri = 0 ; nTri < int( triHold[nVox].vCompoTria[nCompo].size()) ; ++ nTri) { + lstTria.emplace_back( triHold[nVox].vCompoTria[nCompo][nTri]) ; + } + } + } + + return true ; +} + +//---------------------------------------------------------------------------- +// Esegue il flipping dei triangoli contenuti nel TriHolder. +// Se si usa la funzione per la visualizzazione dei triangoli, +// bGraph deve valere true, se si intende calcolare la profondità +// del materiale bGraph deve valere false. +bool +VolZmap::FlipEdgesII( TriHolder& TriHold, bool bGraph) const { // Numero di voxel in cui si presentano sharp feature int nVoxelNum = int( TriHold.size()) ; @@ -2235,6 +2934,24 @@ VolZmap::FlipEdgesII( TriHolder& TriHold) const for ( int n2 = n1 ; n2 < nVoxelNum ; ++ n2) { + // Flag voxel adiacenti o meno + bool bAdj ; + // Se la funzione è chiamata per la grafica i voxel sono ordinati nel TriHolder, + // non serve confrontare i triangoli con quelli dei voxel precedenti. + if ( bGraph) { + bAdj = ( TriHold[n2].i >= TriHold[n1].i && TriHold[n2].i <= TriHold[n1].i + 1) || + ( TriHold[n2].j >= TriHold[n1].j && TriHold[n2].j <= TriHold[n1].j + 1) || + ( TriHold[n2].k >= TriHold[n1].k && TriHold[n2].k <= TriHold[n1].k + 1) ; + } + // Se la funzione è chiamata per il calcolo della profondità del materiale i voxel non + // sono ordinati nel TriHolder, quindi può capitare di dover confrontare i triangoli con + // quelli dei voxel precedenti. + else { + bAdj = abs( TriHold[n2].i - TriHold[n1].i) == 1 || + abs( TriHold[n2].j - TriHold[n1].j) == 1 || + abs( TriHold[n2].k - TriHold[n1].k) == 1 ; + } + // Se i voxel sono adiacenti proseguo if ( ( TriHold[n2].i >= TriHold[n1].i && TriHold[n2].i <= TriHold[n1].i + 1) || ( TriHold[n2].j >= TriHold[n1].j && TriHold[n2].j <= TriHold[n1].j + 1) || @@ -2325,6 +3042,108 @@ VolZmap::FlipEdgesII( TriHolder& TriHold) const return true ; } +//////// SPERIMENTALE /////////////////// +//bool +//VolZmap::FlipEdgesII( TriHolder& TriHold) const +//{ +// // Numero di voxel in cui si presentano sharp feature +// int nVoxelNum = int( TriHold.size()) ; +// +// // Ciclo su tali voxel +// for ( int n1 = 0 ; n1 < nVoxelNum ; ++ n1) { +// +// for ( int n2 = n1 ; n2 < nVoxelNum ; ++ n2) { +// +// // Se i voxel sono adiacenti proseguo +// if ( abs( TriHold[n2].i - TriHold[n1].i) == 1 || +// abs( TriHold[n2].j - TriHold[n1].j) == 1 || +// abs( TriHold[n2].k - TriHold[n1].k) == 1) { +// +// // Numero delle componenti connesse nei due voxel +// int nNumCompo1 = int( TriHold[n1].ptCompoVert.size()) ; +// int nNumCompo2 = int( TriHold[n2].ptCompoVert.size()) ; +// +// int nCompo1 = 0 ; +// +// // Ciclo sulle componenti +// for ( ; nCompo1 < nNumCompo1 ; ++ nCompo1) { +// +// int nCompo2 = ( n1 == n2 ? nCompo1 + 1 : 0) ; +// +// for ( ; nCompo2 < nNumCompo2 ; ++ nCompo2) { +// +// // Numero di triangoli per le componenti connesse +// int nTriNum1 = int( TriHold[n1].vCompoTria[nCompo1].size()) ; +// int nTriNum2 = int( TriHold[n2].vCompoTria[nCompo2].size()) ; +// +// for ( int nTri1 = 0 ; nTri1 < nTriNum1 ; ++ nTri1) { +// +// bool bModified = false ; +// +// for ( int nTri2 = 0 ; nTri2 < nTriNum2 ; ++ nTri2) { +// +// INTVECTOR SharedIndex ; +// +// for ( int nVert1 = 0 ; nVert1 < 3 ; ++ nVert1) { +// +// for ( int nVert2 = 0 ; nVert2 < 3 ; ++ nVert2) { +// +// Point3d ptP1 = TriHold[n1].vCompoTria[nCompo1][nTri1].GetP( nVert1) ; +// Point3d ptP2 = TriHold[n2].vCompoTria[nCompo2][nTri2].GetP( nVert2) ; +// +// if ( AreSamePointEpsilon( ptP1, ptP2, EPS_ZERO)) { +// +// Point3d ptVert1 = TriHold[n1].ptCompoVert[nCompo1] ; +// Point3d ptVert2 = TriHold[n2].ptCompoVert[nCompo2] ; +// +// if ( ! ( AreSamePointEpsilon( ptP1, ptVert1, EPS_ZERO) || +// AreSamePointEpsilon( ptP2, ptVert2, EPS_ZERO))) { +// +// SharedIndex.emplace_back( nVert1) ; +// SharedIndex.emplace_back( nVert2) ; +// } +// } +// +// if ( SharedIndex.size() > 2) +// break ; +// } +// +// if ( SharedIndex.size() > 2) +// break ; +// } +// +// // Si deve operare la modifica dei triangoli +// if ( SharedIndex.size() > 2) { +// +// // verifico che i due lati adiacenti siano controversi +// int nProd = ( SharedIndex[2] - SharedIndex[0]) * ( SharedIndex[3] - SharedIndex[1]) ; +// if ( nProd != 1 && nProd != - 2 && nProd != 4 && nProd != 0) { +// +// TriHold[n1].vCompoTria[nCompo1][nTri1].SetP( SharedIndex[0], +// TriHold[n2].ptCompoVert[nCompo2]) ; +// TriHold[n2].vCompoTria[nCompo2][nTri2].SetP( SharedIndex[3], +// TriHold[n1].ptCompoVert[nCompo1]) ; +// +// TriHold[n1].vCompoTria[nCompo1][nTri1].Validate( true) ; +// TriHold[n2].vCompoTria[nCompo2][nTri2].Validate( true) ; +// +// bModified = true ; +// break ; +// } +// } +// } +// +// if ( bModified) +// break ; +// } +// } +// } +// } +// } +// } +// +// return true ; +//} //---------------------------------------------------------------------------- bool diff --git a/VolZmap.h b/VolZmap.h index ac944e1..94a2b8e 100644 --- a/VolZmap.h +++ b/VolZmap.h @@ -96,8 +96,9 @@ class VolZmap : public IVolZmap, public IGeoObjRW bool SetChiselTool( const std::string& sToolName, double dH, double dW, double dTh) override ; bool MillingStep( const Point3d& ptPs, const Vector3d& vtDs, const Point3d& ptPe, const Vector3d& vtDe) override ; bool MillingStep( const Point3d& ptPs, const Vector3d& vtDs, const Vector3d& vtAs, const Point3d& ptPe, const Vector3d& vtDe, const Vector3d& vtAe) override ; - bool GetDepth( const Point3d& ptP, const Vector3d& vtDir, double& dInLength, double& dOutLength) const override ; - bool GetDepthWithVoxel( const Point3d& ptPLoc, const Vector3d& vtDir, double& dInLength, double& dOutLength) const ; + bool GetDepth( const Point3d& ptP, const Vector3d& vtD, double& dInLength, double& dOutLength, bool bExact) const override ; + bool GetDepthWithDexel( const Point3d& ptP, const Vector3d& vtDir, double& dInLength, double& dOutLength) const ; + bool GetDepthWithVoxel( const Point3d& ptP, const Vector3d& vtDir, double& dInLength, double& dOutLength, bool bEnh) const ; bool AvoidBox( const Frame3d& frBox, const Vector3d& vtDiag) const override ; VolZmap* ClonePart( int nPart) const override ; bool RemovePart( int nPart) override ; @@ -129,6 +130,11 @@ class VolZmap : public IVolZmap, public IGeoObjRW VOX_ON_BOUNDARY = 0, VOX_INNER = -1} ; typedef std::unordered_map VoxelContainer ; + + // Struttura voxel + struct Voxel { + int nI, nJ, nK ; + } ; private : bool CopyFrom( const VolZmap& clSrc) ; @@ -139,9 +145,11 @@ class VolZmap : public IVolZmap, public IGeoObjRW bool AddDexelSideFace( int nPos, int nPosAdj, const Point3d& ptP, const Point3d& ptQ, const Vector3d& vtZ, const Vector3d& vtNorm, TRIA3DLIST& lstTria) const ; bool ProcessCube( int nVoxI, int nVoxJ, int nVoxK, int& nCubeIndex, TRIA3DLIST& lstTria) const ; + bool ProcessCube( int nVoxI, int nVoxJ, int nVoxK, TRIA3DLIST& lstTria, TriHolder& triHold, bool bExt) const ; ///// ora bool MarchingCubes( int nBlock, TRIA3DLIST& lstTria) const ; bool ExtMarchingCubes( int nBlock, TRIA3DLIST& lstTria, TriHolder& triHold) const ; - bool FlipEdgesII( TriHolder& TriHold) const ; + bool ExtMarchingCubes( std::vector& vVox, TRIA3DLIST& lstTria, bool bEnh) const ; + bool FlipEdgesII( TriHolder& TriHold, bool bGraph) const ; bool FlipEdgesBB( TriaMatrix& InterTria) const ; bool IsThereMat( int nI, int nJ, int nK) const ; int CalcIndex( int nI, int nJ, int nK) const ; @@ -239,6 +247,7 @@ class VolZmap : public IVolZmap, public IGeoObjRW bool IntersLineBox( const Point3d& ptP, const Vector3d& vtV, const Point3d& ptMin, const Point3d& ptMax) const ; bool IntersLineBox( const Point3d& ptP, const Vector3d& vtV, const Point3d& ptMin, const Point3d& ptMax, double& dU1, double& dU2) const ; + bool IntersLineZMapLattice( const Point3d& ptP, const Vector3d& vtV, double& dU1, double& dU2) const ; bool IntersLineZMapBBox( const Point3d& ptP, const Vector3d& vtV, unsigned int nGrid, double& dU1, double& dU2) const ; bool IntersLineDexel( const Point3d& ptP, const Vector3d& vtV, unsigned int nGrid, unsigned int nI, unsigned int nJ, double& dU1, double& dU2) const ; diff --git a/VolZmapTool.cpp b/VolZmapTool.cpp index 9dfcbf8..b70ecd7 100644 --- a/VolZmapTool.cpp +++ b/VolZmapTool.cpp @@ -317,7 +317,7 @@ VolZmap::SetMortiserTool( const std::string& sToolName, double dH, double dW, do //---------------------------------------------------------------------------- bool -VolZmap::SetChiselTool( const std::string& sToolName, double dH, double dW, double dTh) { +VolZmap::SetChiselTool( const string& sToolName, double dH, double dW, double dTh) { // Reset nome m_sToolName.clear() ;