//----------------------------------------------------------------------------
//  EgalTech 2024-2024
//----------------------------------------------------------------------------
// File : CAvSilhouetteSurfTm.cpp    Data : 08.06.24    Versione : 2.6f2
// Contenuto : Implementazione della funzione CAvSilhouetteSurfTm.
//
//
//
// Modifiche : 08.06.24 DS Creazione modulo.
//
//
//----------------------------------------------------------------------------

#include "stdafx.h"
#include "CAvSilhouetteSurfTm.h"
#include "GeoConst.h"
#include "/EgtDev/Include/EGkChainCurves.h"
#include "/EgtDev/Include/EGkOffsetCurve.h"
#include "/EgtDev/Include/EGkCurveLine.h"
#include "/EgtDev/Include/EGkCurveArc.h"
#include "/EgtDev/Include/EGkIntersCurves.h"
#include "/EgtDev/Include/EGkIntersLineSurfTm.h"
#include "/EgtDev/Include/EGkSfrCreate.h"
#include "/EgtDev/Include/EGkDistPointLine.h"
#include <thread>
#include <future>

using namespace std ;

// Test di Debug per disegnare le posizioni dei punti per la ricostruzione degli spigoli vivi
#define ENABLE_SHARPED_EDGES_DEBUG 0
#if ENABLE_SHARPED_EDGES_DEBUG
   string Sharped_Edges_Debug_File_Name = "C:\\Temp\\SharpedEdges.nge" ;
#endif
// Test di Debug per disegnare le quote dove sono trovate le intersesioni tra il tool e le trimesh
#define ENABLE_COLORED_GRID_DEBUG 0
#if ENABLE_COLORED_GRID_DEBUG
   string Colored_Grid_Debug_File_Name = "C:\\Temp\\ColoredGrid.nge" ;
#endif
// Test di Debug per disegnare i segmenti ricavati dal marchingSquares
#define ENABLE_PROCESS_SQUARE_DEBUG 0
#if ENABLE_PROCESS_SQUARE_DEBUG
   string Process_Square_Debug_File_Name = "C:\\Temp\\ProcessSquares.nge" ;
#endif
// Test di Debug per disegnare i punti ricavati con il metodo di Bisezione
#define ENABLE_BISECTION_DEBUG 0
#if ENABLE_BISECTION_DEBUG
   string Bisection_Debug_File_Name = "C:\\Temp\\Bisection.nge" ;
#endif
// Inclusione file per salvataggio
#if ENABLE_SHARPED_EDGES_DEBUG  ||  ENABLE_COLORED_GRID_DEBUG  ||  ENABLE_PROCESS_SQUARE_DEBUG  ||  ENABLE_BISECTION_DEBUG
   #include "/EgtDev/Include/EGkGeoPoint3d.h"
   #include "/EgtDev/Include/EGkGeoObjSave.h"
#endif

//----------------------------------------------------------------------------
// Funzioni locali per calcolo isolinee con metodo Marching Squares
//----------------------------------------------------------------------------

//----------------------------------------------------------------------------
// Quadrato (C=corner E=edge) :
// 
//      C3 - E2 - C2
//      |          |
//      E3        E1
//      |          |
//      C0 - E0 - C1
// 
//----------------------------------------------------------------------------

// tipo di punto della PolyLinea
enum PointType {
   VALID = 0 ,             // punto valido da conservare
   INVALID = -1,           // punto da scartare
   SHARPED_EXT = 1,        // punto per spigolo vivo (con angolo esterno)
   SHARPED_INT = 2,        // punto per spigolo vivo (con angolo interno)
} ;

//----------------------------------------------------------------------------
static Point3d
CalcPoint( const Point3d& ptS, const Point3d& ptE, double dLevel, double dOffsR,
           const ICAvToolSurfTm& cavTstm, const Frame3d &frGrid)
{
  // Ricerca con metodo di bisezione (si arresta quando il medio è allineato agli estremi)
   const int MAX_ITER = 10 ;
   int nCount = 0 ;
   Point3d ptP1 = ptS ;
   Point3d ptP2 = ptE ;
   Point3d ptMid = Media( ptP1, ptP2) ;
   while ( nCount < MAX_ITER) {
     // verifica del punto medio
      ptMid.z = 0 ;
      Point3d ptTest = GetToGlob( ptMid + ( cavTstm.GetToolHeight() - dOffsR) * Z_AX, frGrid) ;
      double dMove ;
      cavTstm.TestPosition( ptTest, frGrid.VersZ(), frGrid.VersZ(), dMove) ;
      ptMid.z = dMove ;
     // se estremi coincidenti, allora mi fermo
      if ( SqDistXY( ptP1, ptP2) < 4 * SQ_EPS_SMALL) {
         ptMid.z = dLevel ;
         return ptMid ;
      }
     // altrimenti nuova suddivisione della parte con estremi da parte opposta rispetto al livello
      bool b1On = ( ptP1.z > dLevel) ;
      bool bMidOn = ( ptMid.z > dLevel) ;
      if ( b1On != bMidOn)
         ptP2 = ptMid ;
      else
         ptP1 = ptMid ;
      ptMid = Media( ptP1, ptP2) ;
      ++ nCount ;
   }
   ptMid.z = dLevel ;
   return ptMid ;
}

//----------------------------------------------------------------------------
static int
ProcessSquare( int nFlag, double dLevel, double dQPt0, double dQpt1, double dQpt2, double dQpt3,
               int& nI1s, int& nI1e, int& nI2s, int& nI2e)
{
   static int LineTable[16][5] = { { 0, -1, -1, -1, -1},    // ( 0)
                                   { 1,  0,  3, -1, -1},    // ( 1)
                                   { 1,  1,  0, -1, -1},    // ( 2)
                                   { 1,  1,  3, -1, -1},    // ( 3)
                                   { 1,  2,  1, -1, -1},    // ( 4)
                                   { 2,  0,  1,  2,  3},    // ( 5.0) -> primo caso ambiguo
                                   { 1,  2,  0, -1, -1},    // ( 6)
                                   { 1,  2,  3, -1, -1},    // ( 7)
                                   { 1,  3,  2, -1, -1},    // ( 8)
                                   { 1,  0,  2, -1, -1},    // ( 9)
                                   { 2,  1,  2,  3,  0},    // ( 10.0) -> primo caso ambiguo
                                   { 1,  1,  2, -1, -1},    // ( 11)
                                   { 1,  3,  1, -1, -1},    // ( 12)
                                   { 1,  0,  1, -1, -1},    // ( 13)
                                   { 1,  3,  0, -1, -1},    // ( 14)
                                   { 0, -1, -1, -1, -1}} ;  // ( 15)

   static int LineTableAmbiguos[2][5] = { { 2, 0, 3, 2, 1},       // ( 5.1) -> secondo caso ambiguo
                                          { 2, 1, 0, 3, 2}} ;     // ( 10.1) -> secondo caso ambiguo

  // flag fuori dai limiti
   if ( nFlag < 0  ||  nFlag > 15)
      return -1 ;
  // nessuna linea
   if ( LineTable[nFlag][0] == 0)
      return 0 ;
  // una linea
   else if ( LineTable[nFlag][0] == 1) {
      nI1s = LineTable[nFlag][1] ;
      nI1e = LineTable[nFlag][2] ;
      return 1 ;
   }
  // due linee
   else if ( LineTable[nFlag][0] == 2) {
     // controllo in quale caso sono
      if ( nFlag == 5) {
        // caso 5.0
         if ( ( dQPt0 - dLevel) * ( dQpt2 - dLevel) > ( dQpt1 - dLevel) * ( dQpt3 - dLevel)) {
            nI1s = LineTable[nFlag][1] ;
            nI1e = LineTable[nFlag][2] ;
            nI2s = LineTable[nFlag][3] ;
            nI2e = LineTable[nFlag][4] ;
         }
        // caso 5.1
         else {
            nI1s = LineTableAmbiguos[0][1] ;
            nI1e = LineTableAmbiguos[0][2] ;
            nI2s = LineTableAmbiguos[0][3] ;
            nI2e = LineTableAmbiguos[0][4] ;
         }
      }
      else if ( nFlag == 10) {
        // caso 10.0
         if ( ( dQPt0 - dLevel) * ( dQpt2 - dLevel) < ( dQpt1 - dLevel) * ( dQpt3 - dLevel)) {
            nI1s = LineTable[nFlag][1] ;
            nI1e = LineTable[nFlag][2] ;
            nI2s = LineTable[nFlag][3] ;
            nI2e = LineTable[nFlag][4] ;
         }
        // caso 10.1
         else {
            nI1s = LineTableAmbiguos[1][1] ;
            nI1e = LineTableAmbiguos[1][2] ;
            nI2s = LineTableAmbiguos[1][3] ;
            nI2e = LineTableAmbiguos[1][4] ;
         }
      }
      return 2 ;
   }

   return -1 ;
}

