// Reconstruction - A tool for building cad models from meshes
// Copyright (C) 2010-2026 Eric Bechet, Borhen Louhichi
//
// See the LICENSE file for license information and contributions.
// Please report all bugs and problems to <bechet@cadxfem.org>.

#include "BSplineApproxSurface.h"

void myPlateSurfEval(Standard_Integer * Dimension,
                     // Dimension
                     Standard_Real    * UStartEnd,
                     // StartEnd[2] in U
                     Standard_Real    * VStartEnd,
                     // StartEnd[2] in V
                     Standard_Integer * FavorIso,
                     // Choice of constante, 1 for U, 2 for V
                     Standard_Real    * ConstParam,
                     // Value of constant parameter
                     Standard_Integer * NbParams,
                     // Number of parameters N
                     Standard_Real    * Parameters,
                     // Values of parameters,
                     Standard_Integer * UOrder,
                     // Derivative Request in U
                     Standard_Integer * VOrder,
                     // Derivative Request in V
                     Standard_Real    * Result,
                     // Result[Dimension,N]
                     Standard_Integer * ErrorCode)
// Error Code
{
  *ErrorCode = 0;
  Standard_Integer idim,jpar;
  Standard_Real Upar,Vpar;

// Dimension incorrecte
  if (*Dimension!=3)
  {
    *ErrorCode = 1;
  }

// Parametres incorrects
  if (*FavorIso==1)
  {
    Upar = *ConstParam;
    if ((Upar < UStartEnd[0]) || (Upar > UStartEnd[1]))
    {
      *ErrorCode = 2;
    }
    for (jpar=1;jpar<=*NbParams;jpar++)
    {
      Vpar = Parameters[jpar-1];
      if ((Vpar < VStartEnd[0]) || (Vpar > VStartEnd[1]))
      {
        *ErrorCode = 2;
      }
    }
  }
  else
  {
    Vpar = *ConstParam;
    if ((Vpar < VStartEnd[0]) || (Vpar > VStartEnd[1]))
    {
      *ErrorCode = 2;
    }
    for (jpar=1;jpar<=*NbParams;jpar++)
    {
      Upar = Parameters[jpar-1];
      if ((Upar < UStartEnd[0]) || (Upar > UStartEnd[1]))
      {
        *ErrorCode = 2;
      }
    }
  }

// Initialisation
  for (idim=1;idim<=*Dimension;idim++)
  {
    for (jpar=1;jpar<=*NbParams;jpar++)
    {
      Result[idim-1+ (jpar-1) * (*Dimension)] = 0.;
    }
  }


  Standard_Integer Order = *UOrder + *VOrder;
  gp_Pnt pnt;
// gp_Vec vect, v1, v2, v3, v4, v5, v6, v7, v8, v9;
  gp_Vec v1, v2, v3, v4, v5;

  if (*FavorIso==1)
  {
    Upar = *ConstParam;
    switch (Order)
    {
      case 0 :
        for (jpar=1;jpar<=*NbParams;jpar++)
        {
          Vpar = Parameters[jpar-1];
          pnt = fonct->Value(Upar,Vpar);
          Result[(jpar-1) * (*Dimension)] = pnt.X();
          Result[1+ (jpar-1) * (*Dimension)] = pnt.Y();
          Result[2+ (jpar-1) * (*Dimension)] = pnt.Z();
        }
        break;
      case 1 :
        for (jpar=1;jpar<=*NbParams;jpar++)
        {
          Vpar = Parameters[jpar-1];
          fonct->D1(Upar, Vpar, pnt, v1, v2);
          if (*UOrder==1)
          {
            Result[(jpar-1) * (*Dimension)] = v1.X();
            Result[1+ (jpar-1) * (*Dimension)] = v1.Y();
            Result[2+ (jpar-1) * (*Dimension)] = v1.Z();
          }
          else
          {
            Result[(jpar-1) * (*Dimension)] = v2.X();
            Result[1+ (jpar-1) * (*Dimension)] = v2.Y();
            Result[2+ (jpar-1) * (*Dimension)] = v2.Z();
          }
        }
        break;
      case 2 :
        for (jpar=1;jpar<=*NbParams;jpar++)
        {
          Vpar = Parameters[jpar-1];
          fonct->D2(Upar, Vpar, pnt, v1, v2, v3, v4, v5);
          if (*UOrder==2)
          {
            Result[(jpar-1) * (*Dimension)] = v3.X();
            Result[1+ (jpar-1) * (*Dimension)] = v3.Y();
            Result[2+ (jpar-1) * (*Dimension)] = v3.Z();
          }
          else if (*UOrder==1)
          {
            Result[(jpar-1) * (*Dimension)] = v5.X();
            Result[1+ (jpar-1) * (*Dimension)] = v5.Y();
            Result[2+ (jpar-1) * (*Dimension)] = v5.Z();
          }
          else if (*UOrder==0)
          {
            Result[(jpar-1) * (*Dimension)] = v4.X();
            Result[1+ (jpar-1) * (*Dimension)] = v4.Y();
            Result[2+ (jpar-1) * (*Dimension)] = v4.Z();
          }
        }
        break;
    }
  }
  else
  {
    Vpar = *ConstParam;
    switch (Order)
    {
      case 0 :
        for (jpar=1;jpar<=*NbParams;jpar++)
        {
          Upar = Parameters[jpar-1];
          pnt = fonct->Value(Upar,Vpar);
          Result[(jpar-1) * (*Dimension)] = pnt.X();
          Result[1+ (jpar-1) * (*Dimension)] = pnt.Y();
          Result[2+ (jpar-1) * (*Dimension)] = pnt.Z();
        }
        break;
      case 1 :
        for (jpar=1;jpar<=*NbParams;jpar++)
        {
          Upar = Parameters[jpar-1];
          fonct->D1(Upar, Vpar, pnt, v1, v2);
          if (*UOrder==1)
          {
            Result[(jpar-1) * (*Dimension)] = v1.X();
            Result[1+ (jpar-1) * (*Dimension)] = v1.Y();
            Result[2+ (jpar-1) * (*Dimension)] = v1.Z();
          }
          else
          {
            Result[(jpar-1) * (*Dimension)] = v2.X();
            Result[1+ (jpar-1) * (*Dimension)] = v2.Y();
            Result[2+ (jpar-1) * (*Dimension)] = v2.Z();
          }
        }
        break;
      case 2 :
        for (jpar=1;jpar<=*NbParams;jpar++)
        {
          Upar = Parameters[jpar-1];
          fonct->D2(Upar, Vpar, pnt, v1, v2, v3, v4, v5);
          if (*UOrder==2)
          {
            Result[(jpar-1) * (*Dimension)] = v3.X();
            Result[1+ (jpar-1) * (*Dimension)] = v3.Y();
            Result[2+ (jpar-1) * (*Dimension)] = v3.Z();
          }
          else if (*UOrder==1)
          {
            Result[(jpar-1) * (*Dimension)] = v5.X();
            Result[1+ (jpar-1) * (*Dimension)] = v5.Y();
            Result[2+ (jpar-1) * (*Dimension)] = v5.Z();
          }
          else if (*UOrder==0)
          {
            Result[(jpar-1) * (*Dimension)] = v4.X();
            Result[1+ (jpar-1) * (*Dimension)] = v4.Y();
            Result[2+ (jpar-1) * (*Dimension)] = v4.Z();
          }
        }
        break;
    }
  }

}



