//----------------------------------------------------------------------------
//  EgalTech 2015-2016
//----------------------------------------------------------------------------
// File : VolZmap.cpp         Data : 22.01.15           Versione : 1.6a4
// Contenuto : Implementazione della classe Volume Zmap (tre griglie)
//
//
//
// Modifiche : 22.01.15 DS Creazione modulo.
//
//
//----------------------------------------------------------------------------

//--------------------------- Include ----------------------------------------

#include "stdafx.h"
#include "CurveLine.h"
#include "VolZmap.h"
#include "GeoConst.h"
#include "/EgtDev/Include/EGkStmFromCurves.h"
#include "/EgtDev/Include/EGkIntersLineSurfTm.h"
#include "/EgtDev/Include/EGkIntersLinePlane.h"
#include "/EgtDev/Include/EgtNumUtils.h"
#include <future>

using namespace std ;

// ------------------------- CREAZIONE MAPPA --------------------------------------------------------------------------------------

//----------------------------------------------------------------------------
bool
VolZmap::Create( const Point3d& ptO, double dDimX, double dDimY, double dDimZ, double dStep, bool bTriDex, int* nError)
{
  // Controlli sull'ammissibilità delle dimensioni lineari del grezzo e del passo
   if ( dStep < EPS_SMALL  ||  dDimX < EPS_SMALL  ||  dDimY < EPS_SMALL  ||  dDimZ < EPS_SMALL)
      return false ;

  // Il passo di discretizzazione non può essere inferiore a 100 * EPS_SMALL
   m_dStep = max( dStep, 100 * EPS_SMALL) ;

  // Aggiorno la dimensione della mappa 1 o 3
   m_nMapNum = ( bTriDex ? 3 : 1) ;

  // Disponendo i sistemi di riferimento in una successione, le coordinate x,y,z
  // di uno si ottengono da una permutazione ciclica di quelle del precedente sistema.
  // es: X(n) = Z(n-1), Y(n) = X(n-1), Z(n) = Y(n-1)

  // Definisco il sistema di riferimento intrinseco
   m_MapFrame.Set( ptO, X_AX, Y_AX, Z_AX) ;

  // Definisco i vettori dei limiti su indici
   m_nNx[0] = max( int( ( dDimX + EPS_SMALL) / m_dStep + 0.5), 1) ;
   m_nNy[0] = max( int( ( dDimY + EPS_SMALL) / m_dStep + 0.5), 1) ;

  // Numero di componenti connesse
   m_nConnectedCompoCount = 1 ;

  // Se tridexel
   if ( bTriDex) {
      m_nNx[1] = m_nNy[0] ;
      m_nNy[1] = max( int( ( dDimZ + EPS_SMALL) / m_dStep + 0.5), 1) ;
      m_nNx[2] = m_nNy[1] ;
      m_nNy[2] = m_nNx[0] ;
   }

  // altrimenti mono dexel
   else {
      m_nNx[1] = 0 ;
      m_nNy[1] = 0 ;
      m_nNx[2] = 0 ;
      m_nNy[2] = 0 ;
   }

   // Definisco il numero di blocchi lungo x,y e z
   if ( ! CalcBlockNum())
      return false ;

  // Definizione della mappa

  // Creazione delle mappe
  // Calcolo del numero di celle per ogni mappa
   for ( int i = 0 ; i < m_nMapNum ; ++ i)
      m_nDim[i] = m_nNx[i] * m_nNy[i] ;

  // Se versione 32-bit controllo di non superare il numero di Dexel massimo
   #if !defined(_WIN64)
      for ( int i = 0 ; i < ssize( m_nDim) ; ++ i)
         m_nDexelNbr += m_nDim[i] ;
      if ( m_nDexelNbr >= MAX_DEXEL_32_BIT) {
         Clear() ;
         if ( nError != nullptr)
            *nError = 1 ;
         return false ;
      }
   #endif

  // Creazione delle celle per ogni mappa
   for ( int i = 0 ; i < m_nMapNum ; ++ i)
      m_Values[i].resize( m_nDim[i]) ;

  // Riempimento delle celle
   for ( int i = 0 ; i < m_nMapNum ; ++ i) {
      for ( int j = 0 ; j < m_nDim[i] ; ++ j) {

        // Aggiungo il tratto al dexel vuoto
         m_Values[i][j].resize( 1) ;
         m_Values[i][j][0].dMin = 0 ;
         m_Values[i][j][0].nToolMin = 0 ;
         m_Values[i][j][0].nCompo = 1 ;

         switch ( i) {
         case 0 :
            m_Values[i][j][0].vtMinN = - Z_AX ;
            m_Values[i][j][0].dMax = dDimZ ;
            m_Values[i][j][0].vtMaxN = Z_AX ;
            m_Values[i][j][0].nToolMax = 0 ;
            break ;
         case 1 :
            m_Values[i][j][0].vtMinN = - X_AX ;
            m_Values[i][j][0].dMax = dDimX ;
            m_Values[i][j][0].vtMaxN = X_AX ;
            m_Values[i][j][0].nToolMax = 0 ;
            break ;
         case 2 :
            m_Values[i][j][0].vtMinN = - Y_AX ;
            m_Values[i][j][0].dMax = dDimY ;
            m_Values[i][j][0].vtMaxN = Y_AX ;
            m_Values[i][j][0].nToolMax = 0 ;
            break ;
         }
      }
   }

  // Definizione delle limitazioni iniziali in Z per ogni mappa
   m_dMinZ[0] = 0 ;
   m_dMaxZ[0] = dDimZ ;
   m_dMinZ[1] = 0 ;
   m_dMaxZ[1] = ( bTriDex ? dDimX : 0) ;
   m_dMinZ[2] = 0 ;  
   m_dMaxZ[2] = ( bTriDex ? dDimY : 0) ;

   // Tipologia
   m_nShape = BOX ;

  // Aggiornamento dello stato
   m_nStatus = OK ;

   return true ;
}

//----------------------------------------------------------------------------
bool 
VolZmap::CreateEmpty( const Point3d& ptO, double dDimX, double dDimY, double dDimZ, double dStep, bool bTriDex, int* nError)
{
  // Controlli sull'ammissibilità delle dimensioni lineari del grezzo e del passo
   if ( dStep < EPS_SMALL  ||  dDimX < EPS_SMALL  ||  dDimY < EPS_SMALL  ||  dDimZ < EPS_SMALL)
      return false ;

  // Il passo di discretizzazione non può essere inferiore a 100 * EPS_SMALL
   m_dStep = max( dStep, 100 * EPS_SMALL) ;

  // Aggiorno la dimensione della mappa 1 o 3
   m_nMapNum = ( bTriDex ? 3 : 1) ;

  // Disponendo i sistemi di riferimento in una successione, le coordinate x,y,z
  // di uno si ottengono da una permutazione ciclica di quelle del precedente sistema.
  // es: X(n) = Z(n-1), Y(n) = X(n-1), Z(n) = Y(n-1)

  // Definisco il sistema di riferimento intrinseco
   m_MapFrame.Set( ptO, X_AX, Y_AX, Z_AX) ;

  // Definisco i vettori dei limiti su indici
   m_nNx[0] = max( int( ( dDimX + EPS_SMALL) / m_dStep + 0.5), 1) ;
   m_nNy[0] = max( int( ( dDimY + EPS_SMALL) / m_dStep + 0.5), 1) ;

  // Numero di componenti connesse 
   m_nConnectedCompoCount = 1 ;

  // Se tridexel
   if ( bTriDex) {
      m_nNx[1] = m_nNy[0] ;
      m_nNy[1] = max( int( ( dDimZ + EPS_SMALL) / m_dStep + 0.5), 1) ;
      m_nNx[2] = m_nNy[1] ;
      m_nNy[2] = m_nNx[0] ;
   }

  // altrimenti mono dexel
   else {
      m_nNx[1] = 0 ;
      m_nNy[1] = 0 ;
      m_nNx[2] = 0 ;
      m_nNy[2] = 0 ; 
   }

   // Definisco il numero di blocchi lungo x,y e z
   if ( ! CalcBlockNum())
      return false ;

  // Creazione delle mappe
  // Calcolo del numero di celle per ogni mappa
   for ( int i = 0 ; i < m_nMapNum ; ++ i)
      m_nDim[i] = m_nNx[i] * m_nNy[i] ;

  // Se versione 32-bit controllo di non superare il numero di Dexel massimo
   #if !defined(_WIN64)
      for ( int i = 0 ; i < ssize( m_nDim) ; ++ i)
         m_nDexelNbr += m_nDim[i] ;
      if ( m_nDexelNbr >= MAX_DEXEL_32_BIT) {
         Clear() ;
         if ( nError != nullptr)
            *nError = 1 ;
         return false ;
      }
   #endif

  // Creazione delle celle per ogni mappa
   for ( int i = 0 ; i < m_nMapNum ; ++ i)
      m_Values[i].resize( m_nDim[i]) ;

  // Definizione delle limitazioni iniziali in Z per ogni mappa
   m_dMinZ[0] = 0 ;
   m_dMaxZ[0] = dDimZ ;
   m_dMinZ[1] = 0 ;
   m_dMaxZ[1] = ( bTriDex ? dDimX : 0) ;
   m_dMinZ[2] = 0 ;  
   m_dMaxZ[2] = ( bTriDex ? dDimY : 0) ;

   // Tipologia
   m_nShape = GENERIC ;

  // Aggiornamento dello stato
   m_nStatus = OK ;

   return true ;
} 

