//----------------------------------------------------------------------------
//  EgalTech 2020-2020
//----------------------------------------------------------------------------
// File : SurfBezier.cpp       Data : 11.08.20          Versione : 2.2h2
// Contenuto : Implementazione della classe Superfici Bezier.
//
//
//
// Modifiche : 22.03.20 DS Creazione modulo.
//             11.08.20 DS Trasformata in MultiPatch.
//
//----------------------------------------------------------------------------

//--------------------------- Include ----------------------------------------
#include "stdafx.h"
#include "SurfBezier.h"
#include "GeoObjFactory.h"
#include "NgeWriter.h"
#include "NgeReader.h"
#include "Bernstein.h"
#include "CurveBezier.h"
#include "CurveComposite.h"
#include "Tree.h"
#include "Triangulate.h"
#include "SurfTriMesh.h"
#include "/EgtDev/Include/EGkSfrCreate.h"
#include "/EgtDev/Include/EGkStmFromTriangleSoup.h"
#include "/EgtDev/Include/EGkStringUtils3d.h"
#include "/EgtDev/Include/EGkUiUnits.h"
#include "/EgtDev/Include/EgtNumUtils.h"
#include "/EgtDev/Include/EgtPointerOwner.h"
#include "/EgtDev/Include/EGkIntersPlaneSurfTm.h"
#include "/EgtDev/Include/EGkChainCurves.h"
#include "/EgtDev/Include/EGkIntersLineSurfBez.h"
#include "/EgtDev/Include/EGkDistPointSurfTm.h"
#include "/EgtDev/Include/EGkCurveComposite.h"
#include "/EgtDev/Include/EGkGeoPoint3d.h"
#include "/EgtDev/Include/EGkIntervals.h"
#include "/EgtDev/Extern/Eigen/Dense"
#include <limits>

using namespace std ;

//----------------------------------------------------------------------------
GEOOBJ_REGISTER( SRF_BEZIER, NGE_S_BEZ, SurfBezier) ;

//----------------------------------------------------------------------------
SurfBezier::SurfBezier( void)
   : m_pSTM( nullptr), m_nStatus( TO_VERIFY), m_nDegU(), m_nDegV(), m_nSpanU(), m_nSpanV(), m_bRat( false),
     m_bTrimmed( false), m_bClosedU( false), m_bClosedV( false), m_pTrimReg(nullptr), m_nTempProp{0,0}, m_dTempParam{0.0,0.0}
{
}

//----------------------------------------------------------------------------
SurfBezier::~SurfBezier( void)
{
   ResetTrimRegion() ;
   m_OGrMgr.Reset() ;
   ResetAuxSurf() ;
}

//----------------------------------------------------------------------------
bool
SurfBezier::Init( int nDegU, int nDegV, int nSpanU, int nSpanV, bool bIsRational)
{
  // verifico validit� grado
   if ( nDegU < 1  ||  nDegU > MAXDEG  ||  nDegV < 1  ||  nDegV > MAXDEG  ||  nSpanU < 1  ||  nSpanV < 1)
      return false ;

  // imposto gradi e flag di razionale
   m_nDegU = nDegU ;
   m_nDegV = nDegV ;
   m_nSpanU = nSpanU ;
   m_nSpanV = nSpanV ;
   m_bRat = bIsRational ;
   m_bTrimmed = false ;
   ResetTrimRegion() ;

  // dimensiono i vettori dei punti e dei pesi
   m_vPtCtrl.assign( GetDim(), ORIG) ;
   if ( bIsRational)
      m_vWeCtrl.assign( GetDim(), 1) ;
   else
      m_vWeCtrl.clear() ;
   m_nStatus = TO_VERIFY ;

  // imposto ricalcolo della grafica
   ResetAuxSurf() ;
   m_OGrMgr.Reset() ;

   return Validate() ;
}

//----------------------------------------------------------------------------
bool
SurfBezier::SetControlPoint( int nInd, const Point3d& ptCtrl)
{
  // verifico validit� indice
   if ( m_nStatus != OK  ||  m_bRat  ||  nInd < 0  ||  nInd >= GetDim())
      return false ;

  // assegno il valore
   m_vPtCtrl[nInd] = ptCtrl ;

  // se razionale, metto il peso a 1
   if ( m_bRat)
      m_vWeCtrl[nInd] = 1 ;

  // imposto ricalcolo della grafica
   ResetAuxSurf() ;
   m_OGrMgr.Reset() ;

   return true ;
}

//----------------------------------------------------------------------------
bool
SurfBezier::SetControlPoint( int nInd, const Point3d& ptCtrl, double dW)
{
  // verifico validit�, razionalit� e indice
   if ( m_nStatus != OK  ||  ! m_bRat  ||  nInd < 0  ||  nInd >= GetDim())
      return false ;

  // verifico che il peso non sia nullo o negativo
   if ( dW < EPS_SMALL)
      return false ;

  // assegno il valore e il peso
   m_vPtCtrl[nInd] = ptCtrl ;
   m_vWeCtrl[nInd] = dW ;

  // imposto ricalcolo della grafica
   ResetAuxSurf() ;
   m_OGrMgr.Reset() ;

   return true ;
}

//----------------------------------------------------------------------------
bool
SurfBezier::SetTrimRegion( ISurfFlatRegion& sfrTrimReg, bool bIntersectOrSubtract)
{
   // controllo se aveo trim precedenti ed eventualmente faccio un'intersezione con lo spazio esistente

  // verifico la regione passata
   if ( &sfrTrimReg == nullptr  ||  ! sfrTrimReg.IsValid())
      return false ;
  // se la normale ha z negativa ribalto la superficie, sennò le operazioni di intersect e subtract non funzionano
   if ( sfrTrimReg.GetNormVersor().z < 0)
      sfrTrimReg.Invert() ;
  // limito la regione allo spazio parametrico della superficie
   PtrOwner< ISurfFlatRegion> pSfrTrim( GetSurfFlatRegionRectangle( SBZ_TREG_COEFF * m_nSpanU, SBZ_TREG_COEFF * m_nSpanV)) ;
   // bIntersectOrSubtract == true per ottenere lo spazio parametrico trimmato devo fare l'INTERSEZIONE tra il rettangolo totale e l'area passata
   if ( bIntersectOrSubtract) {
      if ( IsNull( pSfrTrim)  ||  ! pSfrTrim->Intersect( sfrTrimReg)  ||  ! pSfrTrim->IsValid()) {
         // provo a offsettare il rettangolo parametrico ( ingrandendolo) per vedere se risolvo problemi di intersezione
         pSfrTrim->Offset( 10* EPS_SMALL, ICurve::OFF_EXTEND) ;
         if ( ! pSfrTrim->Intersect( sfrTrimReg)  ||  ! pSfrTrim->IsValid())
            return false ;
         pSfrTrim->Offset( -10* EPS_SMALL, ICurve::OFF_EXTEND) ;
      }
   }
   // bIntersectOrSubtract == false per ottenere lo spazio parametrico trimmato devo fare la SOTTRAZIONE tra il rettangolo totale e l'area passata
   else {
      if ( IsNull( pSfrTrim)  ||  ! pSfrTrim->Subtract( sfrTrimReg)  ||  ! pSfrTrim->IsValid()) {
         // provo a offsettare il rettangolo parametrico ( ingrandendolo) per vedere se risolvo problemi di sottrazione
         pSfrTrim->Offset( 10* EPS_SMALL, ICurve::OFF_EXTEND) ;
         if ( ! pSfrTrim->Subtract( sfrTrimReg)  ||  ! pSfrTrim->IsValid())
            return false ;
         pSfrTrim->Offset( -10* EPS_SMALL, ICurve::OFF_EXTEND) ;
      }
   }
   ResetAuxSurf() ;
  // assegno la regione di trim
   if ( m_pTrimReg != nullptr ) {
      if ( ! m_pTrimReg->Intersect( *pSfrTrim))
         return false ;
   }
   else
      m_pTrimReg = GetBasicSurfFlatRegion( Release( pSfrTrim)) ;
   // imposto ricalcolo della grafica
   m_OGrMgr.Reset() ;
   // setto la superficie trimmata
   m_bTrimmed = true ;
   return true ;
}

//----------------------------------------------------------------------------
SurfFlatRegion*
SurfBezier::GetTrimRegion( void) const
{
   if ( ! m_bTrimmed  ||  m_pTrimReg == nullptr )
      return nullptr ;
   return m_pTrimReg ;
}

//----------------------------------------------------------------------------
bool
SurfBezier::GetInfo( int& nDegU, int& nDegV, int& nSpanU, int& nSpanV, bool& bIsRat, bool& bTrimmed) const
{
  // verifico validit� superficie
   if ( m_nStatus != OK)
      return false ;
  // restituisco gradi e flag di razionale
   nDegU = m_nDegU ;
   nDegV = m_nDegV ;
   nSpanU = m_nSpanU ;
   nSpanV = m_nSpanV ;
   bIsRat = m_bRat ;
   bTrimmed = m_bTrimmed ;
   return true ;
}

//----------------------------------------------------------------------------
const Point3d&
SurfBezier::GetControlPoint( int nInd, bool* pbOk) const
{
  // verifico validit� e indice
   if ( m_nStatus != OK  ||  nInd < 0  ||  nInd >= GetDim()) {
      if ( pbOk != NULL)
         *pbOk = false ;
      return ORIG ;
   }
  // ritorno i dati
   if ( pbOk != NULL)
      *pbOk = true ;
   return m_vPtCtrl[nInd] ;
}

//----------------------------------------------------------------------------
double
SurfBezier::GetControlWeight( int nInd, bool* pbOk) const
{
  // verifico validit�, razionalit� e indice
   if ( m_nStatus != OK  ||  ! m_bRat  ||  nInd < 0  ||  nInd >= GetDim()) {
      if ( pbOk != NULL)
         *pbOk = false ;
      return 0 ;
   }
  // ritorno i dati
   if ( pbOk != NULL)
      *pbOk = true ;
   return m_vWeCtrl[nInd] ;
}

//----------------------------------------------------------------------------
bool
SurfBezier::IsAPoint( void) const
{
  // verifico lo stato
   if ( m_nStatus != OK)
      return false ;
  // ciclo sui punti
   for ( int i = 1 ; i < GetDim() ; ++ i) {
      if ( ! AreSamePointApprox( m_vPtCtrl[0], m_vPtCtrl[i]))
         return false ;
   }
   
   return true ;
}

//----------------------------------------------------------------------------
bool
SurfBezier::GetArea( double& dArea) const
{
  // controllo parametro di ritorno
   if ( &dArea == nullptr)
      return false ;
  // inizio con area nulla
   dArea = 0 ;
  // calcolo l'area
   if ( m_pSTM == nullptr)
      if ( ! GetAuxSurf())
         return false ;
   return m_pSTM->GetArea( dArea) ;
}

//----------------------------------------------------------------------------
bool
SurfBezier::GetCentroid( Point3d& ptCen) const
{
  // controllo parametro di ritorno
   if ( &ptCen == nullptr)
      return false ;
  // inizio con centro nell'origine
   ptCen = ORIG ;
  // calcolo il baricentro
   if ( m_pSTM == nullptr)
      if ( ! GetAuxSurf())
         return false ;
   return m_pSTM->GetCentroid( ptCen) ;
}

//----------------------------------------------------------------------------
bool
SurfBezier::GetPointD1D2( double dU, double dV, Side nUs, Side nVs,
                          Point3d& ptPos,
                          Vector3d* pvtDerU, Vector3d* pvtDerV, Vector3d* pvtDerUU, Vector3d* pvtDerVV, Vector3d* pvtDerUV) const
{
  // la curva deve essere validata
   if ( m_nStatus != OK)
      return false ;

  // i parametri U e V devono essere compresi tra 0 e il corrispondente numero di Span
   dU = Clamp( dU, 0., double( m_nSpanU)) ;
   dV = Clamp( dV, 0., double( m_nSpanV)) ;

  // determino gli intervalli di span e riduco i parametri in essi
   int nBsU = min( int( dU), m_nSpanU - 1) ;
   double dLocU = dU - nBsU ;
   if ( abs( dLocU) < 5 * EPS_PARAM  &&  nBsU > 0  &&  nUs == ISurfBezier::FROM_MINUS) {
      -- nBsU ;
      dLocU = 1 ;
   }
   else if ( abs( dLocU) > 1 - 5 * EPS_PARAM  &&  nBsU < m_nSpanU - 1  &&  nUs == ISurfBezier::FROM_PLUS) {
      ++ nBsU ;
      dLocU = 0 ;
   }
   int nOffsU = nBsU * m_nDegU ;
   int nBsV = min( int( dV), m_nSpanV - 1) ;
   double dLocV = dV - nBsV ;
   if ( abs( dLocV) < 5 * EPS_PARAM  &&  nBsV > 0  &&  nVs == ISurfBezier::FROM_MINUS) {
      -- nBsV ;
      dLocV = 1 ;
   }
   else if ( abs( dLocV) > 1 - 5 * EPS_PARAM  &&  nBsV < m_nSpanV - 1  &&  nVs == ISurfBezier::FROM_PLUS) {
      ++ nBsV ;
      dLocV = 0 ;
   }
   int nOffsV = nBsV * m_nDegV ;

  // se forma polinomiale (o integrale)
   if ( ! m_bRat) {

     // calcolo dei polinomi di Bernstein per U di grado opportuno
      DBLVECTOR vBernU( m_nDegU + 1) ;
      GetAllBernstein( dLocU, m_nDegU - 2, vBernU) ;

     //// se richiesto, calcolo della derivata seconda
     // if ( pvtDer2 != nullptr  &&  pvtDer1 != nullptr) {
     //    *pvtDer2 = V_NULL ;
     //    for ( int i = 0 ; i <= m_nDeg - 2 ; ++ i) {
     //       *pvtDer2 += vBern[i] * ( m_vPtCtrl[i+2] + m_vPtCtrl[i] - 2 * m_vPtCtrl[i+1]) ;
     //    }
     //    *pvtDer2 *= m_nDeg * ( m_nDeg - 1) ;
     // }

     // aumento il grado
      IncreaseAllBernsteinOneDegree( dLocU, m_nDegU - 1, vBernU) ;

     // se richiesto, calcolo dei vettori intermedi per la derivata prima rispetto ad U
      VCT3DVECTOR vtTemp1( m_nDegV + 1, V_NULL) ;
      if ( pvtDerU != nullptr) {
         for ( int j = 0 ; j <= m_nDegV ; ++ j) {
            for ( int i = 0 ; i <= m_nDegU - 1 ; ++ i)
               vtTemp1[j] += vBernU[i] * ( m_vPtCtrl[GetInd( nOffsU + i + 1, nOffsV + j)] - m_vPtCtrl[GetInd( nOffsU + i, nOffsV + j)]) ;
            vtTemp1[j] *= m_nDegU ;
         }
      }

     // aumento il grado
      IncreaseAllBernsteinOneDegree( dLocU, m_nDegU, vBernU) ;

     // calcolo dei punti intermedi
      PNTVECTOR ptTemp( m_nDegV + 1, ORIG) ;
      for ( int j = 0 ; j <= m_nDegV ; ++ j) {
         for ( int i = 0 ; i <= m_nDegU ; ++ i)
            ptTemp[j] += vBernU[i] * m_vPtCtrl[GetInd( nOffsU + i, nOffsV + j)] ;
      }

     // calcolo dei polinomi di Bernstein per V di grado opportuno
      DBLVECTOR vBernV( m_nDegV + 1) ;
      GetAllBernstein( dLocV, m_nDegV - 1, vBernV) ;

     // se richiesto, calcolo della derivata prima rispetto a V
      if ( pvtDerV != nullptr) {
         *pvtDerV = V_NULL ;
         for ( int j = 0 ; j <= m_nDegV - 1 ; ++ j)
            *pvtDerV += vBernV[j] * ( ptTemp[j+1] - ptTemp[j]) ;
         *pvtDerV *= m_nDegV ;
      }

     // aumento il grado
      IncreaseAllBernsteinOneDegree( dLocV, m_nDegV, vBernV) ;

     // calcolo del punto
      ptPos = ORIG ;
      for ( int j = 0 ; j <= m_nDegV ; ++ j)
         ptPos += vBernV[j] * ptTemp[j] ;

     // se richiesto, calcolo della derivata prima rispetto a U
      if ( pvtDerU != nullptr) {
         *pvtDerU = V_NULL ;
         for ( int j = 0 ; j <= m_nDegV ; ++ j)
            *pvtDerU += vBernV[j] * vtTemp1[j] ;
      }
   }

  // altrimenti forma razionale
   else {

     // porto i punti in forma omogenea moltiplicandoli per i pesi
      PNTVECTOR vPtWCtrl( GetLocDim()) ;
      DBLVECTOR vWeCtrl( GetLocDim()) ;
      for ( int j = 0 ; j <= m_nDegV ; ++ j) {
         for ( int i = 0 ; i <= m_nDegU ; ++ i) {
            vPtWCtrl[GetLocInd( i, j)] = m_vWeCtrl[GetInd( nOffsU + i, nOffsV + j)] * m_vPtCtrl[GetInd( nOffsU + i, nOffsV + j)] ;
            vWeCtrl[GetLocInd( i, j)] = m_vWeCtrl[GetInd( nOffsU + i, nOffsV + j)] ;
         }
      }

     // calcolo dei polinomi di Bernstein di grado opportuno
      DBLVECTOR vBernU( m_nDegU + 1) ;
      GetAllBernstein( dLocU, m_nDegU - 2, vBernU) ;

     //// se richiesto, calcolo della derivata seconda
     // double dW2 = 0 ;
     // if ( pvtDer2 != nullptr  &&  pvtDer1 != nullptr) {
     //    *pvtDer2 = V_NULL ;
     //    for ( int i = 0 ; i <= m_nDeg - 2 ; ++ i) {
     //       *pvtDer2 += vBern[i] * ( vPtWCtrl[i+2] + vPtWCtrl[i] - 2 * vPtWCtrl[i+1]) ;
     //       dW2 += vBern[i] * ( m_vWeCtrl[i+2] + m_vWeCtrl[i] - 2 * m_vWeCtrl[i+1]) ;
     //    }
     //    *pvtDer2 *= m_nDeg * ( m_nDeg - 1) ;
     //    dW2 *= m_nDeg * ( m_nDeg - 1) ;
     // }

     // aumento il grado
      IncreaseAllBernsteinOneDegree( dLocU, m_nDegU - 1, vBernU) ;

     // se richiesto, calcolo dei vettori intermedi per la derivata prima rispetto ad U
      VCT3DVECTOR vtTemp1( m_nDegV + 1, V_NULL) ;
      DBLVECTOR dTemp1( m_nDegV + 1, 0) ;
      if ( pvtDerU != nullptr) {
         for ( int j = 0 ; j <= m_nDegV ; ++ j) {
            for ( int i = 0 ; i <= m_nDegU - 1 ; ++ i) {
               vtTemp1[j] += vBernU[i] * ( vPtWCtrl[GetLocInd( i + 1, j)] - vPtWCtrl[GetLocInd( i, j)]) ;
               dTemp1[j] += vBernU[i] * ( vWeCtrl[GetLocInd( i + 1, j)] - vWeCtrl[GetLocInd( i, j)]) ;
            }
            vtTemp1[j] *= m_nDegU ;
            dTemp1[j] *= m_nDegU ;
         }
      }

     // aumento il grado
      IncreaseAllBernsteinOneDegree( dLocU, m_nDegU, vBernU) ;

     // calcolo dei punti e pesi intermedi
      PNTVECTOR ptTempW( m_nDegV + 1, ORIG) ;
      DBLVECTOR dTempW( m_nDegV + 1, 0) ;
      for ( int j = 0 ; j <= m_nDegV ; ++ j) {
         for ( int i = 0 ; i <= m_nDegU ; ++ i) {
            ptTempW[j] += vBernU[i] * vPtWCtrl[GetLocInd( i, j)] ;
            dTempW[j] += vBernU[i] * vWeCtrl[GetLocInd( i, j)] ;
         }
      }

     // calcolo dei polinomi di Bernstein per V di grado opportuno
      DBLVECTOR vBernV( m_nDegV + 1) ;
      GetAllBernstein( dLocV, m_nDegV - 1, vBernV) ;

     // se richiesto, calcolo della derivata prima rispetto a V
      double dW1v = 0 ;
      if ( pvtDerV != nullptr) {
         *pvtDerV = V_NULL ;
         for ( int j = 0 ; j <= m_nDegV - 1 ; ++ j) {
            *pvtDerV += vBernV[j] * ( ptTempW[j+1] - ptTempW[j]) ;
            dW1v += vBernV[j] * ( dTempW[j+1] - dTempW[j]) ;
         }
         *pvtDerV *= m_nDegV ;
         dW1v *= m_nDegV ;
      }

     // aumento il grado
      IncreaseAllBernsteinOneDegree( dLocV, m_nDegV, vBernV) ;

     // calcolo del punto
      double dW = 0 ;
      ptPos = ORIG ;
      for ( int j = 0 ; j <= m_nDegV ; ++ j) {
         ptPos += vBernV[j] * ptTempW[j] ;
         dW += vBernV[j] * dTempW[j] ;
      }

     // ritrasformo da forma omogenea a forma standard
      double dInvW = 1 / ( ( dW > EPS_ZERO) ? dW : EPS_ZERO) ;
      ptPos *= dInvW ;

     // se richiesto, calcolo della derivata prima rispetto a U
      if ( pvtDerU != nullptr) {
         *pvtDerU = V_NULL ;
         double dW1u = 0 ;
         for ( int j = 0 ; j <= m_nDegV ; ++ j) {
            *pvtDerU += vBernV[j] * vtTemp1[j] ;
            dW1u += vBernV[j] * dTemp1[j] ;
         }
         Vector3d vtPos( ptPos.x, ptPos.y, ptPos.z) ;
         *pvtDerU = ( *pvtDerU - dW1u * vtPos) * dInvW ; 
      }

     // se richiesto, completo il calcolo della derivata prima rispetto a V
      if ( pvtDerV != nullptr) {
         Vector3d vtPos( ptPos.x, ptPos.y, ptPos.z) ;
         *pvtDerV = ( *pvtDerV - dW1v * vtPos) * dInvW ; 
      }

      //if ( pvtDer1 != nullptr) {
      //   Vector3d vtPos( ptPos.x, ptPos.y, ptPos.z) ;
      //   *pvtDer1 = ( *pvtDer1 - dW1 * vtPos) * dInvW ; 
      //   if ( pvtDer2 != nullptr)
      //      *pvtDer2 = ( *pvtDer2 - 2 * dW1 * *pvtDer1 - dW2 * vtPos) * dInvW ;
      //}
   }

   return true ;
}

