// Copyright (c) 2025 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.

#ifndef _LProp_CurveUtils_HeaderFile
#define _LProp_CurveUtils_HeaderFile

#include "LProp_NotDefined.hxx"
#include "LProp_Status.hxx"
#include "Standard_Handle.hxx"

#include <cmath>

//! Template utility functions for CLProps curve local property computation.
//! Provides access policies (DirectAccess, ToolAccess) and shared algorithms
//! for evaluating derivatives, tangents, curvature, normals, and centres of curvature.
namespace LProp_CurveUtils
{

// ==================== Deref helpers ====================

//! Dereference by-value or reference types (identity).
template <typename T>
T& Deref(T& theObj)
{
  return theObj;
}

//! Dereference occ::handle types.
template <typename T>
T& Deref(occ::handle<T>& theHandle)
{
  return *theHandle;
}

//! Dereference const occ::handle types.
template <typename T>
const T& Deref(const occ::handle<T>& theHandle)
{
  return *theHandle;
}

// ==================== Access Policies ====================

//! Direct access policy: calls D0/D1/D2/D3 methods on the curve object.
//! Works with occ::handle<T> and by-value curve types.
struct DirectAccess
{
  template <typename C, typename P>
  static void D0(C& theCurve, double theU, P& thePnt)
  {
    Deref(theCurve).D0(theU, thePnt);
  }

  template <typename C, typename P, typename V>
  static void D1(C& theCurve, double theU, P& thePnt, V& theV1)
  {
    Deref(theCurve).D1(theU, thePnt, theV1);
  }

  template <typename C, typename P, typename V>
  static void D2(C& theCurve, double theU, P& thePnt, V& theV1, V& theV2)
  {
    Deref(theCurve).D2(theU, thePnt, theV1, theV2);
  }

  template <typename C, typename P, typename V>
  static void D3(C& theCurve, double theU, P& thePnt, V& theV1, V& theV2, V& theV3)
  {
    Deref(theCurve).D3(theU, thePnt, theV1, theV2, theV3);
  }

  template <typename C>
  static double FirstParameter(C& theCurve)
  {
    return Deref(theCurve).FirstParameter();
  }

  template <typename C>
  static double LastParameter(C& theCurve)
  {
    return Deref(theCurve).LastParameter();
  }
};

//! Tool-based access policy: delegates to static Tool methods.
//! Used for HLRBRep types where Tool class provides the interface.
template <typename Tool>
struct ToolAccess
{
  template <typename C, typename P>
  static void D0(C& theCurve, double theU, P& thePnt)
  {
    Tool::Value(theCurve, theU, thePnt);
  }

  template <typename C, typename P, typename V>
  static void D1(C& theCurve, double theU, P& thePnt, V& theV1)
  {
    Tool::D1(theCurve, theU, thePnt, theV1);
  }

  template <typename C, typename P, typename V>
  static void D2(C& theCurve, double theU, P& thePnt, V& theV1, V& theV2)
  {
    Tool::D2(theCurve, theU, thePnt, theV1, theV2);
  }

  template <typename C, typename P, typename V>
  static void D3(C& theCurve, double theU, P& thePnt, V& theV1, V& theV2, V& theV3)
  {
    Tool::D3(theCurve, theU, thePnt, theV1, theV2, theV3);
  }

  template <typename C>
  static double FirstParameter(C& theCurve)
  {
    return Tool::FirstParameter(theCurve);
  }