//----------------------------------------------------------------------------
bool
VolZmap::CreateFromFlatRegion( const ISurfFlatRegion& Surf, double dDimZ, double dStep, bool bTriDex, int* nError)
{
  // Aggiorno la dimensione della mappa 1 o 3
   m_nMapNum = ( bTriDex ? 3 : 1) ;

  // Il passo di discretizzazione non può essere inferiore a 100 * EPS_SMALL
   m_dStep = max( dStep, 100 * EPS_SMALL) ;
 
  // Determino il bounding box della flat region
   BBox3d SurfBBox ;
   Surf.GetLocalBBox( SurfBBox, BBF_EXACT) ;

  // Determino i punti estremi del bounding box
   Point3d ptMapOrig, ptMapEnd ; 
   SurfBBox.GetMinMax( ptMapOrig, ptMapEnd) ;
   SurfBBox.Expand( 100 * EPS_SMALL, 100 * EPS_SMALL, 0) ;

  // Sistema di riferimento intrinseco dello Zmap
   m_MapFrame.Set( ptMapOrig, X_AX, Y_AX, Z_AX) ;

  // Determino le dimensioni lineari X Y della griglia
   double dLengthX = ptMapEnd.x - ptMapOrig.x ;
   double dLengthY = ptMapEnd.y - ptMapOrig.y ;

  // A partire dalle dimensioni di xy del grezzo determino il numero di colonne e righe
  // della griglia Zmap e da questi la dimensione del vettore di dexel
   m_nNx[0] = int( ( dLengthX + EPS_SMALL) / m_dStep + 0.5) ;
   m_nNy[0] = int( ( dLengthY + EPS_SMALL) / m_dStep + 0.5) ;
   m_nDim[0] = m_nNx[0] * m_nNy[0] ;
  // Ridimensiono il vettore di dexel e creo lo Zmap
   m_Values[0].resize( m_nDim[0]) ;

  // Numero di componenti connesse 
   m_nConnectedCompoCount = Surf.GetChunkCount() ;

  // Se Tridexel ridimensiono anche gli altri vettori
   if ( bTriDex) {
      m_nNx[1] = m_nNy[0] ;
      m_nNy[1] = int( ( dDimZ + EPS_SMALL) / m_dStep + 0.5) ;
      m_nDim[1] = m_nNx[1] * m_nNy[1] ;
      m_nNx[2] = m_nNy[1] ;
      m_nNy[2] = m_nNx[0] ;
      m_nDim[2] = m_nNx[2] * m_nNy[2] ;
     // Se versione 32-bit controllo di non superare il numero di Dexel massimo
      #if !defined(_WIN64)
         for ( int i = 0 ; i < ssize( m_nDim) ; ++ i)
            m_nDexelNbr += m_nDim[i] ;
         if ( m_nDexelNbr >= MAX_DEXEL_32_BIT) {
            Clear() ;
            if ( nError != nullptr)
               *nError = 1 ;
            return false ;
         }
      #endif
      m_Values[1].resize( m_nDim[1]) ;
      m_Values[2].resize( m_nDim[2]) ;
   }
  // Se dimensione singola
   else {
     // Se versione 32-bit controllo di non superare il numero di Dexel massimo
      #if !defined(_WIN64)
         m_nDexelNbr += m_nDim[0] ;
         if ( m_nDexelNbr >= MAX_DEXEL_32_BIT) {
            Clear() ;
            if ( nError != nullptr)
               *nError = 1 ;
            return false ;
         }
      #endif
      m_nNx[1] = 0 ;
      m_nNy[1] = 0 ;
      m_nDim[1] = 0 ;
      m_nNx[2] = 0 ;
      m_nNy[2] = 0 ;
      m_nDim[2] = 0 ;
   }

  // Definisco il numero di blocchi lungo x,y e z
   if ( ! CalcBlockNum())
      return false ;

  // Metto in cache le curve di contorno della regione
   ICURVEPOVECTOR vpCrvs ;
   INTVECTOR vnCompo ;
   int nChunkNum = Surf.GetChunkCount() ;
   for ( int nChunk = 0 ; nChunk < nChunkNum ; ++ nChunk) {
      int nLoopNum = Surf.GetLoopCount( nChunk) ;
      for ( int nLoop = 0 ; nLoop < nLoopNum ; ++ nLoop) { 
         ICurve* pCrv = Surf.GetLoop( nChunk, nLoop) ;
         if ( pCrv != nullptr) {
            vpCrvs.emplace_back( pCrv) ;
            vnCompo.emplace_back( nChunk + 1) ;
         }
      }
   }

  // Calcolo griglia 0=XY ( se tridexel anche griglia 2=ZX)
   for ( int i = 0 ; i < m_nNx[0] ; ++ i) {

     // Definisco la retta diretta come Y da intersecare con la regione
      double dX = ( i + 0.5) * m_dStep ;
      Point3d ptP0 = ptMapOrig + Vector3d( dX, 0, 0) ;
      CurveLine GridLine ;
      GridLine.SetPVL( ptP0, Y_AX, dLengthY) ;

     // Determino le intersezioni della retta con la regione
      CRVCVECTOR IntersectionResults ;
      Surf.GetCurveClassification( GridLine, EPS_SMALL, IntersectionResults) ;

     // Analizzo le parti in cui la retta è stata divisa
      int nPart = int( IntersectionResults.size()) ;
      for ( int k = 0 ; k < nPart ; ++ k) {

        // Se la retta è interna alla regione o coincidente con parte della sua frontiera
         int nType = IntersectionResults[k].nClass ;
         if ( nType == CRVC_IN  ||  nType == CRVC_ON_P  ||  nType == CRVC_ON_M) {

           // Lunghezze dei tratti di retta corrente
            double dLen1 = IntersectionResults[k].dParS * dLengthY ;
            double dLen2 = IntersectionResults[k].dParE * dLengthY ;

           // Punti estremi della parte di retta corrente
            Point3d ptP1 = ptP0 + dLen1 * Y_AX ;
            Point3d ptP2 = ptP0 + dLen2 * Y_AX ;

           // Ricerca di questi punti sui contorni della regione, per avere le normali e il numero di componente connesso
            int nFind = 0 ;
            int nCompo = 0 ;
            Vector3d vtN1 = - Y_AX ; Vector3d vtN2 = Y_AX ; 
            for ( int m = 0 ; m < int( vpCrvs.size()) ; ++ m) {
              // recupero la curva
               ICurve* pCurve = vpCrvs[m] ;
              // determino posizione primo punto su curva
               double dP1 ;
               if ( ( nFind & 1) == 0  &&  pCurve->GetParamAtPoint( ptP1, dP1, 10 * EPS_SMALL)) {
                  Point3d ptTemp1 ;
                  Vector3d vtT1 ;
                  pCurve->GetPointTang( dP1, ICurve::FROM_MINUS, ptTemp1, vtT1) ;
                  vtN1 = vtT1 ^ Z_AX ;
                  nFind += 1 ;
               }
              // determino posizione secondo punto su curva
               double dP2 ;
               if ( ( nFind & 2) == 0  &&  pCurve->GetParamAtPoint( ptP2, dP2, 10 * EPS_SMALL)) {
                  Point3d ptTemp2 ;
                  Vector3d vtT2 ;
                  pCurve->GetPointTang( dP2, ICurve::FROM_MINUS, ptTemp2, vtT2) ;
                  vtN2 = vtT2 ^ Z_AX ;
                  nFind += 2 ;
               }
              // Se trovati entrambi gli estremi, esco dal ciclo
               if ( nFind == 3) {
                  nCompo = vnCompo[m] ;
                  break ;
               }
            }

           // Verifico di aver trovato i punti sulle curve
            if ( nFind != 3)
               LOG_ERROR( GetEGkLogger(), "Error in VolZmap::CreateFromFlatRegion : point not on baundary")

           // Ridimensiono e riempio i dexel della griglia 0
            int nStartJ = Clamp( int( floor( dLen1 / m_dStep - EPS_SMALL + 0.5)), 0, m_nNy[0] - 1) ;
            int nEndJ = Clamp( int( floor( dLen2 / m_dStep + EPS_SMALL - 0.5)), 0, m_nNy[0] - 1) ;
            for ( int j = nStartJ ; j <= nEndJ ; ++ j) {
              // Determino il dexel
               int nPos0 = j * m_nNx[0] + i ;
              // Aggiungo il tratto al dexel vuoto
               m_Values[0][nPos0].resize( 1) ;
              // Aggiorno i dati del tratto di dexel
               m_Values[0][nPos0][0].dMin = 0 ;
               m_Values[0][nPos0][0].vtMinN = - Z_AX ;
               m_Values[0][nPos0][0].nToolMin = 0 ;
               m_Values[0][nPos0][0].dMax = dDimZ ;
               m_Values[0][nPos0][0].vtMaxN = Z_AX ;
               m_Values[0][nPos0][0].nToolMax = 0 ;
               m_Values[0][nPos0][0].nCompo = nCompo ;
            }

           // Se tridexel riempio i singoli dexel della griglia 2 con gli intervalli
            if ( bTriDex) { 
               for ( int n = 0 ; n < m_nNx[2] ; ++ n) {
                  int nPos2 = i * m_nNx[2] + n ; 
                  int nCurrSize = int( m_Values[2][nPos2].size()) ;
                 // Aggiungo un tratto al dexel
                  m_Values[2][nPos2].resize( nCurrSize + 1) ;
                 // Aggiorno i dati del tratto di dexel
                  m_Values[2][nPos2][nCurrSize].dMin = dLen1 ;
                  m_Values[2][nPos2][nCurrSize].vtMinN = vtN1 ;
                  m_Values[2][nPos2][nCurrSize].nToolMin = 0 ;
                  m_Values[2][nPos2][nCurrSize].dMax = dLen2 ;
                  m_Values[2][nPos2][nCurrSize].vtMaxN = vtN2 ;
                  m_Values[2][nPos2][nCurrSize].nToolMax = 0 ;
                  m_Values[2][nPos2][nCurrSize].nCompo = nCompo ;
               }
            }
         }
      }
   }

  // Se tridexel calcolo griglia 1=YZ
   if ( bTriDex) {

     // ciclo sul lato orizzontale della griglia
      for ( int i = 0 ; i < m_nNx[1] ; ++ i) {

        // Definisco la retta diretta come X da intersecare con la regione
         double dY = ( i + 0.5) * m_dStep ;
         Point3d ptP0 = ptMapOrig + Vector3d( 0, dY, 0) ; 
         CurveLine GridLine ;
         GridLine.SetPVL( ptP0, X_AX, dLengthX) ;  

        // Determino le intersezioni della retta con la regione
         CRVCVECTOR IntersectionResults ;
         Surf.GetCurveClassification( GridLine, EPS_SMALL, IntersectionResults) ;

        // Analizzo le parti
         int nPart = int( IntersectionResults.size()) ; 
         for ( int k = 0 ; k < nPart ; ++ k) {

           // Se la retta è interna alla regione o coincidente con parte della sua frontiera
            int nType = IntersectionResults[k].nClass ;
            if ( nType == CRVC_IN  ||  nType == CRVC_ON_P  ||  nType == CRVC_ON_M) {

              // Lunghezze dei tratti di retta
               double dLen1 = IntersectionResults[k].dParS * dLengthX ;
               double dLen2 = IntersectionResults[k].dParE * dLengthX ;

              // Punti estremi
               Point3d ptP1 = ptP0 + dLen1 * X_AX ;
               Point3d ptP2 = ptP0 + dLen2 * X_AX ;

              // Ricerca di questi punti sui contorni della regione, per avere le normali e il numero di componente connesso
               int nFind = 0 ;
               int nCompo = 0 ;
               Vector3d vtN1 = -X_AX ; Vector3d vtN2 = X_AX ; 
               for ( int m = 0 ; m < int( vpCrvs.size()) ; ++ m) {
                 // recupero la curva
                  ICurve* pCurve = vpCrvs[m] ;
                 // determino posizione primo punto su curva
                  double dP1 ;
                  if ( ( nFind & 1) == 0  &&  pCurve->GetParamAtPoint( ptP1, dP1, 10 * EPS_SMALL)) {
                     Point3d ptTemp1 ; Vector3d vtT1 ;
                     pCurve->GetPointTang( dP1, ICurve::FROM_MINUS, ptTemp1, vtT1) ;
                     vtN1 = vtT1 ^ Z_AX ;
                     nFind += 1 ;
                  }
                 // determino posizione secondo punto su curva
                  double dP2 ;
                  if ( ( nFind & 2) == 0  &&  pCurve->GetParamAtPoint( ptP2, dP2, 10 * EPS_SMALL)) {
                     Point3d ptTemp2 ; Vector3d vtT2 ;
                     pCurve->GetPointTang( dP2, ICurve::FROM_MINUS, ptTemp2, vtT2) ;
                     vtN2 = vtT2 ^ Z_AX ;
                     nFind += 2 ;
                  } 

                 // Se trovati entrambi gli estremi, esco dal ciclo
                  if ( nFind == 3) {
                     nCompo = vnCompo[m] ;
                     break ;
                  }
               }

              // Verifico di aver trovato i punti sulle curve
               if ( nFind != 3)
                  LOG_ERROR( GetEGkLogger(), "Error in VolZmap::CreateFromFlatRegion : point not on baundary")

              // aggiorno i dexel impilati
               for ( int j = 0 ; j < m_nNy[1] ; ++ j) {
                  int nPos1 = j * m_nNx[1] + i ;
                  int nCurrSize = int( m_Values[1][nPos1].size()) ;
                 // Aggiungo un tratto al dexel
                  m_Values[1][nPos1].resize( nCurrSize + 1) ;
                 // Assegno i dati
                  m_Values[1][nPos1][nCurrSize].dMin = dLen1 ;
                  m_Values[1][nPos1][nCurrSize].vtMinN = vtN1 ; 
                  m_Values[1][nPos1][nCurrSize].nToolMin = 0 ;
                  m_Values[1][nPos1][nCurrSize].dMax = dLen2 ;
                  m_Values[1][nPos1][nCurrSize].vtMaxN = vtN2 ;
                  m_Values[1][nPos1][nCurrSize].nToolMax = 0 ;
                  m_Values[1][nPos1][nCurrSize].nCompo = nCompo ;
               }
            }
         }
      }
   }

  // Definizione delle limitazioni iniziali in Z per ogni mappa
   m_dMinZ[0] = 0 ;
   m_dMaxZ[0] = dDimZ ;
   m_dMinZ[1] = 0 ;
   m_dMaxZ[1] = ( bTriDex ? dLengthX : 0) ;
   m_dMinZ[2] = 0 ;  
   m_dMaxZ[2] = ( bTriDex ? dLengthY : 0) ;

  // Tipologia
   m_nShape = ( IsBox() ? BOX : EXTRUSION) ;

  // Aggiornamento dello stato
   m_nStatus = OK ;

   return true ;
}

