//----------------------------------------------------------------------------
//  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 "/EgtDev/Include/EgtPointerOwner.h"
#include "/EgtDev/Include/EGkSurf.h"
#include "/EgtDev/Include/EGkSurfAux.h"
#include "/EgtDev/Include/EGkSurfBezier.h"

using namespace std ;

//----------------------------------------------------------------------------
ISurf*
NurbsToBezierSurface(const CNurbsSurfData& cnData)
{
   //INTVECTOR vInt_sub( 10) ;
   //INTMATRIX vInt( 10, vInt_sub) ;
   //for ( int i = 0 ; i < 10 ; ++i ) {
   //   for ( int j = 0 ; j < 10 ; ++j ) {
   //      vInt[i][j] = i + 10 * j ;
   //   }
   //}
   //vInt_sub.resize( 20) ;
   //vInt[0].resize( 20) ;





   // la superficie Nurbs deve essere in forma canonica
   if ( cnData.bPeriodicU || cnData.bPeriodicV || cnData.bExtraKnotes )
      return nullptr ;
   // controllo sul numero dei nodi
   int nU = cnData.nCPU + cnData.nDegU - 1 ;
   int nV = cnData.nCPV + cnData.nDegV - 1 ;
   // controllo nodi e punti di controllo
   //if ( nU != int( cnData.vU.size())  ||  nV != int( cnData.vV.size())  ||  cnData.nCPU * cnData.nCPV != int( cnData.vCP.size()))
   if ( nU != int(cnData.vU.size()) || nV != int(cnData.vV.size())) {
      return nullptr ;
   }
      
   //// numero degli intervalli
   //int nInt = nU - 2 * cnData.nDeg + 1 ;
   // 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 < cnData.nDegU ; ++ i) {
      if ( abs( cnData.vU[i] - cnData.vU[0]) >= EPS_ZERO)
         bOk = false ;
   }
   for ( int i = 1 ; i < cnData.nDegU ; ++ i) {
      if ( abs( cnData.vU[nU - 1 - i] - cnData.vU[nU - 1]) >= EPS_ZERO)
         bOk = false ;
   }
   // direzione V
   for ( int i = 1 ; i < cnData.nDegV ; ++ i) {
      if ( abs( cnData.vV[i] - cnData.vV[0]) >= EPS_ZERO)
         bOk = false ;
   }
   for ( int i = 1 ; i < cnData.nDegV ; ++ i) {
      if ( abs( cnData.vV[nV - 1 - i] - cnData.vV[nV - 1]) >= EPS_ZERO)
         bOk = false ;
   }
   if ( ! bOk)
      return nullptr ;

   //// se 1 solo intervallo, la Nurbs è già una curva di Bezier
   //if ( nInt == 1) {
   //   // creo la curva di Bezier
   //   PtrOwner<ICurveBezier> pCrvBez( CreateCurveBezier()) ;
   //   if ( IsNull( pCrvBez))
   //      return nullptr ;
   //   // la inizializzo
   //   if ( ! pCrvBez->Init( cnData.nDeg, cnData.bRat))
   //      return nullptr ;
   //   for ( int i = 0 ; i <= cnData.nDeg ; ++ i) {
   //      if ( ! cnData.bRat) {
   //         if ( ! pCrvBez->SetControlPoint( i, cnData.vCP[i]))
   //            return nullptr ;
   //      }
   //      else {
   //         if ( ! pCrvBez->SetControlPoint( i, cnData.vCP[i], cnData.vW[i]))
   //            return nullptr ;
   //      }
   //   }
   //   // se non è una curva ma un punto, la invalido
   //   if ( pCrvBez->IsAPoint())
   //      pCrvBez->Init( cnData.nDeg, cnData.bRat) ;
   //   // restituisco la curva
   //   return Release( pCrvBez) ;
   //}


   // 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 = cnData.nDegU - 1 ;
   int b = cnData.nDegU ;
   int nb = 0 ;  // numero di strisce in U ( lunghezza con U costante)
   //PNTVECTOR vBC ;
   //vBC.resize( cnData.nCPV * cnData.nDegU) ;
   //for (int row = 0 ; row < cnData.nCPV ; ++row ) {
   //   for ( int i = 0 ; i <= cnData.nDegU ; ++i ) {
   //      vBC[nDegU*row + i] = ( cnData.vCP[nCPU*row + i]) ;
   //   }
   //}
   vector<Point3d> vCPV( cnData.nCPV) ;
   vector< vector<Point3d>> mBC (cnData.nDegU + 1,vCPV ) ;
   vector< vector<Point3d>> mBC_next (cnData.nDegU - 1, vCPV) ;
   vector< vector<Point3d>> mPC_strip(cnData.nDegU + 1, vCPV) ; // matrice che verrà ingrandita e conterrà la superficie metà bezier e metà NURBS
   DBLVECTOR vV_W( cnData.nCPV) ;
   vector<DBLVECTOR> mW( cnData.nDegU + 1, vV_W) ;
   vector<DBLVECTOR> mW_next( cnData.nDegU - 1, vV_W) ;
   vector<DBLVECTOR> mW_strip( cnData.nDegU + 1, vV_W) ;
   DBLVECTOR vAlpha ;
   vAlpha.resize( cnData.nDegU - 1) ;
   if ( ! cnData.bRat ) {
      for ( int i = 0 ; i <= cnData.nDegU ; ++i ) {
         for ( int row = 0 ; row < cnData.nCPV ; ++row ) {
            mBC[i][row] = cnData.mCP[i][row] ;
         }
      }
   }
   else {
      for ( int i = 0 ; i <= cnData.nDegU ; ++i ) {
         for (int row = 0 ; row < cnData.nCPV ; ++ row) {
            mW[i][row] = cnData.mW[i][row] ;
            mBC[i][row] = cnData.mCP[i][row] * cnData.mW[i][row] ;
         }
      }
   }
   
   bool bRef = false ;
   // se la superficie è lineare nel parametro U allora è già in forma di Bezier in questo parametro
   if ( b == nU - 1 ) {
      nb = 1 ;
   }
   while ( b < nU - 1) { // qui correggo un probabile errore, mettendo nU anziché nCPV, come indicato nell'algoritmo
      int i = b ;
      while ( b < nU - 1  &&  abs( cnData.vU[b+1] - cnData.vU[b]) < EPS_ZERO)
         ++ b ;
      int mult = b - i + 1 ;
      if ( mult < cnData.nDegU ) {
         bRef = true ;
         // calcolo numeratore e alpha
         double numer = cnData.vU[b] - cnData.vU[a] ;
         for ( int j = cnData.nDegU ; j > mult ; -- j)
            vAlpha[j-mult-1] = numer / ( cnData.vU[a+j] - cnData.vU[a]) ;
         int r = cnData.nDegU - mult ;
         for ( int j = 1 ; j <= cnData.nDegU - mult ; ++j ) {
            int save = r - j ;
            int s = mult + j ;
            //for ( int row = 0 ; row < cnData.nCPV ; ++row) {
            //   for ( int k = cnData.nDegU ; k >= s ; --k ) {
            //      vBC[nCPU*row + k] = vAlpha[k-s]*vBC[nCPU*row + k] + ( 1 -  vAlfa[k-s]) * vBC[nCPU*row + k - 1]
            //   }
            //}
            if ( ! cnData.bRat ) {
               for ( int k = cnData.nDegU ; k >= s ; --k ) {
                  for ( int row = 0 ; row < cnData.nCPV ; ++row) {
                     mBC[k][row] = vAlpha[k-s] * mBC[k][row] + ( 1 -  vAlpha[k-s]) * mBC[k-1][row] ;
                  }
               }
            }
            else {
               for ( int k = cnData.nDegU ; k >= s ; --k ) {
                  for ( int row = 0 ; row < cnData.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 < cnData.nCPV ; ++row) {
                  mBC_next[save][row] =  mBC[cnData.nDegU][row] ;
               }
               if ( cnData.bRat )
                  for ( int row = 0 ; row < cnData.nCPV ; ++row) { 
                     mW_next[save][row] = mW[cnData.nDegU][row] ;
                  }
            }
         }
         mPC_strip.resize( cnData.nDegU * ( nb  + 1) + 1 , vCPV) ;
         mW_strip.resize( cnData.nDegU * ( nb  + 1) + 1, vV_W) ;
         if ( ! cnData.bRat)
            for ( int i = 0 ; i <= cnData.nDegU ; ++i) {
               for ( int row = 0 ; row < cnData.nCPV ; ++row ) {
                  mPC_strip[i+ nb * cnData.nDegU][row] = mBC[i][row] ;
               }
            }
         else {
            for ( int i = 0 ; i <= cnData.nDegU ; ++i) {
               for ( int row = 0 ; row < cnData.nCPV ; ++row ) {
                  mPC_strip[i+ nb * cnData.nDegU][row] = mBC[i][row]/mW[i][row] ;
                  mW_strip[i+ nb * cnData.nDegU][row] = mW[i][row] ;
               }
            }
         }
      }
      ++ nb ;
      // ho finito di definire la patch di Bezier attuale e passo alla successiva
         
      // aggiorno mBC con i valori della prossima pezza di Bezier // corrisponde a nb = nb + 1 
      if ( ! cnData.bRat){
         for (int i = 0 ; i < cnData.nDegU - 1 ; ++ i) {
            for ( int row = 0 ; row < cnData.nCPV ; ++row) {
               mBC[i][row] = mBC_next[i][row] ;
            }
         }
      }
      else {
         for (int i = 0 ; i < cnData.nDegU - 1 ; ++ i) {
            for ( int row = 0 ; row < cnData.nCPV ; ++row) {
               mBC[i][row] = mBC_next[i][row] ;
               mW[i][row] = mW_next[i][row] ;
            }
         }
      }

      if ( b < nU - 1 ) {
         for ( int i = cnData.nDegU - mult ; i <= cnData.nDegU ; ++ i) {
            for (int row = 0 ; row < cnData.nCPV ; ++ row ) {
               mBC[i][row] = cnData.mCP[b - cnData.nDegU + i + 1][row] ;
            }
         }
         if ( cnData.bRat ) {
            for ( int i = cnData.nDegU - mult ; i <= cnData.nDegU ; ++ i) {
               for (int row = 0 ; row < cnData.nCPV ; ++ row ) {
                  mW[i][row] = cnData.mW[b - cnData.nDegU + i + 1][row] ;
               }
            }
         }
         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 = cnData.nCPU ;
      mPC_strip.resize( cnData.nCPU, vCPV) ;
      mW_strip.resize( cnData.nCPU, vV_W) ;
      if ( ! cnData.bRat) {
         for ( int i = 0 ; i < nCPU_ref ; ++i) {
            for ( int row = 0 ; row < cnData.nCPV ; ++ row) {
               mPC_strip[i][row] = cnData.mCP[i][row] ;
            }
         }
      }
      else {
         for ( int i = 0 ; i < nCPU_ref ; ++i) {
            for ( int row = 0 ; row < cnData.nCPV ; ++ row) {
               mPC_strip[i][row] = cnData.mCP[i][row] ;
               mW_strip[i][row] = cnData.mW[i][row] ;
            }
         }
      }
      // devo vedere quante patch ci stanno prendendo i punti che ci sono
      //nb = (cnData.nCPU - 1) / cnData.nDegU ;
   }
   else
      nCPU_ref = cnData.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 = cnData.nDegV - 1 ;
   b = cnData.nDegV ;
   int nc = 0 ;  // numero di strisce in V ( lunghezza con V costante)
   vector<Point3d> vDegV(cnData.nDegV + 1) ;
   vector<Point3d> vDegV_1(cnData.nDegV - 1) ;
   vector< vector<Point3d>> m_BC1( nCPU_ref, vDegV) ;
   vector< vector<Point3d>> m_BC1_next( nCPU_ref, vDegV_1) ;
   DBLVECTOR vV1_W(cnData.nDegV + 1) ;
   DBLVECTOR vV2_W(cnData.nDegV - 1) ;
   vector<DBLVECTOR> mW1( nCPU_ref, vV1_W) ;
   vector<DBLVECTOR> mW1_next( nCPU_ref, vV2_W) ;
   DBLVECTOR vAlpha1( cnData.nDegV - 1) ;
   vector<vector<Point3d>> mPC_tot( nCPU_ref, vDegV) ;
   vector<DBLVECTOR> mW_tot( nCPU_ref, vV1_W) ;
   if ( ! cnData.bRat ) {
      for ( int i = 0 ; i < nCPU_ref ; ++i ) {
         for ( int row = 0 ; row <= cnData.nDegV ; ++row ) {
            m_BC1[i][row] = mPC_strip[i][row] ;
         }
      }
   }
   else {
      for ( int i = 0 ; i < nCPU_ref ; ++i ) {
         for (int row = 0 ; row <= cnData.nDegV ; ++ row) {
            mW1[i][row] = mW_strip[i][row] ;
            m_BC1[i][row] = mPC_strip[i][row] * mW_strip[i][row] ;
         }
      }
   }

   bRef = false ;
   // se la superficie è lineare nel parametro V allora è già in forma di Bezier in questo parametro
   if ( b == nV - 1 ) {
      nc = 1 ;
   }
   while ( b < nV - 1) {
      int i = b ;
      while ( b < nV - 1  &&  abs( cnData.vV[b+1] - cnData.vV[b]) < EPS_ZERO)
         ++ b ;
      int mult = b - i + 1 ;
      if ( mult < cnData.nDegV ) {
         bRef = true ;
         // calcolo numeratore e alpha
         double numer = cnData.vV[b] - cnData.vV[a] ;
         for ( int j = cnData.nDegV ; j > mult ; -- j)
            vAlpha1[j-mult-1] = numer / ( cnData.vV[a+j] - cnData.vV[a]) ;
         int r = cnData.nDegV - mult ;
         for ( int j = 1 ; j <= cnData.nDegV - mult ; ++j ) {
            int save = r - j ;
            int s = mult + j ;
            //for ( int row = 0 ; row < cnData.nCPV ; ++row) {
            //   for ( int k = cnData.nDegU ; k >= s ; --k ) {
            //      vBC[nCPU*row + k] = vAlpha1[k-s]*vBC[nCPU*row + k] + ( 1 -  vAlpha1[k-s]) * vBC[nCPU*row + k - 1]
            //   }
            //}
            if ( ! cnData.bRat) {
               for ( int k = 0 ; k < nCPU_ref ; ++k) {
                  for ( int row = cnData.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 = cnData.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 ( !cnData.bRat ){
                  for ( int i = 0 ; i < nCPU_ref ; ++i) {
                     m_BC1_next[i][save] =  m_BC1[i][cnData.nDegV] ;
                  }
               }
               else {
                  for ( int i = 0 ; i < nCPU_ref ; ++i) {
                     m_BC1_next[i][save] =  m_BC1[i][cnData.nDegV] ;
                     mW1_next[save] = mW1[cnData.nDegV] ;
                  }
               }
            }
         }
      }
      int nRef = cnData.nDegV * ( nc + 1) + 1 ;
      for ( int k = 0 ; k < nCPU_ref; ++k){
         mPC_tot[k].resize( nRef) ;
         mW_tot[k].resize( nRef) ;
      }
      if ( ! cnData.bRat)
         for ( int i = 0 ; i < nCPU_ref ; ++i) {
            for ( int row = 0 ; row <= cnData.nDegV ; ++row ) {
               mPC_tot[i][row + nc * cnData.nDegV] = m_BC1[i][row] ;
            }
         }
      else {
         for ( int i = 0 ; i < nCPU_ref ; ++i) {
            for ( int row = 0 ; row <= cnData.nDegV ; ++row ) {
               mPC_tot[i][row + nc * cnData.nDegV] = m_BC1[i][row]/mW1[i][row] ;
               mW_tot[i][row + nc * cnData.nDegV] = mW1[i][row] ;
            }
         }
      }
      ++ nc ;
      // ho finito di definire la patch di Bezier attuale e passo alla successiva

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

      if ( b < nV - 1) {
         for (int i = 0 ; i < nCPU_ref ; ++ i ) {
            for ( int row = cnData.nDegV - mult ; row <= cnData.nDegV ; ++ row) {
               m_BC1[i][row] = cnData.mCP[i][b - cnData.nDegV + row + 1] ;
            }
         }
         if ( cnData.bRat ) {
            for (int i = 0 ; i < nCPU_ref ; ++ i ) {
               for ( int row = cnData.nDegV - mult ; row <= cnData.nDegV ; ++ row) {
                  mW1[i][row] = cnData.mW[i][b - cnData.nDegV + row + 1] ;
               }
            }
         }
         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 = cnData.nCPV ;
      for ( int k = 0 ; k < nCPU_ref ; ++k){
         mPC_tot[k].resize( cnData.nCPV) ;
         mW_tot[k].resize( cnData.nCPV) ;
      }
      if ( ! cnData.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] ;
            }
         }
      }
      // devo vedere quante patch ci stanno prendendo i punti che ci sono
      //nc = (cnData.nCPV - 1) / cnData.nDegV ;
   }
   else
      nCPV_ref = cnData.nDegV * nc + 1 ;

   // finalmente setto la superficie di bezier totale divisa in nb patch in U e nc patch in V
   PtrOwner<ISurfBezier> pSrfBz( CreateSurfBezier()) ;
   if ( IsNull( pSrfBz))
      return nullptr ;
   pSrfBz->Init(cnData.nDegU, cnData.nDegV, nb, nc, cnData.bRat) ;
   if ( !cnData.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) ;
}