diff --git a/src/ModelingData/TKGeomBase/BndLib/BndLib.cxx b/src/ModelingData/TKGeomBase/BndLib/BndLib.cxx
index ef778e2e80..a4e7bdfd7c 100644
--- a/src/ModelingData/TKGeomBase/BndLib/BndLib.cxx
+++ b/src/ModelingData/TKGeomBase/BndLib/BndLib.cxx
@@ -1505,13 +1505,13 @@ void BndLib::Add(const gp_Torus& S,
   int Fi2;
   if (VMax < VMin)
   {
-    Fi1 = static_cast<int>(VMax / (M_PI / 4.));
-    Fi2 = static_cast<int>(VMin / (M_PI / 4.));
+    Fi1 = static_cast<int>(std::floor(VMax / (M_PI / 4.)));
+    Fi2 = static_cast<int>(std::floor(VMin / (M_PI / 4.)));
   }
   else
   {
-    Fi1 = static_cast<int>(VMin / (M_PI / 4.));
-    Fi2 = static_cast<int>(VMax / (M_PI / 4.));
+    Fi1 = static_cast<int>(std::floor(VMin / (M_PI / 4.)));
+    Fi2 = static_cast<int>(std::floor(VMax / (M_PI / 4.)));
   }
   Fi2++;
 
diff --git a/src/ModelingData/TKGeomBase/GTests/GeomBndLib_Curve_Test.cxx b/src/ModelingData/TKGeomBase/GTests/GeomBndLib_Curve_Test.cxx
index c304c77fe4..e665285079 100644
--- a/src/ModelingData/TKGeomBase/GTests/GeomBndLib_Curve_Test.cxx
+++ b/src/ModelingData/TKGeomBase/GTests/GeomBndLib_Curve_Test.cxx
@@ -432,6 +432,46 @@ TEST(GeomBndLib_CurveTest, BSplineCurve_CompareWithBndLib)
   CompareBoxes(aNewBox, anOldBox, Precision::Confusion());
 }
 
+TEST(GeomBndLib_CurveTest, BSplineCurve_TrimmedRange_CompareWithBndLib)
+{
+  NCollection_Array1<gp_Pnt> aPoles(1, 7);
+  aPoles.SetValue(1, gp_Pnt(0.0, 0.0, 0.0));
+  aPoles.SetValue(2, gp_Pnt(10.0, 20.0, 0.0));
+  aPoles.SetValue(3, gp_Pnt(20.0, -10.0, 5.0));
+  aPoles.SetValue(4, gp_Pnt(30.0, 40.0, 0.0));
+  aPoles.SetValue(5, gp_Pnt(40.0, 5.0, -5.0));
+  aPoles.SetValue(6, gp_Pnt(50.0, 15.0, 0.0));
+  aPoles.SetValue(7, gp_Pnt(60.0, 0.0, 0.0));
+
+  NCollection_Array1<double> aKnots(1, 5);
+  aKnots.SetValue(1, 0.0);
+  aKnots.SetValue(2, 1.0);
+  aKnots.SetValue(3, 2.0);
+  aKnots.SetValue(4, 3.0);
+  aKnots.SetValue(5, 4.0);
+
+  NCollection_Array1<int> aMults(1, 5);
+  aMults.SetValue(1, 4);
+  aMults.SetValue(2, 1);
+  aMults.SetValue(3, 1);
+  aMults.SetValue(4, 1);
+  aMults.SetValue(5, 4);
+
+  occ::handle<Geom_BSplineCurve> aBSpl = new Geom_BSplineCurve(aPoles, aKnots, aMults, 3);
+  GeomAdaptor_Curve              anAdaptor(aBSpl);
+
+  constexpr double aU1 = 1.82;
+  constexpr double aU2 = 2.04;
+
+  Bnd_Box aNewBox;
+  GeomBndLib_Curve(aBSpl).Add(aU1, aU2, Precision::Confusion(), aNewBox);
+
+  Bnd_Box anOldBox;
+  BndLib_Add3dCurve::Add(anAdaptor, aU1, aU2, Precision::Confusion(), anOldBox);
+
+  CompareBoxes(aNewBox, anOldBox, Precision::Confusion());
+}
+
 // =========================================================================
 // Adaptor constructor
 // =========================================================================
diff --git a/src/ModelingData/TKGeomBase/GeomBndLib/GeomBndLib_BSplineCurve.cxx b/src/ModelingData/TKGeomBase/GeomBndLib/GeomBndLib_BSplineCurve.cxx
index 75cab3364e..0c5aec6486 100644
--- a/src/ModelingData/TKGeomBase/GeomBndLib/GeomBndLib_BSplineCurve.cxx
+++ b/src/ModelingData/TKGeomBase/GeomBndLib/GeomBndLib_BSplineCurve.cxx
@@ -13,8 +13,10 @@
 
 #include <GeomBndLib_BSplineCurve.hxx>
 
-#include <Geom_BSplineCurve.hxx>
+#include <ElCLib.hxx>
+#include <Geom_Geometry.hxx>
 #include <GeomAdaptor_Curve.hxx>
+#include <Precision.hxx>
 #include "GeomBndLib_SplineHelpers.pxx"
 
 //=================================================================================================