//----------------------------------------------------------------------------
bool
VolZmap::CreateMapPart( int nMap, int nInfI, int nSupI, int nInfJ, int nSupJ, const Vector3d& vtLen, const Point3d& ptMapOrig,
                        const ISurfTriMesh& Surf, IntersParLinesSurfTm& intPLSTM, bool bForceClosed)
{
   if ( nMap < 0  ||  nMap > 2  ||
        nInfI < 0  ||  nInfI > m_nNx[nMap]  ||
        nSupI < 0  ||  nSupI > m_nNx[nMap]  ||
        nInfJ < 0  ||  nInfJ > m_nNy[nMap]  ||
        nSupJ < 0  ||  nSupJ > m_nNy[nMap])
      return false ;

   double dCosSmall = sin( EPS_ANG_SMALL * DEGTORAD) ;

   INTTOCHECKVEC vIntToCheck ;

   // Determinazione e ridimensionamento dei dexel interni alla trimesh
   for ( int i = nInfI ; i < nSupI ; ++ i) {
      for ( int j = nInfJ ; j < nSupJ ; ++ j) {

         // Definisco la retta da intersecare con la trimesh
         double dX = ( i + 0.5) * m_dStep ;
         double dY = ( j + 0.5) * m_dStep ;
         Point3d ptP0( dX, dY, 0) ;

         // Determino le intersezioni della retta con la TriMesh
         ILSIVECTOR IntersectionResults ;
         intPLSTM.GetInters( ptP0, vtLen.v[(nMap+2)%3], IntersectionResults) ;

         for ( int nI = 0 ; nI < int( IntersectionResults.size()) - 3 ; ++ nI) {
            int nJ = nI + 1 ;  
            int nK = nJ + 1 ;
            int nT = nK + 1 ;
            int nSgnI = IntersectionResults[nI].dCosDN > EPS_SMALL ? 1 : IntersectionResults[nI].dCosDN > -EPS_SMALL ? 0 : - 1 ;
            int nSgnJ = IntersectionResults[nJ].dCosDN > EPS_SMALL ? 1 : IntersectionResults[nJ].dCosDN > -EPS_SMALL ? 0 : - 1 ;
            int nSgnK = IntersectionResults[nK].dCosDN > EPS_SMALL ? 1 : IntersectionResults[nK].dCosDN > -EPS_SMALL ? 0 : - 1 ;
            int nSgnT = IntersectionResults[nT].dCosDN > EPS_SMALL ? 1 : IntersectionResults[nT].dCosDN > -EPS_SMALL ? 0 : - 1 ;
            double dUJ = IntersectionResults[nJ].dU ;
            double dUK = IntersectionResults[nK].dU ;
            if ( nSgnI != 0  &&  nSgnI == nSgnJ  &&  nSgnK != 0  &&  nSgnK == nSgnT  &&  nSgnI == - nSgnT  &&  abs( dUJ - dUK) < EPS_SMALL) {
               IntersectionResults.erase( IntersectionResults.begin() + nK) ;
               IntersectionResults.erase( IntersectionResults.begin() + nJ) ;
            }
         }

         int nInt = int( IntersectionResults.size()) ;

         if ( nInt == 1) {
            int nT = IntersectionResults[0].nT ;
            int nF = Surf.GetFacetFromTria( nT) ;
            Vector3d vtN ; Surf.GetFacetNormal( nF, vtN) ;
            vIntToCheck.emplace_back( IntersectionResults[0].ptI, vtN, IntersectionResults[0].dU, nMap, i, j) ;
            continue ;
         }

         int nPos = j * m_nNx[nMap] + i ;

         bool bInside = false ;
         Point3d ptIn ;
         Vector3d vtInN ;

         for ( int k = 0 ; k < nInt ; ++ k) {

            if ( k > 0) {
               int z = k - 1 ;
               // controllo se l'intersezione corrente è coincidente con la precedente, per esempio se ILTT == 4 ( intersezione su un lato di un triangolo, quindi con 2 triangoli)
               // controllo anche che le proiezioni delle normali abbiano lo stesso segno rispetto allo spillone
               if ( IntersectionResults[k].dU - IntersectionResults[z].dU < EPS_SMALL  &&
                  IntersectionResults[k].dCosDN - IntersectionResults[z].dCosDN < EPS_SMALL  &&
                  IntersectionResults[k].dCosDN * IntersectionResults[z].dCosDN > 0)
                  continue ;
            }

            int nIntType = IntersectionResults[k].nILTT ;

            // Se c'è intersezione
            if ( nIntType != ILTT_NO) {

               double dCos = IntersectionResults[k].dCosDN ;

               // entro nella superficie trimesh
               if ( dCos < - dCosSmall) {
                  if ( bInside) {
                     // salvo la precedente, perché vuol dire che manca la sua uscita
                     // poi procedo con l'intersezione corrente
                     int z = k - 1 ;
                     vIntToCheck.emplace_back( IntersectionResults[z].ptI, vtInN, IntersectionResults[z].dU, nMap, i, j) ;
                  }
                  ptIn = IntersectionResults[k].ptI ;

                  int nT = IntersectionResults[k].nT ;
                  int nF = Surf.GetFacetFromTria( nT) ;

                  Surf.GetFacetNormal( nF, vtInN) ;

                  bInside = true ;
                  if ( k == nInt - 1) {
                     // se ho un ingresso ma è l'ultima intersezione vuol dire che è spaiata
                     vIntToCheck.emplace_back( IntersectionResults[k].ptI, vtInN, IntersectionResults[k].dU, nMap, i, j) ;
                  }
               }

               // esco dalla superficie trimesh
               else if ( dCos > dCosSmall) {
                  Point3d ptOut = IntersectionResults[k].ptI ;

                  int nT = IntersectionResults[k].nT ;
                  int nF = Surf.GetFacetFromTria( nT) ;

                  Vector3d vtOutN ;
                  Surf.GetFacetNormal( nF, vtOutN) ; 

                  if ( ! bInside) {
                     // se l'intersezione è uscente, ma non ho un'entrata allora devo salvare l'intersezione corrente, perché è andata persa la sua entrata
                     // procedo poi con la prossima intersezione
                     vIntToCheck.emplace_back( IntersectionResults[k].ptI, vtOutN, IntersectionResults[k].dU, nMap, i, j) ;
                     continue ;
                  }

                  int nCurrentSize = int( m_Values[nMap][nPos].size()) ;

                  // Aggiungo un tratto al dexel
                  m_Values[nMap][nPos].resize( nCurrentSize + 1) ;

                  // Aggiorno dati del tratto di dexel 
                  m_Values[nMap][nPos][nCurrentSize].dMin = ptIn.v[(nMap+2)%3] - ptMapOrig.v[(nMap+2)%3] ;
                  m_Values[nMap][nPos][nCurrentSize].dMax = ptOut.v[(nMap+2)%3] - ptMapOrig.v[(nMap+2)%3] ;
                  m_Values[nMap][nPos][nCurrentSize].vtMinN = vtInN ;
                  m_Values[nMap][nPos][nCurrentSize].vtMaxN = vtOutN ;
                  m_Values[nMap][nPos][nCurrentSize].nToolMin = 0 ;
                  m_Values[nMap][nPos][nCurrentSize].nToolMax = 0 ;
                  m_Values[nMap][nPos][nCurrentSize].nCompo = 0 ;

                  bInside = false ;
               }
            }   
         }
      }
   }

   if ( bForceClosed) {
      if ( ! AdjustDexelThroughCracks( nMap, vtLen, vIntToCheck, ptMapOrig, Surf, intPLSTM))
         return false ;
      vIntToCheck.clear() ;
      SearchMissingIntervals( nMap, vIntToCheck) ;
      if ( ! AdjustDexelThroughCracks( nMap, vtLen, vIntToCheck, ptMapOrig, Surf, intPLSTM))
         return false ;
   }

   return true ;
}