//----------------------------------------------------------------------------
bool
SurfBezier::GetPointNrmD1D2( double dU, double dV, Side nUs, Side nVs,
                             Point3d& ptPos, Vector3d& vtN,
                             Vector3d* pvtDerU, Vector3d* pvtDerV, Vector3d* pvtDerUU, Vector3d* pvtDerVV, Vector3d* pvtDerUV) const
{
  // la superficie deve essere validata
   if ( m_nStatus != OK)
      return false ;

  // i parametri U e V devono essere compresi tra 0 e il corrispondente numero di Span
   dU = Clamp( dU, 0., double( m_nSpanU)) ;
   dV = Clamp( dV, 0., double( m_nSpanV)) ;

  // eseguo calcolo del punto con le derivate prime
   Vector3d vtDerU ;
   if ( pvtDerU == nullptr)
      pvtDerU = &vtDerU ;
   Vector3d vtDerV ;
   if ( pvtDerV == nullptr)
      pvtDerV = &vtDerV ;
   GetPointD1D2( dU, dV, nUs, nVs, ptPos, pvtDerU, pvtDerV) ;

  // calcolo la normale e la verifico
   vtN = *pvtDerU ^ *pvtDerV ;
   if ( vtN.Normalize())
      return true ;

  // se solo una delle due derivate � piccola, mi sposto lungo il relativo parametro e uso le tangenti
   if ( pvtDerU->Len() < EPS_SMALL  &&  pvtDerV->Len() > 10 * EPS_SMALL) {
      double dCoeff = ( dU - 1000 * EPS_PARAM < 0. ? 1 : -1) ;
      double dUm = dU + 1000 * EPS_PARAM * dCoeff ;
      Point3d ptTmp ;
      Vector3d vtTmpU, vtTmpV ;
      GetPointD1D2( dUm, dV, nUs, nVs, ptTmp, &vtTmpU, &vtTmpV) ;
      vtN = ( *pvtDerV ^ vtTmpV) * dCoeff ;
      if ( vtN.Normalize())
         return true ;
   }
   if ( pvtDerU->Len() > 10 * EPS_SMALL  &&  pvtDerV->Len() < EPS_SMALL) {
      double dCoeff = ( dV - 1000 * EPS_PARAM < 0. ? 1 : -1) ;
      double dVm = dV + 1000 * EPS_PARAM * dCoeff ;
      Point3d ptTmp ;
      Vector3d vtTmpU, vtTmpV ;
      GetPointD1D2( dU, dVm, nUs, nVs, ptTmp, &vtTmpU, &vtTmpV) ;
      vtN = ( vtTmpU ^ *pvtDerU) * dCoeff ;
      if ( vtN.Normalize())
         return true ;
   }

  // ricalcolo con una piccola variazione di entrambi i parametri 
   double dUm ;
   if ( dU - 100 * EPS_PARAM < 0.)
      dUm = dU + 100 * EPS_PARAM ;
   else
      dUm = dU - 100 * EPS_PARAM ;
   double dVm ;
   if ( dV - 100 * EPS_PARAM < 0.)
      dVm = dV + 100 * EPS_PARAM ;
   else
      dVm = dV - 100 * EPS_PARAM ;
   Point3d ptTmp ;
   GetPointD1D2( dUm, dVm, nUs, nVs, ptTmp, pvtDerU, pvtDerV) ;
   vtN = *pvtDerU ^ *pvtDerV ;
   return vtN.Normalize() ;
}

//----------------------------------------------------------------------------
SurfBezier*
SurfBezier::Clone( void) const
{
   // alloco oggetto
   SurfBezier* pSbz = new( nothrow) SurfBezier ;
   if ( pSbz != nullptr) {
      if ( ! pSbz->CopyFrom( *this)) {
         delete pSbz ;
         return nullptr ;
      }
   }

   return pSbz ;
}

//----------------------------------------------------------------------------
bool
SurfBezier::CopyFrom( const IGeoObj* pGObjSrc)
{
   const SurfBezier* pSbz = GetBasicSurfBezier( pGObjSrc) ;
   if ( pSbz == nullptr)
      return false ;
   return CopyFrom( *pSbz) ;
}

//----------------------------------------------------------------------------
bool
SurfBezier::CopyFrom( const SurfBezier& sbSrc)
{
   if ( &sbSrc == this)
      return true ;
   if ( ! Init( sbSrc.m_nDegU, sbSrc.m_nDegV, sbSrc.m_nSpanU, sbSrc.m_nSpanV, sbSrc.m_bRat))
      return false ;
   m_nStatus = sbSrc.m_nStatus ;
   m_vPtCtrl = sbSrc.m_vPtCtrl ;
   if ( sbSrc.m_bRat)
      m_vWeCtrl = sbSrc.m_vWeCtrl ;
   if ( sbSrc.m_bTrimmed) {
      m_bTrimmed = true ;
      m_pTrimReg = sbSrc.m_pTrimReg->Clone() ;
   }
   m_nTempProp[0] = sbSrc.m_nTempProp[0] ;
   m_nTempProp[1] = sbSrc.m_nTempProp[1] ;
   m_dTempParam[0] = sbSrc.m_dTempParam[0] ;
   m_dTempParam[1] = sbSrc.m_dTempParam[1] ;
   return true ;
}

//----------------------------------------------------------------------------
GeoObjType
SurfBezier::GetType( void) const
{
   return static_cast<GeoObjType>( GEOOBJ_GETTYPE( SurfBezier)) ;
}

//----------------------------------------------------------------------------
const string&
SurfBezier::GetTitle( void) const
{
   static const string sTitle = "SurfBezier" ;
   return sTitle ;
}

//----------------------------------------------------------------------------
bool
SurfBezier::Dump( string& sOut, bool bMM, const char* szNewLine) const
{
  // verifico validit� superficie
   if ( m_nStatus != OK)
      sOut += string( "Status=Invalid") + szNewLine ;
  // area
   double dArea ;
   GetArea( dArea) ;
   sOut += "Area=" + ToString( GetAreaInUiUnits( dArea, bMM),1) + szNewLine ;
  // parametri : flag razionale
   sOut += ( m_bRat ? "Rat" : "Int") ;
  //             flag trimmata
   sOut += ( m_bTrimmed ? " Trim " : " Full") ;
  //             gradi in U e V
   sOut += " DegU=" + ToString( m_nDegU) + " DegV=" + ToString( m_nDegV) ;
  //             pezze in U e V
   sOut += " SpanU=" + ToString( m_nSpanU) + " SpanV=" + ToString( m_nSpanV) + szNewLine ;
  //             ciclo sui punti di controllo ( con pesi se razionale)
   for ( int i = 0 ; i < GetDim() ; ++ i) {
      sOut += "PC(" + ToString( GetInUiUnits( m_vPtCtrl[i], bMM), 3) ;
      if ( m_bRat)
         sOut += "," + ToString( m_vWeCtrl[i], 3) ;
      sOut += string( ")") + szNewLine ;
   }

   return true ;
}

//----------------------------------------------------------------------------
int
SurfBezier::GetNgeId( void) const
{
   return GEOOBJ_GETNGEID( SurfBezier) ;
}

//----------------------------------------------------------------------------
bool
SurfBezier::Save( NgeWriter& ngeOut) const
{
  // flag razionale
   if ( ! ngeOut.WriteBool( m_bRat, ";"))
      return false ;
  // flag trimmata
   if ( ! ngeOut.WriteBool( m_bTrimmed, ";"))
      return false ;
  // gradi
   if ( ! ngeOut.WriteInt( m_nDegU, ",")  ||
        ! ngeOut.WriteInt( m_nDegV, ";", false))
      return false ;
  // pezze
   if ( ! ngeOut.WriteInt( m_nSpanU, ",")  ||
        ! ngeOut.WriteInt( m_nSpanV, ";", true))
      return false ;
  // ciclo sui punti di controllo ( con pesi se razionale)
   for ( int i = 0 ; i < GetDim() ; ++ i) {
      if ( ! m_bRat) {
         if ( ! ngeOut.WritePoint( m_vPtCtrl[i], ";", true))
            return false ;
      }
      else {
         if ( ! ngeOut.WritePointW( m_vPtCtrl[i], m_vWeCtrl[i], ";", true))
            return false ;
      }
   }
  // se trimmata, scrittura della regione
   if ( m_bTrimmed) {
     // recupero il gestore di lettura/scrittura della regione
      const IGeoObjRW* pSFrRW = dynamic_cast<const IGeoObjRW*>( m_pTrimReg) ;
      if ( pSFrRW == nullptr)
         return false ;
     // salvataggio della regione
      if ( ! pSFrRW->Save( ngeOut))
         return false ;
   }

   return true ;
}

//----------------------------------------------------------------------------
bool
SurfBezier::Load( NgeReader& ngeIn)
{
  // imposto ricalcolo della grafica
   ResetAuxSurf() ;
   m_OGrMgr.Reset() ;
  // leggo la prossima linea ( 3 parametri)
  // recupero il flag razionale
   bool bIsRat ;
   if ( ! ngeIn.ReadBool( bIsRat, ";"))
      return false ;
  // recupero il flag trimmata
   bool bTrimmed ;
   if ( ! ngeIn.ReadBool( bTrimmed, ";"))
      return false ;
  // recupero i gradi
   int nDegU, nDegV ;
   if ( ! ngeIn.ReadInt( nDegU, ",")  ||
        ! ngeIn.ReadInt( nDegV, ";"))
      return false ;
  // recupero le pezze
   int nSpanU, nSpanV ;
   if ( ! ngeIn.ReadInt( nSpanU, ",")  ||
        ! ngeIn.ReadInt( nSpanV, ";", true))
      return false ;
  // inizializzo la superficie di Bezier
   if ( ! Init( nDegU, nDegV, nSpanU, nSpanV, bIsRat))
      return false ;
  // se integrale
   if ( ! bIsRat) {
     // recupero e setto punti di controllo
      Point3d ptP ;
      for ( int i = 0 ; i < GetDim() ; ++ i) {
        // leggo la prossima linea ( un punto)
         if ( ! ngeIn.ReadPoint( ptP, ";", true))
            return false ;
        // lo assegno
         if ( ! SetControlPoint( i, ptP))
            return false ;
      }
   }
  // altrimenti razionale
   else {
     // recupero e setto punti di controllo
      Point3d ptP ;
      double  dW ;
      for ( int i = 0 ; i < GetDim() ; ++ i) {
        // leggo la prossima linea ( un punto con peso)
         if ( ! ngeIn.ReadPointW( ptP, dW, ";", true))
            return false ;
        // lo assegno
         if ( ! SetControlPoint( i, ptP, dW))
            return false ;
      }
   }
  // se trimmata, lettura della regione
   if ( bTrimmed) {
     // creo l'oggetto
      ResetTrimRegion() ;
      m_pTrimReg = CreateBasicSurfFlatRegion() ;
      if ( m_pTrimReg == nullptr)
         return false ;
     // ne leggo i dati
      IGeoObjRW* pGObjRW = dynamic_cast<IGeoObjRW*>( m_pTrimReg) ;
      if ( pGObjRW == nullptr  ||  ! pGObjRW->Load( ngeIn))
         return false ;
      m_bTrimmed = true ;
   }

  // eseguo validazione
   return Validate() ;
}

//----------------------------------------------------------------------------
bool
SurfBezier::Validate( void)
{
   if ( m_nStatus == TO_VERIFY)
      m_nStatus = ( ( m_nDegU * m_nDegV > 0  &&  m_vPtCtrl.size() > 0) ? OK : ERR) ;

   return ( m_nStatus == OK) ;
}

//----------------------------------------------------------------------------
bool
SurfBezier::GetLocalBBox( BBox3d& b3Loc, int nFlag) const
{
   // basta approssimato
   if ( ( nFlag & BBF_EXACT) == 0) {
      for ( int i = 0 ; i < GetDim() ; ++ i) {
         Point3d ptTemp = m_vPtCtrl[i] ;
         b3Loc.Add( ptTemp) ;
      }
      return true ;
   }
   // deve essere preciso
   else {
     // verifico esistenza trimesh associata
      if ( m_pSTM == nullptr)
         if ( ! GetAuxSurf())
            return false ;
     // calcolo il box della trimesh
      return m_pSTM->GetLocalBBox( b3Loc, nFlag) ;
   }
}

//----------------------------------------------------------------------------
bool
SurfBezier::GetBBox( const Frame3d& frRef, BBox3d& b3Ref, int nFlag) const
{
   // basta approssimato
   if ( ( nFlag & BBF_EXACT) == 0) {
      for ( int i = 0 ; i < GetDim() ; ++ i) {
         Point3d ptTemp = m_vPtCtrl[i] ;
         ptTemp.ToGlob( frRef) ;
         b3Ref.Add( ptTemp) ;
      }
      return true ;
   }
   // deve essere preciso
   else {
     // verifico esistenza trimesh associata
      if ( m_pSTM == nullptr)
         if ( ! GetAuxSurf())
            return false ;
     // calcolo il box della trimesh
      return m_pSTM->GetBBox( frRef, b3Ref, nFlag) ;
   }
}

//----------------------------------------------------------------------------
bool
SurfBezier::Translate( const Vector3d& vtMove)
{
  // la superficie deve essere validata
   if ( m_nStatus != OK)
      return false ;

  // imposto ricalcolo della grafica
   ResetAuxSurf() ;
   m_OGrMgr.Reset() ;

  // traslo i punti di controllo
   for ( auto& ptP : m_vPtCtrl)
      ptP.Translate( vtMove) ;

   return true ;
}

//----------------------------------------------------------------------------
bool 
SurfBezier::Rotate( const Point3d& ptAx, const Vector3d& vtAx, double dCosAng, double dSinAng)
{
  // la superficie deve essere validata
   if ( m_nStatus != OK)
      return false ;
  // verifico validit� dell'asse di rotazione
   if ( vtAx.IsSmall())
      return false ;

  // imposto ricalcolo della grafica
   ResetAuxSurf() ;
   m_OGrMgr.Reset() ;

  // ruoto i punti di controllo
   for ( auto& ptP : m_vPtCtrl)
      ptP.Rotate( ptAx, vtAx, dCosAng, dSinAng) ;

   return true ;
}

//----------------------------------------------------------------------------
bool 
SurfBezier::Scale( const Frame3d& frRef, double dCoeffX, double dCoeffY, double dCoeffZ)
{
  // la superficie deve essere validata
   if ( m_nStatus != OK)
      return false ;
  // verifico non sia nulla
   if ( abs( dCoeffX) < EPS_ZERO  &&  abs( dCoeffY) < EPS_ZERO  &&  abs( dCoeffZ) < EPS_ZERO)
      return false ;

  // calcolo bbox allineato con riferimento di scalatura e senza tener conto dello spessore
  // lo scalo per verificare se tutto si riduce a un punto o una linea (solo se diretta come assi ref)
   BBox3d b3Ref ;
   if ( ! GetBBox( frRef, b3Ref))
      return false ;
   Vector3d vtDelta = b3Ref.GetMax() - b3Ref.GetMin() ;
   bool bZeroX = ( abs( vtDelta.x * dCoeffX) < EPS_SMALL) ;
   bool bZeroY = ( abs( vtDelta.y * dCoeffY) < EPS_SMALL) ;
   bool bZeroZ = ( abs( vtDelta.z * dCoeffZ) < EPS_SMALL) ;
   if ( ( bZeroX && bZeroY)  ||  ( bZeroX && bZeroZ)  ||  ( bZeroY && bZeroZ))
      return false ;

  // imposto ricalcolo della grafica
   ResetAuxSurf() ;
   m_OGrMgr.Reset() ;

  // scalo i punti di controllo
   for ( auto& ptP : m_vPtCtrl)
      ptP.Scale( frRef, dCoeffX, dCoeffY, dCoeffZ) ;

   return true ;
}

//----------------------------------------------------------------------------
bool 
SurfBezier::Mirror( const Point3d& ptOn, const Vector3d& vtNorm)
{
  // la superficie deve essere validata
   if ( m_nStatus != OK)
      return false ;
  // verifico validit� del piano di specchiatura
   if ( vtNorm.IsSmall())
      return false ;

  // imposto ricalcolo della grafica
   ResetAuxSurf() ;
   m_OGrMgr.Reset() ;

  // specchio i punti di controllo
   for ( auto& ptP : m_vPtCtrl)
      ptP.Mirror( ptOn, vtNorm) ;

  // eseguo invert per mantenere l'orientazione
   return Invert() ;
}