//----------------------------------------------------------------------------
static bool
TestSubEdges( unordered_map<int, Point3d>& umEdgePnt, const INTVECTOR& vEdgeInd, int nFirst, int nLast,
              const DBLVECTOR& vdGrid, int nStepX, double dStep,
              double dLevel, double dOffsR, const ICAvToolSurfTm& cavTstm, const Frame3d &frGrid)
{
   for ( int k = nFirst ; k <= nLast ; ++ k) {
     // recupero i differenti indici
      int nKey = vEdgeInd[k] ;
      int nInd = nKey / 2 ;
      bool bOnX = ( ( nKey % 2) == 0) ;
      int j = nInd / ( nStepX + 1) ;
      int i = nInd % ( nStepX + 1) ;
     // determino i punti estremi dell'edge
      Point3d ptS{ i * dStep, j * dStep, vdGrid[nInd]} ;
      Point3d ptE{ ptS.x + ( bOnX ? dStep : 0), ptS.y + ( bOnX ? 0 : dStep), vdGrid[nInd + ( bOnX ? 1 : nStepX + 1)]} ;
     // calcolo il punto
      Point3d ptQ = CalcPoint( ptS, ptE, dLevel, dOffsR, cavTstm, frGrid) ;
      umEdgePnt[ nKey] = ptQ ;
   }
   return true ;
}

//----------------------------------------------------------------------------
static bool
ModifyPolyLineToSharped( PNTIVECTOR& vAdvPt, const CAvToolSurfTm& cavTstm, const Frame3d& frGrid,
                         double dAngTol, double dStep, double dSegLen)
{
  // dichiaro i versori tangenti consecutivi per il percorso della PolyLinea
   Vector3d vtTanPrev, vtTan, vtTanAfter ;
  // ciclo i punti successivi
   for ( int i = 2 ; i < int( vAdvPt.size()) - 1 ; ++ i) {
     // calcolo il versore tangente del segmento prececente
      vtTanPrev = vAdvPt[i-1].first - vAdvPt[i-2].first ;
      vtTanPrev.Normalize() ;
     // calcolo versore tangente attuale
      vtTan = ( vAdvPt[i].first - vAdvPt[i-1].first) ;
      vtTan.Normalize() ;
     // recupero l'angolo con segno tra la direzione precedente e la direzione attuale
      double dAngDeg ;
      if ( ! vtTanPrev.GetAngleXY( vtTan, dAngDeg))
         continue ;
     // potrei modificare il punto i-esimo nei seguenti casi :
     // 1) l'angolo tra i versori supera la tolleranza
     // 2) la distana tra il punto i-esimo e il primo non valido è compresa in due quadrati adiacenti
      bool bModifyPoint = ( ( abs( dAngDeg) > dAngTol)) ;
      if ( bModifyPoint) {
         for ( int j = i - 1 ; j >= max( 0, i - 2)  &&  bModifyPoint ; -- j) {
            if ( vAdvPt[j].second == PointType::INVALID) {
               bModifyPoint = ( SqDist( vAdvPt[i].first, vAdvPt[j].first) < 8 * dStep * dStep) ;
               break ;
            }
         }
      }
      if ( bModifyPoint) {
        // recupero la direzione successiva
         vtTanAfter = ( vAdvPt[i+1].first - vAdvPt[i].first) ;
         vtTanAfter.Normalize() ;
        // definisco i segmenti di raccordo, lunghi quanto la diagonale principale della griglia
         PtrOwner<ICurveLine> pSegPrev( CreateCurveLine()) ;
         if ( IsNull( pSegPrev)  ||  ! pSegPrev->Set( vAdvPt[i-1].first, vAdvPt[i-1].first + vtTanPrev * dSegLen))
            return false ;
         PtrOwner<ICurveLine> pSegAft( CreateCurveLine()) ;
         if ( IsNull( pSegAft)  ||  ! pSegAft->Set( vAdvPt[i].first, vAdvPt[i].first - vtTanAfter * dSegLen))
            return false ;
        // cerco l'intersezione tra i due segmenti
         IntersCurveCurve ICC( *pSegPrev, *pSegAft) ;
         bModifyPoint = ( ICC.GetIntersCount() > 0) ;
         if ( bModifyPoint) {
            IntCrvCrvInfo aInfo ;
            if ( ! ICC.GetIntCrvCrvInfo( 0, aInfo))
               return false ;
           // recupero il punto di intersezione
            Point3d ptInters = aInfo.IciA->ptI ;
           // controllo la distanza con tale punto di intersezione e i punti precedenti.
           // Questa deve essere contenuta in 5 quadrati della griglia
            bModifyPoint = ( SqDist( ptInters, vAdvPt[i-1].first) < 50 * dStep * dStep  &&
                             SqDist( ptInters, vAdvPt[i].first) < 50 * dStep * dStep) ;
            if ( bModifyPoint) {
              // determino versore di movimento
               Vector3d vtMove = ( vtTanPrev - vtTanAfter) ;
               vtMove.Normalize() ;
              // stabilisco se l'angolo attuale è esterno o interno
               bool bIsExternalAngle = ( dAngDeg > 0.) ;
              // dal punto di intersezione mi allontano leggermente dal materiale (1/16 del passo griglia)
               if ( bIsExternalAngle)
                  vtMove *= ( dStep / 16.) ;
               else
                  vtMove *= - ( dStep / 16.) ;
              // definisco il punto di test
               Point3d ptTest = ptInters + vtMove ;
              // verifico se l'utensile tocca il materiale se posizionato in ptTest
               double dToTDist = 0. ;
               ptTest.z += cavTstm.GetToolHeight() ;
               ptTest.ToGlob( frGrid) ;
               if ( ! cavTstm.TestPosition( ptTest, frGrid.VersZ(), frGrid.VersZ(), dToTDist))
                  return false ;
              // modifico se non ho intersezione con TriMesh, allora aggiungo il punto
               bModifyPoint = ( dToTDist < EPS_SMALL) ;
               if ( bModifyPoint) {
                 // il punto i-esimo diventa il punto di intersezione ( quindi di tipo Sharped)
                  vAdvPt[i].first = ptInters ;
                  vAdvPt[i].second = ( bIsExternalAngle ? PointType::SHARPED_EXT : PointType::SHARPED_INT) ;
                 // il punto (i-1)-esimo essendo in tangenza diventa invalido
                  vAdvPt[i-1].second = PointType::INVALID ;
               }
            }
         }
      }
   }

   return true ;
}