//----------------------------------------------------------------------------
bool
VolZmap::AddMapPart( int nMap, int nInfI, int nSupI, int nInfJ, int nSupJ, const Vector3d& vtLen, const Point3d& ptMapOrig,
                     const ISurfTriMesh& Surf, IntersParLinesSurfTm& intPLSTM)
{
  // controllo sui parametri
   if ( nMap < 0  ||  nMap > 2  ||
        nInfI < 0  ||  nInfI > m_nNx[nMap]  ||
        nSupI < 0  ||  nSupI > m_nNx[nMap]  ||
        nInfJ < 0  ||  nInfJ > m_nNy[nMap]  ||
        nSupJ < 0  ||  nSupJ > m_nNy[nMap])
      return false ;

  // determinazione e ridimensionamento dei dexel interni alla trimesh
   for ( int i = nInfI ; i < nSupI ; ++ i) {
      for ( int j = nInfJ ; j < nSupJ ; ++ j) {

        // definisco la retta da intersecare con la trimesh
         double dX = ( i + 0.5) * m_dStep ;
         double dY = ( j + 0.5) * m_dStep ;
         Point3d ptP0( dX, dY, 0) ;

        // intersezioni della retta con la TriMesh
         ILSIVECTOR IntersectionResults ;
         intPLSTM.GetInters( ptP0, vtLen.v[(nMap+2)%3], IntersectionResults) ;

        // rimuovo le intersezioni in eccesso
         for ( int nI = 0 ; nI < int( IntersectionResults.size()) - 3 ; ++ nI) {
            int nJ = nI + 1 ; // prima successiva
            int nK = nJ + 1 ; // seconda successiva
            int nT = nK + 1 ; // terza successiva
           // determino i segni delle 4 intersezioni tra la linea e il trangolo della TriMesh
            int nSgnI = IntersectionResults[nI].dCosDN > EPS_SMALL ? 1 : IntersectionResults[nI].dCosDN > -EPS_SMALL ? 0 : - 1 ;
            int nSgnJ = IntersectionResults[nJ].dCosDN > EPS_SMALL ? 1 : IntersectionResults[nJ].dCosDN > -EPS_SMALL ? 0 : - 1 ;
            int nSgnK = IntersectionResults[nK].dCosDN > EPS_SMALL ? 1 : IntersectionResults[nK].dCosDN > -EPS_SMALL ? 0 : - 1 ;
            int nSgnT = IntersectionResults[nT].dCosDN > EPS_SMALL ? 1 : IntersectionResults[nT].dCosDN > -EPS_SMALL ? 0 : - 1 ;
           // parametri dell'intersezione sulla linea
            double dUJ = IntersectionResults[nJ].dU ;
            double dUK = IntersectionResults[nK].dU ;
           // controllo coerenza con segni...
            if ( nSgnI != 0  &&  nSgnI == nSgnJ  &&
                 nSgnK != 0  &&  nSgnK == nSgnT  &&
                 nSgnI == - nSgnT  &&
                 abs( dUJ - dUK) < EPS_SMALL) {
              // ... ed elimino le intersezioni in eccesso...
               IntersectionResults.erase( IntersectionResults.begin() + nK) ;
               IntersectionResults.erase( IntersectionResults.begin() + nJ) ;
            }
         }

         int nInt = int( IntersectionResults.size()) ; // numero di intersezioni valide
         bool bInside = false ; // Flag entrata/uscita per tratto di retta
         Point3d ptIn ; Vector3d vtInN ;

        // per ogni intersezione valida trovata...
         for ( int k = 0 ; k < nInt ; ++ k) {
           // ricavo il tipo di intersezione
            int nIntType = IntersectionResults[k].nILTT ;
           // se c'è intersezione
            if ( nIntType != ILTT_NO) {
              // ricavo il cos tra i vettori ( normale del triangolo e tangente alla retta)
               double dCos = IntersectionResults[k].dCosDN ;

              // se entro nella superficie trimesh...
               if ( dCos < - EPS_SMALL) {
                  ptIn = IntersectionResults[k].ptI ; // punto di intersezione
                  int nT = IntersectionResults[k].nT ; // triangolo di interesse
                  int nF = Surf.GetFacetFromTria( nT) ; // faccia di interesse
                  Surf.GetFacetNormal( nF, vtInN) ;
                  bInside = true ; // entrata
               }
              // ...se esco dalla superficie trimesh ( prima sono per forza entrato)
               else if ( dCos > EPS_SMALL  &&  bInside) {
                  Point3d ptOut = IntersectionResults[k].ptI ; // punto di intersezione
                  int nT = IntersectionResults[k].nT ;   // triangolo di interesse
                  int nF = Surf.GetFacetFromTria( nT) ;  // faccia di interesse
                  Vector3d vtOutN ; Surf.GetFacetNormal( nF, vtOutN) ; // vettore d'uscita

                 // Aggiungo un tratto al dexel
                  AddIntervals( nMap, i, j,
                                ptIn.v[(nMap+2)%3] - ptMapOrig.v[(nMap+2)%3],
                                ptOut.v[(nMap+2)%3] - ptMapOrig.v[(nMap+2)%3],
                                vtInN, vtOutN, 0, true) ;
                  bInside = false ; // uscita
               }
            }
         }
      }
   }

   return true ;
}