//----------------------------------------------------------------------------
bool 
SurfBezier::Shear( const Point3d& ptOn, const Vector3d& vtNorm, const Vector3d& vtDir, double dCoeff)
{
  // la superficie deve essere validata
   if ( m_nStatus != OK)
      return false ;
  // verifico validit� dei parametri
   if ( vtNorm.IsSmall()  ||  vtDir.IsSmall())
      return false ;

  // imposto ricalcolo della grafica
   ResetAuxSurf() ;
   m_OGrMgr.Reset() ;

  // eseguo scorrimento dei punti di controllo
   for ( auto& ptP : m_vPtCtrl)
      ptP.Shear( ptOn, vtNorm, vtDir, dCoeff) ;

   return true ;
}

//----------------------------------------------------------------------------
bool 
SurfBezier::ToGlob( const Frame3d& frRef)
{
  // la superficie deve essere validata
   if ( m_nStatus != OK)
      return false ;
  // verifico validit� del frame
   if ( frRef.GetType() == Frame3d::ERR)
      return false ;

  // se frame identit�, non devo fare alcunch�
   if ( IsGlobFrame( frRef))
      return true ;

  // imposto ricalcolo della grafica
   ResetAuxSurf() ;
   m_OGrMgr.Reset() ;

  // trasformo i punti di controllo
   for ( auto& ptP : m_vPtCtrl)
      ptP.ToGlob( frRef) ;

   return true ;
}

//----------------------------------------------------------------------------
bool 
SurfBezier::ToLoc( const Frame3d& frRef)
{
  // la superficie deve essere validata
   if ( m_nStatus != OK)
      return false ;
  // verifico validit� del frame
   if ( frRef.GetType() == Frame3d::ERR)
      return false ;

  // se frame identit�, non devo fare alcunch�
   if ( IsGlobFrame( frRef))
      return true ;

  // imposto ricalcolo della grafica
   ResetAuxSurf() ;
   m_OGrMgr.Reset() ;

  // trasformo i punti di controllo
   for ( auto& ptP : m_vPtCtrl)
      ptP.ToLoc( frRef) ;

   return true ;
}

//----------------------------------------------------------------------------
bool 
SurfBezier::LocToLoc( const Frame3d& frOri, const Frame3d& frDest)
{
  // la superficie deve essere validata
   if ( m_nStatus != OK)
      return false ;
  // verifico validit� dei frame
   if ( frOri.GetType() == Frame3d::ERR  ||  frDest.GetType() == Frame3d::ERR)
      return false ;

  // se i due riferimenti coincidono, non devo fare alcunch�
   if ( AreSameFrame( frOri, frDest))
      return true ;

  // imposto ricalcolo della grafica
   ResetAuxSurf() ;
   m_OGrMgr.Reset() ;

  // trasformo i punti di controllo
   for ( auto& ptP : m_vPtCtrl)
      ptP.LocToLoc( frOri, frDest) ;

   return true ;
}

//----------------------------------------------------------------------------
bool
SurfBezier::Invert( void)
{
  // la superficie deve essere validata
   if ( m_nStatus != OK)
      return false ;

  // imposto ricalcolo della grafica
   ResetAuxSurf() ;
   m_OGrMgr.Reset() ;

  // inverto i punti del parametro U
   PNTVECTOR vPtCtrl_inv( m_vPtCtrl.size()) ;
   
   for ( int j = 0 ; j < m_nDegV * m_nSpanV + 1; ++j) {
      for ( int i = 0 ; i < m_nDegU * m_nSpanU + 1; ++i) {
         vPtCtrl_inv[ i + j * ( m_nDegU * m_nSpanU + 1)] = m_vPtCtrl[ ( m_nDegU * m_nSpanU - i) + j * ( m_nDegU * m_nSpanU + 1)] ;
      }
   }
   m_vPtCtrl = vPtCtrl_inv ;

   if ( m_bTrimmed) {
     // inverto la flat region di trim
     //( la specchio rispetto all'asse verticale)
      Point3d pt( m_nSpanU * SBZ_TREG_COEFF/2, 0, 0) ;
      Vector3d vt( 1, 0, 0) ;
      m_pTrimReg->Mirror( pt, vt) ;
   }

   return true ;
}

//----------------------------------------------------------------------------
int
SurfBezier::GetSteps( int nDeg, int nSpan, double dLen, int nQuality) const
{
   const double dCoeff[] = { 0, 1, 2, 2.4, 2.8, 3, 3, 3, 3, 3, 3, 3} ;
   nDeg = Clamp( nDeg, 0, MAXDEG) ;
   nSpan = max( nSpan, 1) ;
   nQuality = Clamp( nQuality, 1, 10) ;
   double dMult = sqrt( max( dLen / nSpan, 10.)) ;
   return ( nSpan * int( dMult * dCoeff[nDeg] * 2 / nQuality)) ;
}

//----------------------------------------------------------------------------
CurveComposite*
SurfBezier::GetCurveOnU( double dV) const
{
  // controlli
   if ( dV < - EPS_PARAM  ||  dV > m_nSpanV + EPS_PARAM)
      return nullptr ;
   dV = Clamp( dV, 0., double( m_nSpanV)) ;

  // determino l'intervallo di span in V e riduco i parametri in essi
   int nBsV = min( int( dV), m_nSpanV - 1) ;
   double dLocV = dV - nBsV ;
   int nOffsV = nBsV * m_nDegV ;

  // se forma polinomiale (o integrale)
   if ( ! m_bRat) {
     // preparazione della curva composita
      PtrOwner<CurveComposite> pCrvCo( CreateBasicCurveComposite()) ;
     // ciclo sugli intervalli
      for ( int k = 0 ; k < m_nSpanU ; ++ k) {
        // preparazione della curva di Bezier
         PtrOwner<CurveBezier> pCbz( CreateBasicCurveBezier()) ;
         if ( IsNull( pCbz)  ||  ! pCbz->Init( m_nDegU, false))
            return nullptr ;
        // calcolo dei polinomi di Bernstein per V di grado opportuno
         DBLVECTOR vBernV( m_nDegV + 1) ;
         GetAllBernstein( dLocV, m_nDegV, vBernV) ;
        // calcolo offset in U
         int nOffsU = k * m_nDegU ;
        // calcolo dei punti di controllo della curva
         for ( int i = 0 ; i <= m_nDegU ; ++ i) {
            Point3d ptP = ORIG ;
            for ( int j = 0 ; j <= m_nDegV ; ++ j)
               ptP += vBernV[j] * m_vPtCtrl[GetInd( nOffsU + i, nOffsV + j)] ;
            pCbz->SetControlPoint( i, ptP) ;
         }
        // inserisco la curva della pezza in quella complessiva
         if ( ! IsNull( pCbz))
            pCrvCo->AddCurve( Release( pCbz)) ;
      }
      return Release( pCrvCo) ;
   }
  // altrimenti forma razionale
   else {
     // preparazione della curva composita
      PtrOwner<CurveComposite> pCrvCo( CreateBasicCurveComposite()) ;
     // ciclo sugli intervalli
      for ( int k = 0 ; k < m_nSpanU ; ++ k) {
        // preparazione della curva di Bezier
         PtrOwner<CurveBezier> pCbz( CreateBasicCurveBezier()) ;
         if ( IsNull( pCbz)  ||  ! pCbz->Init( m_nDegU, true))
            return nullptr ;
        // calcolo dei polinomi di Bernstein per V di grado opportuno
         DBLVECTOR vBernV( m_nDegV + 1) ;
         GetAllBernstein( dLocV, m_nDegV, vBernV) ;
        // calcolo offset in U
         int nOffsU = k * m_nDegU ;
        // porto i punti in forma omogenea moltiplicandoli per i pesi
         PNTVECTOR vPtWCtrl( GetLocDim()) ;
         DBLVECTOR vWeCtrl( GetLocDim()) ;
         for ( int j = 0 ; j <= m_nDegV ; ++ j) {
            for ( int i = 0 ; i <= m_nDegU ; ++ i) {
               vPtWCtrl[GetLocInd( i, j)] = m_vWeCtrl[GetInd( nOffsU + i, nOffsV + j)] * m_vPtCtrl[GetInd( nOffsU + i, nOffsV + j)] ;
               vWeCtrl[GetLocInd( i, j)] = m_vWeCtrl[GetInd( nOffsU + i, nOffsV + j)] ;
            }
         }
        // calcolo dei punti di controllo della curva
         for ( int i = 0 ; i <= m_nDegU ; ++ i) {
            Point3d ptP = ORIG ;
            double dW = 0 ;
            for ( int j = 0 ; j <= m_nDegV ; ++ j) {
               ptP += vBernV[j] * vPtWCtrl[GetLocInd( i, j)] ;
               dW +=  vBernV[j] * vWeCtrl[GetLocInd( i, j)] ;
            }
            double dInvW = 1 / ( ( dW > EPS_ZERO) ? dW : EPS_ZERO) ;
            pCbz->SetControlPoint( i, ptP * dInvW, dW) ;
         }
        // inserisco la curva della pezza in quella complessiva
         if ( ! IsNull( pCbz))
            pCrvCo->AddCurve( Release( pCbz)) ;
      }
      return Release( pCrvCo) ;
   }
}

//----------------------------------------------------------------------------
CurveComposite*
SurfBezier::GetCurveOnV( double dU) const
{
  // controlli
   if ( dU < - EPS_PARAM  ||  dU > m_nSpanU + EPS_PARAM)
      return nullptr ;
   dU = Clamp( dU, 0., double( m_nSpanU)) ;

  // determino l'intervallo di span in U e riduco i parametri in essi
   int nBsU = min( int( dU), m_nSpanU - 1) ;
   double dLocU = dU - nBsU ;
   int nOffsU = nBsU * m_nDegU ;

  // se forma polinomiale (o integrale)
   if ( ! m_bRat) {
     // preparazione della curva composita
      PtrOwner<CurveComposite> pCrvCo( CreateBasicCurveComposite()) ;
     // ciclo sugli intervalli
      for ( int k = 0 ; k < m_nSpanV ; ++ k) {
        // preparazione della curva di Bezier
         PtrOwner<CurveBezier> pCbz( CreateBasicCurveBezier()) ;
         if ( IsNull( pCbz)  ||  ! pCbz->Init( m_nDegV, false))
            return nullptr ;
        // calcolo dei polinomi di Bernstein per U di grado opportuno
         DBLVECTOR vBernU( m_nDegU + 1) ;
         GetAllBernstein( dLocU, m_nDegU, vBernU) ;
        // calcolo offset in V
         int nOffsV = k * m_nDegV ;
        // calcolo dei punti di controllo della curva
         for ( int j = 0 ; j <= m_nDegV ; ++ j) {
            Point3d ptP = ORIG ;
            for ( int i = 0 ; i <= m_nDegU ; ++ i)
               ptP += vBernU[i] * m_vPtCtrl[GetInd( nOffsU + i, nOffsV + j)] ;
            pCbz->SetControlPoint( j, ptP) ;
         }
        // inserisco la curva della pezza in quella complessiva
         if ( ! IsNull( pCbz))
            pCrvCo->AddCurve( Release( pCbz)) ;
      }
      return Release( pCrvCo) ;
   }
  // altrimenti forma razionale
   else {
     // preparazione della curva composita
      PtrOwner<CurveComposite> pCrvCo( CreateBasicCurveComposite()) ;
     // ciclo sugli intervalli
      for ( int k = 0 ; k < m_nSpanV ; ++ k) {
        // preparazione della curva di Bezier
         PtrOwner<CurveBezier> pCbz( CreateBasicCurveBezier()) ;
         if ( IsNull( pCbz)  ||  ! pCbz->Init( m_nDegV, true))
            return nullptr ;
        // calcolo dei polinomi di Bernstein per U di grado opportuno
         DBLVECTOR vBernU( m_nDegU + 1) ;
         GetAllBernstein( dLocU, m_nDegU, vBernU) ;
        // calcolo offset in V
         int nOffsV = k * m_nDegV ;
        // porto i punti in forma omogenea moltiplicandoli per i pesi
         PNTVECTOR vPtWCtrl( GetLocDim()) ;
         DBLVECTOR vWeCtrl( GetLocDim()) ;
         for ( int j = 0 ; j <= m_nDegV ; ++ j) {
            for ( int i = 0 ; i <= m_nDegU ; ++ i) {
               vPtWCtrl[GetLocInd( i, j)] = m_vWeCtrl[GetInd( nOffsU + i, nOffsV + j)] * m_vPtCtrl[GetInd( nOffsU + i, nOffsV + j)] ;
               vWeCtrl[GetLocInd( i, j)] = m_vWeCtrl[GetInd( nOffsU + i, nOffsV + j)] ;
            }
         }
        // calcolo dei punti di controllo della curva
         for ( int j = 0 ; j <= m_nDegV ; ++ j) {
            Point3d ptP = ORIG ;
            double dW = 0 ;
            for ( int i = 0 ; i <= m_nDegU ; ++ i) {
               ptP += vBernU[i] * vPtWCtrl[GetLocInd( i, j)] ;
               dW +=  vBernU[i] * vWeCtrl[GetLocInd( i, j)] ;
            }
            double dInvW = 1 / ( ( dW > EPS_ZERO) ? dW : EPS_ZERO) ;
            pCbz->SetControlPoint( j, ptP * dInvW, dW) ;
         }
        // inserisco la curva della pezza in quella complessiva
         if ( ! IsNull( pCbz))
            pCrvCo->AddCurve( Release( pCbz)) ;
      }
      return Release( pCrvCo) ;
   }
}

//----------------------------------------------------------------------------
CurveComposite*
SurfBezier::GetLoop( int nLoop) const
{
   // Il primo loop sono le 4 isoparametriche di bordo concatenate
   if ( ! m_bTrimmed ) {
      if ( nLoop != 0 )
         return nullptr ;
      // Loop
      PtrOwner<CurveComposite> pLoop( CreateBasicCurveComposite()) ;
      // prima curva isoparametrica in U con V=0
      PtrOwner<CurveComposite> pCrvCoU0( GetCurveOnU( 0)) ;
      if ( ! IsNull( pCrvCoU0)  &&  ! pCrvCoU0->IsAPoint())
         pLoop->AddCurve( Release( pCrvCoU0)) ;
      // seconda curva isoparametrica in V con U=m_nSpanU
      PtrOwner<CurveComposite> pCrvCoV1( GetCurveOnV( m_nSpanU)) ;
      if ( ! IsNull( pCrvCoV1)  &&  ! pCrvCoV1->IsAPoint())
         pLoop->AddCurve( Release( pCrvCoV1)) ;
      // terza curva isoparametrica in U con V=m_nSpanV invertita
      PtrOwner<CurveComposite> pCrvCoU1( GetCurveOnU( m_nSpanV)) ;
      if ( ! IsNull( pCrvCoU1)  &&  ! pCrvCoU1->IsAPoint()) {
         pCrvCoU1->Invert() ;
         pLoop->AddCurve( Release( pCrvCoU1)) ;
      }
      // quarta curva isoparametrica in V con U=0 invertita
      PtrOwner<CurveComposite> pCrvCoV0( GetCurveOnV( 0)) ;
      if ( ! IsNull( pCrvCoV0)  &&  ! pCrvCoV0->IsAPoint()) {
         pCrvCoV0->Invert() ;
         pLoop->AddCurve( Release( pCrvCoV0)) ;
      }
      // se loop chiuso lo restituisco, altrimenti errore
      return ( pLoop->IsClosed() ? Release( pLoop) : nullptr) ;
   }
   // la superficie � trimmata, quindi devo cercare nei vari chunck il loop corrispondente
   else {
      if ( nLoop > m_pTrimReg->GetChunkCount())
         return nullptr ;
      else {
         int nLoopCount = 0 ;
         int nChunck = 0, nLoopLoc = 0;
         INTVECTOR nLoopCountPerChunck ;
         for ( int i = 0 ; i < m_pTrimReg->GetChunkCount()  &&  nLoopCount != nLoop ; ++ i) {
            int nLoopCountLoc = 0 ;
            for ( int j = 0 ; j < m_pTrimReg->GetLoopCount( i) ; ++ j) {
               ++ nLoopCountLoc ;
               ++ nLoopCount ;
               if ( nLoopCount != nLoop ) {
                  nChunck = i ;
                  nLoopLoc = j ;
                  break ;
               }
            }
            nLoopCountPerChunck.push_back( nLoopCountLoc) ;
         }
         if ( nLoopCount < nLoop )
            return nullptr ;
         PtrOwner<CurveComposite> pLoop( GetBasicCurveComposite( m_pTrimReg->GetLoop( nChunck, nLoopLoc))) ;
         return Release( pLoop) ;
      }
   }
}

//----------------------------------------------------------------------------
bool
SurfBezier::GetCurveOnU( double dV, int nStep, PolyLine& plCrvU) const
{
  // controlli
   plCrvU.Clear() ;
   if ( dV < - EPS_PARAM  ||  dV > m_nSpanV + EPS_PARAM)
      return false ;
   dV = Clamp( dV, 0., double( m_nSpanV)) ;
  // recupero la curva
   PtrOwner<CurveComposite> pCrvCo( GetCurveOnU( dV)) ;
   if ( IsNull( pCrvCo))
      return false ;
  // ciclo sulle curve componenti (tutte curve di Bezier)
   int i = 0 ;
   const ICurve* pSCrv = pCrvCo->GetFirstCurve() ;
   while ( pSCrv != nullptr) {
      const CurveBezier* pCbz = GetBasicCurveBezier( pSCrv) ;
      if ( pCbz == NULL)
         return false ;
      int nCurrSpanStep = nStep / m_nSpanU ;
      if ( nCurrSpanStep <= 0) {
         double dLenU = 0 ;
         pCbz->GetApproxLength( dLenU) ;
         nCurrSpanStep = GetSteps( m_nDegU, 1, dLenU, 2) ;
      }
      PolyLine plCurr ;
      if ( ! pCbz->ApproxWithLines( nCurrSpanStep, plCurr))
         return false ;
      plCrvU.Join( plCurr, i) ;
      ++ i ;
      pSCrv = pCrvCo->GetNextCurve() ;
   }
   return true ;
}

//----------------------------------------------------------------------------
bool
SurfBezier::GetCurveOnV( double dU, int nStep, PolyLine& plCrvV) const
{
  // controlli
   plCrvV.Clear() ;
   if ( dU < - EPS_PARAM  ||  dU > m_nSpanU + EPS_PARAM)
      return false ;
   dU = Clamp( dU, 0., double( m_nSpanU)) ;
  // recupero la curva
   PtrOwner<CurveComposite> pCrvCo( GetCurveOnV( dU)) ;
   if ( IsNull( pCrvCo))
      return false ;
  // ciclo sulle curve componenti (tutte curve di Bezier)
   int i = 0 ;
   const ICurve* pSCrv = pCrvCo->GetFirstCurve() ;
   while ( pSCrv != nullptr) {
      const CurveBezier* pCbz = GetBasicCurveBezier( pSCrv) ;
      if ( pCbz == NULL)
         return false ;
      int nCurrSpanStep = nStep / m_nSpanU ;
      if ( nCurrSpanStep <= 0) {
         double dLenV = 0 ;
         pCbz->GetApproxLength( dLenV) ;
         nCurrSpanStep = GetSteps( m_nDegV, 1, dLenV, 2) ;
      }
      PolyLine plCurr ;
      if ( ! pCbz->ApproxWithLines( nCurrSpanStep, plCurr))
         return false ;
      plCrvV.Join( plCurr, i) ;
      ++ i ;
      pSCrv = pCrvCo->GetNextCurve() ;
   }
   return true ;
}