  template <typename C>
  static double LastParameter(C& theCurve)
  {
    return Tool::LastParameter(theCurve);
  }
};

// ==================== Algorithm Utilities ====================

//! Evaluate curve derivatives at parameter theU up to specified order.
//! @param[in,out] theCurve    curve object
//! @param[in]     theU        parameter value
//! @param[in]     theOrder    maximum derivative order (0-3)
//! @param[out]    thePnt      evaluated point
//! @param[out]    theDerivArr derivative array (size >= theOrder)
template <typename Access, typename Curve, typename Pnt, typename Vec>
void EvalDerivatives(Curve& theCurve, double theU, int theOrder, Pnt& thePnt, Vec* theDerivArr)
{
  switch (theOrder)
  {
    case 0:
      Access::D0(theCurve, theU, thePnt);
      break;
    case 1:
      Access::D1(theCurve, theU, thePnt, theDerivArr[0]);
      break;
    case 2:
      Access::D2(theCurve, theU, thePnt, theDerivArr[0], theDerivArr[1]);
      break;
    case 3:
      Access::D3(theCurve, theU, thePnt, theDerivArr[0], theDerivArr[1], theDerivArr[2]);
      break;
  }
}

//! Compute tangent direction with sign correction for higher-order derivatives.
//! When the first significant derivative has order > 1, the sign of the tangent
//! is determined by comparing with the chord direction near the point.
//! @param[in,out] theCurve    curve object
//! @param[in]     theU        current parameter
//! @param[in]     theDerivArr derivative array
//! @param[in]     theRefPnt   reference point (used only for type deduction)
//! @param[in]     theSigOrder order of first significant derivative
//! @param[out]    theDir      computed tangent direction
template <typename Access, typename Curve, typename Vec, typename Pnt, typename Dir>
void ComputeTangent(Curve&     theCurve,
                    double     theU,
                    const Vec* theDerivArr,
                    const Pnt& /*theRefPnt*/,
                    int  theSigOrder,
                    Dir& theDir)
{
  if (theSigOrder == 1)
  {
    theDir = Dir(theDerivArr[0]);
    return;
  }

  constexpr double THE_DIVISION_FACTOR = 1.0e-3;
  constexpr double THE_MIN_STEP        = 1.0e-7;

  const double anUsupremum = Access::LastParameter(theCurve);
  const double anUinfimum  = Access::FirstParameter(theCurve);

  double aDu;
  if ((anUsupremum >= RealLast()) || (anUinfimum <= RealFirst()))
    aDu = 0.0;
  else
    aDu = anUsupremum - anUinfimum;

  const double aDelta = std::max(aDu * THE_DIVISION_FACTOR, THE_MIN_STEP);

  Vec aV = theDerivArr[theSigOrder - 1];

  double anOtherU;
  if (theU - anUinfimum < aDelta)
    anOtherU = theU + aDelta;
  else
    anOtherU = theU - aDelta;

  Pnt aP1, aP2;
  Access::D0(theCurve, std::min(theU, anOtherU), aP1);
  Access::D0(theCurve, std::max(theU, anOtherU), aP2);

  Vec aChord(aP1, aP2);
  if (aV.Dot(aChord) < 0.0)
    aV = -aV;

  theDir = Dir(aV);
}

//! Compute curvature from first and second derivatives.
//! Returns |D1 x D2| / |D1|^3, or 0 if derivatives are collinear or D2 is null.
//! @param[in]  theD1    first derivative
//! @param[in]  theD2    second derivative
//! @param[in]  theTolSq squared linear tolerance
//! @return curvature value
template <typename Vec>
double ComputeCurvature(const Vec& theD1, const Vec& theD2, double theTolSq)
{
  const double aDD1 = theD1.SquareMagnitude();
  const double aDD2 = theD2.SquareMagnitude();

  if (aDD2 <= theTolSq)
    return 0.0;

  const double aN = theD1.CrossSquareMagnitude(theD2);
  const double aT = aN / aDD1 / aDD2;
  if (aT <= theTolSq)
    return 0.0;

  return sqrt(aN) / aDD1 / sqrt(aDD1);
}

//! Compute normal direction from first and second derivatives.
//! Normal = D2*(D1*D1) - D1*(D1*D2), using the vector triple product identity.
//! @param[in]  theD1   first derivative
//! @param[in]  theD2   second derivative
//! @param[out] theDir  computed normal direction
template <typename Vec, typename Dir>
void ComputeNormal(const Vec& theD1, const Vec& theD2, Dir& theDir)
{
  Vec aNorm = theD2 * (theD1 * theD1) - theD1 * (theD1 * theD2);
  theDir    = Dir(aNorm);
}

//! Compute centre of curvature from point, derivatives, and curvature.
//! Centre = Point + Normal / Curvature.
//! @param[in]  thePnt       current point
//! @param[in]  theD1        first derivative
//! @param[in]  theD2        second derivative
//! @param[in]  theCurvature curvature value (must be non-zero)
//! @param[out] theCentre    computed centre of curvature
template <typename Vec, typename Pnt>
void ComputeCentreOfCurvature(const Pnt& thePnt,
                              const Vec& theD1,
                              const Vec& theD2,
                              double     theCurvature,
                              Pnt&       theCentre)
{
  Vec aNorm = theD2 * (theD1 * theD1) - theD1 * (theD1 * theD2);
  aNorm.Normalize();
  aNorm.Divide(theCurvature);
  theCentre = thePnt.Translated(aNorm);
}

// ==================== Higher-Level Method Wrappers ====================

//! SetParameter wrapper: sets parameter, evaluates derivatives, resets tangent status.
//! @param[in,out] theCurve    curve object
//! @param[in]     theU        parameter value
//! @param[out]    theStoredU  stored parameter field
//! @param[in]     theDerOrder current derivative order
//! @param[out]    thePnt      evaluated point
//! @param[out]    theDerivArr derivative array
//! @param[out]    theTanStatus tangent status field
template <typename Access, typename Curve, typename Pnt, typename Vec>
void SetParameter(Curve&        theCurve,
                  double        theU,
                  double&       theStoredU,
                  int           theDerOrder,
                  Pnt&          thePnt,
                  Vec*          theDerivArr,
                  LProp_Status& theTanStatus)
{
  theStoredU = theU;
  EvalDerivatives<Access>(theCurve, theU, theDerOrder, thePnt, theDerivArr);
  theTanStatus = LProp_Undecided;
}

//! Ensure derivatives up to the required order are computed.
//! @param[in,out] theCurve    curve object
//! @param[in]     theU        parameter value
//! @param[in,out] theDerOrder current derivative order (upgraded if needed)
//! @param[in]     theRequired required derivative order
//! @param[out]    thePnt      evaluated point
//! @param[out]    theDerivArr derivative array
//! @return reference to theDerivArr[theRequired-1]
template <typename Access, typename Curve, typename Pnt, typename Vec>
const Vec& EnsureDeriv(Curve& theCurve,
                       double theU,
                       int&   theDerOrder,
                       int    theRequired,
                       Pnt&   thePnt,
                       Vec*   theDerivArr)
{
  if (theDerOrder < theRequired)
  {
    theDerOrder = theRequired;
    EvalDerivatives<Access>(theCurve, theU, theDerOrder, thePnt, theDerivArr);
  }
  return theDerivArr[theRequired - 1];
}

//! IsTangentDefined: searches for first non-null derivative.
//! Calls theProps.D1(), theProps.D2(), theProps.D3() to upgrade derivatives as needed.
//! @param[in,out] theProps     CLProps object
//! @param[in]     theCN        continuity order
//! @param[in]     theLinTol    linear tolerance
//! @param[out]    theSigOrder  order of first significant derivative
//! @param[in,out] theTanStatus tangent status field
//! @return true if tangent is defined
template <typename Vec, typename Props>
bool IsTangentDefined(Props&        theProps,
                      int           theCN,
                      double        theLinTol,
                      int&          theSigOrder,
                      LProp_Status& theTanStatus)
{
  if (theTanStatus == LProp_Undefined)
    return false;
  if (theTanStatus >= LProp_Defined)
    return true;

  const double aTolSq  = theLinTol * theLinTol;
  int          anOrder = 0;
  while (anOrder++ < 4)
  {
    if (theCN >= anOrder)
    {
      Vec aV;
      switch (anOrder)
      {
        case 1:
          aV = theProps.D1();
          break;
        case 2:
          aV = theProps.D2();
          break;
        case 3:
          aV = theProps.D3();
          break;
        default:
          theTanStatus = LProp_Undefined;
          return false;
      }
      if (aV.SquareMagnitude() > aTolSq)
      {
        theSigOrder  = anOrder;
        theTanStatus = LProp_Defined;
        return true;
      }
    }
    else
    {
      theTanStatus = LProp_Undefined;
      return false;
    }
  }
  return false;
}

//! Tangent: checks IsTangentDefined, then computes tangent direction.
//! @param[in,out] theProps     CLProps object
//! @param[in,out] theCurve     curve object
//! @param[in]     theU         current parameter
//! @param[in]     theDerivArr  derivative array
//! @param[in]     theRefPnt    reference point (for type deduction)
//! @param[in]     theSigOrder  order of first significant derivative
//! @param[out]    theDir       computed tangent direction
template <typename Access, typename Props, typename Curve, typename Vec, typename Pnt, typename Dir>
void Tangent(Props&     theProps,
             Curve&     theCurve,
             double     theU,
             const Vec* theDerivArr,
             const Pnt& theRefPnt,
             const int& theSigOrder,
             Dir&       theDir)
{
  if (!theProps.IsTangentDefined())
    throw LProp_NotDefined();
  ComputeTangent<Access>(theCurve, theU, theDerivArr, theRefPnt, theSigOrder, theDir);
}

//! Curvature: checks IsTangentDefined, returns RealLast if higher-order, else computes.
//! @param[in,out] theProps     CLProps object
//! @param[in]     theD1        first derivative
//! @param[in]     theD2        second derivative
//! @param[in]     theLinTol    linear tolerance
//! @param[in]     theSigOrder  order of first significant derivative
//! @param[out]    theCurvature curvature value field
//! @return computed curvature
template <typename Props, typename Vec>
double Curvature(Props&     theProps,
                 const Vec& theD1,
                 const Vec& theD2,
                 double     theLinTol,
                 const int& theSigOrder,
                 double&    theCurvature)
{
  const bool anIsDefined = theProps.IsTangentDefined();
  (void)anIsDefined;
  LProp_NotDefined_Raise_if(!anIsDefined, "CLProps::Curvature()");
  if (theSigOrder > 1)
    return RealLast();
  theCurvature = ComputeCurvature(theD1, theD2, theLinTol * theLinTol);
  return theCurvature;
}

//! Normal: checks curvature, then computes normal direction.
//! @param[in,out] theProps  CLProps object
//! @param[in]     theD1     first derivative
//! @param[in]     theD2     second derivative
//! @param[in]     theLinTol linear tolerance
//! @param[out]    theDir    computed normal direction
template <typename Props, typename Vec, typename Dir>
void Normal(Props& theProps, const Vec& theD1, const Vec& theD2, double theLinTol, Dir& theDir)
{
  const double aCurvature = theProps.Curvature();
  if (aCurvature == RealLast() || std::abs(aCurvature) <= theLinTol)
    throw LProp_NotDefined("CLProps::Normal(): Curvature is null or infinity");
  ComputeNormal(theD1, theD2, theDir);
}

//! CentreOfCurvature: checks curvature, then computes centre.
//! @param[in,out] theProps     CLProps object
//! @param[in]     thePnt       current point
//! @param[in]     theD1        first derivative
//! @param[in]     theD2        second derivative
//! @param[in]     theLinTol    linear tolerance
//! @param[in,out] theCurvature curvature value field
//! @param[out]    theCentre    computed centre of curvature
template <typename Props, typename Vec, typename Pnt>
void CentreOfCurvature(Props&     theProps,
                       const Pnt& thePnt,
                       const Vec& theD1,
                       const Vec& theD2,
                       double     theLinTol,
                       double&    theCurvature,
                       Pnt&       theCentre)
{
  if (std::abs(theProps.Curvature()) <= theLinTol)
    throw LProp_NotDefined();
  ComputeCentreOfCurvature(thePnt, theD1, theD2, theCurvature, theCentre);
}

} // namespace LProp_CurveUtils

#endif // _LProp_CurveUtils_HeaderFile