//----------------------------------------------------------------------------
bool
VolZmap::SubtractMapPart( int nMap, int nInfI, int nSupI, int nInfJ, int nSupJ, const Vector3d& vtLen, const Point3d& ptMapOrig,
                          const ISurfTriMesh& Surf, IntersParLinesSurfTm& intPLSTM)
{
  // controllo sui parametri
   if ( nMap < 0  ||  nMap > 2  ||
        nInfI < 0  ||  nInfI > m_nNx[nMap]  ||
        nSupI < 0  ||  nSupI > m_nNx[nMap]  ||
        nInfJ < 0  ||  nInfJ > m_nNy[nMap]  ||
        nSupJ < 0  ||  nSupJ > m_nNy[nMap])
      return false ;

  // determinazione e ridimensionamento dei dexel interni alla trimesh
   for ( int i = nInfI ; i < nSupI ; ++ i) {
      for ( int j = nInfJ ; j < nSupJ ; ++ j) {

        // definisco la retta da intersecare con la trimesh
         double dX = ( i + 0.5) * m_dStep ;
         double dY = ( j + 0.5) * m_dStep ;
         Point3d ptP0( dX, dY, 0) ;

        // intersezioni della retta con la TriMesh
         ILSIVECTOR IntersectionResults ;
         intPLSTM.GetInters( ptP0, vtLen.v[(nMap+2)%3], IntersectionResults) ;

        // rimuovo le intersezioni in eccesso
         for ( int nI = 0 ; nI < int( IntersectionResults.size()) - 3 ; ++ nI) {
            int nJ = nI + 1 ; // prima successiva
            int nK = nJ + 1 ; // seconda successiva
            int nT = nK + 1 ; // terza successiva
           // determino i segni delle 4 intersezioni tra la linea e il trangolo della TriMesh
            int nSgnI = IntersectionResults[nI].dCosDN > EPS_SMALL ? 1 : IntersectionResults[nI].dCosDN > -EPS_SMALL ? 0 : - 1 ;
            int nSgnJ = IntersectionResults[nJ].dCosDN > EPS_SMALL ? 1 : IntersectionResults[nJ].dCosDN > -EPS_SMALL ? 0 : - 1 ;
            int nSgnK = IntersectionResults[nK].dCosDN > EPS_SMALL ? 1 : IntersectionResults[nK].dCosDN > -EPS_SMALL ? 0 : - 1 ;
            int nSgnT = IntersectionResults[nT].dCosDN > EPS_SMALL ? 1 : IntersectionResults[nT].dCosDN > -EPS_SMALL ? 0 : - 1 ;
           // parametri dell'intersezione sulla linea
            double dUJ = IntersectionResults[nJ].dU ;
            double dUK = IntersectionResults[nK].dU ;
           // controllo coerenza con segni...
            if ( nSgnI != 0  &&  nSgnI == nSgnJ  &&
                 nSgnK != 0  &&  nSgnK == nSgnT  &&
                 nSgnI == - nSgnT  &&
                 abs( dUJ - dUK) < EPS_SMALL) {
              // ... ed elimino le intersezioni in eccesso...
               IntersectionResults.erase( IntersectionResults.begin() + nK) ;
               IntersectionResults.erase( IntersectionResults.begin() + nJ) ;
            }
         }

         int nInt = int( IntersectionResults.size()) ; // numero di intersezioni valide
         bool bInside = false ; // Flag entrata/uscita per tratto di retta
         Point3d ptIn ; Vector3d vtInN ;

        // per ogni intersezione valida trovata...
         for ( int k = 0 ; k < nInt ; ++ k) {
           // ricavo il tipo di intersezione
            int nIntType = IntersectionResults[k].nILTT ;
           // se c'è intersezione
            if ( nIntType != ILTT_NO) {
              // ricavo il cos tra i vettori ( normale del triangolo e tangente alla retta)
               double dCos = IntersectionResults[k].dCosDN ;

              // se entro nella superficie trimesh...
               if ( dCos < - EPS_SMALL) {
                  ptIn = IntersectionResults[k].ptI ; // punto di intersezione
                  int nT = IntersectionResults[k].nT ; // triangolo di interesse
                  int nF = Surf.GetFacetFromTria( nT) ; // faccia di interesse
                  Surf.GetFacetNormal( nF, vtInN) ;
                  bInside = true ; // entrata
               }
              // ...se esco dalla superficie trimesh ( prima sono per forza entrato)
               else if ( dCos > EPS_SMALL  &&  bInside) {
                  Point3d ptOut = IntersectionResults[k].ptI ; // punto di intersezione
                  int nT = IntersectionResults[k].nT ;   // triangolo di interesse
                  int nF = Surf.GetFacetFromTria( nT) ;  // faccia di interesse
                  Vector3d vtOutN ; Surf.GetFacetNormal( nF, vtOutN) ; // vettore d'uscita

                 // Aggiungo un tratto al dexel
                  SubtractIntervals( nMap, i, j,
                                     ptIn.v[(nMap+2)%3] - ptMapOrig.v[(nMap+2)%3],
                                     ptOut.v[(nMap+2)%3] - ptMapOrig.v[(nMap+2)%3],
                                     - vtInN, - vtOutN, 0, true) ;
                  bInside = false ; // uscita
               }
            }
         }
      }
   }

   return true ;
}