//----------------------------------------------------------------------------
bool
SurfBezier::GetControlCurveOnU( int nIndV, PolyLine& plCtrlU) const
{
   plCtrlU.Clear() ;
   if ( nIndV < 0  ||  nIndV > m_nDegV * m_nSpanV)
      return false ;
   for ( int i = 0 ; i <= m_nDegU * m_nSpanU ; ++ i)
      plCtrlU.AddUPoint( double( i) / m_nDegU, m_vPtCtrl[GetInd( i, nIndV)]) ;

   return true ;
}

//----------------------------------------------------------------------------
bool
SurfBezier::GetControlCurveOnV( int nIndU, PolyLine& plCtrlV) const
{
   plCtrlV.Clear() ;
   if ( nIndU < 0  ||  nIndU > m_nDegU * m_nSpanU)
      return false ;
   for ( int j = 0 ; j <= m_nDegV * m_nSpanV ; ++ j)
      plCtrlV.AddUPoint( double( j) / m_nDegV, m_vPtCtrl[GetInd( nIndU, j)]) ;

   return true ;
}

//----------------------------------------------------------------------------
double
SurfBezier::GetCurveOnUApproxLen( double dV) const
{
   PtrOwner<CurveComposite> pCrvCo( GetCurveOnU( dV)) ;
   double dLen ;
   if ( IsNull( pCrvCo)  ||  ! pCrvCo->GetApproxLength( dLen))
      return 0 ;
   return dLen ;
}

//----------------------------------------------------------------------------
double
SurfBezier::GetCurveOnVApproxLen( double dU) const
{
   PtrOwner<CurveComposite> pCrvCo( GetCurveOnV( dU)) ;
   double dLen ;
   if ( IsNull( pCrvCo)  ||  ! pCrvCo->GetApproxLength( dLen))
      return 0 ;
   return dLen ;
}

//----------------------------------------------------------------------------
bool
SurfBezier::UpdateEdgesFromTree( Tree& tr) const
{
   POLYLINEMATRIX mPlEdges ;
   tr.GetEdges3D( mPlEdges) ;
   for ( int i= 0 ; i < int( mPlEdges.size()); ++i) {
      for ( int j = 0 ; j < int ( mPlEdges[i].size()) ; ++j) {
         m_mCCEdge[i].emplace_back(CreateBasicCurveComposite()) ;
         if ( ! m_mCCEdge[i].back()->FromPolyLine(mPlEdges[i][j]) ) {
            Point3d ptStart ; 
            if ( ! mPlEdges[i][j].GetFirstPoint( ptStart))
               continue ;
            m_mCCEdge[i].back()->FromPoint( ptStart) ;
         }
      }
   }

   if ( m_bTrimmed) {
      POLYLINEVECTOR vPl ;
      tr.GetSplitLoops( vPl) ;

      // recupero i loop nel parametrico
      for( int i = 0 ; i < int( vPl.size()); ++i) {
         m_vCCLoop.emplace_back(CreateBasicCurveComposite()) ;
         m_vCCLoop.back()->FromPolyLine(vPl[i]) ;
      }
   }

   return true ;
}

//----------------------------------------------------------------------------
const SurfTriMesh*
SurfBezier::GetAuxSurf( void) const
{
  // la superficie deve essere validata
   if ( m_nStatus != OK) {
      ResetAuxSurf() ;
      return nullptr ;
   }
  // se già calcolata, la restituisco
   if ( m_pSTM != nullptr)
      return m_pSTM ;
  // eseguo calcolo
   m_pSTM = GetApproxSurf( 50 * EPS_SMALL, 100 * EPS_SMALL) ;
   return m_pSTM ;
}

//----------------------------------------------------------------------------
SurfTriMesh*
SurfBezier::GetApproxSurf( double dTol, double dSideMin) const
{
  // la superficie deve essere validata
   if ( m_nStatus != OK)
      return nullptr ;

  // se c'è ausiliaria e richiesta con gli stessi parametri, ne restituisco una copia
   if ( m_pSTM != nullptr  &&
        abs( dTol - 50 * EPS_SMALL) < EPS_SMALL  &&  abs( dSideMin - 100 * EPS_SMALL) < EPS_SMALL)
      return m_pSTM->Clone() ;
  
  // costruttore della superficie
   POLYLINEMATRIX vvPL ;
   //POLYLINEVECTOR vPL ;    // per usare i polygon basic
   //Tree Tree( this, true) ;
   Tree Tree ;
   if ( ! Tree.SetSurf( this, true))
      return nullptr ;
   BIPNTVECTOR vTrees ;
   Tree.GetIndependentTrees( vTrees) ;
   bool bTest = false ; // per debug
   // resetto il vettore degli edge
   m_mCCEdge.clear() ;
   m_mCCEdge = vector<ICRVCOMPOPOVECTOR>(4) ;
   m_vCCLoop.clear() ;
   for ( int i = 0 ; i < (int) vTrees.size() ; ++ i) {
      Point3d ptMin = get<0>( vTrees[i]) ;
      Point3d ptMax = get<1>( vTrees[i]) ;
      Tree.SetSurf( this, true, ptMin, ptMax) ;
      if ( bTest) {
         Tree.BuildTree_test() ;        // per debug
         //Tree.BuildTree( 5 * LIN_TOL_FINE, 1) ; // per debug
         Tree.SetTestMode() ;
      }
      else {
         //Tree.BuildTree( 5 * LIN_TOL_FINE, 0.1) ;
         Tree.BuildTree( dTol, dSideMin) ;
      }
      if ( ! Tree.GetPolygons( vvPL))
         continue ;
      //Tree.GetPolygonsBasic( vPL) ;  // per usare i polygon basic

     // aggiorno la chiusura della superficie
      m_bClosedU = m_bClosedU  ||  Tree.IsClosedU() ;
      m_bClosedV = m_bClosedV  ||  Tree.IsClosedV() ;
     // salvo i bordi in 3d, che servono in caso si voglia trimmare la superficie DOPO aver costruito la trimesh ausiliaria
      UpdateEdgesFromTree( Tree) ;
   }

  //// per usare i polygon basic//////////////////////
   //for (int k = 0 ; k < (int)vPL.size(); ++k) {
   //   vvPL.emplace_back() ;
   //   vvPL.back().push_back(vPL[k]) ;
   //}
  //// per usare i polygon basic///////////////////

   if ( vvPL.empty())
      LOG_DBG_ERR( GetEGkLogger(), "ERROR : Bezier Surface couldn't be triangulated") ;
   
   StmFromTriangleSoup stmSoup ;
   if ( ! stmSoup.Start())
      return nullptr ;

  // prendo i punti di ogni polyline dell'albero, li triangolo e li porto in 3d
   for ( POLYLINEVECTOR vPL : vvPL) {
      PNTVECTOR vPnt ;
      INTVECTOR vTria ;
      Triangulate Tri ;
      if ( ! Tri.Make( vPL, vPnt, vTria))
         return nullptr ;
     
     // porto i punti in 3d
      PNTVECTOR vPnt3d ;
      for ( int i = 0 ; i < int( vPnt.size()) ; ++ i) {
         Point3d pt3d ;
         if ( ! GetPointD1D2( vPnt[i].x / SBZ_TREG_COEFF, vPnt[i].y / SBZ_TREG_COEFF, ISurfBezier::FROM_MINUS, ISurfBezier::FROM_MINUS, pt3d))
            return nullptr ;
         vPnt3d.push_back( pt3d) ;
      }
      int nTria = int( vTria.size()) / 3 ;
      for ( int i = 0 ; i < nTria ; ++i) {
         if ( ! stmSoup.AddTriangle( vPnt3d[vTria[3*i]], vPnt3d[vTria[3*i+1]], vPnt3d[vTria[3*i+2]],
                                     vPnt[vTria[3*i]].x, vPnt[vTria[3*i]].y,
                                     vPnt[vTria[3*i+1]].x, vPnt[vTria[3*i+1]].y,
                                     vPnt[vTria[3*i+2]].x, vPnt[vTria[3*i+2]].y))
            return nullptr ;
      }
   } 

  // termino
   if ( ! stmSoup.End())
      return nullptr ;
  // restituisco
   return GetBasicSurfTriMesh( stmSoup.GetSurf()) ;
}

//----------------------------------------------------------------------------
bool
SurfBezier::GetLeaves( vector<tuple<int, Point3d, Point3d>>& vLeaves) const
{
   Tree Tree ;
   if ( ! Tree.SetSurf( this, true))
      return false ;
   BIPNTVECTOR vTrees ;
   Tree.GetIndependentTrees( vTrees) ;
   for ( int i = 0 ; i < int( vTrees.size()) ; ++ i) {
      Point3d ptMin = get<0>( vTrees[i]) ;
      Point3d ptMax = get<1>( vTrees[i]) ;
      Tree.SetSurf( this, true, ptMin, ptMax) ;
      bool bTest = false ; // per debug
      if ( bTest) {
         Tree.BuildTree_test() ;                    // per debug
         //Tree.BuildTree( 5 * LIN_TOL_FINE, 1) ;     // per debug
      }
      else {
         Tree.BuildTree( 5 * LIN_TOL_FINE, 0.1) ;
      }
      vector<Cell> vCells ;
      Tree.GetLeaves( vCells) ;
      for ( int k = 0 ; k < int( vCells.size()) ; ++ k) {
         vLeaves.emplace_back( vCells[k].m_nId, vCells[k].GetBottomLeft(), vCells[k].GetTopRight()) ;
      }
   }
   return true ;
}

//----------------------------------------------------------------------------
bool
SurfBezier::GetTriangles2D( vector<tuple<int,Point3d, Point3d, Point3d>>& vTria2D) const
{
   const ISurfTriMesh* pSTM = GetAuxSurf() ;
   for ( int t = 0 ; t < int(pSTM->GetTriangleCount()) ; ++t ) {
      double dU0, dU1, dU2, dV0, dV1, dV2 ;
      int nVert[3] ;
      pSTM->GetTriangle( t, nVert);
      pSTM->GetVertexParam( nVert[0], dU0, dV0) ;
      pSTM->GetVertexParam( nVert[1], dU1, dV1) ;
      pSTM->GetVertexParam( nVert[2], dU2, dV2) ;
      Point3d pt0(dU0,dV0), pt1(dU1,dV1), pt2(dU2,dV2) ;
      vTria2D.emplace_back( tuple<int,Point3d, Point3d, Point3d>(t, pt0, pt1, pt2)) ;
   }
   return true ;
}

//----------------------------------------------------------------------------
void
SurfBezier::ResetAuxSurf( void) const
{
   if ( m_pSTM != nullptr)
      delete( m_pSTM) ;
   m_pSTM = nullptr ;
}

//----------------------------------------------------------------------------
void
SurfBezier::ResetTrimRegion( void)
{
   if ( m_pTrimReg != nullptr)
      delete( m_pTrimReg) ;
   m_pTrimReg = nullptr ;
   ResetAuxSurf() ;
   m_bTrimmed = false ;
   // imposto ricalcolo della grafica
   m_OGrMgr.Reset() ;
}

//----------------------------------------------------------------------------
bool
SurfBezier::IncreaseUV( Point3d& ptUV, Vector3d vtH , Point3d* ptUVCopy, bool bModifyOrig) const
{
   if ( ptUVCopy != nullptr) {
      IncreaseUV( ptUV.x, vtH.x, true, &(*ptUVCopy).x, bModifyOrig) ;
      IncreaseUV( ptUV.y, vtH.y, true, &(*ptUVCopy).y, bModifyOrig) ;
   }
   else {
      IncreaseUV( ptUV.x, vtH.x, true, nullptr, bModifyOrig) ;
      IncreaseUV( ptUV.y, vtH.y, true, nullptr, bModifyOrig) ;
   }
   
   return true ;
}

//----------------------------------------------------------------------------
bool
SurfBezier::IncreaseUV( double& dUV, double dxy, bool bUOrV, double* dUVCopy, bool bModifyOrig) const
{
   double dUVTest ;
   if ( dUVCopy != nullptr) {
      *dUVCopy = dUV + dxy ;
      dUVTest = *dUVCopy ;
   }
   if ( bModifyOrig) {
      dUV += dxy ;
      dUVTest = dUV ;
   }
   
   if ( bUOrV) {
      if ( dUVTest < 0)
            dUVTest = 0 ;
      else if ( dUVTest > m_nSpanU * SBZ_TREG_COEFF )
            dUVTest = m_nSpanU * SBZ_TREG_COEFF ;
   }
   else {
      if ( dUVTest < 0)
            dUVTest = 0 ;
      else if ( dUVTest > m_nSpanV * SBZ_TREG_COEFF )
            dUVTest = m_nSpanV * SBZ_TREG_COEFF ;
   }
   if ( bModifyOrig)
      dUV  = dUVTest ;
   if ( dUVCopy != nullptr)
      *dUVCopy = dUVTest ;

   return true ;
}

//----------------------------------------------------------------------------
bool
SurfBezier::UnprojectCurveFromStm( const ICurveComposite* pCC, ICRVCOMPOPVECTOR& vpCC, const Plane3d* pPlCut) const 
{
   // do per scontato che la compo sia una spezzata, visto che arriva dall'intersezione tra un piano e una trimesh
   // creo la chain dei punti che sto riportando nel parametrico
   ChainCurves chainC ;
   double dToler = EPS_SMALL ;
   chainC.Init( false, dToler, 2) ;

   const ICurve* pCrv0 = pCC->GetCurve( 0) ;
   PolyLine pl ;
   Point3d pt3D, pt2D ; pCrv0->GetStartPoint( pt3D) ;
   Point3d pt3DEnd ; pCrv0->GetEndPoint( pt3DEnd) ;
   bool bThroughEdge = false ;
   BOOLVECTOR vbThroughEdge ;
   if ( ! UnprojectPoint(pt3D, pt2D, pt3DEnd, &bThroughEdge, pPlCut))
      return false ;
   vbThroughEdge.push_back( bThroughEdge) ;

   // aggiungo tutti i successivi
   BIPNTVECTOR vBPnt ;

   bThroughEdge = false ;
   int nRejected = 0 ;
   for ( int i = 0 ; i < int( pCC->GetCurveCount()) ; ++i) {
      const ICurve* pCrv = pCC->GetCurve( i) ;
      Point3d pt3DPrev = pt3D ;
      Point3d pt2DPrev = pt2D ;
      bool bPrevIsPole = false ;
      if ( bThroughEdge) {
         // devo cambiare le coordinate di pt2DPrev per periodicità
         // capisco su quale lato è e lo porto sul lato opposto
         Point3d pt = pt2DPrev ;
         if ( m_bClosedU) {
            if ( pt2DPrev.x < 1)
               pt2DPrev.x = m_nSpanU * SBZ_TREG_COEFF ;
            else if ( (m_nSpanU * SBZ_TREG_COEFF - pt2DPrev.x) < 1)
               pt2DPrev.x = 0 ;
         }
         if ( m_bClosedV) {
            if ( pt2DPrev.y < 1)
               pt2DPrev.y = m_nSpanV * SBZ_TREG_COEFF ;
            else if ( (m_nSpanV * SBZ_TREG_COEFF - pt2DPrev.y) < 1)
               pt2DPrev.y = 0 ;
         }
         bPrevIsPole = AreSamePointApprox( pt, pt2DPrev) ;
      }
      pCrv->GetEndPoint( pt3D) ;
      if ( ! UnprojectPoint( pt3D, pt2D, pt3DPrev, &bThroughEdge, pPlCut))
         return false ;
      if ( bPrevIsPole) {
         // se il punto precedente era di polo allora devo correggere le sue coordinate 2D
         if ( ( m_vbPole[0] || m_vbPole[2] )) {
            if ( pt2DPrev.y < 1  ||  m_nSpanV * SBZ_TREG_COEFF - pt2DPrev.y < 1)
               pt2DPrev.x = pt2D.x ;
         }
         if ( ( m_vbPole[1] || m_vbPole[3] )) {
            if ( pt2DPrev.x < 1  ||  m_nSpanU * SBZ_TREG_COEFF - pt2DPrev.x < 1)
               pt2DPrev.y = pt2D.y ;
         } 
      }
      Vector3d vtDir = pt2D - pt2DPrev ;
      vtDir.Normalize() ;
      // se mi accorgo che sto per tracciare un taglio lungo un bordo posso semplicmente evitarlo
      if ( (1 - abs(vtDir.x) < EPS_SMALL  &&  (pt2D.y < EPS_SMALL  ||  m_nSpanV * SBZ_TREG_COEFF - pt2D.y < 1))  ||  // parallelo agli edge 0 e 2 e su uno di questi
            (1 - abs(vtDir.y) < EPS_SMALL  &&  (pt2D.x < EPS_SMALL  ||  m_nSpanU * SBZ_TREG_COEFF - pt2D.x < 1))) {  // parallello agli edge 1 e 3 e su uno di questi
         ++ nRejected ;
         continue ;
      }
      if ( bThroughEdge && vbThroughEdge.back()) {
         double dParamH, dParamL ;
         dParamH = m_nSpanV * SBZ_TREG_COEFF ;
         dParamL = m_nSpanU * SBZ_TREG_COEFF ;
         // sia questo punto che il precedente sono su un edge, ma il segmento che li unisce non è parallelo ad un edge
         // potrei star tracciando un taglio sul bordo di chiusura
         // controllo se sto tracciando una linea che unisce due lati di chiusura, allora in realtà sdtarei tracciando un taglio sull'edge e quindi posso non tracciarlo
         if ( (abs( vtDir.x) > abs( vtDir.y)  &&  Dist( pt2D, pt2DPrev) > dParamL * 0.5)  ||
              (abs( vtDir.y) > abs( vtDir.x)  &&  Dist( pt2D, pt2DPrev) > dParamH * 0.5)) {
            ++ nRejected ;
            continue ;
         }
      }
      vbThroughEdge.push_back( bThroughEdge) ;
      vBPnt.emplace_back( BIPOINT( pt2DPrev, pt2D)) ;
      if ( ! chainC.AddCurve( i + 1 - nRejected, pt2DPrev, vtDir, pt2D, vtDir))
         return false ;
   }
   // ricostruisco le catene in 2D
   Point3d ptNear ; pCrv0->GetStartPoint( ptNear) ;
   INTVECTOR vId ;
   bool bAdded = true ;
   while( chainC.GetChainFromNear( pt3D, false, vId)) {
      PtrOwner<ICurveComposite> pCC2D ( CreateCurveComposite()) ;
      for ( int i = 0 ; i < int( vId.size()) ; ++ i) {
         // creo un segmento di retta
         ICurveLine* pLine( CreateCurveLine()) ;
         if ( pLine == nullptr)
            return false ;
         // recupero gli estremi (non vanno mai invertiti per opzione di concatenamento)
         int nInd = abs( vId[i]) - 1 ;
         Point3d ptStart = ( bAdded ? vBPnt[nInd].first : ptNear) ;
         Point3d ptEnd = vBPnt[nInd].second ;
         // provo ad accodarlo alla composita
         bAdded = ( Dist( ptStart, ptEnd) > dToler / 2  &&
            pLine->Set( ptStart, ptEnd)) ;
         bAdded = bAdded && pCC2D->AddCurve( pLine, true, dToler) ;
         ptNear = ( bAdded ? ptEnd : ptStart) ;
      }
      if ( pCC2D->IsValid())
         vpCC.emplace_back( Release(pCC2D)) ;
   }
   return true ;
}