//----------------------------------------------------------------------------
/*
BSplineApproxSurface::BSplineApproxSurface(const Handle(GeomPlate_Surface)& SurfPlate,
        const AdvApp2Var_Criterion& PlateCrit,
                                           const Standard_Real Tol3d,
        const Standard_Integer Nbmax,
        const Standard_Integer dgmax,
        const GeomAbs_Shape Continuity,
        const Standard_Real EnlargeCoeff)
{
  myPlate = SurfPlate;
  fonct = myPlate;

  Standard_Real U0=0., U1=0., V0=0., V1=0.;
  myPlate->RealBounds(U0, U1, V0, V1);
  U0 = EnlargeCoeff * U0;
  U1 = EnlargeCoeff * U1;
  V0 = EnlargeCoeff * V0;
  V1 = EnlargeCoeff * V1;

  Standard_Integer nb1 = 0, nb2 = 0, nb3 = 1;
  Handle(TColStd_HArray1OfReal) nul1 =
     new TColStd_HArray1OfReal(1,1);
  nul1->Init(0.);
  Handle(TColStd_HArray2OfReal) nul2 =
     new TColStd_HArray2OfReal(1,1,1,4);
  nul2->Init(0.);
  Handle(TColStd_HArray1OfReal) eps3D =
     new TColStd_HArray1OfReal(1,1);
  eps3D->Init(Tol3d);
  Handle(TColStd_HArray2OfReal) epsfr =
     new TColStd_HArray2OfReal(1,1,1,4);
  epsfr->Init(Tol3d);
  GeomAbs_IsoType myType = GeomAbs_IsoV;
  Standard_Integer myPrec = 0;

  AdvApprox_DichoCutting myDec;

//POP pour WNT
  AdvApp2Var_EvaluatorFunc2Var ev = myPlateSurfEval;
  //AdvApp2Var_ApproxAFunc2Var AppPlate(nb1, nb2, nb3,
  AdvApp2Var_ApproxAFunc2Var AppPlate(nb1, nb2, nb3,
                 nul1,nul1,eps3D,
          nul2,nul2,epsfr,
          U0,U1,V0,V1,
          myType,
          Continuity, Continuity,
          myPrec,
          dgmax,dgmax,Nbmax,ev,
//          dgmax,dgmax,Nbmax,myPlateSurfEval,
          PlateCrit,myDec,myDec);
  mySurface = AppPlate.Surface(1);
  myAppError = AppPlate.MaxError(3,1);
  myCritError = AppPlate.CritError(3,1);
#if DEB
  cout<<"Approximation results"<<endl;
  cout<<"  Approximation error : "<<myAppError<<endl;
  cout<<"  Criterium error : "<<myCritError<<endl;
#endif
}
*/