//----------------------------------------------------------------------------
bool
VolZmap::CreateFromTriMesh( const ISurfTriMesh& Surf, double dStep, bool bTriDex, double dExtraBox, int* nError, bool bForceClosed)
{
  // Se la superficie non è chiusa oppure orientata al contrario non ha senso continuare
   double dVol ;
   if ( ( ! Surf.IsClosed()  ||  ! Surf.GetVolume( dVol)  ||  dVol < 0)  &&  ! bForceClosed)
      return false ;

   if ( bForceClosed  && ! bTriDex)
      return false ;

  // Assegno la dimensione della mappa 1 o 3
   m_nMapNum = ( bTriDex ? 3 : 1) ;
 
  // Determino il bounding box della TriMesh
   BBox3d SurfBBox ;
   Surf.GetLocalBBox( SurfBBox) ;

  // Il dexel se parte da un triangolo della trimesh può non trovare l'intersezione,
  // quindi espandiamo il bounding box per ovviare al problema.
   if ( dExtraBox > EPS_ZERO)
      SurfBBox.Expand( dExtraBox) ;
   else
      dExtraBox = 0 ;
    
  // Determino i punti estremi del bounding box
   Point3d ptMapOrig, ptMapEnd ; 
   SurfBBox.GetMinMax( ptMapOrig, ptMapEnd) ;

  // Sistema di riferimento intrinseco dello Zmap
   m_MapFrame.Set( ptMapOrig, Frame3d::TOP) ;

  // Il passo di discretizzazione non può essere inferiore a 100 * EPS_SMALL
   m_dStep = max( dStep, 100 * EPS_SMALL) ;

  // Determino le dimensioni lineari del BBox
   Vector3d vtLen = ptMapEnd - ptMapOrig ;
 
  // A partire dalle dimensioni di xy del grezzo determino il numero di colonne e righe
  // della griglia Zmap e da questi la dimensione del vettore di dexel
   m_nNx[0] = int( ( vtLen.x + EPS_SMALL) / m_dStep + 0.5) ;
   m_nNy[0] = int( ( vtLen.y + EPS_SMALL) / m_dStep + 0.5) ;
   m_nDim[0] = m_nNx[0] * m_nNy[0] ;
  // Ridimensiono il vettore di dexel e creo lo Zmap
   m_Values[0].resize( m_nDim[0]) ;
  
  // Numero di componenti connesse da calcolare
   m_nConnectedCompoCount = - 1 ;

  // Se Tridexel ridimensiono anche gli altri vettori
   if ( bTriDex) {
      m_nNx[1] = m_nNy[0] ;
      m_nNy[1] = int( ( vtLen.z + EPS_SMALL) / m_dStep + 0.5) ;
      m_nDim[1] = m_nNx[1] * m_nNy[1] ;
      m_nNx[2] = m_nNy[1] ;
      m_nNy[2] = m_nNx[0] ;
      m_nDim[2] = m_nNx[2] * m_nNy[2] ;
      #if !defined(_WIN64)
         for ( int i = 0 ; i < ssize( m_nDim) ; ++ i)
            m_nDexelNbr += m_nDim[i] ;
         if ( m_nDexelNbr >= MAX_DEXEL_32_BIT) {
            Clear() ;
            if ( nError != nullptr)
               *nError = 1 ;
            return false ;
         }
      #endif
      m_Values[1].resize( m_nDim[1]) ;
      m_Values[2].resize( m_nDim[2]) ;
   }
  // Se a dimensione singola
   else {
      #if !defined(_WIN64)
         m_nDexelNbr += m_nDim[0] ;
         if ( m_nDexelNbr >= MAX_DEXEL_32_BIT) {
            Clear() ;
            if ( nError != nullptr)
               *nError = 1 ;
            return false ;
         }
      #endif
      m_nNx[1] = 0 ;
      m_nNy[1] = 0 ;
      m_nDim[1] = 0 ;
      m_nNx[2] = 0 ;
      m_nNy[2] = 0 ;
      m_nDim[2] = 0 ;
   }

   // Definisco il numero di blocchi lungo x,y e z
   if ( ! CalcBlockNum())
      return false ;

  // ciclo sulle griglie
   bool bCompleted = true ;
   for ( int nG = 0 ; nG < m_nMapNum ; ++ nG) {

     // Definisco dei sistemi di riferimento ausiliari
      Frame3d frMapFrame ;
      if ( nG == 0) 
         frMapFrame = m_MapFrame ;
      else if ( nG == 1)
         frMapFrame.Set( ptMapOrig, Y_AX, Z_AX, X_AX) ; 
      else if ( nG == 2)
         frMapFrame.Set( ptMapOrig, Z_AX, X_AX, Y_AX) ;

     // Oggetto per calcolo massivo intersezioni
      IntersParLinesSurfTm intPLSTM( frMapFrame, Surf) ;

     // Standarda è multithread
      constexpr bool MULTITHREAD = false ;
      if ( MULTITHREAD) {

         // Numero massimo di thread
         int nThreadMax = max( 1, int( thread::hardware_concurrency()) - 1) ;
         vector< future<bool>> vRes ;
         vRes.resize( nThreadMax) ;
         if ( m_nNx[nG] > m_nNy[nG]) {
            int nDexNum = m_nNx[nG] / nThreadMax ;
            int nRemainder = m_nNx[nG] % nThreadMax ;
            int nInfI = 0 ;
            int nSupI = 0 ;
            for ( int nThread = 0 ; nThread < nThreadMax ; ++ nThread) {
               nInfI = nSupI ;
               nSupI = nInfI + ( nThread < nRemainder ? nDexNum + 1 : nDexNum) ;
               vRes[nThread] = async( launch::async, &VolZmap::CreateMapPart, this, nG,
                                      nInfI, nSupI, 0, m_nNy[nG], ref( vtLen), ref( ptMapOrig), ref( Surf), ref( intPLSTM), bForceClosed) ;
            }
         }
         else {
            int nDexNum = m_nNy[nG] / nThreadMax ;
            int nRemainder = m_nNy[nG] % nThreadMax ;
            int nInfJ = 0 ;
            int nSupJ = 0 ;
            for ( int nThread = 0 ; nThread < nThreadMax ; ++ nThread) {
               nInfJ = nSupJ ;
               nSupJ = nInfJ + ( nThread < nRemainder ? nDexNum + 1 : nDexNum) ;
               vRes[nThread] = async( launch::async, &VolZmap::CreateMapPart, this, nG,
                                      0, m_nNx[nG], nInfJ, nSupJ, ref( vtLen), ref( ptMapOrig), ref( Surf),ref( intPLSTM), bForceClosed) ;
            }
         }
       
         // Ciclo per attendere che tutti gli async abbiano terminato.
         int nTerminated = 0 ;
         while ( nTerminated < nThreadMax) {
            for ( int nL = 0 ; nL < nThreadMax ; ++ nL) {
               // Async terminato
               if ( vRes[nL].valid()  &&  vRes[nL].wait_for( chrono::microseconds{ 1}) == future_status::ready) {
                  ++ nTerminated ;
                  bCompleted = bCompleted  &&  vRes[nL].get() ;
               }
            }
         }
      }

     // !!!! NON MULTITHREAD : SOLO PER DEBUG !!!!
      else {
         CreateMapPart( nG, 0, m_nNx[nG], 0, m_nNy[nG], vtLen, ptMapOrig, Surf, intPLSTM, bForceClosed) ;
      }
   }

  // Assegno il minimo e massimo valore di Z della mappa
   m_dMinZ[0] = dExtraBox ;
   m_dMaxZ[0] = vtLen.z - dExtraBox ;
   m_dMinZ[1] = ( bTriDex ? dExtraBox : 0) ;
   m_dMaxZ[1] = ( bTriDex ? vtLen.x - dExtraBox : 0) ;
   m_dMinZ[2] = ( bTriDex ? dExtraBox : 0) ;  
   m_dMaxZ[2] = ( bTriDex ? vtLen.y - dExtraBox : 0) ;

  // Tipologia
  // Con espansione non va considerato box (calcolo trimesh va in crash)
   m_nShape = ( dExtraBox <= EPS_ZERO && IsBox() ? BOX : GENERIC) ;

  // Aggiornamento dello stato
   m_nStatus = OK ;

   return bCompleted ;
}

vector<Neigh>
VolZmap::GetNeighbours( int nG, int nI, int nJ, double dPar)
{
   vector<Neigh> vNeigh ;
   // parto da in basso a sinistra e procedo per righe, saltando il punto corrente
   if ( nI > 0 && nJ > 0)
      vNeigh.emplace_back( nI - dPar, nJ - dPar, 0) ;
   if ( nJ > 0) {
      vNeigh.emplace_back( nI, nJ - dPar, 1) ;
      if ( nI < m_nNx[nG] - dPar)
         vNeigh.emplace_back( nI + dPar, nJ - dPar, 2) ;
   }
   // passo alla fila contenente il punto corrente
   if ( nI > 0)
      vNeigh.emplace_back( nI - dPar, nJ, 3) ;
   if ( nI < m_nNx[nG] - dPar)
      vNeigh.emplace_back( nI + dPar, nJ, 5) ;
   // passo alla fila successiva a quella del corrente
   if ( nJ < m_nNy[nG] - dPar) {
      if ( nI > 0)
         vNeigh.emplace_back( nI - dPar, nJ + dPar, 6) ;
      vNeigh.emplace_back( nI, nJ + dPar, 7) ;
      if ( nI < m_nNx[nG] - dPar)
         vNeigh.emplace_back( nI + dPar, nJ + dPar, 8) ;
   }
   return vNeigh ;
}