//----------------------------------------------------------------------------
bool
SurfBezier::AddCurveCompoToCuts( ICurveComposite* pCrvCompo, ICRVCOMPOPOVECTOR& vpCCOpen, ICRVCOMPOPOVECTOR& vpCCClosed, double dToler, const Plane3d* pPlCut) const
{
   // se lunghezza curva inferiore a 5 volte la tolleranza, la ignoro e sposto il punto finale nel punto di fine della curva che sto ignorando
   double dCrvLen ;
   if ( ! pCrvCompo->GetLength( dCrvLen)  ||  dCrvLen < 5. * dToler)
      return true ;
   // se curva chiusa entro 5 volte la tolleranza ma considerata aperta, la chiudo bene
   Point3d ptStart, ptEnd ;
   if ( pCrvCompo->GetStartPoint( ptStart)  &&
      pCrvCompo->GetEndPoint( ptEnd)  &&
      AreSamePointEpsilon( ptStart, ptEnd, 5. * dToler)  &&
      ! AreSamePointApprox( ptStart, ptEnd)) {
      // porto il punto finale a coincidere esattamente con l'inizio
      pCrvCompo->ModifyEnd( ptStart) ;
   }
   // unisco segmenti allineati
   pCrvCompo->MergeCurves( 0.5 * dToler, ANG_TOL_STD_DEG) ;

   // porto la curva nello spazio parametrico
   ICRVCOMPOPVECTOR vCC ;
   if ( ! UnprojectCurveFromStm( pCrvCompo, vCC, pPlCut))
      return false ;
   for ( int i = 0 ; i < int( vCC.size()); ++i) {
      // le curve aperte le tengo da parte per giuntarle alla fine col bordo
      if ( ! vCC[i]->IsClosed() )
         vpCCOpen.emplace_back( vCC[i]) ;
      // le curve chiuse le metto tutte insieme subito
      else
         vpCCClosed.emplace_back( vCC[i]) ;
   }
   return true ;
}

//----------------------------------------------------------------------------
typedef tuple<int,int,int> TRINT ;

template<>
struct hash<TRINT> {
   size_t operator()(const TRINT& t) const
   {
      // Compute individual hash values for first, second and third and combine them using XOR and bit shifting:
      return ((hash<int>()(get<0>(t))) ^ (hash<int>()(get<1>(t)) << 1) >> 1) ^ (hash<int>()(get<2>(t)) << 1) ;
   }
};

//----------------------------------------------------------------------------
bool
SurfBezier::Cut( const Plane3d& plPlane, bool bSaveOnEq)
{
   // faccio l'intersezione della trimesh ausiliaria con il piano posso ottenere: punti, curve 3d e triangoli( coplanari al piano di taglio)
   // i punti li escludo
   // le curve 3d le trasformo in curve 2d e le aggiungo alle curve di trim
   // accorpo eventuali triangoli adiacenti ed estraggo i loop delle regioni ottenute; questi vengono poi portati in 2d e aggiunti alle curve di trim

   // se necessario calcolo i poli
   if ( m_vbPole.empty())
      CalcPoles() ;

   PNTVECTOR vPnt ;
   BIPNTVECTOR vBPnt ;
   TRIA3DVECTOR vTria ;
   IntersPlaneSurfTm( plPlane, *GetAuxSurf(), vPnt, vBPnt, vTria) ;

   // concateno le curve 3d
   ChainCurves chainC ;
   double dToler = EPS_SMALL ;
   chainC.Init( false, dToler, int( vBPnt.size())) ;
   for ( int i = 0 ; i < int( vBPnt.size()) ; ++ i) {
      Vector3d vtDir = vBPnt[i].second - vBPnt[i].first ;
      vtDir.Normalize() ;
      if ( ! chainC.AddCurve( i + 1, vBPnt[i].first, vtDir, vBPnt[i].second, vtDir))
         return false ;
   }
  // GESTIONE DELLE CURVE OTTENUTE DALL'INTERSEZIONE

   // recupero i percorsi concatenati
   Point3d ptNear = ( vBPnt.empty() ? ORIG : vBPnt[0].first) ;
   INTVECTOR vId ;

   // separo tra loop chiusi, interni allo spazio parametrico e loop passanti che tagliano lo spazio intersecando i bordi
   ICRVCOMPOPOVECTOR vpCCOpen ;
   ICRVCOMPOPOVECTOR vpCCClosed ;

   while ( chainC.GetChainFromNear( ptNear, false, vId)) {
      // creo una curva composita
      PtrOwner<ICurveComposite> pCrvCompo( CreateCurveComposite()) ;
      if ( IsNull( pCrvCompo))
         return false ;
      // recupero gli estremi dei segmenti, creo le linee e le inserisco nella composita
      bool bAdded = true ;
      for ( int i = 0 ; i < int( vId.size()) ; ++ i) {
         // creo un segmento di retta
         ICurveLine* pLine( CreateCurveLine()) ;
         if ( pLine == nullptr)
            return false ;
         // recupero gli estremi (non vanno mai invertiti per opzione di concatenamento)
         int nInd = abs( vId[i]) - 1 ;
         Point3d ptStart = ( bAdded ? vBPnt[nInd].first : ptNear) ;
         Point3d ptEnd = vBPnt[nInd].second ;
         // provo ad accodarlo alla composita
         bAdded = ( Dist( ptStart, ptEnd) > dToler / 2  &&
            pLine->Set( ptStart, ptEnd)) ;
         bAdded = bAdded && pCrvCompo->AddCurve( pLine, true, dToler) ;
         ptNear = ( bAdded ? ptEnd : ptStart) ;
      }
      if ( ! AddCurveCompoToCuts( pCrvCompo, vpCCOpen, vpCCClosed, EPS_SMALL, &plPlane))
         return false ;
   }

   //comincio a creare la superficie aggiungendo i tagli aperti ai bordi attualmente esistenti
   SurfFlatRegionByContours sfrContour ;
   if ( int(vpCCOpen.size()) != 0 ) {
      // qui devo aggiungere tutto del codice nuovo per ricostruire in altro modo il nuovo bordo della superficie
      // recupero la regione attuale
      
      
      PtrOwner<ISurfFlatRegion> pNewTrim( CreateBasicSurfFlatRegion()) ;
      if ( m_bTrimmed)
         pNewTrim.Set( GetTrimRegion()->Clone()) ;
      else
         pNewTrim.Set( GetSurfFlatRegionRectangle( SBZ_TREG_COEFF * m_nSpanU, SBZ_TREG_COEFF * m_nSpanV)) ;

      // costruisco la mappa delle intersezioni, trovando tutte le intersezioni tra i trim e i loop dei vari chunk della falr region
      unordered_map<TRINT,ICCIVECTOR> mInters ;
      int nInters = 0 ;
      bool bStartFound = false ;
      bool bEndFound = false ;
      // trim
      for ( int t = 0 ; t < int( vpCCOpen.size()); ++t) {
         nInters = 0 ;
         bStartFound = false ;
         bEndFound = false ;
         //chunk
         for ( int c = 0 ; c < pNewTrim->GetChunkCount() ; ++c) {
            // loop
            for ( int l = 0 ; l < pNewTrim->GetLoopCount( c) ; ++l) {
               PtrOwner<ICurve> pLoop( pNewTrim->GetLoop( c, l)) ;
               // prima curva è il loop, seconda curva è il trim
               IntersCurveCurve icc( *pLoop, *vpCCOpen[t]) ;
               if ( icc.GetIntersCount() != 0) {
                  ICCIVECTOR vICC ;
                  for ( int i = 0 ; i < int( icc.GetIntersCount()); ++i) {
                     IntCrvCrvInfo iccInfo ;
                     icc.GetIntCrvCrvInfo( i, iccInfo) ;
                     vICC.emplace_back( iccInfo) ;
                  }
                  mInters.insert( pair<TRINT,ICCIVECTOR>( TRINT(c,l,t), vICC)) ;
                  if ( int( vICC.size() == 2)) {
                     bStartFound = true ;
                     bEndFound = true ;
                  }
                  else if ( int(vICC.size() == 1) ) {
                     if ( vICC[0].IciB->dU < EPS_SMALL)
                        bStartFound = true ;
                     else
                        bEndFound = true ;
                  }
               }
               nInters += int( icc.GetIntersCount()) ;
            }
         }
         if ( nInters != 2) {
            // se un trim non fa 2 intersezioni allora devo estendere la curva allo start e/o all'end per creare le intersezioni
            Point3d ptStart ; vpCCOpen[t]->GetStartPoint( ptStart) ;
            Point3d ptEnd ; vpCCOpen[t]->GetEndPoint( ptEnd) ;
            PtrOwner<ICurve> pCrv( vpCCOpen[t]->Clone()) ;
            double dExtension = m_nSpanU > m_nSpanV ? m_nSpanU : m_nSpanV ;
            dExtension *= SBZ_TREG_COEFF ;
            if ( ! bStartFound)
               pCrv->ExtendStartByLen( dExtension) ;
            if ( ! bEndFound)
               pCrv->ExtendEndByLen( dExtension) ;
            double dDistStart = 1e6, dDistEnd = 1e6 ;
            // vettore per l'intersezione di start
            ICCIVECTOR vICCStart ;
            vICCStart.emplace_back() ;
            // vettore per l'intersezione di end
            ICCIVECTOR vICCEnd ;
            vICCEnd.emplace_back() ;
            TRINT tStart, tEnd ;
            //chunk
            for ( int c = 0 ; c < pNewTrim->GetChunkCount() ; ++c) {
               // loop
               for ( int l = 0 ; l < pNewTrim->GetLoopCount( c) ; ++l) {
                  PtrOwner<ICurve> pLoop( pNewTrim->GetLoop( c, l)) ;
                  // prima curva è il loop, seconda curva è il trim
                  IntersCurveCurve icc( *pLoop, *pCrv) ;
                  if ( icc.GetIntersCount() != 0) {
                     for ( int i = 0 ; i < int( icc.GetIntersCount()); ++i) {
                        IntCrvCrvInfo iccInfo ; icc.GetIntCrvCrvInfo( i, iccInfo) ;
                        if ( ! bStartFound  &&  Dist( iccInfo.IciA->ptI, ptStart) < dDistStart) {
                           dDistStart = Dist( iccInfo.IciA->ptI, ptStart) ;
                           vICCStart[0] = iccInfo ;
                           tStart = TRINT( c, l, t) ;
                        }
                        if ( ! bEndFound  &&  Dist( iccInfo.IciA->ptI, ptEnd) < dDistEnd) {
                           dDistEnd = Dist( iccInfo.IciA->ptI, ptEnd) ;
                           vICCEnd[0] = iccInfo ;
                           tEnd = TRINT( c, l, t) ;
                        }
                     }
                  }
               }
            }

            // ricostruisco gli elementi per la mappa mInters
            if ( ! bStartFound  &&  ! bEndFound) {
               // ridefinisco il taglio aperto con la sua versione estesa che arriva a toccare i loop dello spazio parametrico
               vpCCOpen[t].Set( GetCurveComposite(pCrv->CopyParamRange( vICCStart[0].IciB->dU, vICCEnd[0].IciB->dU))) ;
               // correggo il parametro dell'intersezione allo start
               vICCStart[0].IciB->dU = 0 ;
               if ( tStart == tEnd) {
                  // se ho intersezione con un loop solo allora accorpo i due vettori delle intersezioni
                  vICCStart.emplace_back( vICCEnd[0]) ;
                  mInters.insert( pair<TRINT,ICCIVECTOR>(tStart, vICCStart)) ;
               }
               else{
                  // se ho intersezione con due loop diverse due entry diverse le inserisco
                  mInters.insert( pair<TRINT,ICCIVECTOR>(tStart, vICCStart)) ;
                  mInters.insert( pair<TRINT,ICCIVECTOR>(tEnd, vICCEnd)) ;
               }
            }
            else {
               // devo verificare se avevo già trovato una delle due intersezioni e se era sullo stesso loop o no
               if ( ! bStartFound) {
                  pCrv->TrimStartAtParam( vICCStart[0].IciB->dU) ;
                  vpCCOpen[t].Set( GetCurveComposite( Release( pCrv))) ;
                  // correggo il parametro dell'intersezione allo start
                  vICCStart[0].IciB->dU = 0 ;
                  if ( mInters.count( tStart) == 1)
                     mInters[tStart].emplace_back( vICCStart[0]) ;
                  else
                     mInters.insert( pair<TRINT,ICCIVECTOR>( tStart, vICCStart)) ;
               }
               if ( ! bEndFound) {
                  if ( mInters.count( tEnd) == 1)
                     mInters[tEnd].emplace_back( vICCEnd[0]) ;
                  else
                     mInters.insert( pair<TRINT,ICCIVECTOR>( tEnd, vICCEnd)) ;
                  pCrv->TrimEndAtParam( vICCEnd[0].IciB->dU) ;
                  vpCCOpen[t].Set( GetCurveComposite( Release( pCrv))) ;
               }
            }
         }
      }

      // vettore di flag che mi indica quali tagli aperti sono stati aggiunti al nuovo bordo
      BOOLVECTOR vbAdded( vpCCOpen.size()) ;
      fill( vbAdded.begin(), vbAdded.end(), false) ;
      PtrOwner<ICurveComposite> pCCNewEdge( CreateCurveComposite()) ;
      PtrOwner<ICurveLine> pCL( CreateCurveLine()) ;
      TRINT tiFirstInters ;
      // parto aggiungendo il primo taglio
      int nNewToAdd = 0 ;
      bool bFirstCurveOfEdge = true ;
      while ( nNewToAdd != -1) {
         // aggiungo il  taglio
         pCCNewEdge->AddCurve( Release( vpCCOpen[nNewToAdd])) ;
         // aggiorno la lista degli aggiunti
         vbAdded[nNewToAdd] = true ;
         // di questo taglio mi salvo il chunk e loop di start e end
         TRINT tiStart, tiEnd ;
         for (const auto& pair : mInters) {
            if ( get<2>(pair.first) == nNewToAdd ) {
               for (int p = 0 ; p < int(pair.second.size()) ; ++p) {
                  if ( pair.second[p].IciB->dU < EPS_SMALL) {
                     tiStart = pair.first ;
                     if ( bFirstCurveOfEdge){
                        // salvo l'inizio del taglio che è la prima curva di questa curva compo
                        tiFirstInters = pair.first ;
                        bFirstCurveOfEdge = false ;
                     }
                  }
                  else
                     tiEnd = pair.first ;
               }
            }
         }
         
         // devo trovare fino a che punto seguire il loop che ho trovato come prosecuzione del taglio corrente
         // devo quindi trovare la prossima intersezione con un taglio
         int nInters = -1 ;
         double dNextCut = numeric_limits<double>::infinity() ;
         double dEndCurrentCut ;
         for ( int i = 0 ; i < int( mInters[tiEnd].size()); ++i) {
            // se ho trovato l'intersezione con la fine del taglio corrente, salvo il parametro sul loop
            if ( mInters[tiEnd][i].IciB->dU > EPS_SMALL)
               dEndCurrentCut = mInters[tiEnd][i].IciA->dU ;
         }
         // se non trovo nessuna altra intersezione prima della fine del loop allora devo ripetere tutto cercando a partire dall'inizio del loop
         for ( const auto& pair : mInters) {
            if ( get<0>(pair.first) == get<0>(tiEnd) && get<1>(pair.first) == get<1>(tiEnd)) {
               for ( int i = 0 ; i < int(pair.second.size()); ++i ) {
                  // se trovo una nuova intersezione che incontro prima di quella che mi ero salvato precedentemente allora
                  // mi salvo questa nuova che ho trovato
                  if ( pair.second[i].IciA->dU < dNextCut  &&  pair.second[i].IciA->dU > dEndCurrentCut) {
                     dNextCut = pair.second[i].IciA->dU ;
                     nInters = get<2>(pair.first) ;
                  }
               }
            }
         }
         if ( nInters == -1) {
            dNextCut = numeric_limits<double>::infinity() ;
            for ( const auto& pair : mInters) {
               if ( get<0>(pair.first) == get<0>(tiEnd) && get<1>(pair.first) == get<1>(tiEnd)) {
                  for ( int i = 0 ; i < int(pair.second.size()); ++i ) {
                     // se trovo una nuova intersezione che incontro prima di quella che mi ero salvato precedentemente allora
                     // mi salvo questa nuova che ho trovato
                     if ( pair.second[i].IciA->dU < dNextCut) {
                        dNextCut = pair.second[i].IciA->dU ;
                        nInters = get<2>(pair.first) ;
                     }
                  }
               }
            }
         }

         // se tutto va bene queste due righe sostituiscono tutto il casino qua sotto
         PtrOwner<ICurve> pLoopTrimmed( pNewTrim->GetLoop( get<0>(tiEnd), get<1>(tiEnd))) ;
         pCCNewEdge->AddCurve(pLoopTrimmed->CopyParamRange( dEndCurrentCut, dNextCut)) ;

         // se il prossimo taglio identificato è quello da cui sono partito allora aggiungo il bordo ricostruito fino a questo momento
         // alla flat region e comincio a costruire un altro bordo
         // altrimenti continuo ad aggiungere curve al bordo corrente
         if ( nInters == get<2>(tiFirstInters) ) {
            pCCNewEdge->Close() ;
            sfrContour.AddCurve( Release( pCCNewEdge)) ;
            pCCNewEdge.Set( CreateBasicCurveComposite()) ;
            bFirstCurveOfEdge = true ;
            // trovo il prossimo taglio ancora da aggiungere
            nNewToAdd = -1 ;
            for ( int b = 0 ; b < int(vbAdded.size()) ; ++b ) {
               if ( ! vbAdded[b]) {
                  nNewToAdd = b ;
                  break ;
               }
            }
         }
         else
            nNewToAdd = nInters;
      }
   }

  //GESTIONE DEI TRIANGOLI RISULTANTI DALL'INTERSEZIONE
   StmFromTriangleSoup StmFts ;
   if ( ! StmFts.Start())
      return GDB_ID_NULL ;
   for ( int i = 0 ; i < int( vTria.size()) ; ++ i)
      // inserisco il triangolo nella nuova superficie
      StmFts.AddTriangle( vTria[i]) ;
   // valido la superficie e calcolo le adiacenze
   if ( ! StmFts.End())
      return GDB_ID_NULL ;
   // se superficie con triangoli
   PtrOwner<ISurfTriMesh> pNewStm( StmFts.GetSurf()) ;
   POLYLINEVECTOR vPLTria ; 
   if ( ! IsNull( pNewStm)  &&  pNewStm->GetTriangleCount() > 0) {
      pNewStm->GetLoops( vPLTria) ;
   }

   // aggiungo i loop chiusi
   for ( int i = 0 ; i < int( vpCCClosed.size()); ++i )
      sfrContour.AddCurve( Release( vpCCClosed[i])) ;
   // aggiungo loop derivati dai triangoli
   for ( int i = 0 ; i < int( vPLTria.size()); ++i ) {
      PtrOwner<ICurveComposite> pCC( CreateCurveComposite()) ;
      pCC->FromPolyLine( vPLTria[i]) ;
      sfrContour.AddCurve( Release( pCC)) ;
   }
   PtrOwner<ISurfFlatRegion> pSFR( sfrContour.GetSurf()) ;
   if ( IsNull( pSFR)  ||  ! pSFR->IsValid())
      return false ;

   // se la superficie ha normale con z negativa la inverto
   if ( pSFR->GetNormVersor().z < 0)
      pSFR->Invert() ;

   // verifico se la superficie che ho ottenuto è corretta o devo prendere il complementare ( rispetto allo spazio parametrico totale)
   // per verificarlo prendo un punto su questa superficie e verifico dove sta il suo corrispettivo 3D rispetto al piano di taglio
   int nChunkMin = 0 , nLoopMin = 0 ;
   // sono nello spazio parametrico, quindi le aree delle curve possono essere molto grandi
   double dAreaMin = 1e30 ;
   bool bPos = false ;
   for ( int c = 0 ; c < int( pSFR->GetChunkCount()); ++c) {
      PtrOwner<ISurfFlatRegion> pSurf( pSFR->CloneChunk( c)) ;
         for ( int l = 0 ; l < pSurf->GetLoopCount( 0); ++l) {
            PtrOwner<ICurve> pCrv( pSurf->GetLoop( 0, l)) ;
            double dArea ; pCrv->GetAreaXY( dArea) ;
            if ( abs( dArea) < dAreaMin) {
               nChunkMin = c ;
               nLoopMin = l ;
               dAreaMin = abs( dArea) ;
               bPos = dArea > 0 ;
            }
         }
   }

   // aggiorno la superficie di trim
   Point3d ptStart ;
   Vector3d vtDir, vtDirS, vtDirE ;
   PtrOwner<ICurve> pCrv( pSFR->GetLoop( nChunkMin, nLoopMin)) ;
   pCrv->GetStartPoint( ptStart) ;
   pCrv->GetStartDir( vtDirS) ;
   pCrv->GetEndDir( vtDirE) ;
   vtDir = vtDirS + vtDirE ;
   PtrOwner<ICurveLine> pCL( CreateCurveLine()) ;
   pCL->SetPVL( ptStart, vtDir, 1e6) ;
   IntersCurveCurve icc( *pCL, *pCrv) ;
   IntCrvCrvInfo iccInfo ;
   // verifico di guardare verso l'interno ( il numero di intersezioni deve essere pari visto che partivo da un punto sulla curva)
   if ( icc.GetIntersCount()%2 != 0) {
      vtDir = vtDirS - vtDirE ;
      PtrOwner<ICurveLine> pCL2( CreateCurveLine()) ; pCL2->SetPVL( ptStart, vtDir, 1e6) ;
      IntersCurveCurve icc2( *pCL2, *pCrv) ;
      if ( icc2.GetIntersCount()%2 != 0) {
         vtDir = vtDirS ;
         vtDir.Rotate( Z_AX, bPos? 90 : -90) ;
         PtrOwner<ICurveLine> pCL3( CreateCurveLine()) ; pCL3->SetPVL( ptStart, vtDir, 1e6) ;
         IntersCurveCurve icc3( *pCL3, *pCrv) ;
         if ( icc3.GetIntersCount()%2 != 0)
            return false ;
         icc3.GetIntCrvCrvInfo( 1, iccInfo) ;
      }
      else
         icc2.GetIntCrvCrvInfo( 1, iccInfo) ;
   }
   else
      icc.GetIntCrvCrvInfo( 1, iccInfo) ;

   Point3d ptI = iccInfo.IciA[0].ptI ;
   Point3d ptToCheck = ( ptStart + ptI) / 2 ;
   Point3d pt3D ; GetPointD1D2( ptToCheck.x / SBZ_TREG_COEFF, ptToCheck.y / SBZ_TREG_COEFF, ISurfBezier::FROM_MINUS, ISurfBezier::FROM_MINUS, pt3D) ;
   double dDist = DistPointPlane( pt3D, plPlane) ;

   // ho due casi in cui devo invertire(prendere il complementare rispetto allo spazio parametrico) la superficie:
   // 1. se il punto è sopra il piano ed era dentro una curva CCW
   // 2. se il punto è sotto il piano ed era interno ad una curva CW

   // la SetTrimRegion controlla se avevo trim precedenti ed eventualmente fa l'intersezione con lo spazio esistente
   if ( ( dDist > 0  &&  bPos) ||  ( dDist < 0  &&  ! bPos)) {
      if ( ! SetTrimRegion( *pSFR, false)  ||  ! m_pTrimReg->IsValid())
         return false ;
   }
   else {
      if ( ! SetTrimRegion( *pSFR)  ||  ! m_pTrimReg->IsValid())
         return false ;
   }

   // imposto ricalcolo della grafica
   m_OGrMgr.Reset() ;

   return true ;
}

