//----------------------------------------------------------------------------
//  EgalTech 2015-2021
//----------------------------------------------------------------------------
// File : VolZmap.cpp         Data : 13.09.21           Versione : 2.3i1
// Contenuto : Implementazione della classe Volume Zmap (tre griglie)
//
//
//
// Modifiche : 22.01.15 DS Creazione modulo.
//             13.09.21 LM Correzione funzione che determina se voxel su frontiera di un blocco.
//
//----------------------------------------------------------------------------

//--------------------------- Include ----------------------------------------
#include "stdafx.h"
#include "CurveLine.h"
#include "VolZmap.h"
#include "GeoConst.h"
#include "MC_Tables.h"
#include "PolygonPlane.h"
#include "IntersLineBox.h"
#include "SurfTriMesh.h"
#include "/EgtDev/Include/EGkIntervals.h"
#include "/EgtDev/Include/EGkStringUtils3d.h"
#include "/EgtDev/Include/EGkChainCurves.h"
#include "/EgtDev/Include/EgtNumUtils.h"
#define EIGEN_NO_IO
#include "/EgtDev/Extern/Eigen/Core"
#include "/EgtDev/Extern/Eigen/SVD"
#include <thread>
#include <future>

using namespace std ;

// ------------------------- Tipi per SVD con Eigen ------------------------------------------------------------------------------
static const int MAX_FAN_BASE_VERTS = 7 ; 
typedef Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor, MAX_FAN_BASE_VERTS, 3> SvdMatrix ;
typedef Eigen::Matrix<double, Eigen::Dynamic, 1, Eigen::ColMajor, MAX_FAN_BASE_VERTS, 1> SvdVector ;
typedef Eigen::JacobiSVD<SvdMatrix, Eigen::QRPreconditioners::ColPivHouseholderQRPreconditioner> SvdDecomposer ;

// ------------------------- FUNZIONE TEST SULLE NORMALI --------------------------------------------------------------------------
enum FatureType { NO_FEATURE = 0, CORNER = 1, EDGE = 2} ; 
enum CanonicDir { X_PLUS = 1, X_MINUS = -1, Y_PLUS = 2, Y_MINUS = -2, Z_PLUS = 3, Z_MINUS = -3} ;

//----------------------------------------------------------------------------
bool
Config2VertOrder( int nInd)
{
   if ( nInd == 111  ||  nInd == 119  ||  nInd == 159  ||  nInd == 187  ||  nInd == 207  ||  nInd == 221  ||
        nInd == 238  ||  nInd == 243  ||  nInd == 246  ||  nInd == 249  ||  nInd == 205)
      return true ;
   return false ;
}

//----------------------------------------------------------------------------
bool
Config3Duality(int nInd)
{
   if (nInd == 5 || nInd == 10 || nInd == 18 || nInd == 24 || nInd == 33 || nInd == 36 ||
      nInd == 66 || nInd == 72 || nInd == 80 || nInd == 129 || nInd == 132 || nInd == 160)
      return true;
   return false;
}


//----------------------------------------------------------------------------
bool
Config6Duality( int nInd)
{
   if ( nInd == 21  ||  nInd == 22  ||  nInd == 28  ||  nInd == 41  ||  nInd == 42  ||  nInd == 44  ||
        nInd == 52  ||  nInd == 56  ||  nInd == 67  ||  nInd == 69  ||  nInd == 73  ||  nInd == 81  ||  
        nInd == 84  ||  nInd == 97  ||  nInd == 104  ||  nInd == 131  ||  nInd == 134  ||  nInd == 138  || 
        nInd == 146  ||  nInd == 148  ||  nInd == 162  ||  nInd == 168  ||  nInd == 193  ||  nInd == 194)
      return true ;
   return false ;
}

//----------------------------------------------------------------------------
bool
Config7Duality( int nInd)
{
   if ( nInd == 26  ||  nInd == 37  ||  nInd == 74  ||  nInd == 82  ||
        nInd == 88  ||  nInd == 133  ||  nInd == 161  ||  nInd == 164)
      return true ;
   return false ;
}

//----------------------------------------------------------------------------
int
TestOnNormal( const AppliedVector CompoField[], int nCompoElem)
{
  // Cerco la massima deviazione tra le normali nei punti della parte connessa
   int nI = -1, nJ = -1 ;
   double dMinCosTheta = 2 ;
   for ( int i = 0 ; i < nCompoElem ; ++ i) {
      for  ( int j = i + 1 ; j < nCompoElem ; ++ j) {
         double dCurrCos = CompoField[i].vtVec * CompoField[j].vtVec ;
         if ( dCurrCos < dMinCosTheta) {
            dMinCosTheta = dCurrCos ;
            nI = i ;
            nJ = j ;
         }
      }
   }

  // Se la massima deviazione non supera il limite non è feature
   const double SHARP_COS_SUP = 0.866 ;
   const double SHARP_COS_INF =  - 0.985 ;
   if ( dMinCosTheta >= SHARP_COS_SUP  ||  dMinCosTheta <= SHARP_COS_INF)
      return NO_FEATURE ;

  // Verifico se Edge o Corner
   // direzione perpendicolare alle normali con massima differenza (potenziale edge)
   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( CompoField[i].vtVec * vtK) ;
      if ( dAbsCurrentCos > dMaxAbsCos)
         dMaxAbsCos = dAbsCurrentCos ;
   }
   // se esiste normale diretta quasi come potenziale edge, allora corner
   const double CORNER_COS = 0.39 ;
   if ( dMaxAbsCos > CORNER_COS) 
      return CORNER ;
   else
      return EDGE ;
}

//----------------------------------------------------------------------------
bool
CanonicPlaneTest( const AppliedVector CompoField[], int nDir, double& dPos, int& nTool)
{
  // Verifico posizione dei punti
   int nI = -1 ;
   switch ( nDir) {
   case X_PLUS : case X_MINUS : nI = 0 ; break ;
   case Y_PLUS : case Y_MINUS : nI = 1 ; break ;
   case Z_PLUS : case Z_MINUS : nI = 2 ; break ;
   }
   double dPos0 = 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)
      return false ;
  // Verifico direzione normali
   Vector3d vtN ;
   switch ( nDir) {
   case X_PLUS : vtN = X_AX ; break ;
   case X_MINUS : vtN = - X_AX ; break ;
   case Y_PLUS : vtN = Y_AX ; break ;
   case Y_MINUS : vtN = - Y_AX ; break ;
   case Z_PLUS : vtN = Z_AX ; break ;
   case Z_MINUS : vtN = - Z_AX ; break ;
   }
   for ( int i = 0 ; i < 4 ; ++ i) {
      if ( CompoField[i].vtVec * vtN < 0.999)
         return false ;
      for ( int j = i + 1 ; j < 4 ; ++ j) {
         if ( CompoField[i].nPropIndex != CompoField[j].nPropIndex)
            return false ;
      }
   }
   nTool = CompoField[0].nPropIndex ;
  // Superati tutti i test
   return true ;
}

//----------------------------------------------------------------------------
bool
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 = CompoField[i].vtVec * CompoField[j].vtVec ;
         if ( dCurrCos < dMinCosTheta) {
            dMinCosTheta = dCurrCos ;
         }
      }
   }

  // se normali sparpagliate oltre limite
   if ( dMinCosTheta < dThreshold)
      return false ;

  // determino media delle normali
   vtAvg = V_NULL ;
   for ( int i = 0 ; i < nCompoElem ; ++ i)
      vtAvg += CompoField[i].vtVec ;
   vtAvg /= nCompoElem ;

   return true ;
}

//----------------------------------------------------------------------------
bool 
CreateBigTriangleXY( double dMinX, double dMaxX, double dMinY, double dMaxY, double dZ,
                     bool bNormZPlus, int nGrade, Triangle3dEx& trTria1, Triangle3dEx& trTria2)
{
  // Punti che definiscono il perimetro del rettangolo 
   Point3d ptP0( dMinX, dMinY, dZ) ;
   Point3d ptP1( dMinX, dMaxY, dZ) ;
   Point3d ptP2( dMaxX, dMaxY, dZ) ;
   Point3d ptP3( dMaxX, dMinY, dZ) ;   

  // Versore normale diretto come Z+ 
   if ( bNormZPlus) {        
      trTria1.Set( ptP0, ptP2, ptP1) ;
      trTria2.Set( ptP0, ptP3, ptP2) ; 
   }
  // Versore normale diretto come Z-
   else {            
      trTria1.Set( ptP0, ptP1, ptP2) ;
      trTria2.Set( ptP0, ptP2, ptP3) ;             
   }
            
   trTria1.SetGrade( nGrade) ;
   trTria2.SetGrade( nGrade) ;
        
   return trTria1.Validate( true)  &&  trTria2.Validate( true)  ; 
}

//----------------------------------------------------------------------------
bool 
CreateBigTriangleXZ( double dMinX, double dMaxX, double dMinZ, double dMaxZ, double dY,
                     bool bNormYPlus, int nGrade, Triangle3dEx& trTria1, Triangle3dEx& trTria2)
{
  // Punti che definiscono il perimetro del rettangolo   
   Point3d ptP0( dMinX, dY, dMinZ) ;
   Point3d ptP1( dMinX, dY, dMaxZ) ;
   Point3d ptP2( dMaxX, dY, dMaxZ) ;
   Point3d ptP3( dMaxX, dY, dMinZ) ;  

  // Versore normale diretto come Y+ 
   if ( bNormYPlus) {    
      trTria1.Set( ptP0, ptP1, ptP2) ;
      trTria2.Set( ptP0, ptP2, ptP3) ;   
   }
  // Versore normale diretto come Y-  
   else {              
      trTria1.Set( ptP0, ptP2, ptP1) ;
      trTria2.Set( ptP0, ptP3, ptP2) ;      
   }            

   trTria1.SetGrade( nGrade) ;
   trTria2.SetGrade( nGrade) ;
        
   return trTria1.Validate( true)  &&  trTria2.Validate( true)  ; 
}

//----------------------------------------------------------------------------
bool 
CreateBigTriangleYZ( double dMinY, double dMaxY, double dMinZ, double dMaxZ, double dX,
                     bool bNormXPlus, int nGrade, Triangle3dEx& trTria1, Triangle3dEx& trTria2)
{
  // Punti che definiscono il perimetro del rettangolo 
   Point3d ptP0( dX, dMinY, dMinZ) ;
   Point3d ptP1( dX, dMinY, dMaxZ) ;
   Point3d ptP2( dX, dMaxY, dMaxZ) ;
   Point3d ptP3( dX, dMaxY, dMinZ) ;  

  // Versore normale diretto come X+ 
   if ( bNormXPlus) {                                
      trTria1.Set( ptP0, ptP2, ptP1) ;
      trTria2.Set( ptP0, ptP3, ptP2) ;     
   }    
  // Versore normale diretto come X-
   else {                         
      trTria1.Set( ptP0, ptP1, ptP2) ;
      trTria2.Set( ptP0, ptP2, ptP3) ;          
   }  

   trTria1.SetGrade( nGrade) ;
   trTria2.SetGrade( nGrade) ;   
        
   return trTria1.Validate( true)  &&  trTria2.Validate( true)  ;       
}

//----------------------------------------------------------------------------
double
GetFanTotalNormalDeviation( Point3d ptP, AppliedVector Vector[7], int nVertNum)
{
   double dSum = 0 ;
   for ( int m = 0 ; m < nVertNum ; ++ m) {
      int l = ( m + 1) % nVertNum ;
      Vector3d vtNlm = ( Vector[l].ptPApp - ptP) ^ ( Vector[m].ptPApp - ptP) ;
      vtNlm.Normalize() ;
      dSum += vtNlm * Vector[l].vtVec ;
      dSum += vtNlm * Vector[m].vtVec ;
   }
   return dSum ;
}