//----------------------------------------------------------------------------
bool
VolZmap::AdjustDexelThroughCracks( int nG, const Vector3d& vtLen, const INTTOCHECKVEC& vIntToCheck, const Point3d& ptMapOrig,
                                   const ISurfTriMesh& Surf, const IntersParLinesSurfTm& intPLSTM)
{
   double dAngSameFace = 45 ;
   double dCosSmall = sin( EPS_ANG_SMALL * DEGTORAD) ;
   // per ogni intersezione segnalata devo cercare sugli spilloni vicini l'intervallo con l'estremo più vicino a questa intersezione
   // e usare l'altro estremo, insieme a quello degli altri spilloni vicini, per estrapolare quello mancante sullo spillone corrente
   Vector3d vtAx ;
   if ( nG == 0)
      vtAx = Z_AX ;
   else if ( nG == 1)
      vtAx = X_AX ;
   else if ( nG == 2)
      vtAx = Y_AX ;
   for ( const auto& itc : vIntToCheck) {
      int nPos = itc.nJ * m_nNx[nG] + itc.nI ;
      double dPosRef = itc.ptInters.v[( nG + 2) % 3] ;
      bool bUseNeighInfo = true ;
      // se il punto di intersezione passato non è valido vuol dire che non ho neanche un'intersezione spaiata da cui partire a cercare
      bool bHaveRef = true ;
      if ( ! isfinite( dPosRef)) {
         bHaveRef = false ;
         bUseNeighInfo = false ;
      }
      PNTVECTOR vNeighInters ;
      Vector3d vtMeanN = V_NULL ;
      if ( bHaveRef) {
         // controllo se esiste già un intervallo con un estremo quasi coincidente a quello segnalato
         // ed entrambi sono entrate o uscite
         bool bNeedToReconstruct = true ;
         for ( const auto& data: m_Values[nG][nPos]) {
            if ( (( abs( dPosRef - data.dMin) < 10 * EPS_SMALL  &&  ( itc.vtN * vtAx) * ( data.vtMinN * vtAx) > 0))  ||
                  (( abs( dPosRef - data.dMax) < 10 * EPS_SMALL  &&  ( itc.vtN * vtAx) * ( data.vtMaxN * vtAx) > 0))) {
               bNeedToReconstruct = false ;
               break ;
            }  
         }
         if ( ! bNeedToReconstruct)
            continue ;
         vector<Neigh> vNeigh = GetNeighbours( nG, itc.nI, itc.nJ, m_dStep) ;
         bool bSearchingInOrOut = itc.vtN * vtAx > 0 ;
         for ( const auto& [i, j, n] : vNeigh) {
            double dMinDist = INFINITO ;
            double dCorrespInters = INFINITO ;
            double dNearestPos = NAN ;
            Vector3d vtNearestN ;
            int nNeighPos = int( j) * m_nNx[nG] + int( i) ;
            for ( const auto& data : m_Values[nG][nNeighPos]) {
               double dDist1 = abs( dPosRef - data.dMin) ;
               double dDist2 = abs( dPosRef - data.dMax) ;
               if ( dDist1 < dMinDist  ||  dDist2 < dMinDist) {
                  dMinDist = min( dDist1, dDist2) ;

                  // se sto cercando un ingresso
                  if ( bSearchingInOrOut) {
                     if ( data.vtMaxN * vtAx < 0) {
                        dNearestPos = data.dMax ;
                        vtNearestN = data.vtMaxN ;
                        dCorrespInters = data.dMin ;
                     }
                     else if ( data.vtMinN * vtAx < 0) {
                        dNearestPos = data.dMin ;
                        vtNearestN = data.vtMinN ;
                        dCorrespInters = data.dMax ;
                     }
                  }
                  // se sto cercando un'uscita
                  else {
                     if ( data.vtMaxN * vtAx > 0) {
                        dNearestPos = data.dMax ;
                        vtNearestN = data.vtMaxN ;
                        dCorrespInters = data.dMin ;
                     }
                     else if ( data.vtMinN * vtAx > 0) {
                        dNearestPos = data.dMin ;
                        vtNearestN = data.vtMinN ;
                        dCorrespInters = data.dMax ;
                     }
                  }
               }
            }
            // verifico che il dNearestPos sia effettivamente associato al tratto che sto cercando di ricostruire
            // e non ad un altro presente sullo spillone corrente
            bool bAddNearest = true ;
            for ( const auto& data : m_Values[nG][nPos]) {
               double dDist1 = abs( dCorrespInters - data.dMin) ;
               double dDist2 = abs( dCorrespInters - data.dMax) ;
               if ( dDist1 < dMinDist  ||  dDist2 < dMinDist) {
                  bAddNearest = false ;
                  break ;
               }
            }
            if ( bAddNearest  &&  isfinite( dNearestPos)) {
               // controllo anche che i punti non siano troppo lontani tra loro, altrimenti
               // il piano calcolato potrebbe essere una media troppo grezza
               double dMaxDist = 3 * m_dStep ;
               const bool bTooFar = any_of( vNeighInters.begin(), vNeighInters.end(),
                                            [dNearestPos, dMaxDist](const Point3d& pt) {
                                               return abs( pt.z - dNearestPos) > dMaxDist ;
                                            }) ;
               double dAng ;
               if ( ! bTooFar) {
                  // controllo anche che gli angoli delle normali non cambino troppo
                  // calcolo l'angolo con la media delle precedenti perché il valore oltre cui
                  // considero le normali diverse è grande
                  vtNearestN.GetAngle( vtMeanN, dAng) ;
               }
               bAddNearest = ! bTooFar  &&  dAng < dAngSameFace ;
               if ( bAddNearest) {
                  double dX = ( i + 0.5) * m_dStep ;
                  double dY = ( j + 0.5) * m_dStep ;
                  vNeighInters.emplace_back( dX, dY, dNearestPos) ;
                  vtMeanN += vtNearestN ;
               }
               else {
                  bUseNeighInfo = false ;
                  break ;
               }
            }
         }
      }
      vector<Data> vNewIntervals ;
      bool bAddAllIntervals = false ;
      Point3d ptInt ;
      // se riesco uso i vicini per calcolare l'intersezione mancante
      if ( bUseNeighInfo) {
         // se non ci sono almeno tre punti, non allineati, non riesco a determinare il piano con cui tagliare lo spillone
         if ( ssize( vNeighInters) < 3)
            continue ;
         else if ( ssize( vNeighInters) == 3){
            // verifico che non siano allineati
            if ( ( vNeighInters[0].x == vNeighInters[1].x  &&  vNeighInters[0].x ==  vNeighInters[1].x)  ||
                  ( vNeighInters[0].y == vNeighInters[1].y  &&  vNeighInters[0].y ==  vNeighInters[1].y))
               continue ;
         }
         // se sono più di tre sicuramente non sono allieanti, visto che sono intorno allo spillone corrente
         vtMeanN.Normalize() ;
         Plane3d plLoc ;
         PolyLine PL ; PL.FromPointVector( vNeighInters) ;
         PL.IsFlat( plLoc, INFINITO) ;
            
         // porto anche il punto dell'intersezione nel frame della griglia
         Point3d ptKnown = itc.ptInters ;
         Frame3d frGrid ;
         if ( nG == 0)
            frGrid.Set( ORIG, Frame3d::TOP) ;
         else if ( nG == 1)
            frGrid.Set( ORIG, Frame3d::RIGHT) ;
         else if ( nG == 2)
            frGrid.Set( ORIG, Z_AX, X_AX, Y_AX) ;
         ptKnown.ToLoc( frGrid) ;
         IntersLinePlane( ptKnown, Z_AX, INFINITO, plLoc, ptInt, false) ;
      }
      // altrimenti sposto di poco lo spillone e ricalcolo le intersezioni con la superficie, finché trovo l'intersezione mancante
      else {
         bool bFound = false ;
         DBLVECTOR vSmallStep ;
         vSmallStep.push_back( Clamp( 0.05 * m_dStep, 50 * EPS_SMALL, 1.)) ;
         vSmallStep.push_back( Clamp( 0.2 * m_dStep, 50 * EPS_SMALL, 1.)) ;
         // ciclo una raggiera di spilloni ( a distanza inferiore allo m_dStep) attorno allo spillone
         // appena trova un'intersezione buona mi fermo
         // se ciclando la raggiera non trovo un'intersezione valida allora potrei essermi spostato troppo poco dallo spillone originale e essere caduto ancora nella fessura
         // quindi ciclo nuovamente aumentando un poco la distanza dallo spillone originale ( comunque sempre meno della distanza m_dStep)
         for ( int c = 0 ; c < ssize( vSmallStep) ; ++c) {
            vector<Neigh> vSmallMoveNeigh = GetNeighbours( nG, itc.nI, itc.nJ, vSmallStep[c]) ;
            for ( int n = 0 ; n < ssize( vSmallMoveNeigh) ; ++n) {
               double dX = ( vSmallMoveNeigh[n].nI + 0.5) * m_dStep ;
               double dY = ( vSmallMoveNeigh[n].nJ + 0.5) * m_dStep ;
               Point3d ptP0( dX, dY, 0) ;
               ILSIVECTOR IntersectionResults ;
               intPLSTM.GetInters( ptP0, vtLen.v[(nG+2)%3], IntersectionResults) ;
               // intervalli dalle intersezioni
               vNewIntervals.clear() ;
               DBLVECTOR vDiscardedInters ;
               bool bInside = false ;
               Vector3d vtInN ;
               Point3d ptIn ;
               for ( int i = 0 ; i < ssize( IntersectionResults) ; ++i) {
                  // Se c'è intersezione
                  if ( IntersectionResults[i].nILTT != ILTT_NO) {
                     if ( i > 0) {
                        int z = i - 1 ;
                        // controllo se l'intersezione corrente è coincidente con la precedente, per esempio se ILTT == 4 ( intersezione su un lato di un triangolo, quindi con 2 triangoli)
                        if ( IntersectionResults[i].dU - IntersectionResults[z].dU < EPS_SMALL  &&
                           IntersectionResults[i].dCosDN - IntersectionResults[z].dCosDN < EPS_SMALL)
                           continue ;
                     }
                     double dCos = IntersectionResults[i].dCosDN ;
                     // entro nella superficie trimesh
                     if ( dCos < - dCosSmall) {
                        if ( bInside) {
                           int k = i - 1 ;
                           vDiscardedInters.push_back( IntersectionResults[k].dU) ;
                        }
                        ptIn = IntersectionResults[i].ptI ;

                        int nT = IntersectionResults[i].nT ;
                        int nF = Surf.GetFacetFromTria( nT) ;

                        Surf.GetFacetNormal( nF, vtInN) ;

                        bInside = true ;
                     }

                     // esco dalla superficie trimesh
                     else if ( dCos > dCosSmall) {
                        if ( ! bInside) {
                           vDiscardedInters.push_back( IntersectionResults[i].dU) ;
                           continue ;
                        }
                        Point3d ptOut = IntersectionResults[i].ptI ;

                        int nT = IntersectionResults[i].nT ;
                        int nF = Surf.GetFacetFromTria( nT) ;

                        Vector3d vtOutN ;
                        Surf.GetFacetNormal( nF, vtOutN) ; 

                        // Aggiungo un tratto al dexel
                        vNewIntervals.emplace_back() ;

                        // Aggiorno dati del tratto di dexel 
                        vNewIntervals.back().dMin = ptIn.v[(nG+2)%3] - ptMapOrig.v[(nG+2)%3] ;
                        vNewIntervals.back().dMax = ptOut.v[(nG+2)%3] - ptMapOrig.v[(nG+2)%3] ;
                        vNewIntervals.back().vtMinN = vtInN ;
                        vNewIntervals.back().vtMaxN = vtOutN ;
                        vNewIntervals.back().nToolMin = 0 ;
                        vNewIntervals.back().nToolMax = 0 ;
                        vNewIntervals.back().nCompo = 0 ;

                        bInside = false ;
                     }
                  } 
               }
               // controllo che abbia almeno un intervallo in più rispetto a quelli di prima
               if ( ssize( vNewIntervals)  >  ssize( m_Values[nG][nPos])) {
                  // se avevo un punto disaccoppiato lo uso come riferimento
                  if ( itc.ptInters.IsValid()) {
                     // devo controllare se l'intervallo in più è quello che cercavo
                     // controllo anche tra le intersezioni scartate
                     double dMinDist = INFINITO ;
                     for ( const auto& data : vNewIntervals) {
                        if ( abs( data.dMin - dPosRef) < dMinDist) {
                           dMinDist = abs( data.dMin - dPosRef) ;
                           ptInt.Set( itc.nI + 0.5, itc.nJ + 0.5, data.dMax) ;
                           vtMeanN = data.vtMaxN ;
                           bFound = true ;
                        }
                        if ( abs( data.dMax - dPosRef) < dMinDist) {
                           dMinDist = abs( data.dMax - dPosRef) ;
                           ptInt.Set( itc.nI + 0.5, itc.nJ + 0.5, data.dMin) ;
                           vtMeanN = data.vtMinN ;
                           bFound = true ;
                        }
                     }
                     for ( const auto& dU : vDiscardedInters) {
                        if ( abs( dU - dPosRef) < dMinDist) {
                           bFound = false ;
                           break ;
                        }
                     }
                  }
                  else {
                     // se non avevo intervalli sul corrente tengo tutto quello trovo
                     if ( ssize( m_Values[nG][nPos]) == 0) {
                        bFound = true ;
                        bAddAllIntervals = true ;
                     }
                     else {
                        // altrimenti mi aspetto di aver attrvarsato due fessure e quindi di avere sullo spillone corrente un
                        // intervallo il cui inizio e fine corrispondono a inizio di un intervallo e fine del successivo
                        for ( int z = 0 ; z < ssize( vNewIntervals) - 1 ; ++z) {
                           Data& data = vNewIntervals[z] ;
                           Data& dataNext = vNewIntervals[z+1] ;
                           bool bSkipNext = false ;
                           for ( int w = 0 ; w < ssize( m_Values[nG][nPos]) ; ++w) {
                              Data& dataCurr = m_Values[nG][nPos][w] ;
                              if ( abs( data.dMin - dataCurr.dMin) < m_dStep  &&  abs( dataNext.dMax - dataCurr.dMax) < m_dStep) {
                                 m_Values[nG][nPos].emplace_back() ;
                                 m_Values[nG][nPos].back().dMin = dataCurr.dMin;
                                 m_Values[nG][nPos].back().dMax = data.dMax ;
                                 m_Values[nG][nPos].back().vtMinN = dataCurr.vtMinN ;
                                 m_Values[nG][nPos].back().vtMaxN = data.vtMaxN ;
                                 m_Values[nG][nPos].back().nToolMin = 0 ;
                                 m_Values[nG][nPos].back().nToolMax = 0 ;
                                 m_Values[nG][nPos].back().nCompo = 0 ;
                                 //secondo
                                 m_Values[nG][nPos].emplace_back() ;
                                 m_Values[nG][nPos].back().dMin = dataNext.dMin;
                                 m_Values[nG][nPos].back().dMax = dataCurr.dMax ;
                                 m_Values[nG][nPos].back().vtMinN = dataNext.vtMinN ;
                                 m_Values[nG][nPos].back().vtMaxN = dataCurr.vtMaxN ;
                                 m_Values[nG][nPos].back().nToolMin = 0 ;
                                 m_Values[nG][nPos].back().nToolMax = 0 ;
                                 m_Values[nG][nPos].back().nCompo = 0 ;
                                 // cancello quello che era
                                 m_Values[nG][nPos].erase( m_Values[nG][nPos].begin() + w) ;
                                 bSkipNext = true ;
                                 break ;
                              }
                           }
                           if ( bSkipNext)
                              ++z ;
                        }
                     }
                  }
               }
               if ( bFound)
                  break ;
            }
            if ( bFound)
               break ;
         }
         // se cercando nei d'intorni dello spillone non ho trovato l'intersezione mancante, mi arrendo
         if ( ! bFound)
            continue ;
      }
      

      if ( ! bAddAllIntervals) {
         double dMin, dMax ;
         Vector3d vtMinN, vtMaxN ;
         if ( itc.dU < ptInt.z) {
            dMin = itc.dU ;
            dMax = ptInt.z ;
            vtMinN = itc.vtN ;
            vtMaxN = vtMeanN ;
         }
         else {
            dMin = ptInt.z ;
            dMax = itc.dU ;
            vtMinN = vtMeanN ;
            vtMaxN = itc.vtN ;
         }
         //// prima di aggiungere il tratto verifico che non si sovrapponga con un tratto già presente
         //for ( const auto& data: m_Values[nG][nPos]) {
         //   if ( ( dMin > data.dMin  &&  dMin < data.dMax)  ||
         //        ( dMax > data.dMin  &&  dMax < data.dMax)  ||
         //        ( dMin < data.dMin  &&  dMax > data.dMax))
         //      continue ;
         //}

         m_Values[nG][nPos].emplace_back() ;
         m_Values[nG][nPos].back().dMin = dMin;
         m_Values[nG][nPos].back().dMax = dMax ;
         m_Values[nG][nPos].back().vtMinN = vtMinN ;
         m_Values[nG][nPos].back().vtMaxN = vtMaxN ;
         m_Values[nG][nPos].back().nToolMin = 0 ;
         m_Values[nG][nPos].back().nToolMax = 0 ;
         m_Values[nG][nPos].back().nCompo = 0 ;
      }
      else {
         for ( const auto& data: vNewIntervals) {
            m_Values[nG][nPos].emplace_back() ;
            m_Values[nG][nPos].back().dMin = data.dMin;
            m_Values[nG][nPos].back().dMax = data.dMax ;
            m_Values[nG][nPos].back().vtMinN = data.vtMinN ;
            m_Values[nG][nPos].back().vtMaxN = data.vtMaxN ;
            m_Values[nG][nPos].back().nToolMin = 0 ;
            m_Values[nG][nPos].back().nToolMax = 0 ;
            m_Values[nG][nPos].back().nCompo = 0 ;
         }
      }
      std::sort( m_Values[nG][nPos].begin(), m_Values[nG][nPos].end(), []( Data& a, Data& b) { return a.dMin < b.dMin ;}) ;
   }

   return true ;
}