//----------------------------------------------------------------------------
static bool
TestEdgesClosedPolyLines( POLYLINEVECTOR& vPL, const CAvToolSurfTm& cavTstm, int nStepX, int nStepY,
                          const Frame3d& frGrid, double dDimZ, double dLevel, double dStep, double dAngTol,
                          bool bAdvCorners)
{
  // se non ho polylinee, allora esco
   if ( vPL.empty())
      return true ;
   dAngTol = max( dAngTol, ANG_TOL_STD_DEG / 5.) ;
   double dSegLen = 50 * dStep ;

  // cerco la PolyLine di area maggiore, determinerà il verso di percorrenza
   double dMaxArea = 0. ;
   int nIndMaxArea = 0 ;
   for ( int i = 0 ; i < int( vPL.size()) ; ++ i) {
     // semplifico rimuovendo i punti allineati
      vPL[i].RemoveAlignedPoints( 10 * EPS_SMALL) ;
     // calcolo l'area
      double dMyArea ;
      if ( vPL[i].GetAreaXY( dMyArea)  &&  dMyArea > dMaxArea) {
         dMaxArea = dMyArea ;
         nIndMaxArea = i ;
      }
   }
   if ( nIndMaxArea > 0)
      swap( vPL[0], vPL[nIndMaxArea]) ;

  // ordino le PolyLines cercando loop interni ed esterni (gli interni vengono invertiti)
   INTMATRIX mIndMat ;
   BOOLVECTOR vbInv ;
   Vector3d vtN ;
   if ( ! CalcRegionPolyLines( vPL, vtN, mIndMat, vbInv))
      return false ;

   #if ENABLE_SHARPED_EDGES_DEBUG
      vector<vector<IGeoObj*>> vvpGObj ; vvpGObj.resize( 7) ;
      vector<vector<Color>> vvCol ; vvCol.resize( 7) ;
     // ombra del tool
      PtrOwner<ICurveArc> pCrvToolShape( CreateCurveArc()) ;
      pCrvToolShape->Set( ORIG, Z_AX, cavTstm.GetToolRadius()) ;
     // 1° gruppo, griglia di punti
      for ( int j = 0 ; j <= nStepY ; ++ j) {
         for ( int i = 0 ; i <= nStepX ; ++ i) {
            Point3d ptP = Point3d( i * dStep, j * dStep, dDimZ) ;
            PtrOwner<IGeoPoint3d> myPtGrid( CreateGeoPoint3d()) ;
            myPtGrid->Set( ptP) ;
            vvpGObj[0].emplace_back( static_cast<IGeoObj*>( Release( myPtGrid))) ;
            vvCol[0].emplace_back( BLUE) ;
         }
      }
     // 2° gruppo, superfici
      CISURFTMPVECTOR myStmVector = cavTstm.GetvStm() ;
      for ( int i = 0 ; i < int( myStmVector.size()) ; ++ i) {
         PtrOwner<ISurfTriMesh> pMyStm( CloneSurfTriMesh( myStmVector[i])) ;
         if ( ! IsNull( pMyStm)) {
            pMyStm->ToLoc( frGrid) ;
            vvpGObj[1].emplace_back( static_cast<IGeoObj*>( Release( pMyStm))) ;
            vvCol[1].emplace_back( Color( 0., 1., 0., .5)) ;
         }
      }
   #endif

  // vettore delle polyline visitate ( viste però come curve composite)
   ICRVCOMPOPOVECTOR vPL_AsCompo ; vPL_AsCompo.reserve( vPL.size()) ;

  // scorro i Chunk della matrice di interi, quindi le righe
   for ( int nChunk = 0 ; nChunk < int( mIndMat.size()) ; ++ nChunk) {
     // scorro le PolyLinee presenti nel Chunk
      for ( int nLoop = 0 ; nLoop < int( mIndMat[nChunk].size()) ; ++ nLoop) {
        // recupero l'indice della PolyLine di riferimento
         int nPol = mIndMat[nChunk][nLoop] ;
         if ( nPol < 0  ||  nPol >= int( vPL.size()))
            return false ;
        // essendo chiusa, se presenti meno di 4 punti, non faccio nulla
         int nPts = vPL[nPol].GetPointNbr() ;
         if ( nPts < 4)
            continue ;
        // creo un vettore contenente tutti i punti ( cerco nel mentre il segmento più lungo)
         PNTVECTOR vPt ; vPt.reserve( nPts) ;
         int nLongEdge = -1 ;
         double dMaxSqLen = 0. ;
         Point3d myPt ;
         if ( ! vPL[nPol].GetFirstPoint( myPt))
            return false ;
         vPt.emplace_back( myPt) ;
         while ( vPL[nPol].GetNextPoint( myPt)) {
            double dCurrSqDist = SqDist( vPt.back(), myPt) ;
            if ( dCurrSqDist > dMaxSqLen) {
               dMaxSqLen = dCurrSqDist ;
               nLongEdge = int( vPt.size()) - 1 ;
            }
            vPt.emplace_back( myPt) ;
         }
         if ( nLongEdge < 0  ||  nLongEdge > nPts - 2)
            return false ;
        // calcolo il punto iniziale come punto medio del segmento più lungo
         Point3d ptStart = Media( vPt[nLongEdge], vPt[nLongEdge + 1]) ;
         ChangePolyLineStart( vPL[nPol], ptStart, 10 * EPS_SMALL) ;
         nPts = vPL[nPol].GetPointNbr() ;

        // creazione vettore di punti avanzati ( tutti validi)
         PNTIVECTOR vAdvPt ; vAdvPt.reserve( nPts) ;
         vPL[nPol].GetFirstPoint( myPt) ;
         vAdvPt.emplace_back( make_pair( myPt, PointType::VALID)) ;
         while ( vPL[nPol].GetNextPoint( myPt))
            vAdvPt.emplace_back( make_pair( myPt, PointType::VALID)) ;

         #if ENABLE_SHARPED_EDGES_DEBUG
           // 3° gruppo, punti della PolyLine
            for ( int i = 0 ; i < int( vAdvPt.size()) ; ++ i) {
               PtrOwner<IGeoPoint3d> myPt( CreateGeoPoint3d()) ;
               myPt->Set( vAdvPt[i].first) ;
               vvpGObj[2].emplace_back( static_cast<IGeoObj*>( Release( myPt))) ;
               vvCol[2].emplace_back( BLACK) ;
            }
         #endif

        // modifico la posizione dei punti
         if ( ! ModifyPolyLineToSharped( vAdvPt, cavTstm, frGrid, dAngTol, dStep, dSegLen))
            return false ;

        // ricostrusico la polyLine a partire da una vuota
         PolyLine PL_New ;
         double dPar = -1 ;
         for ( int j = 0 ; j < int( vAdvPt.size()) ; ++ j) {
           // se punto valido o sharped di angolo interno, lo aggiungo
            if ( j == 0  ||  j == int( vAdvPt.size() - 1)  ||
                 vAdvPt[j].second == PointType::VALID  ||  vAdvPt[j].second == PointType::SHARPED_INT)
               PL_New.AddUPoint( ++ dPar, vAdvPt[j].first) ;
           // se è un punto sharped di un angolo esterno
            else if ( vAdvPt[j].second == PointType::SHARPED_EXT) {
              // se non richiesto il calcolo avanzato del punto a minima distanza dallo spigolo, aggiungo il punto
               if ( ! bAdvCorners)
                  PL_New.AddUPoint( ++ dPar, vAdvPt[j].first) ;
               else {
                 // ricavo il punto a minima distanza dal materiale mediante metodo di bisezione
                  Vector3d vtPrev = ( vAdvPt[j-1].first - vAdvPt[j].first) ;
                  vtPrev.Normalize() ;
                  Vector3d vtAfter = ( vAdvPt[j+1].first - vAdvPt[j].first) ;
                  vtAfter.Normalize() ;
                  Vector3d vtMove = ( vtPrev + vtAfter) ;
                  vtMove.Normalize() ;
                  Point3d ptTest = vAdvPt[j].first + ( 2 * dStep * vtMove) ; // doppio dello step griglia
                  const int MAX_ITER = 20 ;
                  int nCount = 0 ;
                  while ( nCount < MAX_ITER) {
                     double dMove ;
                     ptTest.z += cavTstm.GetToolHeight() ;
                     ptTest.ToGlob( frGrid) ;
                     cavTstm.TestPosition( ptTest, frGrid.VersZ(), frGrid.VersZ(), dMove) ;
                     ptTest.ToLoc( frGrid) ;
                     ptTest.z -= cavTstm.GetToolHeight() ;
                     if ( dMove > EPS_SMALL)
                        ptTest = Media( ptTest, vAdvPt[j].first) ;
                     else
                        break ;
                     ++ nCount ;
                  }
                 // ricavo i punti a minima distanza tra i due segmenti
                  Point3d ptMinDistA, ptMinDistB ;
                  DistPointLine distPtLPrev( ptTest, vAdvPt[j-1].first, vAdvPt[j].first) ;
                  DistPointLine distPtLSucc( ptTest, vAdvPt[j].first, vAdvPt[j+1].first) ;
                  distPtLPrev.GetMinDistPoint( ptMinDistA) ;
                  distPtLSucc.GetMinDistPoint( ptMinDistB) ;
                  PL_New.AddUPoint( ++ dPar, ptMinDistA) ;
                  PL_New.AddUPoint( ++ dPar, ptTest) ;
                  PL_New.AddUPoint( ++ dPar, ptMinDistB) ;

                  #if ENABLE_SHARPED_EDGES_DEBUG
                   // 5° gruppo, punti minDist A e B e ptMid
                    PtrOwner<IGeoPoint3d> myPtA( CreateGeoPoint3d()) ;
                    myPtA->Set( ptMinDistA) ;
                    vvpGObj[4].emplace_back( static_cast<IGeoObj*>( Release( myPtA))) ;
                    vvCol[4].emplace_back( BROWN) ;
                    PtrOwner<IGeoPoint3d> myPtB( CreateGeoPoint3d()) ;
                    myPtB->Set( ptMinDistB) ;
                    vvpGObj[4].emplace_back( static_cast<IGeoObj*>( Release( myPtB))) ;
                    vvCol[4].emplace_back( BROWN) ;
                    PtrOwner<IGeoPoint3d> myPtMid( CreateGeoPoint3d()) ;
                    myPtMid->Set( ptTest) ;
                    vvpGObj[4].emplace_back( static_cast<IGeoObj*>( Release( myPtMid))) ;
                    vvCol[4].emplace_back( BROWN) ;
                  #endif
               }
            }
         }

        // se la polyLinea si autointerseca, allora non la modifico ( essendo a distanza R dalla
        // trimesh rischio di evere Chunk distinti uniti in un unico Chunk
         PtrOwner<ICurveComposite> pCompo( CreateCurveComposite()) ;
         if ( IsNull( pCompo))
            return false ;
         pCompo->FromPolyLine( vPL[nPol]) ;
         SelfIntersCurve SIC( *pCompo) ;
         bool bDiscard = ( SIC.GetCrossOrOverlapIntersCount() > 0) ;
         if ( ! bDiscard) {
           // se la polyLine interseca una delle PolyLine precedenti, non la modifico
            for ( auto& pCompoPL_Prev : vPL_AsCompo) {
               IntersCurveCurve ICC( *pCompo, *pCompoPL_Prev) ;
               if ( ICC.GetCrossOrOverlapIntersCount() > 0) {
                  bDiscard = true ;
                  break ;
               }
            }
         }

        // se nuova approssimazione da conservare, memorizzo i risultati
         if ( ! bDiscard) {
            vPL[nPol] = PL_New ;
            vPL_AsCompo.emplace_back( Release( pCompo)) ;
         }
      }
   }

   #if ENABLE_SHARPED_EDGES_DEBUG
      ICURVEPVECTOR vpCrv ; vpCrv.reserve( vPL.size()) ;
      for ( const auto& PL : vPL) {
        // 6° gruppo, curve composita derivante dalla polyline
         PtrOwner<ICurveComposite> pCrvCompoPoly( CreateCurveComposite()) ;
         pCrvCompoPoly->FromPolyLine( PL) ;
         vvpGObj[5].emplace_back( static_cast<IGeoObj*>( CloneCurveComposite( pCrvCompoPoly))) ;
         vvCol[5].emplace_back( RED) ;
         vpCrv.emplace_back( Release( pCrvCompoPoly)) ;
      }
     // 7° gruppo, controOffset
      ICURVEPOVECTOR vCrvCompoOffs ;
      if ( ! CalcOffsetCurves( vpCrv, vCrvCompoOffs, - ( cavTstm.GetToolRadius() - 2 * EPS_SMALL), ICurve::OFF_EXTEND)) {
         for ( auto& pCrv : vpCrv) { delete( pCrv) ; pCrv = nullptr ;}
         return false ;
      }
      for ( auto& pCrv : vpCrv) { delete( pCrv) ; pCrv = nullptr ;}
      for ( const auto& pCvrOffs : vCrvCompoOffs) {
         vvpGObj[6].emplace_back( static_cast<IGeoObj*>( CloneCurveComposite( pCvrOffs))) ;
         vvCol[6].emplace_back( WHITE) ;
      }
      string myName = Sharped_Edges_Debug_File_Name +
                         to_string( dLevel) + "_" +
                         to_string( dAngTol) + ".nge" ;
       SaveGeoObj( vvpGObj, vvCol, myName) ;
   #endif

   return true ;
}

