From 1cb09f481b94f537db9e93ed18ffaa37b8d58251 Mon Sep 17 00:00:00 2001 From: Pasukhin Dmitry Date: Fri, 7 Aug 2026 19:05:46 +0100 Subject: [PATCH] Modeling - Optimize exact face classification fallback (#1400) - Add a thread-safe lazy grid to CSLib_Class2d for constant-time classification away from polygon boundaries while preserving the exact tolerance path for boundary cells. - Normalize coordinates and explicit tolerances robustly, handle invalid and extreme input values, and preserve cache state correctly across copy and move operations. - Cache face wires, edge occurrences, pcurves, parameter ranges, and optional 2D bounding boxes in BRepClass_FaceExplorer to avoid repeated topology and geometry traversal during exact classification. - Introduce explicit cached-geometry and bounding-box states in BRepClass_Edge, invalidate derived data on topology changes, and isolate bounding-box construction failures per edge. - Reuse a lazily constructed exact face explorer in BRepTopAdaptor_FClass2d, cache periodic surface properties, and serialize access to mutable exact-classification traversal state. - Replace node-based classifier and point sequences with contiguous linear vectors, remove intermediate polygon copies, and make tighter wire discretization transactional. - Consolidate Perform() and TestOnRestriction() through a shared classification path while preserving their boundary and periodic recadrement semantics. --- .../TKMath/CSLib/CSLib_Class2d.cxx | 497 ++++++++++++-- .../TKMath/CSLib/CSLib_Class2d.hxx | 96 +-- .../TKMath/GTests/CSLib_Test.cxx | 314 ++++++++- .../TopOpeBRepBuild_GridSS.cxx | 4 +- .../TopOpeBRepBuild_WireEdgeClassifier.cxx | 4 +- .../TKTopAlgo/BRepClass/BRepClass_Edge.cxx | 93 ++- .../TKTopAlgo/BRepClass/BRepClass_Edge.hxx | 89 ++- .../BRepClass/BRepClass_FaceExplorer.cxx | 300 ++++++++- .../BRepClass/BRepClass_FaceExplorer.hxx | 55 +- .../BRepClass/BRepClass_Intersector.cxx | 43 +- .../BRepTopAdaptor_FClass2d.cxx | 614 ++++++++---------- .../BRepTopAdaptor_FClass2d.hxx | 72 +- .../GTests/BRepClass_FaceExplorer_Test.cxx | 386 +++++++++++ .../GTests/BRepTopAdaptor_FClass2d_Test.cxx | 202 ++++++ .../TKTopAlgo/GTests/FILES.cmake | 2 + 15 files changed, 2164 insertions(+), 607 deletions(-) create mode 100644 src/ModelingAlgorithms/TKTopAlgo/GTests/BRepClass_FaceExplorer_Test.cxx create mode 100644 src/ModelingAlgorithms/TKTopAlgo/GTests/BRepTopAdaptor_FClass2d_Test.cxx diff --git a/src/FoundationClasses/TKMath/CSLib/CSLib_Class2d.cxx b/src/FoundationClasses/TKMath/CSLib/CSLib_Class2d.cxx index 3f8c6097d2..2a066022ea 100644 --- a/src/FoundationClasses/TKMath/CSLib/CSLib_Class2d.cxx +++ b/src/FoundationClasses/TKMath/CSLib/CSLib_Class2d.cxx @@ -16,26 +16,91 @@ #include #include +#include #include +#include +#include #include +#include +#include namespace { +constexpr double THE_MIN_NORMALIZATION_RANGE = 1.0e-10; + +//! Returns true when direct arithmetic cannot reach OCCT's geometric infinity range. +inline bool isSafeForDirectArithmetic(const double theValue) +{ + return !std::isnan(theValue) && !Precision::IsInfinite(theValue); +} + +//! Returns true for representable values, including OCCT's finite infinity sentinels. +inline bool isRepresentableValue(const double theValue) +{ + constexpr double THE_MAX_VALUE = std::numeric_limits::max(); + return theValue >= -THE_MAX_VALUE && theValue <= THE_MAX_VALUE; +} + //! Transforms a coordinate from original space to normalized [0,1] space. //! @param[in] theU Original coordinate value //! @param[in] theUMin Minimum bound of original range -//! @param[in] theURange Range of original domain (theUMax - theUMin) +//! @param[in] theUMax Maximum bound of original range //! @return Normalized coordinate in [0,1], or original value if range is too small -inline double transformToNormalized(const double theU, const double theUMin, const double theURange) +inline double transformToNormalized(const double theU, const double theUMin, const double theUMax) { - constexpr double THE_MIN_RANGE = 1e-10; - if (theURange > THE_MIN_RANGE) + const double aRange = theUMax - theUMin; + if (aRange > THE_MIN_NORMALIZATION_RANGE) { - return (theU - theUMin) / theURange; + const double aDifference = theU - theUMin; + if (isSafeForDirectArithmetic(aDifference) && isSafeForDirectArithmetic(aRange)) + { + return aDifference / aRange; + } + const double aScale = std::max(std::abs(theUMin), std::abs(theUMax)); + const double aScaledMin = theUMin / aScale; + return (theU / aScale - aScaledMin) / (theUMax / aScale - aScaledMin); } return theU; } + +//! Converts an external tolerance to a finite, non-negative value. +inline double sanitizeTolerance(const double theTolerance) +{ + if (std::isnan(theTolerance) || theTolerance <= 0.0) + { + return 0.0; + } + return isRepresentableValue(theTolerance) ? theTolerance : std::numeric_limits::max(); +} + +//! Normalizes a distance without overflowing when the finite bounds span more +//! than the representable double range. +inline double normalizeTolerance(const double theTolerance, + const double theMin, + const double theMax) +{ + const double aTolerance = sanitizeTolerance(theTolerance); + const double aRange = theMax - theMin; + if (!(aRange > THE_MIN_NORMALIZATION_RANGE)) + { + return aTolerance; + } + if (isSafeForDirectArithmetic(aRange)) + { + return sanitizeTolerance(aTolerance / aRange); + } + + const double aScale = std::max(std::abs(theMin), std::abs(theMax)); + const double aScaledRange = theMax / aScale - theMin / aScale; + return sanitizeTolerance((aTolerance / aScale) / aScaledRange); +} + +//! Grid cache for O(1) classification of points far from polygon edges. +//! Construction is delayed until repeated queries amortize its cost. +static constexpr int THE_GRID_SIZE = 32; +static constexpr int THE_GRID_MIN_POINTS = 24; +static constexpr size_t THE_GRID_BUILD_QUERY_COUNT = 64; } // namespace //================================================================================================= @@ -55,46 +120,36 @@ void CSLib_Class2d::init(const TCol_Containers2d& thePnts2d, myVMax = theVMax; // Validate input parameters. - if (theUMax <= theUMin || theVMax <= theVMin || thePnts2d.Length() < 3) + if (!isRepresentableValue(theUMin) || !isRepresentableValue(theVMin) + || !isRepresentableValue(theUMax) || !isRepresentableValue(theVMax) || theUMax <= theUMin + || theVMax <= theVMin || thePnts2d.Length() < 3) { myPointsCount = 0; return; } - myPointsCount = thePnts2d.Length(); - myTolU = theTolU; - myTolV = theTolV; + myPointsCount = thePnts2d.Length(); + myOriginalTolU = sanitizeTolerance(theTolU); + myOriginalTolV = sanitizeTolerance(theTolV); + myTolU = normalizeTolerance(theTolU, theUMin, theUMax); + myTolV = normalizeTolerance(theTolV, theVMin, theVMax); // Allocate arrays with one extra element for closing the polygon. myPnts2dX.Resize(0, myPointsCount, false); myPnts2dY.Resize(0, myPointsCount, false); - const double aDu = theUMax - theUMin; - const double aDv = theVMax - theVMin; - // Transform points to normalized coordinates. const int aLower = thePnts2d.Lower(); for (int i = 0; i < myPointsCount; ++i) { const gp_Pnt2d& aP2D = thePnts2d(i + aLower); - myPnts2dX.ChangeValue(i) = transformToNormalized(aP2D.X(), theUMin, aDu); - myPnts2dY.ChangeValue(i) = transformToNormalized(aP2D.Y(), theVMin, aDv); + myPnts2dX.ChangeValue(i) = transformToNormalized(aP2D.X(), theUMin, theUMax); + myPnts2dY.ChangeValue(i) = transformToNormalized(aP2D.Y(), theVMin, theVMax); } // Close the polygon by copying first point to last position. myPnts2dX.ChangeLast() = myPnts2dX.First(); myPnts2dY.ChangeLast() = myPnts2dY.First(); - - // Normalize tolerances. - constexpr double THE_MIN_RANGE = 1e-10; - if (aDu > THE_MIN_RANGE) - { - myTolU /= aDu; - } - if (aDv > THE_MIN_RANGE) - { - myTolV /= aDv; - } } //================================================================================================= @@ -138,6 +193,285 @@ CSLib_Class2d::CSLib_Class2d(const NCollection_DynamicArray& thePnts2d //================================================================================================= +CSLib_Class2d::CSLib_Class2d(const CSLib_Class2d& theOther) + : myPnts2dX(theOther.myPnts2dX), + myPnts2dY(theOther.myPnts2dY), + myTolU(theOther.myTolU), + myTolV(theOther.myTolV), + myOriginalTolU(theOther.myOriginalTolU), + myOriginalTolV(theOther.myOriginalTolV), + myPointsCount(theOther.myPointsCount), + myUMin(theOther.myUMin), + myVMin(theOther.myVMin), + myUMax(theOther.myUMax), + myVMax(theOther.myVMax), + myQueryCount(theOther.myGridState.load(std::memory_order_acquire) == GridState::Building + ? 0 + : theOther.myQueryCount.load(std::memory_order_relaxed)) +{ + const GridState aState = theOther.myGridState.load(std::memory_order_acquire); + if (aState == GridState::Ready) + { + myGrid = theOther.myGrid; + myGridState.store(GridState::Ready, std::memory_order_relaxed); + } + else if (aState == GridState::Disabled) + { + myGridState.store(GridState::Disabled, std::memory_order_relaxed); + } +} + +//================================================================================================= + +CSLib_Class2d& CSLib_Class2d::operator=(const CSLib_Class2d& theOther) +{ + if (this == &theOther) + { + return *this; + } + CSLib_Class2d aCopy(theOther); + *this = std::move(aCopy); + return *this; +} + +//================================================================================================= + +CSLib_Class2d::CSLib_Class2d(CSLib_Class2d&& theOther) noexcept + : myPnts2dX(std::move(theOther.myPnts2dX)), + myPnts2dY(std::move(theOther.myPnts2dY)), + myTolU(theOther.myTolU), + myTolV(theOther.myTolV), + myOriginalTolU(theOther.myOriginalTolU), + myOriginalTolV(theOther.myOriginalTolV), + myPointsCount(theOther.myPointsCount), + myUMin(theOther.myUMin), + myVMin(theOther.myVMin), + myUMax(theOther.myUMax), + myVMax(theOther.myVMax), + myQueryCount(theOther.myGridState.load(std::memory_order_acquire) == GridState::Building + ? 0 + : theOther.myQueryCount.load(std::memory_order_relaxed)) +{ + const GridState aState = theOther.myGridState.load(std::memory_order_acquire); + if (aState == GridState::Ready) + { + myGrid = std::move(theOther.myGrid); + myGridState.store(GridState::Ready, std::memory_order_relaxed); + } + else if (aState == GridState::Disabled) + { + myGridState.store(GridState::Disabled, std::memory_order_relaxed); + } + theOther.myPointsCount = 0; + theOther.myGridState.store(GridState::Disabled, std::memory_order_release); +} + +//================================================================================================= + +CSLib_Class2d& CSLib_Class2d::operator=(CSLib_Class2d&& theOther) noexcept +{ + if (this == &theOther) + { + return *this; + } + + const GridState aState = theOther.myGridState.load(std::memory_order_acquire); + myPnts2dX = std::move(theOther.myPnts2dX); + myPnts2dY = std::move(theOther.myPnts2dY); + myTolU = theOther.myTolU; + myTolV = theOther.myTolV; + myOriginalTolU = theOther.myOriginalTolU; + myOriginalTolV = theOther.myOriginalTolV; + myPointsCount = theOther.myPointsCount; + myUMin = theOther.myUMin; + myVMin = theOther.myVMin; + myUMax = theOther.myUMax; + myVMax = theOther.myVMax; + myQueryCount.store( + aState == GridState::Building ? 0 : theOther.myQueryCount.load(std::memory_order_relaxed), + std::memory_order_relaxed); + myGrid = NCollection_Array1(); + if (aState == GridState::Ready) + { + myGrid = std::move(theOther.myGrid); + myGridState.store(GridState::Ready, std::memory_order_relaxed); + } + else + { + myGridState.store(aState == GridState::Disabled ? GridState::Disabled : GridState::NotBuilt, + std::memory_order_relaxed); + } + theOther.myPointsCount = 0; + theOther.myGridState.store(GridState::Disabled, std::memory_order_release); + return *this; +} + +//================================================================================================= + +void CSLib_Class2d::buildGridCache() const +{ + if (myPointsCount < THE_GRID_MIN_POINTS) + { + return; + } + + GridState anExpectedState = GridState::NotBuilt; + if (!myGridState.compare_exchange_strong(anExpectedState, + GridState::Building, + std::memory_order_acq_rel, + std::memory_order_acquire)) + { + return; + } + + try + { + const int aTotalCells = THE_GRID_SIZE * THE_GRID_SIZE; + myGrid.Resize(0, aTotalCells - 1, false); + myGrid.Init(GridCell_Unvisited); + + const double aCellSize = 1.0 / THE_GRID_SIZE; + const double anEps = 8.0 * std::numeric_limits::epsilon(); + const double* pX = &myPnts2dX.First(); + const double* pY = &myPnts2dY.First(); + NCollection_LinearVector aCellQueue(static_cast(aTotalCells)); + + // Non-finite polygon data and domain-sized tolerances are valid for the + // exact path but cannot produce a useful, safely indexed cache. + if (myTolU >= 1.0 || myTolV >= 1.0) + { + myGrid = NCollection_Array1(); + myGridState.store(GridState::Disabled, std::memory_order_release); + return; + } + for (int aPointIdx = 0; aPointIdx < myPointsCount; ++aPointIdx) + { + if (!isSafeForDirectArithmetic(pX[aPointIdx]) || !isSafeForDirectArithmetic(pY[aPointIdx])) + { + myGrid = NCollection_Array1(); + myGridState.store(GridState::Disabled, std::memory_order_release); + return; + } + } + + // Rasterize tolerance-expanded edge boxes into the fixed grid. This is + // equivalent to testing every cell box against every edge box, but avoids + // an O(grid cells * polygon edges) scan. + for (int anEdgeIdx = 0; anEdgeIdx < myPointsCount; ++anEdgeIdx) + { + const double aMinX = std::min(pX[anEdgeIdx], pX[anEdgeIdx + 1]) - myTolU - anEps; + const double aMaxX = std::max(pX[anEdgeIdx], pX[anEdgeIdx + 1]) + myTolU + anEps; + const double aMinY = std::min(pY[anEdgeIdx], pY[anEdgeIdx + 1]) - myTolV - anEps; + const double aMaxY = std::max(pY[anEdgeIdx], pY[anEdgeIdx + 1]) + myTolV + anEps; + if (aMaxX < 0.0 || aMinX > 1.0 || aMaxY < 0.0 || aMinY > 1.0) + { + continue; + } + + const double aClippedMinX = std::clamp(aMinX, 0.0, 1.0); + const double aClippedMaxX = std::clamp(aMaxX, 0.0, 1.0); + const double aClippedMinY = std::clamp(aMinY, 0.0, 1.0); + const double aClippedMaxY = std::clamp(aMaxY, 0.0, 1.0); + const int aMinCellX = + std::max(static_cast(std::ceil(aClippedMinX * THE_GRID_SIZE)) - 1, 0); + const int aMaxCellX = + std::min(static_cast(std::floor(aClippedMaxX * THE_GRID_SIZE)), THE_GRID_SIZE - 1); + const int aMinCellY = + std::max(static_cast(std::ceil(aClippedMinY * THE_GRID_SIZE)) - 1, 0); + const int aMaxCellY = + std::min(static_cast(std::floor(aClippedMaxY * THE_GRID_SIZE)), THE_GRID_SIZE - 1); + + for (int aCellY = aMinCellY; aCellY <= aMaxCellY; ++aCellY) + { + for (int aCellX = aMinCellX; aCellX <= aMaxCellX; ++aCellX) + { + myGrid.SetValue(aCellY * THE_GRID_SIZE + aCellX, GridCell_Boundary); + } + } + } + + // A polygon cannot change classification inside a connected set of cells + // that contains no boundary. Classify one center per component and flood + // the result through the remaining cells. + int aUsableCellCount = 0; + for (int aCellIndex = 0; aCellIndex < aTotalCells; ++aCellIndex) + { + if (myGrid.Value(aCellIndex) != GridCell_Unvisited) + { + continue; + } + + const int aSeedX = aCellIndex % THE_GRID_SIZE; + const int aSeedY = aCellIndex / THE_GRID_SIZE; + const double aSeedCenterX = (static_cast(aSeedX) + 0.5) * aCellSize; + const double aSeedCenterY = (static_cast(aSeedY) + 0.5) * aCellSize; + const GridCell aComponentValue = + internalSiDans(aSeedCenterX, aSeedCenterY) ? GridCell_Inside : GridCell_Outside; + + aCellQueue.Clear(); + aCellQueue.Append(aCellIndex); + myGrid.SetValue(aCellIndex, aComponentValue); + for (size_t aQueueIndex = 0; aQueueIndex < aCellQueue.Size(); ++aQueueIndex) + { + const int aCurrentCell = aCellQueue[aQueueIndex]; + ++aUsableCellCount; + + const int aCurrentX = aCurrentCell % THE_GRID_SIZE; + const int aCurrentY = aCurrentCell / THE_GRID_SIZE; + const int aNeighborCells[4] = {aCurrentX > 0 ? aCurrentCell - 1 : -1, + aCurrentX + 1 < THE_GRID_SIZE ? aCurrentCell + 1 : -1, + aCurrentY > 0 ? aCurrentCell - THE_GRID_SIZE : -1, + aCurrentY + 1 < THE_GRID_SIZE ? aCurrentCell + THE_GRID_SIZE + : -1}; + for (size_t aNeighborIndex = 0; aNeighborIndex < 4; ++aNeighborIndex) + { + const int aNeighborCell = aNeighborCells[aNeighborIndex]; + if (aNeighborCell < 0) + { + continue; + } + if (myGrid.Value(aNeighborCell) == GridCell_Unvisited) + { + myGrid.SetValue(aNeighborCell, aComponentValue); + aCellQueue.Append(aNeighborCell); + } + } + } + } + + // Even a partial cache is useful after the sustained-use threshold: a + // classified cell avoids an O(edges) scan, while boundary cells remain exact. + if (aUsableCellCount == 0) + { + myGrid = NCollection_Array1(); + myGridState.store(GridState::Disabled, std::memory_order_release); + return; + } + myGridState.store(GridState::Ready, std::memory_order_release); + } + catch (const Standard_OutOfMemory&) + { + myGrid = NCollection_Array1(); + myQueryCount.store(0, std::memory_order_relaxed); + myGridState.store(GridState::NotBuilt, std::memory_order_release); + } + catch (const std::bad_alloc&) + { + myGrid = NCollection_Array1(); + myQueryCount.store(0, std::memory_order_relaxed); + myGridState.store(GridState::NotBuilt, std::memory_order_release); + } + catch (...) + { + myGrid = NCollection_Array1(); + myQueryCount.store(0, std::memory_order_relaxed); + myGridState.store(GridState::NotBuilt, std::memory_order_release); + throw; + } +} + +//================================================================================================= + CSLib_Class2d::Result CSLib_Class2d::SiDans(const gp_Pnt2d& thePoint) const { if (myPointsCount == 0) @@ -149,8 +483,8 @@ CSLib_Class2d::Result CSLib_Class2d::SiDans(const gp_Pnt2d& thePoint) const double aY = thePoint.Y(); // Compute tolerance in original coordinate space. - const double aTolU = myTolU * (myUMax - myUMin); - const double aTolV = myTolV * (myVMax - myVMin); + const double aTolU = myOriginalTolU; + const double aTolV = myOriginalTolV; // Quick rejection test for points clearly outside the bounding box. if (aX < (myUMin - aTolU) || aX > (myUMax + aTolU) || aY < (myVMin - aTolV) @@ -160,11 +494,44 @@ CSLib_Class2d::Result CSLib_Class2d::SiDans(const gp_Pnt2d& thePoint) const } // Transform to normalized coordinates. - aX = transformToNormalized(aX, myUMin, myUMax - myUMin); - aY = transformToNormalized(aY, myVMin, myVMax - myVMin); + aX = transformToNormalized(aX, myUMin, myUMax); + aY = transformToNormalized(aY, myVMin, myVMax); + + // Build the acceleration grid only for sustained workloads. Short-lived + // classifiers and small polygons remain on the cheaper exact scan. + GridState aGridState = myGridState.load(std::memory_order_acquire); + if (aGridState == GridState::NotBuilt && myPointsCount >= THE_GRID_MIN_POINTS) + { + const size_t aQueryCount = myQueryCount.fetch_add(1, std::memory_order_relaxed) + 1; + if (aQueryCount >= THE_GRID_BUILD_QUERY_COUNT) + { + buildGridCache(); + aGridState = myGridState.load(std::memory_order_acquire); + } + } + + // Fast-path: conservative grid lookup (O(1) away from polygon edges). + if (aGridState == GridState::Ready && isSafeForDirectArithmetic(aX) + && isSafeForDirectArithmetic(aY) && aX >= 0.0 && aX <= 1.0 && aY >= 0.0 && aY <= 1.0) + { + int aIX = static_cast(aX * THE_GRID_SIZE); + int aIY = static_cast(aY * THE_GRID_SIZE); + aIX = std::clamp(aIX, 0, THE_GRID_SIZE - 1); + aIY = std::clamp(aIY, 0, THE_GRID_SIZE - 1); + const GridCell aCell = myGrid.Value(aIY * THE_GRID_SIZE + aIX); + if (aCell == GridCell_Inside) + { + return Result_Inside; + } + if (aCell == GridCell_Outside) + { + return Result_Outside; + } + // Boundary cells fall through to exact classification. + } // Perform classification with ON detection. - const Result aResult = internalSiDansOuOn(aX, aY); + const Result aResult = internalSiDansOuOn(aX, aY, myTolU, myTolV); if (aResult == Result_Uncertain) { return Result_Uncertain; // ON boundary @@ -196,31 +563,34 @@ CSLib_Class2d::Result CSLib_Class2d::SiDans_OnMode(const gp_Pnt2d& thePoint, return Result_Uncertain; } - double aX = thePoint.X(); - double aY = thePoint.Y(); + double aX = thePoint.X(); + double aY = thePoint.Y(); + const double aTolerance = sanitizeTolerance(theTol); + const double aTolU = normalizeTolerance(aTolerance, myUMin, myUMax); + const double aTolV = normalizeTolerance(aTolerance, myVMin, myVMax); // Quick rejection test. - if (aX < (myUMin - theTol) || aX > (myUMax + theTol) || aY < (myVMin - theTol) - || aY > (myVMax + theTol)) + if (aX < (myUMin - aTolerance) || aX > (myUMax + aTolerance) || aY < (myVMin - aTolerance) + || aY > (myVMax + aTolerance)) { return Result_Outside; } // Transform to normalized coordinates. - aX = transformToNormalized(aX, myUMin, myUMax - myUMin); - aY = transformToNormalized(aY, myVMin, myVMax - myVMin); + aX = transformToNormalized(aX, myUMin, myUMax); + aY = transformToNormalized(aY, myVMin, myVMax); // Perform classification with ON detection. - const Result aResult = internalSiDansOuOn(aX, aY); + const Result aResult = internalSiDansOuOn(aX, aY, aTolU, aTolV); // Check corner points with tolerance. - if (theTol > 0.0) + if (aTolU > 0.0 || aTolV > 0.0) { const bool isInside = (aResult == Result_Inside); - if (isInside != internalSiDans(aX - theTol, aY - theTol) - || isInside != internalSiDans(aX + theTol, aY - theTol) - || isInside != internalSiDans(aX - theTol, aY + theTol) - || isInside != internalSiDans(aX + theTol, aY + theTol)) + if (isInside != internalSiDans(aX - aTolU, aY - aTolV) + || isInside != internalSiDans(aX + aTolU, aY - aTolV) + || isInside != internalSiDans(aX - aTolU, aY + aTolV) + || isInside != internalSiDans(aX + aTolU, aY + aTolV)) { return Result_Uncertain; } @@ -234,16 +604,20 @@ CSLib_Class2d::Result CSLib_Class2d::SiDans_OnMode(const gp_Pnt2d& thePoint, bool CSLib_Class2d::internalSiDans(const double thePx, const double thePy) const { // Ray-casting algorithm: count edge crossings with a horizontal ray from (Px, Py) to +infinity. + // Use raw pointers for cache-friendly sequential access and auto-vectorization. + const double* pX = &myPnts2dX.First(); + const double* pY = &myPnts2dY.First(); + int aNbCrossings = 0; - double aPrevDx = myPnts2dX.Value(0) - thePx; - double aPrevDy = myPnts2dY.Value(0) - thePy; + double aPrevDx = pX[0] - thePx; + double aPrevDy = pY[0] - thePy; bool aPrevYIsNegative = (aPrevDy < 0.0); for (int aNextIdx = 1; aNextIdx <= myPointsCount; ++aNextIdx) { - const double aCurrDx = myPnts2dX.Value(aNextIdx) - thePx; - const double aCurrDy = myPnts2dY.Value(aNextIdx) - thePy; + const double aCurrDx = pX[aNextIdx] - thePx; + const double aCurrDy = pY[aNextIdx] - thePy; const bool aCurrYIsNegative = (aCurrDy < 0.0); // Check for edge crossing when Y changes sign. @@ -277,23 +651,29 @@ bool CSLib_Class2d::internalSiDans(const double thePx, const double thePy) const //================================================================================================= CSLib_Class2d::Result CSLib_Class2d::internalSiDansOuOn(const double thePx, - const double thePy) const + const double thePy, + const double theTolU, + const double theTolV) const { // Ray-casting algorithm with ON detection. + // Use raw pointers for cache-friendly sequential access and auto-vectorization. + const double* pX = &myPnts2dX.First(); + const double* pY = &myPnts2dY.First(); + int aNbCrossings = 0; - double aPrevDx = myPnts2dX.Value(0) - thePx; - double aPrevDy = myPnts2dY.Value(0) - thePy; + double aPrevDx = pX[0] - thePx; + double aPrevDy = pY[0] - thePy; bool aPrevYIsNegative = (aPrevDy < 0.0); for (int aNextIdx = 1; aNextIdx <= myPointsCount; ++aNextIdx) { const int aPrevIdx = aNextIdx - 1; - const double aCurrDx = myPnts2dX.Value(aNextIdx) - thePx; - const double aCurrDy = myPnts2dY.Value(aNextIdx) - thePy; + const double aCurrDx = pX[aNextIdx] - thePx; + const double aCurrDy = pY[aNextIdx] - thePy; // Check if point is very close to current vertex. - if (aCurrDx < myTolU && aCurrDx > -myTolU && aCurrDy < myTolV && aCurrDy > -myTolV) + if (aCurrDx < theTolU && aCurrDx > -theTolU && aCurrDy < theTolV && aCurrDy > -theTolV) { return Result_Uncertain; // ON boundary (at vertex) } @@ -301,15 +681,12 @@ CSLib_Class2d::Result CSLib_Class2d::internalSiDansOuOn(const double thePx, // Check if point is ON the edge by computing Y at the test point's X. // Skip interpolation for nearly vertical edges to avoid division instability. // For vertical edges, the ON detection is handled by the tolerance check above. - const double aEdgeDx = myPnts2dX.Value(aNextIdx) - myPnts2dX.Value(aPrevIdx); - if ((myPnts2dX.Value(aPrevIdx) - thePx) * aCurrDx < 0.0 - && std::abs(aEdgeDx) > Precision::PConfusion()) + const double aEdgeDx = pX[aNextIdx] - pX[aPrevIdx]; + if ((pX[aPrevIdx] - thePx) * aCurrDx < 0.0 && std::abs(aEdgeDx) > Precision::PConfusion()) { - const double aInterpY = - myPnts2dY.Value(aNextIdx) - - (myPnts2dY.Value(aNextIdx) - myPnts2dY.Value(aPrevIdx)) / aEdgeDx * aCurrDx; - const double aDeltaY = aInterpY - thePy; - if (aDeltaY >= -myTolV && aDeltaY <= myTolV) + const double aInterpY = pY[aNextIdx] - (pY[aNextIdx] - pY[aPrevIdx]) / aEdgeDx * aCurrDx; + const double aDeltaY = aInterpY - thePy; + if (aDeltaY >= -theTolV && aDeltaY <= theTolV) { return Result_Uncertain; // ON boundary (on edge) } diff --git a/src/FoundationClasses/TKMath/CSLib/CSLib_Class2d.hxx b/src/FoundationClasses/TKMath/CSLib/CSLib_Class2d.hxx index 8e8d8d2ef4..2233990b63 100644 --- a/src/FoundationClasses/TKMath/CSLib/CSLib_Class2d.hxx +++ b/src/FoundationClasses/TKMath/CSLib/CSLib_Class2d.hxx @@ -24,6 +24,8 @@ #include #include +#include + class gp_Pnt2d; //! Low-level algorithm for 2D point-in-polygon classification. @@ -110,37 +112,17 @@ public: double theUMax, double theVMax); + //! Deep-copy constructor. The immutable polygon and a completed grid cache are copied. + Standard_EXPORT CSLib_Class2d(const CSLib_Class2d& theOther); + + //! Deep-copy assignment. A grid under construction is intentionally not copied. + Standard_EXPORT CSLib_Class2d& operator=(const CSLib_Class2d& theOther); + //! Move constructor. - CSLib_Class2d(CSLib_Class2d&& theOther) noexcept - : myPnts2dX(std::move(theOther.myPnts2dX)), - myPnts2dY(std::move(theOther.myPnts2dY)), - myTolU(theOther.myTolU), - myTolV(theOther.myTolV), - myPointsCount(theOther.myPointsCount), - myUMin(theOther.myUMin), - myVMin(theOther.myVMin), - myUMax(theOther.myUMax), - myVMax(theOther.myVMax) - { - } + Standard_EXPORT CSLib_Class2d(CSLib_Class2d&& theOther) noexcept; //! Move assignment operator. - CSLib_Class2d& operator=(CSLib_Class2d&& theOther) noexcept - { - if (this != &theOther) - { - myPnts2dX = std::move(theOther.myPnts2dX); - myPnts2dY = std::move(theOther.myPnts2dY); - myTolU = theOther.myTolU; - myTolV = theOther.myTolV; - myPointsCount = theOther.myPointsCount; - myUMin = theOther.myUMin; - myVMin = theOther.myVMin; - myUMax = theOther.myUMax; - myVMax = theOther.myVMax; - } - return *this; - } + Standard_EXPORT CSLib_Class2d& operator=(CSLib_Class2d&& theOther) noexcept; //! Classifies a point relative to the polygon. //! @@ -173,10 +155,12 @@ private: //! //! Same as internalSiDans() but also detects if the point lies on the boundary. //! - //! @param[in] theX X coordinate in normalized space - //! @param[in] theY Y coordinate in normalized space + //! @param[in] theX X coordinate in normalized space + //! @param[in] theY Y coordinate in normalized space + //! @param[in] theTolU U tolerance in normalized space + //! @param[in] theTolV V tolerance in normalized space //! @return Classification result - Result internalSiDansOuOn(double theX, double theY) const; + Result internalSiDansOuOn(double theX, double theY, double theTolU, double theTolV) const; //! Initializes the classifier with polygon data. //! @tparam TCol_Containers2d Container type (Array1 or Sequence) @@ -189,22 +173,44 @@ private: double theUMax, double theVMax); - //! Copy constructor is deleted. - CSLib_Class2d(const CSLib_Class2d&) = delete; - - //! Copy assignment operator is deleted. - CSLib_Class2d& operator=(const CSLib_Class2d&) = delete; + //! Builds the grid cache for fast point classification on first sustained use. + //! Cells whose box overlaps a tolerance-expanded polygon-edge box remain on + //! the exact path; only provably boundary-free cells are classified/cached. + void buildGridCache() const; private: - NCollection_Array1 myPnts2dX; //!< X coordinates (normalized) - NCollection_Array1 myPnts2dY; //!< Y coordinates (normalized) - double myTolU = 0.0; //!< Tolerance in U direction (normalized) - double myTolV = 0.0; //!< Tolerance in V direction (normalized) - int myPointsCount = 0; //!< Number of polygon vertices - double myUMin = 0.0; //!< Original minimum U bound - double myVMin = 0.0; //!< Original minimum V bound - double myUMax = 0.0; //!< Original maximum U bound - double myVMax = 0.0; //!< Original maximum V bound + //! Grid cell classification for the fast-path cache. + enum GridCell : signed char + { + GridCell_Outside = -1, //!< Cell is fully outside the polygon + GridCell_Boundary = 0, //!< Cell straddles a polygon edge and needs an exact test + GridCell_Inside = 1, //!< Cell is fully inside the polygon + GridCell_Unvisited = 2 //!< Boundary-free cell not yet assigned by flood fill + }; + + //! Lifecycle of the optional immutable grid cache. + enum class GridState : unsigned char + { + NotBuilt, //!< No cache, including after a transient allocation failure + Building, //!< One thread owns construction; other threads use the exact path + Ready, //!< Grid is complete and immutable + Disabled //!< Polygon or tolerance intrinsically cannot use a grid + }; + + NCollection_Array1 myPnts2dX; //!< X coordinates (normalized) + NCollection_Array1 myPnts2dY; //!< Y coordinates (normalized) + double myTolU = 0.0; //!< U tolerance (normalized) + double myTolV = 0.0; //!< V tolerance (normalized) + double myOriginalTolU = 0.0; //!< U tolerance in input coordinates + double myOriginalTolV = 0.0; //!< V tolerance in input coordinates + int myPointsCount = 0; //!< Number of polygon vertices + double myUMin = 0.0; //!< Original minimum U bound + double myVMin = 0.0; //!< Original minimum V bound + double myUMax = 0.0; //!< Original maximum U bound + double myVMax = 0.0; //!< Original maximum V bound + mutable NCollection_Array1 myGrid; //!< Immutable when Ready + mutable std::atomic myGridState{GridState::NotBuilt}; + mutable std::atomic myQueryCount{0}; }; #endif // _CSLib_Class2d_HeaderFile diff --git a/src/FoundationClasses/TKMath/GTests/CSLib_Test.cxx b/src/FoundationClasses/TKMath/GTests/CSLib_Test.cxx index 7540ea57a9..893d2717ca 100644 --- a/src/FoundationClasses/TKMath/GTests/CSLib_Test.cxx +++ b/src/FoundationClasses/TKMath/GTests/CSLib_Test.cxx @@ -28,7 +28,11 @@ #include #include +#include +#include #include +#include +#include namespace { @@ -254,6 +258,28 @@ TEST_F(CSLibClass2dTest, SiDans_PointOnBoundary) EXPECT_EQ(aClassifier.SiDans(aPointOnEdge), 0); } +TEST_F(CSLibClass2dTest, DeepCopyPreservesClassification) +{ + NCollection_Array1 aPnts(1, 4); + aPnts(1) = gp_Pnt2d(0.0, 0.0); + aPnts(2) = gp_Pnt2d(1.0, 0.0); + aPnts(3) = gp_Pnt2d(1.0, 1.0); + aPnts(4) = gp_Pnt2d(0.0, 1.0); + + CSLib_Class2d aSource(aPnts, 0.01, 0.01, 0.0, 0.0, 1.0, 1.0); + const CSLib_Class2d aCopy(aSource); + CSLib_Class2d anAssigned; + anAssigned = aSource; + + const gp_Pnt2d aSamples[] = {gp_Pnt2d(0.5, 0.5), gp_Pnt2d(2.0, 2.0), gp_Pnt2d(0.5, 0.0)}; + for (const gp_Pnt2d& aPoint : aSamples) + { + const CSLib_Class2d::Result aState = aSource.SiDans(aPoint); + EXPECT_EQ(aCopy.SiDans(aPoint), aState); + EXPECT_EQ(anAssigned.SiDans(aPoint), aState); + } +} + TEST_F(CSLibClass2dTest, SiDans_TriangularPolygon) { NCollection_Array1 aPnts(1, 3); @@ -303,6 +329,266 @@ TEST_F(CSLibClass2dTest, InternalSiDans_NormalizedCoordinates) EXPECT_EQ(aClassifier.SiDans(gp_Pnt2d(15.0, 15.0)), CSLib_Class2d::Result_Outside); } +TEST_F(CSLibClass2dTest, LazyGridRemainsExactAwayFromBoundary) +{ + constexpr int THE_POINT_COUNT = 64; + NCollection_Array1 aPnts(1, THE_POINT_COUNT); + for (int anIdx = 0; anIdx < THE_POINT_COUNT; ++anIdx) + { + const double anAngle = + 2.0 * M_PI * static_cast(anIdx) / static_cast(THE_POINT_COUNT); + aPnts(anIdx + 1) = gp_Pnt2d(0.5 + 0.4 * std::cos(anAngle), 0.5 + 0.4 * std::sin(anAngle)); + } + + CSLib_Class2d aClassifier(aPnts, 1.0e-8, 1.0e-8, 0.0, 0.0, 1.0, 1.0); + + // Cross the lazy-build threshold with exact interior queries. + for (int aQueryIdx = 0; aQueryIdx < 64; ++aQueryIdx) + { + EXPECT_EQ(aClassifier.SiDans(gp_Pnt2d(0.5, 0.5)), CSLib_Class2d::Result_Inside); + } + + // Exercise cached cells, excluding an annulus around the polygon boundary. + for (int aY = 0; aY < 16; ++aY) + { + for (int anX = 0; anX < 16; ++anX) + { + const gp_Pnt2d aPoint((static_cast(anX) + 0.5) / 16.0, + (static_cast(aY) + 0.5) / 16.0); + const double aDx = aPoint.X() - 0.5; + const double aDy = aPoint.Y() - 0.5; + const double aRadius = std::sqrt(aDx * aDx + aDy * aDy); + if (aRadius < 0.35) + { + EXPECT_EQ(aClassifier.SiDans(aPoint), CSLib_Class2d::Result_Inside); + } + else if (aRadius > 0.45) + { + EXPECT_EQ(aClassifier.SiDans(aPoint), CSLib_Class2d::Result_Outside); + } + } + } + + // Boundary candidates must stay on the exact/tolerance path after caching. + EXPECT_EQ(aClassifier.SiDans(aPnts(1)), CSLib_Class2d::Result_Uncertain); + + const CSLib_Class2d aCopy(aClassifier); + CSLib_Class2d anAssigned; + anAssigned = aClassifier; + EXPECT_EQ(aCopy.SiDans(gp_Pnt2d(0.5, 0.5)), CSLib_Class2d::Result_Inside); + EXPECT_EQ(anAssigned.SiDans(gp_Pnt2d(0.99, 0.99)), CSLib_Class2d::Result_Outside); + EXPECT_EQ(aCopy.SiDans(aPnts(1)), CSLib_Class2d::Result_Uncertain); +} + +TEST_F(CSLibClass2dTest, LazyGridMatchesExactPathForConcaveDiagonalPolygon) +{ + constexpr int THE_POINT_COUNT = 32; + NCollection_Array1 aPnts(1, THE_POINT_COUNT); + for (int anIdx = 0; anIdx < THE_POINT_COUNT; ++anIdx) + { + const double anAngle = + 2.0 * M_PI * static_cast(anIdx) / static_cast(THE_POINT_COUNT); + const double aRadius = (anIdx % 4 == 1) ? 0.22 : ((anIdx % 2 == 0) ? 0.46 : 0.34); + aPnts(anIdx + 1) = + gp_Pnt2d(0.5 + aRadius * std::cos(anAngle), 0.5 + aRadius * std::sin(anAngle)); + } + + CSLib_Class2d aCached(aPnts, 0.0, 0.0, 0.0, 0.0, 1.0, 1.0); + CSLib_Class2d anExact(aPnts, 0.0, 0.0, 0.0, 0.0, 1.0, 1.0); + for (size_t aQueryIdx = 0; aQueryIdx < 64; ++aQueryIdx) + { + ASSERT_EQ(aCached.SiDans(gp_Pnt2d(0.5, 0.5)), CSLib_Class2d::Result_Inside); + } + + for (int aY = 0; aY < 47; ++aY) + { + for (int anX = 0; anX < 47; ++anX) + { + const gp_Pnt2d aPoint((static_cast(anX) + 0.37) / 47.0, + (static_cast(aY) + 0.61) / 47.0); + EXPECT_EQ(aCached.SiDans(aPoint), anExact.SiDans_OnMode(aPoint, 0.0)) + << "sample (" << anX << ", " << aY << ")"; + } + } +} + +TEST_F(CSLibClass2dTest, LazyGridConcurrentBuildOverlapIsThreadSafe) +{ + // The large polygon keeps the designated build query active while synchronized + // readers exercise the exact fallback. Run this test under TSAN for race detection. + constexpr int THE_POINT_COUNT = 32768; + NCollection_Array1 aPnts(1, THE_POINT_COUNT); + for (int anIdx = 0; anIdx < THE_POINT_COUNT; ++anIdx) + { + const double anAngle = + 2.0 * M_PI * static_cast(anIdx) / static_cast(THE_POINT_COUNT); + aPnts(anIdx + 1) = gp_Pnt2d(0.5 + 0.4 * std::cos(anAngle), 0.5 + 0.4 * std::sin(anAngle)); + } + + CSLib_Class2d aClassifier(aPnts, 1.0e-8, 1.0e-8, 0.0, 0.0, 1.0, 1.0); + for (size_t aQueryIdx = 0; aQueryIdx < 63; ++aQueryIdx) + { + ASSERT_EQ(aClassifier.SiDans(gp_Pnt2d(0.5, 0.5)), CSLib_Class2d::Result_Inside); + } + + std::atomic hasFailure{false}; + std::atomic canStart{false}; + std::atomic hasBuildQueryStarted{false}; + std::atomic hasBuildQueryFinished{false}; + std::atomic anOverlapQueryCount{0}; + std::atomic aReadyCount{0}; + std::thread aBuildThread([&]() { + while (!canStart.load(std::memory_order_acquire)) + { + std::this_thread::yield(); + } + hasBuildQueryStarted.store(true, std::memory_order_release); + if (aClassifier.SiDans(gp_Pnt2d(0.5, 0.5)) != CSLib_Class2d::Result_Inside) + { + hasFailure.store(true, std::memory_order_relaxed); + } + hasBuildQueryFinished.store(true, std::memory_order_release); + }); + + std::array aReaderThreads; + for (size_t aThreadIdx = 0; aThreadIdx < aReaderThreads.size(); ++aThreadIdx) + { + aReaderThreads[aThreadIdx] = std::thread([&]() { + aReadyCount.fetch_add(1, std::memory_order_release); + while (!hasBuildQueryStarted.load(std::memory_order_acquire)) + { + std::this_thread::yield(); + } + for (size_t aQueryIdx = 0; aQueryIdx < 8; ++aQueryIdx) + { + if (!hasBuildQueryFinished.load(std::memory_order_acquire)) + { + anOverlapQueryCount.fetch_add(1, std::memory_order_relaxed); + } + const bool isInside = (aQueryIdx & 1u) == 0u; + const gp_Pnt2d aPoint = isInside ? gp_Pnt2d(0.5, 0.5) : gp_Pnt2d(0.99, 0.99); + const CSLib_Class2d::Result anExpected = + isInside ? CSLib_Class2d::Result_Inside : CSLib_Class2d::Result_Outside; + if (aClassifier.SiDans(aPoint) != anExpected) + { + hasFailure.store(true, std::memory_order_relaxed); + return; + } + } + }); + } + while (aReadyCount.load(std::memory_order_acquire) != aReaderThreads.size()) + { + std::this_thread::yield(); + } + canStart.store(true, std::memory_order_release); + while (!hasBuildQueryStarted.load(std::memory_order_acquire)) + { + std::this_thread::yield(); + } + if (!hasBuildQueryFinished.load(std::memory_order_acquire)) + { + anOverlapQueryCount.fetch_add(1, std::memory_order_relaxed); + if (aClassifier.SiDans(gp_Pnt2d(0.99, 0.99)) != CSLib_Class2d::Result_Outside) + { + hasFailure.store(true, std::memory_order_relaxed); + } + } + aBuildThread.join(); + for (std::thread& aThread : aReaderThreads) + { + aThread.join(); + } + EXPECT_FALSE(hasFailure.load(std::memory_order_relaxed)); + EXPECT_GT(anOverlapQueryCount.load(std::memory_order_relaxed), 0u); +} + +TEST_F(CSLibClass2dTest, SiDansHandlesNegativeNonFiniteAndExtremeInputs) +{ + constexpr int THE_POINT_COUNT = 32; + constexpr int THE_POINTS_PER_SIDE = THE_POINT_COUNT / 4; + NCollection_Array1 aPnts(1, THE_POINT_COUNT); + for (int anIdx = 0; anIdx < THE_POINTS_PER_SIDE; ++anIdx) + { + const double aParameter = static_cast(anIdx) / THE_POINTS_PER_SIDE; + aPnts(1 + anIdx) = gp_Pnt2d(aParameter, 0.0); + aPnts(1 + THE_POINTS_PER_SIDE + anIdx) = gp_Pnt2d(1.0, aParameter); + aPnts(1 + 2 * THE_POINTS_PER_SIDE + anIdx) = gp_Pnt2d(1.0 - aParameter, 1.0); + aPnts(1 + 3 * THE_POINTS_PER_SIDE + anIdx) = gp_Pnt2d(0.0, 1.0 - aParameter); + } + + CSLib_Class2d + aNegative(aPnts, -1.0, -std::numeric_limits::infinity(), 0.0, 0.0, 1.0, 1.0); + EXPECT_EQ(aNegative.SiDans(gp_Pnt2d(0.5, 0.5)), CSLib_Class2d::Result_Inside); + EXPECT_EQ(aNegative.SiDans_OnMode(gp_Pnt2d(1.5, 0.5), -1.0), CSLib_Class2d::Result_Outside); + + CSLib_Class2d anExtremeTolerance(aPnts, + std::numeric_limits::infinity(), + std::numeric_limits::max(), + 0.0, + 0.0, + 1.0, + 1.0); + for (size_t aQueryIdx = 0; aQueryIdx < 80; ++aQueryIdx) + { + EXPECT_EQ(anExtremeTolerance.SiDans(gp_Pnt2d(0.5, 0.5)), CSLib_Class2d::Result_Uncertain); + } + + const double aLimit = Precision::Infinite(); + NCollection_Array1 anExtremePnts(1, 4); + anExtremePnts(1) = gp_Pnt2d(-0.5 * aLimit, -0.5 * aLimit); + anExtremePnts(2) = gp_Pnt2d(0.5 * aLimit, -0.5 * aLimit); + anExtremePnts(3) = gp_Pnt2d(0.5 * aLimit, 0.5 * aLimit); + anExtremePnts(4) = gp_Pnt2d(-0.5 * aLimit, 0.5 * aLimit); + CSLib_Class2d anExtremeCoordinates(anExtremePnts, 0.0, 0.0, -aLimit, -aLimit, aLimit, aLimit); + EXPECT_EQ(anExtremeCoordinates.SiDans(gp_Pnt2d(0.0, 0.0)), CSLib_Class2d::Result_Inside); + EXPECT_EQ(anExtremeCoordinates.SiDans(gp_Pnt2d(0.75 * aLimit, 0.0)), + CSLib_Class2d::Result_Outside); + + const double aLargeTolerance = 0.1 * aLimit; + CSLib_Class2d aTolerantExtreme(anExtremePnts, + aLargeTolerance, + aLargeTolerance, + -aLimit, + -aLimit, + aLimit, + aLimit); + EXPECT_EQ(aTolerantExtreme.SiDans(gp_Pnt2d(0.55 * aLimit, 0.0)), CSLib_Class2d::Result_Uncertain); + EXPECT_EQ(aTolerantExtreme.SiDans(gp_Pnt2d(0.65 * aLimit, 0.0)), CSLib_Class2d::Result_Outside); + EXPECT_EQ(anExtremeCoordinates.SiDans_OnMode(gp_Pnt2d(0.0, 0.55 * aLimit), aLargeTolerance), + CSLib_Class2d::Result_Uncertain); +} + +TEST_F(CSLibClass2dTest, LazyGridPreservesToleranceAcrossCellBoundary) +{ + constexpr int THE_POINTS_PER_SIDE = 8; + constexpr int THE_POINT_COUNT = 4 * THE_POINTS_PER_SIDE; + constexpr double THE_MIN = 0.25; + constexpr double THE_MAX = 0.75; + constexpr double THE_TOLERANCE = 1.0e-4; + NCollection_Array1 aPnts(1, THE_POINT_COUNT); + for (int anIdx = 0; anIdx < THE_POINTS_PER_SIDE; ++anIdx) + { + const double aParameter = static_cast(anIdx) / static_cast(THE_POINTS_PER_SIDE); + aPnts(1 + anIdx) = gp_Pnt2d(THE_MIN + (THE_MAX - THE_MIN) * aParameter, THE_MIN); + aPnts(1 + THE_POINTS_PER_SIDE + anIdx) = + gp_Pnt2d(THE_MAX, THE_MIN + (THE_MAX - THE_MIN) * aParameter); + aPnts(1 + 2 * THE_POINTS_PER_SIDE + anIdx) = + gp_Pnt2d(THE_MAX - (THE_MAX - THE_MIN) * aParameter, THE_MAX); + aPnts(1 + 3 * THE_POINTS_PER_SIDE + anIdx) = + gp_Pnt2d(THE_MIN, THE_MAX - (THE_MAX - THE_MIN) * aParameter); + } + + CSLib_Class2d aClassifier(aPnts, THE_TOLERANCE, THE_TOLERANCE, 0.0, 0.0, 1.0, 1.0); + for (size_t aQueryIdx = 0; aQueryIdx < 64; ++aQueryIdx) + { + ASSERT_EQ(aClassifier.SiDans(gp_Pnt2d(0.5, 0.5)), CSLib_Class2d::Result_Inside); + } + + EXPECT_EQ(aClassifier.SiDans(gp_Pnt2d(THE_MIN - 0.5 * THE_TOLERANCE, 0.5)), + CSLib_Class2d::Result_Uncertain); +} + // Test SiDans_OnMode TEST_F(CSLibClass2dTest, SiDans_OnMode_PointInside) { @@ -318,6 +604,25 @@ TEST_F(CSLibClass2dTest, SiDans_OnMode_PointInside) EXPECT_EQ(aClassifier.SiDans_OnMode(aPointInside, 0.01), 1); } +TEST_F(CSLibClass2dTest, SiDans_OnMode_NormalizesAndReplacesToleranceOnAnisotropicDomain) +{ + NCollection_Array1 aPnts(1, 4); + aPnts(1) = gp_Pnt2d(20.0, -1.0); + aPnts(2) = gp_Pnt2d(80.0, -1.0); + aPnts(3) = gp_Pnt2d(80.0, 1.0); + aPnts(4) = gp_Pnt2d(20.0, 1.0); + + // Constructor tolerances are deliberately larger than the explicit one. + CSLib_Class2d aClassifier(aPnts, 10.0, 1.0, 0.0, -2.0, 100.0, 2.0); + EXPECT_EQ(aClassifier.SiDans_OnMode(gp_Pnt2d(50.0, 1.5), 0.01), CSLib_Class2d::Result_Outside); + + // The same original-space tolerance has different normalized U/V values. + EXPECT_EQ(aClassifier.SiDans_OnMode(gp_Pnt2d(19.75, 0.0), 0.5), CSLib_Class2d::Result_Uncertain); + EXPECT_EQ(aClassifier.SiDans_OnMode(gp_Pnt2d(19.0, 0.0), 0.5), CSLib_Class2d::Result_Outside); + EXPECT_EQ(aClassifier.SiDans_OnMode(gp_Pnt2d(50.0, 1.25), 0.5), CSLib_Class2d::Result_Uncertain); + EXPECT_EQ(aClassifier.SiDans_OnMode(gp_Pnt2d(50.0, 1.75), 0.5), CSLib_Class2d::Result_Outside); +} + // Test with degenerate polygon (less than 3 points effective) TEST_F(CSLibClass2dTest, DegeneratePolygon_InvalidBounds) { @@ -401,10 +706,12 @@ TEST_F(CSLibNormalPolyDefTest, Value_AtSingularPoints) // the tolerance check RealSmall(). // Test that the function doesn't crash at these points. EXPECT_TRUE(aPoly.Value(0.0, aValue)); - EXPECT_TRUE(std::isfinite(aValue)); + EXPECT_FALSE(std::isnan(aValue)); + EXPECT_FALSE(Precision::IsInfinite(aValue)); EXPECT_TRUE(aPoly.Value(M_PI / 2.0, aValue)); - EXPECT_TRUE(std::isfinite(aValue)); + EXPECT_FALSE(std::isnan(aValue)); + EXPECT_FALSE(Precision::IsInfinite(aValue)); } // Test Derivative function @@ -421,7 +728,8 @@ TEST_F(CSLibNormalPolyDefTest, Derivative_AtRegularPoint) double aDeriv; EXPECT_TRUE(aPoly.Derivative(M_PI / 4.0, aDeriv)); // Derivative should be computed without crash - EXPECT_TRUE(std::isfinite(aDeriv)); + EXPECT_FALSE(std::isnan(aDeriv)); + EXPECT_FALSE(Precision::IsInfinite(aDeriv)); } TEST_F(CSLibNormalPolyDefTest, Derivative_AtSingularPoint) diff --git a/src/ModelingAlgorithms/TKBool/TopOpeBRepBuild/TopOpeBRepBuild_GridSS.cxx b/src/ModelingAlgorithms/TKBool/TopOpeBRepBuild/TopOpeBRepBuild_GridSS.cxx index c53b6e6f80..3cc1557876 100644 --- a/src/ModelingAlgorithms/TKBool/TopOpeBRepBuild/TopOpeBRepBuild_GridSS.cxx +++ b/src/ModelingAlgorithms/TKBool/TopOpeBRepBuild/TopOpeBRepBuild_GridSS.cxx @@ -1657,7 +1657,7 @@ static bool AreFacesCoincideInArea(const TopoDS_Shape& theBase double tol2d = Precision::PConfusion(); BRepClass_Intersector anInter; BRepClass_Edge aBCE; - aBCE.Face() = aBaseFace; + aBCE.SetFace(aBaseFace); double maxDist = std::max(BRep_Tool::Tolerance(aBaseFace), BRep_Tool::Tolerance(aFace)); bool isError = false; @@ -1687,7 +1687,7 @@ static bool AreFacesCoincideInArea(const TopoDS_Shape& theBase } BB.UpdateEdge(aE, PC, aBaseFace, tolE); } - aBCE.Edge() = aE; + aBCE.SetEdge(aE); anInter.Perform(aLin, pLinMin, tol2d, aBCE); if (anInter.IsDone()) { diff --git a/src/ModelingAlgorithms/TKBool/TopOpeBRepBuild/TopOpeBRepBuild_WireEdgeClassifier.cxx b/src/ModelingAlgorithms/TKBool/TopOpeBRepBuild/TopOpeBRepBuild_WireEdgeClassifier.cxx index 7212721d05..2d7cc04193 100644 --- a/src/ModelingAlgorithms/TKBool/TopOpeBRepBuild/TopOpeBRepBuild_WireEdgeClassifier.cxx +++ b/src/ModelingAlgorithms/TKBool/TopOpeBRepBuild/TopOpeBRepBuild_WireEdgeClassifier.cxx @@ -54,7 +54,7 @@ TopOpeBRepBuild_WireEdgeClassifier::TopOpeBRepBuild_WireEdgeClassifier( const TopOpeBRepBuild_BlockBuilder& BB) : TopOpeBRepBuild_CompositeClassifier(BB) { - myBCEdge.Face() = TopoDS::Face(F); + myBCEdge.SetFace(TopoDS::Face(F)); } //================================================================================================= @@ -509,7 +509,7 @@ bool TopOpeBRepBuild_WireEdgeClassifier::CompareElement(const TopoDS_Shape& EE) myFirstCompare = false; } - myBCEdge.Edge() = E; + myBCEdge.SetEdge(E); TopAbs_Orientation Eori = E.Orientation(); myFPC.Compare(myBCEdge, Eori); #ifdef OCCT_DEBUG diff --git a/src/ModelingAlgorithms/TKTopAlgo/BRepClass/BRepClass_Edge.cxx b/src/ModelingAlgorithms/TKTopAlgo/BRepClass/BRepClass_Edge.cxx index 3a61a06200..063d071f59 100644 --- a/src/ModelingAlgorithms/TKTopAlgo/BRepClass/BRepClass_Edge.cxx +++ b/src/ModelingAlgorithms/TKTopAlgo/BRepClass/BRepClass_Edge.cxx @@ -15,49 +15,33 @@ // commercial license or contractual agreement. #include -#include #include -#include -#include -#include //================================================================================================= BRepClass_Edge::BRepClass_Edge() - : myMaxTolerance(Precision::Infinite()), + : myFirstParameter(0.0), + myLastParameter(0.0), + myMaxTolerance(Precision::Infinite()), + myBoundingBoxState(BndBoxState::NotBuilt), myUseBndBox(false) { } //================================================================================================= -void BRepClass_Edge::SetNextEdge( - const NCollection_IndexedDataMap, - TopTools_ShapeMapHasher>& theMapVE) +void BRepClass_Edge::SetEdge(const TopoDS_Edge& theEdge) { - if (theMapVE.IsEmpty() || myEdge.IsNull()) - { - return; - } - TopoDS_Vertex aVF, aVL; - TopExp::Vertices(myEdge, aVF, aVL, true); + invalidateDerivedData(); + myEdge = theEdge; +} - if (aVL.IsNull() || aVL.IsSame(aVF)) - { - return; - } - const NCollection_List* aListE = theMapVE.Seek(aVL); - if (aListE->Extent() == 2) - { - for (NCollection_List::Iterator anIt(*aListE); anIt.More(); anIt.Next()) - { - if ((!anIt.Value().IsNull()) && (!anIt.Value().IsSame(myEdge))) - { - myNextEdge = TopoDS::Edge(anIt.Value()); - } - } - } +//================================================================================================= + +void BRepClass_Edge::SetFace(const TopoDS_Face& theFace) +{ + invalidateDerivedData(); + myFace = theFace; } //================================================================================================= @@ -65,7 +49,56 @@ void BRepClass_Edge::SetNextEdge( BRepClass_Edge::BRepClass_Edge(const TopoDS_Edge& E, const TopoDS_Face& F) : myEdge(E), myFace(F), + myFirstParameter(0.0), + myLastParameter(0.0), myMaxTolerance(Precision::Infinite()), + myBoundingBoxState(BndBoxState::NotBuilt), myUseBndBox(false) { } + +//================================================================================================= + +void BRepClass_Edge::SetGeometry(const occ::handle& theCurve, + const double theFirst, + const double theLast) +{ + myCurve = theCurve; + myFirstParameter = theFirst; + myLastParameter = theLast; + myBoundingBox.SetVoid(); + myBoundingBoxState = BndBoxState::NotBuilt; +} + +//================================================================================================= + +void BRepClass_Edge::SetBoundingBox(const Bnd_Box2d& theBox) +{ + if (theBox.IsVoid()) + { + SetBoundingBoxUnavailable(); + return; + } + myBoundingBox = theBox; + myBoundingBoxState = BndBoxState::Ready; +} + +//================================================================================================= + +void BRepClass_Edge::SetBoundingBoxUnavailable() +{ + myBoundingBox.SetVoid(); + myBoundingBoxState = BndBoxState::Unavailable; +} + +//================================================================================================= + +void BRepClass_Edge::invalidateDerivedData() +{ + myNextEdge.Nullify(); + myCurve.Nullify(); + myBoundingBox.SetVoid(); + myFirstParameter = 0.0; + myLastParameter = 0.0; + myBoundingBoxState = BndBoxState::NotBuilt; +} diff --git a/src/ModelingAlgorithms/TKTopAlgo/BRepClass/BRepClass_Edge.hxx b/src/ModelingAlgorithms/TKTopAlgo/BRepClass/BRepClass_Edge.hxx index 1f54e6e781..0d4ddb5d36 100644 --- a/src/ModelingAlgorithms/TKTopAlgo/BRepClass/BRepClass_Edge.hxx +++ b/src/ModelingAlgorithms/TKTopAlgo/BRepClass/BRepClass_Edge.hxx @@ -17,17 +17,14 @@ #ifndef _BRepClass_Edge_HeaderFile #define _BRepClass_Edge_HeaderFile -#include +#include #include #include -#include -#include -#include -#include - #include #include +class Geom2d_Curve; + //! This class is used to send the description of an //! Edge to the classifier. It contains an Edge and a //! Face. So the PCurve of the Edge can be found. @@ -36,28 +33,38 @@ class BRepClass_Edge public: DEFINE_STANDARD_ALLOC + //! State of the cached pcurve bounding box. + enum class BndBoxState + { + NotBuilt, + Ready, + Unavailable + }; + Standard_EXPORT BRepClass_Edge(); Standard_EXPORT BRepClass_Edge(const TopoDS_Edge& E, const TopoDS_Face& F); - //! Returns the current Edge - TopoDS_Edge& Edge() { return myEdge; } - + //! Returns the current edge. const TopoDS_Edge& Edge() const { return myEdge; } - //! Returns the Face for the current Edge - TopoDS_Face& Face() { return myFace; } - + //! Returns the face for the current edge. const TopoDS_Face& Face() const { return myFace; } + //! Sets the current edge and invalidates topology-derived data. + //! @param[in] theEdge new edge + Standard_EXPORT void SetEdge(const TopoDS_Edge& theEdge); + + //! Sets the face and invalidates topology-derived data. + //! @param[in] theFace new face + Standard_EXPORT void SetFace(const TopoDS_Face& theFace); + //! Returns the next Edge const TopoDS_Edge& NextEdge() const { return myNextEdge; } - //! Finds and sets the next Edge for the current - Standard_EXPORT void SetNextEdge( - const NCollection_IndexedDataMap, - TopTools_ShapeMapHasher>& theMapVE); + //! Sets the next edge at the last vertex of the current edge. + //! @param[in] theEdge next edge + void SetNextEdge(const TopoDS_Edge& theEdge) { myNextEdge = theEdge; } //! Returns the maximum tolerance double MaxTolerance() const { return myMaxTolerance; } @@ -74,12 +81,50 @@ public: //! using boxes or not void SetUseBndBox(const bool theValue) { myUseBndBox = theValue; } + //! Sets cached 2D geometry and invalidates its bounding box. + //! @param[in] theCurve pcurve on the associated face + //! @param[in] theFirst first pcurve parameter + //! @param[in] theLast last pcurve parameter + Standard_EXPORT void SetGeometry(const occ::handle& theCurve, + double theFirst, + double theLast); + + //! Sets a successfully computed pcurve bounding box. A void box marks the box unavailable. + //! @param[in] theBox pcurve bounding box + Standard_EXPORT void SetBoundingBox(const Bnd_Box2d& theBox); + + //! Marks the pcurve bounding box as unavailable after a failed build. + Standard_EXPORT void SetBoundingBoxUnavailable(); + + //! Returns cached pcurve, or null when it is unavailable. + const occ::handle& Curve() const { return myCurve; } + + //! Returns cached first pcurve parameter. + double FirstParameter() const { return myFirstParameter; } + + //! Returns cached last pcurve parameter. + double LastParameter() const { return myLastParameter; } + + //! Returns cached pcurve bounding box. + const Bnd_Box2d& BoundingBox() const { return myBoundingBox; } + + //! Returns the state of the cached pcurve bounding box. + BndBoxState BoundingBoxState() const { return myBoundingBoxState; } + private: - TopoDS_Edge myEdge; - TopoDS_Face myFace; - TopoDS_Edge myNextEdge; - double myMaxTolerance; - bool myUseBndBox; + void invalidateDerivedData(); + +private: + TopoDS_Edge myEdge; + TopoDS_Face myFace; + TopoDS_Edge myNextEdge; + occ::handle myCurve; + Bnd_Box2d myBoundingBox; + double myFirstParameter; + double myLastParameter; + double myMaxTolerance; + BndBoxState myBoundingBoxState; + bool myUseBndBox; }; #endif // _BRepClass_Edge_HeaderFile diff --git a/src/ModelingAlgorithms/TKTopAlgo/BRepClass/BRepClass_FaceExplorer.cxx b/src/ModelingAlgorithms/TKTopAlgo/BRepClass/BRepClass_FaceExplorer.cxx index 5d4a4201dd..b318b480da 100644 --- a/src/ModelingAlgorithms/TKTopAlgo/BRepClass/BRepClass_FaceExplorer.cxx +++ b/src/ModelingAlgorithms/TKTopAlgo/BRepClass/BRepClass_FaceExplorer.cxx @@ -18,24 +18,83 @@ // Total rewriting of the method Segment; add the method OtherSegment. #include +#include +#include #include #include #include #include +#include +#include #include +#include +#include #include #include +#include +#include #include static const double Probing_Start = 0.123; static const double Probing_End = 0.7; static const double Probing_Step = 0.2111; +namespace +{ + +constexpr size_t THE_MIN_EDGES_FOR_BOUNDING_BOX = 10; + +void cacheGeometry(BRepClass_Edge& theEdge) +{ + const BRepClass_Edge& anEdge = theEdge; + double aFirst = 0.0; + double aLast = 0.0; + const occ::handle& aCurve = + BRep_Tool::CurveOnSurface(anEdge.Edge(), anEdge.Face(), aFirst, aLast); + if (aCurve.IsNull()) + { + return; + } + + theEdge.SetGeometry(aCurve, aFirst, aLast); +} + +struct VertexEdges +{ + TopoDS_Edge First; + TopoDS_Edge Second; + uint32_t Count = 0; + + void Add(const TopoDS_Edge& theEdge) + { + if (Count == 0) + { + First = theEdge; + } + else if (Count == 1) + { + Second = theEdge; + } + ++Count; + } +}; + +struct EdgeOccurrences +{ + uint32_t First = 0; + uint32_t Last = 0; +}; + +} // namespace + //================================================================================================= BRepClass_FaceExplorer::BRepClass_FaceExplorer(const TopoDS_Face& F) : myFace(F), - myCurEdgeInd(1), + myCurrentWire(0), + myCurrentEdge(0), + myCurrentEdgeEnd(0), + myCurEdgeInd(0), myCurEdgePar(Probing_Start), myMaxTolerance(0.1), myUseBndBox(false), @@ -46,6 +105,113 @@ BRepClass_FaceExplorer::BRepClass_FaceExplorer(const TopoDS_Face& F) { myFace.Orientation(TopAbs_FORWARD); + + constexpr uint32_t THE_NO_EDGE = UINT32_MAX; + NCollection_DataMap anOccurrences; + NCollection_LinearVector aNextOccurrence; + + for (TopExp_Explorer aWireExp(myFace, TopAbs_WIRE); aWireExp.More(); aWireExp.Next()) + { + WireData aWire; + aWire.FirstEdge = static_cast(myEdges.Size()); + + NCollection_DataMap aVertexEdges; + for (TopExp_Explorer anEdgeExp(aWireExp.Current(), TopAbs_EDGE); anEdgeExp.More(); + anEdgeExp.Next()) + { + const TopoDS_Edge& aTopoEdge = TopoDS::Edge(anEdgeExp.Current()); + BRepClass_Edge anEdgeData(aTopoEdge, myFace); + cacheGeometry(anEdgeData); + myEdges.Append(std::move(anEdgeData)); + const uint32_t anEdgeIndex = static_cast(myEdges.Size() - 1); + + aNextOccurrence.Append(THE_NO_EDGE); + EdgeOccurrences* anOccurrence = anOccurrences.ChangeSeek(aTopoEdge); + if (anOccurrence == nullptr) + { + anOccurrences.Bind(aTopoEdge, EdgeOccurrences{anEdgeIndex, anEdgeIndex}); + } + else + { + aNextOccurrence[anOccurrence->Last] = anEdgeIndex; + anOccurrence->Last = anEdgeIndex; + } + + for (TopExp_Explorer aVertexExp(aTopoEdge, TopAbs_VERTEX); aVertexExp.More(); + aVertexExp.Next()) + { + const TopoDS_Shape& aVertex = aVertexExp.Current(); + VertexEdges* anEdges = aVertexEdges.ChangeSeek(aVertex); + if (anEdges == nullptr) + { + VertexEdges aNewEdges; + aNewEdges.Add(aTopoEdge); + aVertexEdges.Bind(aVertex, std::move(aNewEdges)); + } + else + { + anEdges->Add(aTopoEdge); + } + } + ++aWire.NbEdges; + } + + for (uint32_t anEdgeIndex = aWire.FirstEdge; anEdgeIndex < aWire.FirstEdge + aWire.NbEdges; + ++anEdgeIndex) + { + const BRepClass_Edge& anEdgeData = myEdges[anEdgeIndex]; + TopoDS_Vertex aFirstVertex; + TopoDS_Vertex aLastVertex; + TopExp::Vertices(anEdgeData.Edge(), aFirstVertex, aLastVertex, true); + if (aLastVertex.IsNull() || aLastVertex.IsSame(aFirstVertex)) + { + continue; + } + + const VertexEdges* anEdges = aVertexEdges.Seek(aLastVertex); + if (anEdges == nullptr || anEdges->Count != 2) + { + continue; + } + const TopoDS_Edge& aNextEdge = + anEdges->First.IsSame(anEdgeData.Edge()) ? anEdges->Second : anEdges->First; + if (!aNextEdge.IsNull() && !aNextEdge.IsSame(anEdgeData.Edge())) + { + myEdges[anEdgeIndex].SetNextEdge(aNextEdge); + } + } + myWires.Append(aWire); + } + + for (TopExp_Explorer anEdgeExp(myFace, TopAbs_EDGE); anEdgeExp.More(); anEdgeExp.Next()) + { + const TopoDS_Edge& anEdge = TopoDS::Edge(anEdgeExp.Current()); + uint32_t anEdgeIndex = 0; + bool isFound = false; + const EdgeOccurrences* anOccurrence = anOccurrences.Seek(anEdge); + if (anOccurrence != nullptr) + { + for (uint32_t aCandidateIndex = anOccurrence->First; aCandidateIndex != THE_NO_EDGE; + aCandidateIndex = aNextOccurrence[aCandidateIndex]) + { + const BRepClass_Edge& aCandidate = myEdges[aCandidateIndex]; + if (aCandidate.Edge().IsEqual(anEdge)) + { + anEdgeIndex = aCandidateIndex; + isFound = true; + break; + } + } + } + if (!isFound) + { + BRepClass_Edge anEdgeData(anEdge, myFace); + cacheGeometry(anEdgeData); + myEdges.Append(std::move(anEdgeData)); + anEdgeIndex = static_cast(myEdges.Size() - 1); + } + myProbeEdges.Append(anEdgeIndex); + } } //================================================================================================= @@ -110,7 +276,7 @@ bool BRepClass_FaceExplorer::Reject(const gp_Pnt2d&) const bool BRepClass_FaceExplorer::Segment(const gp_Pnt2d& P, gp_Lin2d& L, double& Par) { - myCurEdgeInd = 1; + myCurEdgeInd = 0; myCurEdgePar = Probing_Start; return OtherSegment(P, L, Par); @@ -120,30 +286,19 @@ bool BRepClass_FaceExplorer::Segment(const gp_Pnt2d& P, gp_Lin2d& L, double& Par bool BRepClass_FaceExplorer::OtherSegment(const gp_Pnt2d& P, gp_Lin2d& L, double& Par) { - TopExp_Explorer anExpF(myFace, TopAbs_EDGE); - int i; - double aFPar; - double aLPar; - occ::handle aC2d; - constexpr double aTolParConf2 = Precision::PConfusion() * Precision::PConfusion(); - gp_Pnt2d aPOnC; - double aParamIn; - - for (i = 1; anExpF.More(); anExpF.Next(), i++) + constexpr double aTolParConf2 = Precision::PConfusion() * Precision::PConfusion(); + gp_Pnt2d aPOnC; + while (myCurEdgeInd < myProbeEdges.Size()) { - if (i != myCurEdgeInd) - { - continue; - } - - const TopoDS_Shape& aLocalShape = anExpF.Current(); - const TopAbs_Orientation anOrientation = aLocalShape.Orientation(); + const BRepClass_Edge& anEdgeData = myEdges[myProbeEdges[myCurEdgeInd]]; + const TopoDS_Edge& anEdge = anEdgeData.Edge(); + const TopAbs_Orientation anOrientation = anEdge.Orientation(); if (anOrientation == TopAbs_FORWARD || anOrientation == TopAbs_REVERSED) { - const TopoDS_Edge& anEdge = TopoDS::Edge(aLocalShape); - - aC2d = BRep_Tool::CurveOnSurface(anEdge, myFace, aFPar, aLPar); + const occ::handle& aC2d = anEdgeData.Curve(); + double aFPar = anEdgeData.FirstParameter(); + double aLPar = anEdgeData.LastParameter(); if (!aC2d.IsNull()) { @@ -167,7 +322,7 @@ bool BRepClass_FaceExplorer::OtherSegment(const gp_Pnt2d& P, gp_Lin2d& L, double for (; myCurEdgePar < Probing_End; myCurEdgePar += Probing_Step) { - aParamIn = myCurEdgePar * aFPar + (1. - myCurEdgePar) * aLPar; + const double aParamIn = myCurEdgePar * aFPar + (1. - myCurEdgePar) * aLPar; gp_Vec2d aTanVec; aC2d->D1(aParamIn, aPOnC, aTanVec); @@ -254,7 +409,7 @@ bool BRepClass_FaceExplorer::OtherSegment(const gp_Pnt2d& P, gp_Lin2d& L, double myCurEdgePar += Probing_Step; if (myCurEdgePar >= Probing_End) { - myCurEdgeInd++; + ++myCurEdgeInd; myCurEdgePar = Probing_Start; } @@ -268,7 +423,7 @@ bool BRepClass_FaceExplorer::OtherSegment(const gp_Pnt2d& P, gp_Lin2d& L, double } // if (anOrientation == TopAbs_FORWARD ... // This curve is not valid for line construction. Go to another edge. - myCurEdgeInd++; + ++myCurEdgeInd; myCurEdgePar = Probing_Start; } @@ -283,13 +438,15 @@ bool BRepClass_FaceExplorer::OtherSegment(const gp_Pnt2d& P, gp_Lin2d& L, double void BRepClass_FaceExplorer::InitWires() { - myWExplorer.Init(myFace, TopAbs_WIRE); + myCurrentWire = 0; } //================================================================================================= -bool BRepClass_FaceExplorer::RejectWire(const gp_Lin2d&, const double) const +bool BRepClass_FaceExplorer::RejectWire(const gp_Lin2d& theLine, const double theParameter) const { + (void)theLine; + (void)theParameter; return false; } @@ -297,26 +454,99 @@ bool BRepClass_FaceExplorer::RejectWire(const gp_Lin2d&, const double) const void BRepClass_FaceExplorer::InitEdges() { - myEExplorer.Init(myWExplorer.Current(), TopAbs_EDGE); - myMapVE.Clear(); - TopExp::MapShapesAndAncestors(myWExplorer.Current(), TopAbs_VERTEX, TopAbs_EDGE, myMapVE); + const WireData& aWire = myWires[myCurrentWire]; + myCurrentEdge = aWire.FirstEdge; + myCurrentEdgeEnd = aWire.FirstEdge + aWire.NbEdges; } //================================================================================================= -bool BRepClass_FaceExplorer::RejectEdge(const gp_Lin2d&, const double) const +bool BRepClass_FaceExplorer::RejectEdge(const gp_Lin2d& theLine, const double theParameter) const { + (void)theLine; + (void)theParameter; return false; } //================================================================================================= +void BRepClass_FaceExplorer::SetUseBndBox(const bool theValue) +{ + if (!theValue || myUseBndBox) + { + myUseBndBox = theValue; + return; + } + + for (const WireData& aWire : myWires) + { + for (uint32_t anEdgeIndex = aWire.FirstEdge; anEdgeIndex < aWire.FirstEdge + aWire.NbEdges; + ++anEdgeIndex) + { + BRepClass_Edge& anEdge = myEdges[anEdgeIndex]; + if (anEdge.BoundingBoxState() != BRepClass_Edge::BndBoxState::NotBuilt) + { + continue; + } + if (anEdge.Curve().IsNull()) + { + anEdge.SetBoundingBoxUnavailable(); + continue; + } + + try + { + OCC_CATCH_SIGNALS + Bnd_Box2d aBox; + BndLib_Add2dCurve::Add(anEdge.Curve(), + anEdge.FirstParameter(), + anEdge.LastParameter(), + 0.0, + aBox); + anEdge.SetBoundingBox(aBox); + } + catch (const Standard_Failure&) + { + anEdge.SetBoundingBoxUnavailable(); + } + } + } + myUseBndBox = true; +} + +//================================================================================================= + +bool BRepClass_FaceExplorer::ShouldUseBndBox() const +{ + if (myEdges.Size() <= THE_MIN_EDGES_FOR_BOUNDING_BOX) + { + return false; + } + + size_t aNbSplineEdges = 0; + for (const BRepClass_Edge& anEdge : myEdges) + { + if (anEdge.Curve().IsNull()) + { + continue; + } + const GeomAbs_CurveType aType = Geom2dAdaptor_Curve(anEdge.Curve()).GetType(); + if (aType == GeomAbs_BSplineCurve || aType == GeomAbs_BezierCurve + || aType == GeomAbs_OffsetCurve) + { + ++aNbSplineEdges; + } + } + return aNbSplineEdges * 2 >= myEdges.Size(); +} + +//================================================================================================= + void BRepClass_FaceExplorer::CurrentEdge(BRepClass_Edge& E, TopAbs_Orientation& Or) const { - E.Edge() = TopoDS::Edge(myEExplorer.Current()); - E.Face() = myFace; - Or = E.Edge().Orientation(); - E.SetNextEdge(myMapVE); + E = myEdges[myCurrentEdge]; + const BRepClass_Edge& anEdge = E; + Or = anEdge.Edge().Orientation(); E.SetMaxTolerance(myMaxTolerance); E.SetUseBndBox(myUseBndBox); } diff --git a/src/ModelingAlgorithms/TKTopAlgo/BRepClass/BRepClass_FaceExplorer.hxx b/src/ModelingAlgorithms/TKTopAlgo/BRepClass/BRepClass_FaceExplorer.hxx index 761f815e62..6ae68f1672 100644 --- a/src/ModelingAlgorithms/TKTopAlgo/BRepClass/BRepClass_FaceExplorer.hxx +++ b/src/ModelingAlgorithms/TKTopAlgo/BRepClass/BRepClass_FaceExplorer.hxx @@ -17,23 +17,21 @@ #ifndef _BRepClass_FaceExplorer_HeaderFile #define _BRepClass_FaceExplorer_HeaderFile -#include +#include +#include #include -#include -#include -#include -#include - #include #include -#include -#include + +#include + class gp_Pnt2d; class gp_Lin2d; -class BRepClass_Edge; //! Provide an exploration of a BRep Face for the //! classification. Return UV edges. +//! The explored topology and pcurves form a snapshot. Reconstruct the explorer after modifying the +//! underlying face or its edge representations. class BRepClass_FaceExplorer { public: @@ -65,10 +63,10 @@ public: Standard_EXPORT void InitWires(); //! Returns True if there is a current wire. - bool MoreWires() const { return myWExplorer.More(); } + bool MoreWires() const { return myCurrentWire < myWires.Size(); } //! Sets the explorer to the next wire. - void NextWire() { myWExplorer.Next(); } + void NextWire() { ++myCurrentWire; } //! Returns True if the wire bounding volume does not //! intersect the segment. @@ -79,10 +77,10 @@ public: Standard_EXPORT void InitEdges(); //! Returns True if there is a current edge. - bool MoreEdges() const { return myEExplorer.More(); } + bool MoreEdges() const { return myCurrentEdge < myCurrentEdgeEnd; } //! Sets the explorer to the next edge. - void NextEdge() { myEExplorer.Next(); } + void NextEdge() { ++myCurrentEdge; } //! Returns True if the edge bounding volume does not //! intersect the segment. @@ -104,22 +102,33 @@ public: //! Sets the status of whether we are //! using boxes or not - void SetUseBndBox(const bool theValue) { myUseBndBox = theValue; } + Standard_EXPORT void SetUseBndBox(const bool theValue); + + //! Returns true when cached boxes are expected to benefit this face. + Standard_EXPORT bool ShouldUseBndBox() const; protected: //! Computes UV bounds of a face Standard_EXPORT void ComputeFaceBounds(); private: - TopoDS_Face myFace; - TopExp_Explorer myWExplorer; - TopExp_Explorer myEExplorer; - int myCurEdgeInd; - double myCurEdgePar; - double myMaxTolerance; - bool myUseBndBox; - NCollection_IndexedDataMap, TopTools_ShapeMapHasher> - myMapVE; + struct WireData + { + uint32_t FirstEdge = 0; + uint32_t NbEdges = 0; + }; + + TopoDS_Face myFace; + NCollection_LinearVector myWires; + NCollection_LinearVector myEdges; + NCollection_LinearVector myProbeEdges; + size_t myCurrentWire; + size_t myCurrentEdge; + size_t myCurrentEdgeEnd; + size_t myCurEdgeInd; + double myCurEdgePar; + double myMaxTolerance; + bool myUseBndBox; double myUMin; double myUMax; diff --git a/src/ModelingAlgorithms/TKTopAlgo/BRepClass/BRepClass_Intersector.cxx b/src/ModelingAlgorithms/TKTopAlgo/BRepClass/BRepClass_Intersector.cxx index ac891ef5dd..253bcce3fe 100644 --- a/src/ModelingAlgorithms/TKTopAlgo/BRepClass/BRepClass_Intersector.cxx +++ b/src/ModelingAlgorithms/TKTopAlgo/BRepClass/BRepClass_Intersector.cxx @@ -338,7 +338,16 @@ void BRepClass_Intersector::Perform(const gp_Lin2d& L, const TopoDS_Face& F = E.Face(); // - aC2D = BRep_Tool::CurveOnSurface(EE, F, deb, fin); + aC2D = E.Curve(); + if (!aC2D.IsNull()) + { + deb = E.FirstParameter(); + fin = E.LastParameter(); + } + else + { + aC2D = BRep_Tool::CurveOnSurface(EE, F, deb, fin); + } if (aC2D.IsNull()) { done = false; // !IsDone() @@ -347,12 +356,23 @@ void BRepClass_Intersector::Perform(const gp_Lin2d& L, // Bnd_Box2d aBond; gp_Pnt2d aPntF; - bool anUseBndBox = E.UseBndBox(); - if (anUseBndBox) + bool anUseBndBox = false; + if (E.UseBndBox() && E.BoundingBoxState() != BRepClass_Edge::BndBoxState::Unavailable) { - BndLib_Add2dCurve::Add(aC2D, deb, fin, 0., aBond); - aBond.SetGap(aTolZ); - aPntF = L.Location(); + if (E.BoundingBoxState() == BRepClass_Edge::BndBoxState::Ready) + { + aBond = E.BoundingBox(); + } + else + { + BndLib_Add2dCurve::Add(aC2D, deb, fin, 0.0, aBond); + } + anUseBndBox = !aBond.IsVoid(); + if (anUseBndBox) + { + aBond.SetGap(aTolZ); + aPntF = L.Location(); + } } // Geom2dAdaptor_Curve C(aC2D, deb, fin); @@ -448,9 +468,14 @@ void BRepClass_Intersector::LocalGeometry(const BRepClass_Edge& E, gp_Dir2d& Norm, double& C) const { - double fpar, lpar; - occ::handle aPCurve = BRep_Tool::CurveOnSurface(E.Edge(), E.Face(), fpar, lpar); - GeomLProp_CLProps2d Prop(aPCurve, U, 2, Precision::PConfusion()); + double fpar = E.FirstParameter(); + double lpar = E.LastParameter(); + occ::handle aPCurve = E.Curve(); + if (aPCurve.IsNull()) + { + aPCurve = BRep_Tool::CurveOnSurface(E.Edge(), E.Face(), fpar, lpar); + } + GeomLProp_CLProps2d Prop(aPCurve, U, 2, Precision::PConfusion()); C = 0.; if (Prop.IsTangentDefined()) diff --git a/src/ModelingAlgorithms/TKTopAlgo/BRepTopAdaptor/BRepTopAdaptor_FClass2d.cxx b/src/ModelingAlgorithms/TKTopAlgo/BRepTopAdaptor/BRepTopAdaptor_FClass2d.cxx index e604cc312b..38e8dcf7f5 100644 --- a/src/ModelingAlgorithms/TKTopAlgo/BRepTopAdaptor/BRepTopAdaptor_FClass2d.cxx +++ b/src/ModelingAlgorithms/TKTopAlgo/BRepTopAdaptor/BRepTopAdaptor_FClass2d.cxx @@ -21,6 +21,7 @@ #include #include #include +#include #include #include #include @@ -31,8 +32,7 @@ #include #include #include -#include -#include +#include #include #include #include @@ -40,12 +40,22 @@ #include #include +#include +#include + #ifdef _MSC_VER #include #endif namespace { +//! Returns true for representable values, including OCCT's finite infinity sentinels. +inline bool isRepresentableValue(const double theValue) +{ + constexpr double THE_MAX_VALUE = std::numeric_limits::max(); + return theValue >= -THE_MAX_VALUE && theValue <= THE_MAX_VALUE; +} + // Increments @p theValue by @p theIncrement towards @p theDirection, ensuring that the result is // different from @p theValue. For large values of theValue with small theIncrement the result of // theValue + theIncrement can be equal to theValue due to the limited resolution of double @@ -83,11 +93,59 @@ bool isDegenerated(const BRepAdaptor_Curve& theCurve, return true; } + +//================================================================================================= + +struct PolygonMetrics +{ + double Area = 0.0; + double Perimeter = 0.0; +}; + +PolygonMetrics polygonMetrics(const NCollection_LinearVector& thePoints) +{ + PolygonMetrics aMetrics; + if (thePoints.Size() < 2) + { + return aMetrics; + } + + const size_t aLastUnique = thePoints.Size() - 2; + size_t aPrevious = aLastUnique; + for (size_t aCurrent = 0; aCurrent <= aLastUnique; ++aCurrent) + { + const gp_Pnt2d& aCurrentPoint = thePoints[aCurrent]; + const gp_Pnt2d& aPreviousPoint = thePoints[aPrevious]; + aMetrics.Area += + (aCurrentPoint.X() - aPreviousPoint.X()) * (aCurrentPoint.Y() + aPreviousPoint.Y()) * 0.5; + aMetrics.Perimeter += (aCurrentPoint.XY() - aPreviousPoint.XY()).Modulus(); + aPrevious = aCurrent; + } + return aMetrics; +} + +//================================================================================================= + +bool expectedThickness(const PolygonMetrics& theMetrics, double& theThickness) +{ + if (!isRepresentableValue(theMetrics.Area) || !isRepresentableValue(theMetrics.Perimeter) + || theMetrics.Perimeter <= 0.0) + { + return false; + } + + theThickness = std::max(2.0 * (std::abs(theMetrics.Area) / theMetrics.Perimeter), 1.e-7); + return isRepresentableValue(theThickness); +} } // namespace BRepTopAdaptor_FClass2d::BRepTopAdaptor_FClass2d(const TopoDS_Face& aFace, const double TolUV) : Toluv(TolUV), Face(aFace), + myIsUPeriodic(false), + myIsVPeriodic(false), + myUPeriod(0.0), + myVPeriod(0.0), U1(0.0), V1(0.0), U2(0.0), @@ -96,8 +154,11 @@ BRepTopAdaptor_FClass2d::BRepTopAdaptor_FClass2d(const TopoDS_Face& aFace, const //-- dead end on surfaces defined on more than one period Face.Orientation(TopAbs_FORWARD); - occ::handle surf = new BRepAdaptor_Surface(); - surf->Initialize(aFace, false); + BRepAdaptor_Surface aSurface(aFace, false); + myIsUPeriodic = aSurface.IsUPeriodic(); + myIsVPeriodic = aSurface.IsVPeriodic(); + myUPeriod = myIsUPeriodic ? aSurface.UPeriod() : 0.0; + myVPeriod = myIsVPeriodic ? aSurface.VPeriod() : 0.0; TopoDS_Edge edge; TopAbs_Orientation Or; @@ -112,13 +173,13 @@ BRepTopAdaptor_FClass2d::BRepTopAdaptor_FClass2d(const TopoDS_Face& aFace, const for (TopExp_Explorer aFaceExplorer(Face, TopAbs_WIRE); (aFaceExplorer.More() && !anIsBadWire); aFaceExplorer.Next()) { - int nbpnts = 0; - NCollection_Sequence SeqPnt2d; - int firstpoint = 1; - double FlecheU = 0.0; - double FlecheV = 0.0; - bool WireIsNotEmpty = false; - int NbEdges = 0; + int nbpnts = 0; + NCollection_LinearVector SeqPnt2d; + int firstpoint = 1; + double FlecheU = 0.0; + double FlecheV = 0.0; + bool WireIsNotEmpty = false; + int NbEdges = 0; TopExp_Explorer Explorer; for (Explorer.Init(aFaceExplorer.Current(), TopAbs_EDGE); Explorer.More(); Explorer.Next()) @@ -261,14 +322,14 @@ BRepTopAdaptor_FClass2d::BRepTopAdaptor_FClass2d(const TopoDS_Face& aFace, const //",nbpnts,u,FlecheU,FlecheV,ii,Avant); // if(ii>(Avant+4)) // Modified by Sergey KHROMOV - Fri Apr 19 09:46:12 2002 Begin - if (ii > (Avant + 4) && SeqPnt2d(ii - 2).SquareDistance(SeqPnt2d(ii))) + if (ii > (Avant + 4) && SeqPnt2d[ii - 3].SquareDistance(SeqPnt2d[ii - 1])) // Modified by Sergey KHROMOV - Fri Apr 19 09:46:13 2002 End { - gp_Lin2d Lin(SeqPnt2d(ii - 2), gp_Dir2d(gp_Vec2d(SeqPnt2d(ii - 2), SeqPnt2d(ii)))); - double ul = ElCLib::Parameter(Lin, SeqPnt2d(ii - 1)); + gp_Lin2d Lin(SeqPnt2d[ii - 3], gp_Dir2d(gp_Vec2d(SeqPnt2d[ii - 3], SeqPnt2d[ii - 1]))); + double ul = ElCLib::Parameter(Lin, SeqPnt2d[ii - 2]); gp_Pnt2d Pp = ElCLib::Value(ul, Lin); - double dU = std::abs(Pp.X() - SeqPnt2d(ii - 1).X()); - double dV = std::abs(Pp.Y() - SeqPnt2d(ii - 1).Y()); + double dU = std::abs(Pp.X() - SeqPnt2d[ii - 2].X()); + double dV = std::abs(Pp.Y() - SeqPnt2d[ii - 2].Y()); //-- printf(" (du=%7.5g dv=%7.5g)",dU,dV); if (dU > FlecheU) { @@ -290,19 +351,15 @@ BRepTopAdaptor_FClass2d::BRepTopAdaptor_FClass2d(const TopoDS_Face& aFace, const if (NbEdges) { //-- on compte ++ with a normal explorer and with the Wire Explorer - NCollection_Array1 PClass(1, 2); - //// modified by jgv, 28.04.2009 //// - PClass.Init(gp_Pnt2d(0., 0.)); - ///////////////////////////////////// - TabClass.Append(CSLib_Class2d(PClass, FlecheU, FlecheV, Umin, Vmin, Umax, Vmax)); + NCollection_LinearVector aPoints(2, gp_Pnt2d(0.0, 0.0)); + TabClass.Append(CSLib_Class2d(aPoints.ToArray1(), FlecheU, FlecheV, Umin, Vmin, Umax, Vmax)); anIsBadWire = true; - TabOrien.Append(-1); + TabOrien.Append(WireRole::Invalid); } else if (WireIsNotEmpty) { // double anglep=0,anglem=0; - NCollection_Array1 PClass(1, nbpnts); - double square = 0.0; + double aSignedArea = 0.0; //------------------------------------------------------------------- //-- ** The mode of calculation was somewhat changed @@ -314,45 +371,30 @@ BRepTopAdaptor_FClass2d::BRepTopAdaptor_FClass2d(const TopoDS_Face& aFace, const if (nbpnts > 3) { - // int im2=nbpnts-2; - int im1 = nbpnts - 1; - int im0 = 1; - // PClass(im2)=SeqPnt2d.Value(im2); - PClass(im1) = SeqPnt2d.Value(im1); - PClass(nbpnts) = SeqPnt2d.Value(nbpnts); - - double aPer = 0.; - // for(int ii=1; ii=nbpnts) im2=1; - if (im1 >= nbpnts) - { - im1 = 1; - } - PClass(ii) = SeqPnt2d.Value(ii); - // gp_Vec2d A(PClass(im2),PClass(im1)); - // gp_Vec2d B(PClass(im1),PClass(im0)); - // double N = A.Magnitude() * B.Magnitude(); - - square += (PClass(im0).X() - PClass(im1).X()) * (PClass(im0).Y() + PClass(im1).Y()) * .5; - aPer += (PClass(im0).XY() - PClass(im1).XY()).Modulus(); - - // if(N>1e-16){ double a=A.Angle(B); angle+=a; } + anIsBadWire = true; + TabOrien.Append(WireRole::Invalid); + NCollection_LinearVector aFallbackPoints(2, gp_Pnt2d(0.0, 0.0)); + TabClass.Append( + CSLib_Class2d(aFallbackPoints.ToArray1(), FlecheU, FlecheV, Umin, Vmin, Umax, Vmax)); + continue; } - - double anExpThick = std::max(2. * std::abs(square) / aPer, 1e-7); double aDefl = std::max(FlecheU, FlecheV); double aDiscrDefl = std::min(aDefl * 0.1, anExpThick * 10.); while (aDefl > anExpThick && aDiscrDefl > 1e-7) { // Deflection of the polygon is too much for this ratio of area and perimeter, // and this might lead to self-intersections. - // Discretize the wire more tightly to eliminate the error. - firstpoint = 1; - SeqPnt2d.Clear(); - FlecheU = 0.0; - FlecheV = 0.0; + // Build tighter samples separately so a failed edge leaves the coarse polygon intact. + NCollection_LinearVector aRefinedPoints; + int aRefinedFirstPoint = 1; + double aRefinedFlecheU = 0.0; + double aRefinedFlecheV = 0.0; + bool isRefinementDone = true; for (WireExplorer.Init(TopoDS::Wire(aFaceExplorer.Current()), Face); WireExplorer.More(); WireExplorer.Next()) { @@ -370,6 +412,7 @@ BRepTopAdaptor_FClass2d::BRepTopAdaptor_FClass2d(const TopoDS_Face& aFace, const GCPnts_QuasiUniformDeflection aDiscr(C, aDiscrDefl); if (!aDiscr.IsDone()) { + isRefinementDone = false; break; } int nbp = aDiscr.NbPoints(); @@ -380,69 +423,64 @@ BRepTopAdaptor_FClass2d::BRepTopAdaptor_FClass2d(const TopoDS_Face& aFace, const i = nbp; iEnd = 0; } - if (firstpoint == 2) + if (aRefinedFirstPoint == 2) { i += iStep; } for (; i != iEnd; i += iStep) { - gp_Pnt2d aP2d = C.Value(aDiscr.Parameter(i)); - SeqPnt2d.Append(aP2d); + aRefinedPoints.Append(C.Value(aDiscr.Parameter(i))); } if (nbp > 2) { - int ii = SeqPnt2d.Length(); - gp_Lin2d Lin(SeqPnt2d(ii - 2), gp_Dir2d(gp_Vec2d(SeqPnt2d(ii - 2), SeqPnt2d(ii)))); - double ul = ElCLib::Parameter(Lin, SeqPnt2d(ii - 1)); - gp_Pnt2d Pp = ElCLib::Value(ul, Lin); - double dU = std::abs(Pp.X() - SeqPnt2d(ii - 1).X()); - double dV = std::abs(Pp.Y() - SeqPnt2d(ii - 1).Y()); - if (dU > FlecheU) + const size_t ii = aRefinedPoints.Size(); + gp_Lin2d Lin(aRefinedPoints[ii - 3], + gp_Dir2d(gp_Vec2d(aRefinedPoints[ii - 3], aRefinedPoints[ii - 1]))); + double ul = ElCLib::Parameter(Lin, aRefinedPoints[ii - 2]); + gp_Pnt2d Pp = ElCLib::Value(ul, Lin); + double dU = std::abs(Pp.X() - aRefinedPoints[ii - 2].X()); + double dV = std::abs(Pp.Y() - aRefinedPoints[ii - 2].Y()); + if (dU > aRefinedFlecheU) { - FlecheU = dU; + aRefinedFlecheU = dU; } - if (dV > FlecheV) + if (dV > aRefinedFlecheV) { - FlecheV = dV; + aRefinedFlecheV = dV; } } - firstpoint = 2; + aRefinedFirstPoint = 2; } } - nbpnts = SeqPnt2d.Length(); - PClass.Resize(1, nbpnts, false); - im1 = nbpnts - 1; - im0 = 1; - PClass(im1) = SeqPnt2d.Value(im1); - PClass(nbpnts) = SeqPnt2d.Value(nbpnts); - square = 0.; - aPer = 0.; - for (int ii = 1; ii < nbpnts; ii++, im0++, im1++) - { - if (im1 >= nbpnts) - { - im1 = 1; - } - PClass(ii) = SeqPnt2d.Value(ii); - square += - (PClass(im0).X() - PClass(im1).X()) * (PClass(im0).Y() + PClass(im1).Y()) * .5; - aPer += (PClass(im0).XY() - PClass(im1).XY()).Modulus(); - } - anExpThick = std::max(2. * std::abs(square) / aPer, 1e-7); - aDefl = std::max(FlecheU, FlecheV); - aDiscrDefl = std::min(aDiscrDefl * 0.1, anExpThick * 10.); + const PolygonMetrics aRefinedMetrics = polygonMetrics(aRefinedPoints); + double aRefinedThickness = 0.0; + if (!isRefinementDone || aRefinedPoints.Size() <= 3 + || !expectedThickness(aRefinedMetrics, aRefinedThickness)) + { + break; + } + + SeqPnt2d = std::move(aRefinedPoints); + nbpnts = static_cast(SeqPnt2d.Size()); + FlecheU = aRefinedFlecheU; + FlecheV = aRefinedFlecheV; + aMetrics = aRefinedMetrics; + aSignedArea = aMetrics.Area; + anExpThick = aRefinedThickness; + aDefl = std::max(FlecheU, FlecheV); + aDiscrDefl = std::min(aDiscrDefl * 0.1, anExpThick * 10.); } //-- FlecheU*=10.0; //-- FlecheV*=10.0; - if (aNbE == 1 && FlecheU < eps && FlecheV < eps && std::abs(square) < eps) + if (aNbE == 1 && FlecheU < eps && FlecheV < eps && std::abs(aSignedArea) < eps) { - TabOrien.Append(1); + TabOrien.Append(WireRole::Outer); } else { - TabOrien.Append(((square < 0.0) ? 1 : 0)); + TabOrien.Append(aSignedArea < 0.0 ? WireRole::Outer : WireRole::Inner); } if (FlecheU < Toluv) @@ -453,33 +491,32 @@ BRepTopAdaptor_FClass2d::BRepTopAdaptor_FClass2d(const TopoDS_Face& aFace, const { FlecheV = Toluv; } - TabClass.Append(CSLib_Class2d(PClass, FlecheU, FlecheV, Umin, Vmin, Umax, Vmax)); + TabClass.Append( + CSLib_Class2d(SeqPnt2d.ToArray1(), FlecheU, FlecheV, Umin, Vmin, Umax, Vmax)); } // if(nbpoints>3 else { anIsBadWire = true; - TabOrien.Append(-1); - NCollection_Array1 xPClass(1, 2); - xPClass(1) = SeqPnt2d(1); - xPClass(2) = SeqPnt2d(2); - TabClass.Append(CSLib_Class2d(xPClass, FlecheU, FlecheV, Umin, Vmin, Umax, Vmax)); + TabOrien.Append(WireRole::Invalid); + TabClass.Append( + CSLib_Class2d(SeqPnt2d.ToArray1(), FlecheU, FlecheV, Umin, Vmin, Umax, Vmax)); } } // else if(WareIsNotEmpty } // for(FaceExplorer - int nbtabclass = TabClass.Length(); + const size_t nbtabclass = TabClass.Size(); if (nbtabclass > 0) { //-- If an error was detected on a wire: set all TabOrien to -1 if (anIsBadWire) { - TabOrien(1) = -1; + TabOrien[0] = WireRole::Invalid; } - if (surf->GetType() == GeomAbs_Cone || surf->GetType() == GeomAbs_Cylinder - || surf->GetType() == GeomAbs_Torus || surf->GetType() == GeomAbs_Sphere - || surf->GetType() == GeomAbs_SurfaceOfRevolution) + if (aSurface.GetType() == GeomAbs_Cone || aSurface.GetType() == GeomAbs_Cylinder + || aSurface.GetType() == GeomAbs_Torus || aSurface.GetType() == GeomAbs_Sphere + || aSurface.GetType() == GeomAbs_SurfaceOfRevolution) { double uuu = M_PI + M_PI - (Umax - Umin); @@ -495,7 +532,7 @@ BRepTopAdaptor_FClass2d::BRepTopAdaptor_FClass2d(const TopoDS_Face& aFace, const U1 = U2 = 0.0; } - if (surf->GetType() == GeomAbs_Torus) + if (aSurface.GetType() == GeomAbs_Torus) { double uuu = M_PI + M_PI - (Vmax - Vmin); if (uuu < 0) @@ -512,6 +549,15 @@ BRepTopAdaptor_FClass2d::BRepTopAdaptor_FClass2d(const TopoDS_Face& aFace, const } } +//================================================================================================= + +BRepTopAdaptor_FClass2d::~BRepTopAdaptor_FClass2d() +{ + Destroy(); +} + +//================================================================================================= + TopAbs_State BRepTopAdaptor_FClass2d::PerformInfinitePoint() const { if (Umax == -RealLast() || Vmax == -RealLast() || Umin == RealLast() || Vmin == RealLast()) @@ -522,338 +568,196 @@ TopAbs_State BRepTopAdaptor_FClass2d::PerformInfinitePoint() const return (Perform(P, false)); } -TopAbs_State BRepTopAdaptor_FClass2d::Perform(const gp_Pnt2d& _Puv, - const bool RecadreOnPeriodic) const -{ - int dedans; - int nbtabclass = TabClass.Length(); +//================================================================================================= - if (nbtabclass == 0) +TopAbs_State BRepTopAdaptor_FClass2d::exactState(const gp_Pnt2d& thePoint, + const double theTolerance) const +{ + std::lock_guard aLock(myExactMutex); + if (!myExactExplorer) { - return (TopAbs_IN); + myExactExplorer = std::make_unique(Face); + myExactExplorer->SetUseBndBox(myExactExplorer->ShouldUseBndBox()); + } + BRepClass_FaceClassifier aClassifier(*myExactExplorer, thePoint, theTolerance); + return aClassifier.State(); +} + +//================================================================================================= + +TopAbs_State BRepTopAdaptor_FClass2d::classify(const gp_Pnt2d& thePoint, + const double theTolerance, + const bool theRecadreOnPeriodic, + const ClassificationMode theMode) const +{ + const size_t aClassifierCount = TabClass.Size(); + if (aClassifierCount == 0) + { + return TopAbs_IN; } //-- U1 is the First Param and U2 in this case is U1+Period - double u = _Puv.X(); - double v = _Puv.Y(); - double uu = u, vv = v; + double aU = thePoint.X(); + double aV = thePoint.Y(); + double aReframedU = aU; + double aReframedV = aV; + bool isUReframed = false; + bool isVReframed = false; - occ::handle surf = new BRepAdaptor_Surface(); - surf->Initialize(Face, false); - const bool IsUPer = surf->IsUPeriodic(); - const bool IsVPer = surf->IsVPeriodic(); - const double uperiod = IsUPer ? surf->UPeriod() : 0.0; - const double vperiod = IsVPer ? surf->VPeriod() : 0.0; - TopAbs_State aStatus = TopAbs_UNKNOWN; - bool urecadre = false, vrecadre = false; - - if (RecadreOnPeriodic) + if (theRecadreOnPeriodic) { - if (IsUPer) + if (myIsUPeriodic) { - if (uu < Umin) + if (aReframedU < Umin) { - while (uu < Umin) + while (aReframedU < Umin) { - uu += uperiod; + aReframedU += myUPeriod; } } else { - while (uu >= Umin) + while (aReframedU >= Umin) { - uu -= uperiod; + aReframedU -= myUPeriod; } - uu += uperiod; + aReframedU += myUPeriod; } } - if (IsVPer) + if (myIsVPeriodic) { - if (vv < Vmin) + if (aReframedV < Vmin) { - while (vv < Vmin) + while (aReframedV < Vmin) { - vv += vperiod; + aReframedV += myVPeriod; } } else { - while (vv >= Vmin) + while (aReframedV >= Vmin) { - vv -= vperiod; + aReframedV -= myVPeriod; } - vv += vperiod; + aReframedV += myVPeriod; } } } + TopAbs_State aState = TopAbs_UNKNOWN; for (;;) { - dedans = 1; - gp_Pnt2d Puv(u, v); + const gp_Pnt2d aPoint(aU, aV); - if (TabOrien(1) != -1) + if (TabOrien[0] != WireRole::Invalid) { - for (int n = 1; n <= nbtabclass; n++) + CSLib_Class2d::Result aPolygonResult = CSLib_Class2d::Result_Inside; + for (size_t aWireIndex = 0; aWireIndex < aClassifierCount; ++aWireIndex) { - int cur = TabClass(n).SiDans(Puv); - if (cur == 1) + const int aWireResult = theMode == ClassificationMode::Perform + ? TabClass[aWireIndex].SiDans(aPoint) + : TabClass[aWireIndex].SiDans_OnMode(aPoint, theTolerance); + if (aWireResult == CSLib_Class2d::Result_Inside) { - if (TabOrien(n) == 0) + if (TabOrien[aWireIndex] == WireRole::Inner) { - dedans = -1; + aPolygonResult = CSLib_Class2d::Result_Outside; break; } } - else if (cur == -1) + else if (aWireResult == CSLib_Class2d::Result_Outside) { - if (TabOrien(n) == 1) + if (TabOrien[aWireIndex] == WireRole::Outer) { - dedans = -1; + aPolygonResult = CSLib_Class2d::Result_Outside; break; } } else { - dedans = 0; + aPolygonResult = CSLib_Class2d::Result_Uncertain; break; } } - if (dedans == 0) + + if (aPolygonResult == CSLib_Class2d::Result_Uncertain) { - BRepClass_FaceClassifier aClassifier; - double m_Toluv = (Toluv > 4.0) ? 4.0 : Toluv; - // aClassifier.Perform(Face,Puv,Toluv); - aClassifier.Perform(Face, Puv, m_Toluv); - aStatus = aClassifier.State(); + aState = theMode == ClassificationMode::Perform ? exactState(aPoint, std::min(Toluv, 4.0)) + : TopAbs_ON; } - if (dedans == 1) + else { - aStatus = TopAbs_IN; - } - if (dedans == -1) - { - aStatus = TopAbs_OUT; + aState = aPolygonResult == CSLib_Class2d::Result_Inside ? TopAbs_IN : TopAbs_OUT; } } else - { //-- TabOrien(1)=-1 False Wire - BRepClass_FaceClassifier aClassifier; - aClassifier.Perform(Face, Puv, Toluv); - aStatus = aClassifier.State(); + { + // A malformed wire cannot be classified reliably by its polygon. + aState = exactState(aPoint, theMode == ClassificationMode::Perform ? Toluv : theTolerance); } - if (!RecadreOnPeriodic || (!IsUPer && !IsVPer)) + if (!theRecadreOnPeriodic || (!myIsUPeriodic && !myIsVPeriodic)) { - return aStatus; + return aState; } - if (aStatus == TopAbs_IN || aStatus == TopAbs_ON) + if (aState == TopAbs_IN || aState == TopAbs_ON) { - return aStatus; + return aState; } - if (!urecadre) + if (!isUReframed) { - u = uu; - urecadre = true; + aU = aReframedU; + isUReframed = true; } - else if (IsUPer) + else if (myIsUPeriodic) { - u += uperiod; + aU += myUPeriod; } - if (u > Umax || !IsUPer) + if (aU > Umax || !myIsUPeriodic) { - if (!vrecadre) + if (!isVReframed) { - v = vv; - vrecadre = true; + aV = aReframedV; + isVReframed = true; } - else if (IsVPer) + else if (myIsVPeriodic) { - v += vperiod; + aV += myVPeriod; } - u = uu; + aU = aReframedU; - if (v > Vmax || !IsVPer) + if (aV > Vmax || !myIsVPeriodic) { - return aStatus; + return aState; } } - } // for (;;) + } } -TopAbs_State BRepTopAdaptor_FClass2d::TestOnRestriction(const gp_Pnt2d& _Puv, - const double Tol, - const bool RecadreOnPeriodic) const +//================================================================================================= + +TopAbs_State BRepTopAdaptor_FClass2d::Perform(const gp_Pnt2d& thePoint, + const bool theRecadreOnPeriodic) const { - int dedans; - int nbtabclass = TabClass.Length(); - - if (nbtabclass == 0) - { - return (TopAbs_IN); - } - - //-- U1 is the First Param and U2 in this case is U1+Period - double u = _Puv.X(); - double v = _Puv.Y(); - double uu = u, vv = v; - - occ::handle surf = new BRepAdaptor_Surface(); - surf->Initialize(Face, false); - const bool IsUPer = surf->IsUPeriodic(); - const bool IsVPer = surf->IsVPeriodic(); - const double uperiod = IsUPer ? surf->UPeriod() : 0.0; - const double vperiod = IsVPer ? surf->VPeriod() : 0.0; - TopAbs_State aStatus = TopAbs_UNKNOWN; - bool urecadre = false, vrecadre = false; - - if (RecadreOnPeriodic) - { - if (IsUPer) - { - if (uu < Umin) - { - while (uu < Umin) - { - uu += uperiod; - } - } - else - { - while (uu >= Umin) - { - uu -= uperiod; - } - uu += uperiod; - } - } - if (IsVPer) - { - if (vv < Vmin) - { - while (vv < Vmin) - { - vv += vperiod; - } - } - else - { - while (vv >= Vmin) - { - vv -= vperiod; - } - vv += vperiod; - } - } - } - - for (;;) - { - dedans = 1; - gp_Pnt2d Puv(u, v); - - if (TabOrien(1) != -1) - { - for (int n = 1; n <= nbtabclass; n++) - { - int cur = TabClass(n).SiDans_OnMode(Puv, Tol); - if (cur == 1) - { - if (TabOrien(n) == 0) - { - dedans = -1; - break; - } - } - else if (cur == -1) - { - if (TabOrien(n) == 1) - { - dedans = -1; - break; - } - } - else - { - dedans = 0; - break; - } - } - if (dedans == 0) - { - aStatus = TopAbs_ON; - } - if (dedans == 1) - { - aStatus = TopAbs_IN; - } - if (dedans == -1) - { - aStatus = TopAbs_OUT; - } - } - else - { //-- TabOrien(1)=-1 False Wire - BRepClass_FaceClassifier aClassifier; - aClassifier.Perform(Face, Puv, Tol); - aStatus = aClassifier.State(); - } - - if (!RecadreOnPeriodic || (!IsUPer && !IsVPer)) - { - return aStatus; - } - if (aStatus == TopAbs_IN || aStatus == TopAbs_ON) - { - return aStatus; - } - - if (!urecadre) - { - u = uu; - urecadre = true; - } - else if (IsUPer) - { - u += uperiod; - } - if (u > Umax || !IsUPer) - { - if (!vrecadre) - { - v = vv; - vrecadre = true; - } - else if (IsVPer) - { - v += vperiod; - } - - u = uu; - - if (v > Vmax || !IsVPer) - { - return aStatus; - } - } - } // for (;;) + return classify(thePoint, Toluv, theRecadreOnPeriodic, ClassificationMode::Perform); } +//================================================================================================= + +TopAbs_State BRepTopAdaptor_FClass2d::TestOnRestriction(const gp_Pnt2d& thePoint, + const double theTolerance, + const bool theRecadreOnPeriodic) const +{ + return classify(thePoint, theTolerance, theRecadreOnPeriodic, ClassificationMode::OnRestriction); +} + +//================================================================================================= + void BRepTopAdaptor_FClass2d::Destroy() { - TabClass.Clear(); -} - -#include - -// const BRepTopAdaptor_FClass2d & BRepTopAdaptor_FClass2d::Copy(const BRepTopAdaptor_FClass2d& -// Other) const { -const BRepTopAdaptor_FClass2d& BRepTopAdaptor_FClass2d::Copy(const BRepTopAdaptor_FClass2d&) const -{ -#ifdef OCCT_DEBUG - std::cerr << "Copy not allowed in BRepTopAdaptor_FClass2d" << std::endl; -#endif - throw Standard_ConstructionError(); + TabClass.Clear(true); + TabOrien.Clear(true); + myExactExplorer.reset(); } diff --git a/src/ModelingAlgorithms/TKTopAlgo/BRepTopAdaptor/BRepTopAdaptor_FClass2d.hxx b/src/ModelingAlgorithms/TKTopAlgo/BRepTopAdaptor/BRepTopAdaptor_FClass2d.hxx index 1496b6e65e..253745edcd 100644 --- a/src/ModelingAlgorithms/TKTopAlgo/BRepTopAdaptor/BRepTopAdaptor_FClass2d.hxx +++ b/src/ModelingAlgorithms/TKTopAlgo/BRepTopAdaptor/BRepTopAdaptor_FClass2d.hxx @@ -21,11 +21,16 @@ #include #include -#include +#include #include #include #include +#include +#include +#include + +class BRepClass_FaceExplorer; class gp_Pnt2d; class BRepTopAdaptor_FClass2d @@ -35,6 +40,11 @@ public: Standard_EXPORT BRepTopAdaptor_FClass2d(const TopoDS_Face& F, const double Tol); + BRepTopAdaptor_FClass2d(const BRepTopAdaptor_FClass2d&) = delete; + BRepTopAdaptor_FClass2d(BRepTopAdaptor_FClass2d&&) = delete; + BRepTopAdaptor_FClass2d& operator=(const BRepTopAdaptor_FClass2d&) = delete; + BRepTopAdaptor_FClass2d& operator=(BRepTopAdaptor_FClass2d&&) = delete; + Standard_EXPORT TopAbs_State PerformInfinitePoint() const; Standard_EXPORT TopAbs_State Perform(const gp_Pnt2d& Puv, @@ -42,14 +52,7 @@ public: Standard_EXPORT void Destroy(); - ~BRepTopAdaptor_FClass2d() { Destroy(); } - - Standard_EXPORT const BRepTopAdaptor_FClass2d& Copy(const BRepTopAdaptor_FClass2d& Other) const; - - const BRepTopAdaptor_FClass2d& operator=(const BRepTopAdaptor_FClass2d& Other) const - { - return Copy(Other); - } + Standard_EXPORT ~BRepTopAdaptor_FClass2d(); //! Test a point with +- an offset (Tol) and returns //! On if some points are OUT an some are IN @@ -59,18 +62,45 @@ public: const bool RecadreOnPeriodic = true) const; private: - NCollection_Sequence TabClass; - NCollection_Sequence TabOrien; - double Toluv; - TopoDS_Face Face; - double U1; - double V1; - double U2; - double V2; - double Umin; - double Umax; - double Vmin; - double Vmax; + enum class ClassificationMode : uint8_t + { + Perform, + OnRestriction + }; + + enum class WireRole : int8_t + { + Invalid = -1, + Inner = 0, + Outer = 1 + }; + + TopAbs_State exactState(const gp_Pnt2d& thePoint, const double theTolerance) const; + + TopAbs_State classify(const gp_Pnt2d& thePoint, + const double theTolerance, + const bool theRecadreOnPeriodic, + const ClassificationMode theMode) const; + +private: + NCollection_LinearVector TabClass; + NCollection_LinearVector TabOrien; + double Toluv; + TopoDS_Face Face; + mutable std::unique_ptr myExactExplorer; + mutable std::mutex myExactMutex; + bool myIsUPeriodic; + bool myIsVPeriodic; + double myUPeriod; + double myVPeriod; + double U1; + double V1; + double U2; + double V2; + double Umin; + double Umax; + double Vmin; + double Vmax; }; #endif // _BRepTopAdaptor_FClass2d_HeaderFile diff --git a/src/ModelingAlgorithms/TKTopAlgo/GTests/BRepClass_FaceExplorer_Test.cxx b/src/ModelingAlgorithms/TKTopAlgo/GTests/BRepClass_FaceExplorer_Test.cxx new file mode 100644 index 0000000000..3b41b734b4 --- /dev/null +++ b/src/ModelingAlgorithms/TKTopAlgo/GTests/BRepClass_FaceExplorer_Test.cxx @@ -0,0 +1,386 @@ +// Copyright (c) 2026 OPEN CASCADE SAS +// +// This file is part of Open CASCADE Technology software library. +// +// This library is free software; you can redistribute it and/or modify it under +// the terms of the GNU Lesser General Public License version 2.1 as published +// by the Free Software Foundation, with special exception defined in the file +// OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT +// distribution for complete text of the license and disclaimer of any warranty. +// +// Alternatively, this file may be used under the terms of Open CASCADE +// commercial license or contractual agreement. + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +namespace +{ + +TopoDS_Wire makeRectangle(const double theXMin, + const double theYMin, + const double theXMax, + const double theYMax) +{ + const gp_Pnt aPoints[] = {gp_Pnt(theXMin, theYMin, 0.0), + gp_Pnt(theXMax, theYMin, 0.0), + gp_Pnt(theXMax, theYMax, 0.0), + gp_Pnt(theXMin, theYMax, 0.0)}; + BRepBuilderAPI_MakeWire aWireBuilder; + for (int anIndex = 0; anIndex < 4; ++anIndex) + { + BRepBuilderAPI_MakeEdge anEdgeBuilder(aPoints[anIndex], aPoints[(anIndex + 1) % 4]); + if (!anEdgeBuilder.IsDone()) + { + ADD_FAILURE() << "Edge creation failed"; + return TopoDS_Wire(); + } + aWireBuilder.Add(anEdgeBuilder.Edge()); + } + if (!aWireBuilder.IsDone()) + { + ADD_FAILURE() << "Wire creation failed"; + return TopoDS_Wire(); + } + return aWireBuilder.Wire(); +} + +TopAbs_State classify(const TopoDS_Face& theFace, const gp_Pnt2d& thePoint, const bool theUseBndBox) +{ + BRepClass_FaceExplorer anExplorer(theFace); + anExplorer.SetUseBndBox(theUseBndBox); + BRepClass_FClassifier aClassifier(anExplorer, thePoint, Precision::Confusion()); + return aClassifier.State(); +} + +TopoDS_Edge referenceNextEdge(const TopoDS_Edge& theEdge, const TopoDS_Wire& theWire) +{ + TopoDS_Vertex aFirstVertex; + TopoDS_Vertex aLastVertex; + TopExp::Vertices(theEdge, aFirstVertex, aLastVertex, true); + if (aLastVertex.IsNull() || aLastVertex.IsSame(aFirstVertex)) + { + return TopoDS_Edge(); + } + + NCollection_IndexedDataMap, TopTools_ShapeMapHasher> + aVertexEdges; + TopExp::MapShapesAndAncestors(theWire, TopAbs_VERTEX, TopAbs_EDGE, aVertexEdges); + const NCollection_List* anEdges = aVertexEdges.Seek(aLastVertex); + if (anEdges == nullptr || anEdges->Extent() != 2) + { + return TopoDS_Edge(); + } + + TopoDS_Edge aNextEdge; + for (NCollection_List::Iterator anIt(*anEdges); anIt.More(); anIt.Next()) + { + if (!anIt.Value().IsNull() && !anIt.Value().IsSame(theEdge)) + { + aNextEdge = TopoDS::Edge(anIt.Value()); + } + } + return aNextEdge; +} + +void expectNextEdgesMatchReference(const TopoDS_Face& theFace) +{ + BRepClass_FaceExplorer anExplorer(theFace); + TopExp_Explorer aWireExp(theFace, TopAbs_WIRE); + for (anExplorer.InitWires(); anExplorer.MoreWires(); anExplorer.NextWire(), aWireExp.Next()) + { + ASSERT_TRUE(aWireExp.More()); + const TopoDS_Wire& aWire = TopoDS::Wire(aWireExp.Current()); + TopExp_Explorer anEdgeExp(aWire, TopAbs_EDGE); + for (anExplorer.InitEdges(); anExplorer.MoreEdges(); anExplorer.NextEdge(), anEdgeExp.Next()) + { + ASSERT_TRUE(anEdgeExp.More()); + BRepClass_Edge anEdgeData; + TopAbs_Orientation anOrientation = TopAbs_EXTERNAL; + anExplorer.CurrentEdge(anEdgeData, anOrientation); + const TopoDS_Edge& anExpectedEdge = TopoDS::Edge(anEdgeExp.Current()); + ASSERT_TRUE(anEdgeData.Edge().IsEqual(anExpectedEdge)); + + const TopoDS_Edge aExpectedNext = referenceNextEdge(anExpectedEdge, aWire); + if (aExpectedNext.IsNull()) + { + EXPECT_TRUE(anEdgeData.NextEdge().IsNull()); + } + else + { + EXPECT_TRUE(anEdgeData.NextEdge().IsEqual(aExpectedNext)); + } + } + EXPECT_FALSE(anEdgeExp.More()); + } + EXPECT_FALSE(aWireExp.More()); +} + +} // namespace + +TEST(BRepClassFaceExplorerTest, MultiWireHole_BoundingBoxesPreserveClassification) +{ + const TopoDS_Wire anInnerWire = makeRectangle(3.0, 3.0, 7.0, 7.0); + ASSERT_FALSE(anInnerWire.IsNull()); + + BRepBuilderAPI_MakeFace aFaceBuilder(gp_Pln(gp::XOY()), 0.0, 10.0, 0.0, 10.0); + ASSERT_TRUE(aFaceBuilder.IsDone()); + aFaceBuilder.Add(TopoDS::Wire(anInnerWire.Reversed())); + ASSERT_TRUE(aFaceBuilder.IsDone()); + const TopoDS_Face aFace = aFaceBuilder.Face(); + + const gp_Pnt2d aPoints[] = {gp_Pnt2d(1.0, 1.0), gp_Pnt2d(5.0, 5.0), gp_Pnt2d(12.0, 5.0)}; + const TopAbs_State anExpected[] = {TopAbs_IN, TopAbs_OUT, TopAbs_OUT}; + for (int anIndex = 0; anIndex < 3; ++anIndex) + { + EXPECT_EQ(classify(aFace, aPoints[anIndex], false), anExpected[anIndex]); + EXPECT_EQ(classify(aFace, aPoints[anIndex], true), anExpected[anIndex]); + } +} + +TEST(BRepClassFaceExplorerTest, CylindricalSeam_PreservesWireOccurrenceOrder) +{ + occ::handle aSurface = + new Geom_CylindricalSurface(gp_Ax3(gp::Origin(), gp::DZ()), 2.0); + BRepBuilderAPI_MakeFace aFaceBuilder(aSurface, 0.0, 2.0 * M_PI, 0.0, 5.0, Precision::Confusion()); + ASSERT_TRUE(aFaceBuilder.IsDone()); + const TopoDS_Face aFace = aFaceBuilder.Face(); + + NCollection_Sequence anExpectedEdges; + for (TopExp_Explorer aWireExp(aFace, TopAbs_WIRE); aWireExp.More(); aWireExp.Next()) + { + for (TopExp_Explorer anEdgeExp(aWireExp.Current(), TopAbs_EDGE); anEdgeExp.More(); + anEdgeExp.Next()) + { + anExpectedEdges.Append(TopoDS::Edge(anEdgeExp.Current())); + } + } + + bool hasSeamOccurrences = false; + for (size_t anIndex = 1; anIndex <= anExpectedEdges.Size(); ++anIndex) + { + for (size_t anOther = anIndex + 1; anOther <= anExpectedEdges.Size(); ++anOther) + { + if (anExpectedEdges.Value(anIndex).IsSame(anExpectedEdges.Value(anOther))) + { + hasSeamOccurrences = true; + } + } + } + ASSERT_TRUE(hasSeamOccurrences); + + BRepClass_FaceExplorer anExplorer(aFace); + size_t anExpectedIndex = 1; + for (anExplorer.InitWires(); anExplorer.MoreWires(); anExplorer.NextWire()) + { + for (anExplorer.InitEdges(); anExplorer.MoreEdges(); anExplorer.NextEdge()) + { + ASSERT_LE(anExpectedIndex, anExpectedEdges.Size()); + BRepClass_Edge anEdgeData; + TopAbs_Orientation anOrientation = TopAbs_EXTERNAL; + anExplorer.CurrentEdge(anEdgeData, anOrientation); + const BRepClass_Edge& aConstEdgeData = anEdgeData; + EXPECT_TRUE(aConstEdgeData.Edge().IsEqual(anExpectedEdges.Value(anExpectedIndex))); + EXPECT_EQ(anOrientation, anExpectedEdges.Value(anExpectedIndex).Orientation()); + ++anExpectedIndex; + } + } + EXPECT_EQ(anExpectedIndex, anExpectedEdges.Size() + 1); + expectNextEdgesMatchReference(aFace); + EXPECT_EQ(classify(aFace, gp_Pnt2d(M_PI, 2.5), false), TopAbs_IN); + EXPECT_EQ(classify(aFace, gp_Pnt2d(M_PI, 2.5), true), TopAbs_IN); + EXPECT_EQ(classify(aFace, gp_Pnt2d(M_PI, 6.0), true), TopAbs_OUT); +} + +TEST(BRepClassFaceExplorerTest, ReversedWire_NextEdgesMatchReference) +{ + const TopoDS_Wire aWire = makeRectangle(0.0, 0.0, 2.0, 2.0); + ASSERT_FALSE(aWire.IsNull()); + const TopoDS_Wire aReversedWire = TopoDS::Wire(aWire.Reversed()); + BRepBuilderAPI_MakeFace aFaceBuilder(aReversedWire, true); + ASSERT_TRUE(aFaceBuilder.IsDone()); + + expectNextEdgesMatchReference(aFaceBuilder.Face()); +} + +TEST(BRepClassFaceExplorerTest, SingleClosedEdge_HasNoNextEdge) +{ + occ::handle aCircle = new Geom_Circle(gp_Ax2(gp::Origin(), gp::DZ()), 2.0); + BRepBuilderAPI_MakeEdge anEdgeBuilder(aCircle); + ASSERT_TRUE(anEdgeBuilder.IsDone()); + BRepBuilderAPI_MakeWire aWireBuilder(anEdgeBuilder.Edge()); + ASSERT_TRUE(aWireBuilder.IsDone()); + BRepBuilderAPI_MakeFace aFaceBuilder(aWireBuilder.Wire(), true); + ASSERT_TRUE(aFaceBuilder.IsDone()); + + BRepClass_FaceExplorer anExplorer(aFaceBuilder.Face()); + anExplorer.InitWires(); + ASSERT_TRUE(anExplorer.MoreWires()); + anExplorer.InitEdges(); + ASSERT_TRUE(anExplorer.MoreEdges()); + BRepClass_Edge anEdgeData; + TopAbs_Orientation anOrientation = TopAbs_EXTERNAL; + anExplorer.CurrentEdge(anEdgeData, anOrientation); + EXPECT_TRUE(anEdgeData.NextEdge().IsNull()); + anExplorer.NextEdge(); + EXPECT_FALSE(anExplorer.MoreEdges()); + expectNextEdgesMatchReference(aFaceBuilder.Face()); +} + +TEST(BRepClassFaceExplorerTest, DegenerateEdges_HaveNoNextEdge) +{ + BRepPrimAPI_MakeSphere aSphereBuilder(2.0); + aSphereBuilder.Build(); + ASSERT_TRUE(aSphereBuilder.IsDone()); + const TopoDS_Shape aSphere = aSphereBuilder.Shape(); + ASSERT_FALSE(aSphere.IsNull()); + + bool hasDegenerateEdge = false; + for (TopExp_Explorer aFaceExp(aSphere, TopAbs_FACE); aFaceExp.More(); aFaceExp.Next()) + { + const TopoDS_Face& aFace = TopoDS::Face(aFaceExp.Current()); + BRepClass_FaceExplorer anExplorer(aFace); + for (anExplorer.InitWires(); anExplorer.MoreWires(); anExplorer.NextWire()) + { + for (anExplorer.InitEdges(); anExplorer.MoreEdges(); anExplorer.NextEdge()) + { + BRepClass_Edge anEdgeData; + TopAbs_Orientation anOrientation = TopAbs_EXTERNAL; + anExplorer.CurrentEdge(anEdgeData, anOrientation); + if (BRep_Tool::Degenerated(anEdgeData.Edge())) + { + hasDegenerateEdge = true; + EXPECT_TRUE(anEdgeData.NextEdge().IsNull()); + } + } + } + expectNextEdgesMatchReference(aFace); + } + EXPECT_TRUE(hasDegenerateEdge); +} + +TEST(BRepClassEdgeTest, SetTopology_InvalidatesDerivedData) +{ + const TopoDS_Wire aWire = makeRectangle(0.0, 0.0, 2.0, 2.0); + ASSERT_FALSE(aWire.IsNull()); + BRepBuilderAPI_MakeFace aFaceBuilder(aWire, true); + ASSERT_TRUE(aFaceBuilder.IsDone()); + const TopoDS_Face aFace = aFaceBuilder.Face(); + + TopExp_Explorer anEdgeExp(aFace, TopAbs_EDGE); + ASSERT_TRUE(anEdgeExp.More()); + const TopoDS_Edge anEdge = TopoDS::Edge(anEdgeExp.Current()); + anEdgeExp.Next(); + ASSERT_TRUE(anEdgeExp.More()); + const TopoDS_Edge aNextEdge = TopoDS::Edge(anEdgeExp.Current()); + + double aFirst = 0.0; + double aLast = 0.0; + const occ::handle& aCurve = BRep_Tool::CurveOnSurface(anEdge, aFace, aFirst, aLast); + ASSERT_FALSE(aCurve.IsNull()); + + BRepClass_Edge anEdgeData(anEdge, aFace); + anEdgeData.SetGeometry(aCurve, aFirst, aLast); + Bnd_Box2d aBox; + aBox.Update(0.0, 0.0, 2.0, 2.0); + anEdgeData.SetBoundingBox(aBox); + anEdgeData.SetNextEdge(aNextEdge); + + anEdgeData.SetEdge(TopoDS::Edge(anEdge.Reversed())); + EXPECT_TRUE(anEdgeData.Curve().IsNull()); + EXPECT_TRUE(anEdgeData.NextEdge().IsNull()); + EXPECT_TRUE(anEdgeData.BoundingBox().IsVoid()); + EXPECT_EQ(anEdgeData.BoundingBoxState(), BRepClass_Edge::BndBoxState::NotBuilt); + EXPECT_EQ(anEdgeData.FirstParameter(), 0.0); + EXPECT_EQ(anEdgeData.LastParameter(), 0.0); + + anEdgeData.SetGeometry(aCurve, aFirst, aLast); + anEdgeData.SetBoundingBox(aBox); + anEdgeData.SetFace(TopoDS::Face(aFace.Reversed())); + EXPECT_TRUE(anEdgeData.Curve().IsNull()); + EXPECT_TRUE(anEdgeData.BoundingBox().IsVoid()); + EXPECT_EQ(anEdgeData.BoundingBoxState(), BRepClass_Edge::BndBoxState::NotBuilt); +} + +TEST(BRepClassEdgeTest, VoidBoundingBox_BecomesUnavailable) +{ + BRepClass_Edge anEdgeData; + anEdgeData.SetBoundingBox(Bnd_Box2d()); + EXPECT_TRUE(anEdgeData.BoundingBox().IsVoid()); + EXPECT_EQ(anEdgeData.BoundingBoxState(), BRepClass_Edge::BndBoxState::Unavailable); +} + +TEST(BRepClassIntersectorTest, NotBuiltAndUnavailableBoundingBoxes_IntersectWithoutCacheWrites) +{ + const TopoDS_Wire aWire = makeRectangle(0.0, 0.0, 2.0, 2.0); + ASSERT_FALSE(aWire.IsNull()); + BRepBuilderAPI_MakeFace aFaceBuilder(aWire, true); + ASSERT_TRUE(aFaceBuilder.IsDone()); + const TopoDS_Face aFace = aFaceBuilder.Face(); + + TopExp_Explorer anEdgeExp(aFace, TopAbs_EDGE); + ASSERT_TRUE(anEdgeExp.More()); + const TopoDS_Edge anEdge = TopoDS::Edge(anEdgeExp.Current()); + double aFirst = 0.0; + double aLast = 0.0; + const occ::handle& aCurve = BRep_Tool::CurveOnSurface(anEdge, aFace, aFirst, aLast); + ASSERT_FALSE(aCurve.IsNull()); + + BRepClass_Edge anEdgeData(anEdge, aFace); + anEdgeData.SetGeometry(aCurve, aFirst, aLast); + anEdgeData.SetUseBndBox(true); + + gp_Pnt2d aPoint; + gp_Vec2d aTangent; + aCurve->D1(0.5 * (aFirst + aLast), aPoint, aTangent); + ASSERT_GT(aTangent.SquareMagnitude(), Precision::SquarePConfusion()); + const gp_Dir2d aNormal(-aTangent.Y(), aTangent.X()); + + BRepClass_Intersector aNotBuiltIntersector; + aNotBuiltIntersector.Perform(gp_Lin2d(aPoint, aNormal), + RealLast(), + Precision::Confusion(), + anEdgeData); + EXPECT_TRUE(aNotBuiltIntersector.IsDone()); + EXPECT_EQ(anEdgeData.BoundingBoxState(), BRepClass_Edge::BndBoxState::NotBuilt); + + anEdgeData.SetBoundingBoxUnavailable(); + BRepClass_Intersector anUnavailableIntersector; + anUnavailableIntersector.Perform(gp_Lin2d(aPoint, aNormal), + RealLast(), + Precision::Confusion(), + anEdgeData); + EXPECT_TRUE(anUnavailableIntersector.IsDone()); + EXPECT_EQ(anEdgeData.BoundingBoxState(), BRepClass_Edge::BndBoxState::Unavailable); +} diff --git a/src/ModelingAlgorithms/TKTopAlgo/GTests/BRepTopAdaptor_FClass2d_Test.cxx b/src/ModelingAlgorithms/TKTopAlgo/GTests/BRepTopAdaptor_FClass2d_Test.cxx new file mode 100644 index 0000000000..a1fba5968e --- /dev/null +++ b/src/ModelingAlgorithms/TKTopAlgo/GTests/BRepTopAdaptor_FClass2d_Test.cxx @@ -0,0 +1,202 @@ +// Copyright (c) 2026 OPEN CASCADE SAS +// +// This file is part of Open CASCADE Technology software library. +// +// This library is free software; you can redistribute it and/or modify it under +// the terms of the GNU Lesser General Public License version 2.1 as published +// by the Free Software Foundation, with special exception defined in the file +// OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT +// distribution for complete text of the license and disclaimer of any warranty. +// +// Alternatively, this file may be used under the terms of Open CASCADE +// commercial license or contractual agreement. + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +namespace +{ +TopoDS_Wire makeRectangle(const double theXMin, + const double theYMin, + const double theXMax, + const double theYMax) +{ + const gp_Pnt aPoints[] = {gp_Pnt(theXMin, theYMin, 0.0), + gp_Pnt(theXMax, theYMin, 0.0), + gp_Pnt(theXMax, theYMax, 0.0), + gp_Pnt(theXMin, theYMax, 0.0)}; + BRepBuilderAPI_MakeWire aWireBuilder; + for (int anIndex = 0; anIndex < 4; ++anIndex) + { + BRepBuilderAPI_MakeEdge anEdgeBuilder(aPoints[anIndex], aPoints[(anIndex + 1) % 4]); + if (!anEdgeBuilder.IsDone()) + { + return TopoDS_Wire(); + } + aWireBuilder.Add(anEdgeBuilder.Edge()); + } + return aWireBuilder.IsDone() ? aWireBuilder.Wire() : TopoDS_Wire(); +} + +TopoDS_Face makeFaceWithHole() +{ + const TopoDS_Wire anInnerWire = makeRectangle(3.0, 2.0, 9.0, 6.0); + if (anInnerWire.IsNull()) + { + return TopoDS_Face(); + } + + BRepBuilderAPI_MakeFace aFaceBuilder(gp_Pln(gp::XOY()), 0.0, 12.0, 0.0, 8.0); + if (!aFaceBuilder.IsDone()) + { + return TopoDS_Face(); + } + aFaceBuilder.Add(TopoDS::Wire(anInnerWire.Reversed())); + return aFaceBuilder.IsDone() ? aFaceBuilder.Face() : TopoDS_Face(); +} + +TopoDS_Face makeFaceWithDisconnectedWire() +{ + BRepBuilderAPI_MakeFace aFaceBuilder(gp_Pln(gp::XOY()), 0.0, 10.0, 0.0, 10.0); + if (!aFaceBuilder.IsDone()) + { + return TopoDS_Face(); + } + + TopoDS_Face aFace = aFaceBuilder.Face(); + TopoDS_Edge aFirstEdge; + TopoDS_Edge aThirdEdge; + int anEdgeIndex = 0; + for (TopExp_Explorer anExplorer(aFace, TopAbs_EDGE); anExplorer.More(); anExplorer.Next()) + { + if (anEdgeIndex == 0) + { + aFirstEdge = TopoDS::Edge(anExplorer.Current()); + } + else if (anEdgeIndex == 2) + { + aThirdEdge = TopoDS::Edge(anExplorer.Current()); + } + ++anEdgeIndex; + } + if (aFirstEdge.IsNull() || aThirdEdge.IsNull()) + { + return TopoDS_Face(); + } + + BRep_Builder aBuilder; + TopoDS_Wire aDisconnectedWire; + aBuilder.MakeWire(aDisconnectedWire); + aBuilder.Add(aDisconnectedWire, aFirstEdge); + aBuilder.Add(aDisconnectedWire, aThirdEdge); + aBuilder.Add(aFace, aDisconnectedWire); + return aFace; +} +} // namespace + +TEST(BRepTopAdaptorFClass2dTest, LinearVectorView_ZeroBasedOrderIsClassifiedCorrectly) +{ + NCollection_LinearVector aPoints; + aPoints.Append(gp_Pnt2d(0.0, 0.0)); + aPoints.Append(gp_Pnt2d(7.0, 0.0)); + aPoints.Append(gp_Pnt2d(7.0, 3.0)); + aPoints.Append(gp_Pnt2d(0.0, 3.0)); + aPoints.Append(gp_Pnt2d(0.0, 0.0)); + + NCollection_Array1 aView = aPoints.ToArray1(); + ASSERT_EQ(aView.Lower(), 0); + ASSERT_EQ(aView.Upper(), 4); + EXPECT_TRUE(aView(0).IsEqual(aPoints[0], Precision::PConfusion())); + EXPECT_TRUE(aView(3).IsEqual(aPoints[3], Precision::PConfusion())); + + CSLib_Class2d aClassifier(aPoints.ToArray1(), 1.e-7, 1.e-7, 0.0, 0.0, 7.0, 3.0); + EXPECT_EQ(aClassifier.SiDans(gp_Pnt2d(2.0, 1.0)), CSLib_Class2d::Result_Inside); + EXPECT_EQ(aClassifier.SiDans(gp_Pnt2d(8.0, 1.0)), CSLib_Class2d::Result_Outside); +} + +TEST(BRepTopAdaptorFClass2dTest, PlanarHole_PerformAndRestrictionPreserveWireRoles) +{ + const TopoDS_Face aFace = makeFaceWithHole(); + ASSERT_FALSE(aFace.IsNull()); + BRepTopAdaptor_FClass2d aClassifier(aFace, Precision::PConfusion()); + + const gp_Pnt2d aMaterialPoint(1.0, 1.0); + const gp_Pnt2d aHolePoint(5.0, 4.0); + const gp_Pnt2d anOutsidePoint(14.0, 4.0); + EXPECT_EQ(aClassifier.Perform(aMaterialPoint), TopAbs_IN); + EXPECT_EQ(aClassifier.Perform(aHolePoint), TopAbs_OUT); + EXPECT_EQ(aClassifier.Perform(anOutsidePoint), TopAbs_OUT); + EXPECT_EQ(aClassifier.TestOnRestriction(aMaterialPoint, 1.e-7), TopAbs_IN); + EXPECT_EQ(aClassifier.TestOnRestriction(aHolePoint, 1.e-7), TopAbs_OUT); + EXPECT_EQ(aClassifier.TestOnRestriction(anOutsidePoint, 1.e-7), TopAbs_OUT); +} + +TEST(BRepTopAdaptorFClass2dTest, Boundary_PerformUsesRepeatableExactFallback) +{ + const TopoDS_Face aFace = makeFaceWithHole(); + ASSERT_FALSE(aFace.IsNull()); + BRepTopAdaptor_FClass2d aClassifier(aFace, Precision::PConfusion()); + const gp_Pnt2d aBoundaryPoint(0.0, 3.0); + + BRepClass_FaceClassifier anExact(aFace, aBoundaryPoint, Precision::PConfusion()); + ASSERT_EQ(anExact.State(), TopAbs_ON); + for (int anIteration = 0; anIteration < 20; ++anIteration) + { + EXPECT_EQ(aClassifier.Perform(aBoundaryPoint), anExact.State()); + EXPECT_EQ(aClassifier.TestOnRestriction(aBoundaryPoint, 1.e-6), TopAbs_ON); + } +} + +TEST(BRepTopAdaptorFClass2dTest, PeriodicCylinder_RecadresEquivalentParameters) +{ + occ::handle aSurface = + new Geom_CylindricalSurface(gp_Ax3(gp::Origin(), gp::DZ()), 2.0); + ASSERT_FALSE(aSurface.IsNull()); + BRepBuilderAPI_MakeFace aFaceBuilder(aSurface, 0.0, 2.0 * M_PI, 0.0, 5.0, Precision::Confusion()); + ASSERT_TRUE(aFaceBuilder.IsDone()); + BRepTopAdaptor_FClass2d aClassifier(aFaceBuilder.Face(), Precision::PConfusion()); + + const gp_Pnt2d aCanonicalPoint(0.75, 2.0); + const gp_Pnt2d aShiftedPoint(aCanonicalPoint.X() + 4.0 * M_PI, aCanonicalPoint.Y()); + EXPECT_EQ(aClassifier.Perform(aCanonicalPoint), TopAbs_IN); + EXPECT_EQ(aClassifier.Perform(aShiftedPoint), aClassifier.Perform(aCanonicalPoint)); + EXPECT_EQ(aClassifier.TestOnRestriction(aShiftedPoint, 1.e-7), + aClassifier.TestOnRestriction(aCanonicalPoint, 1.e-7)); +} + +TEST(BRepTopAdaptorFClass2dTest, DisconnectedWire_FallsBackToStableExactClassification) +{ + const TopoDS_Face aFace = makeFaceWithDisconnectedWire(); + ASSERT_FALSE(aFace.IsNull()); + const gp_Pnt2d aPoint(5.0, 5.0); + BRepClass_FaceClassifier anExact(aFace, aPoint, Precision::PConfusion()); + BRepTopAdaptor_FClass2d aClassifier(aFace, Precision::PConfusion()); + + const TopAbs_State anExpected = anExact.State(); + for (int anIteration = 0; anIteration < 10; ++anIteration) + { + EXPECT_EQ(aClassifier.Perform(aPoint, false), anExpected); + EXPECT_EQ(aClassifier.TestOnRestriction(aPoint, Precision::PConfusion(), false), anExpected); + } +} diff --git a/src/ModelingAlgorithms/TKTopAlgo/GTests/FILES.cmake b/src/ModelingAlgorithms/TKTopAlgo/GTests/FILES.cmake index 1702676551..09fd1dc39f 100644 --- a/src/ModelingAlgorithms/TKTopAlgo/GTests/FILES.cmake +++ b/src/ModelingAlgorithms/TKTopAlgo/GTests/FILES.cmake @@ -9,8 +9,10 @@ set(OCCT_TKTopAlgo_GTests_FILES BRepBuilderAPI_Transform_Test.cxx BRepCheck_Face_Test.cxx BRepClass3d_SolidClassifier_Test.cxx + BRepClass_FaceExplorer_Test.cxx BRepExtrema_DistShapeShape_Test.cxx BRepGProp_Test.cxx BRepLib_MakeWire_Test.cxx BRepOffsetAPI_ThruSections_Test.cxx + BRepTopAdaptor_FClass2d_Test.cxx )