// Copyright (c) 1995-1999 Matra Datavision
// Copyright (c) 1999-2014 OPEN CASCADE SAS
//
// This file is part of Open CASCADE Technology software library.
//
// This library is free software; you can redistribute it and/or modify it under
// the terms of the GNU Lesser General Public License version 2.1 as published
// by the Free Software Foundation, with special exception defined in the file
// OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
// distribution for complete text of the license and disclaimer of any warranty.
//
// Alternatively, this file may be used under the terms of Open CASCADE
// commercial license or contractual agreement.

#include <LProp_Status.hxx>
#include <LProp_NotDefined.hxx>
#include <Standard_OutOfRange.hxx>
#include <Standard_DomainError.hxx>
#include <CSLib.hxx>
#include <CSLib_DerivativeStatus.hxx>
#include <CSLib_NormalStatus.hxx>
#include <TColgp_Array2OfVec.hxx>
#include <math_DirectPolynomialRoots.hxx>

static const Standard_Real MinStep   = 1.0e-7;

static Standard_Boolean IsTangentDefined (LProp_SLProps& SProp,
                                          const Standard_Integer cn,
                                          const Standard_Real linTol,
                                          const Standard_Integer  Derivative, 
                                          Standard_Integer&       Order,
                                          LProp_Status&  theStatus)
{
  Standard_Real Tol = linTol * linTol;
  gp_Vec V[2];
  Order = 0;

  while (Order < 3)
  {
    Order++;
    if(cn >= Order)
    {
      switch(Order)
      {
      case 1:
        V[0] = SProp.D1U();
        V[1] = SProp.D1V();
        break;
      case 2:
        V[0] = SProp.D2U();
        V[1] = SProp.D2V();
        break;
      }//switch(Order)

      if(V[Derivative].SquareMagnitude() > Tol)
      {
        theStatus = LProp_Defined;
        return Standard_True;
      }
    }//if(cn >= Order)
    else
    {
      theStatus = LProp_Undefined;
      return Standard_False;
    }
  }

  return Standard_False;
}

LProp_SLProps::LProp_SLProps (const Surface& S, 
                              const Standard_Real U,  
                              const Standard_Real V, 
                              const Standard_Integer N, 
                              const Standard_Real Resolution) 
        : mySurf(S),myDerOrder(N), myCN(4), // (Tool::Continuity(S)),
            myLinTol(Resolution)
{
  Standard_OutOfRange_Raise_if(N < 0 || N > 2, 
                        "LProp_SLProps::LProp_SLProps()");

  SetParameters(U, V);
}

LProp_SLProps::LProp_SLProps (const Surface& S, 
                              const Standard_Integer  N, 
                              const Standard_Real     Resolution) 
        : mySurf(S), myU(RealLast()), myV(RealLast()), myDerOrder(N),
            myCN(4), // (Tool::Continuity(S))
            myLinTol(Resolution), 
            myUTangentStatus (LProp_Undecided),
            myVTangentStatus (LProp_Undecided),
            myNormalStatus   (LProp_Undecided),
            myCurvatureStatus(LProp_Undecided)
{
  Standard_OutOfRange_Raise_if(N < 0 || N > 2, 
                      "LProp_SLProps::LProp_SLProps()");
}

LProp_SLProps::LProp_SLProps (const Standard_Integer  N, 
                              const Standard_Real Resolution) 
        : myU(RealLast()), myV(RealLast()), myDerOrder(N), myCN(0),
        myLinTol(Resolution),
        myUTangentStatus (LProp_Undecided),
        myVTangentStatus (LProp_Undecided),
        myNormalStatus   (LProp_Undecided),
        myCurvatureStatus(LProp_Undecided)
{
  Standard_OutOfRange_Raise_if(N < 0 || N > 2, 
                "LProp_SLProps::LProp_SLProps() bad level");
}

void LProp_SLProps::SetSurface (const Surface& S ) 
{
  mySurf = S;
  myCN = 4; // =Tool::Continuity(S);
}