//-----------------------------------------------------------------------------------

BSplineApproxSurface::BSplineApproxSurface(const Handle(GeomPlate_Surface) & SurfPlate,
    const Standard_Real Tol3d,
    const Standard_Integer Nbmax,
    const Standard_Integer dgmax,
    const Standard_Real dmax,
    const Standard_Integer CritOrder,
    const GeomAbs_Shape Continuity,
    const Standard_Real EnlargeCoeff)
{
  myPlate = SurfPlate;
  fonct = myPlate;

  TColgp_SequenceOfXY Seq2d;
  TColgp_SequenceOfXYZ Seq3d;

  if (CritOrder>=0)
  {

//    contraintes 2d d'ordre 0
    myPlate->Constraints(Seq2d);

//    contraintes 3d correspondantes sur plate
    Standard_Integer i,nbp=Seq2d.Length();
    for (i=1;i<=nbp;i++)
    {
      gp_XY P2d=Seq2d.Value(i);
      gp_Pnt PP;
      gp_Vec v1h,v2h,v3h;
      if (CritOrder==0)
      {
//    a l'ordre 0
        fonct->D0(P2d.X(),P2d.Y(),PP);
        gp_XYZ P3d(PP.X(),PP.Y(),PP.Z());
        Seq3d.Append(P3d);

      }
      else
      {
//    a l'ordre 1
        fonct->D1(P2d.X(),P2d.Y(),PP,v1h,v2h);
        v3h=v1h^v2h;
        gp_XYZ P3d(v3h.X(),v3h.Y(),v3h.Z());
        Seq3d.Append(P3d);
      }
    }
  }


  ///*********************** le 30-09-2010 ******************************//

  std::cout<<"Seq2d.Length() : "<< Seq2d.Length() <<std::endl;
  for (int i=1; i<= Seq2d.Length(); ++i)
  {
    std::cout<<"Seq2d.Value(i) :  X : "<< Seq2d.Value(i).X() <<" Y : "<< Seq2d.Value(i).Y() <<std::endl;

  }

  std::cout<<"Seq3d.Length() : "<< Seq3d.Length() <<std::endl;
  for (int i=1; i<= Seq3d.Length(); ++i)
  {
    std::cout<<"Seq3d.Value( "<< i <<" ) :  X : "<< Seq3d.Value(i).X() <<" Y : "<< Seq3d.Value(i).Y() <<" Z : "<< Seq3d.Value(i).Z() <<std::endl;

  }

  ///********************************************************************//



  Standard_Real U0=0., U1=0., V0=0., V1=0.;
  myPlate->RealBounds(U0, U1, V0, V1);

  ///********** B le 30/09/2010
  std::cout<<" U0 : "<< U0 <<" U1 : "<< U1 <<" V0 : "<< V0 <<" V1 : "<< V1 <<std::endl;
  ///**************************

  U0 = EnlargeCoeff * U0;
  U1 = EnlargeCoeff * U1;
  V0 = EnlargeCoeff * V0;
  V1 = EnlargeCoeff * V1;

  ///********** B le 30/09/2010
  std::cout<<" U0 : "<< U0 <<" U1 : "<< U1 <<" V0 : "<< V0 <<" V1 : "<< V1 <<std::endl;
  ///**************************


  Standard_Real seuil = Tol3d;
  if (CritOrder==0&&Tol3d<10*dmax)
  {
    seuil=10*dmax;
#if DEB
    cout<<"Seuil G0 choisi trop faible par rapport au contour. On prend "<<seuil<<endl;
#endif


  }

  if (CritOrder==1&&Tol3d<10*dmax)
  {
    seuil=10*dmax;
#if DEB
    cout<<"Seuil G1 choisi trop faible par rapport au contour. On prend "<<seuil<<endl;
#endif
  }
  Standard_Integer nb1 = 0, nb2 = 0, nb3 = 1;
  Handle(TColStd_HArray1OfReal) nul1 =
    new TColStd_HArray1OfReal(1,1);
  nul1->Init(0.);
  Handle(TColStd_HArray2OfReal) nul2 =
    new TColStd_HArray2OfReal(1,1,1,4);
  nul2->Init(0.);
  Handle(TColStd_HArray1OfReal) eps3D =
    new TColStd_HArray1OfReal(1,1);
  eps3D->Init(Tol3d);
  Handle(TColStd_HArray2OfReal) epsfr =
    new TColStd_HArray2OfReal(1,1,1,4);
  epsfr->Init(Tol3d);

  GeomAbs_IsoType myType = GeomAbs_IsoV;

  Standard_Integer myPrec = 0;

  AdvApprox_DichoCutting myDec;

  if (CritOrder==-1)
  {
    myPrec = 1;
// POP pour NT
    AdvApp2Var_EvaluatorFunc2Var ev = myPlateSurfEval;
    ApproxAFunc2Var AppPlate(nb1, nb2, nb3,
                             nul1,nul1,eps3D,
                             nul2,nul2,epsfr,
                             U0,U1,V0,V1,
                             myType,
                             Continuity, Continuity,
                             myPrec,
                             dgmax,dgmax,Nbmax,ev,
                             myDec,myDec);
    mySurface = AppPlate.Surface(1);
    myAppError = AppPlate.MaxError(3,1);
    myCritError = 0.;
#if DEB
    cout<<"Approximation results"<<endl;
    cout<<"  Approximation error : "<<myAppError<<endl;
#endif
  }


  else if (CritOrder==0)
  {

    GeomPlate_PlateG0Criterion Crit0(Seq2d,Seq3d,seuil);
// POP pour NT
    AdvApp2Var_EvaluatorFunc2Var ev = myPlateSurfEval;


    ApproxAFunc2Var AppPlate(nb1, nb2, nb3,
                             nul1,nul1,eps3D,
                             nul2,nul2,epsfr,
                             U0,U1,V0,V1,
                             myType,
                             Continuity, Continuity,
                             myPrec,
                             dgmax,dgmax,Nbmax,ev,
//     dgmax,dgmax,Nbmax,myPlateSurfEval,
                             Crit0,myDec,myDec);

    mySurface = AppPlate.Surface(1);
    myAppError = AppPlate.MaxError(3,1);
    myCritError = AppPlate.CritError(3,1);

#if DEB
    cout<<"Approximation results"<<endl;
    cout<<"  Approximation error : "<<myAppError<<endl;
    cout<<"  Criterium error : "<<myCritError<<endl;
#endif
  }
  else if (CritOrder==1)
  {
    GeomPlate_PlateG1Criterion Crit1(Seq2d,Seq3d,seuil);
// POP pour NT
    AdvApp2Var_EvaluatorFunc2Var ev = myPlateSurfEval;
    ApproxAFunc2Var AppPlate(nb1, nb2, nb3,
                             nul1,nul1,eps3D,
                             nul2,nul2,epsfr,
                             U0,U1,V0,V1,
                             myType,
                             Continuity, Continuity,
                             myPrec,
                             dgmax,dgmax,Nbmax,ev,
//     dgmax,dgmax,Nbmax,myPlateSurfEval,
                             Crit1,myDec,myDec);
    mySurface = AppPlate.Surface(1);
    myAppError = AppPlate.MaxError(3,1);
    myCritError = AppPlate.CritError(3,1);
#if DEB
    cout<<"Approximation results"<<endl;
    cout<<"  Approximation error : "<<myAppError<<endl;
    cout<<"  Criterium error : "<<myCritError<<endl;
#endif
  }
}


//===================================================================================


//=======================================================================
//function : Surface
//purpose  :
//=======================================================================

Handle(Geom_BSplineSurface) BSplineApproxSurface::Surface() const
{
  return mySurface;
}


//=======================================================================
//function : ApproxError
//purpose  :
//=======================================================================

Standard_Real BSplineApproxSurface::ApproxError() const
{
  return myAppError;
}


//=======================================================================
//function : CriterionError
//purpose  :
//=======================================================================

Standard_Real BSplineApproxSurface::CriterionError() const
{
  return myCritError;
}