//----------------------------------------------------------------------------
bool
SurfBezier::UnprojectPointFromStm( int nT, const Point3d& ptI, Point3d& ptSP, int nIL) const 
{
   return UnprojectPointFromStm( nT, ptI, ptSP, nIL, P_INVALID) ;
}

//----------------------------------------------------------------------------
bool
SurfBezier::UnprojectPointFromStm( int nT, const Point3d& ptI, Point3d& ptSP, int nIL, const Point3d& ptIPrevOrNext, bool* bThroughEdge) const 
{
   ptSP = ORIG ;
   if ( bThroughEdge != nullptr)
      *bThroughEdge = false ;
   // dato un punto sulla trimesh ausiliaria, ne ricavo le coordinate parametriche
   const ISurfTriMesh* pSurfTm = GetAuxSurf() ;
   int nTriaIndex = nT ;
   if ( nT == -1) {
      DistPointSurfTm distPtStm0( ptI, *pSurfTm) ;
      distPtStm0.GetMinDistTriaIndex( nTriaIndex) ;
   }
   // aggiungo il primo punto
   // devo subito capire se sono in un polo o no
   // se sono in polo e mi hanno passato un punto precedente allora devo prendere il triangolo di quel punto
   bool bIsPole = false ;
   int nInters = 0 ;
   INTVECTOR vInters(4) ;
   fill( vInters.begin(), vInters.end(), 0) ;
   // se il vettore dei poli non è stato riempito vuol dire che quando è stata creata la superficie non è stata chiamata la funzione CalcPoles
   if ( int( m_vbPole.size()) == 0)
      return false ;
   if ( m_vbPole[0] ||  m_vbPole[1]  ||  m_vbPole[2]  ||  m_vbPole[3]  ||  m_bClosedU  || m_bClosedV) {
      // scorro sugli edge
      for ( int c = 0 ; c < 4 ; ++c) {
         // scorro sui tratti che compongono l'edge
         for ( int i = 0 ; i < int( m_mCCEdge[c].size()) ; ++i) {
            if ( ! m_mCCEdge[c][i]->IsValid()) {
               Point3d pt ;
               if ( ! m_mCCEdge[c][i]->GetOnlyPoint(pt)) 
                  return false ;
               vInters[c] = AreSamePointApprox( pt, ptI) ? 1 : 0 ;
               nInters += vInters[c] ;
            }
            else {
               vInters[c] = m_mCCEdge[c][i]->IsPointOn(ptI) ? 1 : 0 ;
               nInters += vInters[c] ;
            }
         }
      }

      // se ho tre intersezioni vuol dire che un lato è collassato in un punto e il punto di cui voglio la controimmagine è esattamente nel polo
      if ( nInters == 3  ||  ( m_bClosedU &&  ( vInters[1] == 1 || vInters[3] == 1))  ||  ( m_bClosedV  &&  ( vInters[1] == 1  ||  vInters[3] == 1))) {
         if ( nInters == 3)
            bIsPole = true ;
         // visto che sono in un polo o su un lato di chiusura devo verificare di aver ricevuto il triangolo giusto
         // se è stato passato il punto successivo o precedente mi sposto verso quello e ricalcolo il triangolo di appartenenza
         if ( ! ptIPrevOrNext.IsValid())
            return false ;
         if ( bThroughEdge != nullptr)
            *bThroughEdge = true ;
         Point3d ptI2 = ptI + ( ptIPrevOrNext - ptI) * EPS_SMALL ;
         // ricalcolo il triangolo di appartenenza
         DistPointSurfTm dPtStm( ptI2, *pSurfTm) ;
         dPtStm.GetMinDistTriaIndex( nTriaIndex) ;
      }
   }

   // recupero i dati dei vertici del triangolo che fa intersezione
   int nVert[3] ;
   pSurfTm->GetTriangle( nTriaIndex, nVert) ;
   PNTVECTOR vPtPa(3) ;
   pSurfTm->GetVertexParam( nVert[0], vPtPa[0].x,vPtPa[0].y) ;
   pSurfTm->GetVertexParam( nVert[1], vPtPa[1].x,vPtPa[1].y) ;
   pSurfTm->GetVertexParam( nVert[2], vPtPa[2].x,vPtPa[2].y) ;
   PNTVECTOR vPT(3) ;
   pSurfTm->GetVertex( nVert[0], vPT[0]) ;
   pSurfTm->GetVertex( nVert[1], vPT[1]) ;
   pSurfTm->GetVertex( nVert[2], vPT[2]) ;
   // se la superficie è chiusa controllo se devo tenere conto della periodicità nel prendere le coordinate parametriche dei vertici
   double dParamH, dParamL ;
   if ( m_bClosedU || m_bClosedV) {
      dParamH = m_nSpanV * SBZ_TREG_COEFF ;
      dParamL = m_nSpanU * SBZ_TREG_COEFF ;
      // devo trovare il lato più lungo e confrontarlo con le dimensioni dello spazio parametrico
      Vector3d vtDir ;
      double dDist ;
      if ( DistXY( vPtPa[0], vPtPa[1]) > DistXY( vPtPa[1], vPtPa[2])  &&  DistXY( vPtPa[0], vPtPa[1]) > Dist( vPtPa[0], vPtPa[2])){
         vtDir = vPtPa[1] - vPtPa[0] ;
         dDist = DistXY( vPtPa[0], vPtPa[1]) ;
      }
      else if ( DistXY( vPtPa[1], vPtPa[2]) > DistXY( vPtPa[0], vPtPa[1])  &&  DistXY( vPtPa[1], vPtPa[2]) > Dist( vPtPa[0], vPtPa[2])){
         vtDir = vPtPa[2] - vPtPa[1] ;
         dDist = DistXY( vPtPa[1], vPtPa[2]) ;
      }
      else if ( DistXY( vPtPa[0], vPtPa[2]) > DistXY( vPtPa[0], vPtPa[1])  &&  DistXY( vPtPa[0], vPtPa[2]) > Dist( vPtPa[1], vPtPa[2])){
         vtDir = vPtPa[2] - vPtPa[0] ;
         dDist = DistXY( vPtPa[0], vPtPa[2]) ;
      }
      vtDir.Normalize() ;
      // se la dimensione maggiore è grande come la dimensione dello spazio parametrico allora potrebbe essere che le coordinate parametriche di un vertice
      // siano da correggere per periodicità
      if ( m_bClosedU  &&  abs(vtDir.x) > abs( vtDir.y) && dDist > dParamL * 0.5 ) {
         // trovo se dei vertici del triangolo sono sul bordo dello spazio parametrico
         INTVECTOR vOn(3) ;
         fill( vOn.begin(), vOn.end(), -1) ;
         int nVertOnPole = -1 ;
         INTVECTOR vEdgesClosed = { 1, 3} ;
         // scorro sui vertici
         for ( int p = 0 ; p < 3; ++p ) {
            // scorro sugli edge
            for ( int ed : vEdgesClosed) {
               // scorro sui tratti che compongono l'edge
               for ( int i = 0 ; i < int( m_mCCEdge[ed].size()) ; ++i) {
                  if ( ! m_mCCEdge[ed][i]->IsValid()) {
                     Point3d pt ;
                     if ( ! m_mCCEdge[ed][i]->GetOnlyPoint(pt))
                        return false ;
                     if ( AreSamePointApprox( pt, vPT[p])) {
                        vOn[p] = ed ;
                        // se un vertice sta su un polo me lo segno
                        nVertOnPole = p ;
                     }
                  }
                  else {
                     if (m_mCCEdge[ed][i]->IsPointOn(vPT[p]) && vOn[p] == -1 )
                        vOn[p] = ed ;
                  }
               }
            }
         }
         // controllo che almeno un vertice sia su un edge e se è l'unico vertice, che non sia su un polo
         if ( vOn[0] > 0  ||  vOn[1] > 0  ||  vOn[2] > 0) {
            // se ho più un vertice sul lato oppure se ne ho solo uno ma non è sul polo allora procedo alla correzione delle coordinate
            if ( vOn[0] * vOn[1] * vOn[2] < 0  ||
               (vOn[0] > 0  &&  vOn[0] != nVertOnPole)  ||
               (vOn[1] > 0  &&  vOn[1] != nVertOnPole)  ||
               (vOn[2] > 0  &&  vOn[2] != nVertOnPole)) {
               double dRightX ;
               // tengo per buone le coordinate dei vertici che NON sono sul bordo dello spazio parametrico
               for ( int p = 0 ; p < 3; ++p) {
                  if ( vOn[p] == -1) {
                     dRightX = vPtPa[p].x ;
                     break ;
                  }
               }
               for ( int p = 0 ; p < 3; ++p) {
                  if ( abs(vPtPa[p].x - dRightX) > EPS_SMALL ) {
                     if ( vPtPa[p].x < EPS_SMALL)
                        vPtPa[p].x = dParamL ;
                     else
                        vPtPa[p].x = 0 ;
                  }
               }
            }
         }
      }
      else if ( m_bClosedV  &&  abs(vtDir.y) > abs(vtDir.x)  &&  dDist > dParamH * 0.5) {
         INTVECTOR vOn(3) ;
         fill( vOn.begin(), vOn.end(), -1) ;
         int nVertOnPole = -1 ;
         INTVECTOR vEdgesClosed = { 0, 2} ;
         // scorro sui vertici
         for ( int p = 0 ; p < 3; ++p ) {
            //scorro sugli edge
            for (int ed : vEdgesClosed) {
               // scorro sui tratti che compongono l'edge
               for ( int i = 0 ; i < int( m_mCCEdge[ed].size()) ; ++i) {
                  if ( ! m_mCCEdge[ed][i]->IsValid()) {
                     Point3d pt ;
                     if ( ! m_mCCEdge[ed][i]->GetOnlyPoint( pt))
                        return false ;
                     if ( AreSamePointApprox( pt, vPT[p])) {
                        vOn[p] = ed ;
                        // se un vertice sta su un polo me lo segno
                        nVertOnPole = p ;
                     }
                  }
                  else {
                     if( m_mCCEdge[ed][i]->IsPointOn( vPT[p])  &&  vOn[p] == -1)
                        vOn[p] = ed ;
                  }
               }
            }
         }
         // controllo che almeno un vertice sia su un edge
         if ( vOn[0] > 0  ||  vOn[1] > 0  ||  vOn[2] > 0) {
            // se ho più un vertice sul lato oppure se ne ho solo uno ma non è sul polo allora procedo alla correzione delle coordinate
            if ( vOn[0] * vOn[1] * vOn[2] < 0  ||
               (vOn[0] > 0  &&  vOn[0] != nVertOnPole)  ||
               (vOn[1] > 0  &&  vOn[1] != nVertOnPole)  ||
               (vOn[2] > 0  &&  vOn[2] != nVertOnPole)) {
               double dRightY ;
               // tengo per buone le coordinate dei vertici che NON sono sul bordo dello spazio parametrico
               for ( int p = 0 ; p < 3; ++p) {
                  if ( vOn[p] == -1) {
                     dRightY = vPtPa[p].y ;
                     break ;
                  }
               }
               for ( int p = 0 ; p < 3; ++p) {
                  if ( abs(vPtPa[p].y - dRightY) > EPS_SMALL) {
                     if ( vPtPa[p].y < EPS_SMALL)
                        vPtPa[p].y = dParamH ;
                     else
                        vPtPa[p].y = 0 ;
                  }
               }
            }
         }
      }
   }
   // devo anche tener conto della possibilità che i lati siano collassati in poli
   if ( bIsPole) {
      int nInters = -1 ;
      for ( int c = 0 ; c < 4 ; ++c) {
         if ( ( c == 0 && vInters[0] == vInters[3])  ||
              ( c != 0 && vInters[c] == vInters[c - 1])){
            nInters = c ;
            break ;
         }
      }
      // se non ho trovato il lato su cui ho il polo
      if ( nInters == -1)
         return false ;
      // trovo quale vertice è sull'edge di polo
      BOOLVECTOR vbOn(3) ;
      fill( vbOn.begin(), vbOn.end(), false) ;
      for ( int p = 0 ; p < 3; ++p ) {
         for ( int c = 0 ; c < 4; ++c) {
            for( int i = 0 ; int( m_mCCEdge[c].size()) ; ++i) {
               if ( ! m_mCCEdge[c][i]->IsValid()) {
                  Point3d pt ;
                  if ( ! m_mCCEdge[c][i]->GetOnlyPoint( pt))
                     return false ;
                  vbOn[p] = vbOn[p]  ||  AreSamePointApprox( pt, vPT[p]) ;
               }
               else
                  vbOn[p] = vbOn[p]  ||  m_mCCEdge[c][i]->IsPointOn( vPT[p]) ;
            }
         }
      }
      // trovo la coordinata giusta da tenere ( x o y a seconda dell'edge)
      double dRightX, dRightY ;
      for ( int p = 0 ; p < 3; ++p) {
         if ( ! vbOn[p]) {
            if ( nInters == 0 || nInters == 2) {
               dRightX = vPtPa[p].x ;
               dRightY = nInters == 0 ? dParamH : 0 ;
            }
            else if ( nInters == 1 || nInters == 3) {
               dRightX = nInters == 1 ? 0 : dParamL ;
               dRightY = vPtPa[p].y ;
            }
         }
      }
      // correggo le coordinate del punto sull'edge di polo
      for ( int p = 0 ; p < 3 ; ++p) {
         if ( vbOn[p]) {
            vPtPa[p].x = dRightX ;
            vPtPa[p].y = dRightY ;
         }
      }
   }

   // se l'intersezione era su un vertice ( NON DI POLO) restituisco le coordinate parametriche del vertice
   if ( nIL == 3  && ! bIsPole) {
      if ( AreSamePointApprox(ptI, vPT[0]))
         ptSP = vPtPa[0] ;
      else if ( AreSamePointApprox(ptI, vPT[1]))
         ptSP = vPtPa[1] ;
      else if ( AreSamePointApprox(ptI, vPT[2]))
         ptSP = vPtPa[2] ;
      return true ;
   }
   // calcolo approssimativamente le coordinate nello spazio parametrico del punto di intersezione
   // quindi prima calcolo la composizione lineare tra i vertici del triangolo per ottenere il punto di intersezione
   Eigen::Matrix3d mA ;
   mA.col(0) << vPT[0].x, vPT[0].y , vPT[0].z ;
   mA.col(1) << vPT[1].x, vPT[1].y , vPT[1].z ;
   mA.col(2) << vPT[2].x, vPT[2].y , vPT[2].z ;
   int nCount = 0 ;
   Vector3d vtMod( 0, 0, 0) ;
   while ( abs( mA.determinant()) < EPS_SMALL) {
      if ( nCount %3 ==  0) {
         mA.row(0) << vPT[0].x + 1, vPT[1].x + 1, vPT[2].x + 1 ;
         vtMod.x += 1 ;
      }
      else if ( nCount %3 == 1) {
         mA.row(1) << vPT[0].y + 1, vPT[1].y + 1, vPT[2].y + 1 ;
         vtMod.y += 1 ;
      }
      else if ( nCount %3 == 2) {
         mA.row(2) << vPT[0].z + 1, vPT[1].z + 1, vPT[2].z + 1 ;
         vtMod.z += 1 ;
      }
      ++nCount ;
      if ( nCount == 10)
         return false ;
   }
   Point3d ptNewI = ptI + vtMod ;
   Eigen::Vector3d b ( ptNewI.x, ptNewI.y, ptNewI.z) ;
   Eigen::Vector3d x = mA.fullPivLu().solve(b) ;
   // applico questa composizione alle loro coordinate parametriche
   Eigen::Matrix3d mB ;
   mB.col(0) << vPtPa[0].x, vPtPa[0].y, 0 ;
   mB.col(1) << vPtPa[1].x, vPtPa[1].y, 0 ;
   mB.col(2) << vPtPa[2].x, vPtPa[2].y, 0 ;
   Eigen::Vector3d ptParam = mB * x ;
   IncreaseUV( ptSP.x, ptParam.x(), true) ;
   IncreaseUV( ptSP.y, ptParam.y(), false) ;
   return true ;
}