//----------------------------------------------------------------------------
bool
VolZmap::FindAdjComp( const vector<VoxelContainer>& 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 
VolZmap::GetDexelLines( int nDir, int nPos1, int nPos2, POLYLINELIST& lstPL) const 
{
  // Se richiesti spilloni ( 0 <= nDir < 3)
   if ( nDir < 3) {
     // Controllo l'ammissibilità della griglia
      if ( nDir < 0 || nDir > 2)
         return false ;
     // Verifiche sugli indici
      if ( nPos1 < 0  ||  nPos1 >= int( m_nNx[nDir])  ||  nPos2 < 0  ||  nPos2 >= int( m_nNy[nDir]))
         return false ;
      int nPos = nPos1 + nPos2 * m_nNx[nDir] ;
      if ( nPos < 0  ||  nPos >= int( m_Values[nDir].size()))
         return false ;
     // Definisco un sistema di riferimento opportuno
      Frame3d frMapFrame ;
      if ( nDir == 0) 
         frMapFrame = m_MapFrame ;
      else if ( nDir == 1)
         frMapFrame.Set( m_MapFrame.Orig(), m_MapFrame.VersY(), m_MapFrame.VersZ(), m_MapFrame.VersX()) ;
      else if ( nDir == 2)
         frMapFrame.Set( m_MapFrame.Orig(), m_MapFrame.VersZ(), m_MapFrame.VersX(), m_MapFrame.VersY()) ;
     // Calcolo coordinate punto
      double dX = m_dStep * ( 0.5 + nPos1) ;
      double dY = m_dStep * ( 0.5 + nPos2) ;
     // Determino il punto di partenza sulla griglia 
      Point3d ptP = frMapFrame.Orig() + dX * frMapFrame.VersX() + dY * frMapFrame.VersY() ;
     // Creo le polilinee
      for ( int j = 0 ; j < int( m_Values[nDir][nPos].size()) ; j += 1) {
        // aggiungo polilinea a lista
         lstPL.emplace_back() ;
        // inserisco punti estremi
         lstPL.back().AddUPoint( 0, ptP + m_Values[nDir][nPos][j].dMin * frMapFrame.VersZ()) ;
         lstPL.back().AddUPoint( 1, ptP + m_Values[nDir][nPos][j].dMax * frMapFrame.VersZ()) ;
      }
      return true ;
   }

  // altrimenti richieste normali ( 3 <= nDir < 6)
   else {
     // riporto a indice griglia
      nDir -= 3 ;
     // Controllo l'ammissibilità della griglia
      if ( nDir < 0 || nDir > 2)
         return false ;
     // Verifiche sugli indici
      if ( nPos1 < 0  ||  nPos1 >= int( m_nNx[nDir])  ||  nPos2 < 0  ||  nPos2 >= int( m_nNy[nDir]))
         return false ;
      int nPos = nPos1 + nPos2 * m_nNx[nDir] ;
      if ( nPos < 0  ||  nPos >= int( m_Values[nDir].size()))
         return false ;
     // Definisco un sistema di riferimento opportuno
      Frame3d frMapFrame ;
      if ( nDir == 0) 
         frMapFrame = m_MapFrame ;
      else if ( nDir == 1)
         frMapFrame.Set( m_MapFrame.Orig(), m_MapFrame.VersY(), m_MapFrame.VersZ(), m_MapFrame.VersX()) ;
      else if ( nDir == 2)
         frMapFrame.Set( m_MapFrame.Orig(), m_MapFrame.VersZ(), m_MapFrame.VersX(), m_MapFrame.VersY()) ;
     // Calcolo coordinate punto
      double dX = m_dStep * ( 0.5 + nPos1) ;
      double dY = m_dStep * ( 0.5 + nPos2) ;
     // Determino il punto di partenza sulla griglia 
      Point3d ptP = frMapFrame.Orig() + dX * frMapFrame.VersX() + dY * frMapFrame.VersY() ;
     // Creo le polilinee
      for ( int j = 0 ; j < int( m_Values[nDir][nPos].size()) ; j += 1) {
        // aggiungo polilinea a lista
         lstPL.emplace_back() ;
        // calcolo e inserisco punto inizio spillone
         Point3d ptQ = ptP + m_Values[nDir][nPos][j].dMin * frMapFrame.VersZ() ;
         lstPL.back().AddUPoint( 0, ptQ) ;
        // calcolo e inserisco punto su termine sua normale
         Vector3d vtV = m_Values[nDir][nPos][j].vtMinN ;
         vtV.ToGlob( m_MapFrame) ;
         lstPL.back().AddUPoint( 1, ptQ + vtV * m_dStep / 4) ;
        // aggiungo polilinea a lista
         lstPL.emplace_back() ;
        // calcolo e inserisco punto fine spillone
         Point3d ptR = ptP + m_Values[nDir][nPos][j].dMax * frMapFrame.VersZ() ;
         lstPL.back().AddUPoint( 0, ptR) ;
        // calcolo e inserisco punto su termine sua normale
         Vector3d vtW = m_Values[nDir][nPos][j].vtMaxN ;
         vtW.ToGlob( m_MapFrame) ;
         lstPL.back().AddUPoint( 1, ptR + vtW * m_dStep / 4) ;
      }
      return true ;
   }
}

//----------------------------------------------------------------------------
bool 
VolZmap::GetBlockTriangles( int nBlock, TRIA3DEXVECTOR& vTria) const
{
  // Dexel (singola mappa)
   if ( m_nMapNum == 1) {
     // se il blocco non esiste, errore
      if ( nBlock < 0  ||  nBlock > m_nNumBlock - 1)
         return false ;
     // lancio eventuale aggiornamento pendente della grafica
      UpdateSingleMapGraphics() ;
     // copio i triangoli
      vTria = m_SingleMapTria[nBlock] ;
      return true ;
   }

  // Tridexel (tre mappe)
   else {
     // se il blocco non esiste, errore
      if ( nBlock < 0  ||  nBlock > ( m_nNumBlock + 1) - 1)
         return false ;
     // lancio eventuale aggiornamento pendente della grafica
      UpdateTripleMapGraphics() ;
     // copio i triangoli del blocco
      if ( nBlock != m_nNumBlock) {
         vTria.clear() ;
         vTria.reserve( 10000) ;
        // triangoli smooth
         for ( int nVx = 0 ; nVx < int( m_BlockSmoothTria[nBlock].size()) ; ++ nVx) {
            for ( int nTr = 0 ; nTr < int( m_BlockSmoothTria[nBlock][nVx].vTria.size()) ; ++ nTr) {
               vTria.emplace_back( m_BlockSmoothTria[nBlock][nVx].vTria[nTr]) ;
               vTria.back().ResetEdgeFlags() ;
            }
         }
        // triangoli grandi piatti
         for ( int tBl = 0 ; tBl < int( m_BlockBigTria[nBlock].size()) ; ++ tBl) {
            vTria.emplace_back( m_BlockBigTria[nBlock][tBl]) ;  
            vTria.back().ResetEdgeFlags() ;
         }
        // triangoli di feature nel blocco (ciclo sui voxel del blocco)
         for ( int t1 = 0 ;  t1 < int( m_BlockSharpTria[nBlock].size()) ; ++ t1) {
           // ciclo sulle componenti connesse del voxel
            for ( int t2 = 0 ; t2 < int( m_BlockSharpTria[nBlock][t1].vCompoTria.size()) ; ++ t2) {
              // ciclo sui triangoli delle componenti connesse 
               for ( int t3 = 0 ; t3 < int( m_BlockSharpTria[nBlock][t1].vCompoTria[t2].size()) ; ++ t3) {
                 // aggiungo triangolo alla lista
                  vTria.emplace_back( m_BlockSharpTria[nBlock][t1].vCompoTria[t2][t3]) ;
                  vTria.back().ResetEdgeFlags() ;
               }
            }
         }
      }
     // blocco speciale aggiunto per i triangoli feature a cavallo dei blocchi
      else {
         vTria.clear() ;
         vTria.reserve( 1000) ;
         for ( int t = 0 ; t < int( m_InterBlockSharpTria.size()) ; ++ t) {
            for ( int t1 = 0 ; t1 < int( m_InterBlockSharpTria[t].size()) ; ++ t1) {
               for ( int t2 = 0 ; t2 < int( m_InterBlockSharpTria[t][t1].vCompoTria.size()) ; ++ t2) {
                  for ( int t3 = 0 ; t3 < int( m_InterBlockSharpTria[t][t1].vCompoTria[t2].size()) ; ++ t3) {
                     if ( m_InterBlockSharpTria[t][t1].vCompoTria[t2][t3].GetArea() > SQ_EPS_SMALL) {
                       // aggiungo triangolo alla lista
                        vTria.emplace_back( m_InterBlockSharpTria[t][t1].vCompoTria[t2][t3]) ;
                        vTria.back().ResetEdgeFlags() ;
                     }
                  }
               }
            }
         }
      }
      return true ;
   }
}

//----------------------------------------------------------------------------
ISurfTriMesh*
VolZmap::GetSurfTriMesh( void) const
{
  // controllo che lo Zmap sia valido
   if ( ! IsValid())
      return nullptr ;

  // inizializzo la superficie
   PtrOwner<SurfTriMesh> pStm( CreateBasicSurfTriMesh()) ;
   if ( IsNull( pStm)  ||  ! pStm->Init( 3, 1))
      return nullptr ;
   PointGrid3d VertGrid ; VertGrid.Init( 50000) ;

  // ciclo lungo i blocchi dello Zmap
   for ( int nB = 0 ; nB < GetBlockCount() ; ++ nB) {
      TRIA3DEXVECTOR vTria ;
      GetBlockTriangles( nB, vTria) ;
      if ( ! pStm->AddTriaFromZMap( vTria, VertGrid))
         return nullptr ;
   }

  // sistemo la topologia
   if ( ! pStm->AdjustTopologyFromZMap())
      return nullptr ;

   return ( Release( pStm)) ;
}

//----------------------------------------------------------------------------
int
VolZmap::GetBlockCount( void) const
{
   return ( m_nNumBlock + ( m_nMapNum == 1 ? 0 : 1)) ;
}

//----------------------------------------------------------------------------
int
VolZmap::GetBlockUpdatingCounter( int nBlock) const
{
  // se il blocco non esiste, errore
   if ( nBlock < 0  ||  nBlock >= int( m_BlockUpdatingCounter.size()))
      return -1 ;
  // Lancio eventuale aggiornamento pendente della grafica
   if ( m_nMapNum == 1)
      UpdateSingleMapGraphics() ;
   else
      UpdateTripleMapGraphics() ;
  // restituisco il suo indice di aggiornamento
   return m_BlockUpdatingCounter[nBlock] ;
}

//----------------------------------------------------------------------------
bool
VolZmap::UpdateSingleMapGraphics( void) const
{
   const int MAX_DIM_CHUNK = 128 ;

  // Ciclo sui blocchi
   for ( int t = 0 ; t < m_nNumBlock ; ++ t) {

     // Se il blocco deve essere aggiornato, eseguo
      if ( m_BlockToUpdate[t]) {

         m_SingleMapTria[t].clear() ;
        
        // Calcolo posizione del blocco nella griglia
         int nIBlock = int( t) % int( m_nFracLin[0]) ;
         int nJBlock = int( t) / int( m_nFracLin[0]) ;

        // Calcolo limiti per l'indice i
         int nStartI = nIBlock * int( m_nVoxNumPerBlock) * m_nDexVoxRatio ;
         int nEndI = ( nIBlock + 1 == int( m_nFracLin[0]) ? 
                       int( m_nNx[0]) : ( nIBlock + 1) * int( m_nVoxNumPerBlock)) * m_nDexVoxRatio ;

        // Calcolo limiti per l'indice j
         int nStartJ = nJBlock * int( m_nVoxNumPerBlock) * m_nDexVoxRatio ;
         int nEndJ = ( nJBlock + 1 == int( m_nFracLin[1]) ? 
                       int( m_nNy[0]) : ( nJBlock + 1) * int( m_nVoxNumPerBlock)) * m_nDexVoxRatio ;

        // Ciclo su i e j 
         for ( int i = nStartI ; i < nEndI ; i += MAX_DIM_CHUNK) {
            int nDimChunkX = min( MAX_DIM_CHUNK, nEndI - i) ;
            for ( int j = nStartJ ; j < nEndJ ; j += MAX_DIM_CHUNK) {
               int nDimChunkY = min( MAX_DIM_CHUNK, nEndJ - j) ;
               GetChunkPrisms( i, j, nDimChunkX, nDimChunkY, MAX_DIM_CHUNK, t) ;
            } 
         }

         m_BlockToUpdate[t] = false ;
         ++ m_BlockUpdatingCounter[t] ;
      }
   }

   return true ;
}

//----------------------------------------------------------------------------
bool
VolZmap::GetChunkPrisms( int nPos1, int nPos2, int nDim1, int nDim2, int nDimChk, int nBlock) const
{
  // determino se è un semplice parallelepipedo
   bool bIsSimple = true ;
   double dBotZ = 0 ;
   double dTopZ = 0 ;
   for ( int i = 0 ; i < nDim1 && bIsSimple ; ++ i) {
      for ( int j = 0 ; j < nDim2 && bIsSimple ; ++ j) {
         int nPos = ( nPos1 + i) + ( nPos2 + j) * m_nNx[0] ;
         if ( nPos > int( m_nDim[0])  ||
              int( m_Values[0][nPos].size()) != 1)
            bIsSimple = false ;
         else if ( i == 0  &&  j == 0) {
            dBotZ = m_Values[0][nPos][0].dMin ;
            dTopZ = m_Values[0][nPos][0].dMax ;
         }
         else if ( abs( m_Values[0][nPos][0].dMin - dBotZ) > EPS_SMALL  ||
                   abs( m_Values[0][nPos][0].dMax - dTopZ) > EPS_SMALL)
            bIsSimple = false ;
      }
   }

  // se semplice parallelepipedo
   if ( bIsSimple) {
      CalcChunkPrisms( nPos1, nPos2, nDim1, nDim2, nBlock) ;
   }
  // se chunk di dimensioni accettabili
   else if ( nDimChk >= 4) {
      int nNewDimChk = nDimChk / 2 ;
      for ( int i = nPos1 ; i < int( nPos1 + nDim1) ; i += nNewDimChk) {
         int nDimChunkX = min( nNewDimChk, int( nPos1 + nDim1) - i) ;
         for ( int j = nPos2 ; j < int( nPos2 + nDim2) ; j += nNewDimChk) {
            int nDimChunkY = min( nNewDimChk, int( nPos2 + nDim2) - j) ;
            GetChunkPrisms( i, j, nDimChunkX, nDimChunkY, nNewDimChk, nBlock) ;
         }
      }
   }
  // altrimenti
   else {
     // elaboro ogni singolo dexel
      for ( int i = 0 ; i < nDim1 ; ++ i) {
         for ( int j = 0 ; j < nDim2 ; ++ j) {
            CalcDexelPrisms( nPos1 + i, nPos2 + j, nBlock) ;
         }
      }
   }
   return true ;
}

//----------------------------------------------------------------------------
bool
VolZmap::CalcChunkPrisms( int nPos1, int nPos2, int nDim1, int nDim2, int nBlock) const
{
  // verifiche sugli indici
   if ( nPos1 < 0  ||  nPos1 + nDim1 > int( m_nNx[0])  ||  nPos2 < 0  ||  nPos2 + nDim2 > int( m_nNy[0]))
      return false ;
   int nPos = nPos1 + nPos2 * m_nNx[0] ;
   if ( nPos < 0  ||  nPos >= int( m_nDim[0]))
      return false ;

  // calcolo coordinate punti
   double dX = m_dStep * nPos1 ;
   double dY = m_dStep * nPos2 ;
   Point3d ptP1 = m_MapFrame.Orig() + dX * m_MapFrame.VersX() + dY * m_MapFrame.VersY() ;
   Point3d ptP2 = ptP1 + nDim1 * m_dStep * m_MapFrame.VersX() ;
   Point3d ptP3 = ptP2 + nDim2 * m_dStep * m_MapFrame.VersY() ;
   Point3d ptP4 = ptP1 + nDim2 * m_dStep * m_MapFrame.VersY() ;

  // creo le facce sopra e sotto
   Vector3d vtDZt = m_Values[0][nPos][0].dMax * m_MapFrame.VersZ() ;
   Vector3d vtDZb = m_Values[0][nPos][0].dMin * m_MapFrame.VersZ() ;
   // faccia superiore P1t->P2t->P3t->P4t : sempre visibile
   m_SingleMapTria[nBlock].emplace_back() ;
   m_SingleMapTria[nBlock].back().Set( ptP1 + vtDZt, ptP2 + vtDZt, ptP3 + vtDZt, m_MapFrame.VersZ()) ;
   m_SingleMapTria[nBlock].emplace_back() ;
   m_SingleMapTria[nBlock].back().Set( ptP3 + vtDZt, ptP4 + vtDZt, ptP1 + vtDZt, m_MapFrame.VersZ()) ;
   // faccia inferiore P1b->P4b->P3b->P2b : sempre visibile
   m_SingleMapTria[nBlock].emplace_back() ;
   m_SingleMapTria[nBlock].back().Set( ptP1 + vtDZb, ptP4 + vtDZb, ptP3 + vtDZb, - m_MapFrame.VersZ()) ;
   m_SingleMapTria[nBlock].emplace_back() ;
   m_SingleMapTria[nBlock].back().Set( ptP3 + vtDZb, ptP2 + vtDZb, ptP1 + vtDZb, - m_MapFrame.VersZ()) ;

  // creo le facce laterali
   for ( int j = 0 ; j < nDim2 ; ++ j) {
      int nPosD = nPos + nDim1 - 1 + j * m_nNx[0] ;
      int nPosEst = ( nPos1 + nDim1 - 1 < int( m_nNx[0] - 1) ? nPosD + 1 : - 1) ;
      Point3d ptP2D = ptP2 + j * m_dStep * m_MapFrame.VersY() ;
      Point3d ptP3D = ptP2D + m_dStep * m_MapFrame.VersY() ;
      AddDexelSideFace( nPosD, nPosEst, ptP2D, ptP3D, m_MapFrame.VersZ(), m_MapFrame.VersX(), nBlock) ;
   }
   for ( int i = 0 ; i < nDim1 ; ++ i) {
      int nPosD = nPos + ( nDim2 - 1) * m_nNx[0] + i ;
      int nPosNord = ( nPos2 + nDim2 - 1 < int( m_nNy[0] - 1) ? nPosD + m_nNx[0] : - 1) ;
      Point3d ptP4D = ptP4 + i * m_dStep * m_MapFrame.VersX() ;
      Point3d ptP3D = ptP4D + m_dStep * m_MapFrame.VersX() ;
      AddDexelSideFace( nPosD, nPosNord, ptP3D, ptP4D, m_MapFrame.VersZ(), m_MapFrame.VersY(), nBlock) ;
   }
   for ( int j = 0 ; j < nDim2 ; ++ j) {
      int nPosD = nPos + j * m_nNx[0] ;
      int nPosWest = ( nPos1 > 0 ? nPosD - 1 : - 1) ;
      Point3d ptP1D = ptP1 + j * m_dStep * m_MapFrame.VersY() ;
      Point3d ptP4D = ptP1D + m_dStep * m_MapFrame.VersY() ;
      AddDexelSideFace( nPosD, nPosWest, ptP4D, ptP1D, m_MapFrame.VersZ(), - m_MapFrame.VersX(), nBlock) ;
   }
   for ( int i = 0 ; i < nDim1 ; ++ i) {
      int nPosD = nPos + i ;
      int nPosSud = ( nPos2 > 0 ? nPosD - m_nNx[0] : - 1) ;
      Point3d ptP1D = ptP1 + i * m_dStep * m_MapFrame.VersX() ;
      Point3d ptP2D = ptP1D + m_dStep * m_MapFrame.VersX() ;
      AddDexelSideFace( nPosD, nPosSud, ptP1D, ptP2D, m_MapFrame.VersZ(), - m_MapFrame.VersY(), nBlock) ;
   }
//
   return true ;
}

//----------------------------------------------------------------------------
bool
VolZmap::CalcDexelPrisms( int nPos1, int nPos2, int nBlock) const
{
  // verifiche sugli indici
   if ( nPos1 < 0  ||  nPos1 >= int( m_nNx[0])  ||  nPos2 < 0  ||  nPos2 >= int( m_nNy[0]))
      return false ;
   int nPos = nPos1 + nPos2 * m_nNx[0] ;
   if ( nPos < 0  ||  nPos >= int( m_nDim[0]))
      return false ;

  // calcolo coordinate punto
   double dX = m_dStep * nPos1 ;
   double dY = m_dStep * nPos2 ;
   Point3d ptP1 = m_MapFrame.Orig() + dX * m_MapFrame.VersX() + dY * m_MapFrame.VersY() ;
   Point3d ptP2 = ptP1 + m_dStep * m_MapFrame.VersX() ;
   Point3d ptP3 = ptP2 + m_dStep * m_MapFrame.VersY() ;
   Point3d ptP4 = ptP1 + m_dStep * m_MapFrame.VersY() ;

  // creo le facce sopra e sotto di ogni intervallo (sempre visibili)
   for ( int i = 0 ; i < int( m_Values[0][nPos].size()) ; i += 1) {
      Vector3d vtDZt = m_Values[0][nPos][i].dMax * m_MapFrame.VersZ() ;
      Vector3d vtDZb = m_Values[0][nPos][i].dMin * m_MapFrame.VersZ() ;
     // faccia superiore P1t->P2t->P3t->P4t : sempre visibile
      m_SingleMapTria[nBlock].emplace_back() ;
      m_SingleMapTria[nBlock].back().Set( ptP1 + vtDZt, ptP2 + vtDZt, ptP3 + vtDZt, m_MapFrame.VersZ()) ;
      m_SingleMapTria[nBlock].emplace_back() ;
      m_SingleMapTria[nBlock].back().Set( ptP3 + vtDZt, ptP4 + vtDZt, ptP1 + vtDZt, m_MapFrame.VersZ()) ;
     // faccia inferiore P1b->P4b->P3b->P2b : sempre visibile
      m_SingleMapTria[nBlock].emplace_back() ;
      m_SingleMapTria[nBlock].back().Set( ptP1 + vtDZb, ptP4 + vtDZb, ptP3 + vtDZb, - m_MapFrame.VersZ()) ;
      m_SingleMapTria[nBlock].emplace_back() ;
      m_SingleMapTria[nBlock].back().Set( ptP3 + vtDZb, ptP2 + vtDZb, ptP1 + vtDZb, - m_MapFrame.VersZ()) ;
   }

  // creo le facce laterali
   int nPosEst = ( nPos1 < int( m_nNx[0] - 1) ? nPos + 1 : - 1) ;
   AddDexelSideFace( nPos, nPosEst, ptP2, ptP3, m_MapFrame.VersZ(), m_MapFrame.VersX(), nBlock) ;
   int nPosNord = ( nPos2 < int( m_nNy[0] - 1) ? nPos + m_nNx[0] : - 1) ;
   AddDexelSideFace( nPos, nPosNord, ptP3, ptP4, m_MapFrame.VersZ(), m_MapFrame.VersY(), nBlock) ;
   int nPosWest = ( nPos1 > 0 ? nPos - 1 : - 1) ;
   AddDexelSideFace( nPos, nPosWest, ptP4, ptP1, m_MapFrame.VersZ(), - m_MapFrame.VersX(), nBlock) ;
   int nPosSud = ( nPos2 > 0 ? nPos - m_nNx[0] : - 1) ;
   AddDexelSideFace( nPos, nPosSud, ptP1, ptP2, m_MapFrame.VersZ(), - m_MapFrame.VersY(), nBlock) ;

   return true ;
}

//----------------------------------------------------------------------------
bool
VolZmap::AddDexelSideFace( int nPos, int nPosAdj, const Point3d& ptP, const Point3d& ptQ,
                           const Vector3d& vtZ, const Vector3d& vtNorm, int nBlock) const
{
   Intervals intFace ;
   for ( int i = 0 ; i < int( m_Values[0][nPos].size()) ; i += 1)
       intFace.Add( m_Values[0][nPos][i].dMin, m_Values[0][nPos][i].dMax) ;
   if ( nPosAdj > 0) {
      for ( int i = 0 ; i < int( m_Values[0][nPosAdj].size()) ; i += 1)
          intFace.Subtract( m_Values[0][nPosAdj][i].dMin, m_Values[0][nPosAdj][i].dMax) ;
   }
   double dMin, dMax ;
   bool bFound = intFace.GetFirst( dMin, dMax) ;
   while ( bFound) {
      Vector3d vtDZt = dMax * vtZ ;
      Vector3d vtDZb = dMin * vtZ ;
      m_SingleMapTria[nBlock].emplace_back() ;
      m_SingleMapTria[nBlock].back().Set( ptP + vtDZb, ptQ + vtDZb, ptQ + vtDZt, vtNorm) ;
      m_SingleMapTria[nBlock].emplace_back() ;
      m_SingleMapTria[nBlock].back().Set( ptQ + vtDZt, ptP + vtDZt, ptP + vtDZb, vtNorm) ;
      bFound = intFace.GetNext( dMin, dMax) ;
   }
   return true ;
}

//----------------------------------------------------------------------------
bool 
VolZmap::UpdateTripleMapGraphics( void) const
{
  // se non ci sono blocchi da aggiornare, esco
   if ( find( begin( m_BlockToUpdate), end( m_BlockToUpdate), true) == end( m_BlockToUpdate))
      return true ;

  // contenitori temporanei
   vector<VoxelContainer> vVoxContainerVec ;
   vVoxContainerVec.resize( m_nNumBlock) ;
   SharpTriaMatrix VecTriHold ;
   VecTriHold.resize( m_nNumBlock) ;

  // Ciclo sui blocchi per eliminare le slice fra blocchi da aggiornare
   for ( int t = 0 ; t < m_nNumBlock ; ++ t) {
      for ( auto it = m_SliceXY[t].begin() ; it != m_SliceXY[t].end() ;) {
         int nSlIJK[3] ;
         if ( GetVoxIJKFromN( it->first, nSlIJK[0], nSlIJK[1], nSlIJK[2])) {
            int nBlockIJK[3] ;
            if ( GetVoxelBlockIJK( nSlIJK, nBlockIJK)) {
               int nLimits[6] ;
               int nDeltaIndex[3] ;
               if ( GetBlockLimitsIJK( nBlockIJK, nLimits)  &&
                    IsAVoxelOnBoundary( nLimits, nSlIJK, nDeltaIndex)) {
                  for ( int nInd = 0 ; nInd < 3 ; ++ nInd)
                     nSlIJK[nInd] += nDeltaIndex[nInd] ;
                  int nAdBlockIJK[3] ;
                  int nAdBlockNum ;
                  if ( GetVoxelBlockIJK( nSlIJK, nAdBlockIJK)  && 
                       GetBlockNFromIJK( nAdBlockIJK, nAdBlockNum)  && 
                       m_BlockToUpdate[nAdBlockNum]) {
                     it = m_SliceXY[t].erase( it) ;
                     continue ;
                  }
               }
            }
         }
         ++ it ;
      }  
      for ( auto it = m_SliceXZ[t].begin() ; it != m_SliceXZ[t].end() ;) {
         int nSlIJK[3] ;
         if ( GetVoxIJKFromN( it->first, nSlIJK[0], nSlIJK[1], nSlIJK[2])) {
            int nBlockIJK[3] ;
            if ( GetVoxelBlockIJK( nSlIJK, nBlockIJK)) { 
               int nLimits[6] ;
               int nDeltaIndex[3] ;
               if ( GetBlockLimitsIJK( nBlockIJK, nLimits)  &&  
                    IsAVoxelOnBoundary( nLimits, nSlIJK, nDeltaIndex)) {
                  for ( int nInd = 0 ; nInd < 3 ; ++ nInd)
                     nSlIJK[nInd] += nDeltaIndex[nInd] ; 
                  int nAdBlockIJK[3] ;
                  int nAdBlockNum ;
                  if ( GetVoxelBlockIJK( nSlIJK, nAdBlockIJK)  && 
                       GetBlockNFromIJK( nAdBlockIJK, nAdBlockNum)  && 
                       m_BlockToUpdate[nAdBlockNum]) {
                     it = m_SliceXZ[t].erase( it) ;
                     continue ;
                  }
               }
            }
         }
         ++ it ;
      }
      for ( auto it = m_SliceYZ[t].begin() ; it != m_SliceYZ[t].end() ;) {
         int nSlIJK[3] ;
         if ( GetVoxIJKFromN( it->first, nSlIJK[0], nSlIJK[1], nSlIJK[2])) {
            int nBlockIJK[3] ;
            if ( GetVoxelBlockIJK( nSlIJK, nBlockIJK)) { 
               int nLimits[6] ;
               int nDeltaIndex[3] ;
               if ( GetBlockLimitsIJK( nBlockIJK, nLimits)  &&  
                    IsAVoxelOnBoundary( nLimits, nSlIJK, nDeltaIndex)) {
                  for ( int nInd = 0 ; nInd < 3 ; ++ nInd)
                     nSlIJK[nInd] += nDeltaIndex[nInd] ; 
                  int nAdBlockIJK[3] ;
                  int nAdBlockNum ;
                  if ( GetVoxelBlockIJK( nSlIJK, nAdBlockIJK)  && 
                       GetBlockNFromIJK( nAdBlockIJK, nAdBlockNum)  && 
                       m_BlockToUpdate[nAdBlockNum]) {
                     it = m_SliceYZ[t].erase( it) ;
                     continue ;
                  }
               }
            }
         }
         ++ it ;
      }
   }

  // Standard è multithread
   constexpr bool MULTITHREAD = true ;
   if ( MULTITHREAD) {

     // Calcolo i triangoli sui blocchi
      int nBlockUpdated = 0 ;
      vector< future<bool>> vRes ;
      vRes.resize( m_nNumBlock) ;
      for ( int i = 0 ; i < m_nNumBlock ; ++ i) {
        // Se il blocco deve essere processato
         if ( m_BlockToUpdate[i]) {
           // processo ...
            ++ nBlockUpdated ;
            vRes[i] = async( launch::async, &VolZmap::ExtMarchingCubes, this, i, ref( vVoxContainerVec[i])) ;
         }
      }
      bool bOk = true ;
      int nTerminated = 0 ;
      while ( nTerminated < nBlockUpdated) {
         for ( int i = 0 ; i < m_nNumBlock ; ++ i) {
            if ( m_BlockToUpdate[i]  &&  vRes[i].valid()  &&  vRes[i].wait_for( chrono::nanoseconds{ 1}) == future_status::ready) {
               bOk = vRes[i].get()  &&  bOk ;
               ++ nTerminated ;
            }
         }
      }
   }
   else {
     // Calcolo i triangoli sui blocchi
      bool bOk = true ;
      for ( int i = 0 ; i < m_nNumBlock ; ++ i) {
        // Se il blocco deve essere processato
         if ( m_BlockToUpdate[i]) {
           // processo ...
            bOk = ExtMarchingCubes( i, vVoxContainerVec[i]) && bOk ;
         }
      }
   }

  // Regolarizzo la catena
   RegulateFeaturesChain( vVoxContainerVec) ;

  // Costruisco i triangoli di feature
   bool bCalcInterBlock = false ;
   for ( int t = 0 ; t < m_nNumBlock ; ++ t) {
     // Se il blocco è stato processato 
      if ( m_BlockToUpdate[t]) {
         CreateSharpFeatureTriangle( t, vVoxContainerVec[t]) ;
        // Flipping fra voxel interni
         FlipEdgesII( t) ;
         bCalcInterBlock = true ;
         m_BlockToUpdate[t] = false ;
         ++ m_BlockUpdatingCounter[t] ;
        // Sistemo le normali ai vertici (ciclo sui voxel)
         for ( int t1 = 0 ;  t1 < int( m_BlockSharpTria[t].size()) ; ++ t1) {
           // ciclo sulle componenti connesse del voxel
            for ( int t2 = 0 ; t2 < int( m_BlockSharpTria[t][t1].vCompoTria.size()) ; ++ t2) {
              // ciclo sui triangoli delle componenti connesse 
               for ( int t3 = 0 ; t3 < int( m_BlockSharpTria[t][t1].vCompoTria[t2].size()) ; ++ t3) {
                 // Controllo normali
                  Vector3d vtN = m_BlockSharpTria[t][t1].vCompoTria[t2][t3].GetN() ;
                  bool bNormN = vtN.IsNormalized() ;
                  for ( int nV = 0 ; nV < 3 ; ++ nV) {
                     Vector3d vtNV = m_BlockSharpTria[t][t1].vCompoTria[t2][t3].GetVertexNorm( nV) ;
                     bool bNormV = vtNV.IsNormalized() ;
                     if ( bNormN  &&  bNormV  &&  vtN * vtNV < 0.7)
                        m_BlockSharpTria[t][t1].vCompoTria[t2][t3].SetVertexNorm( nV, vtN) ;
                  }
               }
            }
         }
      }
   }

  // Calcolo i triangoli di frontiera tra feature di blocchi diversi (lavoro su copia e conservo gli originali)
   SharpTriaMatrix InterBlockTria ;
   if ( bCalcInterBlock) {
     // Eseguo
      m_InterBlockSharpTria = m_InterBlockOriginalSharpTria ;   
      FlipEdgesBB() ;
      ++ m_BlockUpdatingCounter[m_nNumBlock] ;
     // Sistemo le normali ai vertici
      for ( int t = 0 ; t < int( m_InterBlockSharpTria.size()) ; ++ t) {
         for ( int t1 = 0 ; t1 < int( m_InterBlockSharpTria[t].size()) ; ++ t1) {
            for ( int t2 = 0 ; t2 < int( m_InterBlockSharpTria[t][t1].vCompoTria.size()) ; ++ t2) {
               for ( int t3 = 0 ; t3 < int( m_InterBlockSharpTria[t][t1].vCompoTria[t2].size()) ; ++ t3) {
                  if ( m_InterBlockSharpTria[t][t1].vCompoTria[t2][t3].GetArea() > SQ_EPS_SMALL) {
                    // Controllo normali
                     Vector3d vtN = m_InterBlockSharpTria[t][t1].vCompoTria[t2][t3].GetN() ;
                     bool bNormN = vtN.IsNormalized() ;
                     for ( int nV = 0 ; nV < 3 ; ++ nV) {
                        Vector3d vtNV = m_InterBlockSharpTria[t][t1].vCompoTria[t2][t3].GetVertexNorm( nV) ;
                        bool bNormV = vtNV.IsNormalized() ;
                        if ( bNormN  &&  bNormV  &&  vtN * vtNV < 0.7)
                           m_InterBlockSharpTria[t][t1].vCompoTria[t2][t3].SetVertexNorm( nV, vtN) ;
                     }
                  }
               }
            }
         }
      }
   }

   return true ;
}

//----------------------------------------------------------------------------
bool
VolZmap::ExtMarchingCubes( int nBlock, 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 
  // 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() ;
   m_InterBlockOriginalSharpTria[nBlock].clear() ;
   m_BlockSharpTria[nBlock].clear() ;
   m_BlockSmoothTria[nBlock].clear() ;
   m_BlockBigTria[nBlock].clear() ;

  // Unordered Map per la riduzione del numero di triangoli
   int nDim = m_nVoxNumPerBlock * m_nVoxNumPerBlock ;
   FlatVoxelContainer VoxContXZInf( nDim) ;
   FlatVoxelContainer VoxContXZSup( nDim) ;
   FlatVoxelContainer VoxContXYInf( nDim) ;
   FlatVoxelContainer VoxContXYSup( nDim) ;
   FlatVoxelContainer VoxContYZInf( nDim) ;
   FlatVoxelContainer VoxContYZSup( nDim) ;
   
  // Unordered map per la coerenza topologica nel blocco
   InterVoxMatter SliceXY( 200) ;
   InterVoxMatter SliceXZ( 200) ;
   InterVoxMatter SliceYZ( 200) ;
 
  // 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) {
            
            if ( m_nShape == BOX  &&  ! IsVoxelOnBoxEdge( i, j, k))
               continue ;
           // Classificazione dei vertici: interni o esterni al materiale 
            int nIndex = CalcIndex( i, j, k) ;

           // Se vi è qualche intersezione fra segmenti e superficie
           // continuo altrimenti passo al prossimo voxel. 
            if ( EdgeTable[nIndex] == 0) 
               continue ;

           // Indici i,j,k dei vertici
            int IndexCorner[8][3] = {
               { i, j, k},
               { i + 1, j, k},
               { i + 1, j + 1, k},
               { i, j + 1, k},
               { i, j, k + 1},
               { i + 1, j, k + 1},
               { i + 1, j + 1, k + 1},
               { i, j + 1, k + 1}
            } ;

            static int intersections[12][2] = {
               { 0, 1 }, { 1, 2 }, { 3, 2 }, { 0, 3 }, { 4, 5 }, { 5, 6 },
               { 7, 6 }, { 4, 7 }, { 0, 4 }, { 1, 5 }, { 2, 6 }, { 3, 7 }
            } ;

           // Array di strutture punto di intersezione e normale alla superficie in esso.
            AppliedVector VecField[12] ;

           // Flag di regolarità dei campi scalare e vettoriale
            bool bReg = true ;

           // Ciclo sui segmenti
            for ( int EdgeIndex = 0 ; EdgeIndex < 12 ; ++ EdgeIndex) {
              // Se il segmento non attraversa la superficie passo al successivo
               if ( ( EdgeTable[nIndex] & ( 1 << EdgeIndex)) == 0)
                  continue ;
              // Indici per linee di griglia sui vertici
               int n1 = intersections[EdgeIndex][0] ;
               int n2 = intersections[EdgeIndex][1] ;
              // Flag posizione corner
               bool bN1 = ( ( nIndex & ( 1 << n1)) != 0) ;
              // Determino con precisione il punto di intersezione sullo spigolo,
              // se i campi scalare e vettoriale non sono regolari bReg diviene falso.
               if ( ! IntersPos( IndexCorner[n1], IndexCorner[n2], bN1, VecField[EdgeIndex]))
                  bReg = false ;
            }

            
           // 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.
            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[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]] ;
              // 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) {
                  CompoTriVert[nComp][nVert] =   VecField[TriangleTableEn[nIndex][0][nStdTabOff + nVert + 2]] ;
                  CompoTriVert[nComp][nVert+1] = VecField[TriangleTableEn[nIndex][0][nStdTabOff + nVert + 1]] ;
                  CompoTriVert[nComp][nVert+2] = VecField[TriangleTableEn[nIndex][0][nStdTabOff + nVert]] ;
               }

              // Aggiorno gli offsets per raggiungere i vertici della componente successiva.
               nExtTabOff += nVertComp[nComp] ;
               nStdTabOff += 3 * ( nVertComp[nComp] - 2) ;
            }

           // Controllo se il voxel ha una sola faccia che giace in un piano canonico e quindi ha gestione speciale
            if ( m_nShape != BOX  &&  m_nShape != OFFSET) {
              // Faccia XY normale Z+
               if ( nIndex == 15) {
                  int nTool ; double dPos ;
                  if ( CanonicPlaneTest( CompoVert[0], Z_PLUS, dPos, nTool)) {
                     int nN ; GetVoxNFromIJK( i, j, k, nN) ;
                     VoxContXYSup.emplace( nN, HeigthAndColor( nTool, dPos)) ;
                     continue ;
                  }
               }
              // Faccia YZ normale X+
               else if ( nIndex == 153) {
                  int nTool ; double dPos ;
                  if ( CanonicPlaneTest( CompoVert[0], X_PLUS, dPos, nTool)) {
                     int nN ; GetVoxNFromIJK( i, j, k, nN) ;
                     VoxContYZSup.emplace( nN, HeigthAndColor( nTool, dPos)) ;
                     continue ;
                  }
               }
              // Faccia ZX normale Y+
               else if ( nIndex == 51) {
                  int nTool ; double dPos ;
                  if ( CanonicPlaneTest( CompoVert[0], Y_PLUS, dPos, nTool)) {
                     int nN ; GetVoxNFromIJK( i, j, k, nN) ;
                     VoxContXZSup.emplace( nN, HeigthAndColor( nTool, dPos)) ;
                     continue ;
                  }
               }
              // Faccia YX normale Z-
               else if ( nIndex == 240) {
                  int nTool ; double dPos ;
                  if ( CanonicPlaneTest( CompoVert[0], Z_MINUS, dPos, nTool)) {
                     int nN ; GetVoxNFromIJK( i, j, k, nN) ;
                     VoxContXYInf.emplace( nN, HeigthAndColor( nTool, dPos)) ;
                     continue ;
                  }
               }
              // Faccia ZY normale X-
               else if ( nIndex == 102) {
                  int nTool ; double dPos ;
                  if ( CanonicPlaneTest( CompoVert[0], X_MINUS, dPos, nTool)) {
                     int nN ; GetVoxNFromIJK( i, j, k, nN) ;  
                     VoxContYZInf.emplace( nN, HeigthAndColor( nTool, dPos)) ;
                     continue ;
                  }
               }
              // Faccia XZ normale Y-
               else if ( nIndex == 204) {
                  int nTool ; double dPos ;
                  if ( CanonicPlaneTest( CompoVert[0], Y_MINUS, dPos, nTool)) {
                     int nN ; GetVoxNFromIJK( i, j, k, nN) ;
                     VoxContXZInf.emplace( nN, HeigthAndColor( nTool, dPos)) ;
                     continue ;
                  }
               }
            }
               
            // Configurazione 3
            if ( nAllConfig[nIndex] == 3) {
               // Test sulla topologia
               bool bDefTopology = false ;
               bool bMatOnSlice = false ;
               int nCount = 0 ;
               while ( nIndexConfig3[nCount] != nIndex)
                  ++ nCount ;
               // Vedo se la topologia è definita: se sì uso l'informazione
               // passata dall'altro voxel, altrimenti la calcolo
               int nIJKSl[3] = { ( nAdjVox3[nCount] != 1 ? i : i + 1),
                                 ( nAdjVox3[nCount] != 2 ? j : j + 1),
                                 ( nAdjVox3[nCount] != 3 ? k : k + 1) } ;
               int nSliceN = -1 ;
               int nSlBlockN = -1 ;
               if ( GetVoxNFromIJK(nIJKSl[0], nIJKSl[1], nIJKSl[2], nSliceN)) {
                  int nSlBlockIJK[3] ;
                  GetVoxelBlockIJK(nIJKSl, nSlBlockIJK) ;
                  if ( abs( nAdjVox3[nCount]) == 1) {
                     auto it = SliceYZ.find( nSliceN) ;
                     if ( it != SliceYZ.end()) {
                        bMatOnSlice = it->second ;
                        bDefTopology = true ;
                     }
                     if ( GetBlockNFromIJK( nSlBlockIJK, nSlBlockN)) {
                        while ( m_SliceFlag.test_and_set( memory_order_acquire))
                           m_SliceFlag.wait( true, memory_order_relaxed) ;
                        auto it = m_SliceYZ[nSlBlockN].find( nSliceN) ;
                        if ( it != m_SliceYZ[nSlBlockN].end()) {
                           bMatOnSlice = it->second ;
                           bDefTopology = true ;
                        }
                        m_SliceFlag.clear( memory_order_release) ;
                        m_SliceFlag.notify_one() ;
                     }
                  }
                  else if ( abs( nAdjVox3[nCount]) == 2) {
                     auto it = SliceXZ.find( nSliceN) ;
                     if ( it != SliceXZ.end()) {
                        bMatOnSlice = it->second ;
                        bDefTopology = true ;
                     }
                     if ( GetBlockNFromIJK( nSlBlockIJK, nSlBlockN)) {
                        while ( m_SliceFlag.test_and_set( memory_order_acquire))
                           m_SliceFlag.wait( true, memory_order_relaxed) ;
                        auto it = m_SliceXZ[nSlBlockN].find( nSliceN) ;
                        if ( it != m_SliceXZ[nSlBlockN].end()) {
                           bMatOnSlice = it->second ;
                           bDefTopology = true ;
                        }
                        m_SliceFlag.clear( memory_order_release) ;
                        m_SliceFlag.notify_one() ;
                     }
                  }
                  else if ( abs( nAdjVox3[nCount]) == 3) {
                     auto it = SliceXY.find( nSliceN) ;
                     if ( it != SliceXY.end()) {
                        bMatOnSlice = it->second ;
                        bDefTopology = true ;
                     }
                     if ( GetBlockNFromIJK( nSlBlockIJK, nSlBlockN)) {
                        while ( m_SliceFlag.test_and_set( memory_order_acquire))
                           m_SliceFlag.wait( true, memory_order_relaxed) ;
                        auto it = m_SliceXY[nSlBlockN].find( nSliceN) ;
                        if ( it != m_SliceXY[nSlBlockN].end()) {
                           bMatOnSlice = it->second ;
                           bDefTopology = true ;
                        }
                        m_SliceFlag.clear( memory_order_release) ;
                        m_SliceFlag.notify_one() ;
                     }
                  }
               }
              // La topologia è indefinita: la calcolo
               if ( ! bDefTopology  &&  bReg) {

                  Point3d ptFirstBar = ( CompoVert[0][0].ptPApp + CompoVert[0][1].ptPApp + CompoVert[0][2].ptPApp) / 3 ;
                  Point3d ptSecondBar = ( CompoVert[1][0].ptPApp + CompoVert[1][1].ptPApp + CompoVert[1][2].ptPApp) / 3 ;
                  Vector3d vtDiagBar = ptSecondBar - ptFirstBar ;
                  vtDiagBar.Normalize() ;
                  int nSum = 0 ;
                  for ( int nVec = 0 ; nVec < 3 ; ++ nVec) {
                     if ( CompoVert[0][nVec].vtVec * vtDiagBar > 0.7)
                        ++ nSum ;
                     else if ( CompoVert[0][nVec].vtVec * vtDiagBar < -0.7)
                        -- nSum;
                     if ( CompoVert[1][nVec].vtVec * vtDiagBar < -0.7)
                        ++ nSum ;
                     else if ( CompoVert[1][nVec].vtVec * vtDiagBar > 0.7)
                        -- nSum ;
                  }

                  if ( nSum < - 2)
                     bMatOnSlice = true ;
                  if ( nSum <= 2) {
                     Vector3d vtVoxCentre ;
                     if ( nAdjVox3[nCount] == - 1)
                        vtVoxCentre = X_AX ;
                     else if ( nAdjVox3[nCount] == 1)
                        vtVoxCentre = - X_AX ;
                     else if ( nAdjVox3[nCount] == - 2)
                        vtVoxCentre = Y_AX ;
                     else if ( nAdjVox3[nCount] == 2) 
                        vtVoxCentre = - Y_AX ;
                     else if ( nAdjVox3[nCount] == - 3)
                        vtVoxCentre = Z_AX ;
                     else if ( nAdjVox3[nCount] == 3)
                        vtVoxCentre = - Z_AX ;

                     int nPlusNum = 0 ;
                     int nMinusNum = 0 ;
                     for ( int nVec = 0 ; nVec < 3 ; ++ nVec) {
                        if ( CompoVert[0][nVec].vtVec * vtVoxCentre > EPS_SMALL)
                           ++ nPlusNum ;
                        else if ( CompoVert[0][nVec].vtVec * vtVoxCentre < - EPS_SMALL)
                           -- nMinusNum ;
                        if ( CompoVert[1][nVec].vtVec * vtVoxCentre > EPS_SMALL)
                           ++ nPlusNum ;
                        else if ( CompoVert[1][nVec].vtVec * vtVoxCentre < - EPS_SMALL)
                           -- nMinusNum ;
                     }

                     if ( nPlusNum >= nMinusNum)
                        bMatOnSlice = true ;
                  }
                  else
                     ;
               }
               // Conservo l'informazione per i voxel successivi
               if ( GetVoxNFromIJK( nIJKSl[0], nIJKSl[1], nIJKSl[2], nSliceN)) {
                  if ( abs(nAdjVox3[nCount]) == 1) {
                     if ( nSlBlockN == nBlock)
                        SliceYZ.emplace( nSliceN, bMatOnSlice) ;
                     else {
                        while ( m_SliceFlag.test_and_set( memory_order_acquire))
                           m_SliceFlag.wait( true, memory_order_relaxed) ;
                        m_SliceYZ[nSlBlockN].emplace( nSliceN, bMatOnSlice) ;
                        m_SliceFlag.clear( memory_order_release) ;
                        m_SliceFlag.notify_one() ;
                     }
                  }
                  else if ( abs(nAdjVox3[nCount]) == 2) {
                     if ( nSlBlockN == nBlock)
                        SliceXZ.emplace( nSliceN, bMatOnSlice) ;
                     else {
                        while ( m_SliceFlag.test_and_set( memory_order_acquire))
                           m_SliceFlag.wait( true, memory_order_relaxed) ;
                        m_SliceXZ[nSlBlockN].emplace( nSliceN, bMatOnSlice) ;
                        m_SliceFlag.clear( memory_order_release) ;
                        m_SliceFlag.notify_one() ;
                     }
                  }
                  else if ( abs(nAdjVox3[nCount]) == 3) {
                     if ( nSlBlockN == nBlock)
                        SliceXY.emplace(nSliceN, bMatOnSlice) ;
                     else {
                        while ( m_SliceFlag.test_and_set( memory_order_acquire))
                           m_SliceFlag.wait( true, memory_order_relaxed) ;
                        m_SliceXY[nSlBlockN].emplace( nSliceN, bMatOnSlice) ;
                        m_SliceFlag.clear( memory_order_release) ;
                        m_SliceFlag.notify_one() ;
                     }
                  }
               }

               bool bNewTopology = ( bDefTopology ? ! bMatOnSlice : bMatOnSlice) ;

              // Si passa alla seconda topologia
               if ( bNewTopology) {
                 // 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
                     nVertComp[nC - 1] = Cases3Plus[nRotCase][1][nC] ;
                    // 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) ;
                  }
               }
            }

           // Configurazione 6
            else if ( nAllConfig[nIndex] == 6) {
              // Test sulla topologia
               bool bDefTopology = false ;
               bool bMatOnSlice = false ;
               int nCount = 0 ;
               while ( nIndexConfig6[nCount] != nIndex)
                  ++ nCount ;
              // Vedo se la topologia è definita: altrimenti la ricalcolo
               int nIJKSl[3] = { ( nAdjVox6[nCount] != 1 ? i : i + 1),
                                 ( nAdjVox6[nCount] != 2 ? j : j + 1),
                                 ( nAdjVox6[nCount] != 3 ? k : k + 1)} ;
               int nSliceN ;
               int nSlBlockN ;
               if ( GetVoxNFromIJK( nIJKSl[0], nIJKSl[1], nIJKSl[2], nSliceN)) {
                  int nSlBlockIJK[3] ;
                  GetVoxelBlockIJK( nIJKSl, nSlBlockIJK) ;
                  if ( abs( nAdjVox6[nCount]) == 1) {
                     auto it = SliceYZ.find( nSliceN) ;
                     if ( it != SliceYZ.end()) { 
                        bMatOnSlice = it->second ;
                        bDefTopology = true ;
                     }
                     if ( GetBlockNFromIJK( nSlBlockIJK, nSlBlockN)) {
                        auto it = m_SliceYZ[nSlBlockN].find( nSliceN) ;
                        if ( it != m_SliceYZ[nSlBlockN].end()) {
                           bMatOnSlice = it->second ;
                           bDefTopology = true ;
                        }
                     }
                  }
                  else if ( abs( nAdjVox6[nCount]) == 2) {
                     auto it = SliceXZ.find( nSliceN) ;
                     if ( it != SliceXZ.end()) { 
                        bMatOnSlice = it->second ;
                        bDefTopology = true ;
                     }
                     if ( GetBlockNFromIJK( nSlBlockIJK, nSlBlockN)) {
                        auto it = m_SliceXZ[nSlBlockN].find( nSliceN) ;
                        if ( it != m_SliceXZ[nSlBlockN].end()) {
                           bMatOnSlice = it->second ;
                           bDefTopology = true ;
                        }
                     }
                  }
                  else if ( abs( nAdjVox6[nCount]) == 3) {
                     auto it = SliceXY.find( nSliceN) ;
                     if ( it != SliceXY.end()) { 
                        bMatOnSlice = it->second ;
                        bDefTopology = true ;
                     }
                     if ( GetBlockNFromIJK( nSlBlockIJK, nSlBlockN)) {
                        auto it = m_SliceXY[nSlBlockN].find( nSliceN) ;
                        if ( it != m_SliceXY[nSlBlockN].end()) {
                           bMatOnSlice = it->second ;
                           bDefTopology = true ;
                        }
                     }
                  }
               }
              // Topologia indefinita: la calcolo
               if ( ! bDefTopology  &&  bReg) {
                  // Test sulla topologia
                  Point3d ptFirstBar = ( CompoVert[0][0].ptPApp + CompoVert[0][1].ptPApp + 
                                         CompoVert[0][2].ptPApp + CompoVert[0][3].ptPApp) / 4 ;
                  Point3d ptSecondBar = ( CompoVert[1][0].ptPApp + CompoVert[1][1].ptPApp + 
                                          CompoVert[1][2].ptPApp) / 3 ;
                  Vector3d vtDiagBar = ptSecondBar - ptFirstBar ;
                  vtDiagBar.Normalize() ;
                  int nSum = 0 ;
                  for ( int nVec = 0 ; nVec < 4 ; ++ nVec) {
                     if ( CompoVert[0][nVec].vtVec * vtDiagBar > 0.7)
                        ++ nSum ;
                     else if ( CompoVert[0][nVec].vtVec * vtDiagBar < - 0.7)
                        -- nSum ;
                  }
                  for ( int nVec = 0 ; nVec < 3 ; ++ nVec) {
                     if ( CompoVert[1][nVec].vtVec * vtDiagBar < - 0.7)
                        ++ nSum ;
                     else if ( CompoVert[1][nVec].vtVec * vtDiagBar > 0.7)
                        -- nSum ;
                  }

                  if ( nSum < - 3) 
                     bMatOnSlice = true ;
                  else if ( nSum <= 3) {
                     Vector3d vtVoxCentre ;
                     if ( nAdjVox6[nCount] == -1)
                        vtVoxCentre = X_AX ;
                     else if ( nAdjVox6[nCount] == 1)
                        vtVoxCentre = - X_AX ;
                     else if ( nAdjVox6[nCount] == -2)
                        vtVoxCentre = Y_AX ;
                     else if ( nAdjVox6[nCount] == 2)
                        vtVoxCentre = - Y_AX ;
                     else if ( nAdjVox6[nCount] == -3)
                        vtVoxCentre = Z_AX ;
                     else if ( nAdjVox6[nCount] == 3)
                        vtVoxCentre = - Z_AX ;

                     int nPlusNum = 0 ;
                     int nMinusNum = 0 ;
                     for ( int nVec = 0 ; nVec < 4 ; ++ nVec) {
                        if ( CompoVert[0][nVec].vtVec * vtVoxCentre > EPS_SMALL)
                           ++ nPlusNum ;
                        else if ( CompoVert[0][nVec].vtVec * vtVoxCentre < -EPS_SMALL)
                           -- nMinusNum ;
                     }
                     for ( int nVec = 0 ; nVec < 3 ; ++ nVec) {
                        if ( CompoVert[1][nVec].vtVec * vtVoxCentre > EPS_SMALL)
                           ++ nPlusNum ;
                        else if ( CompoVert[1][nVec].vtVec * vtVoxCentre < -EPS_SMALL)
                           -- nMinusNum ;
                     }

                     if ( nPlusNum >= nMinusNum)
                        bMatOnSlice = true ;

                  }
                  else 
                     ;
               }
              // Conservo l'informazione
               if ( GetVoxNFromIJK( nIJKSl[0], nIJKSl[1], nIJKSl[2], nSliceN)) {
                  if ( abs(nAdjVox6[nCount]) == 1) {
                     if ( nSlBlockN == nBlock)
                        SliceYZ.emplace( nSliceN, bMatOnSlice) ;
                     else
                        m_SliceYZ[nSlBlockN].emplace( nSliceN, bMatOnSlice) ;
                  }
                  else if ( abs(nAdjVox6[nCount]) == 2) {
                     if ( nSlBlockN == nBlock)
                        SliceXZ.emplace( nSliceN, bMatOnSlice) ;
                     else
                        m_SliceXZ[nSlBlockN].emplace( nSliceN, bMatOnSlice) ;
                  }
                  else if (abs(nAdjVox6[nCount]) == 3) {
                     if ( nSlBlockN == nBlock)
                        SliceXY.emplace( nSliceN, bMatOnSlice) ;
                     else
                        m_SliceXY[nSlBlockN].emplace( nSliceN, bMatOnSlice) ;
                  }
               }
               bool bNewTopology = ( bDefTopology ? ! bMatOnSlice : bMatOnSlice) ;
              // Si deve passare alla seconda topologia
               if ( bNewTopology) {
                 // Ricerca del caso corrispondente della nuova topologia 
                  int nt = 0 ;
                  while ( nIndexVsIndex6[nt][0] != nIndex)
                     ++ nt ;
                  int nRotCase = nIndexVsIndex6[nt][1] ;
                 // Aggiorno numero di componenti
                  nComponents = Cases6Plus[nRotCase][1][0] ;
                 // Riaggiorno gli offsets
                  nExtTabOff = nComponents ;
                  nStdTabOff = 0 ;
                 // Modifico le matrici
                  for ( int nC = 1 ; nC <= nComponents ; ++ nC) {
                    // Numero vertici per componenti
                     nVertComp[nC - 1] = Cases6Plus[nRotCase][1][nC] ;
                    // Matrice dei vertici della base del fan
                     for ( int nFanVert = 0 ; nFanVert < nVertComp[nC - 1] ; ++ nFanVert)
                        CompoVert[nC - 1][nFanVert] = VecField[Cases6Plus[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[Cases6Plus[nRotCase][0][nStdTabOff + nTriVert+2]] ;
                        CompoTriVert[nC - 1][nTriVert+1] = VecField[Cases6Plus[nRotCase][0][nStdTabOff + nTriVert+1]] ;
                        CompoTriVert[nC - 1][nTriVert+2] =   VecField[Cases6Plus[nRotCase][0][nStdTabOff + nTriVert]] ;
                     }
                    // Aggiorno gli offsets per raggiungere i vertici della componente successiva.
                     nExtTabOff += nVertComp[nC - 1] ;
                     nStdTabOff += 3 * ( nVertComp[nC - 1] - 2) ;
                  }
               }
            }

           // Configurazione 7
            else if ( nAllConfig[nIndex] == 7) {
              // !!! Versione provvisoria, deve essere riveduta e semplificata !!!
              // Test sulla topologia
               bool bMatOnSliceYZ = false ;
               bool bMatOnSliceXZ = false ;
               bool bMatOnSliceXY = false ;
               bool bDefSliceYZ = false ;
               bool bDefSliceXZ = false ;
               bool bDefSliceXY = false ;
               int nCount = 0 ;
               while ( nIndexConfig7[nCount] != nIndex)
                  ++ nCount ;

              // Vedo se la topologia è definita: se sì uso l'informazione già posseduta,
              // altrimenti devo calcolare la topologia

              // Numerazione delle facce del voxel:  0: YZ-   1: XZ-   2: YZ+   3: XZ+   4: XY-   5: XY+
              // - 1 faccia non in gioco, 0 faccia in gioco ma con topologia non definita,
              //  1 faccia in gioco con topologia definita
               int nFace[6] = { -1, -1, -1, -1, -1, -1} ;
              // Facce in gioco
               int nCurFaceYZ = 0 ;
               int nCurFaceXZ = 1 ;
               int nCurFaceXY = 4 ;
               
               int nIJKSlYZ[3] = { i, j, k} ;
               if ( nAdjVox7[nCount][0] == 1) {
                  ++ nIJKSlYZ[0] ;
                  nCurFaceYZ = 2 ;
               }
               int nIJKSlXZ[3] = { i, j, k} ;
               if ( nAdjVox7[nCount][1] == 2) {
                  ++ nIJKSlXZ[1] ;
                  nCurFaceXZ = 3 ;
               }
               int nIJKSlXY[3] = { i, j, k} ;
               if ( nAdjVox7[nCount][2] == 3) {
                  ++ nIJKSlXY[2] ;
                  nCurFaceXY = 5 ;
               }
              // Assegno l'indice zero alle facce in gioco. L'indice zero è per le facce in gioco 
              // su cui non si è certi se ci sia o meno materiale.
               nFace[nCurFaceYZ] = 0 ;
               nFace[nCurFaceXZ] = 0 ;
               nFace[nCurFaceXY] = 0 ;

              // Per ogni faccia in gioco vediamo se c'è materiale: se sì assegnamo l'indice 1
               int nSlYZN = -1 ;
               int nSlYZBlockN =-1 ;
               if ( GetVoxNFromIJK( nIJKSlYZ[0], nIJKSlYZ[1], nIJKSlYZ[2], nSlYZN)) {
                 // Slice interna al blocco
                  auto it = SliceYZ.find( nSlYZN) ;
                 // Topologia definita 
                  if ( it != SliceYZ.end()) { 
                     bMatOnSliceYZ = it->second ;
                     bDefSliceYZ = true ;
                     nFace[nCurFaceYZ] = 1 ;
                  }
                 // Slice sulla frontiera
                  int nSlBlockIJK[3] ;
                  GetVoxelBlockIJK( nIJKSlYZ, nSlBlockIJK) ; 
                  if ( GetBlockNFromIJK( nSlBlockIJK, nSlYZBlockN)) {
                     auto it = m_SliceYZ[nSlYZBlockN].find( nSlYZN) ;
                    // Topologia definita 
                     if ( it != m_SliceYZ[nSlYZBlockN].end()) {
                        bMatOnSliceYZ = it->second ; 
                        bDefSliceYZ = true ;
                        nFace[nCurFaceYZ] = 1 ;
                     }
                  }
               }
               int nSlXZN = -1 ;
               int nSlXZBlockN =-1 ;
               if ( GetVoxNFromIJK( nIJKSlXZ[0], nIJKSlXZ[1], nIJKSlXZ[2], nSlXZN)) {
                 // Slice interna al blocco
                  auto it = SliceXZ.find( nSlXZN) ;
                 // Topologia definita 
                  if ( it != SliceXZ.end()) { 
                     bMatOnSliceXZ = it->second ;
                     bDefSliceXZ = true ;
                     nFace[nCurFaceXZ] = 1 ;
                  }
                 // Slice sulla frontiera
                  int nSlBlockIJK[3] ;
                  GetVoxelBlockIJK( nIJKSlXZ, nSlBlockIJK) ;
                  if ( GetBlockNFromIJK( nSlBlockIJK, nSlXZBlockN)) {
                     auto it = m_SliceYZ[nSlXZBlockN].find( nSlXZN) ;
                    // Topologia definita 
                     if ( it != m_SliceYZ[nSlXZBlockN].end()) {
                        bMatOnSliceXZ = it->second ; 
                        bDefSliceXZ = true ;
                        nFace[nCurFaceXZ] = 1 ;
                     }
                  }
               }
               int nSlXYN =-1 ;
               int nSlXYBlockN =-1 ;
               if ( GetVoxNFromIJK( nIJKSlXY[0], nIJKSlXY[1], nIJKSlXY[2], nSlXYN)) {
                 // Slice interna al blocco
                  auto it = SliceXY.find( nSlXYN) ;
                 // Topologia definita
                  if ( it != SliceXY.end()) { 
                     bMatOnSliceXY = it->second ;
                     bDefSliceXY = true ;
                     nFace[nCurFaceXY] = 1 ;
                  }
                 // Slice sulla frontiera
                  int nSlBlockIJK[3] ;
                  GetVoxelBlockIJK( nIJKSlXY, nSlBlockIJK) ;
                  if ( GetBlockNFromIJK( nSlBlockIJK, nSlXYBlockN)) {
                     auto it = m_SliceYZ[nSlXYBlockN].find( nSlXYN) ;
                    // Topologia definita
                     if ( it != m_SliceYZ[nSlXYBlockN].end()) {
                        bMatOnSliceXY = it->second ; 
                        bDefSliceXY = true ;
                        nFace[nCurFaceXY] = 1 ;
                     }
                  }
               }

              // Numerazione delle facce del voxel 0: YZ-   1: XZ-   2: YZ+   3: XZ+   4: XY-   5: XY+
              // Gli spigoli sono ordinati in senso antiorario dal punto di vista di un osservatore esterno al voxel
               static int nSliceEdges[6][4] = { {  3,  8,  7, 11}, {  0,  9,  4,  8}, {  1, 10,  5,  9}, {  2, 11,  6,  10},
                                                {  0,  3,  2,  1}, {  4,  5,  6,  7} } ;
              // Le facce in gioco hano indici nCurFaceXY, nCurFaceXZ e nCurFaceYZ. Se per almeno una faccia non si sa se c'è materiale,
              // ovvero almeno una faccia ha indice 0, e il campo vettoriale normale è regolare verifico la presenza di materiale.
               if ( ( nFace[nCurFaceXY] == 0  ||  nFace[nCurFaceXZ]  == 0  ||  nFace[nCurFaceYZ] == 0)  &&  bReg) {
                 // Ciclo sulle facce 
                  for ( int nFaceN = 0 ; nFaceN < 6 ; ++ nFaceN) {
                    // Faccia da analizzare
                     if ( nFace[nFaceN] == 0) {
                       // Dalla tabella determino le due componenti connesse che si appoggiano alla faccia
                        int nFaceCompo1 = 0 ;
                        int nFaceCompo2 = 0 ;
                       // Ciclo sulle componenti connesse del voxel
                        for ( int nCurComp = 1 ; nCurComp <= 3 ; ++ nCurComp) {
                           int nMatchEdge = 0 ;
                         // Ciclo sugli edge della componente
                           for ( int nEdge = 0 ; nEdge < 3 ; ++ nEdge) {
                             // Ciclo sugli edge della faccia
                              for ( int nFaceEdge = 0 ; nFaceEdge < 4 ; ++ nFaceEdge) {
                                // Edge della componente coincide con quello della faccia
                                 if ( nSliceEdges[nFaceN][nFaceEdge] == TriangleTableEn[nIndex][1][3 + nCurComp - 1 + nEdge]) {
                                    ++ nMatchEdge ;
                                    break ;
                                 }
                              }
                             // La componente corrente ha due edge sulla faccia, 
                             // interrompiamo la ricerca di edge corrispondenti
                              if ( nMatchEdge == 2)  
                                 break ;
                           } 
                          // Abbiamo trovato una nuova componente con due edge sulla faccia 
                           if ( nMatchEdge == 2) {
                              if ( nFaceCompo1 == 0)
                                 nFaceCompo1 = nCurComp ;
                              else 
                                 nFaceCompo2 = nCurComp ;
                           }
                          // Se le componenti con due edge sulla faccia sono due, interrompiamo la ricerca
                           if ( nFaceCompo2 != 0)
                              break ;
                        }
                        // Valuto la topologia   VecField[EdgeIndex] edgeIndex = 0, 1, ..., 11
                        double dDotSum = 0 ;
                        for ( int nV = 0 ; nV < 3 ; ++ nV) {
                           Vector3d vtV1 = VecField[TriangleTableEn[nIndex][1][3 + nFaceCompo1 - 1 + nV]].vtVec ;
                           Vector3d vtV2 = VecField[TriangleTableEn[nIndex][1][3 + nFaceCompo2 - 1 + nV]].vtVec ;
                           dDotSum += ( vtV1 * vtV2) ; 
                        }
                        if ( nFace[nFaceN] == 0  ||  nFace[nFaceN] == 2) 
                           bMatOnSliceXZ = dDotSum > 0. ;
                        else if ( nFace[nFaceN] == 1  ||  nFace[nFaceN] == 3) 
                           bMatOnSliceYZ = dDotSum > 0. ;
                        else 
                           bMatOnSliceXY = dDotSum > 0. ;
                     }
                  }
               }
              // Uso le informazioni per scegliere la topologia
               int nFaceWithMatNum = 0 ;
               if ( bMatOnSliceXZ  &&  ! bDefSliceXZ)
                  ++ nFaceWithMatNum ;
               if ( bMatOnSliceYZ  &&  ! bDefSliceYZ)
                  ++ nFaceWithMatNum ;
               if ( bMatOnSliceXY  &&  ! bDefSliceXY)
                  ++ nFaceWithMatNum ;
               if ( nFaceWithMatNum == 1) {
                  int nFaceCase = ( bMatOnSliceYZ ? 0 : ( bMatOnSliceXZ ? 1 : 2)) ;
                  // Aggiorno numero di componenti
                  nComponents = Cases7Plus[nCount][nFaceCase][1][0] ;
                  // Riaggiorno gli offsets
                  nExtTabOff = nComponents ;
                  nStdTabOff = 0 ;
                  // Modifico le matrici
                  for ( int nC = 1 ; nC <= nComponents ; ++ nC) {
                     // Numero vertici per componenti
                     nVertComp[nC - 1] = Cases7Plus[nCount][nFaceCase][1][nC] ;
                     // Matrice dei vertici della base del fan
                     for ( int nFanVert = 0 ; nFanVert < nVertComp[nC - 1] ; ++ nFanVert)
                        CompoVert[nC - 1][nFanVert] = VecField[Cases7Plus[nCount][nFaceCase][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[Cases7Plus[nCount][nFaceCase][0][nStdTabOff + nTriVert + 2]] ;
                        CompoTriVert[nC - 1][nTriVert + 1] = VecField[Cases7Plus[nCount][nFaceCase][0][nStdTabOff + nTriVert + 1]] ;
                        CompoTriVert[nC - 1][nTriVert + 2] = VecField[Cases7Plus[nCount][nFaceCase][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 ( nFaceWithMatNum == 2) {
                  ;
               }
              // Conservo l'informazione
               if ( GetVoxNFromIJK( nIJKSlYZ[0], nIJKSlYZ[1], nIJKSlYZ[2], nSlYZN)) {
                  if ( nSlYZBlockN == nBlock) 
                     SliceYZ.emplace( nSlYZN, bMatOnSliceYZ) ;
                  else 
                     m_SliceYZ[nSlYZBlockN].emplace( nSlYZN, bMatOnSliceYZ) ;
               }
               if ( GetVoxNFromIJK( nIJKSlXZ[0], nIJKSlXZ[1], nIJKSlXZ[2], nSlXZN)) {
                  if ( nSlXZBlockN == nBlock) 
                     SliceXZ.emplace( nSlXZN, bMatOnSliceXZ) ;
                  else 
                     m_SliceXZ[nSlXZBlockN].emplace( nSlXZN, bMatOnSliceXZ) ;
               } 
               if ( GetVoxNFromIJK( nIJKSlXY[0], nIJKSlXY[1], nIJKSlXY[2], nSlXYN)) {
                  if ( nSlXYBlockN == nBlock)
                     SliceXY.emplace( nSlXYN, bMatOnSliceXY) ;
                  else 
                     m_SliceXY[nSlXYBlockN].emplace( nSlXYN, bMatOnSliceXY) ;
               }
            }

           // Configurazione 10
            else if ( nAllConfig[nIndex] == 10) {
               // Test sulla topologia
               bool bDefStTopology = false ;
               bool bMatOnSlice = false ;
               int nCount = 0 ;
               while ( nIndexConfig10[nCount] != nIndex)
                  ++ nCount ;
              // Vedo se la topologia è definita: se sì uso l'informazione già posseduta,
              // altrimenti devo calcolare la topologia
               int nIJKSlSt[3] = { i, j, k} ;
               int nIJKSlEn[3] = { ( nAdjVox10[nCount] != 1 ? i : i + 1),
                                   ( nAdjVox10[nCount] != 2 ? j : j + 1),
                                   ( nAdjVox10[nCount] != 3 ? k : k + 1)} ;
               int nSliceStN, nSliceEnN ;
               int nSlBlockEnN ;
               if ( GetVoxNFromIJK( nIJKSlSt[0], nIJKSlSt[1], nIJKSlSt[2], nSliceStN)  && 
                    GetVoxNFromIJK( nIJKSlEn[0], nIJKSlEn[1], nIJKSlEn[2], nSliceEnN)) {
                  if ( abs( nAdjVox10[nCount]) == 1) {
                     auto itSt = SliceYZ.find( nSliceStN) ;
                     if ( itSt != SliceYZ.end()) {
                        bMatOnSlice = itSt->second ;
                        bDefStTopology = true ;
                     }
                  }
                  else if ( abs( nAdjVox10[nCount]) == 2) {
                     auto itSt = SliceXZ.find( nSliceStN) ;
                     if ( itSt != SliceXZ.end()) {
                        bMatOnSlice = itSt->second ;
                        bDefStTopology = true ;
                     }
                  }
                  else if ( abs( nAdjVox10[nCount]) == 3) {
                     auto itSt = SliceXY.find( nSliceStN) ;
                     if ( itSt != SliceXY.end()) {
                        bMatOnSlice = itSt->second ;
                        bDefStTopology = true ;
                     }
                  }
               }
              // La topologia non è definita, la calcolo
               if ( ! bDefStTopology  &&  bReg) {
                  // 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) ;
                  bMatOnSlice = ( ! bTest0  ||  ! bTest1) ;
               }
              // Conservo l'informazioe e la trasmetto al voxel successivo
               if ( GetVoxNFromIJK( nIJKSlSt[0], nIJKSlSt[1], nIJKSlSt[2], nSliceStN)  && 
                    GetVoxNFromIJK( nIJKSlEn[0], nIJKSlEn[1], nIJKSlEn[2], nSliceEnN)) {
                  if ( abs( nAdjVox6[nCount]) == 1) {
                     if ( GetBlockNFromIJK( nIJKSlEn, nSlBlockEnN)) {
                        auto it = m_SliceYZ[nSlBlockEnN].find( nSliceEnN) ;
                        if ( it != m_SliceYZ[nSlBlockEnN].end()) {
                           if ( it->second != bMatOnSlice)
                              m_BlockToUpdate[nSlBlockEnN] = true ;
                           it->second = bMatOnSlice ;
                        }
                        else
                           m_SliceYZ[nSlBlockEnN].emplace( nSliceEnN, bMatOnSlice) ;
                     }
                  }
                  else if ( abs( nAdjVox6[nCount]) == 2) {
                     if ( GetBlockNFromIJK( nIJKSlEn, nSlBlockEnN)) {
                        auto it = m_SliceXZ[nSlBlockEnN].find( nSliceEnN) ;
                        if ( it != m_SliceXZ[nSlBlockEnN].end()) {
                           if ( it->second != bMatOnSlice)
                              m_BlockToUpdate[nSlBlockEnN] = true ;
                           it->second = bMatOnSlice ;
                        }
                        else
                           m_SliceXZ[nSlBlockEnN].emplace( nSliceEnN, bMatOnSlice) ;
                     }
                  }
                  else if ( abs( nAdjVox6[nCount]) == 3) {
                     if ( GetBlockNFromIJK( nIJKSlEn, nSlBlockEnN)) {
                        auto it = m_SliceXY[nSlBlockEnN].find( nSliceEnN) ;
                        if ( it != m_SliceXY[nSlBlockEnN].end()) {
                           if ( it->second != bMatOnSlice)
                              m_BlockToUpdate[nSlBlockEnN] = true ;
                           it->second = bMatOnSlice ;
                        }
                        else
                           m_SliceXY[nSlBlockEnN].emplace( nSliceEnN, bMatOnSlice) ;
                     }
                  }
               }
               bool bNewTopology = ( bDefStTopology ? ! bMatOnSlice : bMatOnSlice) ;
               // Si passa alla seconda topologia
               if ( bNewTopology) {
                  // 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) ;
                  }
               }
            }

            Voxel VoxConf ;
            VoxConf.nNumComp = 0 ;

            SmoothTriaStruct VoxSmoothTria ;
            int nVoxSmootSizePrev = int( VoxSmoothTria.vTria.size()) ;
            int nVoxSmootSize = nVoxSmootSizePrev ;

           // Ciclo sulle componenti
            for ( int nComp = 0 ; nComp < nComponents ; ++ nComp) {

               int nFeatureType = NO_FEATURE ;
              // Se i componenti sono regolari valuto le normali per stabilire se eseguire ExtMC o MC 
               if ( bReg) 
                  nFeatureType = TestOnNormal( CompoVert[nComp], nVertComp[nComp]) ;

              // Extended MC
               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].ptPApp ;
                  ptGravityCenter = ptGravityCenter / nVertComp[nComp] ;

                  Vector3d vtTrasf[12] ;
                  for ( int ni = 0 ; ni < nVertComp[nComp] ; ++ ni) 
                     vtTrasf[ni] = CompoVert[nComp][ni].ptPApp - ptGravityCenter ;

                 // Preparo le matrici per il sistema
                  SvdMatrix dMatrixN( nVertComp[nComp], 3) ;
                  SvdVector dKnownVector( nVertComp[nComp]) ;

                 // 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].vtVec ;
                  for ( int ni = 0 ; ni < nVertComp[nComp] ; ++ ni) {
                     int nj = ( ni + 1) % nVertComp[nComp] ;
                     if ( vtNorm[ni] * vtNorm[nj] > 0.92) {
                        Vector3d vtNI = ( 0.6 * vtNorm[ni] + 0.4 * vtNorm[nj]) ;
                        Vector3d vtNJ = ( 0.4 * vtNorm[ni] + 0.6 * vtNorm[nj]) ; 
                        vtNorm[ni] = vtNI ;
                        vtNorm[ni].Normalize() ;
                        vtNorm[nj] = vtNJ ;
                        vtNorm[nj].Normalize() ;
                        ++ ni ;
                     }
                  }

                 // 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 ;
                     dMatrixN( ni, 2) = vtNorm[ni].z ;
                     dKnownVector( ni) = vtNorm[ni] * vtTrasf[ni] ;
                  }

                 // calcolo SVD 
                  SvdDecomposer svd( dMatrixN, Eigen::ComputeThinU | Eigen::ComputeThinV) ;
                  const auto& dSingularValue = svd.singularValues() ;
                  const auto& dMatrixV = svd.matrixV() ;

                 // 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
                  auto 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 ;

                 // Gestisco i ribaltamenti dei triangoli
                  bool bExtConfirmed = true ;
                  Vector3d vtNullSpace( dMatrixV( 0, 2), dMatrixV( 1, 2), dMatrixV( 2, 2)) ;
                  if ( nFeatureType == EDGE  &&  vtNullSpace.Normalize()  &&  ! IsPointInsideVoxelApprox( i, j, k, ptSol)) {
                     double dBestDisp = 0 ;
                     double dMaxVal = - DBL_MAX ;
                     double dDisplH = 0.1 * m_dStep ;
                     double dDispStop = 0.5 * m_dStep ;
                     double dDisp = - 0.5 * m_dStep;
                     while ( dDisp < dDispStop + EPS_SMALL) {
                        double dCrVal = GetFanTotalNormalDeviation( ptSol + dDisp * vtNullSpace, CompoVert[nComp], nVertComp[nComp]) ;
                        if ( dCrVal > dMaxVal) {
                           dMaxVal = dCrVal ;
                           dBestDisp = dDisp ;
                        }
                        dDisp += dDisplH ;
                     }
                     ptSol += dBestDisp * vtNullSpace ;
                     int nNewAdjVoxIJK[3] ;
                     if ( GetPointVoxel( ptSol, nNewAdjVoxIJK[0], nNewAdjVoxIJK[1], nNewAdjVoxIJK[2])) {
                        int nNewAdjVoxIndex = CalcIndex( nNewAdjVoxIJK[0], nNewAdjVoxIJK[1], nNewAdjVoxIJK[2]) ;
                        if ( nNewAdjVoxIndex != 0) {
                           bool bOverTurning = false ;
                           for ( int m = 0 ; m < nVertComp[nComp] ; ++ m) {
                              int l = ( m + 1) % nVertComp[nComp] ;
                              Vector3d vtNlm = ( CompoVert[nComp][l].ptPApp - ptSol) ^ ( CompoVert[nComp][m].ptPApp - ptSol) ;
                              vtNlm.Normalize() ;
                              double dDotl = vtNlm * CompoVert[nComp][l].vtVec ;
                              double dDotm = vtNlm * CompoVert[nComp][m].vtVec ;
                              if ( ( dDotl < -0.01  &&  dDotm < -0.01)  ||  dDotl < -0.95  ||  dDotm < -0.95) {
                                 bOverTurning = true ;
                                 break ;
                              }
                           }
                           if ( bOverTurning) {
                              Point3d ptMinDiag( ( i * m_nDexVoxRatio + 0.5) * m_dStep,
                                                 ( j * m_nDexVoxRatio + 0.5) * m_dStep,
                                                 ( k * m_nDexVoxRatio + 0.5) * m_dStep) ;
                              Point3d ptMaxDiag( ( ( i + 1) * m_nDexVoxRatio + 0.5) * m_dStep,
                                                 ( ( j + 1) * m_nDexVoxRatio + 0.5) * m_dStep,
                                                 ( ( k + 1) * m_nDexVoxRatio + 0.5) * m_dStep) ;
                              double dNewU1, dNewU2 ;
                              if ( IntersLineBox( ptSol, vtNullSpace, ptMinDiag, ptMaxDiag, dNewU1, dNewU2)) {
                                 if ( dNewU1 > 0.) {
                                    ptSol += ( dNewU1 + EPS_SMALL * ( dNewU2 - dNewU1)) * vtNullSpace ;
                                 }
                                 else if ( dNewU2 < 0.) {
                                    ptSol += ( dNewU2 + EPS_SMALL * ( dNewU1 - dNewU2)) * vtNullSpace ;
                                 }
                              }
                              else {
                           
                                 if ( nVertComp[nComp] == 5) {                                    
                                    int nNumPar = 0 ;
                                    for ( int m = 0 ; m < nVertComp[nComp] - 1 ; ++ m) {
                                       for ( int l = m + 1 ; l < nVertComp[nComp] ; ++ l) {
                                          if ( AreSameVectorExact( CompoVert[nComp][m].vtVec, CompoVert[nComp][l].vtVec))
                                             ++ nNumPar ;
                                       }
                                    }
                                    bExtConfirmed = ( nNumPar == 3) ;
                                 }
                                 else if ( nVertComp[nComp] < 5) {
                                    bExtConfirmed = false ;
                                 }
                              }
                           }
                        }
                     }
                  }

                 // Gestione casi speciali: tutti i vettori quasi equiversi tranne uno
                  if ( bExtConfirmed  &&  ( nVertComp[nComp] == 4  ||  nVertComp[nComp] == 5)) {
                     int nParNum = 0 ;
                     for ( int m = 0 ; m < nVertComp[nComp] - 1 ; ++ m) {
                        for ( int l = m + 1 ; l < nVertComp[nComp] ; ++ l) {
                           if ( CompoVert[nComp][m].vtVec * CompoVert[nComp][l].vtVec > 0.9) {
                              ++ nParNum ;
                           }
                        }
                     }
                     if ( nParNum == ( nVertComp[nComp] - 1) * ( nVertComp[nComp] - 2) / 2) {
                        int nNotParIndex = - 1 ;      
                        for ( int m = 0 ; m < nVertComp[nComp] ; ++ m) {
                           int nNumDiff = 0 ;
                           for ( int l = 0 ; l < nVertComp[nComp] ; ++ l) {
                              if ( CompoVert[nComp][m].vtVec * CompoVert[nComp][l].vtVec <= 0.9) {
                                 ++ nNumDiff ;
                              }
                           }
                           if ( nNumDiff == nVertComp[nComp] - 1) {
                              nNotParIndex = m ;
                              break ;
                           }
                        }
                        int nPrevL = ( nNotParIndex > 0 ? nNotParIndex - 1 : nVertComp[nComp] - 1) ;
                        int nNextL = ( nNotParIndex < nVertComp[nComp] - 1 ? nNotParIndex + 1 : 0) ;
                        Vector3d vtPrevN = ( CompoVert[nComp][nNotParIndex].ptPApp - ptSol) ^ ( CompoVert[nComp][nPrevL].ptPApp - ptSol) ;
                        Vector3d vtNextN = ( CompoVert[nComp][nNextL].ptPApp - ptSol) ^ ( CompoVert[nComp][nNotParIndex].ptPApp - ptSol) ;
                        vtPrevN.Normalize() ;
                        vtNextN.Normalize() ;
                        if ( nNotParIndex >= 0  &&
                             CompoVert[nComp][nNotParIndex].vtVec * vtPrevN < - EPS_SMALL  &&
                             CompoVert[nComp][nNotParIndex].vtVec * vtNextN < - EPS_SMALL) {
                           bExtConfirmed = false ;
                        }
                     }
                  }

                 // ExtMC confermato
                  if ( bExtConfirmed) {
                     int tOldCompo = VoxConf.nNumComp ;
                     ++ VoxConf.nNumComp ;
                     ConComp& CurFan = VoxConf.Compo[tOldCompo] ;
                     CurFan.nVertNum = nVertComp[nComp] ;
                     for ( int nV = 0 ; nV < nVertComp[nComp] ; ++ nV) {
                        CurFan.CompVecField[nV] = CompoVert[nComp][nV] ;
                     }
                      for ( int nV = 0 ; nV < CurFan.nVertNum  &&  CurFan.nVertNum >= 2 ; ++ nV) {
                        int nW = ( nV + 1) % CurFan.nVertNum ;
                        if ( AreSamePointApprox( CurFan.CompVecField[nV].ptPApp, CurFan.CompVecField[nW].ptPApp)  &&
                             CurFan.CompVecField[nV].vtVec * CurFan.CompVecField[nW].vtVec > 0.9) {
                           Vector3d vtAvNorm = 0.5 * ( CurFan.CompVecField[nV].vtVec + CurFan.CompVecField[nW].vtVec) ;
                           vtAvNorm.Normalize() ;
                           CurFan.CompVecField[nV].ptPApp = 0.5 * ( CurFan.CompVecField[nV].ptPApp + CurFan.CompVecField[nW].ptPApp) ;
                           for ( int nL = nW ; nL < CurFan.nVertNum - 1 ; ++ nL) {
                              CurFan.CompVecField[nV].ptPApp = CurFan.CompVecField[nL + 1].ptPApp ;
                              CurFan.CompVecField[nL].vtVec = CurFan.CompVecField[nL + 1].vtVec ;
                              CurFan.CompVecField[nL].nPropIndex = CurFan.CompVecField[nL + 1].nPropIndex;
                           }
                           -- CurFan.nVertNum ;
                        }
                     }
                     CurFan.ptVert = ptSol ;
                     CurFan.vtNullSpace = vtNullSpace ;
                     CurFan.bInside = IsPointInsideVoxelApprox( i, j, k, ptSol, EPS_SMALL) ;
                     CurFan.bCorner = ( nFeatureType == CORNER) ;
                  }
                 // ExtMC non confermato
                  else
                     CreateSmoothTriangle( nIndex, nVertComp[nComp], CompoTriVert[nComp], true, VoxSmoothTria) ;
               }
              // Standard MC
               else if ( m_nShape != BOX) {
                  CreateSmoothTriangle( nIndex, nVertComp[nComp], CompoTriVert[nComp], false, VoxSmoothTria) ;
               }
            }
           // 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) ;
            }
           // Se nel voxel abbiamo trovato componenti smooth
           // aggiorno i contenitori
            nVoxSmootSize = int( VoxSmoothTria.vTria.size()) ;
            if ( nVoxSmootSize > nVoxSmootSizePrev) {
               VoxSmoothTria.i = i ;
               VoxSmoothTria.j = j ;
               VoxSmoothTria.k = k ;
               m_BlockSmoothTria[nBlock].emplace_back( VoxSmoothTria) ;
               nVoxSmootSizePrev = nVoxSmootSize ;
            }
         }
      }
   }

  // Se il solido è un parallelepipedo creo direttamente i triangoli grandi
   if ( m_nShape == BOX) {
      int nBlockIJK[3] ;
      GetBlockIJKFromN( nBlock, nBlockIJK) ;
     // Determino il primo nodo pieno della mappa
      int nFirstVoxI, nFirstVoxJ, nFirstVoxK ;  
      GetFirstVoxIJK( nFirstVoxI, nFirstVoxJ, nFirstVoxK) ;
     // Determino il primo nodo pieno della mappa   
      int nLastVoxI, nLastVoxJ, nLastVoxK ;  
      GetLastVoxIJK( nLastVoxI, nLastVoxJ, nLastVoxK) ; 
     // Costruisco i triangoli paralleli al piano YZ
      if ( nBlockIJK[0] == 0  ||  nBlockIJK[0] + 1 == m_nFracLin[0]) {  
        // Indici dei dexel corrispondenti al confine inferiore fra voxel piatti e feature 
         int nDexMinJ = ( nBlockIJK[1] == 0 ? m_nDexVoxRatio * ( nFirstVoxJ + 1) :
                                              m_nDexVoxRatio * nBlockIJK[1] * m_nVoxNumPerBlock) ;
         int nDexMinK = ( nBlockIJK[2] == 0 ? m_nDexVoxRatio * ( nFirstVoxK + 1) :
                                              m_nDexVoxRatio * nBlockIJK[2] * m_nVoxNumPerBlock) ;    
        // Indici dei dexel corrispondenti al confine superiore fra voxel piatti e feature 
         int nDexMaxJ = ( nBlockIJK[1] + 1 < int( m_nFracLin[1]) ? m_nDexVoxRatio * ( nBlockIJK[1] + 1) * m_nVoxNumPerBlock :
                                                                   m_nDexVoxRatio * nLastVoxJ) ;
         int nDexMaxK = ( nBlockIJK[2] + 1 < int( m_nFracLin[2]) ? m_nDexVoxRatio * ( nBlockIJK[2] + 1) * m_nVoxNumPerBlock :
                                                                   m_nDexVoxRatio * nLastVoxK) ; 
        // Determino coordinate minime e massime dei punti dei triangoli
         double dYMin = ( nDexMinJ + 0.5) * m_dStep ;
         double dZMin = ( nDexMinK + 0.5) * m_dStep ;           
         double dYMax = ( nDexMaxJ + 0.5) * m_dStep ;
         double dZMax = ( nDexMaxK + 0.5) * m_dStep ;                                                                                     
        // Determino colore dei triangoli       
         int nFaceInfGrade = ( ! m_Values[1].empty() ? m_Values[1][nDexMinK * m_nNx[1] + nDexMinJ][0].nToolMin : 0) ;
         int nFaceSupGrade = ( ! m_Values[1].empty() ? m_Values[1][nDexMinK * m_nNx[1] + nDexMinJ][0].nToolMax : 0) ;
        // Piano di coordinata x inferiore: versore normale rivolto come X-                            
         if ( nBlockIJK[0] == 0) {                  
            Triangle3dEx trTria1, trTria2 ;
            CreateBigTriangleYZ( dYMin, dYMax, dZMin, dZMax, m_dMinZ[1], false, nFaceInfGrade, trTria1, trTria2) ;  
            trTria1.ToGlob( m_MapFrame) ;
            trTria2.ToGlob( m_MapFrame) ;
            m_BlockBigTria[nBlock].emplace_back( trTria1) ;
            m_BlockBigTria[nBlock].emplace_back( trTria2) ;                 
         }
        // Piano di coordinata x superiore: versore normale rivolto come X+
         if ( nBlockIJK[0] + 1 == m_nFracLin[0]) {            
            Triangle3dEx trTria1, trTria2 ;
            CreateBigTriangleYZ( dYMin, dYMax, dZMin, dZMax, m_dMaxZ[1], true, nFaceSupGrade, trTria1, trTria2) ;                         
            trTria1.ToGlob( m_MapFrame) ;
            trTria2.ToGlob( m_MapFrame) ;
            m_BlockBigTria[nBlock].emplace_back( trTria1) ;
            m_BlockBigTria[nBlock].emplace_back( trTria2) ;         
         }                          
      }
     // Costruisco i triangoli paralleli al piano XZ
      if ( nBlockIJK[1] == 0  ||  nBlockIJK[1] + 1 == m_nFracLin[1]) {
        // Indici dei dexel corrispondenti al confine inferiore fra voxel piatti e feature   
         int nDexMinI = ( nBlockIJK[0] == 0 ? m_nDexVoxRatio * ( nFirstVoxI + 1) :
                                              m_nDexVoxRatio * nBlockIJK[0] * m_nVoxNumPerBlock) ;
         int nDexMinK = ( nBlockIJK[2] == 0 ? m_nDexVoxRatio * ( nFirstVoxK + 1) :
                                              m_nDexVoxRatio * nBlockIJK[2] * m_nVoxNumPerBlock) ;
        // Indici dei dexel corrispondenti al confine superiore fra voxel piatti e feature 
         int nDexMaxI = ( nBlockIJK[0] + 1 < int( m_nFracLin[0]) ?  m_nDexVoxRatio * ( nBlockIJK[0] + 1) * m_nVoxNumPerBlock :
                                                                    m_nDexVoxRatio * nLastVoxI) ;
         int nDexMaxK = ( nBlockIJK[2] + 1 < int( m_nFracLin[2]) ?  m_nDexVoxRatio * ( nBlockIJK[2] + 1) * m_nVoxNumPerBlock : 
                                                                    m_nDexVoxRatio * nLastVoxK) ; 
        // Determino coordinate minime e massime dei punti dei triangoli
         double dXMin = ( nDexMinI + 0.5) * m_dStep ;
         double dZMin = ( nDexMinK + 0.5) * m_dStep ;
         double dXMax = ( nDexMaxI + 0.5) * m_dStep ;
         double dZMax = ( nDexMaxK + 0.5) * m_dStep ;                                                     
        // Determino colore dei triangoli        
         int nFaceInfGrade = ( ! m_Values[2].empty() ? m_Values[2][nDexMinI * m_nNx[2] + nDexMinK][0].nToolMin : 0) ;
         int nFaceSupGrade = ( ! m_Values[2].empty() ? m_Values[2][nDexMinI * m_nNx[2] + nDexMinK][0].nToolMax : 0) ;              
        // Piano di coordinata y inferiore: versore normale rivolto come Y-       
         if ( nBlockIJK[1] == 0) {               
            Triangle3dEx trTria1, trTria2 ;  
            CreateBigTriangleXZ( dXMin, dXMax, dZMin, dZMax, m_dMinZ[2], false, nFaceInfGrade, trTria1, trTria2) ;      
            trTria1.ToGlob( m_MapFrame) ;
            trTria2.ToGlob( m_MapFrame) ;
            m_BlockBigTria[nBlock].emplace_back( trTria1) ;
            m_BlockBigTria[nBlock].emplace_back( trTria2) ;            
         }         
        // Piano di coordinata y superiore: versore normale rivolto come Y+  
         if ( nBlockIJK[1] + 1 == m_nFracLin[1]) {           
            Triangle3dEx trTria1, trTria2 ;  
            CreateBigTriangleXZ( dXMin, dXMax, dZMin, dZMax, m_dMaxZ[2], true, nFaceSupGrade, trTria1, trTria2) ;   
            trTria1.ToGlob( m_MapFrame) ;
            trTria2.ToGlob( m_MapFrame) ;
            m_BlockBigTria[nBlock].emplace_back( trTria1) ;
            m_BlockBigTria[nBlock].emplace_back( trTria2) ;                       
         }                  
      }
     // Costruisco i triangoli paralleli al piano XY
      if ( nBlockIJK[2] == 0  ||  nBlockIJK[2] + 1 == m_nFracLin[2]) {  
        // Indici dei dexel corrispondenti al confine inferiore fra voxel piatti e feature 
         int nDexMinI = ( nBlockIJK[0] == 0 ? m_nDexVoxRatio * ( nFirstVoxI + 1) : 
                                              m_nDexVoxRatio * nBlockIJK[0] * m_nVoxNumPerBlock) ; 
         int nDexMinJ = ( nBlockIJK[1] == 0 ? m_nDexVoxRatio * ( nFirstVoxJ + 1) : 
                                              m_nDexVoxRatio * nBlockIJK[1] * m_nVoxNumPerBlock) ;  
        // Indici dei dexel corrispondenti al confine superiore fra voxel piatti e feature 
         int nDexMaxI = ( nBlockIJK[0] + 1 < int( m_nFracLin[0]) ? m_nDexVoxRatio * ( nBlockIJK[0] + 1) * m_nVoxNumPerBlock :
                                                                   m_nDexVoxRatio * nLastVoxI) ; 
         int nDexMaxJ = ( nBlockIJK[1] + 1 < int( m_nFracLin[1]) ? m_nDexVoxRatio * ( nBlockIJK[1] + 1) * m_nVoxNumPerBlock :
                                                                   m_nDexVoxRatio * nLastVoxJ) ; 
        // Determino coordinate minime e massime dei punti dei triangoli
         double dXMin = ( nDexMinI + 0.5) * m_dStep ;
         double dYMin = ( nDexMinJ + 0.5) * m_dStep ;
         double dXMax = ( nDexMaxI + 0.5) * m_dStep ;
         double dYMax = ( nDexMaxJ + 0.5) * m_dStep ;            
        // Determino colore dei triangoli         
         int nFaceInfGrade = ( ! m_Values[0].empty() ? m_Values[0][nDexMinJ * m_nNx[0] + nDexMinI][0].nToolMin : 0) ;
         int nFaceSupGrade = ( ! m_Values[0].empty() ? m_Values[0][nDexMinJ * m_nNx[0] + nDexMinI][0].nToolMax : 0) ;     
        // Piano di coordinata Z inferiore: versore normale rivolto come Z-  
         if ( nBlockIJK[2] == 0) {          
            Triangle3dEx trTria1, trTria2 ;
            CreateBigTriangleXY( dXMin, dXMax, dYMin, dYMax, m_dMinZ[0], false, nFaceInfGrade, trTria1, trTria2) ;    
            trTria1.ToGlob( m_MapFrame) ;
            trTria2.ToGlob( m_MapFrame) ;
            m_BlockBigTria[nBlock].emplace_back( trTria1) ;
            m_BlockBigTria[nBlock].emplace_back( trTria2) ;     
         }
        // Piano di coordinata Z superiore: versore normale rivolto come Z+
         if ( nBlockIJK[2] + 1 == m_nFracLin[2]) {           
            Triangle3dEx trTria1, trTria2 ;
            CreateBigTriangleXY( dXMin, dXMax, dYMin, dYMax, m_dMaxZ[0], true, nFaceSupGrade, trTria1, trTria2) ;
            trTria1.ToGlob( m_MapFrame) ;
            trTria2.ToGlob( m_MapFrame) ;
            m_BlockBigTria[nBlock].emplace_back( trTria1) ;
            m_BlockBigTria[nBlock].emplace_back( trTria2) ;    
         }                   
      }         
   }
  // Processo i Voxel con possibile superficie piana
   else {
      ProcessVoxContXY( VoxContXYInf, nBlock, false) ;
      ProcessVoxContXY( VoxContXYSup, nBlock, true) ;
      ProcessVoxContYZ( VoxContYZInf, nBlock, false) ;
      ProcessVoxContYZ( VoxContYZSup, nBlock, true) ;
      ProcessVoxContXZ( VoxContXZInf, nBlock, false) ;
      ProcessVoxContXZ( VoxContXZSup, nBlock, true) ;
   }
   
   return true ;
}

