//----------------------------------------------------------------------------
//  EgalTech 2023-2023
//----------------------------------------------------------------------------
// File : SurfAux.cpp      Data : 09.08.23          Versione :
// Contenuto : Implementazione di alcune funzioni di utilit? per le Superfici.
//
//
//
// Modifiche : 09.08.23 DB Creazione modulo.
//
//
//----------------------------------------------------------------------------

//--------------------------- Include ----------------------------------------
#include "stdafx.h"
#include "CurveAux.h"
#include "GeoConst.h"
#include "CurveLine.h"
#include "CurveArc.h"
#include "CurveBezier.h"
#include "CurveComposite.h"
#include "SurfBezier.h"
#include "/EgtDev/Include/EGkSurf.h"
#include "/EgtDev/Include/EGkSurfAux.h"
#include "/EgtDev/Include/EGkSfrCreate.h"
#include "/EgtDev/Include/EgtPointerOwner.h"
#include "/EgtDev/Include/EGkChainCurves.h"


#define SAVEMKUNIF_CRVS 0
#if SAVEMKUNIF_CRVS
   #include "/EgtDev/Include/EGkGeoObjSave.h"
#endif

using namespace std ;

//----------------------------------------------------------------------------
bool
NurbsSurfaceCanonicalize( SNurbsSurfData& snData)
{
   // per rendere una superficie non periodica devo recuperare i punti di controllo, suddivisi in isoparametriche lungo una direzione
   // e applicare la trasformazione ad ogni isoparametrica, sostituendola poi a quella originale. ( si mantiene il numero di punti)
   if ( snData.bPeriodicU  ||  ! snData.bClampedU) {
      bool bIsRational = snData.bRat ;
      // vettore dei nodi
      DBLVECTOR vU ;
      int nKnot = ssize( snData.vU) ;
      for ( int k = 0 ; k < nKnot ; ++k ) {
         double dKnot = snData.vU[k] ;
         vU.push_back( dKnot) ;
      }
      for ( int j = 0 ; j < snData.nCPV ; ++j) {
         CNurbsData nuCurve ;
         nuCurve.bPeriodic = true ;
         nuCurve.bRat = snData.bRat ;
         nuCurve.nDeg = snData.nDegU ;
         nuCurve.vU = vU ;
         // vettore dei punti di controllo
         PNTVECTOR vPtCtrl ;
         // vettore dei pesi
         DBLVECTOR vWeCtrl ;
         for ( int i = 0 ; i < snData.nCPU ; ++i ) {
            if ( bIsRational) {
               vPtCtrl.push_back( snData.mCP[i][j] / snData.mW[i][j]) ;
               vWeCtrl.push_back( snData.mW[i][j]) ;
            }
            else
               vPtCtrl.push_back( snData.mCP[i][j]) ;
         }
         nuCurve.vCP = vPtCtrl ;
         nuCurve.vW = vWeCtrl ;
         // i punti dell' oggetto nuCurve devono essere in forma non omogenea
         if ( NurbsCurveCanonicalize( nuCurve)) {  // se NurbsCurveCanonicalize ha restituito false (la curva potrebbe esserre un punto di polo) allora non modifico i punti e il vettore dei nodi della superficie
            if ( snData.mCP.size() != nuCurve.vCP.size() ) {
               snData.mCP.resize( nuCurve.vCP.size()) ;
               if ( snData.bRat)
                  snData.mW.resize( nuCurve.vW.size()) ;
            }
            for ( int i = 0 ; i < ssize( nuCurve.vCP) ; ++i) {
               if ( snData.mCP[i].empty())
                  snData.mCP[i].resize( snData.nCPV) ;
               snData.mCP[i][j] = nuCurve.vCP[i] ;
               if ( snData.bRat) {
                  if ( snData.mW[i].empty())
                     snData.mW[i].resize( snData.nCPV) ;
                  snData.mW[i][j] = nuCurve.vW[i] ;
                  snData.mCP[i][j] *= nuCurve.vW[i] ;
               }
            }
            snData.vU = nuCurve.vU ;
         }
      }
      snData.bPeriodicU = false ;
      snData.nCPU = ssize( snData.mCP) ;
   }
   if ( snData.bPeriodicV  ||  ! snData.bClampedV) {
      bool bIsRational = snData.bRat ;
      // vettore dei nodi
      DBLVECTOR vV ;
      int nKnot = ssize( snData.vV) ;
      for ( int k = 0 ; k < nKnot ; ++k ) {
         double dKnot = snData.vV[k] ;
         vV.push_back( dKnot) ;
      }
      for ( int i = 0 ; i < snData.nCPU ; ++i) {
         CNurbsData nuCurve ;
         nuCurve.bPeriodic = true ;
         nuCurve.bRat = snData.bRat ;
         nuCurve.nDeg = snData.nDegV ;
         nuCurve.vU = vV ;
         // vettore dei punti di controllo
         PNTVECTOR vPtCtrl ;
         // vettore dei pesi
         DBLVECTOR vWeCtrl ;
         for ( int j = 0 ; j < snData.nCPV ; ++j ) {
            if ( bIsRational) {
               vPtCtrl.push_back( snData.mCP[i][j] / snData.mW[i][j]) ;
               vWeCtrl.push_back( snData.mW[i][j]) ;
            }
            else
               vPtCtrl.push_back( snData.mCP[i][j]) ;
         }
         nuCurve.vCP = vPtCtrl ;
         nuCurve.vW = vWeCtrl ;
         // i punti dell' oggetto nuCurve devono essere in forma non omogenea
         if ( NurbsCurveCanonicalize( nuCurve)) { // se NurbsCurveCanonicalize ha restituito false (la curva potrebbe esserre un punto di polo) allora non modifico i punti e il vettore dei nodi della superficie
            if ( snData.mCP[i].size() != nuCurve.vCP.size()){
               snData.mCP[i].clear() ;
               snData.mCP[i].resize( nuCurve.vCP.size()) ;
               if ( snData.bRat ) {
                  snData.mW[i].clear() ;
                  snData.mW[i].resize( nuCurve.vW.size()) ;
               }
            }
            for ( int j = 0 ; j < int( nuCurve.vCP.size()) ; ++j ) {
               snData.mCP[i][j] = nuCurve.vCP[j] ;
               if ( snData.bRat ) {
                  snData.mW[i][j] = nuCurve.vW[j] ;
                  snData.mCP[i][j] *= nuCurve.vW[j] ;
               }
            }
            snData.vV = nuCurve.vU ;
         }
      }
      snData.bPeriodicV = false ;
      snData.nCPV = ssize( snData.mCP[0]) ;
   }
   return true ;
}