void LProp_SLProps::SetParameters (const Standard_Real U, 
                                   const Standard_Real V)
{
  myU = U;
  myV = V;
  switch (myDerOrder)
  {
  case 0:
    Tool::Value(mySurf, myU, myV, myPnt);
    break;
  case 1:
    Tool::D1(mySurf, myU, myV, myPnt, myD1u, myD1v);
    break;
  case 2:
    Tool::D2(mySurf, myU, myV, myPnt, myD1u, myD1v, myD2u, myD2v, myDuv);
    break;
  }

  myUTangentStatus  = LProp_Undecided;
  myVTangentStatus  = LProp_Undecided;
  myNormalStatus    = LProp_Undecided;
  myCurvatureStatus = LProp_Undecided;
}

const gp_Pnt& LProp_SLProps::Value() const
{
  return myPnt;
}

const gp_Vec& LProp_SLProps::D1U()
{
  if (myDerOrder < 1)
  {
    myDerOrder =1;
    Tool::D1(mySurf,myU,myV,myPnt,myD1u,myD1v);
  }

  return myD1u;
}

const gp_Vec& LProp_SLProps::D1V()
{
  if (myDerOrder < 1)
  {
    myDerOrder =1;
    Tool::D1(mySurf,myU,myV,myPnt,myD1u,myD1v);
  }

  return myD1v;
}

const gp_Vec& LProp_SLProps::D2U()
{
  if (myDerOrder < 2) 
  {
    myDerOrder =2;
    Tool::D2(mySurf,myU,myV,myPnt,myD1u,myD1v,myD2u,myD2v,myDuv);
  }

  return myD2u;
}

const gp_Vec& LProp_SLProps::D2V()
{
  if (myDerOrder < 2) 
  {
    myDerOrder =2;
    Tool::D2(mySurf,myU,myV,myPnt,myD1u,myD1v,myD2u,myD2v,myDuv);
  }

  return myD2v;
}

const gp_Vec& LProp_SLProps::DUV()
{
  if (myDerOrder < 2) 
  {
    myDerOrder =2;
    Tool::D2(mySurf,myU,myV,myPnt,myD1u,myD1v,myD2u,myD2v,myDuv);
  }

  return myDuv;
}

Standard_Boolean LProp_SLProps::IsTangentUDefined ()
{
  if (myUTangentStatus == LProp_Undefined)
    return Standard_False;
  else if (myUTangentStatus >= LProp_Defined)
    return Standard_True; 

  // uTangentStatus == Lprop_Undecided 
  // we have to calculate the first non null U derivative
  return IsTangentDefined(*this, myCN, myLinTol, 0, 
          mySignificantFirstDerivativeOrderU, myUTangentStatus);
}

void LProp_SLProps::TangentU (gp_Dir& D) 
{
  if(!IsTangentUDefined())
    throw LProp_NotDefined();

  if(mySignificantFirstDerivativeOrderU == 1)
    D = gp_Dir(myD1u);
  else
  {
    const Standard_Real DivisionFactor = 1.e-3;
    Standard_Real anUsupremum, anUinfium;
    Standard_Real anVsupremum, anVinfium;
    Tool::Bounds(mySurf,anUinfium,anVinfium,anUsupremum,anVsupremum);
    Standard_Real du;
    if((anUsupremum >= RealLast()) || (anUinfium <= RealFirst()))
      du = 0.0;
    else
      du = anUsupremum-anUinfium;
    
    const Standard_Real aDeltaU = Max(du*DivisionFactor,MinStep);

    gp_Vec V = myD2u;
    
    Standard_Real u;
          
    if(myU-anUinfium < aDeltaU)
      u = myU+aDeltaU;
    else
      u = myU-aDeltaU;
    
    gp_Pnt P1, P2;
    Tool::Value(mySurf, Min(myU, u),myV,P1);
    Tool::Value(mySurf, Max(myU, u),myV,P2);
    
    gp_Vec V1(P1,P2);
    Standard_Real aDirFactor = V.Dot(V1);
    
    if(aDirFactor < 0.0)
      V = -V;
    
    D = gp_Dir(V);
  }
}

Standard_Boolean LProp_SLProps::IsTangentVDefined ()
{
  if (myVTangentStatus == LProp_Undefined)
    return Standard_False;
  else if (myVTangentStatus >= LProp_Defined)
    return Standard_True; 

  // vTangentStatus == Lprop_Undecided 
  // we have to calculate the first non null V derivative
  return IsTangentDefined(*this, myCN, myLinTol, 1, 
          mySignificantFirstDerivativeOrderV, myVTangentStatus);
}

