OCCT

// 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 GeomBndLib_SplineHelpers_HeaderFile
#define GeomBndLib_SplineHelpers_HeaderFile

#include <Bnd_Box.hxx>
#include <Bnd_Box2d.hxx>
#include <BSplCLib.hxx>
#include <ElCLib.hxx>
#include <NCollection_Array1.hxx>
#include <Precision.hxx>
#include <Standard_Handle.hxx>
#include <gp_Pnt.hxx>
#include <gp_Pnt2d.hxx>

#include <GeomBndLib_OptimizationHelpers.pxx>
#include <GeomBndLib_SamplingHelpers.pxx>

#include <algorithm>
#include <cmath>
#include <type_traits>

namespace GeomBndLib_SplineHelpers
{

template <typename TPoint>
struct PointOps;

template <>
struct PointOps<gp_Pnt>
{
  static gp_Pnt Mid(const gp_Pnt& theP1, const gp_Pnt& theP2)
  {
    return gp_Pnt((theP1.XYZ() + theP2.XYZ()) / 2.0);
  }
};

template <>
struct PointOps<gp_Pnt2d>
{
  static gp_Pnt2d Mid(const gp_Pnt2d& theP1, const gp_Pnt2d& theP2)
  {
    return gp_Pnt2d((theP1.XY() + theP2.XY()) / 2.0);
  }
};

//! Fill box with sampled points and return maximum deflection.
template <typename TBox, typename TAdaptor, typename TPoint>
double FillBox(TBox&           theBox,
               const TAdaptor& theCurve,
               const double    theFirst,
               const double    theLast,
               const int       theN)
{
  theBox.Add(theCurve.EvalD0(theFirst));
  double       aP   = theFirst;
  double       aTol = 0.0;
  const double aDP  = theLast - theFirst;
  if (std::abs(aDP) > Precision::PConfusion())
  {
    TPoint       aP1   = theCurve.EvalD0(theFirst);
    const double aStep = aDP / (2 * theN);
    for (int anIdx = 1; anIdx <= theN; ++anIdx)
    {
      aP += aStep;
      const TPoint aP2 = theCurve.EvalD0(aP);
      theBox.Add(aP2);
      aP += aStep;
      const TPoint aP3 = theCurve.EvalD0(aP);
      theBox.Add(aP3);
      const TPoint aPc = PointOps<TPoint>::Mid(aP1, aP3);
      aTol             = std::max(aTol, aPc.Distance(aP2));
      aP1              = aP3;
    }
  }
  else
  {
    theBox.Add(theCurve.EvalD0(theLast));
  }
  return aTol;
}

//! Reduce box using poles convex hull.
template <typename TPoint, typename TBox>
void ReduceSplineBox(const NCollection_Array1<TPoint>& thePoles,
                     const TBox&                       theOrigBox,
                     TBox&                             theReducedBox)
{
  TBox aPolesBox;
  for (int anIdx = thePoles.Lower(); anIdx <= thePoles.Upper(); ++anIdx)
  {
    aPolesBox.Add(thePoles.Value(anIdx));
  }

  if constexpr (std::is_same_v<TBox, Bnd_Box>)
  {
    const auto [aXMin, aXMax, aYMin, aYMax, aZMin, aZMax] = theOrigBox.Get();
    if (!aPolesBox.IsVoid())
    {
      const auto [aPXMin, aPXMax, aPYMin, aPYMax, aPZMin, aPZMax] = aPolesBox.Get();
      theReducedBox.Update(std::max(aXMin, aPXMin),
                           std::max(aYMin, aPYMin),
                           std::max(aZMin, aPZMin),
                           std::min(aXMax, aPXMax),
                           std::min(aYMax, aPYMax),
                           std::min(aZMax, aPZMax));
    }
    else
    {
      theReducedBox.Update(aXMin, aYMin, aZMin, aXMax, aYMax, aZMax);
    }
  }
  else
  {
    const auto [aXMin, aXMax, aYMin, aYMax] = theOrigBox.Get();
    if (!aPolesBox.IsVoid())
    {
      const auto [aPXMin, aPXMax, aPYMin, aPYMax] = aPolesBox.Get();
      theReducedBox.Update(std::max(aXMin, aPXMin),
                           std::max(aYMin, aPYMin),
                           std::min(aXMax, aPXMax),
                           std::min(aYMax, aPYMax));
    }
    else
    {
      theReducedBox.Update(aXMin, aYMin, aXMax, aYMax);
    }
  }
}

//! Compute the pole index range [theOutMinIdx, theOutMaxIdx] that covers the parameter
//! interval [theMin, theMax] in a B-spline defined by theKnots/theMults/theDegree.
//! For periodic curves theOutMaxIdx may exceed theMaxPoleIdx (wraps into [1, theMaxPoleIdx]).
inline void ComputePoleIndexRange(const NCollection_Array1<double>& theKnots,
                                  const NCollection_Array1<int>&    theMults,
                                  const int                         theDegree,
                                  const double                      theMin,
                                  const double                      theMax,
                                  const int                         theMaxPoleIdx,
                                  const bool                        theIsPeriodic,
                                  int&                              theOutMinIdx,
                                  int&                              theOutMaxIdx)
{
  BSplCLib::Hunt(theKnots, theMin, theOutMinIdx);
  theOutMinIdx = std::clamp(theOutMinIdx, theKnots.Lower(), theKnots.Upper());

  BSplCLib::Hunt(theKnots, theMax, theOutMaxIdx);
  theOutMaxIdx++;
  theOutMaxIdx          = std::clamp(theOutMaxIdx, theKnots.Lower(), theKnots.Upper());
  const int aMultiplier = theMults(theOutMaxIdx);

  theOutMinIdx = BSplCLib::PoleIndex(theDegree, theOutMinIdx, theIsPeriodic, theMults) + 1;
  theOutMinIdx = std::max(theOutMinIdx, 1);
  theOutMaxIdx = BSplCLib::PoleIndex(theDegree, theOutMaxIdx, theIsPeriodic, theMults) + 1;
  theOutMaxIdx += theDegree - aMultiplier;
  if (!theIsPeriodic)
  {
    theOutMaxIdx = std::min(theOutMaxIdx, theMaxPoleIdx);
  }
}

//! Reduce box using a selected range of poles [theMinIdx, theMaxIdx].
//! For periodic curves, indices beyond the array upper bound wrap around.
template <typename TPoint, typename TBox>
void ReduceSplineBox(const NCollection_Array1<TPoint>& thePoles,
                     const int                         theMinIdx,
                     const int                         theMaxIdx,
                     const TBox&                       theOrigBox,
                     TBox&                             theReducedBox)
{
  TBox      aPolesBox;
  const int aNbPoles = thePoles.Size();
  for (int anIdx = theMinIdx; anIdx <= theMaxIdx; ++anIdx)
  {
    int anI = anIdx;
    if (anI > thePoles.Upper())
    {
      anI -= aNbPoles;
    }
    aPolesBox.Add(thePoles.Value(anI));
  }

  if constexpr (std::is_same_v<TBox, Bnd_Box>)
  {
    const auto [aXMin, aXMax, aYMin, aYMax, aZMin, aZMax] = theOrigBox.Get();
    if (!aPolesBox.IsVoid())
    {
      const auto [aPXMin, aPXMax, aPYMin, aPYMax, aPZMin, aPZMax] = aPolesBox.Get();
      theReducedBox.Update(std::max(aXMin, aPXMin),
                           std::max(aYMin, aPYMin),
                           std::max(aZMin, aPZMin),
                           std::min(aXMax, aPXMax),
                           std::min(aYMax, aPYMax),
                           std::min(aZMax, aPZMax));
    }
    else
    {
      theReducedBox.Update(aXMin, aYMin, aZMin, aXMax, aYMax, aZMax);
    }
  }
  else
  {
    const auto [aXMin, aXMax, aYMin, aYMax] = theOrigBox.Get();
    if (!aPolesBox.IsVoid())
    {
      const auto [aPXMin, aPXMax, aPYMin, aPYMax] = aPolesBox.Get();
      theReducedBox.Update(std::max(aXMin, aPXMin),
                           std::max(aYMin, aPYMin),
                           std::min(aXMax, aPXMax),
                           std::min(aYMax, aPYMax));
    }
    else
    {
      theReducedBox.Update(aXMin, aYMin, aXMax, aYMax);
    }
  }
}

//! Compute bounding box for a Bezier curve (3D or 2D) over [theU1, theU2].
//! For sub-ranges, segments the curve and uses the segmented poles' convex hull.
//! For the full range, uses FillBox sampling + ReduceSplineBox.
template <typename TBezier, typename TAdaptor, typename TBox, typename TPoint>
TBox BezierCurveBox(const occ::handle<TBezier>& theGeom,
                    const double                theU1,
                    const double                theU2,
                    const double                theTol)
{
  // For a strict sub-range, segment the curve so ReduceSplineBox uses the exact sub-curve
  // poles. Without segmentation, the full poles' convex hull may exceed the sub-range extent.
  constexpr double aSub = Precision::PConfusion();
  if (theU1 - theGeom->FirstParameter() > aSub || theGeom->LastParameter() - theU2 > aSub)
  {
    occ::handle<TBezier> aBz = occ::down_cast<TBezier>(theGeom->Copy());
    aBz->Segment(theU1, theU2);
    TBox aBox;
    for (int anI = 1; anI <= aBz->NbPoles(); ++anI)
    {
      aBox.Add(aBz->Pole(anI));
    }
    aBox.Enlarge(theTol);
    return aBox;
  }

  constexpr double aWeakness = 1.5;
  TAdaptor         aGACurve(theGeom);
  TBox             aB1;
  const double     aTol =
    FillBox<TBox, TAdaptor, TPoint>(aB1, aGACurve, theU1, theU2, theGeom->Degree());
  aB1.Enlarge(aWeakness * aTol);
  TBox aBox;
  ReduceSplineBox(theGeom->Poles(), aB1, aBox);
  aBox.Enlarge(theTol);
  return aBox;
}

//! Compute bounding box for a BSpline curve (3D or 2D) over [theU1, theU2].
//! Handles periodic curves via Copy+Segment, sub-ranges via Segment+poles,
//! and full-range via knot-span FillBox + ComputePoleIndexRange + ReduceSplineBox.
template <typename TBSpline, typename TAdaptor, typename TBox, typename TPoint>
TBox BSplineCurveBox(const occ::handle<TBSpline>& theGeom,
                     const double                 theU1,
                     const double                 theU2,
                     const double                 theTol)
{
  constexpr double aWeakness = 1.5;
  const int        aDegree   = theGeom->Degree();

  if (theGeom->IsPeriodic())
  {
    // Periodic: use Segment to correctly unwrap the parameter range.
    occ::handle<TBSpline> aBs = occ::down_cast<TBSpline>(theGeom->Copy());
    double                aU1 = theU1, aU2 = theU2;
    ElCLib::AdjustPeriodic(aBs->FirstParameter(),
                           aBs->LastParameter(),
                           Precision::PConfusion(),
                           aU1,
                           aU2);

    const double aPeriod           = aBs->LastParameter() - aBs->FirstParameter();
    const double aDirectDiff       = std::abs(aU2 - aU1);
    const double aCrossPeriodDiff1 = std::abs(aU2 - aPeriod - aU1);
    const double aCrossPeriodDiff2 = std::abs(aU1 - aPeriod - aU2);
    double       aSegmentTol       = 2. * Precision::PConfusion();
    const double aMinDiff = std::min(aDirectDiff, std::min(aCrossPeriodDiff1, aCrossPeriodDiff2));
    if (aMinDiff < aSegmentTol)
    {
      aSegmentTol = aMinDiff * 0.01;
    }
    aBs->Segment(aU1, aU2, aSegmentTol);

    const int                         aK1    = aBs->FirstUKnotIndex();
    const int                         aK2    = aBs->LastUKnotIndex();
    const NCollection_Array1<double>& aKnots = aBs->Knots();
    TAdaptor                          aGACurve(aBs);
    double                            aFirst = aKnots(aK1), aLast;
    double                            aTol   = 0.0;
    TBox                              aB1;
    for (int aK = aK1 + 1; aK <= aK2; aK++)
    {
      aLast = aKnots(aK);
      aTol = std::max(FillBox<TBox, TAdaptor, TPoint>(aB1, aGACurve, aFirst, aLast, aDegree), aTol);
      aFirst = aLast;
    }
    TBox aBox;
    if (!aB1.IsVoid())
    {
      aB1.Enlarge(aWeakness * aTol);
      ReduceSplineBox(aBs->Poles(), aB1, aBox);
      aBox.Enlarge(theTol);
    }
    else
    {
      // Degenerate parameter range: evaluate endpoints directly.
      aBox.Add(aBs->EvalD0(aU1));
      aBox.Add(aBs->EvalD0(aU2));
      aBox.Enlarge(theTol);
    }
    return aBox;
  }

  // Non-periodic: traverse knot spans directly without Copy/Segment.
  const double aU1 = std::max(theU1, theGeom->FirstParameter());
  const double aU2 = std::min(theU2, theGeom->LastParameter());

  // For a strict sub-range, segment the curve so ReduceSplineBox uses the exact sub-curve
  // poles. Without segmentation, ComputePoleIndexRange includes adjacent-span poles whose
  // convex hull may exceed the actual sub-range extent.
  constexpr double aSub = Precision::PConfusion();
  if (aU1 - theGeom->FirstParameter() > aSub || theGeom->LastParameter() - aU2 > aSub)
  {
    occ::handle<TBSpline> aBs = occ::down_cast<TBSpline>(theGeom->Copy());
    aBs->Segment(aU1, aU2, 2.0 * Precision::PConfusion());
    TBox aBox;
    for (int anI = 1; anI <= aBs->NbPoles(); ++anI)
    {
      aBox.Add(aBs->Pole(anI));
    }
    aBox.Enlarge(theTol);
    return aBox;
  }

  const NCollection_Array1<double>& aKnots   = theGeom->Knots();
  const NCollection_Array1<int>&    aMults   = theGeom->Multiplicities();
  const int                         aNbPoles = theGeom->NbPoles();

  int aKMin = aKnots.Lower();
  BSplCLib::Hunt(aKnots, aU1, aKMin);
  aKMin = std::clamp(aKMin, aKnots.Lower(), aKnots.Upper() - 1);

  int aKMax = aKnots.Lower();
  BSplCLib::Hunt(aKnots, aU2, aKMax);
  aKMax = std::clamp(aKMax + 1, aKnots.Lower(), aKnots.Upper());

  TAdaptor aGACurve(theGeom);
  TBox     aB1;
  double   aTol   = 0.0;
  double   aFirst = aU1;
  for (int aK = aKMin + 1; aK <= aKMax; aK++)
  {
    const double aLast = (aK < aKMax) ? std::min(aU2, aKnots(aK)) : aU2;
    if (aLast > aFirst + Precision::PConfusion())
    {
      aTol = std::max(FillBox<TBox, TAdaptor, TPoint>(aB1, aGACurve, aFirst, aLast, aDegree), aTol);
    }
    aFirst = aLast;
    if (aFirst >= aU2 - Precision::PConfusion())
    {
      break;
    }
  }
  TBox aBox;
  if (!aB1.IsVoid())
  {
    aB1.Enlarge(aWeakness * aTol);
    int aPoleMin = 1, aPoleMax = aNbPoles;
    ComputePoleIndexRange(aKnots, aMults, aDegree, aU1, aU2, aNbPoles, false, aPoleMin, aPoleMax);
    ReduceSplineBox(theGeom->Poles(), aPoleMin, aPoleMax, aB1, aBox);
    aBox.Enlarge(theTol);
  }
  else
  {
    // Degenerate parameter range: evaluate endpoints directly.
    aBox.Add(aGACurve.EvalD0(aU1));
    aBox.Add(aGACurve.EvalD0(aU2));
    aBox.Enlarge(theTol);
  }
  return aBox;
}

//! Compute precise bounding box for a curve (3D or 2D) using sampling + numerical optimization.
//! Works for both Adaptor3d_Curve and Adaptor2d_Curve2d subtypes.
template <typename TAdaptor, typename TBox, typename TPoint>
TBox CurveBoxOptimal(TAdaptor&    theAdaptor,
                     const double theU1,
                     const double theU2,
                     const double theTol)
{
  using TCoord         = std::conditional_t<std::is_same_v<TPoint, gp_Pnt>, gp_XYZ, gp_XY>;
  constexpr int NDim   = std::is_same_v<TBox, Bnd_Box> ? 3 : 2;
  const int aNbSamples = GeomBndLib_SamplingHelpers::ComputeNbSamplesT(theAdaptor, theU1, theU2);

  double aCoordMin[NDim], aCoordMax[NDim], aDeflMax[NDim];
  for (int aK = 0; aK < NDim; ++aK)
  {
    aCoordMin[aK] = RealLast();
    aCoordMax[aK] = -RealLast();
    aDeflMax[aK]  = -RealLast();
  }

  const double               aDU  = (theU2 - theU1) / (aNbSamples - 1);
  const double               aDU2 = aDU / 2.;
  NCollection_Array1<TCoord> aPnts(1, aNbSamples);
  for (int anI = 1, aJ = 0; anI <= aNbSamples; anI++, aJ++)
  {
    const double aU = theU1 + aJ * aDU;
    const TPoint aP = theAdaptor.EvalD0(aU);
    if constexpr (std::is_same_v<TPoint, gp_Pnt>)
    {
      aPnts(anI) = aP.XYZ();
    }
    else
    {
      aPnts(anI) = aP.XY();
    }
    for (int aK = 0; aK < NDim; ++aK)
    {
      if (aCoordMin[aK] > aP.Coord(aK + 1))
      {
        aCoordMin[aK] = aP.Coord(aK + 1);
      }
      if (aCoordMax[aK] < aP.Coord(aK + 1))
      {
        aCoordMax[aK] = aP.Coord(aK + 1);
      }
    }
    if (anI > 1)
    {
      const TCoord aPm   = 0.5 * (aPnts(anI - 1) + aPnts(anI));
      const TPoint aPMid = theAdaptor.EvalD0(aU - aDU2);
      TCoord       aD;
      if constexpr (std::is_same_v<TPoint, gp_Pnt>)
      {
        aD = aPMid.XYZ() - aPm;
      }
      else
      {
        aD = aPMid.XY() - aPm;
      }
      for (int aK = 0; aK < NDim; ++aK)
      {
        if (aCoordMin[aK] > aPMid.Coord(aK + 1))
        {
          aCoordMin[aK] = aPMid.Coord(aK + 1);
        }
        if (aCoordMax[aK] < aPMid.Coord(aK + 1))
        {
          aCoordMax[aK] = aPMid.Coord(aK + 1);
        }
        const double aDiff = std::abs(aD.Coord(aK + 1));
        if (aDeflMax[aK] < aDiff)
        {
          aDeflMax[aK] = aDiff;
        }
      }
    }
  }

  const double anEps = std::max(theTol, Precision::Confusion());
  for (int aK = 0; aK < NDim; ++aK)
  {
    const double aDiff = aDeflMax[aK];
    if (aDiff <= anEps)
    {
      continue;
    }
    double aCMin = aCoordMin[aK];
    double aCMax = aCoordMax[aK];
    for (int anI = 1; anI <= aNbSamples; ++anI)
    {
      if (aPnts(anI).Coord(aK + 1) - aCMin < aDiff)
      {
        const double anUMin  = theU1 + std::max(0, anI - 2) * aDU;
        const double anUMax  = theU1 + std::min(aNbSamples - 1, anI) * aDU;
        const double aLocMin = GeomBndLib_OptimizationHelpers::AdjustExtrT(theAdaptor,
                                                                           anUMin,
                                                                           anUMax,
                                                                           aCMin,
                                                                           aK + 1,
                                                                           anEps,
                                                                           true);
        if (aLocMin < aCMin)
        {
          aCMin = aLocMin;
        }
      }
      else if (aCMax - aPnts(anI).Coord(aK + 1) < aDiff)
      {
        const double anUMin  = theU1 + std::max(0, anI - 2) * aDU;
        const double anUMax  = theU1 + std::min(aNbSamples - 1, anI) * aDU;
        const double aLocMax = GeomBndLib_OptimizationHelpers::AdjustExtrT(theAdaptor,
                                                                           anUMin,
                                                                           anUMax,
                                                                           aCMax,
                                                                           aK + 1,
                                                                           anEps,
                                                                           false);
        if (aLocMax > aCMax)
        {
          aCMax = aLocMax;
        }
      }
    }
    aCoordMin[aK] = aCMin;
    aCoordMax[aK] = aCMax;
  }

  TBox aBox;
  if constexpr (std::is_same_v<TBox, Bnd_Box>)
  {
    aBox.Add(gp_Pnt(aCoordMin[0], aCoordMin[1], aCoordMin[2]));
    aBox.Add(gp_Pnt(aCoordMax[0], aCoordMax[1], aCoordMax[2]));
  }
  else
  {
    aBox.Add(gp_Pnt2d(aCoordMin[0], aCoordMin[1]));
    aBox.Add(gp_Pnt2d(aCoordMax[0], aCoordMax[1]));
  }
  aBox.Enlarge(anEps);
  return aBox;
}

} // namespace GeomBndLib_SplineHelpers

#endif // GeomBndLib_SplineHelpers_HeaderFile