//----------------------------------------------------------------------------
bool
VolZmap::RegulateFeaturesChain( vector<VoxelContainer>& 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) {
              // Vertice fuori dal suo voxel
               if ( ! CurVox.Compo[nComp].bInside) {
                  if (  IsPointInsideVoxelApprox( CurVox.i, CurVox.j, CurVox.k, CurVox.Compo[nComp].ptVert, 1.5 * EPS_SMALL))
                     continue ;
                 // Caso corner 
                  if ( CurVox.Compo[nComp].bCorner) {
                    // Cerco i primi vicini
                     INTVECTOR 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) ;  
                        INTVECTOR vVecNFCur = ( nF < nInnSizeF ? vNearFirst : vBordNearFirst) ;
                       // Cerco i secondi vicini
                        INTVECTOR 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) ; 
                           INTVECTOR 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 {
                     INTVECTOR vNearInn, vNearBord ;
                     FindAdjComp( vVecVox, nBlock, nVox, nComp, vNearInn, vNearBord) ;
                     int nSizeInn = int( vNearInn.size()) ;
                     int nSizeBord = int( vNearBord.size() );
                     if ( nSizeInn + nSizeBord == 6) {
                        const Voxel* pVoxSt = nullptr ;
                        const Voxel* pVoxEn = nullptr ;
                        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) ;
                        }
                        Point3d ptPCur = CurVox.Compo[nComp].ptVert ;
                        Point3d ptPSt ;
                        Point3d ptPEn ;
                        if ( nSizeInn == 6) {
                           ptPSt = pVoxSt->Compo[vNearInn[2]].ptVert ;
                           ptPEn = pVoxEn->Compo[vNearInn[5]].ptVert ; 
                        }
                        else if ( nSizeBord ==  6) {
                           ptPSt = pVoxSt->Compo[vNearBord[2]].ptVert ;
                           ptPEn = pVoxEn->Compo[vNearBord[5]].ptVert ; 
                        }
                        else {
                           ptPSt = pVoxSt->Compo[vNearInn[2]].ptVert ;
                           ptPEn = pVoxEn->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 ;
                        if ( dMidU < dCurU) {
                           ptNew = ptPEn + ( ptMid - ptPEn) * vtCurrEn * vtCurrEn ;
                        }
                        else {
                           ptNew = ptPSt + ( ptMid - ptPSt) * vtStCurr * vtStCurr ;
                        }
                        if ( abs( vtStCurr * vtStEn) > 0.95  &&  abs( vtStCurr * vtCurrEn) > 0.95  &&
                             abs( vtStEn * vtCurrEn) > 0.95) {
                           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) {
                  if (  IsPointInsideVoxelApprox( CurVox.i, CurVox.j, CurVox.k, CurVox.Compo[nComp].ptVert, 1.5 * EPS_SMALL))
                     continue ;
                 // Caso feature
                  if ( ! CurVox.Compo[nComp].bCorner) {
                     INTVECTOR vNearInn, vNearBord ;
                     FindAdjComp( vVecVox, nBlock, nVox, nComp, vNearInn, vNearBord) ;
                     int nSizeInn = int( vNearInn.size()) ;
                     int nSizeBord = int( vNearBord.size() );
                     if ( nSizeInn + nSizeBord == 6) {
                        const Voxel* pVoxSt = nullptr ; 
                        const Voxel* pVoxEn = nullptr ;
                        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) ;
                        }
                        Point3d ptPCur = CurVox.Compo[nComp].ptVert ;
                        Point3d ptPSt ;
                        Point3d ptPEn ;
                        if ( nSizeInn == 6) {
                           ptPSt = pVoxSt->Compo[vNearInn[2]].ptVert ;
                           ptPEn = pVoxEn->Compo[vNearInn[5]].ptVert ;
                        }
                        else if ( nSizeBord ==  6) {
                           ptPSt = pVoxSt->Compo[vNearBord[2]].ptVert ;
                           ptPEn = pVoxEn->Compo[vNearBord[5]].ptVert ;
                        }
                        else {
                           ptPSt = pVoxSt->Compo[vNearInn[2]].ptVert ;
                           ptPEn = pVoxEn->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 * m_nDexVoxRatio + 0.5) * m_dStep,
                                           ( CurVox.j * m_nDexVoxRatio + 0.5) * m_dStep,
                                           ( CurVox.k * m_nDexVoxRatio + 0.5) * m_dStep) ;
                        Point3d ptCubeSup( ( ( CurVox.i + 1) * m_nDexVoxRatio + 0.5) * m_dStep,
                                           ( ( CurVox.j + 1) * m_nDexVoxRatio + 0.5) * m_dStep,
                                           ( ( CurVox.k + 1) * m_nDexVoxRatio + 0.5) * m_dStep) ;
                        double dU1, dU2 ;
                        if ( 1 - abs( vtStCurr * vtCurrEn ) < EPS_ZERO  &&  
                             IntersLineBox( ptPLine, vtDLine, ptCubeInf, ptCubeSup, dU1, dU2)) {
                           double dU = ( dU1 + dU2) / 2 ;
                           ptNew = ptPLine + dU * vtDLine ;
                        }
                        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) 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) {
      int tOldSize = int( m_BlockSharpTria[nBlock].size()) ;
      m_BlockSharpTria[nBlock].resize( tOldSize + 1) ;
      m_BlockSharpTria[nBlock][tOldSize].i = it->second.i ;
      m_BlockSharpTria[nBlock][tOldSize].j = it->second.j ;
      m_BlockSharpTria[nBlock][tOldSize].k = it->second.k ;
     // Ciclo sulle componenti connesse del voxel
      for ( int nComp = 0 ; nComp < it->second.nNumComp ; ++ nComp) {
          m_BlockSharpTria[nBlock][tOldSize].ptCompoVert.emplace_back( it->second.Compo[nComp].ptVert) ;
          m_BlockSharpTria[nBlock][tOldSize].ptCompoVert.back().ToGlob( m_MapFrame) ;
          int tOldCompNum = int( m_BlockSharpTria[nBlock][tOldSize].vCompoTria.size()) ;
          m_BlockSharpTria[nBlock][tOldSize].vCompoTria.resize( tOldCompNum + 1) ;
          m_BlockSharpTria[nBlock][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) ;
           // Setto i flag a false
            CurrTri.SetEdgeFlag( 0, false) ;
            CurrTri.SetEdgeFlag( 1, false) ;
            CurrTri.SetEdgeFlag( 2, false) ;
           // Valido il triangolo
            if ( ! CurrTri.Validate( true))
               CurrTri.SetVertexNorm( 0, 0.5 * ( CurrTri.GetVertexNorm( 1) + CurrTri.GetVertexNorm( 2))) ;
            m_BlockSharpTria[nBlock][tOldSize].vCompoTria[tOldCompNum].emplace_back( CurrTri) ;
            m_BlockSharpTria[nBlock][tOldSize].vCompoTria[tOldCompNum].back().ToGlob( m_MapFrame) ;
            int nTri = int( m_BlockSharpTria[nBlock][tOldSize].vbFlipped[tOldCompNum].size()) ;
            m_BlockSharpTria[nBlock][tOldSize].vbFlipped[tOldCompNum].resize( nTri + 1) ;
            m_BlockSharpTria[nBlock][tOldSize].vbFlipped[tOldCompNum][nTri] = false ;
         }
      }
   }

   
  // Pulisco il contenitore dei voxel di frontiera
   m_InterBlockOriginalSharpTria[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( m_BlockSharpTria[nBlock].size()) ;
      m_BlockSharpTria[nBlock].resize( nOldSizeInn + 1) ;
      m_BlockSharpTria[nBlock][nOldSizeInn].i = nVoxIJK[0] ;
      m_BlockSharpTria[nBlock][nOldSizeInn].j = nVoxIJK[1] ;
      m_BlockSharpTria[nBlock][nOldSizeInn].k = nVoxIJK[2] ;
     // Ridimensiono il contenitore dei triangoli di frontiera
      int nOldSizeBor = int( m_InterBlockOriginalSharpTria[nBlock].size()) ;
      m_InterBlockOriginalSharpTria[nBlock].resize( nOldSizeBor + 1) ;
      m_InterBlockOriginalSharpTria[nBlock][nOldSizeBor].i = nVoxIJK[0] ;
      m_InterBlockOriginalSharpTria[nBlock][nOldSizeBor].j = nVoxIJK[1] ;
      m_InterBlockOriginalSharpTria[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) ;
           // Setto i flag a false
            CurrTri.SetEdgeFlag( 0, false) ;
            CurrTri.SetEdgeFlag( 1, false) ;
            CurrTri.SetEdgeFlag( 2, false) ;
           // Valido il triangolo
            if ( ! CurrTri.Validate( true))
               CurrTri.SetVertexNorm( 0, 0.5 * ( CurrTri.GetVertexNorm( 1) + CurrTri.GetVertexNorm( 2)));
           // Smisto i triangoli fra di frontiera e interni
            bool bTriOnBorder = IsTriangleOnBorder( CurrTri, nLimits, nVoxIJK) ;
           // Triangolo di frontiera
            if ( bTriOnBorder) {   
               if ( bNewCompBor) {
                  Point3d ptVert = itVox->second.Compo[nComp].ptVert ;
                  ptVert.ToGlob( m_MapFrame) ;
                  m_InterBlockOriginalSharpTria[nBlock][nOldSizeBor].ptCompoVert.emplace_back( ptVert) ;
                  int tOldComp = int( m_InterBlockOriginalSharpTria[nBlock][nOldSizeBor].vCompoTria.size()) ;
                  m_InterBlockOriginalSharpTria[nBlock][nOldSizeBor].vCompoTria.resize( tOldComp + 1) ;
                  m_InterBlockOriginalSharpTria[nBlock][nOldSizeBor].vbFlipped.resize( tOldComp + 1) ;
                  bNewCompBor = false ;
               }
               int tCurrSz = int( m_InterBlockOriginalSharpTria[nBlock][nOldSizeBor].vCompoTria.size()) ; 
               m_InterBlockOriginalSharpTria[nBlock][nOldSizeBor].vCompoTria[tCurrSz - 1].emplace_back( CurrTri) ;
               m_InterBlockOriginalSharpTria[nBlock][nOldSizeBor].vCompoTria[tCurrSz - 1].back().ToGlob( m_MapFrame) ;
               int tTri =int( m_InterBlockOriginalSharpTria[nBlock][nOldSizeBor].vbFlipped[tCurrSz - 1].size()) ;
               m_InterBlockOriginalSharpTria[nBlock][nOldSizeBor].vbFlipped[tCurrSz - 1].resize( tTri + 1) ;
               m_InterBlockOriginalSharpTria[nBlock][nOldSizeBor].vbFlipped[tCurrSz - 1][tTri] = false ;
            }
           // Triangolo interno
            else {
               if ( bNewCompInn) {
                  Point3d ptVert =  itVox->second.Compo[nComp].ptVert ;
                  ptVert.ToGlob( m_MapFrame) ;
                  m_BlockSharpTria[nBlock][nOldSizeInn].ptCompoVert.emplace_back( ptVert) ;
                  int tOldComp = int( m_BlockSharpTria[nBlock][nOldSizeInn].vCompoTria.size()) ;
                  m_BlockSharpTria[nBlock][nOldSizeInn].vCompoTria.resize( tOldComp + 1) ;
                  m_BlockSharpTria[nBlock][nOldSizeInn].vbFlipped.resize( tOldComp + 1) ;
                  bNewCompInn = false ;
               }
               int tCurrSz = int( m_BlockSharpTria[nBlock][nOldSizeInn].vCompoTria.size()) ;
               m_BlockSharpTria[nBlock][nOldSizeInn].vCompoTria[tCurrSz - 1].emplace_back( CurrTri) ;
               m_BlockSharpTria[nBlock][nOldSizeInn].vCompoTria[tCurrSz - 1].back().ToGlob( m_MapFrame) ;
               int tTri = int( m_BlockSharpTria[nBlock][nOldSizeInn].vbFlipped[tCurrSz - 1].size()) ;
               m_BlockSharpTria[nBlock][nOldSizeInn].vbFlipped[tCurrSz - 1].resize( tTri + 1) ;
               m_BlockSharpTria[nBlock][nOldSizeInn].vbFlipped[tCurrSz - 1][tTri] = false ;
            }
         }
      }
   }

   return true ;
}