void LProp_SLProps::TangentV (gp_Dir& D)
{
  if(!IsTangentVDefined())
    throw LProp_NotDefined();
  
  if(mySignificantFirstDerivativeOrderV == 1)
    D = gp_Dir(myD1v);
  else
  {
    const Standard_Real DivisionFactor = 1.e-3;
    Standard_Real anUsupremum, anUinfium;
    Standard_Real anVsupremum, anVinfium;
    Tool::Bounds(mySurf,anUinfium,anVinfium,anUsupremum,anVsupremum);
    
    Standard_Real dv;
    if((anVsupremum >= RealLast()) || (anVinfium <= RealFirst()))
      dv = 0.0;
    else
      dv = anVsupremum-anVinfium;
    
    const Standard_Real aDeltaV = Max(dv*DivisionFactor,MinStep);

    gp_Vec V = myD2v;
    
    Standard_Real v;
          
    if(myV-anVinfium < aDeltaV)
      v = myV+aDeltaV;
    else
      v = myV-aDeltaV;
    
    gp_Pnt P1, P2;
    Tool::Value(mySurf, myU, Min(myV, v),P1);
    Tool::Value(mySurf, myU, Max(myV, v),P2);
    
    gp_Vec V1(P1,P2);
    Standard_Real aDirFactor = V.Dot(V1);
    
    if(aDirFactor < 0.0)
      V = -V;
    
    D = gp_Dir(V);
  }
}

Standard_Boolean LProp_SLProps::IsNormalDefined()
{
  if (myNormalStatus == LProp_Undefined)
    return Standard_False;
  else if (myNormalStatus >= LProp_Defined)
    return Standard_True;

  // status = UnDecided 
  // first try the standard computation of the normal.
  CSLib_DerivativeStatus aStatus = CSLib_Done;
  CSLib::Normal(myD1u, myD1v, myLinTol, aStatus, myNormal);
  if (aStatus == CSLib_Done)
  {
    myNormalStatus = LProp_Computed;
    return Standard_True;
  }

  // else solve the degenerated case only if continuity >= 2

  myNormalStatus = LProp_Undefined;
  return Standard_False;
}

const gp_Dir& LProp_SLProps::Normal ()
{
  if(!IsNormalDefined())
  {
    throw LProp_NotDefined();
  }

  return myNormal;
}