//----------------------------------------------------------------------------
bool
SurfBezier::UnprojectPoint( const Point3d& pt3D, Point3d& ptParam, const Point3d& ptIPrev, bool* bThroughEdge, const Plane3d* pPlCut) const 
{
   // dato il punto pt3D sulla superficie di Bezier si cercano le coordinate parametriche ( ptParam) , iterativamente con Newton
   // trovato un primo candidato ptParam, ne calcolo l'immagine sulla superficie ( ptBez) e ne calcolo la distanza con il punto pt3D
   // ripeto cercando di avvicinarmi il più possibile
   // per trovare il primo punto trovo il triangolo della trimesh ausiliaria più vicino e il punto più vicino
   DistPointSurfTm dptSurfTm( pt3D, *GetAuxSurf()) ;
   Point3d ptI ; dptSurfTm.GetMinDistPoint( ptI) ;
   if ( ! UnprojectPointFromStm( -1, ptI, ptParam, 5, ptIPrev, bThroughEdge))
      return false ;
   Point3d ptBez ;
   GetPointD1D2( ptParam.x / SBZ_TREG_COEFF, ptParam.y / SBZ_TREG_COEFF, ISurfBezier::FROM_MINUS, ISurfBezier::FROM_MINUS, ptBez) ;
   // usando un algoritmo di newton cerco di avvicinarmi il più possibile al punto

   double dDistNew = pPlCut == nullptr ? Dist( pt3D, ptBez) : abs(DistPointPlane( ptBez, *pPlCut)) ;
   double dDistPre ;
   double dDist0, dDist1;

   int nCount = 0 ;
   double dh = EPS_SMALL ;
   // metodo di newton in più dimensioni
   // vario sia il parametro U che il parametro V e verifico se la distanza dalla retta diminuisce per scostamenti positivi o negativi.
   bool bRetry = false ;
   double dApproach = 0.01 ;
   bool bDesperate = false ;
   double dAng = 1 ;
   double dr = 5 ;
   double dfdU, dfdV ;

   while ( dDistNew > 2 * EPS_SMALL  &&  nCount < 100) {

      if ( ! bRetry) {
         dDistPre = dDistNew ;
         // derivata in U
         Point3d ptIBzNew1 ;
         //GetPointD1D2( ( ptParam.x + dh)  / SBZ_TREG_COEFF, ptParam.y  / SBZ_TREG_COEFF, ISurfBezier::FROM_MINUS, ISurfBezier::FROM_MINUS, ptIBzNew1) ;
         double dUh ; IncreaseUV( ptParam.x, dh, true, &dUh, false) ;
         GetPointD1D2( dUh / SBZ_TREG_COEFF, ptParam.y  / SBZ_TREG_COEFF, ISurfBezier::FROM_MINUS, ISurfBezier::FROM_MINUS, ptIBzNew1) ;
         dDist0 = pPlCut == nullptr ? Dist( pt3D, ptIBzNew1) : abs(DistPointPlane( ptIBzNew1, *pPlCut)) ;
         dfdU = ( dDist0 - dDistPre) / dh ;
         // derivata in V
         Point3d ptIBzNew2 ;
         //GetPointD1D2( ptParam.x  / SBZ_TREG_COEFF, ( ptParam.y + dh) / SBZ_TREG_COEFF, ISurfBezier::FROM_MINUS, ISurfBezier::FROM_MINUS, ptIBzNew2) ;
         double dVh ; IncreaseUV( ptParam.y, dh, false, &dVh, false) ;
         GetPointD1D2( ptParam.x  / SBZ_TREG_COEFF, dVh / SBZ_TREG_COEFF, ISurfBezier::FROM_MINUS, ISurfBezier::FROM_MINUS, ptIBzNew2) ;
         dDist1 = pPlCut == nullptr ? Dist( pt3D, ptIBzNew2) : abs(DistPointPlane( ptIBzNew2, *pPlCut)) ;
         dfdV = ( dDist1 - dDistPre) / dh ;
      }
      // calcolo le nuove coordinate
      Vector3d vtDir ;
      double dASum = abs(dfdU) + abs(dfdV) ;
      double dSSum = sqrt( pow(dfdU,2) + pow( dfdV,2)) ;
      vtDir.Set( - dfdU, - dfdV, 0) ;
      if ( ! vtDir.Normalize() )
         vtDir.Set( - dfdU / dSSum, - dfdV / dSSum, 0) ;
      dr = dDistPre / dASum ;
      // in modalità Retry riduco lo spostamento
      vtDir *= dr * ( bRetry ? 0.1 : 0.5) ;
      if ( bDesperate) {
         // in depserate mode riduco lo spostamento e comincio a cambiare progressivamente la direzione oscillando tra destra e sinistra della direzione "naturale"
         vtDir *= 0.5 ;
         dAng *= -1.25 ;
         // riduco ulteriormente lo spostamento dopo qualche giro in desperate mode
         if ( abs( dAng) > 5)
            vtDir *= 0.5 ;
         vtDir.Rotate( Z_AX, dAng) ;
      }
      Point3d ptCopy = ptParam ;
      IncreaseUV( ptParam, vtDir, nullptr, true) ;
      // calcolo la nuova distanza tra il punto di partenza e quello che sto trovando con le coordinate attuali
      GetPointD1D2( ptParam.x / SBZ_TREG_COEFF, ptParam.y / SBZ_TREG_COEFF, ISurfBezier::FROM_MINUS, ISurfBezier::FROM_MINUS, ptBez) ;
      dDistNew = pPlCut == nullptr ? Dist( pt3D, ptBez) : abs(DistPointPlane( ptBez, *pPlCut)) ;

      dApproach = dDistPre - dDistNew ;
      // se ho peggiorato la situazione rispetto allo step precedente torno indietro e vado in Retry mode
      if ( dApproach < EPS_ZERO) {
         if ( bRetry)
            bDesperate = true ;  // entro in desperate mode
         ptParam = ptCopy ;
         bRetry = true ;
         // se ho già fatto 9 giri in desperate mode allora esco
         if ( abs( dAng) > 8)
            break ;
      }
      else {
         bRetry = false ;
         bDesperate = false ;
         dAng = 1 ;
      }
      ++nCount ;
   }

   // se il punto era su un edge allora verifico che sia ancora su un edge, sennò ce lo riporto
   // guardo a quale dei due lati sono più vicino
   // devo distinguere il caso di un triangolo a metà dello spazio, con un vertice su un lato di polo, ma senza vertici su lati di chiusura
   if ( bThroughEdge != nullptr  &&  *bThroughEdge) {
      if ( m_bClosedU) {
         if ( ptParam.x < 1)
            ptParam.x = 0 ;
         else if ( abs( m_nSpanU * SBZ_TREG_COEFF - ptParam.x) < 1)
            ptParam.x = m_nSpanU * SBZ_TREG_COEFF ;
      }
      if ( m_bClosedV) {
         if ( ptParam.y < 1)
            ptParam.y = 0 ;
         else if ( abs( m_nSpanV * SBZ_TREG_COEFF - ptParam.y) < 1)
            ptParam.y = m_nSpanV * SBZ_TREG_COEFF ;
      }
   }

   return nCount != 100  ||  (dDistNew < dDistPre ? dDistNew : dDistPre) < 10 * EPS_SMALL ;
}

//----------------------------------------------------------------------------
bool
SurfBezier::CalcPoles( void)
{
   // la funzione identifica se degli edge della superficie non trimmata sono in realtà dei poli
   for ( int i = 0 ; i < 4 ; ++i)
      m_vbPole.emplace_back( true) ;
   // scorro i punti di controllo e vedo subito
   bool bOk = false ;
   bool bPole0 = true, bPole1 = true ;
   Point3d ptU0, ptU1 ;
   // controllo l'edge 0 e 2 per vedere se tutti i punti dell'edge sono coincidenti
   Point3d ptP00 = GetControlPoint( 0, &bOk) ;
   Point3d ptP10 = GetControlPoint( m_nDegU * m_nSpanU, &bOk) ;
   for ( int i = 1 ; i < m_nDegV * m_nSpanV + 1 ; ++ i) {
      ptU0 = GetControlPoint( i * ( m_nDegU * m_nSpanU + 1), &bOk) ;
      bPole0 = bPole0  &&  AreSamePointApprox( ptP00, ptU0) ;
      ptU1 = GetControlPoint( ( i + 1) * ( m_nDegU * m_nSpanU + 1) - 1, &bOk) ;
      bPole1 = bPole1  &&  AreSamePointApprox( ptP10, ptU1) ;
      if ( ! bPole0 && ! bPole1)
         break ;
   }
   m_vbPole[1] = bPole0 ;  // u = 0 corrisponde all'edge 1
   m_vbPole[3] = bPole1 ;  // u = 1 corrisponde all'edge 3
   // controllo l'edge 1 e 3 per vedere se tutti i punti dell'edge sono coincidenti
   Point3d ptV0, ptV1 ;
   Point3d ptP01 = GetControlPoint( ( m_nDegU * m_nSpanU + 1) * ( m_nDegV * m_nSpanV), &bOk) ;
   bPole0 = true ;
   bPole1 = true ;
   for ( int i = 1 ; i < m_nDegU * m_nSpanU + 1 ; ++ i) {
      ptV0 = GetControlPoint( i, &bOk) ;
      bPole0 = bPole0  &&  AreSamePointApprox( ptP00, ptV0) ;
      ptV1 = GetControlPoint( i + ( m_nDegU * m_nSpanU + 1) * ( m_nDegV * m_nSpanV), &bOk) ;
      bPole1 = bPole1  &&  AreSamePointApprox( ptP01, ptV1) ;
      if ( ! bPole0 && ! bPole1)
         break ;
   }
   m_vbPole[0] = bPole1 ;  // v = 1 corrisponde all'edge 0
   m_vbPole[2] = bPole0 ;  // v = 0 corrisponde all'edge 2

   return true ;
}