//----------------------------------------------------------------------------
static bool
MarchingSquares( const DBLVECTOR& vdGrid, int nStepX, int nStepY, double dStep, double dRad,
                 double dOffsR, double dLevel, const CAvToolSurfTm& cavTstm, const Frame3d &frGrid,
                 double dDimZ, double dCalcLevel, bool bTool, POLYLINEVECTOR& vPL)
{
   #if ENABLE_PROCESS_SQUARE_DEBUG  ||  ENABLE_BISECTION_DEBUG
      vector<IGeoObj*> VT_GO ;
      vector<Color> VT_CO ;
   #endif

  // Analizzo gli edge da cui passano le curve cercate
   unordered_map<int, Point3d> umEdgePnt( 4 * ( nStepX + nStepY)) ;
   INTVECTOR vEdgeInd ;
   vEdgeInd.reserve( 4 * ( nStepX + nStepY)) ;
   for ( int j = 0 ; j < nStepY ; ++ j) {
      for ( int i = 0 ; i < nStepX ; ++ i) {
        // indici dei vertici nella griglia
         int nInd0 = i + j * ( nStepX + 1) ;
         int nInd1 = ( i + 1) + j * ( nStepX + 1) ;
         int nInd2 = ( i + 1) + ( j + 1) * ( nStepX + 1) ;
         int nInd3 = i + ( j + 1) * ( nStepX + 1) ;
        // flag del quadrato
         bool bUp0 = ( vdGrid[nInd0] > dLevel) ;
         bool bUp1 = ( vdGrid[nInd1] > dLevel) ;
         bool bUp2 = ( vdGrid[nInd2] > dLevel) ;
         bool bUp3 = ( vdGrid[nInd3] > dLevel) ;

         #if ENABLE_BISECTION_DEBUG
            PtrOwner<IGeoPoint3d> myPt( CreateGeoPoint3d()) ;
            myPt->Set( Point3d( i * dStep, j * dStep, dDimZ)) ;
            VT_GO.emplace_back( static_cast<IGeoObj*>( Release( myPt))) ;
            VT_CO.emplace_back( bUp0 ? BLUE : RED) ;
            myPt.Set( CreateGeoPoint3d()) ;
            myPt->Set( Point3d( ( i + 1) * dStep, j * dStep, dDimZ)) ;
            VT_GO.emplace_back( static_cast<IGeoObj*>( Release( myPt))) ;
            VT_CO.emplace_back( bUp1 ? BLUE : RED) ;
            myPt.Set( CreateGeoPoint3d()) ;
            myPt->Set( Point3d( ( i + 1) * dStep, ( j + 1) * dStep, dDimZ)) ;
            VT_GO.emplace_back( static_cast<IGeoObj*>( Release( myPt))) ;
            VT_CO.emplace_back( bUp2 ? BLUE : RED) ;
            myPt.Set( CreateGeoPoint3d()) ;
            myPt->Set( Point3d( i * dStep, ( j + 1) * dStep, dDimZ)) ;
            VT_GO.emplace_back( static_cast<IGeoObj*>( Release( myPt))) ;
            VT_CO.emplace_back( bUp3 ? BLUE : RED) ;
         #endif

        // se tutti uguali, passo al successivo
         if ( bUp0 == bUp1  &&  bUp0 == bUp2  &&  bUp0 == bUp3)
            continue ;

         #if ENABLE_PROCESS_SQUARE_DEBUG
            PtrOwner<IGeoPoint3d> myPt( CreateGeoPoint3d()) ;
            myPt->Set( Point3d( i * dStep, j * dStep, dDimZ)) ;
            VT_GO.emplace_back( static_cast<IGeoObj*>( Release( myPt))) ;
            VT_CO.emplace_back( bUp0 ? BLUE : RED) ;
            myPt.Set( CreateGeoPoint3d()) ;
            myPt->Set( Point3d( ( i + 1) * dStep, j * dStep, dDimZ)) ;
            VT_GO.emplace_back( static_cast<IGeoObj*>( Release( myPt))) ;
            VT_CO.emplace_back( bUp1 ? BLUE : RED) ;
            myPt.Set( CreateGeoPoint3d()) ;
            myPt->Set( Point3d( ( i + 1) * dStep, ( j + 1) * dStep, dDimZ)) ;
            VT_GO.emplace_back( static_cast<IGeoObj*>( Release( myPt))) ;
            VT_CO.emplace_back( bUp2 ? BLUE : RED) ;
            myPt.Set( CreateGeoPoint3d()) ;
            myPt->Set( Point3d( i * dStep, ( j + 1) * dStep, dDimZ)) ;
            VT_GO.emplace_back( static_cast<IGeoObj*>( Release( myPt))) ;
            VT_CO.emplace_back( bUp3 ? BLUE : RED) ;
         #endif

        // verifico quali calcolare
         if ( bUp0 != bUp1) {
            int nKey = 2 * nInd0 ;
            if ( umEdgePnt.find( nKey) == umEdgePnt.end()) {
               umEdgePnt[ nKey] = P_INVALID ;
               vEdgeInd.emplace_back( nKey) ;
            }
         }
         if ( bUp3 != bUp2) {
            int nKey = 2 * nInd3 ;
            if ( umEdgePnt.find( nKey) == umEdgePnt.end()) {
               umEdgePnt[ nKey] = P_INVALID ;
               vEdgeInd.emplace_back( nKey) ;
            }
         }
         if ( bUp0 != bUp3) {
            int nKey = 2 * nInd0 + 1 ;
            if ( umEdgePnt.find( nKey) == umEdgePnt.end()) {
               umEdgePnt[ nKey] = P_INVALID ;
               vEdgeInd.emplace_back( nKey) ;
            }
         }
         if ( bUp1 != bUp2) {
            int nKey = 2 * nInd1 + 1 ;
            if ( umEdgePnt.find( nKey) == umEdgePnt.end()) {
               umEdgePnt[ nKey] = P_INVALID ;
               vEdgeInd.emplace_back( nKey) ;
            }
         }
      }
   }
  // Numero di edge da valutare
   int nEdgeCnt = int( vEdgeInd.size()) ;
  // Recupero il numero massimo di thread concorrenti
   int nThreadMax = thread::hardware_concurrency() ;
   bool bOk = true ;
  // Se un solo thread o pochi punti
   if ( nThreadMax <= 1  ||  nEdgeCnt < 50) {
      TestSubEdges( umEdgePnt, vEdgeInd, 0, nEdgeCnt - 1, vdGrid, nStepX, dStep, dLevel, dOffsR, cavTstm, frGrid) ;
   }
  // altrimenti
   else {
      const int MAX_PARTS = 32 ;
      INTINTVECTOR vFstLst( MAX_PARTS) ;
     // calcolo le parti del vettore
      int nPartCnt = min( nThreadMax, MAX_PARTS) ;
      int nPartDim = nEdgeCnt / nPartCnt + 1 ;
      for ( int i = 0 ; i < nPartCnt ; ++ i) {
         vFstLst[i].first = i * nPartDim ;
         vFstLst[i].second = min( ( i + 1) * nPartDim, nEdgeCnt) - 1 ;
      }
     // processo le parti
      future<bool> vRes[MAX_PARTS] ;
      for ( int i = 0 ; i < nPartCnt ; ++ i)
         vRes[i] = async( launch::async, &TestSubEdges, ref( umEdgePnt), cref( vEdgeInd), vFstLst[i].first, vFstLst[i].second,
                          cref( vdGrid), nStepX, dStep, dLevel, dOffsR, cref( cavTstm), cref( frGrid)) ;
     // attendo i risultati
      int nFin = 0 ;
      while ( nFin < nPartCnt) {
         for ( int i = 0 ; i < nPartCnt ; ++ i) {
            if ( vRes[i].valid()  &&  vRes[i].wait_for( chrono::nanoseconds{ 1}) == future_status::ready) {
               bOk = vRes[i].get() && bOk ;
               ++ nFin ;
            }
         }
      }
   }

  // Predispongo il concatenamento
   BIPNTVECTOR vBiPnt ;
   vBiPnt.reserve( 2 * ( nStepX + nStepY)) ;
   ChainCurves chainC ;
   chainC.Init( false, EPS_SMALL, 2 * ( nStepX + nStepY)) ;
  // Ciclo sui quadrati da analizzare
   for ( int j = 0 ; j < nStepY ; ++ j) {
      for ( int i = 0 ; i < nStepX ; ++ i) {
        // indici dei vertici nella griglia
         int nInd0 = i + j * ( nStepX + 1) ;
         int nInd1 = ( i + 1) + j * ( nStepX + 1) ;
         int nInd2 = ( i + 1) + ( j + 1) * ( nStepX + 1) ;
         int nInd3 = i + ( j + 1) * ( nStepX + 1) ;
        // flag del quadrato
         int nFlag = ( vdGrid[nInd0] > dLevel ? 1 : 0) +
                     ( vdGrid[nInd1] > dLevel ? 2 : 0) +
                     ( vdGrid[nInd2] > dLevel ? 4 : 0) +
                     ( vdGrid[nInd3] > dLevel ? 8 : 0) ;
        // se quadrato con vertici tutti dello stesso tipo, passo al successivo
         if ( nFlag == 0  ||  nFlag == 15)
            continue ;
        // chiavi
         int vKey[4] = { 2 * nInd0, 2 * nInd1 + 1, 2 * nInd3, 2 * nInd0 + 1} ;
        // calcolo segmenti da inserire
         int nI1s, nI1e, nI2s, nI2e ;
         int nSegCnt = ProcessSquare( nFlag, dLevel, vdGrid[nInd0], vdGrid[nInd1], vdGrid[nInd2],
                                      vdGrid[nInd3], nI1s, nI1e, nI2s, nI2e) ;
         if ( nSegCnt == -1)
            return false ;
         else if ( nSegCnt == 1) {
            Point3d ptL1s = umEdgePnt.find( vKey[nI1s])->second ;
            Point3d ptL1e = umEdgePnt.find( vKey[nI1e])->second ;

            #if ENABLE_BISECTION_DEBUG
               if ( SqDist( Point3d( i * dStep, j * dStep, dDimZ), Point3d( 445, 190, dDimZ)) < 10) {
                  PtrOwner<IGeoPoint3d> myPt( CreateGeoPoint3d()) ;
                  myPt->Set( Point3d( ptL1s + frGrid.VersZ() * 10)) ;
                  VT_GO.emplace_back( static_cast<IGeoObj*>( Release( myPt))) ;
                  VT_CO.emplace_back( WHITE) ;
                  myPt.Set( CreateGeoPoint3d()) ;
                  myPt->Set( Point3d( ptL1e + frGrid.VersZ() * 10)) ;
                  VT_GO.emplace_back( static_cast<IGeoObj*>( Release( myPt))) ;
                  VT_CO.emplace_back( WHITE) ;
               }
            #endif

            vBiPnt.emplace_back( ptL1s, ptL1e) ;
            Vector3d vtDir1 = ptL1e - ptL1s ; vtDir1.Normalize() ;
            chainC.AddCurve( int( vBiPnt.size()), ptL1s, vtDir1, ptL1e, vtDir1) ;

            #if ENABLE_PROCESS_SQUARE_DEBUG
               PtrOwner<ICurveLine> pMyLine( CreateCurveLine()) ;
               pMyLine->Set( ptL1s, ptL1e) ;
               VT_GO.emplace_back( static_cast<IGeoObj*>( Release( pMyLine))) ;
               VT_CO.emplace_back( Color( double( rand()) / RAND_MAX, double( rand()) / RAND_MAX, double( rand()) / RAND_MAX, 1.)) ;
            #endif
         }
         else if ( nSegCnt == 2) {
            Point3d ptL1s = umEdgePnt.find( vKey[nI1s])->second ;
            Point3d ptL1e = umEdgePnt.find( vKey[nI1e])->second ;
            vBiPnt.emplace_back( ptL1s, ptL1e) ;
            Vector3d vtDir1 = ptL1e - ptL1s ; vtDir1.Normalize() ;
            chainC.AddCurve( int( vBiPnt.size()), ptL1s, vtDir1, ptL1e, vtDir1) ;
            Point3d ptL2s = umEdgePnt.find( vKey[nI2s])->second ;
            Point3d ptL2e = umEdgePnt.find( vKey[nI2e])->second ;
            vBiPnt.emplace_back( ptL2s, ptL2e) ;
            Vector3d vtDir2 = ptL2e - ptL2s ; vtDir2.Normalize() ;
            chainC.AddCurve( int( vBiPnt.size()), ptL2s, vtDir2, ptL2e, vtDir2) ;
            #if ENABLE_PROCESS_SQUARE_DEBUG
               PtrOwner<ICurveLine> pMyLine( CreateCurveLine()) ;
               pMyLine->Set( ptL1s, ptL1e) ;
               VT_GO.emplace_back( static_cast<IGeoObj*>( Release( pMyLine))) ;
               VT_CO.emplace_back( Color( double( rand()) / RAND_MAX, double( rand()) / RAND_MAX, double( rand()) / RAND_MAX, 1.)) ;
               pMyLine.Set( CreateCurveLine()) ;
               pMyLine->Set( ptL2s, ptL2e) ;
               VT_GO.emplace_back( static_cast<IGeoObj*>( Release( pMyLine))) ;
               VT_CO.emplace_back( Color( double( rand()) / RAND_MAX, double( rand()) / RAND_MAX, double( rand()) / RAND_MAX, 1.)) ;
            #endif
         }
      }
   }

   #if ENABLE_PROCESS_SQUARE_DEBUG
      VT_GO.emplace_back( static_cast<IGeoObj*>( cavTstm.GetvStm()[0]->Clone())) ;
      VT_GO.back()->ToLoc( frGrid) ;
      VT_CO.emplace_back( Color( 0., 255., 0., 1.)) ;
      SaveGeoObj( VT_GO, VT_CO, Process_Square_Debug_File_Name) ;
   #endif
   #if ENABLE_BISECTION_DEBUG
      VT_GO.emplace_back( static_cast<IGeoObj*>( cavTstm.GetvStm()[0]->Clone())) ;
      VT_GO.back()->ToLoc( frGrid) ;
      VT_CO.emplace_back( Color( 0., 255., 0., 1.)) ;
      SaveGeoObj( VT_GO, VT_CO, Bisection_Debug_File_Name) ;
   #endif

  // Recupero i contorni
   INTVECTOR vnId ;
   while ( chainC.GetChainFromNear( ORIG, false, vnId)) {
     // creo una composita
      CurveComposite crvCompo ;
      crvCompo.AddPoint( vBiPnt[vnId[0]-1].first) ;
      for ( int i = 0 ; i < int( vnId.size()) ; ++ i) {
         Point3d ptCurr = vBiPnt[vnId[i]-1].second ;
         crvCompo.AddLine( ptCurr) ;
      }
      crvCompo.Close() ;
     // elimino le parti allineate
      crvCompo.MergeCurves( 10 * EPS_SMALL, ANG_TOL_STD_DEG) ;
     // salto i contorni orari con area inferiore al doppio del quadrato
      double dCmpArea ;
      if ( ! crvCompo.GetAreaXY( dCmpArea)  ||  ( dCmpArea < 0  &&  abs( dCmpArea) < 2 * dStep * dStep)) 
         continue ;
     // memorizzo i risultati ottenuti
      vPL.emplace_back( PolyLine()) ;
      crvCompo.ApproxWithLines( 0, 0, ICurve::APL_SPECIAL, vPL.back()) ;
   }

  // le polyline ricavate, definendo dei contorni di regioni, sono chiuse ; cerco di ricostruire gli spigoli vivi
   POLYLINEVECTOR vPL_Sharped( vPL.size()) ;
   for ( int i = 0 ; i < int( vPL_Sharped.size()) ; ++ i)
      vPL_Sharped[i] = vPL[i] ;
   if ( TestEdgesClosedPolyLines( vPL_Sharped, cavTstm, nStepX, nStepY, frGrid, dDimZ, dLevel, dStep, ANG_TOL_STD_DEG, false))
      swap( vPL_Sharped, vPL) ;

  // se vettore di PolyLine vuoto, non faccio nulla
   if ( vPL.empty())
      return true ;

  // se non ho impostato un utensile, devo effettuare un contro-offset
   if ( ! bTool) {
      ICURVEPOVECTOR vpOwCrv ; vpOwCrv.reserve( vPL.size()) ;
      ICURVEPVECTOR vpCrv ; vpCrv.reserve( vPL.size()) ;
      for ( const auto& PL : vPL) {
         PtrOwner<ICurveComposite> pCrvCompoPoly( CreateCurveComposite()) ;
         pCrvCompoPoly->FromPolyLine( PL) ;
         vpOwCrv.emplace_back( Release( pCrvCompoPoly)) ;
         vpCrv.emplace_back( vpOwCrv.back()) ;
      }
     // calcolo l'Offset
      ICURVEPOVECTOR vpCrvOffs ;
      if ( ! CalcOffsetCurves( vpCrv, vpCrvOffs, - ( dRad - 2 * EPS_SMALL), ICurve::OFF_EXTEND))
         return false ;
     // resituisco le PolyLine
      POLYLINEVECTOR vPL_Offs ; vPL_Offs.reserve( vpCrvOffs.size()) ;
      for ( auto& pCrv : vpCrvOffs) {
         if ( pCrv != nullptr  &&  pCrv->IsValid()) {
            PolyLine myPolyLine ;
            pCrv->ApproxWithLines( 10 * EPS_SMALL, ANG_TOL_STD_DEG, ICurve::APL_STD, myPolyLine) ;
            vPL_Offs.emplace_back( myPolyLine) ;
         }
      }
      swap( vPL_Offs, vPL) ;
   }

   #if ENABLE_BISECTION_DEBUG
      VT_GO.emplace_back( static_cast<IGeoObj*>( cavTstm.GetvStm()[0]->Clone())) ;
      VT_GO.back()->ToLoc( frGrid) ;
      VT_CO.emplace_back( GRAY) ;
      for ( int i = 0 ; i < int( vPL.size()) ; ++ i) {
         PtrOwner<ICurveComposite> pCompo( CreateCurveComposite()) ;
         pCompo->FromPolyLine( vPL[i]) ;
         VT_GO.emplace_back( static_cast<IGeoObj*>( pCompo->Clone())) ;
         VT_CO.emplace_back( WHITE) ;
      }
      SaveGeoObj( VT_GO, VT_CO, Bisection_Debug_File_Name) ;
   #endif

   return true ;
}