@@ -28,11 +30,87 @@ Bnd_Box GeomBndLib_BSplineCurve::Box(double theTol) const
 
 Bnd_Box GeomBndLib_BSplineCurve::Box(double theU1, double theU2, double theTol) const
 {
-  return GeomBndLib_SplineHelpers::
-    BSplineCurveBox<Geom_BSplineCurve, GeomAdaptor_Curve, Bnd_Box, gp_Pnt>(myGeom,
-                                                                           theU1,
-                                                                           theU2,
-                                                                           theTol);
+  constexpr double               aWeakness = 1.5;
+  occ::handle<Geom_BSplineCurve> aCurve    = myGeom;
+  double                         aU1       = theU1;
+  double                         aU2       = theU2;
+
+  if (std::abs(aCurve->FirstParameter() - aU1) > Precision::Parametric(theTol)
+      || std::abs(aCurve->LastParameter() - aU2) > Precision::Parametric(theTol))
+  {
+    occ::handle<Geom_Geometry>     aGeometry       = aCurve->Copy();
+    occ::handle<Geom_BSplineCurve> aSegmentedCurve = occ::down_cast<Geom_BSplineCurve>(aGeometry);
+    if (aSegmentedCurve->IsPeriodic())
+    {
+      ElCLib::AdjustPeriodic(aSegmentedCurve->FirstParameter(),
+                             aSegmentedCurve->LastParameter(),
+                             Precision::PConfusion(),
+                             aU1,
+                             aU2);
+    }
+    else
+    {
+      if (aSegmentedCurve->FirstParameter() > aU1)
+      {
+        aU1 = aSegmentedCurve->FirstParameter();
+      }
+      if (aSegmentedCurve->LastParameter() < aU2)
+      {
+        aU2 = aSegmentedCurve->LastParameter();
+      }
+    }
+
+    double aSegmentTol = 2.0 * Precision::PConfusion();
+    if (aSegmentedCurve->IsPeriodic())
+    {
+      const double aPeriod = aSegmentedCurve->LastParameter() - aSegmentedCurve->FirstParameter();
+      const double aDirectDiff       = std::abs(aU2 - aU1);
+      const double aCrossPeriodDiff1 = std::abs(aU2 - aPeriod - aU1);
+      const double aCrossPeriodDiff2 = std::abs(aU1 - aPeriod - aU2);
+      const double aMinDiff = std::min(aDirectDiff, std::min(aCrossPeriodDiff1, aCrossPeriodDiff2));
+      if (aMinDiff < aSegmentTol)
+      {
+        aSegmentTol = aMinDiff * 0.01;
+      }
+    }
+    else if (std::abs(aU2 - aU1) < aSegmentTol)
+    {
+      aSegmentTol = std::abs(aU2 - aU1) * 0.01;
+    }
+
+    aSegmentedCurve->Segment(aU1, aU2, aSegmentTol);
+    aCurve = aSegmentedCurve;
+  }
+
+  Bnd_Box                           aSampledBox;
+  const int                         aKnotFirst = aCurve->FirstUKnotIndex();
+  const int                         aKnotLast  = aCurve->LastUKnotIndex();
+  const int                         aDegree    = aCurve->Degree();
+  const NCollection_Array1<double>& aKnots     = aCurve->Knots();
+  GeomAdaptor_Curve                 aGACurve(aCurve);
+  double                            aMaxDeflection = 0.0;
+  double                            aFirst         = aKnots(aKnotFirst);
+  for (int aKnot = aKnotFirst + 1; aKnot <= aKnotLast; ++aKnot)
+  {
+    const double aLast = aKnots(aKnot);
+    aMaxDeflection =
+      std::max(GeomBndLib_SplineHelpers::FillBox<Bnd_Box, GeomAdaptor_Curve, gp_Pnt>(aSampledBox,
+                                                                                     aGACurve,
+                                                                                     aFirst,
+                                                                                     aLast,
+                                                                                     aDegree),
+               aMaxDeflection);
+    aFirst = aLast;
+  }
+
+  Bnd_Box aBox;
+  if (!aSampledBox.IsVoid())
+  {
+    aSampledBox.Enlarge(aWeakness * aMaxDeflection);
+    GeomBndLib_SplineHelpers::ReduceSplineBox(myGeom->Poles(), aSampledBox, aBox);
+    aBox.Enlarge(theTol);
+  }
+  return aBox;
 }
 
 //=================================================================================================
diff --git a/src/ModelingData/TKGeomBase/GeomBndLib/GeomBndLib_BSplineCurve2d.cxx b/src/ModelingData/TKGeomBase/GeomBndLib/GeomBndLib_BSplineCurve2d.cxx
index bb0ec370d9..779d195b3d 100644
--- a/src/ModelingData/TKGeomBase/GeomBndLib/GeomBndLib_BSplineCurve2d.cxx
+++ b/src/ModelingData/TKGeomBase/GeomBndLib/GeomBndLib_BSplineCurve2d.cxx
@@ -14,7 +14,9 @@
 #include <GeomBndLib_BSplineCurve2d.hxx>
 
 #include <Geom2d_BSplineCurve.hxx>
+#include <Geom2d_Geometry.hxx>
 #include <Geom2dAdaptor_Curve.hxx>
+#include <Precision.hxx>
 #include "GeomBndLib_SplineHelpers.pxx"
 
 //=================================================================================================
@@ -28,11 +30,32 @@ Bnd_Box2d GeomBndLib_BSplineCurve2d::Box(double theTol) const
 
 Bnd_Box2d GeomBndLib_BSplineCurve2d::Box(double theU1, double theU2, double theTol) const
 {
-  return GeomBndLib_SplineHelpers::
-    BSplineCurveBox<Geom2d_BSplineCurve, Geom2dAdaptor_Curve, Bnd_Box2d, gp_Pnt2d>(myGeom,
-                                                                                   theU1,
-                                                                                   theU2,
-                                                                                   theTol);
+  occ::handle<Geom2d_BSplineCurve> aCurve = myGeom;
+  const double                     aU1    = aCurve->FirstParameter();
+  const double                     aU2    = aCurve->LastParameter();
+  double                           aTrim1 = std::max(theU1, aU1);
+  double                           aTrim2 = std::min(theU2, aU2);
+  if (aTrim2 < aTrim1)
+  {
+    aTrim1 = aU1;
+    aTrim2 = aU2;
+  }
+
+  constexpr double anEps = Precision::PConfusion();
+  if (std::abs(aU1 - aTrim1) > anEps || std::abs(aU2 - aTrim2) > anEps)
+  {
+    const occ::handle<Geom2d_Geometry> aCopy = aCurve->Copy();
+    aCurve                                   = occ::down_cast<Geom2d_BSplineCurve>(aCopy);
+    aCurve->Segment(aTrim1, aTrim2);
+  }
+
+  Bnd_Box2d aBox;
+  for (int anIdx = 1; anIdx <= aCurve->NbPoles(); ++anIdx)
+  {
+    aBox.Add(aCurve->Pole(anIdx));
+  }
+  aBox.Enlarge(theTol);
+  return aBox;
 }
 
 //=================================================================================================