//----------------------------------------------------------------------------
bool
VolZmap::SearchMissingIntervals( int nMap, INTTOCHECKVEC& vIntToCheck)
{
   
   for ( int nI = 0 ; nI < m_nNx[nMap] ; ++nI) {
      for ( int nJ = 0 ; nJ < m_nNy[nMap] ; ++nJ) {
         int nPos = int( nJ) * m_nNx[nMap] + int( nI) ;
         vector<Neigh> vNeigh = GetNeighbours( nMap, nI, nJ, 1) ;
         for ( const auto& [i,j,n] : vNeigh) {
            int nNeighPos = int( j) * m_nNx[nMap] + int( i) ;
            if ( ssize( m_Values[nMap][nNeighPos]) > ssize( m_Values[nMap][nPos])) {
               int nNeigh = 8 - n ;
               auto itOpp = find_if( vNeigh.begin(), vNeigh.end(), [nNeigh]( const Neigh& a) { return a.nN == nNeigh ;}) ;
               int nIndOpp = distance( begin( vNeigh), itOpp) ;
               if ( itOpp != vNeigh.end()) {
                  int nOppNeighPos = int( vNeigh[nIndOpp].nJ) * m_nNx[nMap] + int( vNeigh[nIndOpp].nI) ;
                  if ( ssize( m_Values[nMap][nOppNeighPos]) == ssize( m_Values[nMap][nNeighPos])) {
                     vIntToCheck.emplace_back( P_INVALID, V_INVALID, double( NAN), nMap, nI, nJ) ;
                     break ;
                  }
               }
            }
         }
      }
   }
   return true ;
}