//----------------------------------------------------------------------------
//  Silhouette rispetto ad una direzione e sopra una quota di una superficie TriMesh
//  ( dTol è il passo della griglia di campionamento e il raggio del cilindro di prova)
//----------------------------------------------------------------------------
bool
CAvSilhouetteSurfTm( const ISurfTriMesh& Stm, const Plane3d& plPlane, double dTol, POLYLINEVECTOR& vPL)
{
  // verifico superficie
   if ( &Stm == nullptr  ||  ! Stm.IsValid())
      return false ;
  // verifico piano
   if ( &plPlane == nullptr  ||  plPlane.GetVersN().IsSmall())
      return false ;
  // verifico tolleranza
   dTol = max( dTol, 100 * EPS_SMALL) ;
  // verifico parametri di ritorno
   if ( &vPL == nullptr)
      return false ;
   vPL.clear() ;

  // sistema di riferimento nel piano
   Frame3d frGrid ;
   if ( ! frGrid.Set( plPlane.GetPoint(), plPlane.GetVersN()))
      return false ;

  // bounding box della superficie nel riferimento del piano
   BBox3d b3Surf ;
   if ( ! Stm.GetBBox( GetInvert( frGrid), b3Surf, BBF_STANDARD))
      return false ;
  // calcolo dati della griglia
   const double EXTRA_XY = 1.5 * dTol ;
   const double EXTRA_Z = 10 * dTol ;
   b3Surf.Expand( EXTRA_XY, EXTRA_XY, EXTRA_Z) ;
   int nStepX = int( ceil( b3Surf.GetDimX() / dTol)) ;
   int nStepY = int( ceil( b3Surf.GetDimY() / dTol)) ;
   frGrid.ChangeOrig( GetToGlob( b3Surf.GetMin(), frGrid)) ;
   double dDimZ = b3Surf.GetDimZ() ;
   double dLevelOffs = - b3Surf.GetMin().z ;

  // calcolo dei punti della griglia (sul top del cilindro)
   PNTUVECTOR vPntM( ( nStepX + 1) * ( nStepY + 1)) ;
   for ( int j = 0 ; j <= nStepY ; ++ j) {
      for ( int i = 0 ; i <= nStepX ; ++ i) {
         int nInd = i + j * ( nStepX + 1) ;
         Point3d ptP = GetToGlob( Point3d( i * dTol, j * dTol, dDimZ), frGrid) ;
         vPntM[nInd] = { ptP, 0.} ;
      }
   }

  // esecuzione della verifica
   double dRad = SQRT1_2 * dTol ;
   CAvToolSurfTm cavTstm ;
   cavTstm.SetStdTool( dDimZ, dRad, 0) ;
   cavTstm.SetSurfTm( Stm) ;
   if ( ! cavTstm.TestSeries( vPntM, frGrid.VersZ(), frGrid.VersZ(), -1))
      return false ;

  // griglia degli spostamenti
   DBLVECTOR vdGrid( ( nStepX + 1) * ( nStepY + 1)) ;
   for ( int j = 0 ; j <= nStepY ; ++ j) {
      for ( int i = 0 ; i <= nStepX ; ++ i) {
         int nInd = i + j * ( nStepX + 1) ;
         vdGrid[nInd] = vPntM[nInd].second ;
      }
   }

  // calcolo della silhouette con il metodo MarchingSquares
   double dLevel = dLevelOffs ;
   if ( ! MarchingSquares( vdGrid, nStepX, nStepY, dTol, dRad, 0., dLevel, cavTstm, frGrid, -1., -1., false, vPL))
      return false ;
  // riporto nella corretta posizione le curve trovate
   for ( auto& PL : vPL)
      PL.ToGlob( frGrid) ;

   return true ;
}