Standard_Boolean LProp_SLProps::IsCurvatureDefined ()
{
  if (myCurvatureStatus == LProp_Undefined)
    return Standard_False;
  else if (myCurvatureStatus >= LProp_Defined)
    return Standard_True;
  
  if(myCN < 2)
  {
    myCurvatureStatus = LProp_Undefined;
    return Standard_False;
  }

  // status = UnDecided 
  if (!IsNormalDefined())
  {
    myCurvatureStatus = LProp_Undefined;
    return Standard_False;
  }

  // pour eviter un plantage dans le cas du caro pointu
  // en fait on doit pouvoir calculer les courbure
  // avoir
  if(!IsTangentUDefined() || !IsTangentVDefined())
  {
    myCurvatureStatus = LProp_Undefined;
    return Standard_False;
  }

  // here we compute the curvature features of the surface

  gp_Vec Norm (myNormal);

  Standard_Real E = myD1u.SquareMagnitude();
  Standard_Real F = myD1u.Dot(myD1v);
  Standard_Real G = myD1v.SquareMagnitude();

  if(myDerOrder < 2)
    this->D2U();

  Standard_Real L = Norm.Dot(myD2u);
  Standard_Real M = Norm.Dot(myDuv);
  Standard_Real N = Norm.Dot(myD2v);

  Standard_Real A = E * M - F * L;
  Standard_Real B = E * N - G * L;
  Standard_Real C = F * N - G * M;

  Standard_Real MaxABC = Max(Max(Abs(A),Abs(B)),Abs(C));
  if (MaxABC < RealEpsilon())    // ombilic
  {
    myMinCurv = N / G;
    myMaxCurv = myMinCurv;
    myDirMinCurv = gp_Dir (myD1u);
    myDirMaxCurv = gp_Dir (myD1u.Crossed(Norm));
    myMeanCurv = myMinCurv;            // (Cmin + Cmax) / 2.
    myGausCurv = myMinCurv * myMinCurv;  // (Cmin * Cmax)
    myCurvatureStatus = LProp_Computed;
    return Standard_True;
  }

  A = A / MaxABC;
  B = B / MaxABC;
  C = C / MaxABC;
  Standard_Real Curv1, Curv2, Root1, Root2;
  gp_Vec VectCurv1, VectCurv2;

  if (Abs(A) > RealEpsilon())
  {
    math_DirectPolynomialRoots Root (A, B, C);
    if(Root.NbSolutions() != 2)
    {
      myCurvatureStatus = LProp_Undefined;
      return Standard_False;
    }
    else
    {
      Root1 = Root.Value(1);
      Root2 = Root.Value(2);
      Curv1 = ((L * Root1 + 2. * M) * Root1 + N) /
	      ((E * Root1 + 2. * F) * Root1 + G);
      Curv2 = ((L * Root2 + 2. * M) * Root2 + N) /
	      ((E * Root2 + 2. * F) * Root2 + G);
      VectCurv1 = Root1 * myD1u + myD1v;
      VectCurv2 = Root2 * myD1u + myD1v;
    }
  }
  else if (Abs(C) > RealEpsilon())
  {
    math_DirectPolynomialRoots Root(C, B, A);

    if((Root.NbSolutions() != 2))
    {
      myCurvatureStatus = LProp_Undefined;
      return Standard_False;
    }
    else
    {
      Root1 = Root.Value(1);
      Root2 = Root.Value(2);
      Curv1 = ((N * Root1 + 2. * M) * Root1 + L) /
	      ((G * Root1 + 2. * F) * Root1 + E);
      Curv2 = ((N * Root2 + 2. * M) * Root2 + L) /
	      ((G * Root2 + 2. * F) * Root2 + E);
      VectCurv1 = myD1u + Root1 * myD1v;
      VectCurv2 = myD1u + Root2 * myD1v;
    }
  }
  else
  {
    Curv1 = L / E;
    Curv2 = N / G;
    VectCurv1 = myD1u;
    VectCurv2 = myD1v;
  }

  if (Curv1 < Curv2)
  {
    myMinCurv = Curv1;
    myMaxCurv = Curv2;
    myDirMinCurv = gp_Dir (VectCurv1);
    myDirMaxCurv = gp_Dir (VectCurv2);
  }
  else
  {
    myMinCurv = Curv2;
    myMaxCurv = Curv1;
    myDirMinCurv = gp_Dir (VectCurv2);
    myDirMaxCurv = gp_Dir (VectCurv1);
  }

  myMeanCurv = ((N * E) - (2. * M * F) + (L * G)) // voir Farin p.282
               / (2. * ((E * G) - (F * F)));
  myGausCurv = ((L * N) - (M * M)) 
               / ((E * G) - (F * F));
  myCurvatureStatus = LProp_Computed;
  return Standard_True;
}

Standard_Boolean LProp_SLProps::IsUmbilic ()
{
  if(!IsCurvatureDefined())
    throw LProp_NotDefined();
  
  return Abs(myMaxCurv - myMinCurv) < Abs(Epsilon(myMaxCurv));
}

Standard_Real LProp_SLProps::MaxCurvature ()
{
  if(!IsCurvatureDefined())
    throw LProp_NotDefined();

  return myMaxCurv;
}

Standard_Real LProp_SLProps::MinCurvature ()
{
  if(!IsCurvatureDefined())
    throw LProp_NotDefined();
  
  return myMinCurv;
}

void LProp_SLProps::CurvatureDirections(gp_Dir& Max, gp_Dir& Min)
{
  if(!IsCurvatureDefined())
    throw LProp_NotDefined();

  Max = myDirMaxCurv;
  Min = myDirMinCurv;
}

Standard_Real LProp_SLProps::MeanCurvature ()
{
  if(!IsCurvatureDefined())
    throw LProp_NotDefined();

  return myMeanCurv;
}

Standard_Real LProp_SLProps::GaussianCurvature ()
{
  if(!IsCurvatureDefined())
    throw LProp_NotDefined();

  return myGausCurv;
}