//----------------------------------------------------------------------------
ISurf*
NurbsToBezierSurface(const SNurbsSurfData& snData)
{
  // la superficie Nurbs deve essere in forma canonica
   if ( snData.bPeriodicU || snData.bPeriodicV || snData.bExtraKnotes )
      return nullptr ;
  // controllo sul numero dei nodi
   int nU = snData.nCPU + snData.nDegU - 1 ;
   int nV = snData.nCPV + snData.nDegV - 1 ;
  // controllo nodi e punti di controllo
   if ( nU != int(snData.vU.size()) || nV != int(snData.vV.size())) {
      return nullptr ;
   }

  // verifico le condizioni agli estremi sui nodi (i primi nDeg nodi e gli ultimi nDeg nodi devono essere uguali tra loro)
   bool bOk = true ;
  // direzione U
   for ( int i = 1 ; i < snData.nDegU ; ++ i) {
      if ( abs( snData.vU[i] - snData.vU[0]) >= EPS_ZERO)
         bOk = false ;
   }
   for ( int i = 1 ; i < snData.nDegU ; ++ i) {
      if ( abs( snData.vU[nU - 1 - i] - snData.vU[nU - 1]) >= EPS_ZERO)
         bOk = false ;
   }
  // direzione V
   for ( int i = 1 ; i < snData.nDegV ; ++ i) {
      if ( abs( snData.vV[i] - snData.vV[0]) >= EPS_ZERO)
         bOk = false ;
   }
   for ( int i = 1 ; i < snData.nDegV ; ++ i) {
      if ( abs( snData.vV[nV - 1 - i] - snData.vV[nV - 1]) >= EPS_ZERO)
         bOk = false ;
   }
   if ( ! bOk)
      return nullptr ;

  // algoritmo 5.7 del libro "The NURBS book"
  // creazione delle strips nella direzione U ( trasformo le curve iso con U costante in bezier)
   int a = snData.nDegU - 1 ;
   int b = snData.nDegU ;
   int nb = 0 ;  // numero di strisce in U ( lunghezza con U costante)
   vector<Point3d> vCPV( snData.nCPV) ;
   vector< vector<Point3d>> mBC (snData.nDegU + 1,vCPV ) ;
   vector< vector<Point3d>> mBC_next (snData.nDegU - 1, vCPV) ;
   vector< vector<Point3d>> mPC_strip(snData.nDegU + 1, vCPV) ; // matrice che verrà ingrandita e conterrà la superficie metà bezier e metà NURBS
   DBLVECTOR vV_W( snData.nCPV) ;
   vector<DBLVECTOR> mW( snData.nDegU + 1, vV_W) ;
   vector<DBLVECTOR> mW_next( snData.nDegU - 1, vV_W) ;
   vector<DBLVECTOR> mW_strip( snData.nDegU + 1, vV_W) ;
   DBLVECTOR vAlpha ;
   vAlpha.resize( snData.nDegU - 1) ;
   if ( ! snData.bRat ) {
      for ( int i = 0 ; i <= snData.nDegU ; ++i) {
         for ( int row = 0 ; row < snData.nCPV ; ++row ) {
            mBC[i][row] = snData.mCP[i][row] ;
         }
      }
   }
   else {
      for ( int i = 0 ; i <= snData.nDegU ; ++i) {
         for (int row = 0 ; row < snData.nCPV ; ++ row) {
            mW[i][row] = snData.mW[i][row] ;
            mBC[i][row] = snData.mCP[i][row] * snData.mW[i][row] ;
         }
      }
   }

   bool bRef = false ;
   while ( snData.nDegU != 1 ? b < nU - 1 : b < nU) { // qui correggo un probabile errore, mettendo nU anziché nCPV, come indicato nell'algoritmo
      int i = b ;
      while ( b < nU - 1  &&  abs( snData.vU[b+1] - snData.vU[b]) < EPS_ZERO)
         ++ b ;
      int mult = b - i + 1 ;
      if ( mult < snData.nDegU) {
         bRef = true ;
         // calcolo numeratore e alpha
         double numer = snData.vU[b] - snData.vU[a] ;
         for ( int j = snData.nDegU ; j > mult ; -- j)
            vAlpha[j-mult-1] = numer / ( snData.vU[a+j] - snData.vU[a]) ;
         int r = snData.nDegU - mult ;
         for ( int j = 1 ; j <= snData.nDegU - mult ; ++j) {
            int save = r - j ;
            int s = mult + j ;
            if ( ! snData.bRat ) {
               for ( int k = snData.nDegU ; k >= s ; --k) {
                  for ( int row = 0 ; row < snData.nCPV ; ++row) {
                     mBC[k][row] = vAlpha[k-s] * mBC[k][row] + ( 1 -  vAlpha[k-s]) * mBC[k-1][row] ;
                  }
               }
            }
            else {
               for ( int k = snData.nDegU ; k >= s ; --k) {
                  for ( int row = 0 ; row < snData.nCPV ; ++row) {
                     mBC[k][row] = vAlpha[k-s] * mBC[k][row] + ( 1 -  vAlpha[k-s]) * mBC[k-1][row] ;
                     mW[k][row] = vAlpha[k-s] * mW[k][row] + ( 1 - vAlpha[k-s]) * mW[k-1][row] ;
                  }
               }
            }

            if ( b < nU - 1 ) {
               for ( int row = 0 ; row < snData.nCPV ; ++row) {
                  mBC_next[save][row] =  mBC[snData.nDegU][row] ;
               }
               if ( snData.bRat )
                  for ( int row = 0 ; row < snData.nCPV ; ++row) { 
                     mW_next[save][row] = mW[snData.nDegU][row] ;
                  }
            }
         }
      }
      mPC_strip.resize( snData.nDegU * ( nb  + 1) + 1, vCPV) ;
      mW_strip.resize( snData.nDegU * ( nb  + 1) + 1, vV_W) ;
      if ( ! snData.bRat)
         for ( int i = 0 ; i <= snData.nDegU ; ++i) {
            for ( int row = 0 ; row < snData.nCPV ; ++row) {
               mPC_strip[i+ nb * snData.nDegU][row] = mBC[i][row] ;
            }
         }
      else {
         for ( int i = 0 ; i <= snData.nDegU ; ++i) {
            for ( int row = 0 ; row < snData.nCPV ; ++row) {
               mPC_strip[i+ nb * snData.nDegU][row] = mBC[i][row]/mW[i][row] ;
               mW_strip[i+ nb * snData.nDegU][row] = mW[i][row] ;
            }
         }
      }
     // ho finito di definire la patch di Bezier attuale e passo alla successiva
      ++ nb ;

     // aggiorno mBC con i valori della prossima pezza di Bezier // corrisponde a nb = nb + 1 
      if ( ! snData.bRat){
         for (int i = 0 ; i < snData.nDegU - 1 ; ++ i) {
            for ( int row = 0 ; row < snData.nCPV ; ++row) {
               mBC[i][row] = mBC_next[i][row] ;
            }
         }
      }
      else {
         for (int i = 0 ; i < snData.nDegU - 1 ; ++ i) {
            for ( int row = 0 ; row < snData.nCPV ; ++row) {
               mBC[i][row] = mBC_next[i][row] ;
               mW[i][row] = mW_next[i][row] ;
            }
         }
      }

      if ( b < nU - 1 ) {
         if ( ! snData.bRat ) {
            for ( int i = snData.nDegU - mult ; i <= snData.nDegU ; ++ i) {
               for (int row = 0 ; row < snData.nCPV ; ++ row ) {
                  mBC[i][row] = snData.mCP[b - snData.nDegU + i + 1][row] ;
               }
            }
         }
         else {
            for ( int i = snData.nDegU - mult ; i <= snData.nDegU ; ++ i) {
               for (int row = 0 ; row < snData.nCPV ; ++ row ) {
                  mBC[i][row] = snData.mCP[b - snData.nDegU + i + 1][row] * snData.mW[b - snData.nDegU + i + 1][row] ;
                  mW[i][row] = snData.mW[b - snData.nDegU + i + 1][row] ;
               }
            }
         }
         a = b ;
         ++b ;
      }
      else if ( snData.nDegU == 1  &&  b < nU) {
         a = b ;
         ++b ;
      }
   }

  // se non ho raffinato allora tutti i nodi avevano gi? molteplicit? massima. Converto direttamente in Bezier la dir U
   int nCPU_ref ; // numero dei punti di controllo in U dopo il raffinamento
   if ( ! bRef ) {
      nCPU_ref = snData.nCPU ;
      mPC_strip.resize( snData.nCPU, vCPV) ;
      mW_strip.resize( snData.nCPU, vV_W) ;
      if ( ! snData.bRat) {
         for ( int i = 0 ; i < nCPU_ref ; ++i) {
            for ( int row = 0 ; row < snData.nCPV ; ++ row) {
               mPC_strip[i][row] = snData.mCP[i][row] ;
            }
         }
      }
      else {
         for ( int i = 0 ; i < nCPU_ref ; ++i) {
            for ( int row = 0 ; row < snData.nCPV ; ++ row) {
               mPC_strip[i][row] = snData.mCP[i][row] ;
               mW_strip[i][row] = snData.mW[i][row] ;
            }
         }
      }
   }
   else
      nCPU_ref = snData.nDegU * nb + 1 ; // numero dei punti di controllo in U dopo il raffinamento

  // ora ho ottenuto le strisce nDegU x nCPV
  // devo ripetere la procedura, sulla dir V, per ottenere le patch nDegU x nDegV
   a = snData.nDegV - 1 ;
   b = snData.nDegV ;
   int nc = 0 ;  // numero di strisce in V ( lunghezza con V costante)
   vector<Point3d> vDegV(snData.nDegV + 1) ;
   vector<Point3d> vDegV_1(snData.nDegV - 1) ;
   vector< vector<Point3d>> m_BC1( nCPU_ref, vDegV) ;
   vector< vector<Point3d>> m_BC1_next( nCPU_ref, vDegV_1) ;
   DBLVECTOR vV1_W(snData.nDegV + 1) ;
   DBLVECTOR vV2_W(snData.nDegV - 1) ;
   vector<DBLVECTOR> mW1( nCPU_ref, vV1_W) ;
   vector<DBLVECTOR> mW1_next( nCPU_ref, vV2_W) ;
   DBLVECTOR vAlpha1( snData.nDegV - 1) ;
   vector<vector<Point3d>> mPC_tot( nCPU_ref, vDegV) ;
   vector<DBLVECTOR> mW_tot( nCPU_ref, vV1_W) ;
   if ( ! snData.bRat ) {
      for ( int i = 0 ; i < nCPU_ref ; ++i) {
         for ( int row = 0 ; row <= snData.nDegV ; ++row ) {
            m_BC1[i][row] = mPC_strip[i][row] ;
         }
      }
   }
   else {
      for ( int i = 0 ; i < nCPU_ref ; ++i) {
         for (int row = 0 ; row <= snData.nDegV ; ++ row) {
            mW1[i][row] = mW_strip[i][row] ;
            m_BC1[i][row] = mPC_strip[i][row] * mW_strip[i][row] ;
         }
      }
   }

   bRef = false ;
   while ( snData.nDegV != 1 ? b < nV - 1 : b < nV) { // qui correggo un probabile errore, mettendo nU anziché nCPV, come indicato nell'algoritmo
      int i = b ;
      while ( b < nV - 1  &&  abs( snData.vV[b+1] - snData.vV[b]) < EPS_ZERO)
         ++ b ;
      int mult = b - i + 1 ;
      if ( mult < snData.nDegV ) {
         bRef = true ;
        // calcolo numeratore e alpha
         double numer = snData.vV[b] - snData.vV[a] ;
         for ( int j = snData.nDegV ; j > mult ; -- j)
            vAlpha1[j-mult-1] = numer / ( snData.vV[a+j] - snData.vV[a]) ;
         int r = snData.nDegV - mult ;
         for ( int j = 1 ; j <= snData.nDegV - mult ; ++j) {
            int save = r - j ;
            int s = mult + j ;
            if ( ! snData.bRat) {
               for ( int k = 0 ; k < nCPU_ref ; ++k) {
                  for ( int row = snData.nDegV ; row >= s ; --row) {
                     m_BC1[k][row] = vAlpha1[row-s] * m_BC1[k][row] + ( 1 -  vAlpha1[row-s]) * m_BC1[k][row-1] ;
                  }
               }
            }
            else {
               for ( int k = 0 ; k < nCPU_ref ; ++k) {
                  for ( int row = snData.nDegV ; row >= s ; --row) {
                     m_BC1[k][row] = vAlpha1[row-s] * m_BC1[k][row] + ( 1 -  vAlpha1[row-s]) * m_BC1[k][row-1] ;
                     mW1[k][row] = vAlpha1[row-s] * mW1[k][row] + ( 1 -  vAlpha1[row-s]) * mW1[k][row-1] ;
                  }
               }
            }

            if ( b < nV - 1 ) {
               if ( !snData.bRat ){
                  for ( int i = 0 ; i < nCPU_ref ; ++i) {
                     m_BC1_next[i][save] =  m_BC1[i][snData.nDegV] ;
                  }
               }
               else {
                  for ( int i = 0 ; i < nCPU_ref ; ++i) {
                     m_BC1_next[i][save] =  m_BC1[i][snData.nDegV] ;
                     mW1_next[i][save] = mW1[i][snData.nDegV] ;
                  }
               }
            }
         }
      }
      int nRef = snData.nDegV * ( nc + 1) + 1 ;
      for ( int k = 0 ; k < nCPU_ref; ++k) {
         mPC_tot[k].resize( nRef) ;
         mW_tot[k].resize( nRef) ;
      }
      if ( ! snData.bRat)
         for ( int i = 0 ; i < nCPU_ref ; ++i) {
            for ( int row = 0 ; row <= snData.nDegV ; ++row) {
               mPC_tot[i][row + nc * snData.nDegV] = m_BC1[i][row] ;
            }
         }
      else {
         for ( int i = 0 ; i < nCPU_ref ; ++i) {
            for ( int row = 0 ; row <= snData.nDegV ; ++row) {
               mPC_tot[i][row + nc * snData.nDegV] = m_BC1[i][row] / mW1[i][row] ;
               mW_tot[i][row + nc * snData.nDegV] = mW1[i][row] ;
            }
         }
      }
     // ho finito di definire la patch di Bezier attuale e passo alla successiva
      ++ nc ;

     // aggiorno mBC con i valori della prossima pezza di Bezier // corrisponde a nc = nc + 1 
      if ( ! snData.bRat){
         for (int i = 0 ; i < nCPU_ref ; ++ i) {
            for ( int row = 0 ; row < snData.nDegV - 1 ; ++row) {
               m_BC1[i][row] = m_BC1_next[i][row] ;
            }
         }
      }
      else {
         for (int i = 0 ; i < nCPU_ref ; ++ i) {
            for ( int row = 0 ; row < snData.nDegV - 1 ; ++row) {
               m_BC1[i][row] = m_BC1_next[i][row] ;
               mW1[i][row] = mW1_next[i][row] ;
            }
         }
      }

      if ( b < nV - 1) {
         if ( ! snData.bRat) {
            for (int i = 0 ; i < nCPU_ref ; ++ i ) {
               for ( int row = snData.nDegV - mult ; row <= snData.nDegV ; ++ row) {
                  m_BC1[i][row] = mPC_strip[i][b - snData.nDegV + row + 1] ;
               }
            }
         }
         else {
            for (int i = 0 ; i < nCPU_ref ; ++ i ) {
               for ( int row = snData.nDegV - mult ; row <= snData.nDegV ; ++ row) {
                  m_BC1[i][row] = mPC_strip[i][b - snData.nDegV + row + 1] * mW_strip[i][b - snData.nDegV + row + 1] ;
                  mW1[i][row] = mW_strip[i][b - snData.nDegV + row + 1] ;
               }
            }
         }
         a = b ;
         ++b ;
      }
      else if ( snData.nDegV == 1  &&  b < nV) {
         a = b ;
         ++b ;
      }
   }
  // se non ho raffinato allora aggiungo direttamente alle matrici della superficie totale
   int nCPV_ref ; // numero dei punti di controllo in V dopo il raffinamento
   if ( ! bRef) {
      nCPV_ref = snData.nCPV ;
      for ( int k = 0 ; k < nCPU_ref ; ++k){
         mPC_tot[k].resize( snData.nCPV) ;
         mW_tot[k].resize( snData.nCPV) ;
      }
      if ( ! snData.bRat) {
         for ( int i = 0 ; i < nCPU_ref ; ++i) {
            for ( int row = 0 ; row < nCPV_ref ; ++ row) {
               mPC_tot[i][row] = mPC_strip[i][row] ;
            }
         }
      }
      else {
         for ( int i = 0 ; i < nCPU_ref ; ++i) {
            for ( int row = 0 ; row < nCPV_ref ; ++ row) {
               mPC_tot[i][row] = mPC_strip[i][row] ;
               mW_tot[i][row] = mW_strip[i][row] ;
            }
         }
      }
   }
   else
      nCPV_ref = snData.nDegV * nc + 1 ;

  // finalmente setto la superficie di bezier totale divisa in nb patch in U e nc patch in V
   PtrOwner<SurfBezier> pSrfBz( CreateBasicSurfBezier()) ;
   if ( IsNull( pSrfBz))
      return nullptr ;
   pSrfBz->Init(snData.nDegU, snData.nDegV, nb, nc, snData.bRat) ;
   if ( !snData.bRat ) {
      for ( int i = 0 ; i < nCPU_ref; ++ i) {
         for (int j = 0 ; j < nCPV_ref; ++j) {
            pSrfBz->SetControlPoint( i + nCPU_ref * j, mPC_tot[i][j]) ;
         }
      }
   }
   else {
      for ( int i = 0 ; i < nCPU_ref; ++ i) {
         for (int j = 0 ; j < nCPV_ref; ++j) {
            pSrfBz->SetControlPoint( i + nCPU_ref * j, mPC_tot[i][j] /*/ mW_tot[i][j]*/, mW_tot[i][j]) ;
         }
      }
   }
   return Release( pSrfBz) ;
}