//----------------------------------------------------------------------------
//  Silhouette rispetto ad una direzione e sopra diverse quote di una superficie TriMesh
//----------------------------------------------------------------------------
ICAvParSilhouettesSurfTm*
CreateCAvParSilhouettesSurfTm( void)
{
   return static_cast<ICAvParSilhouettesSurfTm*> ( new(nothrow) CAvParSilhouettesSurfTm) ;
}

//----------------------------------------------------------------------------
CAvParSilhouettesSurfTm::CAvParSilhouettesSurfTm( void)
   : m_dTol( 100 * EPS_SMALL), m_nStepX( 0), m_nStepY( 0), m_dRad( m_dTol), m_dOffsR( 0), m_dLevelOffs( 0), m_bGridOk( false)
{
}

//----------------------------------------------------------------------------
bool
CAvParSilhouettesSurfTm::SetData( const CISURFTMPVECTOR& vpStm, const Frame3d& frPlanes, double dTol)
{
   m_vpStm = vpStm ;
   m_frGrid = frPlanes ;
   m_dTol = max( dTol, 100 * EPS_SMALL) ;
   m_dRad = SQRT1_2 * m_dTol ;
   m_dSharpedTol = ANG_TOL_STD_DEG ;
   m_dOffsR = 0. ;
   m_nStepX = 0 ;
   m_nStepY = 0 ;
   m_dDimZ = 0 ;
   m_dLevelOffs = 0 ;
   m_dCornRad = 0. ;
   m_dMaxMat = INFINITO ;
   m_dOffsR = 0. ;
   m_dSideAng = 0. ;
   m_dMaxDepth = 0. ;
   m_bGridOk = false ;
   m_bTool = false ;
   return true ;
}

