//----------------------------------------------------------------------------
//  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/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)
{
   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) ;

   // 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()) ;

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

         bool bInside = false ;
         Point3d ptIn ;
         Vector3d vtInN ;

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

            int nIntType = IntersectionResults[k].nILTT ;

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

               double dCos = IntersectionResults[k].dCosDN ;

               // entro nella superficie trimesh
               if ( dCos < - dCosSmall) {

                  ptIn = IntersectionResults[k].ptI ;

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

                  Surf.GetFacetNormal( nF, vtInN) ;

                  bInside = true ;
               }

               // esco dalla superficie trimesh
               else if ( dCos > dCosSmall  &&  bInside) {

                  Point3d ptOut = IntersectionResults[k].ptI ;

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

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

                  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 ;
               }
            }   
         }
      }
   }
   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)
{
  // Se la superficie non è chiusa oppure orientata al contrario non ha senso continuare
   double dVol ;
   if ( ! Surf.IsClosed()  ||  ! Surf.GetVolume( dVol)  ||  dVol < 0)
      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 = true ;
      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)) ;
            }
         }
         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)) ;
            }
         }
       
         // 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) ;
      }
   }

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