//----------------------------------------------------------------------------
ICurveComposite*
GetRectangleCurve( const Point3d& ptStart, double dWidth, double dHeight)
{
   PolyLine PL ;
   PL.AddUPoint( 0, ptStart) ;
   PL.AddUPoint( 1, Point3d( ptStart.x + dWidth, ptStart.y)) ;
   PL.AddUPoint( 2, Point3d( ptStart.x + dWidth, ptStart.y + dHeight, 0)) ;
   PL.AddUPoint( 3, Point3d( ptStart.x , ptStart.y + dHeight, 0)) ;
   PL.AddUPoint( 4, ptStart) ;
   // creo la curva e la inserisco nel GDB
   PtrOwner<ICurveComposite> pCrvCompo( CreateCurveComposite()) ;
   bool bOk = true ;
   bOk = bOk && ! IsNull( pCrvCompo) ;
   // inserisco i segmenti che uniscono i punti
   bOk = bOk && pCrvCompo->FromPolyLine( PL) ;
   return Release( pCrvCompo) ;
}

//----------------------------------------------------------------------------
bool
MakeUniform( ISurfFlatRegion*& pSfr, bool& bRescaled, const DBLVECTOR& vU0, const DBLVECTOR& vV0,
             int nDegU, int nDegV, double dScaleU, double dScaleV, bool bRetry) 
{
  // la superficie in input arriva già scalata
   bool bRescaledU = false ;
   bool bRescaledV = false ;
   int nSpanU = 1, nSpanV = 1 ;
   ICRVCOMPOPOVECTOR vUniformedCurves ;
   BOOLVECTOR vbUniform(2) ;
   fill( vbUniform.begin(), vbUniform.end(), true) ;
   DBLMATRIX mKnots(2) ;
   for ( int nDir = 0 ; nDir <= 1 ; ++ nDir) {
     // vettore dei nodi
      DBLVECTOR& vU = mKnots[nDir] ;
      int nExtraKnots = 0 ;
     // controllo in U
      if ( nDir == 0) {
         for ( int i = nExtraKnots ; i < ssize( vU0) - nExtraKnots ; ++i ) {
            double dKnot = vU0[i] * SBZ_TREG_COEFF ;
            // lo aggiungo solo se è diverso dal precedente
            if ( i == nExtraKnots  ||  dKnot > vU.back() + EPS_SMALL)
               vU.push_back( dKnot) ;
         }
         nSpanU = ssize( vU) - 1 ;
      }
     // controllo in V
      else if ( nDir == 1 ) {
         for ( int i = nExtraKnots ; i < ssize( vV0) - nExtraKnots ; ++i ) {
            double dKnot = vV0[i] * SBZ_TREG_COEFF ;
            // lo aggiungo solo se è diverso dal precedente
            if ( i == nExtraKnots  ||  dKnot > vU.back() + EPS_SMALL)
               vU.push_back( dKnot) ;
         }
         nSpanV = ssize( vU) - 1 ;
      }

     // controllo se il vettore dei nodi è uniforme
      int a = 0, b = 1 ;
      double d0 = abs( vU[b] - vU[a]), d1 = d0 ;
      // il vettore è uniforme quando la distanza tra nodi consecutivi è sempre zero o un valore costante
      while ( b < (int)vU.size()  &&  ( ( d1 < d0 + EPS_SMALL  &&  d1 > d0 - EPS_SMALL)  ||  d1 < EPS_SMALL)) {
         a = b ;
         ++b ;
         if ( b < (int)vU.size())
            d1 = abs( vU[b] - vU[a]) ;
      }
      if ( b != (int)vU.size())
         vbUniform[nDir] = false ;
   }
   // vettore delle curve di loop della regione di trim
   ICRVCOMPOPOVECTOR vLoop ;
   if ( ! vbUniform[0]  ||  ! vbUniform[1]) {
      for ( int c = 0 ;  c < pSfr->GetChunkCount() ; ++c) {
         for ( int l = 0 ; l < pSfr->GetLoopCount( c); ++l)
            vLoop.emplace_back( ConvertCurveToComposite( pSfr->GetLoop( c, l))) ;
      }
   }

#if SAVEMKUNIF_CRVS
   //debug
   vector<IGeoObj*> vGeo ;
   for( int i = 0 ; i < ssize( vLoop); ++i){
      vGeo.push_back(vLoop[i]->Clone()) ;
   }
   SaveGeoObj( vGeo, "D:\\Temp\\bezier\\import3dm\\trim_error\\failed_trim_crv_loops.nge") ;
   //debug
#endif

   for ( int nDir = 0 ; nDir <= 1 ; ++ nDir) {
      DBLVECTOR& vU = mKnots[nDir] ;
      if ( ! vbUniform[nDir]) {
         nDir == 0 ? bRescaledU = true : bRescaledV = true ;
         // creo il vettore delle curve all'interno di una striscia
         ICRVCOMPOPOVECTOR vCrvStrip ; 
         for ( int p = 0 ; p < (int)vU.size() - 1 ; ++p) {
            double dLenStrip = abs( vU[p+1] - vU[p]) ;
            if ( dLenStrip < EPS_SMALL)
               continue ;
            // creo la maschera per tagliare la superficie originale e ottenere una striscia
            PtrOwner<ICurveComposite> pTrimMask ;

            // ricavo la maschera del trim, con cui poi farò l'intersezione con la sfr iniziale
            if ( nDir == 0) {
               Point3d ptStart( abs(vU[p] - vU.front()), - 1, 0) ;
               pTrimMask.Set( GetRectangleCurve( ptStart, dLenStrip, dScaleV + 2)) ;
            }
            else{
               Point3d ptStart( - 1, abs(vU[p] - vU.front()), 0) ;
               pTrimMask.Set( GetRectangleCurve( ptStart, dScaleU + 2, dLenStrip)) ;
            }

            // qui potrei decidere di eliminare tutti i tratti dei loop paralleli alla direzione del parametro che sto analizzando ( e di valore coincidente a quello di un nodo)

            for ( int l = 0 ; l < ssize( vLoop); ++l) {
               #if SAVEMKUNIF_CRVS
                  //debug
                  vector<IGeoObj*> vGeo ;
                  vGeo.push_back(pTrimMask->Clone()) ;
                  vGeo.push_back(vLoop[l]->Clone()) ;
                  SaveGeoObj( vGeo, "D:\\Temp\\bezier\\import3dm\\trim_error\\failed_trim_crv_inters.nge") ;
                  //debug
               #endif
               IntersCurveCurve icc( *pTrimMask, *vLoop[l]) ;
               CRVCVECTOR vCurveClass ;
               icc.GetCurveClassification( 1, 10 * EPS_SMALL, vCurveClass) ;
               for( int i = 0 ; i < ssize( vCurveClass); ++i) {
                  if( vCurveClass[i].nClass == CRVC_IN  ||  vCurveClass[i].nClass == CRVC_ON_P)
                     vCrvStrip.emplace_back( ConvertCurveToComposite( vLoop[l]->CopyParamRange( vCurveClass[i].dParS, vCurveClass[i].dParE))) ;
               }
            }

            #if SAVEMKUNIF_CRVS
               //debug
               vector<IGeoObj*> vGeo ;
               for( int i = 0 ; i < ssize( vCrvStrip); ++i){
                  vGeo.push_back(vCrvStrip[i]->Clone()) ;
               }
               SaveGeoObj( vGeo, "D:\\Temp\\bezier\\import3dm\\trim_error\\failed_trim_crv_strip.nge") ;
               //debug
            #endif

            // riscalo le curve nella striscia
            if ( ! vCrvStrip.empty()) {
               double dCoeffX = 1, dCoeffY = 1 ;
               
               if ( nDir == 0)
                  dCoeffX = SBZ_TREG_COEFF / dLenStrip ;
               else 
                  dCoeffY = SBZ_TREG_COEFF / dLenStrip ;
               for( int i = 0 ; i < ssize( vCrvStrip); ++i) {
                  if( ! IsNull( vCrvStrip[i]))
                     vCrvStrip[i]->Scale( GLOB_FRM, dCoeffX, dCoeffY, 1) ;
                  else
                     return false ;
               }

               // prima di riunire le curve al resto devo traslarle sul bordo destro della superificie che sto ricostruendo
               Point3d pt ;
               nDir == 0 ? pt.Set( abs(vU[p] - vU.front()), 0, 0) : pt.Set( 0,abs(vU[p] - vU.front()), 0) ;
               if ( nDir == 0)
                  pt.Scale( GLOB_FRM, SBZ_TREG_COEFF / dLenStrip, 1, 1) ;
               else
                  pt.Scale( GLOB_FRM, 1, SBZ_TREG_COEFF / dLenStrip, 1) ;

               Vector3d vtJoin ;
               if ( nDir == 0)
                  vtJoin.Set( p * SBZ_TREG_COEFF - pt.x, 0, 0) ;
               else
                  vtJoin.Set( 0, p  * SBZ_TREG_COEFF - pt.y, 0) ;
               for( int i = 0 ; i < ssize( vCrvStrip); ++i)
                  vCrvStrip[i]->Translate( vtJoin) ;

            #if SAVEMKUNIF_CRVS
               //debug
               vector<IGeoObj*> vGeo ;
               for( int i = 0 ; i < ssize( vCrvStrip); ++i){
                  vGeo.push_back(vCrvStrip[i]->Clone()) ;
               }
               SaveGeoObj( vGeo, "D:\\Temp\\bezier\\import3dm\\trim_error\\failed_trim_crv_strip.nge") ;
               //debug
            #endif
               
               // faccio la chain con le curve delle striscie precedenti
               
               if ( ! vUniformedCurves.empty()  ||  ! vCrvStrip.empty()) {
                  #if SAVEMKUNIF_CRVS
                     //debug
                     vector<IGeoObj*> vGeo ;
                     for( int i = 0 ; i < ssize( vUniformedCurves); ++i){
                        vGeo.push_back(vUniformedCurves[i]->Clone()) ;
                     }
                     SaveGeoObj( vGeo, "D:\\Temp\\bezier\\import3dm\\trim_error\\failed_trim_crv_unif.nge") ;
                     //debug
                  #endif

                  ChainCurves chainCrv ;
                  double dChainTol = 5 * EPS_SMALL ;
                  chainCrv.Init( false, dChainTol, max(ssize( vUniformedCurves), ssize(vCrvStrip))) ;
                  for ( int c = 0 ; c < ssize( vUniformedCurves) ; ++c) {
                     Point3d ptStart, ptEnd ;
                     Vector3d vtStart, vtEnd ;
                     vUniformedCurves[c]->GetStartPoint( ptStart) ;
                     vUniformedCurves[c]->GetEndPoint( ptEnd) ;
                     vUniformedCurves[c]->GetStartDir( vtStart) ;
                     vUniformedCurves[c]->GetEndDir( vtEnd) ;
                     chainCrv.AddCurve( 1 + c, ptStart, vtStart, ptEnd, vtEnd) ;
                  }
                  for ( int c = 0 ; c < ssize( vCrvStrip); ++c) {
                     Point3d ptStart, ptEnd ;
                     Vector3d vtStart, vtEnd ;
                     vCrvStrip[c]->GetStartPoint( ptStart) ;
                     vCrvStrip[c]->GetEndPoint( ptEnd) ;
                     vCrvStrip[c]->GetStartDir( vtStart) ;
                     vCrvStrip[c]->GetEndDir( vtEnd) ;
                     chainCrv.AddCurve( 1 + ssize( vUniformedCurves) + c, ptStart, vtStart, ptEnd, vtEnd) ;
                  }
                  INTVECTOR vIds ;
                  ICRVCOMPOPOVECTOR vNewCrv ;
                  int nCrvPrec = ssize( vUniformedCurves) ;
                  while ( chainCrv.GetChainFromNear( ORIG, false, vIds)) {
                     vNewCrv.emplace_back( CreateBasicCurveComposite()) ;
                     for ( int nId : vIds) {
                        nId -= 1 ;
                        if ( nId < nCrvPrec)
                           vNewCrv.back()->AddCurve( Release( vUniformedCurves[nId]), true, dChainTol) ;
                        else
                           vNewCrv.back()->AddCurve( Release( vCrvStrip[nId - ssize( vUniformedCurves)]), true, dChainTol) ;
                     }
                  }
                  // aggiorno le curve 
                  vUniformedCurves.clear() ;
                  vUniformedCurves.swap( vNewCrv)  ;
               }
            }
            vCrvStrip.clear() ;
         }
         if ( nDir == 0) {
            dScaleU = ((int)vU.size() - 1) * SBZ_TREG_COEFF ;
            if( ! vbUniform[1]) {
               vLoop.swap( vUniformedCurves) ;
               vUniformedCurves.clear() ;
            }
         }
         else
            dScaleV = ((int)vU.size() - 1) * SBZ_TREG_COEFF ;
      }
   }

   if ( ! vUniformedCurves.empty()) {
      #if SAVEMKUNIF_CRVS
         //debug
         vector<IGeoObj*> vGeo ;
         for( int i = 0 ; i < ssize( vUniformedCurves); ++i){
            vGeo.push_back(vUniformedCurves[i]->Clone()) ;
         }
         SaveGeoObj( vGeo, "D:\\Temp\\bezier\\import3dm\\trim_error\\failed_trim_crv_unif.nge") ;
         //debug
      #endif

      // controllo che tutte le curve siano chiuse, sennò vuol dire che ho perso qualche pezzo durante le intersezioni
      for ( int i = 0 ; i < ssize( vUniformedCurves); ++i) {
         if ( ! vUniformedCurves[i]->IsClosed())
            return false ;
      }

      // creo una regione dalle curve riscalate
      SurfFlatRegionByContours sfrRescaled ;
      for ( int i = 0 ; i < ssize( vUniformedCurves); ++i)
         sfrRescaled.AddCurve( Release( vUniformedCurves[i])) ;

      PtrOwner<ISurfFlatRegion> pRescaledSfr( sfrRescaled.GetSurf()) ;
      if ( ! IsNull( pRescaledSfr)  &&  pRescaledSfr->IsValid()) {
         delete pSfr ;
         pSfr = Release( pRescaledSfr) ;
      }
      else
         return false ;
   }
   else if( ! vbUniform[0]  ||  ! vbUniform[1])
      return false ;

   if ( ! bRescaledU  &&  ! bRescaledV)
      pSfr->Scale( GLOB_FRM, nSpanU / dScaleU * SBZ_TREG_COEFF, nSpanV / dScaleV * SBZ_TREG_COEFF, 1) ;
   else if ( bRescaledU  &&  ! bRescaledV)
      pSfr->Scale( GLOB_FRM, 1, nSpanV / dScaleV * SBZ_TREG_COEFF, 1) ;
   else if ( ! bRescaledU  &&  bRescaledV)
      pSfr->Scale( GLOB_FRM, nSpanU / dScaleU * SBZ_TREG_COEFF, 1, 1) ;

   bRescaled = ( bRescaledU || bRescaledV) ;
   return true ;
}