//----------------------------------------------------------------------------
bool
CAvParSilhouettesSurfTm::SetData( const CISURFTMPVECTOR& vpStm, const Frame3d& frPlanes, double dTol,
                                  double dSideAng, double dDiam, double dCornRad, double dMaxMat, double dOffsR,
                                  double dMaxDepth)
{
   m_vpStm = vpStm ;
   m_frGrid = frPlanes ;
   m_dTol = max( dTol, 100 * EPS_SMALL) ;
   m_nStepX = 0 ;
   m_nStepY = 0 ;
   m_dDimZ = 0 ;
   m_dLevelOffs = 0 ;
   m_dRad = dDiam / 2. ;
   m_dCornRad = dCornRad ;
   m_dMaxMat = dMaxMat ;
   m_dOffsR = dOffsR ;
   m_dSideAng = dSideAng ;
   m_dMaxDepth = dMaxDepth ;
   m_bGridOk = false ;
   m_bTool = true ;
   return true ;
}

//----------------------------------------------------------------------------
bool
CAvParSilhouettesSurfTm::Prepare( void)
{
  // ingombro delle superfici nel riferimento dei piani
   BBox3d b3All ;
   Frame3d frInv = GetInvert( m_frGrid) ;
   for ( auto pStm : m_vpStm) {
      BBox3d b3Surf ;
      if ( ! pStm->GetBBox( frInv, b3Surf, BBF_STANDARD))
         return false ;
      b3All.Add( b3Surf) ;
   }
  // espansione in Z del Box
   const double EXTRA_Z = 10 * m_dTol + ( m_dMaxDepth > EPS_SMALL ? max( 0., m_dMaxDepth - b3All.GetDimZ()) : 0) ;
   b3All.Expand( 0., 0., EXTRA_Z) ;
   m_dDimZ = b3All.GetDimZ() ;

  // set utensile corrente per calcolo delle collisioni
   double dExtraXY = m_dRad + m_dOffsR ;
   if ( m_dSideAng < EPS_ANG_SMALL)
      m_cavTstm.SetStdTool( m_dDimZ + m_dOffsR, m_dRad + m_dOffsR, m_dCornRad + m_dOffsR) ;
   else {
      double dDeltaRad ;
      double dSideAngRad = m_dSideAng * DEGTORAD ;
      if ( m_dSideAng > 0) {
         if ( m_dCornRad < EPS_SMALL)
            dDeltaRad = m_dMaxMat * tan( dSideAngRad) ;
         else
            dDeltaRad = ( m_dCornRad * cos( dSideAngRad) +
                        ( m_dMaxMat + m_dCornRad * ( sin( dSideAngRad) - 1)) * tan( dSideAngRad)) ;
      }
      else
         dDeltaRad = tan( dSideAngRad) * m_dMaxMat ;

      double dStemRad = m_dRad + dDeltaRad ;
      double dTipRad = m_dRad ;
      dExtraXY = dStemRad + m_dOffsR ;
      m_cavTstm.SetAdvTool( m_dDimZ + m_dOffsR, dStemRad + m_dOffsR, m_dMaxMat, dTipRad + m_dOffsR, m_dCornRad + m_dOffsR) ;
   }

  // calcolo dati della griglia
   const double EXTRA_XY = ( m_bTool ? dExtraXY + 2. : 1.5 * m_dTol) ;
   b3All.Expand( EXTRA_XY, EXTRA_XY, 0.) ;
   m_frGrid.ChangeOrig( GetToGlob( b3All.GetMin(), m_frGrid)) ;
   m_dLevelOffs = - b3All.GetMin().z ;
   m_nStepX = int( ceil( b3All.GetDimX() / m_dTol)) ;
   m_nStepY = int( ceil( b3All.GetDimY() / m_dTol)) ;

  // calcolo dei punti della griglia (sul top del cilindro)
   PNTUVECTOR vPntM( ( m_nStepX + 1) * ( m_nStepY + 1)) ;
   for ( int j = 0 ; j <= m_nStepY ; ++ j) {
      for ( int i = 0 ; i <= m_nStepX ; ++ i) {
         int nInd = i + j * ( m_nStepX + 1) ;
         Point3d ptP = GetToGlob( Point3d( i * m_dTol, j * m_dTol, m_dDimZ), m_frGrid) ;
         vPntM[nInd] = { ptP, 0.} ;
      }
   }

  // esecuzione verifica della collisione
   if ( m_vpStm.empty()  ||  ! m_cavTstm.SetSurfTm( *( m_vpStm[0])))
      return false ;
   for ( int k = 1 ; k < int( m_vpStm.size()) ; ++ k)
      m_cavTstm.AddSurfTm( *( m_vpStm[k])) ;
   if ( ! m_cavTstm.TestSeries( vPntM, m_frGrid.VersZ(), m_frGrid.VersZ(), -1))
      return false ;

  // griglia degli spostamenti
   m_vdGrid.clear() ;
   m_vdGrid.resize( ( m_nStepX + 1) * ( m_nStepY + 1)) ;
   for ( int j = 0 ; j <= m_nStepY ; ++ j) {
      for ( int i = 0 ; i <= m_nStepX ; ++ i) {
         int nInd = i + j * ( m_nStepX + 1) ;
         m_vdGrid[nInd] = vPntM[nInd].second ;
      }
   }

  // disegno la griglia colorando i punti in base alle altezze trovate
   #if ENABLE_COLORED_GRID_DEBUG
      vector<IGeoObj*> VT_GO ;
      vector<Color> VT_CO ;
      for ( int i = 0 ; i < int( vPntM.size()) ; ++ i) {
         PtrOwner<IGeoPoint3d> myPt( CreateGeoPoint3d()) ;
         myPt->Set( vPntM[i].first) ;
         double myAngle = 240 + 120 * ( vPntM[i].second / m_dDimZ) ;
         VT_GO.emplace_back( static_cast<IGeoObj*>( Release( myPt))) ;
         VT_CO.emplace_back( Color( GetColorFromHSV( HSV( myAngle, 1., 1.)))) ;
      }
      SaveGeoObj( VT_GO, VT_CO, Colored_Grid_Debug_File_Name) ;
   #endif

   return true ;
}

//----------------------------------------------------------------------------
bool
CAvParSilhouettesSurfTm::GetSilhouette( double dLevel, POLYLINEVECTOR& vPL)
{
  // reset risultato
   vPL.clear() ;

  // se necessario eseguo i calcoli preparatori
   if ( ! m_bGridOk  &&  ! Prepare())
      return false ;
   m_bGridOk = true ;

   dLevel += m_dLevelOffs ;
  // controllo se ho impostato un utensile, in caso negativo, la silhouette richiede un contro-offset
   double dCalcLevel = max( dLevel, 0.) ;
  // calcolo della silhouette con il metodo MarchingSquare
   if ( ! MarchingSquares( m_vdGrid, m_nStepX, m_nStepY, m_dTol, m_dRad, m_dOffsR, dCalcLevel, m_cavTstm, m_frGrid,
                           m_dDimZ, dCalcLevel, m_bTool, vPL))
      return false ;

  // riporto nella corretta posizione le curve trovate
   for ( auto& PL : vPL) {
      if ( dLevel < dCalcLevel - EPS_SMALL)
         PL.Translate( Vector3d{ 0, 0, dLevel - dCalcLevel}) ;
      PL.ToGlob( m_frGrid) ;
   }
   return true ;
}