Modeling - Optimize BndLib and add GTests (#856)

- Fix negative modulo in torus bounding box computation
- Fix hyperbola extrema loop early break condition
- Fix OpenMin/OpenMax direction sign for infinite bounds
- Remove dead code in cone bounding computation
- Replace sqrt/log with more efficient computations in hyperbola bounds
- General code cleanup and modernization
This commit is contained in:
Pasukhin Dmitry
2025-11-27 09:22:44 +00:00
committed by GitHub
parent 1663c65625
commit c04f5e0b4c
7 changed files with 1733 additions and 428 deletions

View File

@@ -42,6 +42,12 @@ static Standard_Integer ComputeBox(const gp_Hypr& aHypr,
namespace namespace
{ {
//! Cosine of M_PI/8 (22.5 degrees) - used for 8-point polygon approximation.
constexpr Standard_Real THE_COS_PI8 = 0.92387953251128674;
//! Cosine (and sine) of M_PI/4 (45 degrees) - used for diagonal points.
constexpr Standard_Real THE_COS_PI4 = 0.70710678118654746;
//! Compute method //! Compute method
template <class PointType, class BndBoxType> template <class PointType, class BndBoxType>
void Compute(const Standard_Real theP1, void Compute(const Standard_Real theP1,
@@ -74,19 +80,11 @@ void Compute(const Standard_Real theP1,
} }
else else
{ {
// Normalize aTeta1 to [0, 2*PI) range
aTeta1 = std::fmod(aTeta1, 2. * M_PI);
if (aTeta1 < 0.) if (aTeta1 < 0.)
{ {
do aTeta1 += 2. * M_PI;
{
aTeta1 += 2. * M_PI;
} while (aTeta1 < 0.);
}
else if (aTeta1 > 2. * M_PI)
{
do
{
aTeta1 -= 2. * M_PI;
} while (aTeta1 > 2. * M_PI);
} }
aTeta2 = aTeta1 + aDelta; aTeta2 = aTeta1 + aDelta;
} }
@@ -104,149 +102,76 @@ void Compute(const Standard_Real theP1,
if (aDelta > M_PI / 8.) if (aDelta > M_PI / 8.)
{ {
// Main radiuses to take into account only 8 points (/cos(Pi/8.)) // Main radiuses to take into account only 8 points (/cos(Pi/8.))
aRam = theRa / 0.92387953251128674; aRam = theRa / THE_COS_PI8;
aRbm = theRb / 0.92387953251128674; aRbm = theRb / THE_COS_PI8;
} }
else else
{ {
// Main radiuses to take into account the arrow // Main radiuses to take into account the arrow
Standard_Real aTc = cos(aDelta / 2); Standard_Real aTc = std::cos(aDelta / 2);
aRam = theRa / aTc; aRam = theRa / aTc;
aRbm = theRb / aTc; aRbm = theRb / aTc;
} }
theB.Add(PointType(theO.Coord() + aRam * aCn1 * theXd.Coord() + aRbm * aSn1 * theYd.Coord())); theB.Add(PointType(theO.Coord() + aRam * aCn1 * theXd.Coord() + aRbm * aSn1 * theYd.Coord()));
theB.Add(PointType(theO.Coord() + aRam * aCn2 * theXd.Coord() + aRbm * aSn2 * theYd.Coord())); theB.Add(PointType(theO.Coord() + aRam * aCn2 * theXd.Coord() + aRbm * aSn2 * theYd.Coord()));
// cos or sin M_PI/4. // X and Y multipliers for 8 polygon points at 45-degree intervals (0, 45, 90, ..., 315 degrees).
#define PI4 0.70710678118654746 // Point i corresponds to angle i * 45 degrees.
constexpr Standard_Real aXMult[8] =
{1., THE_COS_PI4, 0., -THE_COS_PI4, -1., -THE_COS_PI4, 0., THE_COS_PI4};
constexpr Standard_Real aYMult[8] =
{0., THE_COS_PI4, 1., THE_COS_PI4, 0., -THE_COS_PI4, -1., -THE_COS_PI4};
// 8 points of the polygon // Lambda to add polygon point by index (0-7).
#define addPoint0 theB.Add(PointType(theO.Coord() + aRam * theXd.Coord())) const auto addPoint = [&](Standard_Integer theIdx) {
#define addPoint1 \ theB.Add(PointType(theO.Coord() + aRam * aXMult[theIdx] * theXd.Coord()
theB.Add(PointType(theO.Coord() + aRam * PI4 * theXd.Coord() + aRbm * PI4 * theYd.Coord())) + aRbm * aYMult[theIdx] * theYd.Coord()));
#define addPoint2 theB.Add(PointType(theO.Coord() + aRbm * theYd.Coord())) };
#define addPoint3 \
theB.Add(PointType(theO.Coord() - aRam * PI4 * theXd.Coord() + aRbm * PI4 * theYd.Coord()))
#define addPoint4 theB.Add(PointType(theO.Coord() - aRam * theXd.Coord()))
#define addPoint5 \
theB.Add(PointType(theO.Coord() - aRam * PI4 * theXd.Coord() - aRbm * PI4 * theYd.Coord()))
#define addPoint6 theB.Add(PointType(theO.Coord() - aRbm * theYd.Coord()))
#define addPoint7 \
theB.Add(PointType(theO.Coord() + aRam * PI4 * theXd.Coord() - aRbm * PI4 * theYd.Coord()))
Standard_Integer aDeb = (Standard_Integer)(aTeta1 / (M_PI / 4.)); Standard_Integer aDeb = static_cast<Standard_Integer>(aTeta1 / (M_PI / 4.));
Standard_Integer aFin = (Standard_Integer)(aTeta2 / (M_PI / 4.)); Standard_Integer aFin = static_cast<Standard_Integer>(aTeta2 / (M_PI / 4.));
aDeb++; aDeb++;
if (aDeb > aFin) if (aDeb > aFin)
return;
switch (aDeb)
{ {
case 1: { return;
addPoint1; }
if (aFin <= 1)
break; // Add polygon points from aDeb to aFin, wrapping around using modulo 8.
} for (Standard_Integer i = aDeb; i <= aFin; ++i)
Standard_FALLTHROUGH {
case 2: { addPoint(i % 8);
addPoint2;
if (aFin <= 2)
break;
}
Standard_FALLTHROUGH
case 3: {
addPoint3;
if (aFin <= 3)
break;
}
Standard_FALLTHROUGH
case 4: {
addPoint4;
if (aFin <= 4)
break;
}
Standard_FALLTHROUGH
case 5: {
addPoint5;
if (aFin <= 5)
break;
}
Standard_FALLTHROUGH
case 6: {
addPoint6;
if (aFin <= 6)
break;
}
Standard_FALLTHROUGH
case 7: {
addPoint7;
if (aFin <= 7)
break;
}
Standard_FALLTHROUGH
case 8: {
addPoint0;
if (aFin <= 8)
break;
}
Standard_FALLTHROUGH
case 9: {
addPoint1;
if (aFin <= 9)
break;
}
Standard_FALLTHROUGH
case 10: {
addPoint2;
if (aFin <= 10)
break;
}
Standard_FALLTHROUGH
case 11: {
addPoint3;
if (aFin <= 11)
break;
}
Standard_FALLTHROUGH
case 12: {
addPoint4;
if (aFin <= 12)
break;
}
Standard_FALLTHROUGH
case 13: {
addPoint5;
if (aFin <= 13)
break;
}
Standard_FALLTHROUGH
case 14: {
addPoint6;
if (aFin <= 14)
break;
}
Standard_FALLTHROUGH
case 15: {
addPoint7;
if (aFin <= 15)
break;
}
} }
} }
} // end namespace } // end namespace
static void OpenMin(const gp_Dir& V, Bnd_Box& B) static void OpenMin(const gp_Dir& V, Bnd_Box& B)
{ {
gp_Dir OX(gp_Dir::D::X); // OpenMin opens the box in the direction of decreasing parameter.
gp_Dir OY(gp_Dir::D::Y); // For a line P(t) = Origin + t*V, as t -> -Inf:
gp_Dir OZ(gp_Dir::D::Z); // - If V.X() > 0, x-coordinate -> -Inf, so open Xmin
if (V.IsParallel(OX, Precision::Angular())) // - If V.X() < 0, x-coordinate -> +Inf, so open Xmax
B.OpenXmin(); if (V.IsParallel(gp::DX(), Precision::Angular()))
else if (V.IsParallel(OY, Precision::Angular())) {
B.OpenYmin(); if (V.X() > 0.)
else if (V.IsParallel(OZ, Precision::Angular())) B.OpenXmin();
B.OpenZmin(); else
B.OpenXmax();
}
else if (V.IsParallel(gp::DY(), Precision::Angular()))
{
if (V.Y() > 0.)
B.OpenYmin();
else
B.OpenYmax();
}
else if (V.IsParallel(gp::DZ(), Precision::Angular()))
{
if (V.Z() > 0.)
B.OpenZmin();
else
B.OpenZmax();
}
else else
{ {
B.OpenXmin(); B.OpenXmin();
@@ -257,15 +182,31 @@ static void OpenMin(const gp_Dir& V, Bnd_Box& B)
static void OpenMax(const gp_Dir& V, Bnd_Box& B) static void OpenMax(const gp_Dir& V, Bnd_Box& B)
{ {
gp_Dir OX(gp_Dir::D::X); // OpenMax opens the box in the direction of increasing parameter.
gp_Dir OY(gp_Dir::D::Y); // For a line P(t) = Origin + t*V, as t -> +Inf:
gp_Dir OZ(gp_Dir::D::Z); // - If V.X() > 0, x-coordinate -> +Inf, so open Xmax
if (V.IsParallel(OX, Precision::Angular())) // - If V.X() < 0, x-coordinate -> -Inf, so open Xmin
B.OpenXmax(); if (V.IsParallel(gp::DX(), Precision::Angular()))
else if (V.IsParallel(OY, Precision::Angular())) {
B.OpenYmax(); if (V.X() > 0.)
else if (V.IsParallel(OZ, Precision::Angular())) B.OpenXmax();
B.OpenZmax(); else
B.OpenXmin();
}
else if (V.IsParallel(gp::DY(), Precision::Angular()))
{
if (V.Y() > 0.)
B.OpenYmax();
else
B.OpenYmin();
}
else if (V.IsParallel(gp::DZ(), Precision::Angular()))
{
if (V.Z() > 0.)
B.OpenZmax();
else
B.OpenZmin();
}
else else
{ {
B.OpenXmax(); B.OpenXmax();
@@ -276,20 +217,17 @@ static void OpenMax(const gp_Dir& V, Bnd_Box& B)
static void OpenMinMax(const gp_Dir& V, Bnd_Box& B) static void OpenMinMax(const gp_Dir& V, Bnd_Box& B)
{ {
gp_Dir OX(gp_Dir::D::X); if (V.IsParallel(gp::DX(), Precision::Angular()))
gp_Dir OY(gp_Dir::D::Y);
gp_Dir OZ(gp_Dir::D::Z);
if (V.IsParallel(OX, Precision::Angular()))
{ {
B.OpenXmax(); B.OpenXmax();
B.OpenXmin(); B.OpenXmin();
} }
else if (V.IsParallel(OY, Precision::Angular())) else if (V.IsParallel(gp::DY(), Precision::Angular()))
{ {
B.OpenYmax(); B.OpenYmax();
B.OpenYmin(); B.OpenYmin();
} }
else if (V.IsParallel(OZ, Precision::Angular())) else if (V.IsParallel(gp::DZ(), Precision::Angular()))
{ {
B.OpenZmax(); B.OpenZmax();
B.OpenZmin(); B.OpenZmin();
@@ -307,12 +245,24 @@ static void OpenMinMax(const gp_Dir& V, Bnd_Box& B)
static void OpenMin(const gp_Dir2d& V, Bnd_Box2d& B) static void OpenMin(const gp_Dir2d& V, Bnd_Box2d& B)
{ {
gp_Dir2d OX(gp_Dir2d::D::X); // OpenMin opens the box in the direction of decreasing parameter.
gp_Dir2d OY(gp_Dir2d::D::Y); // For a line P(t) = Origin + t*V, as t -> -Inf:
if (V.IsParallel(OX, Precision::Angular())) // - If V.X() > 0, x-coordinate -> -Inf, so open Xmin
B.OpenXmin(); // - If V.X() < 0, x-coordinate -> +Inf, so open Xmax
else if (V.IsParallel(OY, Precision::Angular())) if (V.IsParallel(gp::DX2d(), Precision::Angular()))
B.OpenYmin(); {
if (V.X() > 0.)
B.OpenXmin();
else
B.OpenXmax();
}
else if (V.IsParallel(gp::DY2d(), Precision::Angular()))
{
if (V.Y() > 0.)
B.OpenYmin();
else
B.OpenYmax();
}
else else
{ {
B.OpenXmin(); B.OpenXmin();
@@ -322,12 +272,24 @@ static void OpenMin(const gp_Dir2d& V, Bnd_Box2d& B)
static void OpenMax(const gp_Dir2d& V, Bnd_Box2d& B) static void OpenMax(const gp_Dir2d& V, Bnd_Box2d& B)
{ {
gp_Dir2d OX(gp_Dir2d::D::X); // OpenMax opens the box in the direction of increasing parameter.
gp_Dir2d OY(gp_Dir2d::D::Y); // For a line P(t) = Origin + t*V, as t -> +Inf:
if (V.IsParallel(OX, Precision::Angular())) // - If V.X() > 0, x-coordinate -> +Inf, so open Xmax
B.OpenXmax(); // - If V.X() < 0, x-coordinate -> -Inf, so open Xmin
else if (V.IsParallel(OY, Precision::Angular())) if (V.IsParallel(gp::DX2d(), Precision::Angular()))
B.OpenYmax(); {
if (V.X() > 0.)
B.OpenXmax();
else
B.OpenXmin();
}
else if (V.IsParallel(gp::DY2d(), Precision::Angular()))
{
if (V.Y() > 0.)
B.OpenYmax();
else
B.OpenYmin();
}
else else
{ {
B.OpenXmax(); B.OpenXmax();
@@ -337,14 +299,12 @@ static void OpenMax(const gp_Dir2d& V, Bnd_Box2d& B)
static void OpenMinMax(const gp_Dir2d& V, Bnd_Box2d& B) static void OpenMinMax(const gp_Dir2d& V, Bnd_Box2d& B)
{ {
gp_Dir2d OX(gp_Dir2d::D::X); if (V.IsParallel(gp::DX2d(), Precision::Angular()))
gp_Dir2d OY(gp_Dir2d::D::Y);
if (V.IsParallel(OX, Precision::Angular()))
{ {
B.OpenXmax(); B.OpenXmax();
B.OpenXmin(); B.OpenXmin();
} }
else if (V.IsParallel(OY, Precision::Angular())) else if (V.IsParallel(gp::DY2d(), Precision::Angular()))
{ {
B.OpenYmax(); B.OpenYmax();
B.OpenYmin(); B.OpenYmin();
@@ -1162,6 +1122,9 @@ void BndLib::Add(const gp_Cylinder& S,
const Standard_Real Tol, const Standard_Real Tol,
Bnd_Box& B) Bnd_Box& B)
{ {
// Cache axis direction for infinite cases.
const gp_Dir& aDir = S.Axis().Direction();
if (Precision::IsNegativeInfinite(VMin)) if (Precision::IsNegativeInfinite(VMin))
{ {
if (Precision::IsNegativeInfinite(VMax)) if (Precision::IsNegativeInfinite(VMax))
@@ -1170,19 +1133,19 @@ void BndLib::Add(const gp_Cylinder& S,
} }
else if (Precision::IsPositiveInfinite(VMax)) else if (Precision::IsPositiveInfinite(VMax))
{ {
OpenMinMax(S.Axis().Direction(), B); OpenMinMax(aDir, B);
} }
else else
{ {
ComputeCyl(S, UMin, UMax, 0., VMax, B); ComputeCyl(S, UMin, UMax, 0., VMax, B);
OpenMin(S.Axis().Direction(), B); OpenMin(aDir, B);
} }
} }
else if (Precision::IsPositiveInfinite(VMin)) else if (Precision::IsPositiveInfinite(VMin))
{ {
if (Precision::IsNegativeInfinite(VMax)) if (Precision::IsNegativeInfinite(VMax))
{ {
OpenMinMax(S.Axis().Direction(), B); OpenMinMax(aDir, B);
} }
else if (Precision::IsPositiveInfinite(VMax)) else if (Precision::IsPositiveInfinite(VMax))
{ {
@@ -1191,7 +1154,7 @@ void BndLib::Add(const gp_Cylinder& S,
else else
{ {
ComputeCyl(S, UMin, UMax, 0., VMax, B); ComputeCyl(S, UMin, UMax, 0., VMax, B);
OpenMax(S.Axis().Direction(), B); OpenMax(aDir, B);
} }
} }
else else
@@ -1199,12 +1162,12 @@ void BndLib::Add(const gp_Cylinder& S,
if (Precision::IsNegativeInfinite(VMax)) if (Precision::IsNegativeInfinite(VMax))
{ {
ComputeCyl(S, UMin, UMax, VMin, 0., B); ComputeCyl(S, UMin, UMax, VMin, 0., B);
OpenMin(S.Axis().Direction(), B); OpenMin(aDir, B);
} }
else if (Precision::IsPositiveInfinite(VMax)) else if (Precision::IsPositiveInfinite(VMax))
{ {
ComputeCyl(S, UMin, UMax, VMin, 0., B); ComputeCyl(S, UMin, UMax, VMin, 0., B);
OpenMax(S.Axis().Direction(), B); OpenMax(aDir, B);
} }
else else
{ {
@@ -1264,8 +1227,9 @@ void BndLib::Add(const gp_Cone& S,
const Standard_Real Tol, const Standard_Real Tol,
Bnd_Box& B) Bnd_Box& B)
{ {
// Cache axis direction for infinite cases.
const gp_Dir& aD = S.Axis().Direction();
Standard_Real A = S.SemiAngle();
if (Precision::IsNegativeInfinite(VMin)) if (Precision::IsNegativeInfinite(VMin))
{ {
if (Precision::IsNegativeInfinite(VMax)) if (Precision::IsNegativeInfinite(VMax))
@@ -1274,22 +1238,19 @@ void BndLib::Add(const gp_Cone& S,
} }
else if (Precision::IsPositiveInfinite(VMax)) else if (Precision::IsPositiveInfinite(VMax))
{ {
gp_Dir D(std::cos(A) * S.Axis().Direction()); OpenMinMax(aD, B);
OpenMinMax(D, B);
} }
else else
{ {
ComputeCone(S, UMin, UMax, 0., VMax, B); ComputeCone(S, UMin, UMax, 0., VMax, B);
gp_Dir D(std::cos(A) * S.Axis().Direction()); OpenMin(aD, B);
OpenMin(D, B);
} }
} }
else if (Precision::IsPositiveInfinite(VMin)) else if (Precision::IsPositiveInfinite(VMin))
{ {
if (Precision::IsNegativeInfinite(VMax)) if (Precision::IsNegativeInfinite(VMax))
{ {
gp_Dir D(std::cos(A) * S.Axis().Direction()); OpenMinMax(aD, B);
OpenMinMax(D, B);
} }
else if (Precision::IsPositiveInfinite(VMax)) else if (Precision::IsPositiveInfinite(VMax))
{ {
@@ -1298,8 +1259,7 @@ void BndLib::Add(const gp_Cone& S,
else else
{ {
ComputeCone(S, UMin, UMax, 0., VMax, B); ComputeCone(S, UMin, UMax, 0., VMax, B);
gp_Dir D(std::cos(A) * S.Axis().Direction()); OpenMax(aD, B);
OpenMax(D, B);
} }
} }
else else
@@ -1307,14 +1267,12 @@ void BndLib::Add(const gp_Cone& S,
if (Precision::IsNegativeInfinite(VMax)) if (Precision::IsNegativeInfinite(VMax))
{ {
ComputeCone(S, UMin, UMax, VMin, 0., B); ComputeCone(S, UMin, UMax, VMin, 0., B);
gp_Dir D(std::cos(A) * S.Axis().Direction()); OpenMin(aD, B);
OpenMin(D, B);
} }
else if (Precision::IsPositiveInfinite(VMax)) else if (Precision::IsPositiveInfinite(VMax))
{ {
ComputeCone(S, UMin, UMax, VMin, 0., B); ComputeCone(S, UMin, UMax, VMin, 0., B);
gp_Dir D(std::cos(A) * S.Axis().Direction()); OpenMax(aD, B);
OpenMax(D, B);
} }
else else
{ {
@@ -1570,18 +1528,18 @@ void BndLib::Add(const gp_Torus& S,
Standard_Integer Fi2; Standard_Integer Fi2;
if (VMax < VMin) if (VMax < VMin)
{ {
Fi1 = (Standard_Integer)(VMax / (M_PI / 4.)); Fi1 = static_cast<Standard_Integer>(VMax / (M_PI / 4.));
Fi2 = (Standard_Integer)(VMin / (M_PI / 4.)); Fi2 = static_cast<Standard_Integer>(VMin / (M_PI / 4.));
} }
else else
{ {
Fi1 = (Standard_Integer)(VMin / (M_PI / 4.)); Fi1 = static_cast<Standard_Integer>(VMin / (M_PI / 4.));
Fi2 = (Standard_Integer)(VMax / (M_PI / 4.)); Fi2 = static_cast<Standard_Integer>(VMax / (M_PI / 4.));
} }
Fi2++; Fi2++;
Standard_Real Ra = S.MajorRadius(); const Standard_Real Ra = S.MajorRadius();
Standard_Real Ri = S.MinorRadius(); const Standard_Real Ri = S.MinorRadius();
if (Fi2 < Fi1) if (Fi2 < Fi1)
return; return;
@@ -1593,181 +1551,61 @@ void BndLib::Add(const gp_Torus& S,
return; return;
} }
#define SC 0.71 // Cache direction vectors.
#define addP0 \ const gp_XYZ aZDir = S.Axis().Direction().XYZ();
(Compute(UMin, \ const gp_XYZ aLocXYZ = S.Location().XYZ();
UMax, \ const gp_Pnt aXd(S.XAxis().Direction().XYZ());
Ra + Ri, \ const gp_Pnt aYd(S.YAxis().Direction().XYZ());
Ra + Ri, \
gp_Pnt(S.XAxis().Direction().XYZ()), \
gp_Pnt(S.YAxis().Direction().XYZ()), \
S.Location(), \
B))
#define addP1 \
(Compute(UMin, \
UMax, \
Ra + Ri * SC, \
Ra + Ri * SC, \
gp_Pnt(S.XAxis().Direction().XYZ()), \
gp_Pnt(S.YAxis().Direction().XYZ()), \
gp_Pnt(S.Location().XYZ() + (Ri * SC) * S.Axis().Direction().XYZ()), \
B))
#define addP2 \
(Compute(UMin, \
UMax, \
Ra, \
Ra, \
gp_Pnt(S.XAxis().Direction().XYZ()), \
gp_Pnt(S.YAxis().Direction().XYZ()), \
gp_Pnt(S.Location().XYZ() + Ri * S.Axis().Direction().XYZ()), \
B))
#define addP3 \
(Compute(UMin, \
UMax, \
Ra - Ri * SC, \
Ra - Ri * SC, \
gp_Pnt(S.XAxis().Direction().XYZ()), \
gp_Pnt(S.YAxis().Direction().XYZ()), \
gp_Pnt(S.Location().XYZ() + (Ri * SC) * S.Axis().Direction().XYZ()), \
B))
#define addP4 \
(Compute(UMin, \
UMax, \
Ra - Ri, \
Ra - Ri, \
gp_Pnt(S.XAxis().Direction().XYZ()), \
gp_Pnt(S.YAxis().Direction().XYZ()), \
S.Location(), \
B))
#define addP5 \
(Compute(UMin, \
UMax, \
Ra - Ri * SC, \
Ra - Ri * SC, \
gp_Pnt(S.XAxis().Direction().XYZ()), \
gp_Pnt(S.YAxis().Direction().XYZ()), \
gp_Pnt(S.Location().XYZ() - (Ri * SC) * S.Axis().Direction().XYZ()), \
B))
#define addP6 \
(Compute(UMin, \
UMax, \
Ra, \
Ra, \
gp_Pnt(S.XAxis().Direction().XYZ()), \
gp_Pnt(S.YAxis().Direction().XYZ()), \
gp_Pnt(S.Location().XYZ() - Ri * S.Axis().Direction().XYZ()), \
B))
#define addP7 \
(Compute(UMin, \
UMax, \
Ra + Ri * SC, \
Ra + Ri * SC, \
gp_Pnt(S.XAxis().Direction().XYZ()), \
gp_Pnt(S.YAxis().Direction().XYZ()), \
gp_Pnt(S.Location().XYZ() - (Ri * SC) * S.Axis().Direction().XYZ()), \
B))
switch (Fi1) // Multipliers for torus cross-section points at 45-degree intervals.
// radiusMult[i]: multiplier for Ri in radius calculation (Ra + Ri * radiusMult[i])
// zMult[i]: multiplier for Ri in Z offset calculation (Ri * zMult[i])
// THE_COS_PI4 = cos(45 deg) = sin(45 deg) = 0.707...
constexpr Standard_Real aRadiusMult[8] =
{1., THE_COS_PI4, 0., -THE_COS_PI4, -1., -THE_COS_PI4, 0., THE_COS_PI4};
constexpr Standard_Real aZMult[8] =
{0., THE_COS_PI4, 1., THE_COS_PI4, 0., -THE_COS_PI4, -1., -THE_COS_PI4};
// Lambda to add torus cross-section point by index (0-7).
const auto addTorusPoint = [&](Standard_Integer theIdx) {
const Standard_Real aRadius = Ra + Ri * aRadiusMult[theIdx];
const gp_Pnt aCenter(aLocXYZ + (Ri * aZMult[theIdx]) * aZDir);
Compute(UMin, UMax, aRadius, aRadius, aXd, aYd, aCenter, B);
};
// Add points from Fi1 to Fi2, handling wrap-around for indices.
// Use ((i % 8) + 8) % 8 to handle negative indices correctly
// (C++ modulo can return negative values for negative dividends).
for (Standard_Integer i = Fi1; i <= Fi2; ++i)
{ {
case 0: { addTorusPoint(((i % 8) + 8) % 8);
addP0;
if (Fi2 <= 0)
break;
}
Standard_FALLTHROUGH
case 1: {
addP1;
if (Fi2 <= 1)
break;
}
Standard_FALLTHROUGH
case 2: {
addP2;
if (Fi2 <= 2)
break;
}
Standard_FALLTHROUGH
case 3: {
addP3;
if (Fi2 <= 3)
break;
}
Standard_FALLTHROUGH
case 4: {
addP4;
if (Fi2 <= 4)
break;
}
Standard_FALLTHROUGH
case 5: {
addP5;
if (Fi2 <= 5)
break;
}
Standard_FALLTHROUGH
case 6: {
addP6;
if (Fi2 <= 6)
break;
}
Standard_FALLTHROUGH
case 7: {
addP7;
if (Fi2 <= 7)
break;
}
Standard_FALLTHROUGH
case 8:
default: {
addP0;
switch (Fi2)
{
case 15:
addP7;
Standard_FALLTHROUGH
case 14:
addP6;
Standard_FALLTHROUGH
case 13:
addP5;
Standard_FALLTHROUGH
case 12:
addP4;
Standard_FALLTHROUGH
case 11:
addP3;
Standard_FALLTHROUGH
case 10:
addP2;
Standard_FALLTHROUGH
case 9:
addP1;
Standard_FALLTHROUGH
case 8:
break;
}
}
} }
B.Enlarge(Tol); B.Enlarge(Tol);
} }
void BndLib::Add(const gp_Torus& S, const Standard_Real Tol, Bnd_Box& B) void BndLib::Add(const gp_Torus& S, const Standard_Real Tol, Bnd_Box& B)
{ {
const Standard_Real aRMa = S.MajorRadius();
Standard_Real RMa = S.MajorRadius(); const Standard_Real aRmi = S.MinorRadius();
Standard_Real Rmi = S.MinorRadius(); const Standard_Real aR = aRMa + aRmi;
gp_XYZ O = S.Location().XYZ(); const gp_XYZ aO = S.Location().XYZ();
gp_XYZ Xd = S.XAxis().Direction().XYZ(); const gp_XYZ aXd = S.XAxis().Direction().XYZ();
gp_XYZ Yd = S.YAxis().Direction().XYZ(); const gp_XYZ aYd = S.YAxis().Direction().XYZ();
gp_XYZ Zd = S.Axis().Direction().XYZ(); const gp_XYZ aZd = S.Axis().Direction().XYZ();
B.Add(gp_Pnt(O - (RMa + Rmi) * Xd - (RMa + Rmi) * Yd + Rmi * Zd)); // Precompute scaled direction vectors.
B.Add(gp_Pnt(O - (RMa + Rmi) * Xd - (RMa + Rmi) * Yd - Rmi * Zd)); const gp_XYZ aRXd = aR * aXd;
B.Add(gp_Pnt(O + (RMa + Rmi) * Xd - (RMa + Rmi) * Yd + Rmi * Zd)); const gp_XYZ aRYd = aR * aYd;
B.Add(gp_Pnt(O + (RMa + Rmi) * Xd - (RMa + Rmi) * Yd - Rmi * Zd)); const gp_XYZ aRiZd = aRmi * aZd;
B.Add(gp_Pnt(O - (RMa + Rmi) * Xd + (RMa + Rmi) * Yd + Rmi * Zd)); // Add 8 corner points of torus bounding box.
B.Add(gp_Pnt(O - (RMa + Rmi) * Xd + (RMa + Rmi) * Yd - Rmi * Zd)); B.Add(gp_Pnt(aO - aRXd - aRYd + aRiZd));
B.Add(gp_Pnt(O + (RMa + Rmi) * Xd + (RMa + Rmi) * Yd + Rmi * Zd)); B.Add(gp_Pnt(aO - aRXd - aRYd - aRiZd));
B.Add(gp_Pnt(O + (RMa + Rmi) * Xd + (RMa + Rmi) * Yd - Rmi * Zd)); B.Add(gp_Pnt(aO + aRXd - aRYd + aRiZd));
B.Add(gp_Pnt(aO + aRXd - aRYd - aRiZd));
B.Add(gp_Pnt(aO - aRXd + aRYd + aRiZd));
B.Add(gp_Pnt(aO - aRXd + aRYd - aRiZd));
B.Add(gp_Pnt(aO + aRXd + aRYd + aRiZd));
B.Add(gp_Pnt(aO + aRXd + aRYd - aRiZd));
B.Enlarge(Tol); B.Enlarge(Tol);
} }
@@ -1804,7 +1642,8 @@ Standard_Integer ComputeBox(const gp_Hypr& aHypr,
aRmaj = aHypr.MajorRadius(); aRmaj = aHypr.MajorRadius();
aRmin = aHypr.MinorRadius(); aRmin = aHypr.MinorRadius();
// //
aT3 = 0; // Find extrema for each coordinate (X, Y, Z) independently.
// Each coordinate can have its extremum at a different parameter value.
for (i = 1; i <= 3; ++i) for (i = 1; i <= 3; ++i)
{ {
aA = aRmin * aYDir.Coord(i); aA = aRmin * aYDir.Coord(i);
@@ -1813,8 +1652,8 @@ Standard_Integer ComputeBox(const gp_Hypr& aHypr,
aABP = aA + aB; aABP = aA + aB;
aBAM = aB - aA; aBAM = aB - aA;
// //
aABP = fabs(aABP); aABP = std::abs(aABP);
aBAM = fabs(aBAM); aBAM = std::abs(aBAM);
// //
if (aABP < aEps || aBAM < aEps) if (aABP < aEps || aBAM < aEps)
{ {
@@ -1822,23 +1661,17 @@ Standard_Integer ComputeBox(const gp_Hypr& aHypr,
} }
// //
aCf = aBAM / aABP; aCf = aBAM / aABP;
aT3 = log(sqrt(aCf)); aT3 = 0.5 * std::log(aCf);
// //
if (aT3 < aT1 || aT3 > aT2) if (aT3 < aT1 || aT3 > aT2)
{ {
continue; continue;
} }
// Add extremum point for this coordinate.
aP3 = ElCLib::Value(aT3, aHypr);
aBox.Add(aP3);
iErr = 0; iErr = 0;
break;
} }
// //
if (iErr)
{
return iErr;
}
//
aP3 = ElCLib::Value(aT3, aHypr);
aBox.Add(aP3);
//
return iErr; return iErr;
} }

View File

@@ -143,7 +143,7 @@ public:
Standard_Integer NbVariables() const { return 1; } Standard_Integer NbVariables() const { return 1; }
private: private:
Curv2dMaxMinCoordMVar& operator=(const Curv2dMaxMinCoordMVar& theOther); Curv2dMaxMinCoordMVar& operator=(const Curv2dMaxMinCoordMVar&) = delete;
Standard_Boolean CheckInputData(Standard_Real theParam) Standard_Boolean CheckInputData(Standard_Real theParam)
{ {
@@ -190,7 +190,7 @@ public:
} }
private: private:
Curv2dMaxMinCoord& operator=(const Curv2dMaxMinCoord& theOther); Curv2dMaxMinCoord& operator=(const Curv2dMaxMinCoord&) = delete;
Standard_Boolean CheckInputData(Standard_Real theParam) Standard_Boolean CheckInputData(Standard_Real theParam)
{ {
@@ -418,7 +418,8 @@ void BndLib_Box2dCurve::PerformBezier()
aTb[1] = aT2; aTb[1] = aT2;
} }
// //
if (!(aT1 == aTb[0] && aT2 == aTb[1])) constexpr Standard_Real anEps = Precision::PConfusion();
if (std::abs(aT1 - aTb[0]) > anEps || std::abs(aT2 - aTb[1]) > anEps)
{ {
aG = aCBz->Copy(); aG = aCBz->Copy();
// //
@@ -477,7 +478,7 @@ void BndLib_Box2dCurve::PerformBSpline()
// //
constexpr Standard_Real eps = Precision::PConfusion(); constexpr Standard_Real eps = Precision::PConfusion();
if (fabs(aT1 - aTb[0]) > eps || fabs(aT2 - aTb[1]) > eps) if (std::abs(aT1 - aTb[0]) > eps || std::abs(aT2 - aTb[1]) > eps)
{ {
aG = aCBS->Copy(); aG = aCBS->Copy();
// //
@@ -738,7 +739,7 @@ void BndLib_Box2dCurve::D0(const Standard_Real aU, gp_Pnt2d& aP2D)
// //
aA = aV1.Y(); aA = aV1.Y();
aB = -aV1.X(); aB = -aV1.X();
aR = sqrt(aA * aA + aB * aB); aR = std::sqrt(aA * aA + aB * aB);
if (aR <= aRes) if (aR <= aRes)
{ {
myErrorStatus = 13; myErrorStatus = 13;
@@ -948,18 +949,18 @@ void BndLib_Box2dCurve::Compute(const Handle(Geom2d_Conic)& aConic2D,
} }
// //
// aType==GeomAbs_Circle || aType==GeomAbs_Ellipse // aType==GeomAbs_Circle || aType==GeomAbs_Ellipse
aEps = 1.e-14; aEps = Precision::Angular();
aTwoPI = 2. * M_PI; aTwoPI = 2. * M_PI;
dT = aT2 - aT1; dT = aT2 - aT1;
// //
Standard_Real aT1z = AdjustToPeriod(aT1, aTwoPI); Standard_Real aT1z = AdjustToPeriod(aT1, aTwoPI);
if (fabs(aT1z) < aEps) if (std::abs(aT1z) < aEps)
{ {
aT1z = 0.; aT1z = 0.;
} }
// //
Standard_Real aT2z = aT1z + dT; Standard_Real aT2z = aT1z + dT;
if (fabs(aT2z - aTwoPI) < aEps) if (std::abs(aT2z - aTwoPI) < aEps)
{ {
aT2z = aTwoPI; aT2z = aTwoPI;
} }
@@ -1033,12 +1034,12 @@ Standard_Integer BndLib_Box2dCurve::Compute(const Handle(Geom2d_Conic)& aConic2D
aLy = (!i) ? 1. : 0.; aLy = (!i) ? 1. : 0.;
aBx = aLx * aA21 - aLy * aA11; aBx = aLx * aA21 - aLy * aA11;
aBy = aLx * aA22 - aLy * aA12; aBy = aLx * aA22 - aLy * aA12;
aB = sqrt(aBx * aBx + aBy * aBy); aB = std::sqrt(aBx * aBx + aBy * aBy);
// //
aCosFi = aBx / aB; aCosFi = aBx / aB;
aSinFi = aBy / aB; aSinFi = aBy / aB;
// //
aFi = acos(aCosFi); aFi = std::acos(aCosFi);
if (aSinFi < 0.) if (aSinFi < 0.)
{ {
aFi = aTwoPI - aFi; aFi = aTwoPI - aFi;
@@ -1060,7 +1061,7 @@ Standard_Integer BndLib_Box2dCurve::Compute(const Handle(Geom2d_Conic)& aConic2D
Standard_Real aA1, aA2; Standard_Real aA1, aA2;
Handle(Geom2d_Parabola) aPR2D; Handle(Geom2d_Parabola) aPR2D;
// //
aEps = 1.e-12; aEps = Precision::Angular();
// //
aPR2D = Handle(Geom2d_Parabola)::DownCast(aConic2D); aPR2D = Handle(Geom2d_Parabola)::DownCast(aConic2D);
aFc = aPR2D->Focal(); aFc = aPR2D->Focal();
@@ -1072,7 +1073,7 @@ Standard_Integer BndLib_Box2dCurve::Compute(const Handle(Geom2d_Conic)& aConic2D
aLy = (!i) ? 1. : 0.; aLy = (!i) ? 1. : 0.;
// //
aA2 = aLx * aSinBt - aLy * aCosBt; aA2 = aLx * aSinBt - aLy * aCosBt;
if (fabs(aA2) < aEps) if (std::abs(aA2) < aEps)
{ {
continue; continue;
} }
@@ -1092,7 +1093,7 @@ Standard_Integer BndLib_Box2dCurve::Compute(const Handle(Geom2d_Conic)& aConic2D
Standard_Real aEps, aB1, aB2, aB12, aB22, aZ, aD; Standard_Real aEps, aB1, aB2, aB12, aB22, aZ, aD;
Handle(Geom2d_Hyperbola) aHP2D; Handle(Geom2d_Hyperbola) aHP2D;
// //
aEps = 1.e-12; aEps = Precision::Angular();
// //
aHP2D = Handle(Geom2d_Hyperbola)::DownCast(aConic2D); aHP2D = Handle(Geom2d_Hyperbola)::DownCast(aConic2D);
aR1 = aHP2D->MajorRadius(); aR1 = aHP2D->MajorRadius();
@@ -1107,12 +1108,12 @@ Standard_Integer BndLib_Box2dCurve::Compute(const Handle(Geom2d_Conic)& aConic2D
aB1 = aR1 * (aLx * aSinBt - aLy * aCosBt); aB1 = aR1 * (aLx * aSinBt - aLy * aCosBt);
aB2 = aR2 * (aLx * aSinGm - aLy * aCosGm); aB2 = aR2 * (aLx * aSinGm - aLy * aCosGm);
// //
if (fabs(aB1) < aEps) if (std::abs(aB1) < aEps)
{ {
continue; continue;
} }
// //
if (fabs(aB2) < aEps) if (std::abs(aB2) < aEps)
{ {
pT[j] = 0.; pT[j] = 0.;
++j; ++j;
@@ -1126,14 +1127,14 @@ Standard_Integer BndLib_Box2dCurve::Compute(const Handle(Geom2d_Conic)& aConic2D
continue; continue;
} }
// //
aD = sqrt(aB12 - aB22); aD = std::sqrt(aB12 - aB22);
//------------- //-------------
for (k = -1; k < 2; k += 2) for (k = -1; k < 2; k += 2)
{ {
aZ = (aB1 + k * aD) / aB2; aZ = (aB1 + k * aD) / aB2;
if (fabs(aZ) < 1.) if (std::abs(aZ) < 1.)
{ {
pT[j] = -log((1. + aZ) / (1. - aZ)); pT[j] = -std::log((1. + aZ) / (1. - aZ));
++j; ++j;
} }
} }
@@ -1155,12 +1156,12 @@ Standard_Real BndLib_Box2dCurve::AdjustToPeriod(const Standard_Real aT, const St
aTRet = aT; aTRet = aT;
if (aT < 0.) if (aT < 0.)
{ {
k = 1 + (Standard_Integer)(-aT / aPeriod); k = 1 + static_cast<Standard_Integer>(-aT / aPeriod);
aTRet = aT + k * aPeriod; aTRet = aT + k * aPeriod;
} }
else if (aT > aPeriod) else if (aT > aPeriod)
{ {
k = (Standard_Integer)(aT / aPeriod); k = static_cast<Standard_Integer>(aT / aPeriod);
aTRet = aT - k * aPeriod; aTRet = aT - k * aPeriod;
} }
if (aTRet == aPeriod) if (aTRet == aPeriod)
@@ -1197,8 +1198,7 @@ void BndLib_Add2dCurve::Add(const Adaptor2d_Curve2d& aC,
const Standard_Real aTol, const Standard_Real aTol,
Bnd_Box2d& aBox2D) Bnd_Box2d& aBox2D)
{ {
Adaptor2d_Curve2d* pC = (Adaptor2d_Curve2d*)&aC; const Geom2dAdaptor_Curve* pA = dynamic_cast<const Geom2dAdaptor_Curve*>(&aC);
Geom2dAdaptor_Curve* pA = dynamic_cast<Geom2dAdaptor_Curve*>(pC);
if (!pA) if (!pA)
{ {
Standard_Real U, DU; Standard_Real U, DU;

View File

@@ -48,53 +48,46 @@ static void reduceSplineBox(const Adaptor3d_Curve& theCurve,
const Bnd_Box& theOrigBox, const Bnd_Box& theOrigBox,
Bnd_Box& theReducedBox) Bnd_Box& theReducedBox)
{ {
// Guaranteed bounding box based on poles of bspline. // Guaranteed bounding box based on poles of bspline/bezier.
Bnd_Box aPolesBox; Bnd_Box aPolesBox;
Standard_Real aPolesXMin, aPolesYMin, aPolesZMin, aPolesXMax, aPolesYMax, aPolesZMax;
if (theCurve.GetType() == GeomAbs_BSplineCurve) if (theCurve.GetType() == GeomAbs_BSplineCurve)
{ {
Handle(Geom_BSplineCurve) aC = theCurve.BSpline(); Handle(Geom_BSplineCurve) aC = theCurve.BSpline();
const TColgp_Array1OfPnt& aPoles = aC->Poles(); const TColgp_Array1OfPnt& aPoles = aC->Poles();
for (Standard_Integer anIdx = aPoles.Lower(); anIdx <= aPoles.Upper(); ++anIdx) for (Standard_Integer anIdx = aPoles.Lower(); anIdx <= aPoles.Upper(); ++anIdx)
{ {
aPolesBox.Add(aPoles.Value(anIdx)); aPolesBox.Add(aPoles.Value(anIdx));
} }
} }
if (theCurve.GetType() == GeomAbs_BezierCurve) else if (theCurve.GetType() == GeomAbs_BezierCurve)
{ {
Handle(Geom_BezierCurve) aC = theCurve.Bezier(); Handle(Geom_BezierCurve) aC = theCurve.Bezier();
const TColgp_Array1OfPnt& aPoles = aC->Poles(); const TColgp_Array1OfPnt& aPoles = aC->Poles();
for (Standard_Integer anIdx = aPoles.Lower(); anIdx <= aPoles.Upper(); ++anIdx) for (Standard_Integer anIdx = aPoles.Lower(); anIdx <= aPoles.Upper(); ++anIdx)
{ {
aPolesBox.Add(aPoles.Value(anIdx)); aPolesBox.Add(aPoles.Value(anIdx));
} }
} }
aPolesBox.Get(aPolesXMin, aPolesYMin, aPolesZMin, aPolesXMax, aPolesYMax, aPolesZMax); // Get original box limits using C++17 structured bindings.
const auto [aXMin, aXMax, aYMin, aYMax, aZMin, aZMax] = theOrigBox.Get();
Standard_Real x, y, z, X, Y, Z; // If poles box is valid, intersect with original box.
theOrigBox.Get(x, y, z, X, Y, Z); if (!aPolesBox.IsVoid())
{
// Left bound. const auto [aPXMin, aPXMax, aPYMin, aPYMax, aPZMin, aPZMax] = aPolesBox.Get();
if (aPolesXMin > x) theReducedBox.Update(std::max(aXMin, aPXMin),
x = aPolesXMin; std::max(aYMin, aPYMin),
if (aPolesYMin > y) std::max(aZMin, aPZMin),
y = aPolesYMin; std::min(aXMax, aPXMax),
if (aPolesZMin > z) std::min(aYMax, aPYMax),
z = aPolesZMin; std::min(aZMax, aPZMax));
}
// Right bound. else
if (aPolesXMax < X) {
X = aPolesXMax; theReducedBox.Update(aXMin, aYMin, aZMin, aXMax, aYMax, aZMax);
if (aPolesYMax < Y) }
Y = aPolesYMax;
if (aPolesZMax < Z)
Z = aPolesZMax;
theReducedBox.Update(x, y, z, X, Y, Z);
} }
//================================================================================================= //=================================================================================================
@@ -152,8 +145,8 @@ void BndLib_Add3dCurve::Add(const Adaptor3d_Curve& C,
const Standard_Real Tol, const Standard_Real Tol,
Bnd_Box& B) Bnd_Box& B)
{ {
static Standard_Real weakness = 1.5; // OCC566(apo) constexpr Standard_Real weakness = 1.5; // OCC566(apo)
Standard_Real tol = 0.0; Standard_Real tol = 0.0;
switch (C.GetType()) switch (C.GetType())
{ {
@@ -270,9 +263,9 @@ void BndLib_Add3dCurve::Add(const Adaptor3d_Curve& C,
break; break;
} }
default: { default: {
Bnd_Box B1; Bnd_Box B1;
static Standard_Integer N = 33; constexpr Standard_Integer N = 33;
tol = FillBox(B1, C, U1, U2, N); tol = FillBox(B1, C, U1, U2, N);
B1.Enlarge(weakness * tol); B1.Enlarge(weakness * tol);
Standard_Real x, y, z, X, Y, Z; Standard_Real x, y, z, X, Y, Z;
B1.Get(x, y, z, X, Y, Z); B1.Get(x, y, z, X, Y, Z);
@@ -464,7 +457,7 @@ public:
Standard_Integer NbVariables() const { return 1; } Standard_Integer NbVariables() const { return 1; }
private: private:
CurvMaxMinCoordMVar& operator=(const CurvMaxMinCoordMVar& theOther); CurvMaxMinCoordMVar& operator=(const CurvMaxMinCoordMVar&) = delete;
Standard_Boolean CheckInputData(Standard_Real theParam) Standard_Boolean CheckInputData(Standard_Real theParam)
{ {
@@ -511,7 +504,7 @@ public:
} }
private: private:
CurvMaxMinCoord& operator=(const CurvMaxMinCoord& theOther); CurvMaxMinCoord& operator=(const CurvMaxMinCoord&) = delete;
Standard_Boolean CheckInputData(Standard_Real theParam) Standard_Boolean CheckInputData(Standard_Real theParam)
{ {

View File

@@ -447,10 +447,10 @@ void BndLib_AddSurface::Add(const Adaptor3d_Surface& S,
gp_Pnt P; gp_Pnt P;
for (Standard_Integer i = 1; i <= Nu; i++) for (Standard_Integer i = 1; i <= Nu; i++)
{ {
Standard_Real U = UMin + ((UMax - UMin) * (i - 1) / (Nu - 1)); const Standard_Real U = UMin + ((UMax - UMin) * (i - 1) / (Nu - 1));
for (Standard_Integer j = 1; j <= Nv; j++) for (Standard_Integer j = 1; j <= Nv; j++)
{ {
Standard_Real V = VMin + ((VMax - VMin) * (j - 1) / (Nv - 1)); const Standard_Real V = VMin + ((VMax - VMin) * (j - 1) / (Nv - 1));
S.D0(U, V, P); S.D0(U, V, P);
B.Add(P); B.Add(P);
} }
@@ -769,7 +769,7 @@ public:
Standard_Integer NbVariables() const { return 2; } Standard_Integer NbVariables() const { return 2; }
private: private:
SurfMaxMinCoord& operator=(const SurfMaxMinCoord& theOther); SurfMaxMinCoord& operator=(const SurfMaxMinCoord&) = delete;
Standard_Boolean CheckInputData(const math_Vector& theParams) Standard_Boolean CheckInputData(const math_Vector& theParams)
{ {

File diff suppressed because it is too large Load Diff

View File

@@ -2,6 +2,7 @@
set(OCCT_TKGeomBase_GTests_FILES_LOCATION "${CMAKE_CURRENT_LIST_DIR}") set(OCCT_TKGeomBase_GTests_FILES_LOCATION "${CMAKE_CURRENT_LIST_DIR}")
set(OCCT_TKGeomBase_GTests_FILES set(OCCT_TKGeomBase_GTests_FILES
BndLib_Test.cxx
Extrema_ExtPC_Test.cxx Extrema_ExtPC_Test.cxx
IntAna_IntQuadQuad_Test.cxx IntAna_IntQuadQuad_Test.cxx
) )

View File

@@ -18,7 +18,7 @@ checkshape b_1
set xmin0 -62.992562093513783 set xmin0 -62.992562093513783
set ymin0 -51.45371311501097 set ymin0 -51.45371311501097
set zmin0 -33.223762093513777 set zmin0 -33.223762093513777
set xmax0 43.263112093513783 set xmax0 43.257600511674156
set ymax0 33.211062093513789 set ymax0 33.211062093513789
set zmax0 37.744962093513792 set zmax0 37.744962093513792
set bb [bounding b_1] set bb [bounding b_1]