//----------------------------------------------------------------------------
bool
VolZmap::CreateSmoothTriangle(int nIndex, int nVertNum, AppliedVector TriVert[], bool bWasSharp, SmoothTriaStruct& VoxSmoothTria) const
{
   TRIA3DEXVECTOR vTria ;
  // Costruzione dei triangoli
   for ( int TriIndex = 0 ; TriIndex < (nVertNum - 2) * 3 ; TriIndex += 3) {
     // Il triangolo è pronto
      Triangle3dEx CurrentTriangle ;
      CurrentTriangle.Set( TriVert[TriIndex].ptPApp,
                           TriVert[TriIndex + 1].ptPApp,
                           TriVert[TriIndex + 2].ptPApp) ;
     // Setto il numero di utensile (conta solo positivo, nullo o negativo)
      int nTool0 = Clamp( TriVert[TriIndex].nPropIndex, - 1, 1) ;
      int nTool1 = Clamp( TriVert[TriIndex + 1].nPropIndex, - 1, 1) ;
      int nTool2 = Clamp( TriVert[TriIndex + 2].nPropIndex, - 1, 1) ;
      int nGrade = 0 ;
      if ( nTool0 == -1  ||  nTool1 == -1  ||  nTool2 == -1)
         nGrade = -1 ;
      else if ( nTool0 == 1  ||  nTool1 == 1  ||  nTool2 == 1)
         nGrade = 1 ;
      CurrentTriangle.SetGrade( nGrade) ;
      CurrentTriangle.SetAttrib( 0, nTool0) ;
      CurrentTriangle.SetAttrib( 1, nTool1) ;
      CurrentTriangle.SetAttrib( 2, nTool2) ;
     // Valido il triangolo e setto le normali del campo vettoriale ai corrispondenti vertici
      if ( CurrentTriangle.Validate(  true)) {
         double dCosAngThreshold = ( bWasSharp ? 0.7 : 0.5) ;
         for ( int nV = 0 ; nV < 3 ; ++ nV) {
            const Vector3d& vtVertNorm = TriVert[TriIndex + nV].vtVec ;
            if ( CurrentTriangle.GetN() * vtVertNorm > dCosAngThreshold)
               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) {
      int nFst1 = - 1 ;
      int nFst2 = - 1 ;
      int nSec1 = - 1 ;
      int nSec2 = - 1 ;
      for ( int nF = 0 ; nF < 3 ; ++ nF) {
         for ( int nS = 0 ; nS < 3 ; ++ nS) {
            if ( AreSamePointExact( vTria[0].GetP( nF), vTria[1].GetP( nS))) {
               if ( nFst1 < 0)
                  nFst1 = nF ;
               else
                  nFst2 = nF ;
               if ( nSec1 < 0)
                  nSec1 = nS ;
               else
                  nSec2 = nS ;
            }
         }
      }
      int nFst3 = -1 ;
      int nSec3 = -1 ;
      for ( int n = 0 ; n < 3 ; ++n) {
         if ( n != nFst1  &&  n != nFst2)
            nFst3 = n ;
         if ( n != nSec1  &&  n != nSec2)
            nSec3 = n ;
      }
      if ( vTria[0].GetAttrib( nFst1) != vTria[0].GetAttrib( nFst2)  && 
           vTria[0].GetAttrib( nFst3) == vTria[1].GetAttrib( nSec3)) {
         vTria[0].SetP( nFst1, vTria[1].GetP( nSec3)) ;
         vTria[0].SetAttrib( nFst1, vTria[1].GetAttrib( nSec3)) ;
         vTria[0].SetVertexNorm( nFst1, vTria[1].GetVertexNorm( nSec3)) ;
         int nFirstGrade = 0 ;
         if ( vTria[0].GetAttrib(0) == -1  ||  vTria[0].GetAttrib(1) == -1  ||  vTria[0].GetAttrib(2) == -1)
            nFirstGrade = -1 ;
         else if ( vTria[0].GetAttrib(0) == 1  ||  vTria[0].GetAttrib(1) == 1  ||  vTria[0].GetAttrib(2) == 1)
            nFirstGrade = 1 ;
         vTria[0].SetGrade( nFirstGrade) ;
         vTria[0].Validate( true) ;
         vTria[1].SetP( nSec2, vTria[0].GetP( nFst3)) ;
         vTria[1].SetAttrib( nSec2, vTria[0].GetAttrib( nFst3)) ;
         vTria[1].SetVertexNorm( nSec2, vTria[0].GetVertexNorm( nFst3)) ;
         int nSecondGrade = 0 ;
         if ( vTria[1].GetAttrib(0) == -1  ||  vTria[1].GetAttrib(1) == -1  ||  vTria[1].GetAttrib(2) == -1)
            nSecondGrade = -1 ;
         else if ( vTria[1].GetAttrib(0) == 1  ||  vTria[1].GetAttrib(1) == 1  ||  vTria[1].GetAttrib(2) == 1)
            nSecondGrade = 1 ;
         vTria[1].SetGrade( nSecondGrade) ;
         vTria[1].Validate( true) ;
      }
      else if ( vTria[0].GetAttrib(nFst1) == vTria[0].GetAttrib(nFst2)  &&
                vTria[0].GetAttrib(nFst3) != vTria[1].GetAttrib(nSec3)) {
         int nFirstGrade = 0 ;
         if ( vTria[0].GetAttrib(0) == -1  ||  vTria[0].GetAttrib(1) == -1  ||  vTria[0].GetAttrib(2) == -1)
            nFirstGrade = -1 ;
         else if ( vTria[0].GetAttrib(0) == 1  ||  vTria[0].GetAttrib(1) == 1  ||  vTria[0].GetAttrib(2) == 1)
            nFirstGrade = 1 ;
         vTria[0].SetGrade( nFirstGrade) ;
         int nSecondGrade = 0 ;
         if ( vTria[1].GetAttrib(0) == -1  ||  vTria[1].GetAttrib(1) == -1  ||  vTria[1].GetAttrib(2) == -1)
            nSecondGrade = -1 ;
         else if ( vTria[1].GetAttrib(0) == 1  ||  vTria[1].GetAttrib(1) == 1  ||  vTria[1].GetAttrib(2) == 1)
            nSecondGrade = 1 ;
         vTria[1].SetGrade( nSecondGrade) ;              
      }
   }
   for ( int nT = 0 ; nT < int( vTria.size()) ; ++nT)
      VoxSmoothTria.vTria.emplace_back( vTria[nT]) ;
   return true ;
}

//----------------------------------------------------------------------------
// Esegue il flipping dei triangoli contenuti nel TriHolder,
// bGraph indica se chiamata per grafica o per calcolo di profondità.
bool
VolZmap::FlipEdgesII( int nBlock) const
{
  // Numero di voxel in cui si presentano sharp feature
   int nVoxelNum = int( m_BlockSharpTria[nBlock].size()) ;

  // Ciclo sui voxel con sharp feature
   for ( int n1 = 0 ; n1 < nVoxelNum ; ++ n1) {
      SharpTriaStruct& SharpTria1 = m_BlockSharpTria[nBlock][n1] ;
      for ( int n2 = n1 ; n2 < nVoxelNum ; ++ n2) {
         SharpTriaStruct& SharpTria2 = m_BlockSharpTria[nBlock][n2] ;
        // Se non adiacenti o coincidenti vado oltre
         if ( abs( SharpTria2.i - SharpTria1.i) > 1  ||
              abs( SharpTria2.j - SharpTria1.j) > 1  ||
              abs( SharpTria2.k - SharpTria1.k) > 1 )
            continue ;  
        // Ciclo sulle componenti connesse del primo voxel
         int nNumCompo1 = int( SharpTria1.ptCompoVert.size()) ;
         for ( int nCompo1 = 0 ; nCompo1 < nNumCompo1 ; ++ nCompo1) {
            TRIA3DEXVECTOR& vTria1 = SharpTria1.vCompoTria[nCompo1] ;
           // Numero di triangoli della componente connessa 
            int nTriNum1 = int( vTria1.size()) ;
           // Ciclo sulle componenti connesse del secondo voxel
            int nNumCompo2 = int( SharpTria2.ptCompoVert.size()) ;
            int nCompo2 = ( n1 == n2 ? nCompo1 + 1 : 0) ;
            for ( ; nCompo2 < nNumCompo2 ; ++ nCompo2) {
               TRIA3DEXVECTOR& vTria2 = SharpTria2.vCompoTria[nCompo2] ;
              // Numero di triangoli della componente connessa 
               int nTriNum2 = int( vTria2.size()) ;
               for ( int nTri1 = 0 ; nTri1 < nTriNum1 ; ++ nTri1) { 
                  bool bModified = false ;
                  for ( int nTri2 = 0 ; nTri2 < nTriNum2 ; ++ nTri2) {
                    // Punti che devono essere in comune fra i due triangoli 
                     const Point3d& ptP11 = vTria1[nTri1].GetP( 1) ;
                     const Point3d& ptP12 = vTria1[nTri1].GetP( 2) ;
                     const Point3d& ptP21 = vTria2[nTri2].GetP( 1) ;
                     const Point3d& ptP22 = vTria2[nTri2].GetP( 2) ;
                    // I triangoli sono da flippare
                     if ( AreSamePointEpsilon( ptP11, ptP22, EPS_ZERO)  &&
                          AreSamePointEpsilon( ptP12, ptP21, EPS_ZERO)  && 
                          ! ( SharpTria1.vbFlipped[nCompo1][nTri1]  ||
                              SharpTria2.vbFlipped[nCompo2][nTri2])) {
                       // Assegno l'array dei punti di contorno
                        Point3d vPnt[4] ;
                        vPnt[0] = vTria1[nTri1].GetP( 1) ;
                        vPnt[1] = SharpTria1.ptCompoVert[nCompo1] ;
                        vPnt[2] = vTria1[nTri1].GetP( 2) ;
                        vPnt[3] = SharpTria2.ptCompoVert[nCompo2] ;
                       // Valuto se i triangoli giacciono su un piano
                        PolygonPlane Polygon ;
                        for ( int i = 0 ; i < 4 ; ++ i)
                           Polygon.AddPoint( vPnt[i]) ;
                        Plane3d plPlane ;
                        bool bOnPlane = Polygon.GetPlane( plPlane) ;
                        for ( int i = 0 ; i < 4 && bOnPlane ; ++ i)
                           bOnPlane = PointInPlaneApprox( vPnt[i], plPlane) ;
                       // Se sono su un piano controllo se avviene inversione
                        bool bInv = false ;
                        if ( bOnPlane) {
                           Triangle3d trT1 = vTria1[nTri1] ;
                           Triangle3d trT2 = vTria2[nTri2] ;
                           int nVert1 = -1, nVert2 = -1 ;
                          // Determino gli indici dei punti sharp-feature
                           for ( int nP = 0 ; nP < 3 ; ++ nP) {  
                              if ( nP != 1  &&  nP != 2)
                                 nVert1 = nP ;
                              if ( nP != 2  &&  nP != 1)
                                 nVert2 = nP ;
                           }
                           trT1.SetP( 1, trT2.GetP( nVert2)) ;
                           trT2.SetP( 1, trT1.GetP( nVert1)) ;
                           trT1.Validate( true) ;
                           trT2.Validate( true) ;
                           bInv = ( trT1.GetN() * trT2.GetN() < 0) ;
                        }
                       // Se non vi è inversione eseguo il flipping
                        if ( ! bInv) {
                          // Vertice condiviso fra nTri1 e quello del suo fan
                           int nCol1 = vTria1[nTri1].GetAttrib( 2) ;
                           int nCol2 = vTria2[nTri2].GetAttrib( 2) ; 
                          // Modifico i punti e gli indici 
                           vTria1[nTri1].SetP( 1, SharpTria2.ptCompoVert[nCompo2]) ;
                           vTria2[nTri2].SetP( 1, SharpTria1.ptCompoVert[nCompo1]) ;
                           vTria1[nTri1].SetGrade( nCol1) ;
                           vTria2[nTri2].SetGrade( nCol2) ;
                          // Setto le normali 
                           vTria1[nTri1].SetVertexNorm( 0, V_NULL) ;
                           vTria1[nTri1].SetVertexNorm( 1, V_NULL) ;
                           vTria2[nTri2].SetVertexNorm( 0, V_NULL) ;
                           vTria2[nTri2].SetVertexNorm( 1, V_NULL) ;
                          // Valido i triangoli
                           vTria1[nTri1].Validate( true) ;
                           vTria2[nTri2].Validate( true) ;
                          // Setto i bFlags
                           if ( ! AreSameVectorApprox( vTria1[nTri1].GetN(), vTria2[nTri2].GetN())) {
                              vTria1[nTri1].SetEdgeFlag( 0, true) ;                            
                           }
                          // Setto i triangoli come flippati
                           SharpTria1.vbFlipped[nCompo1][nTri1] = true ;
                           SharpTria2.vbFlipped[nCompo2][nTri2] = true ;
                          // Avvenuto flipping
                           bModified = true ;
                           break ;
                        }
                        else {
                          // In questo caso i due triangoli sono necessariamente su un piano, 
                          // quindi hanno normali concordi. A entrambi assegno il colore del vertice con
                          // normale più concorde a quella del triangolo
                           double dDotVec[4] ;
                           dDotVec[0] = vTria1[nTri1].GetN() * vTria1[nTri1].GetVertexNorm( 1) ;
                           dDotVec[1] = vTria2[nTri2].GetN() * vTria2[nTri2].GetVertexNorm( 2) ;
                           dDotVec[2] = vTria1[nTri1].GetN() * vTria1[nTri1].GetVertexNorm( 2) ;
                           dDotVec[3] = vTria2[nTri2].GetN() * vTria2[nTri2].GetVertexNorm( 1) ;
                          // Cerco il massimo dei prodotti scalari
                           int nMaxPos = 0 ;
                           double dMaxDot = - 1 ;
                           for ( int nPos = 0 ; nPos < 4  &&  dMaxDot < 1 ; ++ nPos) {
                              if ( dDotVec[nPos] > dMaxDot) {
                                 dMaxDot = dDotVec[nPos] ; 
                                 nMaxPos = nPos ;
                              }
                           }
                          // Trovo il colore associato al vertice di massimo prodotto scalare
                           int nCol = -1 ;
                           switch ( nMaxPos) {
                           case 0 :
                              nCol = vTria1[nTri1].GetAttrib( 1) ; 
                              break ;
                           case 1 :
                              nCol = vTria2[nTri2].GetAttrib( 2) ; 
                              break ;
                           case 2 :
                              nCol = vTria1[nTri1].GetAttrib( 2) ; 
                              break ;
                           case 3 :
                              nCol = vTria2[nTri2].GetAttrib( 1) ;
                              break ;
                           }
                          // Assegno il colore ai triangoli
                           vTria1[nTri1].SetGrade( nCol) ;
                           vTria2[nTri2].SetGrade( nCol) ;
                          // Setto i flag dei vertici
                           double dDist02 = ( vTria1[nTri1].GetP( 0) - vTria1[nTri1].GetP( 2)).Len() +
                                            ( vTria2[nTri2].GetP( 1) - vTria2[nTri2].GetP( 0)).Len() ;
                           double dDist20 = ( vTria1[nTri1].GetP( 1) - vTria1[nTri1].GetP( 0)).Len() +
                                            ( vTria2[nTri2].GetP( 0) - vTria2[nTri2].GetP( 2)).Len() ;
                           if ( dDist02 > dDist20) {
                              vTria1[nTri1].SetEdgeFlag( 0, true) ;
                              vTria2[nTri2].SetEdgeFlag( 2, true) ;
                           }
                           else {
                              vTria1[nTri1].SetEdgeFlag( 2, true) ;
                              vTria2[nTri2].SetEdgeFlag( 0, true) ;
                           }
                        }
                     }
                  }
                  if ( bModified)
                     break ;
               }
            }
         }
      }
   }

  // Ciclo sui voxel con sharp feature
   for ( int nVox = 0 ; nVox < nVoxelNum ; ++nVox) {
      SharpTriaStruct& SharpTria = m_BlockSharpTria[nBlock][nVox] ;
     // Ciclo sulle componenti connesse del voxel
      int nCompNum = int( SharpTria.ptCompoVert.size()) ;
      for ( int nComp = 0 ; nComp < nCompNum ; ++ nComp) {
         TRIA3DEXVECTOR& vTria = SharpTria.vCompoTria[nComp] ;
        // Ciclo sui triangoli della componente connessa
         int nTriaNum = int( vTria.size()) ;
         for ( int nTria = 0 ; nTria < nTriaNum ; ++ nTria) {
           // Lunghezza del segmento opposto al vertice
            double dDist = Dist( vTria[nTria].GetP( 1), vTria[nTria].GetP( 2)) ;
           // Se è quasi nullo
            if ( dDist < 2 * EPS_SMALL  &&  ! SharpTria.vbFlipped[nComp][nTria]) {
              // Cerco triangolo adiacente con normale più vicina
               int nPrevTria = ( nTria - 1 >= 0 ? nTria - 1 : nTriaNum - 1) ;
               int nNextTria = ( nTria + 1 < nTriaNum ? nTria + 1 : 0) ;
               double dDotPrev = vTria[nTria].GetN() * vTria[nPrevTria].GetN() ;
               double dDotNext = vTria[nTria].GetN() * vTria[nNextTria].GetN() ;
               if ( dDotPrev > dDotNext) {
                  double dAdjDist = Dist( vTria[nPrevTria].GetP( 1), vTria[nPrevTria].GetP( 2)) ;
                 // Se normale sufficientemente vicina e lato opposto al vertice abbastanza lungo assegno normale al vertice del triangolo corrente
                  if ( dDotPrev > 0.5  &&  dAdjDist > 2 * EPS_SMALL)
                     vTria[nTria].SetVertexNorm( 0, vTria[nPrevTria].GetVertexNorm( 0)) ;
               }
               else {
                  double dAdjDist = Dist( vTria[nNextTria].GetP( 1), vTria[nNextTria].GetP( 2)) ;
                 // Se normale sufficientemente vicina e lato opposto al vertice abbastanza lungo assegno normale al vertice del triangolo corrente
                  if ( dDotNext > 0.5  &&  dAdjDist > 2 * EPS_SMALL)
                     vTria[nTria].SetVertexNorm( 0, vTria[nNextTria].GetVertexNorm( 0)) ;
               }
            }
         }
      }
   }

  // Ciclo sui voxel con sharp feature
   for ( int nVox = 0 ; nVox < nVoxelNum ; ++ nVox) {
      SharpTriaStruct& SharpTria = m_BlockSharpTria[nBlock][nVox] ;
     // Ciclo sulle componenti connesse del voxel
      int nCompNum = int( SharpTria.ptCompoVert.size()) ;
      for ( int nComp = 0 ; nComp < nCompNum ; ++ nComp) {
         TRIA3DEXVECTOR& vTria = SharpTria.vCompoTria[nComp] ;
        // Ciclo sui triangoli della componente connessa
         int nTriaNum = int( vTria.size()) ;
         for ( int nTria = 0 ; nTria < nTriaNum ; ++ nTria) {
           // Se il triangolo ha partecipato a un flipping
            if ( vTria[nTria].GetVertexNorm( 0).IsSmall()  &&
                 vTria[nTria].GetVertexNorm( 1).IsSmall()) {  
               const Point3d& ptP0 = vTria[nTria].GetP( 0) ;
               const Point3d& ptP1 = vTria[nTria].GetP( 1) ;
               const Point3d& ptP2 = vTria[nTria].GetP( 2) ;
               int nPrevTria = ( nTria - 1 >= 0 ? nTria - 1 : nTriaNum - 1) ;
               const Point3d& ptPrevP0 = vTria[nPrevTria].GetP( 0) ;
               const Point3d& ptPrevP1 = vTria[nPrevTria].GetP( 1) ;
               if ( AreSamePointApprox( ptP0, ptPrevP0)  &&  AreSamePointApprox( ptP2, ptPrevP1)  &&
                    vTria[nPrevTria].GetVertexNorm( 0).IsNormalized())
                  vTria[nTria].SetVertexNorm( 0, vTria[nPrevTria].GetVertexNorm( 0)) ;
               else {
                  vTria[nTria].Validate( true) ;
                  vTria[nTria].SetVertexNorm( 1, V_NULL) ;
               }
              // Ciclo su voxel adiacenti
               bool bFound = false ;
               for ( int nAdjVox = 0 ; nAdjVox < nVoxelNum ; ++ nAdjVox) {
                  SharpTriaStruct& AdjSharpTria = m_BlockSharpTria[nBlock][nAdjVox] ;
                 // Se coincidenti o non adiacenti vado oltre
                  if ( nAdjVox == nVox  ||
                       abs( AdjSharpTria.i - SharpTria.i) > 1  ||
                       abs( AdjSharpTria.j - SharpTria.j) > 1  ||
                       abs( AdjSharpTria.k - SharpTria.k) > 1)
                     continue ;
                 // Ciclo sulle componenti connesse del voxel adiacente
                  int nAdjCompNum = int( AdjSharpTria.ptCompoVert.size()) ;
                  for ( int nAdjComp = 0 ; nAdjComp < nAdjCompNum ; ++ nAdjComp) {
                     TRIA3DEXVECTOR& vAdjTria = AdjSharpTria.vCompoTria[nAdjComp] ;
                    // Ciclo sui triangoli della componente connessa
                     int nAdjTriaNum = int( vAdjTria.size()) ;
                     for ( int nAdjTria = 0 ; nAdjTria < nAdjTriaNum ; ++ nAdjTria) {
                        const Point3d& ptAdjP0 = vAdjTria[nAdjTria].GetP( 0) ;
                        const Point3d& ptAdjP2 = vAdjTria[nAdjTria].GetP( 2) ;
                        if ( AreSamePointApprox( ptP1, ptAdjP0)  &&  AreSamePointApprox( ptP2, ptAdjP2)) {
                           if ( vAdjTria[nAdjTria].GetVertexNorm( 0).IsNormalized())
                              vTria[nTria].SetVertexNorm( 1, vAdjTria[nAdjTria].GetVertexNorm( 0)) ;
                           else 
                              vTria[nTria].Validate( true) ;
                           bFound = true ;
                           break ;
                        }  
                     }
                     if ( bFound)
                        break ;
                  }
                  if ( bFound)
                     break ;
               }        
               vTria[nTria].Validate( true) ;
            }
         }
      }
   }

   return true ;
}

//----------------------------------------------------------------------------
bool
VolZmap::FlipEdgesBB( void) const
{
  // Numero di blocchi
   int nBlocksNum = int( m_InterBlockSharpTria.size()) ;
  // ciclo sui blocchi 
   for ( int tFB = 0 ; tFB < nBlocksNum ; ++ tFB)  {
      int nFBijk[3] ;
      GetBlockIJKFromN( int( tFB), nFBijk) ;
      for ( int tLB = tFB ; tLB < nBlocksNum ; ++ tLB) {
         int nLBijk[3] ;
         GetBlockIJKFromN( int( tLB), nLBijk) ;
        // Se i blocchi non sono adiacenti salto l'iterazione
         if ( ! ( abs( nFBijk[0] - nLBijk[0]) <= 1  &&
                  abs( nFBijk[1] - nLBijk[1]) <= 1  &&
                  abs( nFBijk[2] - nLBijk[2]) <= 1))
            continue ;
        // Numero di voxel nei blocchi correnti
         int nVoxelNumFB = int( m_InterBlockSharpTria[tFB].size()) ;
         int nVoxelNumLB = int( m_InterBlockSharpTria[tLB].size()) ;
        // Ciclo sui voxel dei due blocchi
         for ( int tVFB = 0 ; tVFB < nVoxelNumFB ; ++ tVFB) {
            for ( int tVLB = 0 ; tVLB < nVoxelNumLB ; ++ tVLB) {
              // Se i voxel non sono adiacenti salto l'iterazione
               if ( ! ( abs( m_InterBlockSharpTria[tFB][tVFB].i - m_InterBlockSharpTria[tLB][tVLB].i) <= 1  &&
                        abs( m_InterBlockSharpTria[tFB][tVFB].j - m_InterBlockSharpTria[tLB][tVLB].j) <= 1  &&
                        abs( m_InterBlockSharpTria[tFB][tVFB].k - m_InterBlockSharpTria[tLB][tVLB].k) <= 1))
                  continue ;
              // Numero di componenti connesse dei voxel
               int nCompoVFBNum = int( m_InterBlockSharpTria[tFB][tVFB].ptCompoVert.size()) ;
               int nCompoVLBNum = int( m_InterBlockSharpTria[tLB][tVLB].ptCompoVert.size()) ;
              // Ciclo sulle componenti connesse
               for ( int tCmpF = 0 ; tCmpF < nCompoVFBNum ; ++ tCmpF) {
                  for ( int tCmpL = 0 ; tCmpL < nCompoVLBNum ; ++ tCmpL) {
                    // Numero di triangoli delle componenti connesse
                     int nTriFBNum = int( m_InterBlockSharpTria[tFB][tVFB].vCompoTria[tCmpF].size()) ;
                     int nTriLBNum = int( m_InterBlockSharpTria[tLB][tVLB].vCompoTria[tCmpL].size()) ;
                    // Ciclo sui triangoli
                     for ( int tTriFB = 0 ; tTriFB < nTriFBNum ; ++ tTriFB) {
                        bool bModified = false ;
                        for ( int tTriLB = 0 ; tTriLB < nTriLBNum ; ++ tTriLB) {
                          // Punti che devono essere in comune fra i due triangoli 
                           Point3d ptPF1 = m_InterBlockSharpTria[tFB][tVFB].vCompoTria[tCmpF][tTriFB].GetP( 1) ;
                           Point3d ptPF2 = m_InterBlockSharpTria[tFB][tVFB].vCompoTria[tCmpF][tTriFB].GetP( 2) ;
                           Point3d ptPL1 = m_InterBlockSharpTria[tLB][tVLB].vCompoTria[tCmpL][tTriLB].GetP( 1) ;
                           Point3d ptPL2 = m_InterBlockSharpTria[tLB][tVLB].vCompoTria[tCmpL][tTriLB].GetP( 2) ;
                          // Si deve operare la modifica dei triangoli
                           if ( AreSamePointEpsilon( ptPF1, ptPL2, EPS_ZERO)  &&
                                AreSamePointEpsilon( ptPF2, ptPL1, EPS_ZERO)  &&  
                                ! ( m_InterBlockSharpTria[tFB][tVFB].vbFlipped[tCmpF][tTriFB]  ||
                                    m_InterBlockSharpTria[tLB][tVLB].vbFlipped[tCmpL][tTriLB])) {
                             // Assegno l'array dei punti di contorno
                              Point3d vPnt[4] ;
                              vPnt[0] = m_InterBlockSharpTria[tFB][tVFB].vCompoTria[tCmpF][tTriFB].GetP( 1) ;
                              vPnt[1] = m_InterBlockSharpTria[tFB][tVFB].ptCompoVert[tCmpF] ;
                              vPnt[2] = m_InterBlockSharpTria[tFB][tVFB].vCompoTria[tCmpF][tTriFB].GetP( 2) ;
                              vPnt[3] = m_InterBlockSharpTria[tLB][tVLB].ptCompoVert[tCmpL] ;
                             // Valuto se i triangoli giacciono su un piano
                              PolygonPlane Polygon ;
                              for ( int i = 0 ; i < 4 ; ++ i)
                                 Polygon.AddPoint( vPnt[i]) ;
                              Plane3d plPlane ;
                              bool bOnPlane = Polygon.GetPlane( plPlane) ;
                              for ( int i = 0 ; i < 4 && bOnPlane ; ++ i)
                                 bOnPlane = PointInPlaneApprox( vPnt[i], plPlane) ;
                             // Se sono su un piano controllo se avviene inversione
                              bool bInv = false ;
                              if ( bOnPlane) {
                                 Triangle3dEx trTF = m_InterBlockSharpTria[tFB][tVFB].vCompoTria[tCmpF][tTriFB] ;
                                 Triangle3dEx trTL = m_InterBlockSharpTria[tLB][tVLB].vCompoTria[tCmpL][tTriLB] ;
                                 int nVertF = -1, nVertL =-1 ;
                                // Determino gli indici dei punti sharp-feature
                                 for ( int nP = 0 ; nP < 3 ; ++ nP) {  
                                    if ( nP != 1  &&  nP != 2)
                                       nVertF = nP ;
                                    if ( nP != 2  &&  nP != 1)
                                       nVertL = nP ;
                                 }
                                 trTF.SetP( 1, trTL.GetP( nVertL)) ;
                                 trTF.Validate( true) ;
                                 trTL.SetP( 1, trTF.GetP( nVertF)) ;
                                 trTL.Validate( true) ;
                                 bInv = ( trTF.GetN() * trTL.GetN() < 0) ;
                              }
                             // Se non vi è inversione eseguo il flipping
                              if ( ! bInv) {
                                // Vertice condiviso fra nTri1 e quello del suo fan
                                 int nColF = m_InterBlockSharpTria[tFB][tVFB].vCompoTria[tCmpF][tTriFB].GetAttrib( 2) ;
                                 int nColL = m_InterBlockSharpTria[tLB][tVLB].vCompoTria[tCmpL][tTriLB].GetAttrib( 2) ; 
                                // modifico punti e colori
                                 m_InterBlockSharpTria[tFB][tVFB].vCompoTria[tCmpF][tTriFB].SetP( 1, m_InterBlockSharpTria[tLB][tVLB].ptCompoVert[tCmpL]) ;
                                 m_InterBlockSharpTria[tLB][tVLB].vCompoTria[tCmpL][tTriLB].SetP( 1, m_InterBlockSharpTria[tFB][tVFB].ptCompoVert[tCmpF]) ;
                                 m_InterBlockSharpTria[tFB][tVFB].vCompoTria[tCmpF][tTriFB].SetGrade( nColF) ;
                                 m_InterBlockSharpTria[tLB][tVLB].vCompoTria[tCmpL][tTriLB].SetGrade( nColL) ;
                                // Valido i triangoli
                                 m_InterBlockSharpTria[tFB][tVFB].vCompoTria[tCmpF][tTriFB].Validate( true) ;
                                 m_InterBlockSharpTria[tLB][tVLB].vCompoTria[tCmpL][tTriLB].Validate( true) ;
                                // Setto le normali a ogni vertice                                   
                                 Vector3d vtNormF = m_InterBlockSharpTria[tFB][tVFB].vCompoTria[tCmpF][tTriFB].GetVertexNorm( 2) ;
                                 Vector3d vtNormL = m_InterBlockSharpTria[tLB][tVLB].vCompoTria[tCmpL][tTriLB].GetVertexNorm( 2) ;                                
                                 m_InterBlockSharpTria[tFB][tVFB].vCompoTria[tCmpF][tTriFB].SetVertexNorm( 1, vtNormF) ;
                                 m_InterBlockSharpTria[tFB][tVFB].vCompoTria[tCmpF][tTriFB].SetVertexNorm( 0, vtNormF) ;
                                 m_InterBlockSharpTria[tLB][tVLB].vCompoTria[tCmpL][tTriLB].SetVertexNorm( 1, vtNormL) ;
                                 m_InterBlockSharpTria[tLB][tVLB].vCompoTria[tCmpL][tTriLB].SetVertexNorm( 0, vtNormL) ;
                                // Valido i triangoli
                                 m_InterBlockSharpTria[tFB][tVFB].vCompoTria[tCmpF][tTriFB].Validate( true) ;
                                 m_InterBlockSharpTria[tLB][tVLB].vCompoTria[tCmpL][tTriLB].Validate( true) ;
                                // Setto i bFlags
                                 if ( ! AreSameVectorApprox( m_InterBlockSharpTria[tFB][tVFB].vCompoTria[tCmpF][tTriFB].GetN(),
                                                             m_InterBlockSharpTria[tLB][tVLB].vCompoTria[tCmpL][tTriLB].GetN())) {
                                    m_InterBlockSharpTria[tFB][tVFB].vCompoTria[tCmpF][tTriFB].SetEdgeFlag( 0, true) ;                                    
                                 }
                                // Setto i triangoli come flippati 
                                 m_InterBlockSharpTria[tFB][tVFB].vbFlipped[tCmpF][tTriFB] = true ; 
                                 m_InterBlockSharpTria[tLB][tVLB].vbFlipped[tCmpL][tTriLB] = true ;
                                 bModified = true ;
                                 break ;
                              }
                              else {
                                // In questo caso i due triangoli sono necessariamente su un piano, 
                                // quindi hanno normali concordi. A entrambi assegno il colore del vertice con
                                // normale più concorde a quella del triangolo
                                 double dDotVec[4] ;
                                 dDotVec[0] = m_InterBlockSharpTria[tFB][tVFB].vCompoTria[tCmpF][tTriFB].GetN() *
                                              m_InterBlockSharpTria[tFB][tVFB].vCompoTria[tCmpF][tTriFB].GetVertexNorm( 1) ;
                                 dDotVec[1] = m_InterBlockSharpTria[tLB][tVLB].vCompoTria[tCmpL][tTriLB].GetN() *
                                              m_InterBlockSharpTria[tLB][tVLB].vCompoTria[tCmpL][tTriLB].GetVertexNorm( 2) ; 
                                 dDotVec[2] = m_InterBlockSharpTria[tFB][tVFB].vCompoTria[tCmpF][tTriFB].GetN() *
                                              m_InterBlockSharpTria[tFB][tVFB].vCompoTria[tCmpF][tTriFB].GetVertexNorm( 2) ;
                                 dDotVec[3] = m_InterBlockSharpTria[tLB][tVLB].vCompoTria[tCmpL][tTriLB].GetN() *
                                              m_InterBlockSharpTria[tLB][tVLB].vCompoTria[tCmpL][tTriLB].GetVertexNorm( 1) ;
                                // Cerco il massimo dei prodotti scalari
                                 int nMaxPos = 0 ;
                                 double dMaxDot = - 1 ;
                                 for ( int nPos = 0 ; nPos < 4  &&  dMaxDot < 1 ; ++ nPos) {
                                    if ( dDotVec[nPos] > dMaxDot) {
                                       dMaxDot = dDotVec[nPos] ;
                                       nMaxPos = nPos ;
                                    }
                                 }
                                // Trovo il colore associato al vertice di massimo prodotto scalare
                                 int nCol = -1 ;
                                 switch ( nMaxPos) {
                                 case 0 :
                                    nCol = m_InterBlockSharpTria[tFB][tVFB].vCompoTria[tCmpF][tTriFB].GetAttrib( 1) ; 
                                    break ;
                                 case 1 :
                                    nCol = m_InterBlockSharpTria[tLB][tVLB].vCompoTria[tCmpL][tTriLB].GetAttrib( 2) ; 
                                    break ;
                                 case 2 :
                                    nCol = m_InterBlockSharpTria[tFB][tVFB].vCompoTria[tCmpF][tTriFB].GetAttrib( 2) ; 
                                    break ;
                                 case 3 :
                                    nCol = m_InterBlockSharpTria[tLB][tVLB].vCompoTria[tCmpL][tTriLB].GetAttrib( 1) ;
                                    break ;
                                 }
                                // Assegno il colore ai triangoli
                                 m_InterBlockSharpTria[tFB][tVFB].vCompoTria[tCmpF][tTriFB].SetGrade( nCol) ;
                                 m_InterBlockSharpTria[tLB][tVLB].vCompoTria[tCmpL][tTriLB].SetGrade( nCol) ;
                                // Setto i flag dei vertici
                                 double dDist02 = ( m_InterBlockSharpTria[tFB][tVFB].vCompoTria[tCmpF][tTriFB].GetP( 0) -
                                                    m_InterBlockSharpTria[tFB][tVFB].vCompoTria[tCmpF][tTriFB].GetP( 2)).Len() +
                                                  ( m_InterBlockSharpTria[tLB][tVLB].vCompoTria[tCmpL][tTriLB].GetP( 1) -
                                                    m_InterBlockSharpTria[tLB][tVLB].vCompoTria[tCmpL][tTriLB].GetP( 0)).Len() ;
                                 double dDist20 = ( m_InterBlockSharpTria[tFB][tVFB].vCompoTria[tCmpF][tTriFB].GetP( 1) -
                                                    m_InterBlockSharpTria[tFB][tVFB].vCompoTria[tCmpF][tTriFB].GetP( 0)).Len() +
                                                  ( m_InterBlockSharpTria[tLB][tVLB].vCompoTria[tCmpL][tTriLB].GetP( 0) -
                                                    m_InterBlockSharpTria[tLB][tVLB].vCompoTria[tCmpL][tTriLB].GetP( 2)).Len() ;
                                 if ( dDist02 > dDist20) {
                                    m_InterBlockSharpTria[tFB][tVFB].vCompoTria[tCmpF][tTriFB].SetEdgeFlag( 0, true) ;
                                    m_InterBlockSharpTria[tLB][tVLB].vCompoTria[tCmpL][tTriLB].SetEdgeFlag( 2, true) ;
                                 }
                                 else {
                                    m_InterBlockSharpTria[tFB][tVFB].vCompoTria[tCmpF][tTriFB].SetEdgeFlag( 2, true) ;
                                    m_InterBlockSharpTria[tLB][tVLB].vCompoTria[tCmpL][tTriLB].SetEdgeFlag( 0, true) ;
                                 }
                              }
                           }
                        }
                        if ( bModified)
                           break ;
                     }
                  }
               }
            }
         }
      }
   }
   return true ;
}

//----------------------------------------------------------------------------
bool
VolZmap::IsThereMat( int nI, int nJ, int nK) const
{
  // Trasformo gli indici della griglia voxel in quelli della griglia dexel
   nI *= m_nDexVoxRatio ; 
   nJ *= m_nDexVoxRatio ; 
   nK *= m_nDexVoxRatio ; 

  // Se l'indice è alla frontiera del reticolo non vi è materiale
   if ( nI <= - 1  ||  nI >= int( m_nNx[0])  ||
        nJ <= - 1  ||  nJ >= int( m_nNy[0])  ||
        nK <= - 1  ||  nK >= int( m_nNy[1]))
      return false ;

  // ciclo sulle griglie
   int nCount = 0 ;
   int nMinPos[3] = { -1, -1, -1} ;
   int nMinIndex[3] ;
   double dZ[3] ;
   int nDexSize[3] ;
   bool bInterOnNode[3] = { false, false, false} ;
   for ( int nGrid = 0 ; nGrid < int( m_nMapNum) ; ++ nGrid) {
     // assegnazione dati vertice dipendenti dalla griglia
      int nGrI = -1, nGrJ = -1 ;
      switch ( nGrid) {
      case 0 :
         nGrI = nI ;
         nGrJ = nJ ; 
         dZ[nGrid] = ( nK + 0.5) * m_dStep ;
         break ;
      case 1 :
         nGrI = nJ ;
         nGrJ = nK ;
         dZ[nGrid] = ( nI + 0.5) * m_dStep ;
         break ;
      case 2 :
         nGrI = nK ;
         nGrJ = nI ;
         dZ[nGrid] = ( nJ + 0.5) * m_dStep ;
         break ;
      }
     // verifica spillone su vertice
      double dMinDist = INFINITO ;
      int nIndex = 0 ;
      int nPos = nGrJ * m_nNx[nGrid] + nGrI ;
      nDexSize[nGrid] = int( m_Values[nGrid][nPos].size()) ;
     // scorro i sotto-interalli dello spillone
      while ( nIndex < nDexSize[nGrid]) {
        // distanza tra la "Z" attuale e il parametro minimo e massimo dell'intervallo nIndex-esimo
         double dDistInf = dZ[nGrid] - m_Values[nGrid][nPos][nIndex].dMin + 2 * EPS_SMALL ;
         double dDistSup = dZ[nGrid] - m_Values[nGrid][nPos][nIndex].dMax - 2 * EPS_SMALL ;
        // se "Z" attuale a cavallo tra queste due distanze...
         if ( dDistInf > 0.  &&  dDistSup < 0.) {
            nMinIndex[nGrid] = nIndex ;      // aggiorno l'indice
            ++ nCount ;
            bInterOnNode[nGrid] = true ;    // flag T per griglia zero
            break ;
         }
        // se "Z" attuale tutta sopra o tutta sotto
         else {
            double dDist = min( abs( dDistInf), abs( dDistSup)) ;
            if ( dDist < dMinDist) {
               nMinPos[nGrid] = nPos ;
               nMinIndex[nGrid] = nIndex ;
               dMinDist = dDist ;
            }
         }
        // sotto-intervallo successivo
         nIndex += 1 ;
      }
   } // fine ciclo sulle griglie

   if ( nCount == 3) // ... se interno a tutte e 3 le griglie, allora c'è materiale ...
      return true ;
   else if ( nCount == 2) { // ... se interno solo a 2 griglie ...
     // recupero la griglia sulla quale è esterno
      int nGrid = ( bInterOnNode[0] ? ( bInterOnNode[1] ? 2 : 1) : 0) ;
     // se tale griglia non ha sotto-intervalli allora non c'è materiale
      if ( nDexSize[nGrid] == 0)
         return false ;
     // se il valore in "Z" dello spillone è vicino al punto ( 1/10 dello step), aggiorno il
     // parametro minimo e massimo e considero la presenza di materiale
      if ( dZ[nGrid] > m_Values[nGrid][nMinPos[nGrid]][nMinIndex[nGrid]].dMin - 0.1 * m_dStep  &&
           dZ[nGrid] < m_Values[nGrid][nMinPos[nGrid]][nMinIndex[nGrid]].dMax + 0.1 * m_dStep) {
         double dDistInf = abs( dZ[nGrid] - m_Values[nGrid][nMinPos[nGrid]][nMinIndex[nGrid]].dMin) ;
         double dDistSup = abs( dZ[nGrid] - m_Values[nGrid][nMinPos[nGrid]][nMinIndex[nGrid]].dMax) ;
         if ( dDistInf < dDistSup)
            const_cast<double&>( m_Values[nGrid][nMinPos[nGrid]][nMinIndex[nGrid]].dMin) = dZ[nGrid] ;
         else
            const_cast<double&>( m_Values[nGrid][nMinPos[nGrid]][nMinIndex[nGrid]].dMax) = dZ[nGrid] ;
         return true ;
      }
      else
         return false ;
   }
   else // ... se invece interno a 1 o a nessuna delle griglie, allora non c'è materiale
      return false ;
}

//----------------------------------------------------------------------------
int
VolZmap::CalcIndex( int nI, int nJ, int nK) const
{
   int nIndex = 0 ;
   if ( IsThereMat( nI, nJ, nK))
      nIndex |= ( 1 << 0) ;
   if ( IsThereMat( nI + 1, nJ, nK))
      nIndex |= ( 1 << 1) ;
   if ( IsThereMat( nI + 1, nJ + 1, nK))
      nIndex |= ( 1 << 2) ;
   if ( IsThereMat( nI, nJ + 1, nK))
      nIndex |= ( 1 << 3) ;
   if ( IsThereMat( nI, nJ, nK + 1))
      nIndex |= ( 1 << 4) ;
   if ( IsThereMat( nI + 1, nJ, nK + 1))
      nIndex |= ( 1 << 5) ;
   if ( IsThereMat( nI + 1, nJ + 1, nK + 1))
      nIndex |= ( 1 << 6) ;
   if ( IsThereMat( nI, nJ + 1, nK + 1))
      nIndex |= ( 1 << 7) ;
   return nIndex ;
}

//----------------------------------------------------------------------------
bool
VolZmap::IntersPos( int nVec1[], int nVec2[], bool bFirstCorner, AppliedVector& vfField) const
{
   const double dEps = 2 * EPS_SMALL ;
   const double dEpsClamp = EPS_SMALL ; // EPS_SMALL ; 0

   bool bFound = false ;

   if ( nVec1[0] != nVec2[0]) {

      int nMinI = min( nVec1[0], nVec2[0]) * m_nDexVoxRatio ;
      int nMaxI = max( nVec1[0], nVec2[0]) * m_nDexVoxRatio ;

      double dMinX = ( nMinI + 0.5) * m_dStep ;
      double dMaxX = ( nMaxI + 0.5) * m_dStep ;

      vfField.ptPApp.y = ( nVec1[1] * m_nDexVoxRatio + 0.5) * m_dStep ;
      vfField.ptPApp.z = ( nVec1[2] * m_nDexVoxRatio + 0.5) * m_dStep ;

      int nDexel = ( nVec1[2] * m_nNx[1] + nVec1[1]) * m_nDexVoxRatio ;
      int nSize = int( m_Values[1][nDexel].size()) ;

      if ( bFirstCorner) {
         int n = nSize - 1 ;
         double dX = m_Values[1][nDexel][n].dMax ;
         while ( n >= 0  &&  dX > dMinX - dEps) {
            if ( dX < dMaxX + dEps) {
               vfField.ptPApp.x = Clamp( dX, dMinX + dEpsClamp, dMaxX - dEpsClamp) ;
               vfField.vtVec = m_Values[1][nDexel][n].vtMaxN ;
               vfField.nPropIndex = m_Values[1][nDexel][n].nToolMax ;
               bFound = true ;
               break ;
            }
            n -= 1 ;
            if ( n >= 0) 
               dX = m_Values[1][nDexel][n].dMax ;
         }
      }

      else {
         int n = 0 ;
         double dX = m_Values[1][nDexel][n].dMin ;
         while ( n < nSize  &&  dX < dMaxX + dEps) {
            if ( dX > dMinX - dEps) {
               vfField.ptPApp.x = Clamp( dX, dMinX + dEpsClamp, dMaxX - dEpsClamp) ;
               vfField.vtVec = m_Values[1][nDexel][n].vtMinN ;
               vfField.nPropIndex = m_Values[1][nDexel][n].nToolMin ;
               bFound = true ;
               break ;
            }
            n += 1 ;
            if ( n < nSize)
               dX = m_Values[1][nDexel][n].dMin ;
         }
      }

      if ( ! bFound)
         vfField.ptPApp.x = 0.5 * ( dMinX + dMaxX) ;
   }

   else if ( nVec1[1] != nVec2[1]) {

      int nMinJ = min( nVec1[1], nVec2[1]) * m_nDexVoxRatio ;
      int nMaxJ = max( nVec1[1], nVec2[1]) * m_nDexVoxRatio ;
      
      double dMinY = ( nMinJ + 0.5) * m_dStep ;
      double dMaxY = ( nMaxJ + 0.5) * m_dStep ;

      vfField.ptPApp.x = ( nVec1[0] * m_nDexVoxRatio + 0.5) * m_dStep ;
      vfField.ptPApp.z = ( nVec1[2] * m_nDexVoxRatio + 0.5) * m_dStep ;

      int nDexel = ( nVec1[0] * m_nNx[2] + nVec1[2]) * m_nDexVoxRatio ;
      int nSize = int( m_Values[2][nDexel].size()) ;

      if ( bFirstCorner) {
         int n = nSize - 1 ;
         double dY = m_Values[2][nDexel][n].dMax ;
         while ( n >= 0  &&  dY > dMinY - dEps) {
            if ( dY < dMaxY + dEps) {
               vfField.ptPApp.y = Clamp( dY, dMinY + dEpsClamp, dMaxY - dEpsClamp) ;
               vfField.vtVec = m_Values[2][nDexel][n].vtMaxN ;
               vfField.nPropIndex = m_Values[2][nDexel][n].nToolMax;
               bFound = true ;
               break ;
            }
            n -= 1 ;
            if ( n >= 0)
               dY = m_Values[2][nDexel][n].dMax ;
         }
      }

      else {
         int n = 0 ;
         double dY = m_Values[2][nDexel][n].dMin ;
         while ( n < nSize  &&  dY < dMaxY + dEps) {
            if ( dY > dMinY - dEps) {
               vfField.ptPApp.y = Clamp( dY, dMinY + dEpsClamp, dMaxY - dEpsClamp) ;
               vfField.vtVec = m_Values[2][nDexel][n].vtMinN ;
               vfField.nPropIndex = m_Values[2][nDexel][n].nToolMin ;
               bFound = true ;
               break ;
            }
            n += 1 ;
            if ( n < nSize)
               dY = m_Values[2][nDexel][n].dMin ;
         }
      }

      if ( ! bFound)
         vfField.ptPApp.y = 0.5 * ( dMinY + dMaxY) ;
   }

   else if ( nVec1[2] != nVec2[2]) {

      int nMinK = min( nVec1[2], nVec2[2]) * m_nDexVoxRatio ;
      int nMaxK = max( nVec1[2], nVec2[2]) * m_nDexVoxRatio ;

      double dMinZ = ( nMinK + 0.5) * m_dStep ;
      double dMaxZ = ( nMaxK + 0.5) * m_dStep ;

      vfField.ptPApp.x = ( nVec1[0] * m_nDexVoxRatio + 0.5) * m_dStep ;
      vfField.ptPApp.y = ( nVec1[1] * m_nDexVoxRatio + 0.5) * m_dStep ;

      int nDexel = ( nVec1[1] * m_nNx[0] + nVec1[0]) * m_nDexVoxRatio ;
      int nSize = int( m_Values[0][nDexel].size()) ;

      if ( bFirstCorner) {
         int n = nSize - 1 ;
         double dZ = m_Values[0][nDexel][n].dMax ;
         while ( n >= 0  &&  dZ > dMinZ - dEps) {
            if ( dZ < dMaxZ + dEps) {
               vfField.ptPApp.z = Clamp( dZ, dMinZ + dEpsClamp, dMaxZ - dEpsClamp) ;
               vfField.vtVec = m_Values[0][nDexel][n].vtMaxN ;
               vfField.nPropIndex = m_Values[0][nDexel][n].nToolMax ;
               bFound = true ;
               break ;
            }
            n -= 1 ;
            if ( n >= 0) 
               dZ = m_Values[0][nDexel][n].dMax ;
         }
      }

      else {
         int n = 0 ;
         double dZ = m_Values[0][nDexel][n].dMin ;
         while ( n < nSize  &&  dZ < dMaxZ + dEps) {
            if ( dZ > dMinZ - dEps) {
               vfField.ptPApp.z = Clamp( dZ, dMinZ + dEpsClamp, dMaxZ - dEpsClamp) ;
               vfField.vtVec = m_Values[0][nDexel][n].vtMinN ;
               vfField.nPropIndex = m_Values[0][nDexel][n].nToolMin ;
               bFound = true ;
               break ;
            }
            n += 1 ;
            if ( n < nSize)
               dZ = m_Values[0][nDexel][n].dMin ;
         }
      }

      if ( ! bFound)
         vfField.ptPApp.z = 0.5 * ( dMinZ + dMaxZ) ;
   }

   return bFound ;
}

//----------------------------------------------------------------------------
bool
VolZmap::IsValidVoxel( int nN) const
{
  // Calcolo il numero di voxel lungo X,Y e Z
   int nVoxNumX = m_nNx[0] / m_nDexVoxRatio + ( m_nNx[0] % m_nDexVoxRatio == 0 ? 1 : 2) ;  
   int nVoxNumY = m_nNy[0] / m_nDexVoxRatio + ( m_nNy[0] % m_nDexVoxRatio == 0 ? 1 : 2) ; 
   int nVoxNumZ = m_nNy[1] / m_nDexVoxRatio + ( m_nNy[1] % m_nDexVoxRatio == 0 ? 1 : 2) ;  
  // Verifico la validità del voxel
   return ( nN >= 0  &&  nN < ( nVoxNumX * nVoxNumY * nVoxNumZ)) ;
}

//----------------------------------------------------------------------------
bool 
VolZmap::IsValidVoxel( int nI, int nJ, int nK) const 
{
  // Calcolo il numero di voxel lungo X,Y e Z
   int nVoxNumX = int( m_nNx[0] / m_nDexVoxRatio + ( m_nNx[0] % m_nDexVoxRatio == 0 ? 1 : 2)) ;  
   int nVoxNumY = int( m_nNy[0] / m_nDexVoxRatio + ( m_nNy[0] % m_nDexVoxRatio == 0 ? 1 : 2)) ; 
   int nVoxNumZ = int( m_nNy[1] / m_nDexVoxRatio + ( m_nNy[1] % m_nDexVoxRatio == 0 ? 1 : 2)) ;  
  // Verifico la validità del voxel
   return ( nI >= -1  &&  nI <= nVoxNumX - 2  &&
            nJ >= -1  &&  nJ <= nVoxNumY - 2  &&
            nK >= -1  &&  nK <= nVoxNumZ - 2) ;
}

//----------------------------------------------------------------------------
bool
VolZmap::GetVoxIJKFromN( int nN, int& nI, int& nJ, int& nK) const
{
  // Calcolo il numero di voxel lungo X,Y e Z
   int nVoxNumX = int( m_nNx[0] / m_nDexVoxRatio + ( m_nNx[0] % m_nDexVoxRatio == 0 ? 1 : 2)) ;
   int nVoxNumY = int( m_nNy[0] / m_nDexVoxRatio + ( m_nNy[0] % m_nDexVoxRatio == 0 ? 1 : 2)) ;
   int nVoxNumZ = int( m_nNy[1] / m_nDexVoxRatio + ( m_nNy[1] % m_nDexVoxRatio == 0 ? 1 : 2)) ;

  // Controllo sulla validità del voxel
   if ( nN < 0 || nN >= nVoxNumX * nVoxNumY * nVoxNumZ)
      return false ;

  // Calcolo gli indici
   nI = ( nN % ( nVoxNumX * nVoxNumY)) % nVoxNumX - 1 ; 
   nJ = ( nN % ( nVoxNumX * nVoxNumY)) / nVoxNumX - 1 ;
   nK = ( nN / ( nVoxNumX * nVoxNumY)) - 1 ;

  // Controllo la sensatezza del risultato ottenuto
   return ( nI >= -1  &&  nI <= nVoxNumX - 2  &&
            nJ >= -1  &&  nJ <= nVoxNumY - 2  &&
            nK >= -1  &&  nK <= nVoxNumZ - 2) ;
}

//----------------------------------------------------------------------------
bool
VolZmap::GetVoxNFromIJK( int nI, int nJ, int nK, int& nN) const
{
  // Calcolo il numero di voxel lungo X,Y e Z
   int nVoxNumX = int( m_nNx[0] / m_nDexVoxRatio + ( m_nNx[0] % m_nDexVoxRatio == 0 ? 1 : 2)) ;
   int nVoxNumY = int( m_nNy[0] / m_nDexVoxRatio + ( m_nNy[0] % m_nDexVoxRatio == 0 ? 1 : 2)) ;
   int nVoxNumZ = int( m_nNy[1] / m_nDexVoxRatio + ( m_nNy[1] % m_nDexVoxRatio == 0 ? 1 : 2)) ;

  // Controllo la validità del voxel
   if ( nI < -1  ||  nI > nVoxNumX - 2  ||
        nJ < -1  ||  nJ > nVoxNumY - 2  ||
        nK < -1  ||  nK > nVoxNumZ - 2)
      return false ;

  // Calcolo il numero di Voxel
   nN = nVoxNumX * nVoxNumY * ( nK + 1) + nVoxNumX * ( nJ + 1) + nI + 1 ;

  // controllo sulla sensatezza del numero ottenuto
   return ( nN >= 0  &&  nN < nVoxNumX * nVoxNumY * nVoxNumZ) ;
}

//----------------------------------------------------------------------------
bool 
VolZmap::GetVoxelBox( int nI, int nJ, int nK, BBox3d& b3VoxBox) const
{
  // Calcolo il numero di voxel lungo X,Y e Z
   int nVoxNumX = int( m_nNx[0] / m_nDexVoxRatio + ( m_nNx[0] % m_nDexVoxRatio == 0 ? 1 : 2)) ;
   int nVoxNumY = int( m_nNy[0] / m_nDexVoxRatio + ( m_nNy[0] % m_nDexVoxRatio == 0 ? 1 : 2)) ;
   int nVoxNumZ = int( m_nNy[1] / m_nDexVoxRatio + ( m_nNy[1] % m_nDexVoxRatio == 0 ? 1 : 2)) ;

  // Controllo la validità del voxel
   if ( nI < -1  ||  nI > nVoxNumX - 2  ||
        nJ < -1  ||  nJ > nVoxNumY - 2  ||
        nK < -1  ||  nK > nVoxNumZ - 2)
      return false ;

  // Punti della diagonale del voxel
   Point3d ptCubeInf( ( nI * m_nDexVoxRatio + 0.5) * m_dStep,
                      ( nJ * m_nDexVoxRatio + 0.5) * m_dStep,
                      ( nK * m_nDexVoxRatio + 0.5) * m_dStep) ;
   Point3d ptCubeSup( ( ( nI + 1) * m_nDexVoxRatio + 0.5) * m_dStep,
                      ( ( nJ + 1) * m_nDexVoxRatio + 0.5) * m_dStep,
                      ( ( nK + 1) * m_nDexVoxRatio + 0.5) * m_dStep) ;
   b3VoxBox.Set( ptCubeInf, ptCubeSup) ;
   return true ;
}

//----------------------------------------------------------------------------
bool
VolZmap::GetBlockIJKFromN( int nBlock, int nIJK[]) const 
{
  // Controllo sulla validità del blocco
   if ( nBlock < 0  ||  nBlock >= int( m_nNumBlock))
      return false ;

  // Calcolo posizione del blocco nel reticolo 
   nIJK[0] = ( nBlock % int( m_nFracLin[0] * m_nFracLin[1])) % int( m_nFracLin[0]) ;
   nIJK[1] = ( nBlock % int( m_nFracLin[0] * m_nFracLin[1])) / int( m_nFracLin[0]) ;
   nIJK[2] = ( nBlock / int( m_nFracLin[0] * m_nFracLin[1])) ;

   return true ; 
}   

//----------------------------------------------------------------------------
bool
VolZmap::GetBlockNFromIJK( int nIJK[], int& nBlock) const
{
  // Controllo sulla validità degli indici i, j, k del blocco
   if ( nIJK[0] < 0 || nIJK[0] >= int( m_nFracLin[0]) ||
        nIJK[1] < 0 || nIJK[1] >= int( m_nFracLin[1]) || 
        nIJK[2] < 0 || nIJK[2] >= int( m_nFracLin[2]))
      return false ;

  // Determino il blocco
   nBlock = m_nFracLin[0] * m_nFracLin[1] * nIJK[2] + m_nFracLin[0] * nIJK[1] + nIJK[0] ;

   return true ;
}

//----------------------------------------------------------------------------
bool
VolZmap::GetBlockLimitsIJK( const int nIJK[], int nLimits[]) const
{
  // Controllo sulla validità degli indici i, j, k del blocco
   if ( nIJK[0] < 0  ||  nIJK[0] >= m_nFracLin[0] ||
        nIJK[1] < 0  ||  nIJK[1] >= m_nFracLin[1] || 
        nIJK[2] < 0  ||  nIJK[2] >= m_nFracLin[2])
      return false ;

   // Calcolo il numero di voxel lungo X,Y e Z
   int nVoxNumX = int( m_nNx[0] / m_nDexVoxRatio + ( m_nNx[0] % m_nDexVoxRatio == 0 ? 1 : 2)) ;
   int nVoxNumY = int( m_nNy[0] / m_nDexVoxRatio + ( m_nNy[0] % m_nDexVoxRatio == 0 ? 1 : 2)) ;
   int nVoxNumZ = int( m_nNy[1] / m_nDexVoxRatio + ( m_nNy[1] % m_nDexVoxRatio == 0 ? 1 : 2)) ;

   // Calcolo limiti per l'indice i
   nLimits[0] = ( nIJK[0] == 0 ? - 1 : nIJK[0] * m_nVoxNumPerBlock) ;
   nLimits[1] = ( nIJK[0] + 1 == m_nFracLin[0] ? 
                  nVoxNumX - 1 : ( nIJK[0] + 1) * m_nVoxNumPerBlock) ;

  // Calcolo limiti per l'indice j
   nLimits[2] = ( nIJK[1] == 0 ? - 1 : nIJK[1] * m_nVoxNumPerBlock) ;
   nLimits[3] = ( nIJK[1] + 1 == m_nFracLin[1] ? 
                  nVoxNumY - 1 : ( nIJK[1] + 1) * m_nVoxNumPerBlock) ;

  // Calcolo limiti per l'indice k
   nLimits[4] = ( nIJK[2] == 0 ? - 1 : nIJK[2] * m_nVoxNumPerBlock) ;
   nLimits[5] = ( nIJK[2] + 1 == m_nFracLin[2] ? 
                  nVoxNumZ - 1 : ( nIJK[2] + 1) * m_nVoxNumPerBlock) ;

   return true ;
}

//----------------------------------------------------------------------------
bool
VolZmap::GetBlockBox( const int nIJK[], BBox3d& b3VoxBox) const
{
  // Calcolo limiti sugli indici del blocco; se il blocco non è valido, errore
   int nLimits[6] ;
   if ( ! GetBlockLimitsIJK( nIJK, nLimits))
      return false ;
  // Costruisco il bounding-bx del blocco
   Point3d ptMinDiag( ( nLimits[0] * m_nDexVoxRatio + 0.5) * m_dStep,
                      ( nLimits[2] * m_nDexVoxRatio + 0.5) * m_dStep,
                      ( nLimits[4] * m_nDexVoxRatio + 0.5) * m_dStep) ;
   Point3d ptMaxDiag( ( nLimits[1] * m_nDexVoxRatio + 0.5) * m_dStep,
                      ( nLimits[3] * m_nDexVoxRatio + 0.5) * m_dStep,
                      ( nLimits[5] * m_nDexVoxRatio + 0.5) * m_dStep) ;
   b3VoxBox.Set( ptMinDiag, ptMaxDiag) ;
   return true ;
}

//----------------------------------------------------------------------------
bool
VolZmap::IsPointInsideVoxelApprox( int nI, int nJ, int nK, const Point3d& ptP, double dPrec) const 
{
  // Calcolo il numero di voxel lungo X,Y e Z
   int nVoxNumX = int( m_nNx[0] / m_nDexVoxRatio + ( m_nNx[0] % m_nDexVoxRatio == 0 ? 1 : 2)) ;
   int nVoxNumY = int( m_nNy[0] / m_nDexVoxRatio + ( m_nNy[0] % m_nDexVoxRatio == 0 ? 1 : 2)) ;
   int nVoxNumZ = int( m_nNy[1] / m_nDexVoxRatio + ( m_nNy[1] % m_nDexVoxRatio == 0 ? 1 : 2)) ;
  // Controllo sulla validità del voxel
   if ( nI < -1  ||  nI > nVoxNumX - 2  ||
        nJ < -1  ||  nJ > nVoxNumY - 2  ||
        nK < -1  ||  nK > nVoxNumZ - 2)
      return false ;
  // Se la tolleranza è minore di EPS_SMALL, la considero nulla
   if ( dPrec < EPS_ZERO)
      dPrec = 0 ;
  // Controllo che ogni coordinata stia dentro le relative limitazioni del voxel
   bool bI = ptP.x > ( nI * m_nDexVoxRatio + 0.5) * m_dStep - dPrec  &&
             ptP.x < ( ( nI + 1) * m_nDexVoxRatio + 0.5) * m_dStep + dPrec ;
   bool bJ = ptP.y > ( nJ * m_nDexVoxRatio + 0.5) * m_dStep - dPrec  &&
             ptP.y < ( ( nJ + 1) * m_nDexVoxRatio + 0.5) * m_dStep + dPrec ;
   bool bK = ptP.z > ( nK * m_nDexVoxRatio + 0.5) * m_dStep - dPrec  &&
             ptP.z < ( ( nK + 1) * m_nDexVoxRatio + 0.5) * m_dStep + dPrec ;

   return ( bI  &&  bJ  &&  bK) ;
}

//----------------------------------------------------------------------------
bool
VolZmap::GetPointVoxel( const Point3d& ptP, int& nVoxI, int& nVoxJ, int& nVoxK) const 
{
  // Calcolo il numero di voxel lungo X,Y e Z
   int nVoxNumX = int( m_nNx[0] / m_nDexVoxRatio + ( m_nNx[0] % m_nDexVoxRatio == 0 ? 1 : 2)) ;
   int nVoxNumY = int( m_nNy[0] / m_nDexVoxRatio + ( m_nNy[0] % m_nDexVoxRatio == 0 ? 1 : 2)) ;
   int nVoxNumZ = int( m_nNy[1] / m_nDexVoxRatio + ( m_nNy[1] % m_nDexVoxRatio == 0 ? 1 : 2)) ;
  // Calcolo gli indici del voxel
   nVoxI = int( floor( ( ptP.x - 0.5 * m_dStep) / ( m_dStep * m_nDexVoxRatio))) ;
   nVoxJ = int( floor( ( ptP.y - 0.5 * m_dStep) / ( m_dStep * m_nDexVoxRatio))) ;
   nVoxK = int( floor( ( ptP.z - 0.5 * m_dStep) / ( m_dStep * m_nDexVoxRatio))) ;
  // Controllo la validità del voxel
   return ( nVoxI >= - 1  &&  nVoxI <= nVoxNumX - 2  &&
            nVoxJ >= - 1  &&  nVoxJ <= nVoxNumY - 2  &&
            nVoxK >= - 1  &&  nVoxK <= nVoxNumZ - 2) ;
}

//----------------------------------------------------------------------------
bool
VolZmap::GetVoxelBlockIJK( const int nVoxIJK[], int nBlockIJK[]) const 
{
  // Calcolo il numero di voxel lungo X,Y e Z
   int nVoxNumX = int( m_nNx[0] / m_nDexVoxRatio + ( m_nNx[0] % m_nDexVoxRatio == 0 ? 1 : 2)) ;
   int nVoxNumY = int( m_nNy[0] / m_nDexVoxRatio + ( m_nNy[0] % m_nDexVoxRatio == 0 ? 1 : 2)) ;
   int nVoxNumZ = int( m_nNy[1] / m_nDexVoxRatio + ( m_nNy[1] % m_nDexVoxRatio == 0 ? 1 : 2)) ;
  // Controllo sull'ammissibilità del voxel
   if ( nVoxIJK[0] < -1  ||  nVoxIJK[0] > nVoxNumX - 2  ||
        nVoxIJK[1] < -1  ||  nVoxIJK[1] > nVoxNumY - 2  ||
        nVoxIJK[2] < -1  ||  nVoxIJK[2] > nVoxNumZ - 2)
      return false ;
  // Divisioni intere
   int nIntRatio0 = nVoxIJK[0] / m_nVoxNumPerBlock ;
   int nIntRatio1 = nVoxIJK[1] / m_nVoxNumPerBlock ;
   int nIntRatio2 = nVoxIJK[2] / m_nVoxNumPerBlock ;
  // Calcolo indici del blocco
   nBlockIJK[0] = ( nVoxIJK[0] == -1 ? 0 : max( 0, nIntRatio0 - ( nIntRatio0 == m_nFracLin[0] ? 1 : 0))) ;
   nBlockIJK[1] = ( nVoxIJK[1] == -1 ? 0 : max( 0, nIntRatio1 - ( nIntRatio1 == m_nFracLin[1] ? 1 : 0))) ;
   nBlockIJK[2] = ( nVoxIJK[2] == -1 ? 0 : max( 0, nIntRatio2 - ( nIntRatio2 == m_nFracLin[2] ? 1 : 0))) ;

   return true ;
}

//----------------------------------------------------------------------------
bool
VolZmap::GetVoxelBlockN( int nVox, int& nBlock) const
{
   int nVoxIJK[3] ;
   int nBlockIJK[3] ;
   return ( GetVoxIJKFromN( nVox, nVoxIJK[0], nVoxIJK[1], nVoxIJK[2])  &&
            GetVoxelBlockIJK( nVoxIJK, nBlockIJK)  &&
            GetBlockNFromIJK( nBlockIJK, nBlock)) ;
}

//----------------------------------------------------------------------------
bool
VolZmap::GetAdjBlockToBlock( int nBlockN, int nDeltaI, int nDeltaJ, int nDeltaK, int& nAdjBlockN) const
{
  // Test sulla validità degli incrementi su i,j,k
   if ( nDeltaI < -1  ||  nDeltaI > 1  ||
        nDeltaJ < -1  ||  nDeltaJ > 1  ||
        nDeltaK < -1  ||  nDeltaK > 1)
      return false ;

  // Determino blocco adiacente
   nAdjBlockN = nBlockN ;
   nAdjBlockN += nDeltaI ;
   nAdjBlockN += nDeltaJ * m_nFracLin[0] ;
   nAdjBlockN += nDeltaK * m_nFracLin[0] * m_nFracLin[1] ;

  // Se il blocco adiacente esiste restituisco vero, altrimenti falso.
   return ( nAdjBlockN > -1  &&  nAdjBlockN < int( m_nNumBlock)) ;
}

//----------------------------------------------------------------------------
bool 
VolZmap::IsAVoxelOnBoundary( const int nLimits[], const int nIJK[], bool bType) const
{
  // Calcolo il numero di voxel lungo X,Y e Z
   int nVoxNumX = int( m_nNx[0] / m_nDexVoxRatio + ( m_nNx[0] % m_nDexVoxRatio == 0 ? 1 : 2)) ;
   int nVoxNumY = int( m_nNy[0] / m_nDexVoxRatio + ( m_nNy[0] % m_nDexVoxRatio == 0 ? 1 : 2)) ;
   int nVoxNumZ = int( m_nNy[1] / m_nDexVoxRatio + ( m_nNy[1] % m_nDexVoxRatio == 0 ? 1 : 2)) ;

  // Test sulla validità dei limiti
   if ( nLimits[0] < -1  ||  nLimits[0] > nVoxNumX - 1  ||
        nLimits[1] < -1  ||  nLimits[1] > nVoxNumX - 1  ||
        nLimits[2] < -1  ||  nLimits[2] > nVoxNumY - 1  ||
        nLimits[3] < -1  ||  nLimits[3] > nVoxNumY - 1  ||
        nLimits[4] < -1  ||  nLimits[4] > nVoxNumZ - 1  ||
        nLimits[5] < -1  ||  nLimits[5] > nVoxNumZ - 1  )
      return false ;

  // Controllo sull'ammissibilità del voxel
   if ( nIJK[0] < -1  ||  nIJK[0] > nVoxNumX - 2  ||
        nIJK[1] < -1  ||  nIJK[1] > nVoxNumY - 2  ||
        nIJK[2] < -1  ||  nIJK[2] > nVoxNumZ - 2)
      return false ;

  // Se cerchiamo i voxel che sono sulla frontiera del blocco
   if ( ! bType) {
      return ( nIJK[0] == nLimits[0]  ||  nIJK[0] == nLimits[1] - 1  ||
               nIJK[1] == nLimits[2]  ||  nIJK[1] == nLimits[3] - 1  ||
               nIJK[2] == nLimits[4]  ||  nIJK[2] == nLimits[5] - 1) ;
   }

  // Altrimenti cerchiamo i voxel che confinano con quelli di altri blocchi
  // Condizione necessaria è che il voxel sia sulla frontiera
   if ( nIJK[0] == nLimits[0]  ||  nIJK[0] == nLimits[1] - 1  ||
        nIJK[1] == nLimits[2]  ||  nIJK[1] == nLimits[3] - 1  ||
        nIJK[2] == nLimits[4]  ||  nIJK[2] == nLimits[5] - 1) {
     // Indici del blocco
      int nCurrentBlockIJK[3] ;
      GetVoxelBlockIJK( nIJK, nCurrentBlockIJK) ;
     // Ciclo sulle possibili posizioni dei voxel adiacenti
      for ( int nDeltaI = -1 ; nDeltaI <= 1 ; ++ nDeltaI) {
         for ( int nDeltaJ = -1 ; nDeltaJ <= 1 ; ++ nDeltaJ) {
            for ( int nDeltaK = -1 ; nDeltaK <= 1 ; ++ nDeltaK) {
              // Evito controllo con se stesso
               if ( nDeltaI == 0  &&  nDeltaJ == 0  &&  nDeltaK == 0)
                  continue ;
              // Indici del voxel adiacente
               int nAdjIJK[3] = { nIJK[0] + nDeltaI, nIJK[1] + nDeltaJ, nIJK[2] + nDeltaK} ;
              // Determino (se esiste) il blocco in cui cade il voxel adiacente.
               int nAdjBlockIJK[3] ;
               if ( GetVoxelBlockIJK( nAdjIJK, nAdjBlockIJK)) {
                  if ( nAdjBlockIJK[0] > -1  &&  nAdjBlockIJK[0] < m_nFracLin[0]  &&
                       nAdjBlockIJK[1] > -1  &&  nAdjBlockIJK[1] < m_nFracLin[1]  &&
                       nAdjBlockIJK[2] > -1  &&  nAdjBlockIJK[2] < m_nFracLin[2]  && 
                       ( nCurrentBlockIJK[0] != nAdjBlockIJK[0]  ||
                         nCurrentBlockIJK[1] != nAdjBlockIJK[1]  ||
                         nCurrentBlockIJK[2] != nAdjBlockIJK[2])) {
                     return true ;
                  }
               }
            }
         }
      }
   }
   return false ;
}

//----------------------------------------------------------------------------
bool
VolZmap::IsAVoxelOnBoundary( const int nLimits[], const int nIJK[], int nDeltaIndex[]) const
{
  // Calcolo il numero di voxel lungo X,Y e Z
   int nVoxNumX = int( m_nNx[0] / m_nDexVoxRatio + ( m_nNx[0] % m_nDexVoxRatio == 0 ? 1 : 2)) ;
   int nVoxNumY = int( m_nNy[0] / m_nDexVoxRatio + ( m_nNy[0] % m_nDexVoxRatio == 0 ? 1 : 2)) ;
   int nVoxNumZ = int( m_nNy[1] / m_nDexVoxRatio + ( m_nNy[1] % m_nDexVoxRatio == 0 ? 1 : 2)) ;
  // Test sulla validità dei limiti
   if ( nLimits[0] < -1  ||  nLimits[0] > nVoxNumX - 1  ||
        nLimits[1] < -1  ||  nLimits[1] > nVoxNumX - 1  ||
        nLimits[2] < -1  ||  nLimits[2] > nVoxNumY - 1  ||
        nLimits[3] < -1  ||  nLimits[3] > nVoxNumY - 1  ||
        nLimits[4] < -1  ||  nLimits[4] > nVoxNumZ - 1  ||
        nLimits[5] < -1  ||  nLimits[5] > nVoxNumZ - 1  )
      return false ; 
  // Controllo sull'ammissibilità del voxel
   if ( nIJK[0] < -1  ||  nIJK[0] > nVoxNumX - 2  ||
        nIJK[1] < -1  ||  nIJK[1] > nVoxNumY - 2  ||
        nIJK[2] < -1  ||  nIJK[2] > nVoxNumZ - 2)
      return false ;
   nDeltaIndex[0] = 0 ;
   nDeltaIndex[1] = 0 ;
   nDeltaIndex[2] = 0 ;
   if ( nIJK[0] == nLimits[0])
      -- nDeltaIndex[0] ;
   else if ( nIJK[0] == nLimits[1] - 1)
      ++ nDeltaIndex[0] ;   
   if ( nIJK[1] == nLimits[2])
      -- nDeltaIndex[1];
   else if ( nIJK[1] == nLimits[3] - 1)
      ++ nDeltaIndex[1];   
   if ( nIJK[2] == nLimits[4])
      -- nDeltaIndex[2];
   else if ( nIJK[2] == nLimits[5] - 1)
      ++ nDeltaIndex[2];   
   return ( nDeltaIndex[0] != 0  ||  nDeltaIndex[1] != 0  ||  nDeltaIndex[2] != 0) ;
}

//----------------------------------------------------------------------------
// Per riferimento restituisce i minimi degli indici ijk riferiti ai Voxel.
// Si assume che il volume sia un parallelepipedo.
bool
VolZmap::GetFirstVoxIJK( int& i, int& j, int& k) const 
{
  // Se non è un parallelepipedo, errore
   if ( m_nShape != BOX)
      return false ;      
  // Calcolo  indici 
   int ni, nj ;
   for ( ni = 0 ; ni < 2 ; ++ ni) {
      bool bNotEmpty = false ;
      for ( nj = 0 ; nj < 2 ; ++ nj) {
         int nDex = nj * int( m_nNx[0]) + ni ;      
         if ( ! m_Values[0].empty()  &&  ! m_Values[0][nDex].empty()) {
            bNotEmpty = true ;
            break ;    
         }            
      }
      if ( bNotEmpty)
            break ;   
   }    
   int mi, mj ;
   for ( mi = 0 ; mi < 2 ; ++ mi) {
      bool bNotEmpty = false ;
      for ( mj = 0 ; mj < 2 ; ++ mj) {
         int nDex = mj * int( m_nNx[1]) + mi ;
         if ( ! m_Values[1].empty()  &&  ! m_Values[1][nDex].empty()) {
            bNotEmpty = true ;
            break ;
         }           
      }
      if ( bNotEmpty)
         break ;
   }

   i = ( ni % m_nDexVoxRatio != 0 ? (ni + m_nDexVoxRatio - 1) / m_nDexVoxRatio - 1 : ni / m_nDexVoxRatio - 1) ;
   j = ( nj % m_nDexVoxRatio != 0 ? (nj + m_nDexVoxRatio - 1) / m_nDexVoxRatio - 1 : nj / m_nDexVoxRatio - 1) ;
   k = ( mj % m_nDexVoxRatio != 0 ? (mj + m_nDexVoxRatio - 1) / m_nDexVoxRatio - 1 : mj / m_nDexVoxRatio - 1) ;
   return true ;
}

//----------------------------------------------------------------------------
// Per riferimento restituisce i massimi degli indici ijk riferiti ai Voxel.
// Si assume che il volume sia un parallelepipedo.
bool
VolZmap::GetLastVoxIJK( int& i, int& j, int& k) const
{
  // Se non è un parallelepipedo, errore
   if ( m_nShape != BOX)
      return false ;   
  // Calcolo indici     
   int ni = -1, nj = -1 ;
   for ( ni = int( m_nNx[0]) - 1 ; ni > int( m_nNx[0]) - 3 ; -- ni) {
      bool bNotEmpty = false ;
      for ( nj = int( m_nNy[0]) - 1 ; nj > int( m_nNy[0]) - 3 ; -- nj) {
         int nDex = nj * int( m_nNx[0]) + ni ;
         if ( ! m_Values[0].empty()  &&  ! m_Values[0][nDex].empty()) {
            bNotEmpty = true ;
            break ; 
         }        
      }
      if ( bNotEmpty)
         break ;  
   } 
   int mi = -1, mj = -1 ;
   for ( mi = int( m_nNx[1]) - 1 ; mi > int( m_nNx[1]) - 3 ; -- mi) {
      bool bNotEmpty = false ;
      for ( mj = int( m_nNy[1]) - 1 ; mj > int( m_nNy[1]) - 3 ; -- mj) {
         int nDex = mj * int( m_nNx[1]) + mi ;        
         if ( ! m_Values[1].empty()  &&  ! m_Values[1][nDex].empty()) {
            bNotEmpty = true ;
            break ;    
         }
      }
      if ( bNotEmpty)
         break ;  
   }
   i = ni / m_nDexVoxRatio ;
   j = nj / m_nDexVoxRatio ;
   k = mj / m_nDexVoxRatio ; 
   return true ;
}

//----------------------------------------------------------------------------
// Il volume deve essere un parallelepipedo.
// Verifica se il voxel contiene uno dei dodici spigoli del parallelepipedo.
bool 
VolZmap::IsVoxelOnBoxEdge( int i, int j, int k) const
{
  // Se non è un parallelepipedo, errore
   if ( m_nShape != BOX)
      return false ;     
  // Determino il primo nodo pieno della mappa
   int nFirstVoxI, nFirstVoxJ, nFirstVoxK ;  
   GetFirstVoxIJK( nFirstVoxI, nFirstVoxJ, nFirstVoxK) ;
  // Determino il primo nodo pieno della mappa   
   int nLastVoxI, nLastVoxJ, nLastVoxK ;  
   GetLastVoxIJK( nLastVoxI, nLastVoxJ, nLastVoxK) ;

  // Determino se il voxel è su un edge del box
   int nIndexOnLimitNum = 0 ; 
   if ( i == nFirstVoxI  ||  i == nLastVoxI) 
      ++ nIndexOnLimitNum ; 
   if ( j == nFirstVoxJ  ||  j == nLastVoxJ)
      ++ nIndexOnLimitNum ;  
   if ( k == nFirstVoxK  ||  k == nLastVoxK)
      ++ nIndexOnLimitNum ;  
   return ( nIndexOnLimitNum >= 2) ;  
}

//----------------------------------------------------------------------------
bool
VolZmap::IsTriangleOnBorder( const Triangle3dEx& trTria, const int nBlockLimits[], const int nVoxIJK[]) const
{
  // Punti del triangolo sulla griglia 
   Point3d ptFirstGrPt = trTria.GetP( 1) ;
   Point3d ptSecondGrPt = trTria.GetP( 2) ;

  // Verifico se tali punti sono sulla griglia 
   for ( int nC = 0 ; nC < 3 ; ++ nC) {
      if ( nVoxIJK[nC] == nBlockLimits[2*nC]) {
         double dGrid = ( nBlockLimits[2*nC] * m_nDexVoxRatio + 0.5) * m_dStep ;
         if ( abs( ptFirstGrPt.v[nC] - dGrid) < EPS_SMALL  &&
              abs( ptSecondGrPt.v[nC] - dGrid) < EPS_SMALL) 
            return true ;
      }
      if ( nVoxIJK[nC] + 1 == nBlockLimits[2*nC+1]) {
         double dGrid = ( nBlockLimits[2*nC+1] * m_nDexVoxRatio + 0.5) * m_dStep ;
         if ( abs( ptFirstGrPt.v[nC] - dGrid) < EPS_SMALL  &&
              abs( ptSecondGrPt.v[nC] - dGrid) < EPS_SMALL) 
            return true ;
      }
   }

   return false ;
}

//----------------------------------------------------------------------------
bool
VolZmap::CheckForFanNodeInterferance( AppliedVector BaseVecField[], Point3d& ptFanVert,
                                      int nBasePointNum, int nVoxConfig, int nVoxI, int nVoxJ, int nVoxK) const
{
   bool bInterferance = false ;
   int nFilterIndex = 1 ;
   for ( int n = 0 ; n < 7 ; ++ n) {
      bool bNodeInside = true ;
      if ( ( nVoxConfig & nFilterIndex) != 0) {
         int nNodeI = nVoxI ;
         int nNodeJ = nVoxJ ;
         int nNodeK = nVoxK ;
         if ( ( n + 1) % 4 >= 2)
            ++ nNodeI ;
         if ( n % 4 >= 2)
            ++ nNodeJ ;
         if ( n / 4 == 1)
            ++ nNodeK ;
         Point3d ptNodeP( ( nNodeI * m_nDexVoxRatio + 0.5) * m_dStep,
                          ( nNodeJ * m_nDexVoxRatio + 0.5) * m_dStep,
                          ( nNodeK * m_nDexVoxRatio + 0.5) * m_dStep) ;
         for ( int nV = 0 ; nV < nBasePointNum ; ++ nV) {
            int nW = ( nV + 1) % nBasePointNum ;
            Plane3d plPlane ;
            Vector3d vtN = ( BaseVecField[nW].ptPApp - ptFanVert) ^ ( BaseVecField[nV].ptPApp - ptFanVert) ;
            if ( plPlane.Set( ptFanVert, vtN)  &&  DistPointPlane( ptNodeP, plPlane) > -EPS_SMALL) {
               bNodeInside = false ;
               break ;
            }
         }
         if ( bNodeInside) {
            for ( int nV = 2 ; nV < nBasePointNum ; ++ nV) {
               Plane3d plPlane ;
               Vector3d vtN = ( BaseVecField[nV - 1].ptPApp - BaseVecField[0].ptPApp) ^ ( BaseVecField[nV].ptPApp - BaseVecField[0].ptPApp) ;
               if ( plPlane.Set( BaseVecField[0].ptPApp, vtN)  &&  DistPointPlane( ptNodeP, plPlane) > -EPS_SMALL) {
                  bNodeInside = false ;
                  break ;
               }
            }
         }
      }
      else
         bNodeInside = false ;

      bInterferance = bInterferance  ||  bNodeInside ;

      nFilterIndex *= 2 ;
   }

   return bInterferance ;
}

//----------------------------------------------------------------------------
bool
VolZmap::ProcessVoxContXY( FlatVoxelContainer& VoxContXY, int nBlock, bool bPlus) const
{
   for ( auto it = VoxContXY.cbegin() ;
         it != VoxContXY.end() ;
         it = VoxContXY.begin()) {

      int nN = it->first ;
      int nI, nJ, nK ;
      GetVoxIJKFromN( nN, nI, nJ, nK) ;

     // Costruzione del primo rettangolo: un singolo voxel
      int nMinI = nI ;
      int nMinJ = nJ ;
      int nMaxI = nI ;
      int nMaxJ = nJ ;
      int nToolNum = it->second.nTool ; 
      double dCordZ = it->second.dHeigth ;

     // Flag sul ritrovamento di un rettangolo più grande.
      bool bOkI = true ;
      bool bOkJ = true ;
      while ( bOkI  ||  bOkJ) {
        // Se precedente espansione ok e lato I più lungo o non ok J 
         if ( bOkI  &&  ( nMaxI - nMinI >= nMaxJ - nMinJ  ||  ! bOkJ)) {
           // Analizzo linea superiore
            bool bSupJ = true ;
            for ( int i = nMinI ; bSupJ && i <= nMaxI ; ++ i) {
               if ( ! Find( VoxContXY, i, nMaxJ + 1, nK, dCordZ, nToolNum))
                  bSupJ = false ;
            }
           // Se linea superiore accettata, elimino i voxel
            if ( bSupJ) {
               for ( int i = nMinI ; i <= nMaxI ; ++ i)
                  Remove( VoxContXY, i, nMaxJ + 1, nK) ;
               ++ nMaxJ ;
            }
           // Analizzo linea inferiore
            bool bInfJ = true ;
            for ( int i = nMinI ; bInfJ && i <= nMaxI ; ++ i) {
               if ( ! Find( VoxContXY, i, nMinJ - 1, nK, dCordZ, nToolNum))
                  bInfJ = false ;
            }
           // Se linea inferiore accettata, elimino i voxel
            if ( bInfJ) {
               for ( int i = nMinI ; i <= nMaxI ; ++ i)
                  Remove( VoxContXY, i, nMinJ - 1, nK) ;
               -- nMinJ ;
            }
            bOkI = bSupJ || bInfJ ;
         }
        // altrimenti se precedente espansione J ok
         else if ( bOkJ) {
           // Analizzo linea destra
            bool bSupI = true ;
            for ( int j = nMinJ ; bSupI && j <= nMaxJ ; ++ j) { 
               if ( ! Find( VoxContXY, nMaxI + 1, j, nK, dCordZ, nToolNum))
                  bSupI = false ;
            }
           // Se linea destra accettata, elimino i voxel
            if ( bSupI) {
               for ( int j = nMinJ ; j <= nMaxJ ; ++ j)
                  Remove( VoxContXY, nMaxI + 1, j, nK) ;
               ++ nMaxI ;
            }
           // Analizzo linea sinistra
            bool bInfI = true ;
            for ( int j = nMinJ ; bInfI && j <= nMaxJ ; ++ j) {
               if ( ! Find( VoxContXY, nMinI - 1, j, nK, dCordZ, nToolNum))
                  bInfI = false ;
            }
           // Se linea sinistra accettata, elimino i voxel
            if ( bInfI) {
               for ( int j = nMinJ ; j <= nMaxJ ; ++ j)
                  Remove( VoxContXY, nMinI - 1, j, nK) ;
               -- nMinI ;
            }
            bOkJ = bSupI || bInfI ;
         }
      }

      Point3d ptT0( ( nMinI * m_nDexVoxRatio + 0.5) * m_dStep,
                    ( nMinJ * m_nDexVoxRatio + 0.5) * m_dStep, dCordZ) ;
      Point3d ptT1( ( ( nMaxI + 1) * m_nDexVoxRatio + 0.5) * m_dStep, 
                    ( ( nMinJ * m_nDexVoxRatio + 0.5) * m_dStep), dCordZ) ;
      Point3d ptT2( ( ( nMaxI + 1) * m_nDexVoxRatio + 0.5) * m_dStep, 
                    ( ( ( nMaxJ + 1) * m_nDexVoxRatio + 0.5) * m_dStep), dCordZ) ;
      Point3d ptT3( ( nMinI * m_nDexVoxRatio + 0.5) * m_dStep,
                    ( ( ( nMaxJ + 1) * m_nDexVoxRatio + 0.5) * m_dStep), dCordZ) ;

      ptT0.ToGlob( m_MapFrame) ;
      ptT1.ToGlob( m_MapFrame) ;
      ptT2.ToGlob( m_MapFrame) ;
      ptT3.ToGlob( m_MapFrame) ;

      Triangle3dEx Tria0, Tria1 ;
      if ( bPlus) {
         Tria0.Set( ptT0, ptT1, ptT3) ;
         Tria1.Set( ptT1, ptT2, ptT3) ;
      }
      else {
         Tria0.Set( ptT0, ptT3, ptT1) ;
         Tria1.Set( ptT1, ptT3, ptT2) ;
      }
      Tria0.SetGrade( nToolNum) ;
      Tria1.SetGrade( nToolNum) ;
      Tria0.Validate( true) ;
      Tria1.Validate( true) ;
     // Aggiungo alla lista
      m_BlockBigTria[nBlock].emplace_back( Tria0) ;
      m_BlockBigTria[nBlock].emplace_back( Tria1) ;

     // Elimino il voxel da cui sono partito a ingrandire.
      VoxContXY.erase( nN) ;
   }

   return true ;
}

//----------------------------------------------------------------------------
bool
VolZmap::ProcessVoxContYZ( FlatVoxelContainer& VoxContYZ, int nBlock, bool bPlus) const
{
   for ( auto it = VoxContYZ.begin() ;
         it != VoxContYZ.end() ;
         it = VoxContYZ.begin()) {

      int nN = it->first ;
      int nI, nJ, nK ;
      GetVoxIJKFromN( nN, nI, nJ, nK) ;

     // Costruzione del primo rettangolo: un singolo voxel
      int nMinJ = nJ ;
      int nMinK = nK ;
      int nMaxJ = nJ ;
      int nMaxK = nK ;
      int nToolNum = it->second.nTool ; 
      double dCordX = it->second.dHeigth ;

     // Flag sul ritrovamento di un rettangolo più grande.
      bool bOkJ = true ;
      bool bOkK = true ;
      while ( bOkJ  ||  bOkK) {
        // Se precedente espansione ok e lato J più lungo o non ok K
         if ( bOkJ  &&  ( nMaxJ - nMinJ >= nMaxK - nMinK  ||  ! bOkK)) {
           // Analizzo linea superiore
            bool bSupK = true ;
            for ( int j = nMinJ ; bSupK && j <= nMaxJ ; ++ j) {
               if ( ! Find( VoxContYZ, nI, j, nMaxK + 1, dCordX, nToolNum))
                  bSupK = false ;
            }
           // Se linea superiore accettata, elimino i voxel
            if ( bSupK) {
               for ( int j = nMinJ ; j <= nMaxJ ; ++ j)
                  Remove( VoxContYZ, nI, j, nMaxK + 1) ;
               ++ nMaxK ;
            }
           // Analizzo linea inferiore
            bool bInfK = true ;
            for ( int j = nMinJ ; bInfK && j <= nMaxJ ; ++ j) {
               if ( ! Find( VoxContYZ, nI, j, nMinK - 1, dCordX, nToolNum))
                  bInfK = false ;
            }
           // Se linea inferiore accettata, elimino i voxel
            if ( bInfK) {
               for ( int j = nMinJ ; j <= nMaxJ ; ++ j)
                  Remove( VoxContYZ, nI, j, nMinK - 1) ;
               -- nMinK ;
            }
            bOkJ = bSupK || bInfK ;
         }
        // altrimenti se precedente espansione K ok
         else if ( bOkK) {     
           // Analizzo linea destra
            bool bSupJ = true ;
            for ( int k = nMinK ; bSupJ && k <= nMaxK ; ++ k) { 
               if ( ! Find( VoxContYZ, nI, nMaxJ + 1, k, dCordX, nToolNum))
                  bSupJ = false ;
            }
           // Se linea destra accettata, elimino i voxel
            if ( bSupJ) {
               for ( int k = nMinK ; k <= nMaxK ; ++ k)
                  Remove( VoxContYZ, nI, nMaxJ + 1, k) ;
               ++ nMaxJ ;
            }
           // Analizzo linea sinistra
            bool bInfJ = true ;
            for ( int k = nMinK ; bInfJ && k <= nMaxK ; ++ k) {
               if ( ! Find( VoxContYZ, nI, nMinJ - 1, k, dCordX, nToolNum))
                  bInfJ = false ;
            }
           // Se linea sinistra accettata, elimino i voxel
            if ( bInfJ) {
               for ( int k = nMinK ; k <= nMaxK ; ++ k)
                  Remove( VoxContYZ, nI, nMinJ - 1, k) ;
               -- nMinJ ;
            }
            bOkK = bSupJ || bInfJ ;
         }
      }

      Point3d ptT0( dCordX, ( nMinJ * m_nDexVoxRatio + 0.5) * m_dStep, 
                            ( ( nMinK * m_nDexVoxRatio + 0.5) * m_dStep)) ;
      Point3d ptT1( dCordX, ( ( nMaxJ + 1) * m_nDexVoxRatio + 0.5) * m_dStep, 
                            ( ( nMinK * m_nDexVoxRatio + 0.5) * m_dStep)) ;
      Point3d ptT2( dCordX, ( ( nMaxJ + 1) * m_nDexVoxRatio + 0.5) * m_dStep, 
                            ( ( ( nMaxK + 1) * m_nDexVoxRatio + 0.5) * m_dStep)) ;
      Point3d ptT3( dCordX, ( nMinJ * m_nDexVoxRatio + 0.5) * m_dStep, 
                            ( ( ( nMaxK + 1) * m_nDexVoxRatio + 0.5) * m_dStep)) ;

      ptT0.ToGlob( m_MapFrame) ;
      ptT1.ToGlob( m_MapFrame) ;
      ptT2.ToGlob( m_MapFrame) ;
      ptT3.ToGlob( m_MapFrame) ;

      Triangle3dEx Tria0, Tria1 ;
      if ( bPlus) {
         Tria0.Set( ptT0, ptT1, ptT3) ;
         Tria1.Set( ptT1, ptT2, ptT3) ;
      }
      else {
         Tria0.Set( ptT0, ptT3, ptT1) ;
         Tria1.Set( ptT1, ptT3, ptT2) ;
      }
      Tria0.SetGrade( nToolNum) ;
      Tria1.SetGrade( nToolNum) ;
      Tria0.Validate( true) ;
      Tria1.Validate( true) ;
     // Aggiungo alla lista
      m_BlockBigTria[nBlock].emplace_back( Tria0) ;
      m_BlockBigTria[nBlock].emplace_back( Tria1) ;

     // Elimino il voxel da cui sono partito a ingrandire.
      VoxContYZ.erase( nN) ;
   }

   return true ;
}

//----------------------------------------------------------------------------
bool
VolZmap::ProcessVoxContXZ( FlatVoxelContainer& VoxContXZ, int nBlock, bool bPlus) const
{
   for ( auto it = VoxContXZ.begin() ;
         it != VoxContXZ.end() ;
         it = VoxContXZ.begin()) {

      int nN = it->first ;
      int nI, nJ, nK ;
      GetVoxIJKFromN( nN, nI, nJ, nK) ;

     // Costruzione del primo rettangolo: un singolo voxel
      int nMinI = nI ;
      int nMinK = nK ;
      int nMaxI = nI ;
      int nMaxK = nK ;
      int nToolNum = it->second.nTool ; 
      double dCordY = it->second.dHeigth ;

     // Flag sul ritrovamento di un rettangolo più grande.
      bool bOkI = true ;
      bool bOkK = true ;
      while ( bOkI  ||  bOkK) {
        // Se lato I più lungo e precedente espansione ok
         if ( bOkI  &&  ( nMaxI - nMinI >= nMaxK - nMinK  ||  ! bOkK)) {
           // Analizzo linea superiore
            bool bSupK = true ;
            for ( int i = nMinI ; bSupK && i <= nMaxI ; ++ i) {
               if ( ! Find( VoxContXZ, i, nJ, nMaxK + 1, dCordY, nToolNum))
                  bSupK = false ;
            }
           // Se linea superiore accettata, elimino i voxel
            if ( bSupK) {
               for ( int i = nMinI ; i <= nMaxI ; ++ i)
                  Remove( VoxContXZ, i, nJ, nMaxK + 1) ;
               ++ nMaxK ;
            }
           // Analizzo linea inferiore
            bool bInfK = true ;
            for ( int i = nMinI ; bInfK && i <= nMaxI ; ++ i) {
               if ( ! Find( VoxContXZ, i, nJ, nMinK - 1, dCordY, nToolNum))
                  bInfK = false ;
            }
           // Se linea inferiore accettata, elimino i voxel
            if ( bInfK) {
               for ( int i = nMinI ; i <= nMaxI ; ++ i)
                  Remove( VoxContXZ, i, nJ, nMinK - 1) ;
               -- nMinK ;
            }
            bOkI = bSupK || bInfK ;
         }
        // altrimenti se precedente espansione K ok
         else if ( bOkK) {     
           // Analizzo linea destra
            bool bSupI = true ;
            for ( int k = nMinK ; bSupI && k <= nMaxK ; ++ k) {
               if ( ! Find( VoxContXZ, nMaxI + 1, nJ, k, dCordY, nToolNum))
                  bSupI = false ;
            }
           // Se linea destra accettata, elimino i voxel
            if ( bSupI) {
               for ( int k = nMinK ; k <= nMaxK ; ++ k)
                  Remove( VoxContXZ, nMaxI + 1, nJ, k) ;
               ++ nMaxI ;
            }
           // Analizzo linea sinistra
            bool bInfI = true ;
            for ( int k = nMinK ; bInfI && k <= nMaxK ; ++ k) {
               if ( ! Find( VoxContXZ, nMinI - 1, nJ, k, dCordY, nToolNum))
                  bInfI = false ;
            }
           // Se linea sinistra accettata, elimino i voxel
            if ( bInfI) {
               for ( int k = nMinK ; k <= nMaxK ; ++ k)
                  Remove( VoxContXZ, nMinI - 1, nJ, k) ;
               -- nMinI ;
            }
            bOkK = bSupI || bInfI ;
         }
      }

      Point3d ptT0( ( nMinI * m_nDexVoxRatio + 0.5) * m_dStep, dCordY,
                    ( ( nMinK * m_nDexVoxRatio + 0.5) * m_dStep)) ;
      Point3d ptT1( ( ( nMaxI + 1) * m_nDexVoxRatio + 0.5) * m_dStep, dCordY,
                    ( ( nMinK * m_nDexVoxRatio + 0.5) * m_dStep)) ;
      Point3d ptT2( ( ( nMaxI + 1) * m_nDexVoxRatio + 0.5) * m_dStep, dCordY, 
                    ( ( ( nMaxK + 1) * m_nDexVoxRatio + 0.5) * m_dStep)) ;
      Point3d ptT3( ( nMinI * m_nDexVoxRatio + 0.5) * m_dStep, dCordY,
                    ( ( ( nMaxK + 1) * m_nDexVoxRatio + 0.5) * m_dStep)) ;

      ptT0.ToGlob( m_MapFrame) ;
      ptT1.ToGlob( m_MapFrame) ;
      ptT2.ToGlob( m_MapFrame) ;
      ptT3.ToGlob( m_MapFrame) ;

      Triangle3dEx Tria0, Tria1 ;
      if ( bPlus) {
         Tria0.Set( ptT0, ptT3, ptT1) ;
         Tria1.Set( ptT1, ptT3, ptT2) ;
      }
      else {
         Tria0.Set( ptT0, ptT1, ptT3) ;
         Tria1.Set( ptT1, ptT2, ptT3) ;
      }
      Tria0.SetGrade( nToolNum) ;
      Tria1.SetGrade( nToolNum) ;
      Tria0.Validate( true) ;
      Tria1.Validate( true) ;
     // Aggiungo alla lista
      m_BlockBigTria[nBlock].emplace_back( Tria0) ;
      m_BlockBigTria[nBlock].emplace_back( Tria1) ;

     // Elimino il voxel da cui sono partito a ingrandire.
      VoxContXZ.erase( nN) ;
   }

   return true ;
}

//----------------------------------------------------------------------------
bool
VolZmap::Find( const FlatVoxelContainer& VoxCont, int nI, int nJ, int nK, double dPos, int nTool) const
{
  // indice globale del voxel
   int nN ;
   if ( ! GetVoxNFromIJK( nI, nJ, nK, nN))
      return false ;
  // cerco il voxel nel contenitore
   auto iter = VoxCont.find( nN) ;
   return ( iter != VoxCont.end()  &&  
            abs( dPos - iter->second.dHeigth) < EPS_SMALL  &&
            nTool == iter->second.nTool) ;
}

//----------------------------------------------------------------------------
bool
VolZmap::Remove( FlatVoxelContainer& VoxCont, int nI, int nJ, int nK) const
{
   int nN ;
   if ( ! GetVoxNFromIJK( nI, nJ, nK, nN))
      return false ;
   return ( VoxCont.erase( nN) > 0) ;
}

//----------------------------------------------------------------------------
bool
VolZmap::GetEdges( ICURVEPOVECTOR& vpCurve) const
{
  // Se mappa singola, non calcolabili
   if ( m_nMapNum == 1)
      return false ;
   
  // Garantisco grafica aggiornata
   UpdateTripleMapGraphics() ;

  // Vettore dei segmenti di feature (coppie di punti)
   BIPNTVECTOR vBpt ; 

  // Ciclo sui triangoli feature all'interno dei blocchi
   for ( int nBlock = 0 ; nBlock < m_nNumBlock ; ++ nBlock) {
     // Numero di voxel in cui si presentano sharp feature
      int nVoxelNum = int( m_BlockSharpTria[nBlock].size()) ;
     // Ciclo sui voxel con sharp feature
      for ( int n = 0 ; n < nVoxelNum ; ++ n) {
         const SharpTriaStruct& SharpTria = m_BlockSharpTria[nBlock][n] ;        
        // Ciclo sulle componenti connesse del voxel
         int nNumCompo = int( SharpTria.ptCompoVert.size()) ;
         for ( int nCompo = 0 ; nCompo < nNumCompo ; ++ nCompo) {
            const TRIA3DEXVECTOR& vTria = SharpTria.vCompoTria[nCompo] ;
           // Numero di triangoli della componente connessa
            int nTriNum = int( vTria.size()) ;                                                 
            for ( int nTri = 0 ; nTri < nTriNum ; ++ nTri) {
              // Il triangolo ha un lato sulla sharp-feature
               if ( vTria[nTri].GetEdgeFlag( 0)) {
                  const Point3d& ptPS = vTria[nTri].GetP( 0) ;
                  const Point3d& ptPE = vTria[nTri].GetP( 1) ;
                 // Segmento sharp-feature
                  vBpt.emplace_back( ptPS, ptPE) ;
               }
               else if ( vTria[nTri].GetEdgeFlag( 1)) {
                  const Point3d& ptPS = vTria[nTri].GetP( 1) ;
                  const Point3d& ptPE = vTria[nTri].GetP( 2) ;
                 // Segmento sharp-feature
                  vBpt.emplace_back( ptPS, ptPE) ;
               }
               else if ( vTria[nTri].GetEdgeFlag( 2)) {
                  const Point3d& ptPS = vTria[nTri].GetP( 2) ;
                  const Point3d& ptPE = vTria[nTri].GetP( 0) ;
                 // Segmento sharp-feature
                  vBpt.emplace_back( ptPS, ptPE) ;
               }
            }         
         }     
      }
   }

  // Ciclo sui triangoli feature al confine tra blocchi
   for ( int tB = 0 ; tB < m_nNumBlock ; ++ tB) {     
     // Numero di voxel nel blocco corrente
      int nVoxelNum = int( m_InterBlockSharpTria[tB].size()) ;   
     // Ciclo sui voxel del blocco
      for ( int tV = 0 ; tV < nVoxelNum ; ++ tV) {             
        // Numero di componenti connesse del voxel
         int nCompoNum = int( m_InterBlockSharpTria[tB][tV].ptCompoVert.size()) ;      
        // Ciclo sulle componenti connesse
         for ( int tC = 0 ; tC < nCompoNum ; ++ tC) {          
           // Numero di triangoli della componente connessa
            int nTriNum = int( m_InterBlockSharpTria[tB][tV].vCompoTria[tC].size()) ;              
           // Ciclo sui triangoli
            for ( int tT = 0 ; tT < nTriNum ; ++ tT) {                                  
              // Il triangolo ha un lato sulla sharp-feature
               if ( m_InterBlockSharpTria[tB][tV].vCompoTria[tC][tT].GetEdgeFlag( 0)) {
                  const Point3d& ptPS = m_InterBlockSharpTria[tB][tV].vCompoTria[tC][tT].GetP( 0) ;
                  const Point3d& ptPE = m_InterBlockSharpTria[tB][tV].vCompoTria[tC][tT].GetP( 1) ;
                 // Segmento sharp-feature
                  vBpt.emplace_back( ptPS, ptPE) ;
               }
               else if ( m_InterBlockSharpTria[tB][tV].vCompoTria[tC][tT].GetEdgeFlag( 1)) {
                  const Point3d& ptPS = m_InterBlockSharpTria[tB][tV].vCompoTria[tC][tT].GetP( 1) ;
                  const Point3d& ptPE = m_InterBlockSharpTria[tB][tV].vCompoTria[tC][tT].GetP( 2) ;
                 // Segmento sharp-feature
                  vBpt.emplace_back( ptPS, ptPE) ;
               }
               else if ( m_InterBlockSharpTria[tB][tV].vCompoTria[tC][tT].GetEdgeFlag( 2)) {
                  const Point3d& ptPS = m_InterBlockSharpTria[tB][tV].vCompoTria[tC][tT].GetP( 2) ;
                  const Point3d& ptPE = m_InterBlockSharpTria[tB][tV].vCompoTria[tC][tT].GetP( 0) ;
                 // Segmento sharp-feature
                  vBpt.emplace_back( ptPS, ptPE) ;
               }
            }     
         }
      }
   }

  // Concateno le curve rispettando le biforcazioni
   const double dToler = 20 * EPS_SMALL ;
   ChainCurves chainC ;
   chainC.Init( true, dToler, int( vBpt.size())) ;
   for ( int i = 0 ; i < int( vBpt.size()) ; ++ i) {
      Vector3d vtDir = vBpt[i].second - vBpt[i].first ;
      vtDir.Normalize() ;
      if ( ! chainC.AddCurve( i + 1, vBpt[i].first, vtDir, vBpt[i].second, vtDir))
         return false ;
   }
  // recupero i percorsi concatenati
   Point3d ptNear = ( vBpt.empty() ? ORIG : vBpt[0].first) ;
   INTVECTOR vId ;
   while ( chainC.GetChainFromNear( ptNear, true, vId)) {
     // creo una curva composita
      PtrOwner<CurveComposite> pCrvCompo( CreateBasicCurveComposite()) ;
      if ( IsNull( pCrvCompo))
         return false ;
     // recupero gli estremi dei segmenti e li inserisco nella composita
      for ( int i = 0 ; i < int( vId.size()) ; ++ i) {
        // recupero indice e flag di inversione
         int nInd = abs( vId[i]) - 1 ;
         bool bInvert = ( vId[i] < 0) ;
        // se primo segmento, inserisco il punto iniziale
         if ( i == 0)
            pCrvCompo->AddPoint( ! bInvert ? vBpt[nInd].first : vBpt[nInd].second) ;
        // aggiungo il punto finale
         pCrvCompo->AddLine( ! bInvert ? vBpt[nInd].second : vBpt[nInd].first) ;
      }
     // aggiorno il prossimo near
      pCrvCompo->GetEndPoint( ptNear) ;
     // se lunghezza curva inferiore a 5 volte la tolleranza, la salto
      double dCrvLen ;
      if ( ! pCrvCompo->GetLength( dCrvLen)  ||  dCrvLen < 5 * dToler)
         continue ;
     // unisco segmenti allineati
      pCrvCompo->MergeCurves( dToler, ANG_TOL_STD_DEG) ;
     // la salvo
      vpCurve.emplace_back( Release( pCrvCompo)) ;
   }

   return true ;
}