//----------------------------------------------------------------------------
bool
SurfBezier::GetLoops( ICRVCOMPOPOVECTOR& vCC, bool bLineOrBezier, int nEdge) const
{
   // se nEdge non è definito ( == -1) allora restituisco tutti gli edge in ordine
   // se bOpenOrAll è true allora restituisco solo gli edge aperti e che non sono di polo
   if ( m_pSTM == nullptr)
      GetAuxSurf() ;
   
   // se la superficie non è trimmata mi basta recuperare gli edge della superficie
   if ( ! m_bTrimmed) {
      // se decidessi di non restituire gli edge chiusi e i poli posso discriminare qui
      if ( nEdge == -1 ) {
            vCC = ICRVCOMPOPOVECTOR(4) ;
            if ( ! m_bClosedV ) {
               if ( ! m_vbPole[0])
                  vCC[0].Set( GetSingleEdge3D( bLineOrBezier, 0)) ;
               if ( ! m_vbPole[2])
                  vCC[2].Set( GetSingleEdge3D( bLineOrBezier, 2)) ;
            }
            if ( ! m_bClosedU ) {
               if ( ! m_vbPole[1])
                  vCC[1].Set( GetSingleEdge3D( bLineOrBezier, 1)) ;
               if ( ! m_vbPole[3])
                  vCC[3].Set( GetSingleEdge3D( bLineOrBezier, 3)) ;
            }
            //// se li volessi restituire tutti
            //vCC.emplace_back( GetSingleEdge3D( bLineOrBezier, 0)) ;
            //vCC.emplace_back( GetSingleEdge3D( bLineOrBezier, 1)) ;
            //vCC.emplace_back( GetSingleEdge3D( bLineOrBezier, 2)) ;
            //vCC.emplace_back( GetSingleEdge3D( bLineOrBezier, 3)) ;
      }
      else {
         if ( (((nEdge == 0 ||  nEdge == 2) &&  ! m_bClosedV)  ||  ((nEdge == 1 ||  nEdge == 3) &&  ! m_bClosedU))  &&  ! m_vbPole[nEdge])
            vCC.emplace_back( GetSingleEdge3D( bLineOrBezier, nEdge)) ;
      }
   }
   // se la superficie è trimmata devo recuperare i loop dello spazio parametrico
   else {
      // devo ricostruire i bordi aperti in caso di superficie chiusa
      // costruisco gli edge
      ICurveLine* pCLU0( CreateCurveLine()) ;
      ICurveLine* pCLU1( CreateCurveLine()) ;
      ICurveLine* pCLV0( CreateCurveLine()) ;
      ICurveLine* pCLV1( CreateCurveLine()) ;
      if ( m_bClosedU ) {
         pCLU0->Set( Point3d( 0, m_nSpanV * SBZ_TREG_COEFF, 0), ORIG) ;
         pCLU1->Set( Point3d( m_nSpanU * SBZ_TREG_COEFF, 0, 0), Point3d( m_nSpanU * SBZ_TREG_COEFF, m_nSpanV * SBZ_TREG_COEFF, 0)) ;
      }
      if ( m_bClosedV ) {
         pCLV0->Set( ORIG, Point3d( m_nSpanU * SBZ_TREG_COEFF, 0, 0)) ;
         pCLV1->Set( Point3d( m_nSpanU * SBZ_TREG_COEFF, m_nSpanV * SBZ_TREG_COEFF, 0), Point3d( 0, m_nSpanV * SBZ_TREG_COEFF, 0)) ;
      }
      ICRVLINEPOVECTOR vEdge ;
      // li metto con l'ordine degli edge
      vEdge.emplace_back( pCLV1) ;
      vEdge.emplace_back( pCLU0) ;
      vEdge.emplace_back( pCLV0) ;
      vEdge.emplace_back( pCLU1) ;
      // costruisco la mappa delle intersezioni, trovando tutte le intersezioni tra i trim e il loop esterno dello spazio parametrico
      // i 4 elementi fanno riferimento ai 4 edge. per ogni loop viene salvato il vettore delle intersezioni con tale edege
      vector<unordered_map<int,ICCIVECTOR>> vmInters(4) ;
      // comincio anche a creare il vettore di vettori associati ad ogni loop che fa almeno un'intersezione con un edge.
      // in ogni vettore verranno salvati l'indice del loop all'interno del vettore dei loop e gli indici di tutti i loop che deriveranno da sue divisioni in più curve
      unordered_map<int,INTVECTOR> mSplitLoop ;
      int nOriginalLoops = int( m_vCCLoop.size()) ;
      for ( int l = 0 ; l < nOriginalLoops ; ++ l) {
         bool bInters = false ;
         // devo controllare se i loop si appoggiano a bordi chiusi
         // se ne ho più di uno che si appoggia allora devo vedere se trovo dei gruppi di loop che in realtà in 3d sono un edge unico, ma nel 2d risultano a cavallo dell'edge di chiusura
         PtrOwner<ICurveComposite> pLoop( m_vCCLoop[l]->Clone()) ;
         for ( int j = 0 ; j < 4 ; ++j) {
            // il primo è un loop di trim, il secondo è un edge dello spazio parametrico
            IntersCurveCurve icc( *pLoop, *vEdge[j]) ;
            ICCIVECTOR vICC ;
            for ( int i = 0 ; i < int(icc.GetIntersCount()) ; ++i ) {
               IntCrvCrvInfo iccInfo ;
               icc.GetIntCrvCrvInfo( i, iccInfo) ;
               vICC.push_back( iccInfo) ;
            }
            if ( int(vICC.size()) > 0) {
               vmInters[j].insert(pair<int, ICCIVECTOR>(l, vICC)) ;
               bInters = true ;
            }
         }
         if ( bInters)
            mSplitLoop.insert( pair<int,INTVECTOR>( l, INTVECTOR({l}))) ;
      }

      if ( m_bClosedV) {
         // scorro i loop che fanno intersezioni con l'edge 0 ( V=1) e le riporto sull'edge opposto
         for ( pair<int, ICCIVECTOR> vInters0 : vmInters[0]) {
            // scorro le intersezioni del loop
            for ( IntCrvCrvInfo inters0 : vInters0.second) {  
               if ( ! inters0.bOverlap)
                  continue ;
               Point3d ptStart0 = inters0.IciA[0].ptI ;
               Point3d ptEnd0 = inters0.IciA[1].ptI ;
               // porto i punti sull'edge opposto
               ptStart0.y = 0 ;
               ptEnd0.y = 0 ;
               // scorro le intersezioni sull'edge opposto e quando trovo un loop che tra le sue intersezioni contiene uno dei punti che ho appena portato
               // su qusto edge allora cancello la parte comune
               for ( pair<int, ICCIVECTOR> vInters2 : vmInters[2] ) {
                  // scorro le intersezioni di questo loop
                  for ( IntCrvCrvInfo inters2 : vInters2.second) {
                     // se lo start o end del loop corrente è compreso tra lo start e l'end di un loop che fa intesezione sull'edge opposto allora 
                     // devo cancellare la parte comune
                     Point3d ptStart2 = inters2.IciA[0].ptI ;
                     Point3d ptEnd2 = inters2.IciA[1].ptI ;
                     // come ptStart2 prendo quello con la x maggiore
                     if ( ptStart2.x < ptEnd2.x)
                        swap( ptStart2, ptEnd2) ;
                     if ( (( ptEnd0.x - EPS_SMALL< ptStart2.x  &&  ptStart2.x < ptStart0.x + EPS_SMALL)  ||  ( ptEnd0.x - EPS_SMALL < ptEnd2.x  &&  ptEnd2.x < ptStart0.x + EPS_SMALL))  ||
                          (( ptEnd2.x - EPS_SMALL < ptStart0.x  &&  ptStart0.x < ptStart2.x + EPS_SMALL)  ||  ( ptEnd2.x - EPS_SMALL < ptEnd0.x  &&  ptEnd0.x < ptStart2.x + EPS_SMALL))) {
                        PtrOwner<ICurveLine> pCL( CreateBasicCurveLine()) ;
                        pCL->Set( ptStart0, ptEnd0) ;
                        // devo scorrere su tutte le curve che sono state ottenute dallo split del loop originale
                        for ( int w = 0 ; w < int( mSplitLoop[vInters2.first].size()) ; ++w) {
                           int nIndex2 = mSplitLoop[vInters2.first][w] ;
                           IntersCurveCurve icc( *m_vCCLoop[nIndex2], *pCL) ;
                           IntCrvCrvInfo iccInfo ; icc.GetIntCrvCrvInfo( 0, iccInfo) ;
                           for ( int k = 0 ; k < icc.GetIntersCount() ; ++k) {
                              icc.GetIntCrvCrvInfo( k, iccInfo) ;
                              if ( iccInfo.bOverlap)
                                 break ;
                           }
                           if ( ! iccInfo.bOverlap)
                              continue ;
                           // se parto da una curva chiusa semplicemente tolgo un pezzo
                           if ( m_vCCLoop[nIndex2]->IsClosed()) {
                              ICurveComposite* pCC2( GetCurveComposite( m_vCCLoop[nIndex2]->CopyParamRange( iccInfo.IciA[1].dU, iccInfo.IciA[0].dU))) ;
                              m_vCCLoop[nIndex2].Set( pCC2) ;
                           }
                           else {
                              // se la curva era già aperta allora otterrò due curve separate
                              ICurveComposite* pCC2a( GetCurveComposite( m_vCCLoop[nIndex2]->Clone())) ;
                              bool bDoneA = pCC2a->TrimEndAtParam( iccInfo.IciA[0].dU) ;
                              ICurveComposite* pCC2b( GetCurveComposite( m_vCCLoop[nIndex2]->Clone())) ;
                              bool bDoneB = pCC2b->TrimStartAtParam( iccInfo.IciA[1].dU) ;
                              if ( bDoneA) {
                                 m_vCCLoop[nIndex2].Set( pCC2a) ;
                                 if ( bDoneB) {
                                    m_vCCLoop.emplace_back( pCC2b) ;
                                    mSplitLoop[vInters2.first].push_back( m_vCCLoop.size() - 1) ;
                                 }
                              }
                              else if ( bDoneB)
                                 m_vCCLoop[nIndex2].Set( pCC2b) ;
                           }
                           
                           // per togliere la parte comune al loop corrente devo riportare i punti di intersezione sull'edge di partenza
                           Point3d ptOverlapS = iccInfo.IciB[0].ptI ;
                           Point3d ptOverlapE = iccInfo.IciB[1].ptI ;
                           ptOverlapS.y = m_nSpanV * SBZ_TREG_COEFF ;
                           ptOverlapE.y = m_nSpanV * SBZ_TREG_COEFF ;
                           double dStart0, dEnd0 ;
                           // scorro tutte le curve in cui è stato splittato il loop corrente
                           for ( int j = 0 ; j < int( mSplitLoop[vInters0.first].size()) ; ++j ) {
                              int nIndex0 = mSplitLoop[vInters0.first][j] ;
                              if ( ! m_vCCLoop[nIndex0]->GetParamAtPoint( ptOverlapS, dStart0) ||  ! m_vCCLoop[nIndex0]->GetParamAtPoint( ptOverlapE, dEnd0))
                                 continue ;
                              // se parto da una curva chiusa semplicemente tolgo un pezzo
                              if ( m_vCCLoop[nIndex0]->IsClosed()) {
                                 ICurveComposite* pCC0( GetCurveComposite( m_vCCLoop[nIndex0]->CopyParamRange( dStart0, dEnd0))) ;
                                 m_vCCLoop[nIndex0].Set( pCC0) ;
                              }
                              else {
                                 // se la curva era già aperta allora otterrò due curve separate
                                 ICurveComposite* pCC0a( GetCurveComposite( m_vCCLoop[nIndex0]->Clone())) ;
                                 bool bDoneA = pCC0a->TrimEndAtParam( dEnd0) ;
                                 ICurveComposite* pCC0b( GetCurveComposite( m_vCCLoop[nIndex0]->Clone())) ;
                                 bool bDoneB = pCC0b->TrimStartAtParam( dStart0) ;
                                 if ( bDoneA) {
                                    m_vCCLoop[nIndex0].Set( pCC0a) ;
                                    if ( bDoneB) {
                                       m_vCCLoop.emplace_back( pCC0b) ;
                                       mSplitLoop[vInters0.first].push_back( m_vCCLoop.size() - 1) ;
                                    }
                                 }
                                 else if ( bDoneB)
                                    m_vCCLoop[nIndex0].Set( pCC0b) ;
                              }
                              //break ;
                           }
                           //break ;
                        }
                     }
                  }
               }
            }
         }
      }

      if ( m_bClosedU) {
         // scorro i loop che fanno intersezioni con l'edge 1 ( U=0) e le riporto sull'edge opposto
         for ( pair<int, ICCIVECTOR> vInters1 : vmInters[1]) {
            // scorro le intersezioni del loop
            for ( IntCrvCrvInfo inters1 : vInters1.second) {
               if ( ! inters1.bOverlap)
                  continue ;
               Point3d ptStart1 = inters1.IciA[0].ptI ;
               Point3d ptEnd1 = inters1.IciA[1].ptI ;
               // porto i punti sull'edge opposto
               ptStart1.x = m_nSpanU * SBZ_TREG_COEFF ;
               ptEnd1.x = m_nSpanU * SBZ_TREG_COEFF ;
               // scorro le intersezioni sull'edge opposto e quando trovo un loop che tra le sue intersezioni contiene uno dei punti che ho appena portato
               // su qusto edge allora cancello la parte comune
               for ( pair<int, ICCIVECTOR> vInters3 : vmInters[3] ) {
                  // scorro le intersezioni di questo loop
                  for ( IntCrvCrvInfo inters3 : vInters3.second) {
                     // se lo start o end del loop corrente è compreso tra lo start e l'end di un loop che fa intesezione sull'edge opposto allora 
                     // devo cancellare la parte comune
                     Point3d ptStart3 = inters3.IciA[0].ptI ;
                     Point3d ptEnd3 = inters3.IciA[1].ptI ;
                     // come ptStart3 prendo quello con la y maggiore
                     if ( ptStart3.y < ptEnd3.y)
                        swap( ptStart3, ptEnd3) ;
                     if ( (( ptEnd1.y - EPS_SMALL< ptStart3.y  &&  ptStart3.y < ptStart1.y + EPS_SMALL)  ||  ( ptEnd1.y - EPS_SMALL < ptEnd3.y  &&  ptEnd3.y < ptStart1.y + EPS_SMALL))  ||
                        (( ptEnd3.y - EPS_SMALL < ptStart1.y  &&  ptStart1.y < ptStart3.y + EPS_SMALL)  ||  ( ptEnd3.y - EPS_SMALL < ptEnd1.y  &&  ptEnd1.y < ptStart3.y + EPS_SMALL))) {
                        PtrOwner<ICurveLine> pCL( CreateBasicCurveLine()) ;
                        pCL->Set( ptStart1, ptEnd1) ;
                        // devo scorrere su tutte le curve che sono state ottenute dallo split del loop originale
                        for ( int w = 0 ; w < int( mSplitLoop[vInters3.first].size()) ; ++w) {
                           int nIndex3 = mSplitLoop[vInters3.first][w] ;
                           IntersCurveCurve icc( *m_vCCLoop[nIndex3], *pCL) ;
                           IntCrvCrvInfo iccInfo ;
                           for ( int k = 0 ; k < icc.GetIntersCount() ; ++k) {
                              icc.GetIntCrvCrvInfo( k, iccInfo) ;
                              if ( iccInfo.bOverlap)
                                 break ;
                           }
                           if ( ! iccInfo.bOverlap)
                              continue ;
                           // se parto da una curva chiusa semplicemente tolgo un pezzo
                           if ( m_vCCLoop[nIndex3]->IsClosed()) {
                              ICurveComposite* pCC3( GetCurveComposite( m_vCCLoop[nIndex3]->CopyParamRange( iccInfo.IciA[1].dU, iccInfo.IciA[0].dU))) ;
                              m_vCCLoop[nIndex3].Set( pCC3) ;
                           }
                           else {
                              // se la curva era già aperta allora otterrò due curve separate
                              ICurveComposite* pCC3a( GetCurveComposite( m_vCCLoop[nIndex3]->Clone())) ;
                              bool bDoneA = pCC3a->TrimEndAtParam( iccInfo.IciA[0].dU) ;
                              ICurveComposite* pCC3b( GetCurveComposite( m_vCCLoop[nIndex3]->Clone())) ;
                              bool bDoneB = pCC3b->TrimStartAtParam( iccInfo.IciA[1].dU) ;
                              if ( bDoneA) {
                                 m_vCCLoop[nIndex3].Set( pCC3a) ;
                                 if ( bDoneB) {
                                    m_vCCLoop.emplace_back( pCC3b) ;
                                    mSplitLoop[vInters3.first].push_back( m_vCCLoop.size() - 1) ;
                                 }
                              }
                              else if ( bDoneB )
                                 m_vCCLoop[nIndex3].Set( pCC3b) ;
                           }

                           // per togliere la parte comune al loop corrente devo riportare i punti di intersezione sull'edge di partenza
                           Point3d ptOverlapS = iccInfo.IciB[0].ptI ;
                           Point3d ptOverlapE = iccInfo.IciB[1].ptI ;
                           ptOverlapS.x = 0 ;
                           ptOverlapE.x = 0 ;
                           double dStart1, dEnd1 ;
                           // scorro tutte le curve in cui è stato splittato il loop iniziale
                           for ( int j = 0 ; j < int( mSplitLoop[vInters1.first].size()) ; ++j ) {
                              int nIndex1 = mSplitLoop[vInters1.first][j] ;
                              if ( ! m_vCCLoop[nIndex1]->GetParamAtPoint( ptOverlapS, dStart1)  ||  ! m_vCCLoop[nIndex1]->GetParamAtPoint( ptOverlapE, dEnd1))
                                 continue ;
                              // se parto da una curva chiusa semplicemente tolgo un pezzo
                              if ( m_vCCLoop[nIndex1]->IsClosed()) {
                                 ICurveComposite* pCC1( GetCurveComposite( m_vCCLoop[nIndex1]->CopyParamRange( dStart1, dEnd1))) ;
                                 m_vCCLoop[nIndex1].Set( pCC1) ;
                              }
                              else {
                                 // se la curva era già aperta allora otterrò due curve separate
                                 ICurveComposite* pCC1a( GetCurveComposite( m_vCCLoop[nIndex1]->Clone())) ;
                                 bool bDoneA = pCC1a->TrimEndAtParam( dEnd1) ;
                                 ICurveComposite* pCC1b( GetCurveComposite( m_vCCLoop[nIndex1]->Clone())) ;
                                 bool bDoneB = pCC1b->TrimStartAtParam( dStart1) ;
                                 if ( bDoneA) {
                                 m_vCCLoop[nIndex1].Set( pCC1a) ;
                                    if ( bDoneB) {
                                       m_vCCLoop.emplace_back( pCC1b) ;
                                       mSplitLoop[vInters1.first].push_back( m_vCCLoop.size() - 1) ;
                                    }
                                 }
                                 else if ( bDoneB)
                                    m_vCCLoop[nIndex1].Set( pCC1b) ;
                              }
                              //break ;
                           }
                           //break ;
                        }
                     }
                  }
               }
            }
         }
      }

      // qui portei estrarre una funzione che proietta curve dal parametrico al 3D
      // scorro i gruppi di loop 2D formati da loop che partecipano alla formazione dello stesso loop nel 3D
      ICRVCOMPOPOVECTOR vpCCOpen ;
      for( int i = 0 ; i < int( m_vCCLoop.size()); ++i) {
         vpCCOpen.emplace_back(CreateBasicCurveComposite()) ;
         PolyLine pl3D ;
         // la composita è una spezzata composta da linee, quindi la ricostruisco come polyline in 3D
         Point3d pt ; m_vCCLoop[i]->GetStartPoint( pt) ;
         Point3d pt3D ; GetPointD1D2( pt.x / SBZ_TREG_COEFF, pt.y / SBZ_TREG_COEFF, ISurfBezier::FROM_MINUS, ISurfBezier::FROM_MINUS, pt3D) ;
         int nCount = 0 ;
         pl3D.AddUPoint( nCount, pt3D) ;
         ++ nCount ;
         // scorro le curve singole della composita
         for ( int k = 0 ; k < m_vCCLoop[i]->GetCurveCount() ; ++k){
            // recupero l'end point della curva e lo porto in 3D
            m_vCCLoop[i]->GetCurve( k)->GetEndPoint( pt) ;
            GetPointD1D2( pt.x / SBZ_TREG_COEFF, pt.y / SBZ_TREG_COEFF, ISurfBezier::FROM_MINUS, ISurfBezier::FROM_MINUS, pt3D) ;
            pl3D.AddUPoint( nCount, pt3D) ;
            ++ nCount ;
         }
         vpCCOpen.back()->FromPolyLine(pl3D) ;
      }
      // creo la chain a mano con le compo
      BOOLVECTOR vbAdded( vpCCOpen.size()) ;
      fill( vbAdded.begin(), vbAdded.end(), false) ;
      for ( int k = 0 ; k < int( vpCCOpen.size()); ++k ) {
         if ( vbAdded[k])
            continue ;
         PtrOwner<ICurveComposite> pCC( vpCCOpen[k]->Clone()) ;
         vbAdded[k] = true ;
         bool bAddedOne = true ;
         while( bAddedOne) {
            bAddedOne = false ;
            for ( int t = k ; t < int(vpCCOpen.size()); ++t ) {
               if ( vbAdded[t])
                  continue ;
               if ( pCC->AddCurve( vpCCOpen[t]->Clone())) {
                  vbAdded[t] = true ;
                  bAddedOne = true ;
               }
            }
         }
         // aggiungo la curva agli edge 3d da restituire
         vCC.emplace_back( Release( pCC)) ;
      }
   }
   return true ;
}

//----------------------------------------------------------------------------
ICurveComposite*
SurfBezier::GetSingleEdge3D( bool bLineOrBezier, int nEdge) const
{
   if ( m_mCCEdge.size() == 0  &&  bLineOrBezier)
      GetAuxSurf() ;
   // questa funzione dà per scontato che la superficie NON sia trimmata
   if ( nEdge < 0  || nEdge > 3  ||  m_bTrimmed)
      return nullptr ;
   ICurveComposite* pCrvCompo( CreateBasicCurveComposite()) ;
   switch ( nEdge) {
   case 0 : {
      // se il bool è true allora restituisco gli edge con la loro approssimazione in forma di linea spezzata 3D
      if ( bLineOrBezier)
         pCrvCompo = m_mCCEdge[nEdge][0]->Clone() ;
      // se il bool è falso restituisco le curve Bezier di edge
      else {
         //edge 0, scorro sulle patch in U
         for ( int i = 0 ; i < m_nSpanU ; ++i) {
            PtrOwner<ICurveBezier> pCrvBz0( CreateBasicCurveBezier()) ;
            if ( IsNull(pCrvBz0)  ||  ! pCrvBz0->Init( m_nDegU, m_bRat))
               return nullptr ;
            for ( int p = 0 ; p < m_nDegU + 1 ; ++p ) {
               int nIndex = ( m_nDegU * m_nSpanU + 1) * ( m_nDegV * m_nSpanV) + m_nDegU * i + p ;
               if ( ! m_bRat)
                  pCrvBz0->SetControlPoint( p, GetControlPoint( nIndex, nullptr)) ;
               else
                  pCrvBz0->SetControlPoint( p, GetControlPoint( nIndex, nullptr), GetControlWeight( nIndex, nullptr)) ;
            }
            pCrvCompo->AddCurve( Release( pCrvBz0)) ;
         }
      }
      break ;
   }
   case 1 : {
      if ( bLineOrBezier)
         pCrvCompo = m_mCCEdge[nEdge][0]->Clone() ;
      else {
         //edge 1, scorro sulle patch in V
         for ( int i = 0 ; i < m_nSpanV ; ++i) {
            PtrOwner<ICurveBezier> pCrvBz1( CreateBasicCurveBezier()) ;
            if ( IsNull(pCrvBz1)  ||  ! pCrvBz1->Init( m_nDegV, m_bRat))
               return nullptr ;
            for ( int p = 0 ; p < m_nDegU + 1 ; ++p ) {
               int nIndex = ( m_nDegU * m_nSpanU + 1) * p ;
               if ( ! m_bRat)
                  pCrvBz1->SetControlPoint( p, GetControlPoint( nIndex, nullptr)) ;
               else
                  pCrvBz1->SetControlPoint( p, GetControlPoint( nIndex, nullptr), GetControlWeight( nIndex, nullptr)) ;
            }
            pCrvCompo->AddCurve( Release( pCrvBz1)) ;
         }
      }
      break ;
   }
   case 2 : {
      if ( bLineOrBezier)
         pCrvCompo = m_mCCEdge[nEdge][0]->Clone() ;
      else {
         // edge 2, scorro sulle patch in U
         for ( int i = 0 ; i < m_nSpanU ; ++i) {
            PtrOwner<ICurveBezier> pCrvBz2( CreateBasicCurveBezier()) ;
            if ( IsNull(pCrvBz2)  ||  ! pCrvBz2->Init( m_nDegU, m_bRat))
               return nullptr ;
            for ( int p = 0 ; p < m_nDegU + 1 ; ++p ) {
               if ( ! m_bRat)
                  pCrvBz2->SetControlPoint( p, GetControlPoint(m_nDegU * i + p, nullptr)) ;
               else
                  pCrvBz2->SetControlPoint( p, GetControlPoint(m_nDegU * i + p, nullptr), GetControlWeight(m_nDegU * i + p, nullptr)) ;
            }
            pCrvCompo->AddCurve( Release( pCrvBz2)) ;
         }
      }
      break ;
   }
   case 3 : {
      if ( bLineOrBezier)
         pCrvCompo = m_mCCEdge[nEdge][0]->Clone() ;
      else {
         // edge 3, scorro sulle patch in V
         for ( int i = 0 ; i < m_nSpanV ; ++i) {
            PtrOwner<ICurveBezier> pCrvBz3( CreateBasicCurveBezier()) ;
            if ( IsNull(pCrvBz3)  ||  ! pCrvBz3->Init( m_nDegV, m_bRat))
               return nullptr ;
            for ( int p = 0 ; p < m_nDegV + 1 ; ++p ) {
               int nIndex = ( m_nDegU * m_nSpanU + 1) * ( i + p + 1) - 1 ;
               if ( ! m_bRat)
                  pCrvBz3->SetControlPoint( p, GetControlPoint( nIndex, nullptr)) ;
               else
                  pCrvBz3->SetControlPoint( p, GetControlPoint( nIndex, nullptr), GetControlWeight( nIndex, nullptr)) ;
            }
            pCrvCompo->AddCurve( Release( pCrvBz3)) ;
         }
      }
      break ;
   }
   }
   return pCrvCompo ;
}

//----------------------------------------------------------------------------
bool
SurfBezier::IsPlanar( void) const
{
   // costruisco il contorno della superficie unendo gli edge e chiedo se la polyline è piana.
   PtrOwner<ICurveComposite> pCCEdge( GetSingleEdge3D( false, 0)) ;
   pCCEdge->AddCurve( GetSingleEdge3D( false, 1)) ;
   pCCEdge->AddCurve( GetSingleEdge3D( false, 2)) ;
   pCCEdge->AddCurve( GetSingleEdge3D( false, 3)) ;
   PolyLine plApprox ;
   pCCEdge->ApproxWithLines( 0.01, 15, 0, plApprox) ;
   Plane3d plPlane ;
   if ( ! plApprox.IsFlat( plPlane, 2 * EPS_SMALL))
      return false ;
   // in questo caso se è grado 1 in U e V e ho un unica Patch allora sono sicuro sia piana
   if ( m_nDegU == 1  &&  m_nSpanU == 1  &&  m_nDegV == 1  &&  m_nSpanV == 1)  // questa condizione da sola non è sufficiente ( posso avere superfici torte anche se i lati sono segmenti)
      return true ;

   double dULast ; plApprox.GetLastU( dULast) ;
   ++ dULast ;
   // altrimenti devo verificare anche all'interno della superficie, prendendo dei punti campione
   DBLVECTOR vSampling { 0.2, 0.4, 0.6, 0.8} ;
   for ( double i : vSampling) {
      for ( double j : vSampling) {
         Point3d ptBez ;
         GetPointD1D2( i * m_nSpanU, j * m_nSpanV, ISurfBezier::FROM_MINUS, ISurfBezier::FROM_MINUS, ptBez) ;
         if ( plApprox.AddUPoint( dULast, ptBez))
            ++ dULast ;
      }
   }
   plPlane.Reset() ;
   if ( plApprox.IsFlat( plPlane, 2 * EPS_SMALL))
      return true ;

   // nel dubbio restituisco false
   return false ;
}