//----------------------------------------------------------------------------
bool
OnWhichEdge( double u0, double u1, double v0, double v1, const Point3d& ptToAssign, int& nEdge)
{
   Point3d ptTR( u1, v1) ;
   Point3d ptTl( u0, v1) ;
   Point3d ptBL( u0, v0) ;
   Point3d ptBr( u1, v0) ;

   double dEps = 0.1 ;
   if ( AreSamePointEpsilon( ptToAssign, ptTR, dEps))
      nEdge = 7 ;
   else if ( AreSamePointEpsilon( ptToAssign, ptTl, dEps))
      nEdge = 4 ;
   else if ( AreSamePointEpsilon( ptToAssign, ptBL, dEps))
      nEdge = 5 ;
   else if ( AreSamePointEpsilon( ptToAssign, ptBr, dEps))
      nEdge = 6 ;
   else if ( ptToAssign.x > ptBL.x - dEps  &&  ptToAssign.x < ptTR.x + dEps  &&  abs( ptToAssign.y - ptTR.y) < dEps)
      nEdge = 0 ;
   else if ( ptToAssign.y > ptBL.y - dEps  &&  ptToAssign.y < ptTR.y + dEps  &&  abs( ptToAssign.x - ptBL.x) < dEps)
      nEdge = 1 ;
   else if ( ptToAssign.x > ptBL.x - dEps  &&  ptToAssign.x < ptTR.x + dEps  &&  abs( ptToAssign.y - ptBL.y) < dEps)
      nEdge = 2 ;
   else if ( ptToAssign.y > ptBL.y - dEps  &&  ptToAssign.y < ptTR.y + dEps  &&  abs( ptToAssign.x - ptTR.x) < dEps)
      nEdge = 3 ;
   else
      return false ;
   return true ;
}