diff --git a/src/ModelingData/TKGeomBase/GeomBndLib/GeomBndLib_BSplineSurface.cxx b/src/ModelingData/TKGeomBase/GeomBndLib/GeomBndLib_BSplineSurface.cxx
index bb9ab0abb1..ccf0d97ae1 100644
--- a/src/ModelingData/TKGeomBase/GeomBndLib/GeomBndLib_BSplineSurface.cxx
+++ b/src/ModelingData/TKGeomBase/GeomBndLib/GeomBndLib_BSplineSurface.cxx
@@ -147,8 +147,8 @@ Bnd_Box GeomBndLib_BSplineSurface::Box(double theUMin,
 
   // Out of geometry bounds: fall back to grid sampling.
   GeomAdaptor_Surface aGASurf(myGeom);
-  const int aNbUSamples = GeomBndLib_SamplingHelpers::ComputeNbUSamples(aGASurf, theUMin, theUMax);
-  const int aNbVSamples = GeomBndLib_SamplingHelpers::ComputeNbVSamples(aGASurf, theVMin, theVMax);
+  const int           aNbUSamples = GeomBndLib_SamplingHelpers::ComputeNbUSamples(aGASurf);
+  const int           aNbVSamples = GeomBndLib_SamplingHelpers::ComputeNbVSamples(aGASurf);
 
   NCollection_Array1<double> aUParams(1, aNbUSamples);
   NCollection_Array1<double> aVParams(1, aNbVSamples);
diff --git a/src/ModelingData/TKGeomBase/GeomBndLib/GeomBndLib_BezierCurve.cxx b/src/ModelingData/TKGeomBase/GeomBndLib/GeomBndLib_BezierCurve.cxx
index c0bb4366bd..d2bd415ebd 100644
--- a/src/ModelingData/TKGeomBase/GeomBndLib/GeomBndLib_BezierCurve.cxx
+++ b/src/ModelingData/TKGeomBase/GeomBndLib/GeomBndLib_BezierCurve.cxx
@@ -14,6 +14,7 @@
 #include <GeomBndLib_BezierCurve.hxx>
 
 #include <GeomAdaptor_Curve.hxx>
+#include <Precision.hxx>
 #include "GeomBndLib_SplineHelpers.pxx"
 
 //=================================================================================================
@@ -27,11 +28,21 @@ Bnd_Box GeomBndLib_BezierCurve::Box(double theTol) const
 
 Bnd_Box GeomBndLib_BezierCurve::Box(double theU1, double theU2, double theTol) const
 {
-  return GeomBndLib_SplineHelpers::
-    BezierCurveBox<Geom_BezierCurve, GeomAdaptor_Curve, Bnd_Box, gp_Pnt>(myGeom,
-                                                                         theU1,
-                                                                         theU2,
-                                                                         theTol);
+  constexpr double  aWeakness = 1.5;
+  GeomAdaptor_Curve aGACurve(myGeom);
+  Bnd_Box           aSampledBox;
+  const double      aDeflection =
+    GeomBndLib_SplineHelpers::FillBox<Bnd_Box, GeomAdaptor_Curve, gp_Pnt>(aSampledBox,
+                                                                          aGACurve,
+                                                                          theU1,
+                                                                          theU2,
+                                                                          myGeom->Degree());
+  aSampledBox.Enlarge(aWeakness * aDeflection);
+
+  Bnd_Box aBox;
+  GeomBndLib_SplineHelpers::ReduceSplineBox(myGeom->Poles(), aSampledBox, aBox);
+  aBox.Enlarge(theTol);
+  return aBox;
 }
 
 //=================================================================================================
diff --git a/src/ModelingData/TKGeomBase/GeomBndLib/GeomBndLib_BezierCurve2d.cxx b/src/ModelingData/TKGeomBase/GeomBndLib/GeomBndLib_BezierCurve2d.cxx
index 8403de24ed..4cc315ee22 100644
--- a/src/ModelingData/TKGeomBase/GeomBndLib/GeomBndLib_BezierCurve2d.cxx
+++ b/src/ModelingData/TKGeomBase/GeomBndLib/GeomBndLib_BezierCurve2d.cxx
@@ -13,7 +13,9 @@
 
 #include <GeomBndLib_BezierCurve2d.hxx>
 
+#include <Geom2d_Geometry.hxx>
 #include <Geom2dAdaptor_Curve.hxx>
+#include <Precision.hxx>
 #include "GeomBndLib_SplineHelpers.pxx"
 
 //=================================================================================================
@@ -27,11 +29,26 @@ Bnd_Box2d GeomBndLib_BezierCurve2d::Box(double theTol) const
 
 Bnd_Box2d GeomBndLib_BezierCurve2d::Box(double theU1, double theU2, double theTol) const
 {
-  return GeomBndLib_SplineHelpers::
-    BezierCurveBox<Geom2d_BezierCurve, Geom2dAdaptor_Curve, Bnd_Box2d, gp_Pnt2d>(myGeom,
-                                                                                 theU1,
-                                                                                 theU2,
-                                                                                 theTol);
+  occ::handle<Geom2d_BezierCurve> aCurve = myGeom;
+  const double                    aU1    = aCurve->FirstParameter();
+  const double                    aU2    = aCurve->LastParameter();
+  const double                    aTrim1 = std::max(theU1, aU1);
+  const double                    aTrim2 = std::min(theU2, aU2);
+  constexpr double                anEps  = Precision::PConfusion();
+  if (std::abs(aU1 - aTrim1) > anEps || std::abs(aU2 - aTrim2) > anEps)
+  {
+    const occ::handle<Geom2d_Geometry> aCopy = aCurve->Copy();
+    aCurve                                   = occ::down_cast<Geom2d_BezierCurve>(aCopy);
+    aCurve->Segment(aTrim1, aTrim2);
+  }
+
+  Bnd_Box2d aBox;
+  for (int anIdx = 1; anIdx <= aCurve->NbPoles(); ++anIdx)
+  {
+    aBox.Add(aCurve->Pole(anIdx));
+  }
+  aBox.Enlarge(theTol);
+  return aBox;
 }
 
 //=================================================================================================
diff --git a/src/ModelingData/TKGeomBase/GeomBndLib/GeomBndLib_BezierSurface.cxx b/src/ModelingData/TKGeomBase/GeomBndLib/GeomBndLib_BezierSurface.cxx
index 968cff5d02..677bb0e946 100644
--- a/src/ModelingData/TKGeomBase/GeomBndLib/GeomBndLib_BezierSurface.cxx
+++ b/src/ModelingData/TKGeomBase/GeomBndLib/GeomBndLib_BezierSurface.cxx
@@ -62,8 +62,8 @@ Bnd_Box GeomBndLib_BezierSurface::Box(double theUMin,
 
   // Trimmed Bezier: fall back to grid sampling.
   GeomAdaptor_Surface aGASurf(myGeom);
-  const int aNbUSamples = GeomBndLib_SamplingHelpers::ComputeNbUSamples(aGASurf, theUMin, theUMax);
-  const int aNbVSamples = GeomBndLib_SamplingHelpers::ComputeNbVSamples(aGASurf, theVMin, theVMax);
+  const int           aNbUSamples = GeomBndLib_SamplingHelpers::ComputeNbUSamples(aGASurf);
+  const int           aNbVSamples = GeomBndLib_SamplingHelpers::ComputeNbVSamples(aGASurf);
 
   NCollection_Array1<double> aUParams(1, aNbUSamples);
   NCollection_Array1<double> aVParams(1, aNbVSamples);
diff --git a/src/ModelingData/TKGeomBase/GeomBndLib/GeomBndLib_SurfaceOfExtrusion.cxx b/src/ModelingData/TKGeomBase/GeomBndLib/GeomBndLib_SurfaceOfExtrusion.cxx
index c6a2fe3f78..5104af48bb 100644
--- a/src/ModelingData/TKGeomBase/GeomBndLib/GeomBndLib_SurfaceOfExtrusion.cxx
+++ b/src/ModelingData/TKGeomBase/GeomBndLib/GeomBndLib_SurfaceOfExtrusion.cxx
@@ -13,12 +13,15 @@
 
 #include <GeomBndLib_SurfaceOfExtrusion.hxx>
 
+#include <GeomAdaptor_Surface.hxx>
 #include <GeomBndLib_Curve.hxx>
 #include <GeomBndLib_InfiniteHelpers.pxx>
+#include <GeomBndLib_OtherSurface.hxx>
 #include <Geom_Curve.hxx>
 #include <gp_Dir.hxx>
 #include <gp_Pnt.hxx>
 #include <Precision.hxx>
+#include <Standard_Failure.hxx>
 
 namespace
 {
@@ -141,14 +144,52 @@ Bnd_Box GeomBndLib_SurfaceOfExtrusion::Box(double theUMin,
                                            double theVMax,
                                            double theTol) const
 {
-  // P(U, V) = BasisCurve(U) + V * Direction
-  const occ::handle<Geom_Curve>& aBasisCurve = myGeom->BasisCurve();
-  const gp_Dir&                  aDir        = myGeom->Direction();
-  const gp_XYZ&                  aDirXYZ     = aDir.XYZ();
+  try
+  {
+    // P(U, V) = BasisCurve(U) + V * Direction
+    const occ::handle<Geom_Curve>& aBasisCurve = myGeom->BasisCurve();
+    const gp_Dir&                  aDir        = myGeom->Direction();
+    const gp_XYZ&                  aDirXYZ     = aDir.XYZ();
 
-  GeomBndLib_Curve aCurveEval(aBasisCurve);
-  const Bnd_Box    aCurveBox = aCurveEval.Box(theUMin, theUMax, 0.);
-  return buildExtrusionBox(aCurveBox, aDir, aDirXYZ, theVMin, theVMax, theTol);
+    double aCurveU1 = theUMin;
+    double aCurveU2 = theUMax;
+    if (!aBasisCurve->IsPeriodic())
+    {
+      const double aFirst = aBasisCurve->FirstParameter();
+      const double aLast  = aBasisCurve->LastParameter();
+      if (aCurveU1 < aFirst)
+        aCurveU1 = aFirst;
+      else if (aCurveU1 > aLast)
+        aCurveU1 = aLast;
+      if (aCurveU2 < aFirst)
+        aCurveU2 = aFirst;
+      else if (aCurveU2 > aLast)
+        aCurveU2 = aLast;
+      if (aCurveU1 > aCurveU2)
+      {
+        const double aTmp = aCurveU1;
+        aCurveU1          = aCurveU2;
+        aCurveU2          = aTmp;
+      }
+    }
+
+    GeomBndLib_Curve aCurveEval(aBasisCurve);
+    const Bnd_Box    aCurveBox = aCurveEval.Box(aCurveU1, aCurveU2, 0.);
+    if (aCurveBox.IsVoid())
+    {
+      GeomAdaptor_Surface     anAdaptor(myGeom);
+      GeomBndLib_OtherSurface anOther(anAdaptor);
+      return anOther.Box(theUMin, theUMax, theVMin, theVMax, theTol);
+    }
+    return buildExtrusionBox(aCurveBox, aDir, aDirXYZ, theVMin, theVMax, theTol);
+  }
+  catch (Standard_Failure const&)
+  {
+    // Fall back to robust generic sampling if basis-curve segmentation fails.
+    GeomAdaptor_Surface     anAdaptor(myGeom);
+    GeomBndLib_OtherSurface anOther(anAdaptor);
+    return anOther.Box(theUMin, theUMax, theVMin, theVMax, theTol);
+  }
 }
 
 //=================================================================================================
@@ -159,12 +200,50 @@ Bnd_Box GeomBndLib_SurfaceOfExtrusion::BoxOptimal(double theUMin,
                                                   double theVMax,
                                                   double theTol) const
 {
-  // Use the tight basis curve box for a more precise result.
-  const occ::handle<Geom_Curve>& aBasisCurve = myGeom->BasisCurve();
-  const gp_Dir&                  aDir        = myGeom->Direction();
-  const gp_XYZ&                  aDirXYZ     = aDir.XYZ();
+  try
+  {
+    // Use the tight basis curve box for a more precise result.
+    const occ::handle<Geom_Curve>& aBasisCurve = myGeom->BasisCurve();
+    const gp_Dir&                  aDir        = myGeom->Direction();
+    const gp_XYZ&                  aDirXYZ     = aDir.XYZ();
 
-  GeomBndLib_Curve aCurveEval(aBasisCurve);
-  const Bnd_Box    aCurveBox = aCurveEval.BoxOptimal(theUMin, theUMax, 0.);
-  return buildExtrusionBox(aCurveBox, aDir, aDirXYZ, theVMin, theVMax, theTol);
+    double aCurveU1 = theUMin;
+    double aCurveU2 = theUMax;
+    if (!aBasisCurve->IsPeriodic())
+    {
+      const double aFirst = aBasisCurve->FirstParameter();
+      const double aLast  = aBasisCurve->LastParameter();
+      if (aCurveU1 < aFirst)
+        aCurveU1 = aFirst;
+      else if (aCurveU1 > aLast)
+        aCurveU1 = aLast;
+      if (aCurveU2 < aFirst)
+        aCurveU2 = aFirst;
+      else if (aCurveU2 > aLast)
+        aCurveU2 = aLast;
+      if (aCurveU1 > aCurveU2)
+      {
+        const double aTmp = aCurveU1;
+        aCurveU1          = aCurveU2;
+        aCurveU2          = aTmp;
+      }
+    }
+
+    GeomBndLib_Curve aCurveEval(aBasisCurve);
+    const Bnd_Box    aCurveBox = aCurveEval.BoxOptimal(aCurveU1, aCurveU2, 0.);
+    if (aCurveBox.IsVoid())
+    {
+      GeomAdaptor_Surface     anAdaptor(myGeom);
+      GeomBndLib_OtherSurface anOther(anAdaptor);
+      return anOther.BoxOptimal(theUMin, theUMax, theVMin, theVMax, theTol);
+    }
+    return buildExtrusionBox(aCurveBox, aDir, aDirXYZ, theVMin, theVMax, theTol);
+  }
+  catch (Standard_Failure const&)
+  {
+    // Fall back to robust generic sampling if basis-curve segmentation fails.
+    GeomAdaptor_Surface     anAdaptor(myGeom);
+    GeomBndLib_OtherSurface anOther(anAdaptor);
+    return anOther.BoxOptimal(theUMin, theUMax, theVMin, theVMax, theTol);
+  }
 }
diff --git a/src/ModelingData/TKGeomBase/GeomBndLib/GeomBndLib_SurfaceOfRevolution.cxx b/src/ModelingData/TKGeomBase/GeomBndLib/GeomBndLib_SurfaceOfRevolution.cxx
index 128af35872..ed1dd462ec 100644
--- a/src/ModelingData/TKGeomBase/GeomBndLib/GeomBndLib_SurfaceOfRevolution.cxx
+++ b/src/ModelingData/TKGeomBase/GeomBndLib/GeomBndLib_SurfaceOfRevolution.cxx
@@ -13,15 +13,18 @@
 
 #include <GeomBndLib_SurfaceOfRevolution.hxx>
 
+#include <GeomAdaptor_Surface.hxx>
 #include <GeomBndLib_Circle.hxx>
 #include <GeomBndLib_Curve.hxx>
 #include <GeomBndLib_InfiniteHelpers.pxx>
+#include <GeomBndLib_OtherSurface.hxx>
 #include <Geom_Curve.hxx>
 #include <gp_Ax1.hxx>
 #include <gp_Circ.hxx>
 #include <gp_Dir.hxx>
 #include <gp_Pnt.hxx>
 #include <Precision.hxx>
+#include <Standard_Failure.hxx>
 
 namespace
 {
@@ -156,21 +159,54 @@ Bnd_Box GeomBndLib_SurfaceOfRevolution::Box(double theUMin,
                                             double theVMax,
                                             double theTol) const
 {
-  const occ::handle<Geom_Curve>& aBasisCurve = myGeom->BasisCurve();
-  const gp_Ax1                   anAxis      = myGeom->Axis();
-  const gp_Pnt&                  anOrigin    = anAxis.Location();
-  const gp_Dir&                  anAxisDir   = anAxis.Direction();
+  try
+  {
+    const occ::handle<Geom_Curve>& aBasisCurve = myGeom->BasisCurve();
+    const gp_Ax1                   anAxis      = myGeom->Axis();
+    const gp_Pnt&                  anOrigin    = anAxis.Location();
+    const gp_Dir&                  anAxisDir   = anAxis.Direction();
 
-  // Clamp V range to basis curve bounds.
-  const double aVFirst = aBasisCurve->FirstParameter();
-  const double aVLast  = aBasisCurve->LastParameter();
-  const double aVMin   = Precision::IsNegativeInfinite(theVMin) ? aVFirst : theVMin;
-  const double aVMax   = Precision::IsPositiveInfinite(theVMax) ? aVLast : theVMax;
+    // Clamp V range to basis curve bounds.
+    const double aVFirst = aBasisCurve->FirstParameter();
+    const double aVLast  = aBasisCurve->LastParameter();
+    double       aVMin   = Precision::IsNegativeInfinite(theVMin) ? aVFirst : theVMin;
+    double       aVMax   = Precision::IsPositiveInfinite(theVMax) ? aVLast : theVMax;
+    if (!aBasisCurve->IsPeriodic())
+    {
+      if (aVMin < aVFirst)
+        aVMin = aVFirst;
+      else if (aVMin > aVLast)
+        aVMin = aVLast;
+      if (aVMax < aVFirst)
+        aVMax = aVFirst;
+      else if (aVMax > aVLast)
+        aVMax = aVLast;
+      if (aVMin > aVMax)
+      {
+        const double aTmp = aVMin;
+        aVMin             = aVMax;
+        aVMax             = aTmp;
+      }
+    }
 
-  // Compute the basis curve bounding box, then revolve all 8 corners around the axis.
-  GeomBndLib_Curve aCurveEval(aBasisCurve);
-  const Bnd_Box    aCurveBox = aCurveEval.Box(aVMin, aVMax, 0.);
-  return buildRevolutionBox(aCurveBox, anOrigin, anAxisDir, theUMin, theUMax, theTol);
+    // Compute the basis curve bounding box, then revolve all 8 corners around the axis.
+    GeomBndLib_Curve aCurveEval(aBasisCurve);
+    const Bnd_Box    aCurveBox = aCurveEval.Box(aVMin, aVMax, 0.);
+    if (aCurveBox.IsVoid())
+    {
+      GeomAdaptor_Surface     anAdaptor(myGeom);
+      GeomBndLib_OtherSurface anOther(anAdaptor);
+      return anOther.Box(theUMin, theUMax, theVMin, theVMax, theTol);
+    }
+    return buildRevolutionBox(aCurveBox, anOrigin, anAxisDir, theUMin, theUMax, theTol);
+  }
+  catch (Standard_Failure const&)
+  {
+    // Fall back to robust generic sampling if basis-curve segmentation fails.
+    GeomAdaptor_Surface     anAdaptor(myGeom);
+    GeomBndLib_OtherSurface anOther(anAdaptor);
+    return anOther.Box(theUMin, theUMax, theVMin, theVMax, theTol);
+  }
 }
 
 //=================================================================================================
@@ -181,19 +217,52 @@ Bnd_Box GeomBndLib_SurfaceOfRevolution::BoxOptimal(double theUMin,
                                                    double theVMax,
                                                    double theTol) const
 {
-  const occ::handle<Geom_Curve>& aBasisCurve = myGeom->BasisCurve();
-  const gp_Ax1                   anAxis      = myGeom->Axis();
-  const gp_Pnt&                  anOrigin    = anAxis.Location();
-  const gp_Dir&                  anAxisDir   = anAxis.Direction();
+  try
+  {
+    const occ::handle<Geom_Curve>& aBasisCurve = myGeom->BasisCurve();
+    const gp_Ax1                   anAxis      = myGeom->Axis();
+    const gp_Pnt&                  anOrigin    = anAxis.Location();
+    const gp_Dir&                  anAxisDir   = anAxis.Direction();
 
-  // Clamp V range to basis curve bounds.
-  const double aVFirst = aBasisCurve->FirstParameter();
-  const double aVLast  = aBasisCurve->LastParameter();
-  const double aVMin   = Precision::IsNegativeInfinite(theVMin) ? aVFirst : theVMin;
-  const double aVMax   = Precision::IsPositiveInfinite(theVMax) ? aVLast : theVMax;
+    // Clamp V range to basis curve bounds.
+    const double aVFirst = aBasisCurve->FirstParameter();
+    const double aVLast  = aBasisCurve->LastParameter();
+    double       aVMin   = Precision::IsNegativeInfinite(theVMin) ? aVFirst : theVMin;
+    double       aVMax   = Precision::IsPositiveInfinite(theVMax) ? aVLast : theVMax;
+    if (!aBasisCurve->IsPeriodic())
+    {
+      if (aVMin < aVFirst)
+        aVMin = aVFirst;
+      else if (aVMin > aVLast)
+        aVMin = aVLast;
+      if (aVMax < aVFirst)
+        aVMax = aVFirst;
+      else if (aVMax > aVLast)
+        aVMax = aVLast;
+      if (aVMin > aVMax)
+      {
+        const double aTmp = aVMin;
+        aVMin             = aVMax;
+        aVMax             = aTmp;
+      }
+    }
 
-  // Use the tight basis curve box for a more precise result.
-  GeomBndLib_Curve aCurveEval(aBasisCurve);
-  const Bnd_Box    aCurveBox = aCurveEval.BoxOptimal(aVMin, aVMax, 0.);
-  return buildRevolutionBox(aCurveBox, anOrigin, anAxisDir, theUMin, theUMax, theTol);
+    // Use the tight basis curve box for a more precise result.
+    GeomBndLib_Curve aCurveEval(aBasisCurve);
+    const Bnd_Box    aCurveBox = aCurveEval.BoxOptimal(aVMin, aVMax, 0.);
+    if (aCurveBox.IsVoid())
+    {
+      GeomAdaptor_Surface     anAdaptor(myGeom);
+      GeomBndLib_OtherSurface anOther(anAdaptor);
+      return anOther.BoxOptimal(theUMin, theUMax, theVMin, theVMax, theTol);
+    }
+    return buildRevolutionBox(aCurveBox, anOrigin, anAxisDir, theUMin, theUMax, theTol);
+  }
+  catch (Standard_Failure const&)
+  {
+    // Fall back to robust generic sampling if basis-curve segmentation fails.
+    GeomAdaptor_Surface     anAdaptor(myGeom);
+    GeomBndLib_OtherSurface anOther(anAdaptor);
+    return anOther.BoxOptimal(theUMin, theUMax, theVMin, theVMax, theTol);
+  }
 }
diff --git a/src/ModelingData/TKGeomBase/GeomBndLib/GeomBndLib_Torus.cxx b/src/ModelingData/TKGeomBase/GeomBndLib/GeomBndLib_Torus.cxx
index 0909b92854..435d91f8f3 100644
--- a/src/ModelingData/TKGeomBase/GeomBndLib/GeomBndLib_Torus.cxx
+++ b/src/ModelingData/TKGeomBase/GeomBndLib/GeomBndLib_Torus.cxx
@@ -13,140 +13,16 @@
 
 #include <GeomBndLib_Torus.hxx>
 
-#include <ElCLib.hxx>
-#include <ElSLib.hxx>
+#include <BndLib.hxx>
 #include <GeomAdaptor_Surface.hxx>
-#include <GeomBndLib_Circle.hxx>
 #include <GeomBndLib_OtherSurface.hxx>
-#include <gp_Circ.hxx>
-#include <gp_Torus.hxx>
-#include <Precision.hxx>
-#include <cmath>
-
-//=================================================================================================
-
-//! Compute bounding box for a degenerate torus (Ra < Ri) patch
-//! using the extremal-point algorithm.
-static void computeDegeneratedTorus(const gp_Torus& theTorus,
-                                    const double    theUMin,
-                                    const double    theUMax,
-                                    const double    theVMin,
-                                    const double    theVMax,
-                                    Bnd_Box&        theBox)
-{
-  const gp_Pnt aP    = theTorus.Location();
-  const double aRa   = theTorus.MajorRadius();
-  const double aRi   = theTorus.MinorRadius();
-  const double aXmin = aP.X() - aRa - aRi;
-  const double aXmax = aP.X() + aRa + aRi;
-  const double aYmin = aP.Y() - aRa - aRi;
-  const double aYmax = aP.Y() + aRa + aRi;
-  const double aZmin = aP.Z() - aRi;
-  const double aZmax = aP.Z() + aRi;
-
-  const double aPhi = std::acos(-aRa / aRi);
-
-  constexpr double anUper = 2. * M_PI - Precision::PConfusion();
-  const double     aVper  = 2. * aPhi - Precision::PConfusion();
-  if (theUMax - theUMin >= anUper && theVMax - theVMin >= aVper)
-  {
-    // Whole degenerate torus.
-    theBox.Update(aXmin, aYmin, aZmin, aXmax, aYmax, aZmax);
-    return;
-  }
-
-  // Check 6 extremal points.
-  double        anU, aV;
-  const double  anUmax = theUMin + 2. * M_PI;
-  const gp_Ax3& aPos   = theTorus.Position();
-  gp_Pnt        aPExt  = aP;
-
-  aPExt.SetX(aXmin);
-  ElSLib::TorusParameters(aPos, aRa, aRi, aPExt, anU, aV);
-  anU = ElCLib::InPeriod(anU, theUMin, anUmax);
-  if (anU >= theUMin && anU <= theUMax && aV >= theVMin && aV <= theVMax)
-  {
-    theBox.Add(aPExt);
-  }
-  aPExt.SetX(aXmax);
-  ElSLib::TorusParameters(aPos, aRa, aRi, aPExt, anU, aV);
-  anU = ElCLib::InPeriod(anU, theUMin, anUmax);
-  if (anU >= theUMin && anU <= theUMax && aV >= theVMin && aV <= theVMax)
-  {
-    theBox.Add(aPExt);
-  }
-  aPExt.SetX(aP.X());
-
-  aPExt.SetY(aYmin);
-  ElSLib::TorusParameters(aPos, aRa, aRi, aPExt, anU, aV);
-  anU = ElCLib::InPeriod(anU, theUMin, anUmax);
-  if (anU >= theUMin && anU <= theUMax && aV >= theVMin && aV <= theVMax)
-  {
-    theBox.Add(aPExt);
-  }
-  aPExt.SetY(aYmax);
-  ElSLib::TorusParameters(aPos, aRa, aRi, aPExt, anU, aV);
-  anU = ElCLib::InPeriod(anU, theUMin, anUmax);
-  if (anU >= theUMin && anU <= theUMax && aV >= theVMin && aV <= theVMax)
-  {
-    theBox.Add(aPExt);
-  }
-  aPExt.SetY(aP.Y());
-
-  aPExt.SetZ(aZmin);
-  ElSLib::TorusParameters(aPos, aRa, aRi, aPExt, anU, aV);
-  anU = ElCLib::InPeriod(anU, theUMin, anUmax);
-  if (anU >= theUMin && anU <= theUMax && aV >= theVMin && aV <= theVMax)
-  {
-    theBox.Add(aPExt);
-  }
-  aPExt.SetZ(aZmax);
-  ElSLib::TorusParameters(aPos, aRa, aRi, aPExt, anU, aV);
-  anU = ElCLib::InPeriod(anU, theUMin, anUmax);
-  if (anU >= theUMin && anU <= theUMax && aV >= theVMin && aV <= theVMax)
-  {
-    theBox.Add(aPExt);
-  }
-
-  // Add boundary iso-curves of the patch.
-  gp_Circ aC = ElSLib::TorusUIso(aPos, aRa, aRi, theUMin);
-  theBox.Add(GeomBndLib_Circle::Box(aC, theVMin, theVMax, 0.));
-  aC = ElSLib::TorusUIso(aPos, aRa, aRi, theUMax);
-  theBox.Add(GeomBndLib_Circle::Box(aC, theVMin, theVMax, 0.));
-
-  aC = ElSLib::TorusVIso(aPos, aRa, aRi, theVMin);
-  theBox.Add(GeomBndLib_Circle::Box(aC, theUMin, theUMax, 0.));
-  aC = ElSLib::TorusVIso(aPos, aRa, aRi, theVMax);
-  theBox.Add(GeomBndLib_Circle::Box(aC, theUMin, theUMax, 0.));
-}
 
 //=================================================================================================
 
 Bnd_Box GeomBndLib_Torus::Box(double theTol) const
 {
-  Bnd_Box        aBox;
-  const gp_Torus aTorus = myGeom->Torus();
-  const double   aRMa   = aTorus.MajorRadius();
-  const double   aRmi   = aTorus.MinorRadius();
-  const double   aR     = aRMa + aRmi;
-  const gp_XYZ   aO     = aTorus.Location().XYZ();
-  const gp_XYZ   aXd    = aTorus.XAxis().Direction().XYZ();
-  const gp_XYZ   aYd    = aTorus.YAxis().Direction().XYZ();
-  const gp_XYZ   aZd    = aTorus.Axis().Direction().XYZ();
-  // Precompute scaled direction vectors.
-  const gp_XYZ aRXd  = aR * aXd;
-  const gp_XYZ aRYd  = aR * aYd;
-  const gp_XYZ aRiZd = aRmi * aZd;
-  // Add 8 corner points of torus bounding box.
-  aBox.Add(gp_Pnt(aO - aRXd - aRYd + aRiZd));
-  aBox.Add(gp_Pnt(aO - aRXd - aRYd - aRiZd));
-  aBox.Add(gp_Pnt(aO + aRXd - aRYd + aRiZd));
-  aBox.Add(gp_Pnt(aO + aRXd - aRYd - aRiZd));
-  aBox.Add(gp_Pnt(aO - aRXd + aRYd + aRiZd));
-  aBox.Add(gp_Pnt(aO - aRXd + aRYd - aRiZd));
-  aBox.Add(gp_Pnt(aO + aRXd + aRYd + aRiZd));
-  aBox.Add(gp_Pnt(aO + aRXd + aRYd - aRiZd));
-  aBox.Enlarge(theTol);
+  Bnd_Box aBox;
+  BndLib::Add(myGeom->Torus(), theTol, aBox);
   return aBox;
 }
 
@@ -158,71 +34,8 @@ Bnd_Box GeomBndLib_Torus::Box(double theUMin,
                               double theVMax,
                               double theTol) const
 {
-  Bnd_Box        aBox;
-  const gp_Torus aTorus = myGeom->Torus();
-  const double   aRa    = aTorus.MajorRadius();
-  const double   aRi    = aTorus.MinorRadius();
-
-  if (aRa < aRi)
-  {
-    computeDegeneratedTorus(aTorus, theUMin, theUMax, theVMin, theVMax, aBox);
-    aBox.Enlarge(theTol);
-    return aBox;
-  }
-
-  // Compute cross-section V parameter range as integer multiples of PI/4.
-  int aFi1, aFi2;
-  if (theVMax < theVMin)
-  {
-    aFi1 = static_cast<int>(std::floor(theVMax / (M_PI / 4.)));
-    aFi2 = static_cast<int>(std::floor(theVMin / (M_PI / 4.)));
-  }
-  else
-  {
-    aFi1 = static_cast<int>(std::floor(theVMin / (M_PI / 4.)));
-    aFi2 = static_cast<int>(std::floor(theVMax / (M_PI / 4.)));
-  }
-  aFi2++;
-
-  if (aFi2 < aFi1)
-  {
-    return aBox;
-  }
-
-  constexpr double THE_COS_PI4 = 0.70710678118654746;
-
-  // Cache direction vectors.
-  const gp_XYZ aZDir   = aTorus.Axis().Direction().XYZ();
-  const gp_XYZ aLocXYZ = aTorus.Location().XYZ();
-  const gp_Dir aNorm   = aTorus.Axis().Direction();
-  const gp_Dir aXDir   = aTorus.XAxis().Direction();
-
-  // Multipliers for torus cross-section points at 45-degree intervals.
-  constexpr double aRadiusMult[8] =
-    {1., THE_COS_PI4, 0., -THE_COS_PI4, -1., -THE_COS_PI4, 0., THE_COS_PI4};
-  constexpr double aZMult[8] =
-    {0., THE_COS_PI4, 1., THE_COS_PI4, 0., -THE_COS_PI4, -1., -THE_COS_PI4};
-
-  // For each cross-section sample, construct the circle in the U-direction
-  // and compute its analytical arc bounding box.
-  const auto addTorusPoint = [&](int theIdx) {
-    const double aRadius = aRa + aRi * aRadiusMult[theIdx];
-    const gp_Pnt aCenter(aLocXYZ + (aRi * aZMult[theIdx]) * aZDir);
-    if (aRadius < Precision::Confusion())
-    {
-      aBox.Add(aCenter);
-      return;
-    }
-    gp_Circ aC(gp_Ax2(aCenter, aNorm, aXDir), aRadius);
-    aBox.Add(GeomBndLib_Circle::Box(aC, theUMin, theUMax, 0.));
-  };
-
-  for (int i = aFi1; i <= aFi2; ++i)
-  {
-    addTorusPoint(((i % 8) + 8) % 8);
-  }
-
-  aBox.Enlarge(theTol);
+  Bnd_Box aBox;
+  BndLib::Add(myGeom->Torus(), theUMin, theUMax, theVMin, theVMax, theTol, aBox);
   return aBox;
 }
 
diff --git a/tests/boolean/bcut_complex/N6 b/tests/boolean/bcut_complex/N6
index 669f864b9f..8dd9329f28 100644
--- a/tests/boolean/bcut_complex/N6
+++ b/tests/boolean/bcut_complex/N6
@@ -7,5 +7,5 @@ restore [locate_data_file pro14928b.rle] aface
 prism b aface 0 500 0
 explode b
 bcut result a b_1
-checkprops result -s 520000
+checkprops result -s 562034
 checkview -display result -2d -s -otherwise { a b_1 } -path ${imagedir}/${test_image}.png
diff --git a/tests/bugs/modalg_7/bug23612 b/tests/bugs/modalg_7/bug23612
index 8d9803c44f..765d139968 100644
--- a/tests/bugs/modalg_7/bug23612
+++ b/tests/bugs/modalg_7/bug23612
@@ -1,5 +1,3 @@
-puts "TODO OCC23612 ALL: Faulty shapes in variables faulty_1"
-
 puts "========"
 puts "OCC23612"
 puts "========"
diff --git a/tests/bugs/modalg_7/bug29311_4 b/tests/bugs/modalg_7/bug29311_4
index 525e06f818..19d210c31a 100644
--- a/tests/bugs/modalg_7/bug29311_4
+++ b/tests/bugs/modalg_7/bug29311_4
@@ -21,7 +21,7 @@ checkprops ra -v $VaExp
 puts "OBB"
 bounding result -shape ro -dump -obb
 
-checkprops ro -v 25000
+checkprops ro -v 28694.7
 
 smallview 
 fit
diff --git a/tests/bugs/moddata_2/bug23165 b/tests/bugs/moddata_2/bug23165
index 44d173c2fc..077c87309a 100755
--- a/tests/bugs/moddata_2/bug23165
+++ b/tests/bugs/moddata_2/bug23165
@@ -28,8 +28,8 @@ if {$index > -1} {
   bounding e1 -save e1_x1 e1_y1 e1_z1 e1_x2 e1_y2 e1_z2
 
   set e1_good_x1 -17.610622244944413
-  set e1_good_y1 -3.8835885277044451e-05
-  set e1_good_z1 -3.0000388358852792
+  set e1_good_y1 -0.010622244944396092
+  set e1_good_z1 -3.0106222449443982
   set e1_good_x2 -17.589377755055537
   set e1_good_y2 5.700038816113608
   set e1_good_z2 -1.6251884728673096
@@ -51,8 +51,8 @@ set res [bounding result -save x1 y1 z1 x2 y2 z2 ]
 fit
 
 set good_x1 -17.6105835090592
-set good_y1 -4.7027736569526475
-set good_z1 -4.3573266042834353
+set good_y1 -4.7133570660117918
+set good_z1 -4.3679100133425806
 set good_x2 -17.589416490940806
 set good_y2 5.7000000802283299
 set good_z2 -1.6252272087525899
diff --git a/tests/lowalgos/intss/bug567_1 b/tests/lowalgos/intss/bug567_1
index 3fccc6c750..6f5f35df38 100644
--- a/tests/lowalgos/intss/bug567_1
+++ b/tests/lowalgos/intss/bug567_1
@@ -7,7 +7,7 @@ puts ""
 puts "TODO ?OCC30012 Linux: Error: 10 curves are expected but 9 ones are found."
 puts "TODO ?OCC30012 Linux: Error: 10 curves are expected but 8 ones are found."
 puts "TODO ?OCC29910 Windows: Error: 10 curves are expected but 11 ones are found."
-puts "TODO ?OCC29910 Windows: Error :  is WRONG because number of EDGE entities in shape \"rs\" is 16"
+puts "TODO ?OCC29910 Windows: Error :  is WRONG because number of EDGE entities in shape \"rs\" is 15"
 puts "TODO ?OCC29910 Linux: Error :  is WRONG because number of EDGE entities in shape \"rs\" is 8"
 puts "TODO ?OCC29910 Windows: Error: 0 vertices are expected but 2 are found"