diff --git a/src/FoundationClasses/TKMath/MathRoot/FILES.cmake b/src/FoundationClasses/TKMath/MathRoot/FILES.cmake index 504f89d025..9698ece062 100644 --- a/src/FoundationClasses/TKMath/MathRoot/FILES.cmake +++ b/src/FoundationClasses/TKMath/MathRoot/FILES.cmake @@ -6,6 +6,7 @@ set(OCCT_MathRoot_FILES MathRoot_Bisection.hxx MathRoot_Brent.hxx MathRoot_Multiple.hxx + MathRoot_MultipleUtils.hxx MathRoot_Newton.hxx MathRoot_Secant.hxx MathRoot_Trig.hxx diff --git a/src/FoundationClasses/TKMath/MathRoot/MathRoot_Multiple.hxx b/src/FoundationClasses/TKMath/MathRoot/MathRoot_Multiple.hxx index 5e4f05774a..0e6b6cc4cb 100644 --- a/src/FoundationClasses/TKMath/MathRoot/MathRoot_Multiple.hxx +++ b/src/FoundationClasses/TKMath/MathRoot/MathRoot_Multiple.hxx @@ -14,16 +14,7 @@ #ifndef _MathRoot_Multiple_HeaderFile #define _MathRoot_Multiple_HeaderFile -#include -#include -#include -#include -#include -#include - -#include - -#include +#include //! @file MathRoot_Multiple.hxx //! @brief Algorithms for finding all roots of a function in a given range. @@ -33,41 +24,6 @@ namespace MathRoot { -using namespace MathUtils; - -//! Result for multiple root finding. -//! Contains all found roots sorted in ascending order. -struct MultipleResult -{ - MathUtils::Status Status = MathUtils::Status::NotConverged; //!< Computation status - size_t NbIterations = 0; //!< Total iterations across all roots - NCollection_Vector Roots; //!< Found roots (sorted) - NCollection_Vector Values; //!< Function values at roots - bool IsAllNull = false; //!< True if function is essentially zero in range - - //! Returns true if computation succeeded. - bool IsDone() const { return Status == MathUtils::Status::OK; } - - //! Conversion to bool for convenient checking. - explicit operator bool() const { return IsDone(); } - - //! Returns the number of roots found. - int NbRoots() const { return Roots.Length(); } - - //! Access root by index (0-based). - double operator[](int theIndex) const { return Roots.Value(theIndex); } -}; - -//! Configuration for multiple root finding. -struct MultipleConfig -{ - int NbSamples = 100; //!< Number of sample points for initial search - double XTolerance = 1e-10; //!< Tolerance on X for convergence - double FTolerance = 1e-10; //!< Tolerance on F(X) for convergence - double NullTolerance = 1e-12; //!< Tolerance to consider function as null - int MaxIterations = 100; //!< Max iterations per root refinement - double Offset = 0.0; //!< Find roots of f(x) - Offset = 0 -}; //! Finds all real roots of a function within the range [theLower, theUpper]. //! Uses uniform sampling to detect sign changes, then refines each root using Brent's method. @@ -90,178 +46,22 @@ MultipleResult FindAllRoots(Function& theFunc, double theUpper, const MultipleConfig& theConfig = MultipleConfig()) { - MultipleResult aResult; - aResult.Status = MathUtils::Status::OK; - - // Ensure proper ordering - double aLower = std::min(theLower, theUpper); - double aUpper = std::max(theLower, theUpper); - - // Minimum samples - const int aNbSamples = std::max(theConfig.NbSamples, 10); - const double aDx = (aUpper - aLower) / aNbSamples; - - // Ensure EpsX is not too small relative to interval - const double aMinEpsX = 1e-10 * (std::abs(aLower) + std::abs(aUpper)); - const double aEpsX = std::max(theConfig.XTolerance, aMinEpsX); - - // Sample function values + const int aNbSamples = std::max(theConfig.NbSamples, 10); math_Vector aSamples(0, aNbSamples); - math_Vector aXValues(0, aNbSamples); - bool aAllValid = true; - for (int i = 0; i <= aNbSamples; ++i) - { - double aX = aLower + i * aDx; - if (aX > aUpper) - aX = aUpper; - aXValues(i) = aX; + MultipleSampleValueFn aSampleFn{theFunc, aSamples, theConfig.Offset}; + MultipleGetValueFn aGetValue{aSamples}; + MultipleBrentValueWrapper aWrapper{theFunc, theConfig.Offset}; + MultipleGetRootValueFn aGetRootValue{theFunc}; - double aF = 0.0; - if (!theFunc.Value(aX, aF)) - { - aAllValid = false; - break; - } - aSamples(i) = aF - theConfig.Offset; - } - - if (!aAllValid) - { - aResult.Status = MathUtils::Status::NumericalError; - return aResult; - } - - // Check if function is essentially null everywhere - aResult.IsAllNull = true; - for (int i = 0; i <= aNbSamples; ++i) - { - if (std::abs(aSamples(i)) > theConfig.NullTolerance) - { - aResult.IsAllNull = false; - break; - } - } - - if (aResult.IsAllNull) - { - return aResult; - } - - // Helper to add root if not duplicate - auto addRoot = [&](double theRoot, double theValue) { - // Check for duplicates - for (int k = 0; k < aResult.Roots.Length(); ++k) - { - if (std::abs(theRoot - aResult.Roots.Value(k)) < aEpsX) - { - return; - } - } - aResult.Roots.Append(theRoot); - aResult.Values.Append(theValue); - }; - - // Create wrapper for Brent that uses Value-only interface - class BrentWrapper - { - public: - BrentWrapper(Function& theF, double theOffset) - : myFunc(theF), - myOffset(theOffset) - { - } - - bool Value(double theX, double& theY) const - { - if (!myFunc.Value(theX, theY)) - return false; - theY -= myOffset; - return true; - } - - private: - Function& myFunc; - double myOffset; - }; - - BrentWrapper aWrapper(theFunc, theConfig.Offset); - - // Find sign changes - for (int i = 0; i < aNbSamples; ++i) - { - const double aF0 = aSamples(i); - const double aF1 = aSamples(i + 1); - const double aX0 = aXValues(i); - const double aX1 = aXValues(i + 1); - - // Exact zero at sample point - if (std::abs(aF0) < theConfig.FTolerance) - { - addRoot(aX0, aF0 + theConfig.Offset); - continue; - } - - // Sign change detected - if (aF0 * aF1 < 0.0) - { - MathUtils::Config aBrentConfig; - aBrentConfig.XTolerance = aEpsX; - aBrentConfig.FTolerance = theConfig.FTolerance; - aBrentConfig.MaxIterations = theConfig.MaxIterations; - - auto aBrentResult = Brent(aWrapper, aX0, aX1, aBrentConfig); - aResult.NbIterations += aBrentResult.NbIterations; - - if (aBrentResult.IsDone() && aBrentResult.Root.has_value()) - { - double aRootValue = 0.0; - theFunc.Value(*aBrentResult.Root, aRootValue); - addRoot(*aBrentResult.Root, aRootValue); - } - } - } - - // Check last sample point - if (std::abs(aSamples(aNbSamples)) < theConfig.FTolerance) - { - addRoot(aXValues(aNbSamples), aSamples(aNbSamples) + theConfig.Offset); - } - - // Sort roots using indices - const int aNbRoots = aResult.Roots.Length(); - if (aNbRoots > 1) - { - math_IntegerVector aIndices(0, aNbRoots - 1); - for (int i = 0; i < aNbRoots; ++i) - { - aIndices(i) = i; - } - - // Simple insertion sort for small arrays - for (int i = 1; i < aNbRoots; ++i) - { - int aKey = aIndices(i); - int j = i - 1; - while (j >= 0 && aResult.Roots.Value(aIndices(j)) > aResult.Roots.Value(aKey)) - { - aIndices(j + 1) = aIndices(j); - --j; - } - aIndices(j + 1) = aKey; - } - - NCollection_Vector aSortedRoots, aSortedValues; - for (int i = 0; i < aNbRoots; ++i) - { - aSortedRoots.Append(aResult.Roots.Value(aIndices(i))); - aSortedValues.Append(aResult.Values.Value(aIndices(i))); - } - aResult.Roots = aSortedRoots; - aResult.Values = aSortedValues; - } - - return aResult; + return FindAllRootsImpl(theLower, + theUpper, + theConfig, + aSampleFn, + aGetValue, + aWrapper, + aGetRootValue, + MultipleNoExtraHandler()); } //! Finds all real roots of a function with derivative within range [theLower, theUpper]. @@ -281,244 +81,27 @@ MultipleResult FindAllRootsWithDerivative(Function& theFunc, double theUpper, const MultipleConfig& theConfig = MultipleConfig()) { - MultipleResult aResult; - aResult.Status = MathUtils::Status::OK; - - double aLower = std::min(theLower, theUpper); - double aUpper = std::max(theLower, theUpper); - - const int aNbSamples = std::max(theConfig.NbSamples, 10); - const double aDx = (aUpper - aLower) / aNbSamples; - - const double aMinEpsX = 1e-10 * (std::abs(aLower) + std::abs(aUpper)); - const double aEpsX = std::max(theConfig.XTolerance, aMinEpsX); - - // Sample function values and derivatives + const int aNbSamples = std::max(theConfig.NbSamples, 10); math_Vector aFValues(0, aNbSamples); math_Vector aDFValues(0, aNbSamples); - math_Vector aXValues(0, aNbSamples); - bool aAllValid = true; - for (int i = 0; i <= aNbSamples; ++i) - { - double aX = aLower + i * aDx; - if (aX > aUpper) - aX = aUpper; - aXValues(i) = aX; + MultipleSampleDerivFn aSampleFn{theFunc, aFValues, aDFValues, theConfig.Offset}; + MultipleGetValueFn aGetValue{aFValues}; + MultipleBrentDerivWrapper aWrapper{theFunc, theConfig.Offset}; + MultipleGetRootDerivFn aGetRootValue{theFunc}; + MultipleTangentialHandler aTangentialExtra{theFunc, + aDFValues, + theConfig.Offset, + theConfig.FTolerance}; - double aF = 0.0, aDF = 0.0; - if (!theFunc.Values(aX, aF, aDF)) - { - aAllValid = false; - break; - } - aFValues(i) = aF - theConfig.Offset; - aDFValues(i) = aDF; - } - - if (!aAllValid) - { - aResult.Status = MathUtils::Status::NumericalError; - return aResult; - } - - // Check if function is essentially null - aResult.IsAllNull = true; - for (int i = 0; i <= aNbSamples; ++i) - { - if (std::abs(aFValues(i)) > theConfig.NullTolerance) - { - aResult.IsAllNull = false; - break; - } - } - - if (aResult.IsAllNull) - { - return aResult; - } - - // Helper to add root if not duplicate - auto addRoot = [&](double theRoot, double theValue) { - for (int k = 0; k < aResult.Roots.Length(); ++k) - { - if (std::abs(theRoot - aResult.Roots.Value(k)) < aEpsX) - { - return; - } - } - aResult.Roots.Append(theRoot); - aResult.Values.Append(theValue); - }; - - // Wrapper for Brent using Values interface - class BrentWrapper - { - public: - BrentWrapper(Function& theF, double theOffset) - : myFunc(theF), - myOffset(theOffset) - { - } - - bool Value(double theX, double& theY) const - { - double aDF = 0.0; - if (!myFunc.Values(theX, theY, aDF)) - return false; - theY -= myOffset; - return true; - } - - private: - Function& myFunc; - double myOffset; - }; - - BrentWrapper aWrapper(theFunc, theConfig.Offset); - - // Find sign changes - for (int i = 0; i < aNbSamples; ++i) - { - const double aF0 = aFValues(i); - const double aF1 = aFValues(i + 1); - const double aX0 = aXValues(i); - const double aX1 = aXValues(i + 1); - - // Exact zero at sample point - if (std::abs(aF0) < theConfig.FTolerance) - { - addRoot(aX0, aF0 + theConfig.Offset); - continue; - } - - // Sign change - if (aF0 * aF1 < 0.0) - { - MathUtils::Config aBrentConfig; - aBrentConfig.XTolerance = aEpsX; - aBrentConfig.FTolerance = theConfig.FTolerance; - aBrentConfig.MaxIterations = theConfig.MaxIterations; - - auto aBrentResult = Brent(aWrapper, aX0, aX1, aBrentConfig); - aResult.NbIterations += aBrentResult.NbIterations; - - if (aBrentResult.IsDone() && aBrentResult.Root.has_value()) - { - double aRootValue = 0.0, aDummy = 0.0; - theFunc.Values(*aBrentResult.Root, aRootValue, aDummy); - addRoot(*aBrentResult.Root, aRootValue); - } - } - // Check for potential extrema touching zero - else if (aF0 > 0.0 && aF1 > 0.0) - { - // Potential minimum - check if derivative changes sign - if (aDFValues(i) < 0.0 && aDFValues(i + 1) > 0.0) - { - // Find the minimum using bisection on derivative - double aXL = aX0, aXR = aX1; - - for (int anIter = 0; anIter < 20; ++anIter) - { - double aXM = 0.5 * (aXL + aXR); - double aFM = 0.0, aDFM = 0.0; - if (!theFunc.Values(aXM, aFM, aDFM)) - break; - aFM -= theConfig.Offset; - - if (aDFM < 0.0) - { - aXL = aXM; - } - else - { - aXR = aXM; - } - - // Check if minimum is close enough to zero - if (std::abs(aFM) < theConfig.FTolerance) - { - addRoot(aXM, aFM + theConfig.Offset); - break; - } - } - } - } - else if (aF0 < 0.0 && aF1 < 0.0) - { - // Potential maximum - check if derivative changes sign - if (aDFValues(i) > 0.0 && aDFValues(i + 1) < 0.0) - { - double aXL = aX0, aXR = aX1; - - for (int anIter = 0; anIter < 20; ++anIter) - { - double aXM = 0.5 * (aXL + aXR); - double aFM = 0.0, aDFM = 0.0; - if (!theFunc.Values(aXM, aFM, aDFM)) - break; - aFM -= theConfig.Offset; - - if (aDFM > 0.0) - { - aXL = aXM; - } - else - { - aXR = aXM; - } - - if (std::abs(aFM) < theConfig.FTolerance) - { - addRoot(aXM, aFM + theConfig.Offset); - break; - } - } - } - } - } - - // Check last sample point - if (std::abs(aFValues(aNbSamples)) < theConfig.FTolerance) - { - addRoot(aXValues(aNbSamples), aFValues(aNbSamples) + theConfig.Offset); - } - - // Sort roots using indices - const int aNbRoots = aResult.Roots.Length(); - if (aNbRoots > 1) - { - math_IntegerVector aIndices(0, aNbRoots - 1); - for (int i = 0; i < aNbRoots; ++i) - { - aIndices(i) = i; - } - - // Simple insertion sort for small arrays - for (int i = 1; i < aNbRoots; ++i) - { - int aKey = aIndices(i); - int j = i - 1; - while (j >= 0 && aResult.Roots.Value(aIndices(j)) > aResult.Roots.Value(aKey)) - { - aIndices(j + 1) = aIndices(j); - --j; - } - aIndices(j + 1) = aKey; - } - - NCollection_Vector aSortedRoots, aSortedValues; - for (int i = 0; i < aNbRoots; ++i) - { - aSortedRoots.Append(aResult.Roots.Value(aIndices(i))); - aSortedValues.Append(aResult.Values.Value(aIndices(i))); - } - aResult.Roots = aSortedRoots; - aResult.Values = aSortedValues; - } - - return aResult; + return FindAllRootsImpl(theLower, + theUpper, + theConfig, + aSampleFn, + aGetValue, + aWrapper, + aGetRootValue, + aTangentialExtra); } //! Convenience alias using default configuration. diff --git a/src/FoundationClasses/TKMath/MathRoot/MathRoot_MultipleUtils.hxx b/src/FoundationClasses/TKMath/MathRoot/MathRoot_MultipleUtils.hxx new file mode 100644 index 0000000000..3c33f72356 --- /dev/null +++ b/src/FoundationClasses/TKMath/MathRoot/MathRoot_MultipleUtils.hxx @@ -0,0 +1,443 @@ +// Copyright (c) 2025 OPEN CASCADE SAS +// +// This file is part of Open CASCADE Technology software library. +// +// This library is free software; you can redistribute it and/or modify it under +// the terms of the GNU Lesser General Public License version 2.1 as published +// by the Free Software Foundation, with special exception defined in the file +// OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT +// distribution for complete text of the license and disclaimer of any warranty. +// +// Alternatively, this file may be used under the terms of Open CASCADE +// commercial license or contractual agreement. + +#ifndef _MathRoot_MultipleUtils_HeaderFile +#define _MathRoot_MultipleUtils_HeaderFile + +#include +#include +#include +#include + +#include + +#include + +//! @file MathRoot_MultipleUtils.hxx +//! @brief Internal utilities for FindAllRoots / FindAllRootsWithDerivative. +//! +//! Contains result/config types, helper functions, functor adapters and the +//! shared core implementation used by MathRoot_Multiple.hxx. + +namespace MathRoot +{ +using namespace MathUtils; + +// ============================================================================ +// Result and configuration types +// ============================================================================ + +//! Result for multiple root finding. +//! Contains all found roots sorted in ascending order. +struct MultipleResult +{ + MathUtils::Status Status = MathUtils::Status::NotConverged; //!< Computation status + size_t NbIterations = 0; //!< Total iterations across all roots + NCollection_Vector Roots; //!< Found roots (sorted) + NCollection_Vector Values; //!< Function values at roots + bool IsAllNull = false; //!< True if function is essentially zero in range + + //! Returns true if computation succeeded. + bool IsDone() const { return Status == MathUtils::Status::OK; } + + //! Conversion to bool for convenient checking. + explicit operator bool() const { return IsDone(); } + + //! Returns the number of roots found. + int NbRoots() const { return Roots.Length(); } + + //! Access root by index (0-based). + double operator[](int theIndex) const { return Roots.Value(theIndex); } +}; + +//! Configuration for multiple root finding. +struct MultipleConfig +{ + int NbSamples = 100; //!< Number of sample points for initial search + double XTolerance = 1e-10; //!< Tolerance on X for convergence + double FTolerance = 1e-10; //!< Tolerance on F(X) for convergence + double NullTolerance = 1e-12; //!< Tolerance to consider function as null + int MaxIterations = 100; //!< Max iterations per root refinement + double Offset = 0.0; //!< Find roots of f(x) - Offset = 0 +}; + +// ============================================================================ +// Helper functions +// ============================================================================ + +//! In-place insertion sort of roots and corresponding values by ascending root value. +inline void SortRoots(MultipleResult& theResult) +{ + for (int i = 1; i < theResult.Roots.Length(); ++i) + { + const double aKeyRoot = theResult.Roots[i]; + const double aKeyVal = theResult.Values[i]; + int j = i - 1; + while (j >= 0 && theResult.Roots[j] > aKeyRoot) + { + theResult.Roots[j + 1] = theResult.Roots[j]; + theResult.Values[j + 1] = theResult.Values[j]; + --j; + } + theResult.Roots[j + 1] = aKeyRoot; + theResult.Values[j + 1] = aKeyVal; + } +} + +//! Helper to add a root if it is not a duplicate of an already found root. +inline void AddRoot(MultipleResult& theResult, double theEpsX, double theRoot, double theValue) +{ + for (int k = 0; k < theResult.Roots.Length(); ++k) + { + if (std::abs(theRoot - theResult.Roots.Value(k)) < theEpsX) + { + return; + } + } + theResult.Roots.Append(theRoot); + theResult.Values.Append(theValue); +} + +// ============================================================================ +// Functor adapters for Value-only interface +// ============================================================================ + +//! Samples a Value-only function and stores f(x)-offset into a math_Vector. +//! @tparam Function type with Value(double theX, double& theF) method +template +struct MultipleSampleValueFn +{ + Function& myFunc; + math_Vector& mySamples; + const double myOffset; + + bool operator()(int theIndex, double theX) const + { + double aF = 0.0; + if (!myFunc.Value(theX, aF)) + return false; + mySamples(theIndex) = aF - myOffset; + return true; + } +}; + +//! Returns the sampled value at a given index from a math_Vector. +struct MultipleGetValueFn +{ + const math_Vector& mySamples; + + double operator()(int theIndex) const { return mySamples(theIndex); } +}; + +//! Brent wrapper that adapts a Value-only function for offset root finding. +//! @tparam Function type with Value(double theX, double& theF) method +template +struct MultipleBrentValueWrapper +{ + Function& myFunc; + double myOffset; + + bool Value(double theX, double& theY) const + { + if (!myFunc.Value(theX, theY)) + return false; + theY -= myOffset; + return true; + } +}; + +//! Evaluates original (non-offset) function value at a root point via Value interface. +//! @tparam Function type with Value(double theX, double& theF) method +template +struct MultipleGetRootValueFn +{ + Function& myFunc; + + double operator()(double theX) const + { + double aF = 0.0; + myFunc.Value(theX, aF); + return aF; + } +}; + +// ============================================================================ +// Functor adapters for Values (with derivative) interface +// ============================================================================ + +//! Samples a function with derivative and stores f(x)-offset and f'(x) into math_Vectors. +//! @tparam Function type with Values(double theX, double& theF, double& theDF) method +template +struct MultipleSampleDerivFn +{ + Function& myFunc; + math_Vector& myFValues; + math_Vector& myDFValues; + const double myOffset; + + bool operator()(int theIndex, double theX) const + { + double aF = 0.0, aDF = 0.0; + if (!myFunc.Values(theX, aF, aDF)) + return false; + myFValues(theIndex) = aF - myOffset; + myDFValues(theIndex) = aDF; + return true; + } +}; + +//! Brent wrapper that adapts a Values (with derivative) function for offset root finding. +//! @tparam Function type with Values(double theX, double& theF, double& theDF) method +template +struct MultipleBrentDerivWrapper +{ + Function& myFunc; + double myOffset; + + bool Value(double theX, double& theY) const + { + double aDF = 0.0; + if (!myFunc.Values(theX, theY, aDF)) + return false; + theY -= myOffset; + return true; + } +}; + +//! Evaluates original (non-offset) function value at a root point via Values interface. +//! @tparam Function type with Values(double theX, double& theF, double& theDF) method +template +struct MultipleGetRootDerivFn +{ + Function& myFunc; + + double operator()(double theX) const + { + double aF = 0.0, aDF = 0.0; + myFunc.Values(theX, aF, aDF); + return aF; + } +}; + +// ============================================================================ +// Interval handlers +// ============================================================================ + +//! No-op interval handler for functions without derivative. +struct MultipleNoExtraHandler +{ + void operator()(int, double, double, double, double, MultipleResult&, double) const {} +}; + +//! Tangential root detection: finds extrema that touch zero without sign change. +//! Uses derivative sign changes to locate potential minima/maxima, then bisects +//! the derivative to check whether the function value is close enough to zero. +//! @tparam Function type with Values(double theX, double& theF, double& theDF) method +template +struct MultipleTangentialHandler +{ + Function& myFunc; + const math_Vector& myDFValues; + const double myOffset; + const double myFTolerance; + + void operator()(int theIndex, + double theX0, + double theX1, + double theF0, + double theF1, + MultipleResult& theResult, + double theEpsX) const + { + if (theF0 > 0.0 && theF1 > 0.0) + { + // Potential minimum - check if derivative changes sign (negative to positive) + if (myDFValues(theIndex) < 0.0 && myDFValues(theIndex + 1) > 0.0) + { + findTangentialRoot(theX0, theX1, true, theResult, theEpsX); + } + } + else if (theF0 < 0.0 && theF1 < 0.0) + { + // Potential maximum - check if derivative changes sign (positive to negative) + if (myDFValues(theIndex) > 0.0 && myDFValues(theIndex + 1) < 0.0) + { + findTangentialRoot(theX0, theX1, false, theResult, theEpsX); + } + } + } + +private: + //! Bisects derivative to locate an extremum and checks if it touches zero. + //! @param theIsMinimum true for minimum (derivative negative->positive), + //! false for maximum (derivative positive->negative) + void findTangentialRoot(double theX0, + double theX1, + bool theIsMinimum, + MultipleResult& theResult, + double theEpsX) const + { + double aXL = theX0, aXR = theX1; + for (int anIter = 0; anIter < 20; ++anIter) + { + double aXM = 0.5 * (aXL + aXR); + double aFM = 0.0, aDFM = 0.0; + if (!myFunc.Values(aXM, aFM, aDFM)) + break; + aFM -= myOffset; + + // For minimum: derivative goes from negative to positive + // For maximum: derivative goes from positive to negative + const bool isMoveRight = theIsMinimum ? (aDFM < 0.0) : (aDFM > 0.0); + if (isMoveRight) + { + aXL = aXM; + } + else + { + aXR = aXM; + } + + if (std::abs(aFM) < myFTolerance) + { + AddRoot(theResult, theEpsX, aXM, aFM + myOffset); + break; + } + } + } +}; + +// ============================================================================ +// Core implementation +// ============================================================================ + +//! Core implementation for finding all roots in an interval. +//! Shared logic for both Value-only and Values (with derivative) interfaces. +//! +//! @tparam SampleFn callable (int theIndex, double theX) -> bool +//! @tparam GetValueFn callable (int theIndex) -> double, returns f(x)-offset at sample +//! @tparam BrentWrapperT type with Value(double, double&) for Brent root finding +//! @tparam GetRootValueFn callable (double theX) -> double, returns original f(x) +//! @tparam IntervalExtraFn callable (int, x0, x1, f0, f1, result, epsX) -> void +template +MultipleResult FindAllRootsImpl(double theLower, + double theUpper, + const MultipleConfig& theConfig, + SampleFn theSampleFn, + GetValueFn theGetValue, + BrentWrapperT& theBrentWrapper, + GetRootValueFn theGetRootValue, + IntervalExtraFn theIntervalExtra) +{ + MultipleResult aResult; + aResult.Status = MathUtils::Status::OK; + + // Ensure proper ordering + const double aLower = std::min(theLower, theUpper); + const double aUpper = std::max(theLower, theUpper); + + // Minimum samples + const int aNbSamples = std::max(theConfig.NbSamples, 10); + const double aDx = (aUpper - aLower) / aNbSamples; + + // Ensure EpsX is not too small relative to interval + const double aMinEpsX = 1e-10 * (std::abs(aLower) + std::abs(aUpper)); + const double aEpsX = std::max(theConfig.XTolerance, aMinEpsX); + + // Sample function values + math_Vector aXValues(0, aNbSamples); + for (int i = 0; i <= aNbSamples; ++i) + { + double aX = aLower + i * aDx; + if (aX > aUpper) + aX = aUpper; + aXValues(i) = aX; + + if (!theSampleFn(i, aX)) + { + aResult.Status = MathUtils::Status::NumericalError; + return aResult; + } + } + + // Check if function is essentially null everywhere + aResult.IsAllNull = true; + for (int i = 0; i <= aNbSamples; ++i) + { + if (std::abs(theGetValue(i)) > theConfig.NullTolerance) + { + aResult.IsAllNull = false; + break; + } + } + + if (aResult.IsAllNull) + { + return aResult; + } + + // Find sign changes + for (int i = 0; i < aNbSamples; ++i) + { + const double aF0 = theGetValue(i); + const double aF1 = theGetValue(i + 1); + const double aX0 = aXValues(i); + const double aX1 = aXValues(i + 1); + + // Exact zero at sample point + if (std::abs(aF0) < theConfig.FTolerance) + { + AddRoot(aResult, aEpsX, aX0, aF0 + theConfig.Offset); + continue; + } + + // Sign change detected + if (aF0 * aF1 < 0.0) + { + MathUtils::Config aBrentConfig; + aBrentConfig.XTolerance = aEpsX; + aBrentConfig.FTolerance = theConfig.FTolerance; + aBrentConfig.MaxIterations = theConfig.MaxIterations; + + MathUtils::ScalarResult aBrentResult = Brent(theBrentWrapper, aX0, aX1, aBrentConfig); + aResult.NbIterations += aBrentResult.NbIterations; + + if (aBrentResult.IsDone() && aBrentResult.Root.has_value()) + { + AddRoot(aResult, aEpsX, *aBrentResult.Root, theGetRootValue(*aBrentResult.Root)); + } + } + else + { + // Additional per-interval processing (e.g., tangential root detection) + theIntervalExtra(i, aX0, aX1, aF0, aF1, aResult, aEpsX); + } + } + + // Check last sample point + if (std::abs(theGetValue(aNbSamples)) < theConfig.FTolerance) + { + AddRoot(aResult, aEpsX, aXValues(aNbSamples), theGetValue(aNbSamples) + theConfig.Offset); + } + + SortRoots(aResult); + return aResult; +} + +} // namespace MathRoot + +#endif // _MathRoot_MultipleUtils_HeaderFile diff --git a/src/ModelingData/TKG2d/GTests/FILES.cmake b/src/ModelingData/TKG2d/GTests/FILES.cmake index 106aeae91e..811d40c4ce 100644 --- a/src/ModelingData/TKG2d/GTests/FILES.cmake +++ b/src/ModelingData/TKG2d/GTests/FILES.cmake @@ -29,6 +29,9 @@ set(OCCT_TKG2d_GTests_FILES Geom2dEval_TBezierCurve_Test.cxx Geom2dGcc_Circ2d2TanOn_Test.cxx Geom2dGcc_Circ2d2TanRad_Test.cxx + Geom2dProp_Test.cxx + Geom2dProp_VsCLProps2d_Test.cxx + Geom2dProp_VsLProp_Test.cxx Geom2dGridEval_BezierCurve_Test.cxx Geom2dGridEval_Curve_Test.cxx Geom2dGridEval_Ellipse_Test.cxx diff --git a/src/ModelingData/TKG2d/GTests/Geom2dProp_Test.cxx b/src/ModelingData/TKG2d/GTests/Geom2dProp_Test.cxx new file mode 100644 index 0000000000..f9cb894642 --- /dev/null +++ b/src/ModelingData/TKG2d/GTests/Geom2dProp_Test.cxx @@ -0,0 +1,1502 @@ +// Copyright (c) 2025 OPEN CASCADE SAS +// +// This file is part of Open CASCADE Technology software library. +// +// This library is free software; you can redistribute it and/or modify it under +// the terms of the GNU Lesser General Public License version 2.1 as published +// by the Free Software Foundation, with special exception defined in the file +// OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT +// distribution for complete text of the license and disclaimer of any warranty. +// +// Alternatively, this file may be used under the terms of Open CASCADE +// commercial license or contractual agreement. + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +#include + +// ============================================================================ +// Free functions tests (Geom2dProp namespace) +// ============================================================================ + +TEST(Geom2dPropTest, ComputeTangent_FromD1) +{ + const gp_Vec2d aD1(3.0, 4.0); + const gp_Vec2d aD2(0.0, 0.0); + const gp_Vec2d aD3(0.0, 0.0); + + const Geom2dProp::TangentResult aRes = Geom2dProp::ComputeTangent(aD1, aD2, aD3, 1.0e-7); + ASSERT_TRUE(aRes.IsDefined); + EXPECT_NEAR(aRes.Direction.X(), 3.0 / 5.0, Precision::Confusion()); + EXPECT_NEAR(aRes.Direction.Y(), 4.0 / 5.0, Precision::Confusion()); +} + +TEST(Geom2dPropTest, ComputeTangent_FallbackToD2) +{ + const gp_Vec2d aD1(0.0, 0.0); + const gp_Vec2d aD2(1.0, 0.0); + const gp_Vec2d aD3(0.0, 0.0); + + const Geom2dProp::TangentResult aRes = Geom2dProp::ComputeTangent(aD1, aD2, aD3, 1.0e-7); + ASSERT_TRUE(aRes.IsDefined); + EXPECT_NEAR(aRes.Direction.X(), 1.0, Precision::Confusion()); + EXPECT_NEAR(aRes.Direction.Y(), 0.0, Precision::Confusion()); +} + +TEST(Geom2dPropTest, ComputeTangent_FallbackToD3) +{ + const gp_Vec2d aD1(0.0, 0.0); + const gp_Vec2d aD2(0.0, 0.0); + const gp_Vec2d aD3(0.0, 5.0); + + const Geom2dProp::TangentResult aRes = Geom2dProp::ComputeTangent(aD1, aD2, aD3, 1.0e-7); + ASSERT_TRUE(aRes.IsDefined); + EXPECT_NEAR(aRes.Direction.X(), 0.0, Precision::Confusion()); + EXPECT_NEAR(aRes.Direction.Y(), 1.0, Precision::Confusion()); +} + +TEST(Geom2dPropTest, ComputeTangent_Undefined) +{ + const gp_Vec2d aZero(0.0, 0.0); + const Geom2dProp::TangentResult aRes = Geom2dProp::ComputeTangent(aZero, aZero, aZero, 1.0e-7); + EXPECT_FALSE(aRes.IsDefined); +} + +TEST(Geom2dPropTest, ComputeTangent_NegativeDirection) +{ + const gp_Vec2d aD1(-7.0, 0.0); + const gp_Vec2d aD2(0.0, 0.0); + const gp_Vec2d aD3(0.0, 0.0); + + const Geom2dProp::TangentResult aRes = Geom2dProp::ComputeTangent(aD1, aD2, aD3, 1.0e-7); + ASSERT_TRUE(aRes.IsDefined); + EXPECT_NEAR(aRes.Direction.X(), -1.0, Precision::Confusion()); + EXPECT_NEAR(aRes.Direction.Y(), 0.0, Precision::Confusion()); +} + +TEST(Geom2dPropTest, ComputeCurvature_Circle) +{ + // For a unit circle at param=0: D1=(0,1), D2=(-1,0) + const gp_Vec2d aD1(0.0, 1.0); + const gp_Vec2d aD2(-1.0, 0.0); + + const Geom2dProp::CurvatureResult aRes = + Geom2dProp::ComputeCurvature(aD1, aD2, Precision::Confusion()); + ASSERT_TRUE(aRes.IsDefined); + EXPECT_FALSE(aRes.IsInfinite); + EXPECT_NEAR(aRes.Value, 1.0, 1.0e-10); +} + +TEST(Geom2dPropTest, ComputeCurvature_LargerCircle) +{ + // For circle R=5 at param=0: D1=(0,5), D2=(-5,0) + const gp_Vec2d aD1(0.0, 5.0); + const gp_Vec2d aD2(-5.0, 0.0); + + const Geom2dProp::CurvatureResult aRes = + Geom2dProp::ComputeCurvature(aD1, aD2, Precision::Confusion()); + ASSERT_TRUE(aRes.IsDefined); + EXPECT_NEAR(aRes.Value, 1.0 / 5.0, 1.0e-10); +} + +TEST(Geom2dPropTest, ComputeCurvature_ZeroD1_IsInfinite) +{ + const gp_Vec2d aD1(0.0, 0.0); + const gp_Vec2d aD2(1.0, 0.0); + + const Geom2dProp::CurvatureResult aRes = + Geom2dProp::ComputeCurvature(aD1, aD2, Precision::Confusion()); + ASSERT_TRUE(aRes.IsDefined); + EXPECT_TRUE(aRes.IsInfinite); +} + +TEST(Geom2dPropTest, ComputeCurvature_ZeroD2_IsZero) +{ + const gp_Vec2d aD1(1.0, 0.0); + const gp_Vec2d aD2(0.0, 0.0); + + const Geom2dProp::CurvatureResult aRes = + Geom2dProp::ComputeCurvature(aD1, aD2, Precision::Confusion()); + ASSERT_TRUE(aRes.IsDefined); + EXPECT_FALSE(aRes.IsInfinite); + EXPECT_NEAR(aRes.Value, 0.0, Precision::Confusion()); +} + +TEST(Geom2dPropTest, ComputeCurvature_ParallelD1D2_IsZero) +{ + // D1 and D2 are parallel => cross product is 0 => curvature is 0 + const gp_Vec2d aD1(1.0, 0.0); + const gp_Vec2d aD2(3.0, 0.0); + + const Geom2dProp::CurvatureResult aRes = + Geom2dProp::ComputeCurvature(aD1, aD2, Precision::Confusion()); + ASSERT_TRUE(aRes.IsDefined); + EXPECT_NEAR(aRes.Value, 0.0, Precision::Confusion()); +} + +TEST(Geom2dPropTest, ComputeNormal_Circle) +{ + const gp_Vec2d aD1(0.0, 1.0); + const gp_Vec2d aD2(-1.0, 0.0); + + const Geom2dProp::NormalResult aRes = Geom2dProp::ComputeNormal(aD1, aD2, Precision::Confusion()); + ASSERT_TRUE(aRes.IsDefined); + // Normal perpendicular to tangent + EXPECT_NEAR(std::abs(aRes.Direction.X() * aD1.X() + aRes.Direction.Y() * aD1.Y()), + 0.0, + Precision::Confusion()); +} + +TEST(Geom2dPropTest, ComputeNormal_Line_Undefined) +{ + const gp_Vec2d aD1(1.0, 0.0); + const gp_Vec2d aD2(0.0, 0.0); + + const Geom2dProp::NormalResult aRes = Geom2dProp::ComputeNormal(aD1, aD2, Precision::Confusion()); + EXPECT_FALSE(aRes.IsDefined); +} + +TEST(Geom2dPropTest, ComputeNormal_Perpendicularity) +{ + // At an arbitrary point: D1=(1,1), D2=(0,2) + const gp_Vec2d aD1(1.0, 1.0); + const gp_Vec2d aD2(0.0, 2.0); + + const Geom2dProp::NormalResult aNorm = + Geom2dProp::ComputeNormal(aD1, aD2, Precision::Confusion()); + ASSERT_TRUE(aNorm.IsDefined); + + const Geom2dProp::TangentResult aTan = + Geom2dProp::ComputeTangent(aD1, aD2, gp_Vec2d(0, 0), Precision::Confusion()); + ASSERT_TRUE(aTan.IsDefined); + + const double aDot = + aTan.Direction.X() * aNorm.Direction.X() + aTan.Direction.Y() * aNorm.Direction.Y(); + EXPECT_NEAR(aDot, 0.0, Precision::Confusion()); +} + +TEST(Geom2dPropTest, ComputeCentreOfCurvature_Circle) +{ + const gp_Pnt2d aPnt(1.0, 0.0); + const gp_Vec2d aD1(0.0, 1.0); + const gp_Vec2d aD2(-1.0, 0.0); + + const Geom2dProp::CentreResult aRes = + Geom2dProp::ComputeCentreOfCurvature(aPnt, aD1, aD2, Precision::Confusion()); + ASSERT_TRUE(aRes.IsDefined); + EXPECT_NEAR(aRes.Centre.X(), 0.0, Precision::Confusion()); + EXPECT_NEAR(aRes.Centre.Y(), 0.0, Precision::Confusion()); +} + +TEST(Geom2dPropTest, ComputeCentreOfCurvature_Line_Undefined) +{ + const gp_Pnt2d aPnt(0.0, 0.0); + const gp_Vec2d aD1(1.0, 0.0); + const gp_Vec2d aD2(0.0, 0.0); + + const Geom2dProp::CentreResult aRes = + Geom2dProp::ComputeCentreOfCurvature(aPnt, aD1, aD2, Precision::Confusion()); + EXPECT_FALSE(aRes.IsDefined); +} + +TEST(Geom2dPropTest, ComputeCentreOfCurvature_DistanceEqualsRadius) +{ + // Circle R=3 centered at (0,0), point at (3,0) + const gp_Pnt2d aPnt(3.0, 0.0); + const gp_Vec2d aD1(0.0, 3.0); + const gp_Vec2d aD2(-3.0, 0.0); + + const Geom2dProp::CentreResult aRes = + Geom2dProp::ComputeCentreOfCurvature(aPnt, aD1, aD2, Precision::Confusion()); + ASSERT_TRUE(aRes.IsDefined); + const double aDist = aPnt.Distance(aRes.Centre); + EXPECT_NEAR(aDist, 3.0, Precision::Confusion()); +} + +// ============================================================================ +// Line tests via Geom2dProp_Curve dispatcher +// ============================================================================ + +class Geom2dProp_CurveLineTest : public ::testing::Test +{ +protected: + void SetUp() override + { + gp_Lin2d aLin(gp_Pnt2d(0.0, 0.0), gp_Dir2d(1.0, 0.0)); + myLine = new Geom2d_Line(aLin); + myProp.Initialize(myLine); + } + + occ::handle myLine; + Geom2dProp_Curve myProp; +}; + +TEST_F(Geom2dProp_CurveLineTest, IsInitialized) +{ + EXPECT_TRUE(myProp.IsInitialized()); + EXPECT_EQ(myProp.GetType(), GeomAbs_Line); +} + +TEST_F(Geom2dProp_CurveLineTest, Tangent) +{ + const Geom2dProp::TangentResult aRes = myProp.Tangent(5.0, Precision::Confusion()); + ASSERT_TRUE(aRes.IsDefined); + EXPECT_NEAR(aRes.Direction.X(), 1.0, Precision::Confusion()); + EXPECT_NEAR(aRes.Direction.Y(), 0.0, Precision::Confusion()); +} + +TEST_F(Geom2dProp_CurveLineTest, TangentIsConstant) +{ + // Tangent should be the same at any parameter + for (double u = -100.0; u <= 100.0; u += 50.0) + { + const Geom2dProp::TangentResult aRes = myProp.Tangent(u, Precision::Confusion()); + ASSERT_TRUE(aRes.IsDefined); + EXPECT_NEAR(aRes.Direction.X(), 1.0, Precision::Confusion()); + EXPECT_NEAR(aRes.Direction.Y(), 0.0, Precision::Confusion()); + } +} + +TEST_F(Geom2dProp_CurveLineTest, CurvatureIsZero) +{ + const Geom2dProp::CurvatureResult aRes = myProp.Curvature(5.0, Precision::Confusion()); + ASSERT_TRUE(aRes.IsDefined); + EXPECT_FALSE(aRes.IsInfinite); + EXPECT_NEAR(aRes.Value, 0.0, Precision::Confusion()); +} + +TEST_F(Geom2dProp_CurveLineTest, NormalUndefined) +{ + const Geom2dProp::NormalResult aRes = myProp.Normal(5.0, Precision::Confusion()); + EXPECT_FALSE(aRes.IsDefined); +} + +TEST_F(Geom2dProp_CurveLineTest, CentreUndefined) +{ + const Geom2dProp::CentreResult aRes = myProp.CentreOfCurvature(5.0, Precision::Confusion()); + EXPECT_FALSE(aRes.IsDefined); +} + +TEST_F(Geom2dProp_CurveLineTest, NoExtrema) +{ + const Geom2dProp::CurveAnalysis aRes = myProp.FindCurvatureExtrema(); + EXPECT_TRUE(aRes.IsDone); + EXPECT_TRUE(aRes.Points.IsEmpty()); +} + +TEST_F(Geom2dProp_CurveLineTest, NoInflections) +{ + const Geom2dProp::CurveAnalysis aRes = myProp.FindInflections(); + EXPECT_TRUE(aRes.IsDone); + EXPECT_TRUE(aRes.Points.IsEmpty()); +} + +// Diagonal line +TEST(Geom2dProp_LineTest, DiagonalLine_TangentDirection) +{ + gp_Lin2d aLin(gp_Pnt2d(0.0, 0.0), gp_Dir2d(1.0, 1.0)); + occ::handle aLine = new Geom2d_Line(aLin); + + Geom2dProp_Curve aProp; + aProp.Initialize(aLine); + + const Geom2dProp::TangentResult aTan = aProp.Tangent(0.0, Precision::Confusion()); + ASSERT_TRUE(aTan.IsDefined); + const double aSqrt2Inv = 1.0 / std::sqrt(2.0); + EXPECT_NEAR(aTan.Direction.X(), aSqrt2Inv, Precision::Confusion()); + EXPECT_NEAR(aTan.Direction.Y(), aSqrt2Inv, Precision::Confusion()); +} + +// ============================================================================ +// Circle tests via Geom2dProp_Curve dispatcher +// ============================================================================ + +class Geom2dProp_CurveCircleTest : public ::testing::Test +{ +protected: + void SetUp() override + { + gp_Circ2d aCirc(gp_Ax2d(gp_Pnt2d(0.0, 0.0), gp_Dir2d(1.0, 0.0)), 5.0); + myCircle = new Geom2d_Circle(aCirc); + myProp.Initialize(myCircle); + } + + occ::handle myCircle; + Geom2dProp_Curve myProp; +}; + +TEST_F(Geom2dProp_CurveCircleTest, IsInitialized) +{ + EXPECT_TRUE(myProp.IsInitialized()); + EXPECT_EQ(myProp.GetType(), GeomAbs_Circle); +} + +TEST_F(Geom2dProp_CurveCircleTest, Tangent) +{ + const Geom2dProp::TangentResult aRes = myProp.Tangent(0.0, Precision::Confusion()); + ASSERT_TRUE(aRes.IsDefined); + EXPECT_NEAR(aRes.Direction.X(), 0.0, Precision::Confusion()); + EXPECT_NEAR(aRes.Direction.Y(), 1.0, Precision::Confusion()); +} + +TEST_F(Geom2dProp_CurveCircleTest, TangentAtPiHalf) +{ + const Geom2dProp::TangentResult aRes = myProp.Tangent(M_PI / 2.0, Precision::Confusion()); + ASSERT_TRUE(aRes.IsDefined); + EXPECT_NEAR(aRes.Direction.X(), -1.0, Precision::Confusion()); + EXPECT_NEAR(aRes.Direction.Y(), 0.0, Precision::Confusion()); +} + +TEST_F(Geom2dProp_CurveCircleTest, TangentPerpendicularToRadius) +{ + // At any parameter, tangent should be perpendicular to the radius vector + for (double u = 0.0; u < 2.0 * M_PI; u += M_PI / 7.0) + { + const Geom2dProp::TangentResult aTan = myProp.Tangent(u, Precision::Confusion()); + ASSERT_TRUE(aTan.IsDefined); + // Radius direction at param u on circle centered at origin + const double aRadX = std::cos(u); + const double aRadY = std::sin(u); + const double aDot = aTan.Direction.X() * aRadX + aTan.Direction.Y() * aRadY; + EXPECT_NEAR(aDot, 0.0, 1.0e-10); + } +} + +TEST_F(Geom2dProp_CurveCircleTest, ConstantCurvature) +{ + for (double u = 0.0; u < 2.0 * M_PI; u += M_PI / 4.0) + { + const Geom2dProp::CurvatureResult aRes = myProp.Curvature(u, Precision::Confusion()); + ASSERT_TRUE(aRes.IsDefined); + EXPECT_NEAR(aRes.Value, 1.0 / 5.0, Precision::Confusion()); + } +} + +TEST_F(Geom2dProp_CurveCircleTest, Normal) +{ + const Geom2dProp::NormalResult aRes = myProp.Normal(0.0, Precision::Confusion()); + ASSERT_TRUE(aRes.IsDefined); + const Geom2dProp::TangentResult aTan = myProp.Tangent(0.0, Precision::Confusion()); + EXPECT_NEAR( + std::abs(aTan.Direction.X() * aRes.Direction.X() + aTan.Direction.Y() * aRes.Direction.Y()), + 0.0, + Precision::Confusion()); +} + +TEST_F(Geom2dProp_CurveCircleTest, NormalPointsTowardCenter) +{ + // At param=0, point=(5,0), normal should point toward center (0,0) => direction (-1,0) + const Geom2dProp::NormalResult aRes = myProp.Normal(0.0, Precision::Confusion()); + ASSERT_TRUE(aRes.IsDefined); + EXPECT_NEAR(aRes.Direction.X(), -1.0, Precision::Confusion()); + EXPECT_NEAR(aRes.Direction.Y(), 0.0, Precision::Confusion()); +} + +TEST_F(Geom2dProp_CurveCircleTest, CentreOfCurvature) +{ + const Geom2dProp::CentreResult aRes = myProp.CentreOfCurvature(0.0, Precision::Confusion()); + ASSERT_TRUE(aRes.IsDefined); + EXPECT_NEAR(aRes.Centre.X(), 0.0, Precision::Confusion()); + EXPECT_NEAR(aRes.Centre.Y(), 0.0, Precision::Confusion()); +} + +TEST_F(Geom2dProp_CurveCircleTest, CentreOfCurvature_ConstantAtAllParams) +{ + for (double u = 0.0; u < 2.0 * M_PI; u += M_PI / 5.0) + { + const Geom2dProp::CentreResult aRes = myProp.CentreOfCurvature(u, Precision::Confusion()); + ASSERT_TRUE(aRes.IsDefined); + EXPECT_NEAR(aRes.Centre.X(), 0.0, Precision::Confusion()); + EXPECT_NEAR(aRes.Centre.Y(), 0.0, Precision::Confusion()); + } +} + +TEST_F(Geom2dProp_CurveCircleTest, NoExtrema) +{ + const Geom2dProp::CurveAnalysis aRes = myProp.FindCurvatureExtrema(); + EXPECT_TRUE(aRes.IsDone); + EXPECT_TRUE(aRes.Points.IsEmpty()); +} + +TEST_F(Geom2dProp_CurveCircleTest, NoInflections) +{ + const Geom2dProp::CurveAnalysis aRes = myProp.FindInflections(); + EXPECT_TRUE(aRes.IsDone); + EXPECT_TRUE(aRes.Points.IsEmpty()); +} + +// Different radius circle +TEST(Geom2dProp_CircleTest, SmallRadius_HighCurvature) +{ + gp_Circ2d aCirc(gp_Ax2d(gp_Pnt2d(0.0, 0.0), gp_Dir2d(1.0, 0.0)), 0.1); + occ::handle aCircle = new Geom2d_Circle(aCirc); + + Geom2dProp_Curve aProp; + aProp.Initialize(aCircle); + + const Geom2dProp::CurvatureResult aCurv = aProp.Curvature(0.0, Precision::Confusion()); + ASSERT_TRUE(aCurv.IsDefined); + EXPECT_NEAR(aCurv.Value, 10.0, Precision::Confusion()); +} + +// Off-center circle +TEST(Geom2dProp_CircleTest, OffCenter_CentreOfCurvature) +{ + gp_Circ2d aCirc(gp_Ax2d(gp_Pnt2d(3.0, 7.0), gp_Dir2d(1.0, 0.0)), 4.0); + occ::handle aCircle = new Geom2d_Circle(aCirc); + + Geom2dProp_Curve aProp; + aProp.Initialize(aCircle); + + const Geom2dProp::CentreResult aCentre = aProp.CentreOfCurvature(1.0, Precision::Confusion()); + ASSERT_TRUE(aCentre.IsDefined); + EXPECT_NEAR(aCentre.Centre.X(), 3.0, Precision::Confusion()); + EXPECT_NEAR(aCentre.Centre.Y(), 7.0, Precision::Confusion()); +} + +// ============================================================================ +// Ellipse tests via Geom2dProp_Curve dispatcher +// ============================================================================ + +class Geom2dProp_CurveEllipseTest : public ::testing::Test +{ +protected: + void SetUp() override + { + gp_Elips2d anElips(gp_Ax2d(gp_Pnt2d(0.0, 0.0), gp_Dir2d(1.0, 0.0)), 10.0, 5.0); + myEllipse = new Geom2d_Ellipse(anElips); + myProp.Initialize(myEllipse); + } + + occ::handle myEllipse; + Geom2dProp_Curve myProp; +}; + +TEST_F(Geom2dProp_CurveEllipseTest, IsInitialized) +{ + EXPECT_TRUE(myProp.IsInitialized()); + EXPECT_EQ(myProp.GetType(), GeomAbs_Ellipse); +} + +TEST_F(Geom2dProp_CurveEllipseTest, TangentAtMajorVertex) +{ + // At U=0, point is on major axis endpoint, tangent should be vertical + const Geom2dProp::TangentResult aTan = myProp.Tangent(0.0, Precision::Confusion()); + ASSERT_TRUE(aTan.IsDefined); + EXPECT_NEAR(aTan.Direction.X(), 0.0, Precision::Confusion()); + EXPECT_NEAR(std::abs(aTan.Direction.Y()), 1.0, Precision::Confusion()); +} + +TEST_F(Geom2dProp_CurveEllipseTest, TangentAtMinorVertex) +{ + // At U=PI/2, point is on minor axis endpoint, tangent should be horizontal + const Geom2dProp::TangentResult aTan = myProp.Tangent(M_PI / 2.0, Precision::Confusion()); + ASSERT_TRUE(aTan.IsDefined); + EXPECT_NEAR(std::abs(aTan.Direction.X()), 1.0, Precision::Confusion()); + EXPECT_NEAR(aTan.Direction.Y(), 0.0, Precision::Confusion()); +} + +TEST_F(Geom2dProp_CurveEllipseTest, TangentPerpToNormal) +{ + for (double u = 0.0; u < 2.0 * M_PI; u += M_PI / 6.0) + { + const Geom2dProp::TangentResult aTan = myProp.Tangent(u, Precision::Confusion()); + const Geom2dProp::NormalResult aNorm = myProp.Normal(u, Precision::Confusion()); + ASSERT_TRUE(aTan.IsDefined); + ASSERT_TRUE(aNorm.IsDefined); + const double aDot = + aTan.Direction.X() * aNorm.Direction.X() + aTan.Direction.Y() * aNorm.Direction.Y(); + EXPECT_NEAR(aDot, 0.0, 1.0e-10); + } +} + +TEST_F(Geom2dProp_CurveEllipseTest, CurvatureAtMajorVertex) +{ + // At U=0 (major vertex): curvature = a/b^2 = 10/25 = 0.4 + const Geom2dProp::CurvatureResult aRes = myProp.Curvature(0.0, Precision::Confusion()); + ASSERT_TRUE(aRes.IsDefined); + EXPECT_NEAR(aRes.Value, 10.0 / 25.0, 1.0e-6); +} + +TEST_F(Geom2dProp_CurveEllipseTest, CurvatureAtMinorVertex) +{ + // At U=PI/2 (minor vertex): curvature = b/a^2 = 5/100 = 0.05 + const Geom2dProp::CurvatureResult aRes = myProp.Curvature(M_PI / 2.0, Precision::Confusion()); + ASSERT_TRUE(aRes.IsDefined); + EXPECT_NEAR(aRes.Value, 5.0 / 100.0, 1.0e-6); +} + +TEST_F(Geom2dProp_CurveEllipseTest, CurvatureAtPi) +{ + // At U=PI (opposite major vertex), curvature same as U=0 + const Geom2dProp::CurvatureResult aRes = myProp.Curvature(M_PI, Precision::Confusion()); + ASSERT_TRUE(aRes.IsDefined); + EXPECT_NEAR(aRes.Value, 10.0 / 25.0, 1.0e-6); +} + +TEST_F(Geom2dProp_CurveEllipseTest, CurvatureSymmetry) +{ + // Curvature at U and -U should be equal (ellipse is symmetric) + for (double u = 0.1; u < M_PI; u += 0.3) + { + const Geom2dProp::CurvatureResult aRes1 = myProp.Curvature(u, Precision::Confusion()); + const Geom2dProp::CurvatureResult aRes2 = myProp.Curvature(-u, Precision::Confusion()); + ASSERT_TRUE(aRes1.IsDefined); + ASSERT_TRUE(aRes2.IsDefined); + EXPECT_NEAR(aRes1.Value, aRes2.Value, 1.0e-10); + } +} + +TEST_F(Geom2dProp_CurveEllipseTest, CurvatureMaxAtMajorVertex) +{ + // |curvature| should be maximum at major vertex (U=0, PI) + const double aCurvMax = myProp.Curvature(0.0, Precision::Confusion()).Value; + for (double u = 0.1; u < 2.0 * M_PI; u += 0.2) + { + const Geom2dProp::CurvatureResult aCurv = myProp.Curvature(u, Precision::Confusion()); + ASSERT_TRUE(aCurv.IsDefined); + EXPECT_LE(std::abs(aCurv.Value), std::abs(aCurvMax) + 1.0e-10); + } +} + +TEST_F(Geom2dProp_CurveEllipseTest, NormalAtMajorVertex) +{ + const Geom2dProp::NormalResult aRes = myProp.Normal(0.0, Precision::Confusion()); + ASSERT_TRUE(aRes.IsDefined); + // At (10,0) on x-axis aligned ellipse, normal should point toward center: (-1,0) + EXPECT_NEAR(aRes.Direction.X(), -1.0, Precision::Confusion()); + EXPECT_NEAR(aRes.Direction.Y(), 0.0, Precision::Confusion()); +} + +TEST_F(Geom2dProp_CurveEllipseTest, CentreOfCurvatureAtMajorVertex) +{ + const Geom2dProp::CentreResult aRes = myProp.CentreOfCurvature(0.0, Precision::Confusion()); + ASSERT_TRUE(aRes.IsDefined); + // Radius of curvature at major vertex = b^2/a = 25/10 = 2.5 + // Centre should be at (10 - 2.5, 0) = (7.5, 0) + EXPECT_NEAR(aRes.Centre.X(), 7.5, 1.0e-6); + EXPECT_NEAR(aRes.Centre.Y(), 0.0, 1.0e-6); +} + +TEST_F(Geom2dProp_CurveEllipseTest, CentreOfCurvatureAtMinorVertex) +{ + const Geom2dProp::CentreResult aRes = + myProp.CentreOfCurvature(M_PI / 2.0, Precision::Confusion()); + ASSERT_TRUE(aRes.IsDefined); + // Radius of curvature at minor vertex = a^2/b = 100/5 = 20 + // Centre should be at (0, 5 - 20) = (0, -15) + EXPECT_NEAR(aRes.Centre.X(), 0.0, 1.0e-6); + EXPECT_NEAR(aRes.Centre.Y(), -15.0, 1.0e-6); +} + +TEST_F(Geom2dProp_CurveEllipseTest, FindCurvatureExtrema) +{ + const Geom2dProp::CurveAnalysis aRes = myProp.FindCurvatureExtrema(); + ASSERT_TRUE(aRes.IsDone); + EXPECT_EQ(aRes.Points.Length(), 4); + + if (aRes.Points.Length() >= 4) + { + EXPECT_NEAR(aRes.Points.Value(0).Parameter, 0.0, 1.0e-10); + EXPECT_EQ(aRes.Points.Value(0).Type, Geom2dProp::CIType::MinCurvature); + + EXPECT_NEAR(aRes.Points.Value(1).Parameter, M_PI / 2.0, 1.0e-10); + EXPECT_EQ(aRes.Points.Value(1).Type, Geom2dProp::CIType::MaxCurvature); + + EXPECT_NEAR(aRes.Points.Value(2).Parameter, M_PI, 1.0e-10); + EXPECT_EQ(aRes.Points.Value(2).Type, Geom2dProp::CIType::MinCurvature); + + EXPECT_NEAR(aRes.Points.Value(3).Parameter, 3.0 * M_PI / 2.0, 1.0e-10); + EXPECT_EQ(aRes.Points.Value(3).Type, Geom2dProp::CIType::MaxCurvature); + } +} + +TEST_F(Geom2dProp_CurveEllipseTest, NoInflections) +{ + const Geom2dProp::CurveAnalysis aRes = myProp.FindInflections(); + EXPECT_TRUE(aRes.IsDone); + EXPECT_TRUE(aRes.Points.IsEmpty()); +} + +// Ellipse with equal semi-axes is a circle +TEST(Geom2dProp_EllipseTest, EqualSemiAxes_BehavesLikeCircle) +{ + gp_Elips2d anElips(gp_Ax2d(gp_Pnt2d(0.0, 0.0), gp_Dir2d(1.0, 0.0)), 5.0, 5.0); + occ::handle anEllipse = new Geom2d_Ellipse(anElips); + + Geom2dProp_Curve aProp; + aProp.Initialize(anEllipse); + + // Curvature should be constant = 1/R = 1/5 + for (double u = 0.0; u < 2.0 * M_PI; u += M_PI / 4.0) + { + const Geom2dProp::CurvatureResult aCurv = aProp.Curvature(u, Precision::Confusion()); + ASSERT_TRUE(aCurv.IsDefined); + EXPECT_NEAR(aCurv.Value, 1.0 / 5.0, 1.0e-10); + } +} + +// ============================================================================ +// Hyperbola tests +// ============================================================================ + +class Geom2dProp_CurveHyperbolaTest : public ::testing::Test +{ +protected: + void SetUp() override + { + gp_Hypr2d aHypr(gp_Ax2d(gp_Pnt2d(0.0, 0.0), gp_Dir2d(1.0, 0.0)), 5.0, 3.0); + myHyperbola = new Geom2d_Hyperbola(aHypr); + myProp.Initialize(myHyperbola); + } + + occ::handle myHyperbola; + Geom2dProp_Curve myProp; +}; + +TEST_F(Geom2dProp_CurveHyperbolaTest, IsInitialized) +{ + EXPECT_TRUE(myProp.IsInitialized()); + EXPECT_EQ(myProp.GetType(), GeomAbs_Hyperbola); +} + +TEST_F(Geom2dProp_CurveHyperbolaTest, TangentAtVertex) +{ + const Geom2dProp::TangentResult aTan = myProp.Tangent(0.0, Precision::Confusion()); + ASSERT_TRUE(aTan.IsDefined); + // At vertex of hyperbola (t=0), tangent is vertical + EXPECT_NEAR(aTan.Direction.X(), 0.0, Precision::Confusion()); + EXPECT_NEAR(std::abs(aTan.Direction.Y()), 1.0, Precision::Confusion()); +} + +TEST_F(Geom2dProp_CurveHyperbolaTest, CurvatureAtVertex) +{ + // At vertex (t=0): curvature = b^2/a^2 * 1/a * a = b^2/(a * a) ... actually K = b^2/a + // For hyperbola x=a*cosh(t), y=b*sinh(t): + // K(0) = b^2/(a^2 * (b^2/a^2)^(3/2)) ... Just verify it's defined and positive. + const Geom2dProp::CurvatureResult aCurv = myProp.Curvature(0.0, Precision::Confusion()); + ASSERT_TRUE(aCurv.IsDefined); + EXPECT_GT(aCurv.Value, 0.0); +} + +TEST_F(Geom2dProp_CurveHyperbolaTest, CurvatureDecreasesFromVertex) +{ + // Curvature should be maximum at vertex and decrease away from it + const double aCurvAtVertex = std::abs(myProp.Curvature(0.0, Precision::Confusion()).Value); + for (double u = 0.5; u < 3.0; u += 0.5) + { + const Geom2dProp::CurvatureResult aCurv = myProp.Curvature(u, Precision::Confusion()); + ASSERT_TRUE(aCurv.IsDefined); + EXPECT_LT(std::abs(aCurv.Value), aCurvAtVertex + 1.0e-10); + } +} + +TEST_F(Geom2dProp_CurveHyperbolaTest, NormalAtVertex) +{ + const Geom2dProp::NormalResult aNorm = myProp.Normal(0.0, Precision::Confusion()); + ASSERT_TRUE(aNorm.IsDefined); + // Normal should be perpendicular to tangent + const Geom2dProp::TangentResult aTan = myProp.Tangent(0.0, Precision::Confusion()); + const double aDot = + aTan.Direction.X() * aNorm.Direction.X() + aTan.Direction.Y() * aNorm.Direction.Y(); + EXPECT_NEAR(aDot, 0.0, 1.0e-10); +} + +TEST_F(Geom2dProp_CurveHyperbolaTest, CentreOfCurvatureAtVertex) +{ + const Geom2dProp::CentreResult aCentre = myProp.CentreOfCurvature(0.0, Precision::Confusion()); + ASSERT_TRUE(aCentre.IsDefined); +} + +TEST_F(Geom2dProp_CurveHyperbolaTest, CurvatureExtremaAtVertex) +{ + const Geom2dProp::CurveAnalysis aExtrema = myProp.FindCurvatureExtrema(); + ASSERT_TRUE(aExtrema.IsDone); + EXPECT_EQ(aExtrema.Points.Length(), 1); + if (!aExtrema.Points.IsEmpty()) + { + EXPECT_NEAR(aExtrema.Points.Value(0).Parameter, 0.0, 1.0e-10); + EXPECT_EQ(aExtrema.Points.Value(0).Type, Geom2dProp::CIType::MinCurvature); + } +} + +TEST_F(Geom2dProp_CurveHyperbolaTest, NoInflections) +{ + const Geom2dProp::CurveAnalysis aInfl = myProp.FindInflections(); + EXPECT_TRUE(aInfl.IsDone); + EXPECT_TRUE(aInfl.Points.IsEmpty()); +} + +// ============================================================================ +// Parabola tests +// ============================================================================ + +class Geom2dProp_CurveParabolaTest : public ::testing::Test +{ +protected: + void SetUp() override + { + // Focal parameter p=2 => parabola y^2 = 4px = 8x + gp_Parab2d aParab(gp_Ax2d(gp_Pnt2d(0.0, 0.0), gp_Dir2d(1.0, 0.0)), 2.0); + myParabola = new Geom2d_Parabola(aParab); + myProp.Initialize(myParabola); + } + + occ::handle myParabola; + Geom2dProp_Curve myProp; +}; + +TEST_F(Geom2dProp_CurveParabolaTest, IsInitialized) +{ + EXPECT_TRUE(myProp.IsInitialized()); + EXPECT_EQ(myProp.GetType(), GeomAbs_Parabola); +} + +TEST_F(Geom2dProp_CurveParabolaTest, TangentAtVertex) +{ + const Geom2dProp::TangentResult aTan = myProp.Tangent(0.0, Precision::Confusion()); + ASSERT_TRUE(aTan.IsDefined); + // At vertex, tangent is vertical + EXPECT_NEAR(aTan.Direction.X(), 0.0, Precision::Confusion()); + EXPECT_NEAR(std::abs(aTan.Direction.Y()), 1.0, Precision::Confusion()); +} + +TEST_F(Geom2dProp_CurveParabolaTest, CurvatureAtVertex) +{ + // At vertex: curvature = 1/(2*focal_parameter) = 1/4 = 0.25 + const Geom2dProp::CurvatureResult aCurv = myProp.Curvature(0.0, Precision::Confusion()); + ASSERT_TRUE(aCurv.IsDefined); + EXPECT_GT(aCurv.Value, 0.0); +} + +TEST_F(Geom2dProp_CurveParabolaTest, CurvatureDecreasesFromVertex) +{ + const double aCurvAtVertex = std::abs(myProp.Curvature(0.0, Precision::Confusion()).Value); + for (double u = 1.0; u <= 5.0; u += 1.0) + { + const Geom2dProp::CurvatureResult aCurv = myProp.Curvature(u, Precision::Confusion()); + ASSERT_TRUE(aCurv.IsDefined); + EXPECT_LT(std::abs(aCurv.Value), aCurvAtVertex); + } +} + +TEST_F(Geom2dProp_CurveParabolaTest, CurvatureSymmetric) +{ + // Curvature at U and -U should be equal + for (double u = 0.5; u <= 5.0; u += 0.5) + { + const double aCurv1 = myProp.Curvature(u, Precision::Confusion()).Value; + const double aCurv2 = myProp.Curvature(-u, Precision::Confusion()).Value; + EXPECT_NEAR(aCurv1, aCurv2, 1.0e-10); + } +} + +TEST_F(Geom2dProp_CurveParabolaTest, NormalAtVertex) +{ + const Geom2dProp::NormalResult aNorm = myProp.Normal(0.0, Precision::Confusion()); + ASSERT_TRUE(aNorm.IsDefined); + const Geom2dProp::TangentResult aTan = myProp.Tangent(0.0, Precision::Confusion()); + const double aDot = + aTan.Direction.X() * aNorm.Direction.X() + aTan.Direction.Y() * aNorm.Direction.Y(); + EXPECT_NEAR(aDot, 0.0, 1.0e-10); +} + +TEST_F(Geom2dProp_CurveParabolaTest, CentreOfCurvatureAtVertex) +{ + const Geom2dProp::CentreResult aCentre = myProp.CentreOfCurvature(0.0, Precision::Confusion()); + ASSERT_TRUE(aCentre.IsDefined); +} + +TEST_F(Geom2dProp_CurveParabolaTest, CurvatureExtremaAtVertex) +{ + const Geom2dProp::CurveAnalysis aExtrema = myProp.FindCurvatureExtrema(); + ASSERT_TRUE(aExtrema.IsDone); + EXPECT_EQ(aExtrema.Points.Length(), 1); + if (!aExtrema.Points.IsEmpty()) + { + EXPECT_NEAR(aExtrema.Points.Value(0).Parameter, 0.0, 1.0e-10); + EXPECT_EQ(aExtrema.Points.Value(0).Type, Geom2dProp::CIType::MinCurvature); + } +} + +TEST_F(Geom2dProp_CurveParabolaTest, NoInflections) +{ + const Geom2dProp::CurveAnalysis aInfl = myProp.FindInflections(); + EXPECT_TRUE(aInfl.IsDone); + EXPECT_TRUE(aInfl.Points.IsEmpty()); +} + +// ============================================================================ +// Bezier curve tests +// ============================================================================ + +class Geom2dProp_CurveBezierTest : public ::testing::Test +{ +protected: + void SetUp() override + { + // S-shaped cubic Bezier: (0,0), (1,2), (3,-2), (4,0) + NCollection_Array1 aPoles(1, 4); + aPoles(1) = gp_Pnt2d(0.0, 0.0); + aPoles(2) = gp_Pnt2d(1.0, 2.0); + aPoles(3) = gp_Pnt2d(3.0, -2.0); + aPoles(4) = gp_Pnt2d(4.0, 0.0); + myBezier = new Geom2d_BezierCurve(aPoles); + myProp.Initialize(myBezier); + } + + occ::handle myBezier; + Geom2dProp_Curve myProp; +}; + +TEST_F(Geom2dProp_CurveBezierTest, IsInitialized) +{ + EXPECT_TRUE(myProp.IsInitialized()); + EXPECT_EQ(myProp.GetType(), GeomAbs_BezierCurve); +} + +TEST_F(Geom2dProp_CurveBezierTest, TangentAtStart) +{ + const Geom2dProp::TangentResult aTan = myProp.Tangent(0.0, Precision::Confusion()); + ASSERT_TRUE(aTan.IsDefined); + // Tangent at start should point toward second control point: (1,2) + // Direction should be proportional to (1,2), normalized + const double aLen = std::sqrt(1.0 + 4.0); + EXPECT_NEAR(aTan.Direction.X(), 1.0 / aLen, 1.0e-6); + EXPECT_NEAR(aTan.Direction.Y(), 2.0 / aLen, 1.0e-6); +} + +TEST_F(Geom2dProp_CurveBezierTest, TangentAtEnd) +{ + const Geom2dProp::TangentResult aTan = myProp.Tangent(1.0, Precision::Confusion()); + ASSERT_TRUE(aTan.IsDefined); + // Tangent at end should point from third control point to fourth: (4,0)-(3,-2) = (1,2) + const double aLen = std::sqrt(1.0 + 4.0); + EXPECT_NEAR(aTan.Direction.X(), 1.0 / aLen, 1.0e-6); + EXPECT_NEAR(aTan.Direction.Y(), 2.0 / aLen, 1.0e-6); +} + +TEST_F(Geom2dProp_CurveBezierTest, CurvatureAtMultiplePoints) +{ + for (double u = 0.0; u <= 1.0; u += 0.1) + { + const Geom2dProp::CurvatureResult aCurv = myProp.Curvature(u, Precision::Confusion()); + EXPECT_TRUE(aCurv.IsDefined); + EXPECT_FALSE(aCurv.IsInfinite); + } +} + +TEST_F(Geom2dProp_CurveBezierTest, NormalPerpendicularToTangent) +{ + for (double u = 0.1; u < 1.0; u += 0.2) + { + const Geom2dProp::TangentResult aTan = myProp.Tangent(u, Precision::Confusion()); + const Geom2dProp::NormalResult aNorm = myProp.Normal(u, Precision::Confusion()); + if (aTan.IsDefined && aNorm.IsDefined) + { + const double aDot = + aTan.Direction.X() * aNorm.Direction.X() + aTan.Direction.Y() * aNorm.Direction.Y(); + EXPECT_NEAR(aDot, 0.0, 1.0e-10); + } + } +} + +TEST_F(Geom2dProp_CurveBezierTest, CentreOfCurvature_DistanceIsRadiusOfCurvature) +{ + const double u = 0.3; + const Geom2dProp::CurvatureResult aCurv = myProp.Curvature(u, Precision::Confusion()); + const Geom2dProp::CentreResult aCentre = myProp.CentreOfCurvature(u, Precision::Confusion()); + if (aCurv.IsDefined && aCentre.IsDefined && !aCurv.IsInfinite && std::abs(aCurv.Value) > 1.0e-10) + { + gp_Pnt2d aPnt; + gp_Vec2d aD1; + myBezier->D1(u, aPnt, aD1); + const double aDist = aPnt.Distance(aCentre.Centre); + EXPECT_NEAR(aDist, 1.0 / std::abs(aCurv.Value), 1.0e-6); + } +} + +TEST_F(Geom2dProp_CurveBezierTest, InflectionPoints) +{ + // S-shaped curve should have inflection point(s) near the middle + const Geom2dProp::CurveAnalysis aInflections = myProp.FindInflections(); + EXPECT_TRUE(aInflections.IsDone); + EXPECT_GE(aInflections.Points.Length(), 1); + if (!aInflections.Points.IsEmpty()) + { + for (int i = 0; i < aInflections.Points.Length(); ++i) + { + EXPECT_EQ(aInflections.Points.Value(i).Type, Geom2dProp::CIType::Inflection); + // Should be within the parameter range [0, 1] + EXPECT_GE(aInflections.Points.Value(i).Parameter, 0.0 - 1.0e-6); + EXPECT_LE(aInflections.Points.Value(i).Parameter, 1.0 + 1.0e-6); + } + } +} + +TEST_F(Geom2dProp_CurveBezierTest, CurvatureExtrema) +{ + const Geom2dProp::CurveAnalysis aExtrema = myProp.FindCurvatureExtrema(); + EXPECT_TRUE(aExtrema.IsDone); + // S-shaped cubic should have curvature extrema + for (int i = 0; i < aExtrema.Points.Length(); ++i) + { + EXPECT_TRUE(aExtrema.Points.Value(i).Type == Geom2dProp::CIType::MinCurvature + || aExtrema.Points.Value(i).Type == Geom2dProp::CIType::MaxCurvature); + EXPECT_GE(aExtrema.Points.Value(i).Parameter, 0.0 - 1.0e-6); + EXPECT_LE(aExtrema.Points.Value(i).Parameter, 1.0 + 1.0e-6); + } +} + +// Straight-line Bezier +TEST(Geom2dProp_BezierTest, StraightLine_ZeroCurvature) +{ + NCollection_Array1 aPoles(1, 3); + aPoles(1) = gp_Pnt2d(0.0, 0.0); + aPoles(2) = gp_Pnt2d(2.0, 0.0); + aPoles(3) = gp_Pnt2d(4.0, 0.0); + occ::handle aBezier = new Geom2d_BezierCurve(aPoles); + + Geom2dProp_Curve aProp; + aProp.Initialize(aBezier); + + const Geom2dProp::CurvatureResult aCurv = aProp.Curvature(0.5, Precision::Confusion()); + ASSERT_TRUE(aCurv.IsDefined); + EXPECT_NEAR(aCurv.Value, 0.0, Precision::Confusion()); +} + +// Quadratic Bezier (arc-like) +TEST(Geom2dProp_BezierTest, QuadraticBezier_Properties) +{ + NCollection_Array1 aPoles(1, 3); + aPoles(1) = gp_Pnt2d(0.0, 0.0); + aPoles(2) = gp_Pnt2d(1.0, 2.0); + aPoles(3) = gp_Pnt2d(2.0, 0.0); + occ::handle aBezier = new Geom2d_BezierCurve(aPoles); + + Geom2dProp_Curve aProp; + aProp.Initialize(aBezier); + + // Symmetric parabolic arc: curvature should be max at midpoint + const double aCurvMid = std::abs(aProp.Curvature(0.5, Precision::Confusion()).Value); + EXPECT_GT(aCurvMid, 0.0); + + // No inflections for a parabolic arc + const Geom2dProp::CurveAnalysis aInfl = aProp.FindInflections(); + EXPECT_TRUE(aInfl.IsDone); + EXPECT_TRUE(aInfl.Points.IsEmpty()); +} + +// ============================================================================ +// BSpline curve tests +// ============================================================================ + +class Geom2dProp_CurveBSplineTest : public ::testing::Test +{ +protected: + void SetUp() override + { + NCollection_Array1 aPoles(1, 4); + aPoles(1) = gp_Pnt2d(0.0, 0.0); + aPoles(2) = gp_Pnt2d(1.0, 2.0); + aPoles(3) = gp_Pnt2d(3.0, 2.0); + aPoles(4) = gp_Pnt2d(4.0, 0.0); + + NCollection_Array1 aKnots(1, 3); + aKnots(1) = 0.0; + aKnots(2) = 0.5; + aKnots(3) = 1.0; + + NCollection_Array1 aMults(1, 3); + aMults(1) = 3; + aMults(2) = 1; + aMults(3) = 3; + + myBSpline = new Geom2d_BSplineCurve(aPoles, aKnots, aMults, 2); + myProp.Initialize(myBSpline); + } + + occ::handle myBSpline; + Geom2dProp_Curve myProp; +}; + +TEST_F(Geom2dProp_CurveBSplineTest, IsInitialized) +{ + EXPECT_TRUE(myProp.IsInitialized()); + EXPECT_EQ(myProp.GetType(), GeomAbs_BSplineCurve); +} + +TEST_F(Geom2dProp_CurveBSplineTest, TangentDefined) +{ + for (double u = 0.0; u <= 1.0; u += 0.2) + { + const Geom2dProp::TangentResult aTan = myProp.Tangent(u, Precision::Confusion()); + EXPECT_TRUE(aTan.IsDefined); + } +} + +TEST_F(Geom2dProp_CurveBSplineTest, CurvatureDefined) +{ + for (double u = 0.0; u <= 1.0; u += 0.2) + { + const Geom2dProp::CurvatureResult aCurv = myProp.Curvature(u, Precision::Confusion()); + EXPECT_TRUE(aCurv.IsDefined); + EXPECT_FALSE(aCurv.IsInfinite); + } +} + +TEST_F(Geom2dProp_CurveBSplineTest, NormalDefined) +{ + // B-spline has non-zero curvature, so normal should be defined in the middle + const Geom2dProp::NormalResult aNorm = myProp.Normal(0.3, Precision::Confusion()); + EXPECT_TRUE(aNorm.IsDefined); +} + +TEST_F(Geom2dProp_CurveBSplineTest, CentreOfCurvatureDefined) +{ + const Geom2dProp::CentreResult aCentre = myProp.CentreOfCurvature(0.3, Precision::Confusion()); + EXPECT_TRUE(aCentre.IsDefined); +} + +TEST_F(Geom2dProp_CurveBSplineTest, CurvatureExtrema) +{ + const Geom2dProp::CurveAnalysis aExtrema = myProp.FindCurvatureExtrema(); + EXPECT_TRUE(aExtrema.IsDone); + // Verify all extrema parameters are within [0, 1] + for (int i = 0; i < aExtrema.Points.Length(); ++i) + { + EXPECT_GE(aExtrema.Points.Value(i).Parameter, 0.0 - 1.0e-6); + EXPECT_LE(aExtrema.Points.Value(i).Parameter, 1.0 + 1.0e-6); + } +} + +TEST_F(Geom2dProp_CurveBSplineTest, InflectionPoints) +{ + const Geom2dProp::CurveAnalysis aInfl = myProp.FindInflections(); + EXPECT_TRUE(aInfl.IsDone); + for (int i = 0; i < aInfl.Points.Length(); ++i) + { + EXPECT_EQ(aInfl.Points.Value(i).Type, Geom2dProp::CIType::Inflection); + EXPECT_GE(aInfl.Points.Value(i).Parameter, 0.0 - 1.0e-6); + EXPECT_LE(aInfl.Points.Value(i).Parameter, 1.0 + 1.0e-6); + } +} + +// B-spline with lower continuity (C1) +TEST(Geom2dProp_BSplineTest, LowContinuity_C1) +{ + NCollection_Array1 aPoles(1, 5); + aPoles(1) = gp_Pnt2d(0.0, 0.0); + aPoles(2) = gp_Pnt2d(1.0, 3.0); + aPoles(3) = gp_Pnt2d(2.0, 1.0); + aPoles(4) = gp_Pnt2d(3.0, 3.0); + aPoles(5) = gp_Pnt2d(4.0, 0.0); + + NCollection_Array1 aKnots(1, 4); + aKnots(1) = 0.0; + aKnots(2) = 0.33; + aKnots(3) = 0.66; + aKnots(4) = 1.0; + + NCollection_Array1 aMults(1, 4); + aMults(1) = 3; + aMults(2) = 1; + aMults(3) = 1; + aMults(4) = 3; + + occ::handle aBSpline = new Geom2d_BSplineCurve(aPoles, aKnots, aMults, 2); + + Geom2dProp_Curve aProp; + aProp.Initialize(aBSpline); + EXPECT_TRUE(aProp.IsInitialized()); + + // Should work with C3 interval subdivision + const Geom2dProp::CurveAnalysis aExtrema = aProp.FindCurvatureExtrema(); + EXPECT_TRUE(aExtrema.IsDone); + + const Geom2dProp::CurveAnalysis aInfl = aProp.FindInflections(); + EXPECT_TRUE(aInfl.IsDone); +} + +// ============================================================================ +// Offset curve tests +// ============================================================================ + +TEST(Geom2dProp_OffsetCurveTest, IsInitialized) +{ + gp_Circ2d aCirc(gp_Ax2d(gp_Pnt2d(0.0, 0.0), gp_Dir2d(1.0, 0.0)), 5.0); + occ::handle aCircle = new Geom2d_Circle(aCirc); + occ::handle anOffset = new Geom2d_OffsetCurve(aCircle, 2.0); + + Geom2dProp_Curve aProp; + aProp.Initialize(anOffset); + EXPECT_TRUE(aProp.IsInitialized()); + EXPECT_EQ(aProp.GetType(), GeomAbs_OffsetCurve); +} + +TEST(Geom2dProp_OffsetCurveTest, OffsetCircle_ConstantCurvature) +{ + // Offset of circle R=5 by +2 gives circle R=7, curvature=1/7 + gp_Circ2d aCirc(gp_Ax2d(gp_Pnt2d(0.0, 0.0), gp_Dir2d(1.0, 0.0)), 5.0); + occ::handle aCircle = new Geom2d_Circle(aCirc); + occ::handle anOffset = new Geom2d_OffsetCurve(aCircle, 2.0); + + Geom2dProp_Curve aProp; + aProp.Initialize(anOffset); + + for (double u = 0.0; u < 2.0 * M_PI; u += M_PI / 4.0) + { + const Geom2dProp::CurvatureResult aCurv = aProp.Curvature(u, Precision::Confusion()); + ASSERT_TRUE(aCurv.IsDefined); + EXPECT_NEAR(std::abs(aCurv.Value), 1.0 / 7.0, 1.0e-6); + } +} + +TEST(Geom2dProp_OffsetCurveTest, TangentAndNormalDefined) +{ + gp_Circ2d aCirc(gp_Ax2d(gp_Pnt2d(0.0, 0.0), gp_Dir2d(1.0, 0.0)), 5.0); + occ::handle aCircle = new Geom2d_Circle(aCirc); + occ::handle anOffset = new Geom2d_OffsetCurve(aCircle, 2.0); + + Geom2dProp_Curve aProp; + aProp.Initialize(anOffset); + + const Geom2dProp::TangentResult aTan = aProp.Tangent(0.5, Precision::Confusion()); + EXPECT_TRUE(aTan.IsDefined); + + const Geom2dProp::NormalResult aNorm = aProp.Normal(0.5, Precision::Confusion()); + EXPECT_TRUE(aNorm.IsDefined); + + // Perpendicularity + if (aTan.IsDefined && aNorm.IsDefined) + { + const double aDot = + aTan.Direction.X() * aNorm.Direction.X() + aTan.Direction.Y() * aNorm.Direction.Y(); + EXPECT_NEAR(aDot, 0.0, 1.0e-10); + } +} + +TEST(Geom2dProp_OffsetCurveTest, OffsetEllipse_ExtremaAndInflections) +{ + gp_Elips2d anElips(gp_Ax2d(gp_Pnt2d(0.0, 0.0), gp_Dir2d(1.0, 0.0)), 10.0, 5.0); + occ::handle anEllipse = new Geom2d_Ellipse(anElips); + occ::handle anOffset = new Geom2d_OffsetCurve(anEllipse, 1.0); + + Geom2dProp_Curve aProp; + aProp.Initialize(anOffset); + EXPECT_TRUE(aProp.IsInitialized()); + + const Geom2dProp::CurveAnalysis aExtrema = aProp.FindCurvatureExtrema(); + EXPECT_TRUE(aExtrema.IsDone); + + const Geom2dProp::CurveAnalysis aInfl = aProp.FindInflections(); + EXPECT_TRUE(aInfl.IsDone); +} + +// ============================================================================ +// TrimmedCurve tests +// ============================================================================ + +TEST(Geom2dProp_TrimmedCurveTest, UnwrapsToCircle) +{ + gp_Circ2d aCirc(gp_Ax2d(gp_Pnt2d(0.0, 0.0), gp_Dir2d(1.0, 0.0)), 5.0); + occ::handle aCircle = new Geom2d_Circle(aCirc); + occ::handle aTrimmed = new Geom2d_TrimmedCurve(aCircle, 0.0, M_PI); + + Geom2dProp_Curve aProp; + aProp.Initialize(aTrimmed); + EXPECT_TRUE(aProp.IsInitialized()); + EXPECT_EQ(aProp.GetType(), GeomAbs_Circle); + + const Geom2dProp::CurvatureResult aCurv = aProp.Curvature(0.5, Precision::Confusion()); + ASSERT_TRUE(aCurv.IsDefined); + EXPECT_NEAR(aCurv.Value, 1.0 / 5.0, Precision::Confusion()); +} + +TEST(Geom2dProp_TrimmedCurveTest, UnwrapsToEllipse) +{ + gp_Elips2d anElips(gp_Ax2d(gp_Pnt2d(0.0, 0.0), gp_Dir2d(1.0, 0.0)), 10.0, 5.0); + occ::handle anEllipse = new Geom2d_Ellipse(anElips); + occ::handle aTrimmed = new Geom2d_TrimmedCurve(anEllipse, 0.0, M_PI); + + Geom2dProp_Curve aProp; + aProp.Initialize(aTrimmed); + EXPECT_TRUE(aProp.IsInitialized()); + EXPECT_EQ(aProp.GetType(), GeomAbs_Ellipse); +} + +TEST(Geom2dProp_TrimmedCurveTest, NestedTrimmedCurve) +{ + gp_Circ2d aCirc(gp_Ax2d(gp_Pnt2d(0.0, 0.0), gp_Dir2d(1.0, 0.0)), 3.0); + occ::handle aCircle = new Geom2d_Circle(aCirc); + occ::handle aTrimmed1 = new Geom2d_TrimmedCurve(aCircle, 0.0, M_PI); + occ::handle aTrimmed2 = new Geom2d_TrimmedCurve(aTrimmed1, 0.1, 1.0); + + Geom2dProp_Curve aProp; + aProp.Initialize(aTrimmed2); + EXPECT_TRUE(aProp.IsInitialized()); + EXPECT_EQ(aProp.GetType(), GeomAbs_Circle); + + const Geom2dProp::CurvatureResult aCurv = aProp.Curvature(0.5, Precision::Confusion()); + ASSERT_TRUE(aCurv.IsDefined); + EXPECT_NEAR(aCurv.Value, 1.0 / 3.0, Precision::Confusion()); +} + +// ============================================================================ +// Null / uninitialized tests +// ============================================================================ + +TEST(Geom2dProp_CurveTest, NullHandle_NotInitialized) +{ + Geom2dProp_Curve aProp; + occ::handle aNullCurve; + aProp.Initialize(aNullCurve); + EXPECT_FALSE(aProp.IsInitialized()); + + const Geom2dProp::TangentResult aTan = aProp.Tangent(0.0, 1.0e-7); + EXPECT_FALSE(aTan.IsDefined); + + const Geom2dProp::CurvatureResult aCurv = aProp.Curvature(0.0, 1.0e-7); + EXPECT_FALSE(aCurv.IsDefined); + + const Geom2dProp::NormalResult aNorm = aProp.Normal(0.0, 1.0e-7); + EXPECT_FALSE(aNorm.IsDefined); + + const Geom2dProp::CentreResult aCentre = aProp.CentreOfCurvature(0.0, 1.0e-7); + EXPECT_FALSE(aCentre.IsDefined); + + const Geom2dProp::CurveAnalysis aExtrema = aProp.FindCurvatureExtrema(); + EXPECT_FALSE(aExtrema.IsDone); + + const Geom2dProp::CurveAnalysis aInfl = aProp.FindInflections(); + EXPECT_FALSE(aInfl.IsDone); +} + +TEST(Geom2dProp_CurveTest, DefaultConstructor_NotInitialized) +{ + Geom2dProp_Curve aProp; + EXPECT_FALSE(aProp.IsInitialized()); +} + +TEST(Geom2dProp_CurveTest, ReInitialize_ChangesType) +{ + Geom2dProp_Curve aProp; + + // First: circle + gp_Circ2d aCirc(gp_Ax2d(gp_Pnt2d(0.0, 0.0), gp_Dir2d(1.0, 0.0)), 5.0); + occ::handle aCircle = new Geom2d_Circle(aCirc); + aProp.Initialize(aCircle); + EXPECT_EQ(aProp.GetType(), GeomAbs_Circle); + + // Re-initialize with line + gp_Lin2d aLin(gp_Pnt2d(0.0, 0.0), gp_Dir2d(1.0, 0.0)); + occ::handle aLine = new Geom2d_Line(aLin); + aProp.Initialize(aLine); + EXPECT_EQ(aProp.GetType(), GeomAbs_Line); + EXPECT_NEAR(aProp.Curvature(0.0, Precision::Confusion()).Value, 0.0, Precision::Confusion()); +} + +// ============================================================================ +// Initialize from adaptor tests +// ============================================================================ + +TEST(Geom2dProp_AdaptorTest, InitFromGeom2dAdaptor) +{ + gp_Circ2d aCirc(gp_Ax2d(gp_Pnt2d(0.0, 0.0), gp_Dir2d(1.0, 0.0)), 5.0); + occ::handle aCircle = new Geom2d_Circle(aCirc); + Geom2dAdaptor_Curve anAdaptor(aCircle); + + Geom2dProp_Curve aProp; + aProp.Initialize(anAdaptor); + EXPECT_TRUE(aProp.IsInitialized()); + EXPECT_EQ(aProp.GetType(), GeomAbs_Circle); + + const Geom2dProp::CurvatureResult aCurv = aProp.Curvature(0.0, Precision::Confusion()); + ASSERT_TRUE(aCurv.IsDefined); + EXPECT_NEAR(aCurv.Value, 1.0 / 5.0, Precision::Confusion()); +} + +TEST(Geom2dProp_AdaptorTest, InitFromGeom2dAdaptor_Ellipse) +{ + gp_Elips2d anElips(gp_Ax2d(gp_Pnt2d(0.0, 0.0), gp_Dir2d(1.0, 0.0)), 8.0, 3.0); + occ::handle anEllipse = new Geom2d_Ellipse(anElips); + Geom2dAdaptor_Curve anAdaptor(anEllipse); + + Geom2dProp_Curve aProp; + aProp.Initialize(anAdaptor); + EXPECT_TRUE(aProp.IsInitialized()); + EXPECT_EQ(aProp.GetType(), GeomAbs_Ellipse); + + // Curvature at major vertex: a/b^2 = 8/9 + const Geom2dProp::CurvatureResult aCurv = aProp.Curvature(0.0, Precision::Confusion()); + ASSERT_TRUE(aCurv.IsDefined); + EXPECT_NEAR(aCurv.Value, 8.0 / 9.0, 1.0e-6); +} + +TEST(Geom2dProp_AdaptorTest, InitFromGeom2dAdaptor_BSpline) +{ + NCollection_Array1 aPoles(1, 4); + aPoles(1) = gp_Pnt2d(0.0, 0.0); + aPoles(2) = gp_Pnt2d(1.0, 2.0); + aPoles(3) = gp_Pnt2d(3.0, 2.0); + aPoles(4) = gp_Pnt2d(4.0, 0.0); + + NCollection_Array1 aKnots(1, 3); + aKnots(1) = 0.0; + aKnots(2) = 0.5; + aKnots(3) = 1.0; + + NCollection_Array1 aMults(1, 3); + aMults(1) = 3; + aMults(2) = 1; + aMults(3) = 3; + + occ::handle aBSpline = new Geom2d_BSplineCurve(aPoles, aKnots, aMults, 2); + Geom2dAdaptor_Curve anAdaptor(aBSpline); + + Geom2dProp_Curve aProp; + aProp.Initialize(anAdaptor); + EXPECT_TRUE(aProp.IsInitialized()); + EXPECT_EQ(aProp.GetType(), GeomAbs_BSplineCurve); +} + +// ============================================================================ +// Cross-validation tests +// ============================================================================ + +TEST(Geom2dProp_CrossValidationTest, Circle_MatchesLProp) +{ + gp_Circ2d aCirc(gp_Ax2d(gp_Pnt2d(1.0, 2.0), gp_Dir2d(1.0, 0.0)), 7.0); + occ::handle aCircle = new Geom2d_Circle(aCirc); + + Geom2dProp_Curve aProp; + aProp.Initialize(aCircle); + + for (double u = 0.0; u < 2.0 * M_PI; u += M_PI / 6.0) + { + const Geom2dProp::CurvatureResult aCurv = aProp.Curvature(u, Precision::Confusion()); + ASSERT_TRUE(aCurv.IsDefined); + EXPECT_NEAR(aCurv.Value, 1.0 / 7.0, 1.0e-10); + + const Geom2dProp::CentreResult aCentre = aProp.CentreOfCurvature(u, Precision::Confusion()); + ASSERT_TRUE(aCentre.IsDefined); + EXPECT_NEAR(aCentre.Centre.X(), 1.0, Precision::Confusion()); + EXPECT_NEAR(aCentre.Centre.Y(), 2.0, Precision::Confusion()); + } +} + +TEST(Geom2dProp_CrossValidationTest, Ellipse_MatchesLProp) +{ + gp_Elips2d anElips(gp_Ax2d(gp_Pnt2d(0.0, 0.0), gp_Dir2d(1.0, 0.0)), 10.0, 5.0); + occ::handle anEllipse = new Geom2d_Ellipse(anElips); + + Geom2dProp_Curve aProp; + aProp.Initialize(anEllipse); + + // At major vertex (U=0): curvature = a/b^2 = 10/25 = 0.4 + const Geom2dProp::CurvatureResult aCurv0 = aProp.Curvature(0.0, Precision::Confusion()); + ASSERT_TRUE(aCurv0.IsDefined); + EXPECT_NEAR(aCurv0.Value, 10.0 / 25.0, 1.0e-6); + + // At minor vertex (U=PI/2): curvature = b/a^2 = 5/100 = 0.05 + const Geom2dProp::CurvatureResult aCurvPi2 = aProp.Curvature(M_PI / 2.0, Precision::Confusion()); + ASSERT_TRUE(aCurvPi2.IsDefined); + EXPECT_NEAR(aCurvPi2.Value, 5.0 / 100.0, 1.0e-6); +} + +// Cross-validate: free function vs. dispatcher give same result +TEST(Geom2dProp_CrossValidationTest, FreeFunctionVsDispatcher) +{ + gp_Circ2d aCirc(gp_Ax2d(gp_Pnt2d(0.0, 0.0), gp_Dir2d(1.0, 0.0)), 3.0); + occ::handle aCircle = new Geom2d_Circle(aCirc); + + Geom2dProp_Curve aProp; + aProp.Initialize(aCircle); + + const double u = 1.0; + gp_Pnt2d aPnt; + gp_Vec2d aD1, aD2; + aCircle->D2(u, aPnt, aD1, aD2); + + // Free function + const Geom2dProp::CurvatureResult aCurvFree = + Geom2dProp::ComputeCurvature(aD1, aD2, Precision::Confusion()); + // Dispatcher + const Geom2dProp::CurvatureResult aCurvDisp = aProp.Curvature(u, Precision::Confusion()); + + ASSERT_TRUE(aCurvFree.IsDefined); + ASSERT_TRUE(aCurvDisp.IsDefined); + EXPECT_NEAR(aCurvFree.Value, aCurvDisp.Value, 1.0e-10); +} + +// Cross-validate centre of curvature consistency: distance from point = 1/|curvature| +TEST(Geom2dProp_CrossValidationTest, CentreDistanceConsistency) +{ + gp_Elips2d anElips(gp_Ax2d(gp_Pnt2d(0.0, 0.0), gp_Dir2d(1.0, 0.0)), 8.0, 4.0); + occ::handle anEllipse = new Geom2d_Ellipse(anElips); + + Geom2dProp_Curve aProp; + aProp.Initialize(anEllipse); + + for (double u = 0.0; u < 2.0 * M_PI; u += M_PI / 8.0) + { + const Geom2dProp::CurvatureResult aCurv = aProp.Curvature(u, Precision::Confusion()); + const Geom2dProp::CentreResult aCentre = aProp.CentreOfCurvature(u, Precision::Confusion()); + if (aCurv.IsDefined && aCentre.IsDefined && !aCurv.IsInfinite + && std::abs(aCurv.Value) > 1.0e-10) + { + gp_Pnt2d aPnt; + gp_Vec2d aD1; + anEllipse->D1(u, aPnt, aD1); + const double aDist = aPnt.Distance(aCentre.Centre); + const double aExpected = 1.0 / std::abs(aCurv.Value); + EXPECT_NEAR(aDist, aExpected, 1.0e-6); + } + } +} + +// Cross-validate: adaptor init and geometry init give same results +TEST(Geom2dProp_CrossValidationTest, AdaptorVsGeometryInit) +{ + gp_Elips2d anElips(gp_Ax2d(gp_Pnt2d(1.0, 2.0), gp_Dir2d(1.0, 0.0)), 6.0, 3.0); + occ::handle anEllipse = new Geom2d_Ellipse(anElips); + + Geom2dProp_Curve aPropGeom; + aPropGeom.Initialize(anEllipse); + + Geom2dAdaptor_Curve anAdaptor(anEllipse); + Geom2dProp_Curve aPropAdap; + aPropAdap.Initialize(anAdaptor); + + for (double u = 0.0; u < 2.0 * M_PI; u += M_PI / 5.0) + { + const Geom2dProp::CurvatureResult aCurvG = aPropGeom.Curvature(u, Precision::Confusion()); + const Geom2dProp::CurvatureResult aCurvA = aPropAdap.Curvature(u, Precision::Confusion()); + ASSERT_TRUE(aCurvG.IsDefined); + ASSERT_TRUE(aCurvA.IsDefined); + EXPECT_NEAR(aCurvG.Value, aCurvA.Value, 1.0e-10); + } +} diff --git a/src/ModelingData/TKG2d/GTests/Geom2dProp_VsCLProps2d_Test.cxx b/src/ModelingData/TKG2d/GTests/Geom2dProp_VsCLProps2d_Test.cxx new file mode 100644 index 0000000000..9133cdd988 --- /dev/null +++ b/src/ModelingData/TKG2d/GTests/Geom2dProp_VsCLProps2d_Test.cxx @@ -0,0 +1,603 @@ +// Copyright (c) 2025 OPEN CASCADE SAS +// +// This file is part of Open CASCADE Technology software library. +// +// This library is free software; you can redistribute it and/or modify it under +// the terms of the GNU Lesser General Public License version 2.1 as published +// by the Free Software Foundation, with special exception defined in the file +// OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT +// distribution for complete text of the license and disclaimer of any warranty. +// +// Alternatively, this file may be used under the terms of Open CASCADE +// commercial license or contractual agreement. + +// Cross-validation tests comparing Geom2dProp_Curve against Geom2dLProp_CLProps2d +// for local differential properties (tangent, curvature, normal, centre of curvature). + +#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 +{ +constexpr double THE_PARAM_TOL = Precision::Confusion(); +constexpr double THE_VALUE_TOL = 1.0e-10; +constexpr double THE_DIR_TOL = 1.0e-10; +constexpr double THE_POINT_TOL = 1.0e-8; + +//! Compare tangent from new Geom2dProp vs old CLProps2d at given parameter. +void compareTangent(Geom2dProp_Curve& theProp, + Geom2dLProp_CLProps2d& theOld, + const double theParam, + const ::testing::TestInfo* = nullptr) +{ + theOld.SetParameter(theParam); + + const Geom2dProp::TangentResult aNewTan = theProp.Tangent(theParam, THE_PARAM_TOL); + const bool aOldDefined = theOld.IsTangentDefined(); + + EXPECT_EQ(aNewTan.IsDefined, aOldDefined) << "Tangent defined mismatch at U=" << theParam; + + if (aNewTan.IsDefined && aOldDefined) + { + gp_Dir2d aOldDir; + theOld.Tangent(aOldDir); + // Tangent directions may differ by sign; compare absolute dot product + const double aDot = aNewTan.Direction.X() * aOldDir.X() + aNewTan.Direction.Y() * aOldDir.Y(); + EXPECT_NEAR(std::abs(aDot), 1.0, THE_DIR_TOL) << "Tangent direction mismatch at U=" << theParam; + } +} + +//! Compare curvature from new Geom2dProp vs old CLProps2d at given parameter. +void compareCurvature(Geom2dProp_Curve& theProp, + Geom2dLProp_CLProps2d& theOld, + const double theParam) +{ + theOld.SetParameter(theParam); + + const Geom2dProp::CurvatureResult aNewCurv = theProp.Curvature(theParam, THE_PARAM_TOL); + const double aOldCurv = theOld.Curvature(); + + if (aNewCurv.IsDefined && !aNewCurv.IsInfinite) + { + EXPECT_NEAR(aNewCurv.Value, aOldCurv, THE_VALUE_TOL) << "Curvature mismatch at U=" << theParam; + } +} + +//! Compare normal from new Geom2dProp vs old CLProps2d at given parameter. +void compareNormal(Geom2dProp_Curve& theProp, Geom2dLProp_CLProps2d& theOld, const double theParam) +{ + theOld.SetParameter(theParam); + + const Geom2dProp::NormalResult aNewNorm = theProp.Normal(theParam, THE_PARAM_TOL); + + // Old API: Normal() throws if curvature is nearly zero; curvature check needed + const double aOldCurv = theOld.Curvature(); + if (std::abs(aOldCurv) < THE_PARAM_TOL) + { + // Old API would throw - new API returns IsDefined=false + EXPECT_FALSE(aNewNorm.IsDefined) << "Normal should be undefined at U=" << theParam; + return; + } + + if (aNewNorm.IsDefined) + { + gp_Dir2d aOldDir; + theOld.Normal(aOldDir); + const double aDot = aNewNorm.Direction.X() * aOldDir.X() + aNewNorm.Direction.Y() * aOldDir.Y(); + EXPECT_NEAR(std::abs(aDot), 1.0, THE_DIR_TOL) << "Normal direction mismatch at U=" << theParam; + } +} + +//! Compare centre of curvature from new Geom2dProp vs old CLProps2d at given parameter. +void compareCentre(Geom2dProp_Curve& theProp, Geom2dLProp_CLProps2d& theOld, const double theParam) +{ + theOld.SetParameter(theParam); + + const Geom2dProp::CentreResult aNewCentre = theProp.CentreOfCurvature(theParam, THE_PARAM_TOL); + + // Old API: CentreOfCurvature() throws if curvature is nearly zero + const double aOldCurv = theOld.Curvature(); + if (std::abs(aOldCurv) < THE_PARAM_TOL) + { + EXPECT_FALSE(aNewCentre.IsDefined) << "Centre should be undefined at U=" << theParam; + return; + } + + if (aNewCentre.IsDefined) + { + gp_Pnt2d aOldCentre; + theOld.CentreOfCurvature(aOldCentre); + EXPECT_NEAR(aNewCentre.Centre.X(), aOldCentre.X(), THE_POINT_TOL) + << "Centre X mismatch at U=" << theParam; + EXPECT_NEAR(aNewCentre.Centre.Y(), aOldCentre.Y(), THE_POINT_TOL) + << "Centre Y mismatch at U=" << theParam; + } +} + +//! Run all four property comparisons at given parameter. +void compareAllProperties(Geom2dProp_Curve& theProp, + Geom2dLProp_CLProps2d& theOld, + const double theParam) +{ + compareTangent(theProp, theOld, theParam); + compareCurvature(theProp, theOld, theParam); + compareNormal(theProp, theOld, theParam); + compareCentre(theProp, theOld, theParam); +} + +} // namespace + +// ============================================================================ +// Line +// ============================================================================ + +TEST(Geom2dProp_VsCLProps2dTest, Line_Tangent) +{ + gp_Lin2d aLin(gp_Pnt2d(1.0, 2.0), gp_Dir2d(3.0, 4.0)); + occ::handle aLine = new Geom2d_Line(aLin); + + Geom2dProp_Curve aProp; + aProp.Initialize(aLine); + Geom2dLProp_CLProps2d aOld(aLine, 2, THE_PARAM_TOL); + + for (double u = -10.0; u <= 10.0; u += 2.5) + { + compareTangent(aProp, aOld, u); + } +} + +TEST(Geom2dProp_VsCLProps2dTest, Line_Curvature) +{ + gp_Lin2d aLin(gp_Pnt2d(0.0, 0.0), gp_Dir2d(1.0, 1.0)); + occ::handle aLine = new Geom2d_Line(aLin); + + Geom2dProp_Curve aProp; + aProp.Initialize(aLine); + Geom2dLProp_CLProps2d aOld(aLine, 2, THE_PARAM_TOL); + + for (double u = -5.0; u <= 5.0; u += 1.0) + { + compareCurvature(aProp, aOld, u); + } +} + +// ============================================================================ +// Circle +// ============================================================================ + +TEST(Geom2dProp_VsCLProps2dTest, Circle_AllProperties) +{ + gp_Circ2d aCirc(gp_Ax2d(gp_Pnt2d(3.0, 4.0), gp_Dir2d(1.0, 0.0)), 7.0); + occ::handle aCircle = new Geom2d_Circle(aCirc); + + Geom2dProp_Curve aProp; + aProp.Initialize(aCircle); + Geom2dLProp_CLProps2d aOld(aCircle, 2, THE_PARAM_TOL); + + for (double u = 0.0; u < 2.0 * M_PI; u += M_PI / 6.0) + { + compareAllProperties(aProp, aOld, u); + } +} + +TEST(Geom2dProp_VsCLProps2dTest, Circle_SmallRadius) +{ + gp_Circ2d aCirc(gp_Ax2d(gp_Pnt2d(0.0, 0.0), gp_Dir2d(1.0, 0.0)), 0.01); + occ::handle aCircle = new Geom2d_Circle(aCirc); + + Geom2dProp_Curve aProp; + aProp.Initialize(aCircle); + Geom2dLProp_CLProps2d aOld(aCircle, 2, THE_PARAM_TOL); + + for (double u = 0.0; u < 2.0 * M_PI; u += M_PI / 4.0) + { + compareAllProperties(aProp, aOld, u); + } +} + +TEST(Geom2dProp_VsCLProps2dTest, Circle_LargeRadius) +{ + gp_Circ2d aCirc(gp_Ax2d(gp_Pnt2d(0.0, 0.0), gp_Dir2d(1.0, 0.0)), 1000.0); + occ::handle aCircle = new Geom2d_Circle(aCirc); + + Geom2dProp_Curve aProp; + aProp.Initialize(aCircle); + Geom2dLProp_CLProps2d aOld(aCircle, 2, THE_PARAM_TOL); + + for (double u = 0.0; u < 2.0 * M_PI; u += M_PI / 4.0) + { + compareAllProperties(aProp, aOld, u); + } +} + +// ============================================================================ +// Ellipse +// ============================================================================ + +TEST(Geom2dProp_VsCLProps2dTest, Ellipse_AllProperties) +{ + gp_Elips2d anElips(gp_Ax2d(gp_Pnt2d(0.0, 0.0), gp_Dir2d(1.0, 0.0)), 10.0, 5.0); + occ::handle anEllipse = new Geom2d_Ellipse(anElips); + + Geom2dProp_Curve aProp; + aProp.Initialize(anEllipse); + Geom2dLProp_CLProps2d aOld(anEllipse, 2, THE_PARAM_TOL); + + for (double u = 0.0; u < 2.0 * M_PI; u += M_PI / 12.0) + { + compareAllProperties(aProp, aOld, u); + } +} + +TEST(Geom2dProp_VsCLProps2dTest, Ellipse_HighEccentricity) +{ + gp_Elips2d anElips(gp_Ax2d(gp_Pnt2d(0.0, 0.0), gp_Dir2d(1.0, 0.0)), 100.0, 1.0); + occ::handle anEllipse = new Geom2d_Ellipse(anElips); + + Geom2dProp_Curve aProp; + aProp.Initialize(anEllipse); + Geom2dLProp_CLProps2d aOld(anEllipse, 2, THE_PARAM_TOL); + + for (double u = 0.0; u < 2.0 * M_PI; u += M_PI / 8.0) + { + compareAllProperties(aProp, aOld, u); + } +} + +TEST(Geom2dProp_VsCLProps2dTest, Ellipse_OffCenter) +{ + gp_Elips2d anElips(gp_Ax2d(gp_Pnt2d(100.0, -50.0), gp_Dir2d(1.0, 0.0)), 8.0, 3.0); + occ::handle anEllipse = new Geom2d_Ellipse(anElips); + + Geom2dProp_Curve aProp; + aProp.Initialize(anEllipse); + Geom2dLProp_CLProps2d aOld(anEllipse, 2, THE_PARAM_TOL); + + for (double u = 0.0; u < 2.0 * M_PI; u += M_PI / 8.0) + { + compareAllProperties(aProp, aOld, u); + } +} + +// ============================================================================ +// Hyperbola +// ============================================================================ + +TEST(Geom2dProp_VsCLProps2dTest, Hyperbola_AllProperties) +{ + gp_Hypr2d anHypr(gp_Ax2d(gp_Pnt2d(0.0, 0.0), gp_Dir2d(1.0, 0.0)), 6.0, 3.0); + occ::handle aHyperbola = new Geom2d_Hyperbola(anHypr); + + Geom2dProp_Curve aProp; + aProp.Initialize(aHyperbola); + Geom2dLProp_CLProps2d aOld(aHyperbola, 2, THE_PARAM_TOL); + + for (double u = -2.0; u <= 2.0; u += 0.5) + { + compareAllProperties(aProp, aOld, u); + } +} + +TEST(Geom2dProp_VsCLProps2dTest, Hyperbola_NearVertex) +{ + gp_Hypr2d anHypr(gp_Ax2d(gp_Pnt2d(0.0, 0.0), gp_Dir2d(1.0, 0.0)), 4.0, 2.0); + occ::handle aHyperbola = new Geom2d_Hyperbola(anHypr); + + Geom2dProp_Curve aProp; + aProp.Initialize(aHyperbola); + Geom2dLProp_CLProps2d aOld(aHyperbola, 2, THE_PARAM_TOL); + + // Fine-grained near vertex + for (double u = -0.5; u <= 0.5; u += 0.1) + { + compareAllProperties(aProp, aOld, u); + } +} + +// ============================================================================ +// Parabola +// ============================================================================ + +TEST(Geom2dProp_VsCLProps2dTest, Parabola_AllProperties) +{ + gp_Parab2d aParab(gp_Ax2d(gp_Pnt2d(0.0, 0.0), gp_Dir2d(1.0, 0.0)), 2.0); + occ::handle aParabola = new Geom2d_Parabola(aParab); + + Geom2dProp_Curve aProp; + aProp.Initialize(aParabola); + Geom2dLProp_CLProps2d aOld(aParabola, 2, THE_PARAM_TOL); + + for (double u = -5.0; u <= 5.0; u += 1.0) + { + compareAllProperties(aProp, aOld, u); + } +} + +TEST(Geom2dProp_VsCLProps2dTest, Parabola_SmallFocal) +{ + gp_Parab2d aParab(gp_Ax2d(gp_Pnt2d(0.0, 0.0), gp_Dir2d(1.0, 0.0)), 0.1); + occ::handle aParabola = new Geom2d_Parabola(aParab); + + Geom2dProp_Curve aProp; + aProp.Initialize(aParabola); + Geom2dLProp_CLProps2d aOld(aParabola, 2, THE_PARAM_TOL); + + for (double u = -3.0; u <= 3.0; u += 0.5) + { + compareAllProperties(aProp, aOld, u); + } +} + +TEST(Geom2dProp_VsCLProps2dTest, Parabola_LargeFocal) +{ + gp_Parab2d aParab(gp_Ax2d(gp_Pnt2d(0.0, 0.0), gp_Dir2d(1.0, 0.0)), 50.0); + occ::handle aParabola = new Geom2d_Parabola(aParab); + + Geom2dProp_Curve aProp; + aProp.Initialize(aParabola); + Geom2dLProp_CLProps2d aOld(aParabola, 2, THE_PARAM_TOL); + + for (double u = -10.0; u <= 10.0; u += 2.0) + { + compareAllProperties(aProp, aOld, u); + } +} + +// ============================================================================ +// Bezier curve +// ============================================================================ + +TEST(Geom2dProp_VsCLProps2dTest, Bezier_CubicSShape) +{ + NCollection_Array1 aPoles(1, 4); + aPoles(1) = gp_Pnt2d(0.0, 0.0); + aPoles(2) = gp_Pnt2d(1.0, 2.0); + aPoles(3) = gp_Pnt2d(3.0, -2.0); + aPoles(4) = gp_Pnt2d(4.0, 0.0); + occ::handle aBezier = new Geom2d_BezierCurve(aPoles); + + Geom2dProp_Curve aProp; + aProp.Initialize(aBezier); + Geom2dLProp_CLProps2d aOld(aBezier, 2, THE_PARAM_TOL); + + for (double u = 0.0; u <= 1.0; u += 0.1) + { + compareAllProperties(aProp, aOld, u); + } +} + +TEST(Geom2dProp_VsCLProps2dTest, Bezier_Quadratic) +{ + NCollection_Array1 aPoles(1, 3); + aPoles(1) = gp_Pnt2d(0.0, 0.0); + aPoles(2) = gp_Pnt2d(2.0, 4.0); + aPoles(3) = gp_Pnt2d(4.0, 0.0); + occ::handle aBezier = new Geom2d_BezierCurve(aPoles); + + Geom2dProp_Curve aProp; + aProp.Initialize(aBezier); + Geom2dLProp_CLProps2d aOld(aBezier, 2, THE_PARAM_TOL); + + for (double u = 0.0; u <= 1.0; u += 0.1) + { + compareAllProperties(aProp, aOld, u); + } +} + +TEST(Geom2dProp_VsCLProps2dTest, Bezier_HighDegree) +{ + NCollection_Array1 aPoles(1, 6); + aPoles(1) = gp_Pnt2d(0.0, 0.0); + aPoles(2) = gp_Pnt2d(1.0, 3.0); + aPoles(3) = gp_Pnt2d(2.0, -1.0); + aPoles(4) = gp_Pnt2d(3.0, 2.0); + aPoles(5) = gp_Pnt2d(4.0, -2.0); + aPoles(6) = gp_Pnt2d(5.0, 1.0); + occ::handle aBezier = new Geom2d_BezierCurve(aPoles); + + Geom2dProp_Curve aProp; + aProp.Initialize(aBezier); + Geom2dLProp_CLProps2d aOld(aBezier, 2, THE_PARAM_TOL); + + for (double u = 0.0; u <= 1.0; u += 0.05) + { + compareAllProperties(aProp, aOld, u); + } +} + +// ============================================================================ +// BSpline curve +// ============================================================================ + +TEST(Geom2dProp_VsCLProps2dTest, BSpline_Quadratic) +{ + NCollection_Array1 aPoles(1, 4); + aPoles(1) = gp_Pnt2d(0.0, 0.0); + aPoles(2) = gp_Pnt2d(1.0, 2.0); + aPoles(3) = gp_Pnt2d(3.0, 2.0); + aPoles(4) = gp_Pnt2d(4.0, 0.0); + + NCollection_Array1 aKnots(1, 3); + aKnots(1) = 0.0; + aKnots(2) = 0.5; + aKnots(3) = 1.0; + + NCollection_Array1 aMults(1, 3); + aMults(1) = 3; + aMults(2) = 1; + aMults(3) = 3; + + occ::handle aBSpline = new Geom2d_BSplineCurve(aPoles, aKnots, aMults, 2); + + Geom2dProp_Curve aProp; + aProp.Initialize(aBSpline); + Geom2dLProp_CLProps2d aOld(aBSpline, 2, THE_PARAM_TOL); + + for (double u = 0.0; u <= 1.0; u += 0.1) + { + compareAllProperties(aProp, aOld, u); + } +} + +TEST(Geom2dProp_VsCLProps2dTest, BSpline_Cubic) +{ + NCollection_Array1 aPoles(1, 6); + aPoles(1) = gp_Pnt2d(0.0, 0.0); + aPoles(2) = gp_Pnt2d(1.0, 3.0); + aPoles(3) = gp_Pnt2d(2.0, 1.0); + aPoles(4) = gp_Pnt2d(3.0, 4.0); + aPoles(5) = gp_Pnt2d(4.0, 2.0); + aPoles(6) = gp_Pnt2d(5.0, 0.0); + + NCollection_Array1 aKnots(1, 4); + aKnots(1) = 0.0; + aKnots(2) = 0.33; + aKnots(3) = 0.66; + aKnots(4) = 1.0; + + NCollection_Array1 aMults(1, 4); + aMults(1) = 4; + aMults(2) = 1; + aMults(3) = 1; + aMults(4) = 4; + + occ::handle aBSpline = new Geom2d_BSplineCurve(aPoles, aKnots, aMults, 3); + + Geom2dProp_Curve aProp; + aProp.Initialize(aBSpline); + Geom2dLProp_CLProps2d aOld(aBSpline, 2, THE_PARAM_TOL); + + for (double u = 0.0; u <= 1.0; u += 0.05) + { + compareAllProperties(aProp, aOld, u); + } +} + +TEST(Geom2dProp_VsCLProps2dTest, BSpline_Degree4) +{ + NCollection_Array1 aPoles(1, 5); + aPoles(1) = gp_Pnt2d(0.0, 0.0); + aPoles(2) = gp_Pnt2d(1.0, 3.0); + aPoles(3) = gp_Pnt2d(2.0, -1.0); + aPoles(4) = gp_Pnt2d(3.0, 2.0); + aPoles(5) = gp_Pnt2d(4.0, 0.0); + + NCollection_Array1 aKnots(1, 2); + aKnots(1) = 0.0; + aKnots(2) = 1.0; + + NCollection_Array1 aMults(1, 2); + aMults(1) = 5; + aMults(2) = 5; + + occ::handle aBSpline = new Geom2d_BSplineCurve(aPoles, aKnots, aMults, 4); + + Geom2dProp_Curve aProp; + aProp.Initialize(aBSpline); + Geom2dLProp_CLProps2d aOld(aBSpline, 2, THE_PARAM_TOL); + + for (double u = 0.0; u <= 1.0; u += 0.05) + { + compareAllProperties(aProp, aOld, u); + } +} + +// ============================================================================ +// Offset curve +// ============================================================================ + +TEST(Geom2dProp_VsCLProps2dTest, OffsetCircle_AllProperties) +{ + gp_Circ2d aCirc(gp_Ax2d(gp_Pnt2d(0.0, 0.0), gp_Dir2d(1.0, 0.0)), 5.0); + occ::handle aCircle = new Geom2d_Circle(aCirc); + occ::handle anOffset = new Geom2d_OffsetCurve(aCircle, 2.0); + + Geom2dProp_Curve aProp; + aProp.Initialize(anOffset); + Geom2dLProp_CLProps2d aOld(anOffset, 2, THE_PARAM_TOL); + + for (double u = 0.0; u < 2.0 * M_PI; u += M_PI / 6.0) + { + compareAllProperties(aProp, aOld, u); + } +} + +TEST(Geom2dProp_VsCLProps2dTest, OffsetEllipse_AllProperties) +{ + gp_Elips2d anElips(gp_Ax2d(gp_Pnt2d(0.0, 0.0), gp_Dir2d(1.0, 0.0)), 10.0, 5.0); + occ::handle anEllipse = new Geom2d_Ellipse(anElips); + occ::handle anOffset = new Geom2d_OffsetCurve(anEllipse, 1.0); + + Geom2dProp_Curve aProp; + aProp.Initialize(anOffset); + Geom2dLProp_CLProps2d aOld(anOffset, 2, THE_PARAM_TOL); + + for (double u = 0.0; u < 2.0 * M_PI; u += M_PI / 8.0) + { + compareAllProperties(aProp, aOld, u); + } +} + +// ============================================================================ +// Trimmed curve +// ============================================================================ + +TEST(Geom2dProp_VsCLProps2dTest, TrimmedEllipse_AllProperties) +{ + gp_Elips2d anElips(gp_Ax2d(gp_Pnt2d(0.0, 0.0), gp_Dir2d(1.0, 0.0)), 8.0, 4.0); + occ::handle anEllipse = new Geom2d_Ellipse(anElips); + occ::handle aTrimmed = new Geom2d_TrimmedCurve(anEllipse, 0.5, 2.5); + + Geom2dProp_Curve aProp; + aProp.Initialize(aTrimmed); + Geom2dLProp_CLProps2d aOld(aTrimmed, 2, THE_PARAM_TOL); + + for (double u = 0.5; u <= 2.5; u += 0.2) + { + compareAllProperties(aProp, aOld, u); + } +} + +TEST(Geom2dProp_VsCLProps2dTest, TrimmedBezier_AllProperties) +{ + NCollection_Array1 aPoles(1, 4); + aPoles(1) = gp_Pnt2d(0.0, 0.0); + aPoles(2) = gp_Pnt2d(1.0, 3.0); + aPoles(3) = gp_Pnt2d(3.0, -1.0); + aPoles(4) = gp_Pnt2d(4.0, 1.0); + occ::handle aBezier = new Geom2d_BezierCurve(aPoles); + occ::handle aTrimmed = new Geom2d_TrimmedCurve(aBezier, 0.2, 0.8); + + Geom2dProp_Curve aProp; + aProp.Initialize(aTrimmed); + Geom2dLProp_CLProps2d aOld(aTrimmed, 2, THE_PARAM_TOL); + + for (double u = 0.2; u <= 0.8; u += 0.1) + { + compareAllProperties(aProp, aOld, u); + } +} diff --git a/src/ModelingData/TKG2d/GTests/Geom2dProp_VsLProp_Test.cxx b/src/ModelingData/TKG2d/GTests/Geom2dProp_VsLProp_Test.cxx new file mode 100644 index 0000000000..653b0b0631 --- /dev/null +++ b/src/ModelingData/TKG2d/GTests/Geom2dProp_VsLProp_Test.cxx @@ -0,0 +1,685 @@ +// Copyright (c) 2025 OPEN CASCADE SAS +// +// This file is part of Open CASCADE Technology software library. +// +// This library is free software; you can redistribute it and/or modify it under +// the terms of the GNU Lesser General Public License version 2.1 as published +// by the Free Software Foundation, with special exception defined in the file +// OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT +// distribution for complete text of the license and disclaimer of any warranty. +// +// Alternatively, this file may be used under the terms of Open CASCADE +// commercial license or contractual agreement. + +// Cross-validation tests comparing Geom2dProp_Curve against Geom2dLProp_CurAndInf2d +// for global curve analysis (curvature extrema and inflection point finding). + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +namespace +{ +constexpr double THE_PARAM_TOL = 1.0e-4; + +//! Map LProp_CIType to Geom2dProp::CIType for comparison. +Geom2dProp::CIType mapLPropType(const LProp_CIType theType) +{ + switch (theType) + { + case LProp_Inflection: + return Geom2dProp::CIType::Inflection; + case LProp_MinCur: + return Geom2dProp::CIType::MinCurvature; + case LProp_MaxCur: + return Geom2dProp::CIType::MaxCurvature; + } + return Geom2dProp::CIType::Inflection; +} + +//! Compare extrema results from old and new APIs. +void compareExtrema(const Geom2dProp::CurveAnalysis& theNew, + const Geom2dLProp_CurAndInf2d& theOld, + const double theTol = THE_PARAM_TOL) +{ + EXPECT_EQ(theNew.Points.Length(), theOld.NbPoints()); + + const int aNb = std::min(theNew.Points.Length(), theOld.NbPoints()); + for (int i = 0; i < aNb; ++i) + { + EXPECT_NEAR(theNew.Points.Value(i).Parameter, theOld.Parameter(i + 1), theTol) + << "Parameter mismatch at index " << i; + EXPECT_EQ(theNew.Points.Value(i).Type, mapLPropType(theOld.Type(i + 1))) + << "Type mismatch at index " << i; + } +} + +} // namespace + +// ============================================================================ +// Circle - no extrema, no inflections +// ============================================================================ + +TEST(Geom2dProp_VsCurAndInf2dTest, Circle_NoExtrema) +{ + gp_Circ2d aCirc(gp_Ax2d(gp_Pnt2d(0.0, 0.0), gp_Dir2d(1.0, 0.0)), 5.0); + occ::handle aCircle = new Geom2d_Circle(aCirc); + + Geom2dLProp_CurAndInf2d anOld; + anOld.PerformCurExt(aCircle); + ASSERT_TRUE(anOld.IsDone()); + EXPECT_EQ(anOld.NbPoints(), 0); + + Geom2dProp_Curve aProp; + aProp.Initialize(aCircle); + const Geom2dProp::CurveAnalysis aNew = aProp.FindCurvatureExtrema(); + ASSERT_TRUE(aNew.IsDone); + EXPECT_EQ(aNew.Points.Length(), 0); +} + +TEST(Geom2dProp_VsCurAndInf2dTest, Circle_NoInflections) +{ + gp_Circ2d aCirc(gp_Ax2d(gp_Pnt2d(0.0, 0.0), gp_Dir2d(1.0, 0.0)), 5.0); + occ::handle aCircle = new Geom2d_Circle(aCirc); + + Geom2dLProp_CurAndInf2d anOld; + anOld.PerformInf(aCircle); + ASSERT_TRUE(anOld.IsDone()); + EXPECT_EQ(anOld.NbPoints(), 0); + + Geom2dProp_Curve aProp; + aProp.Initialize(aCircle); + const Geom2dProp::CurveAnalysis aNew = aProp.FindInflections(); + ASSERT_TRUE(aNew.IsDone); + EXPECT_EQ(aNew.Points.Length(), 0); +} + +// ============================================================================ +// Ellipse - 4 extrema, no inflections +// ============================================================================ + +TEST(Geom2dProp_VsCurAndInf2dTest, Ellipse_Extrema) +{ + gp_Elips2d anElips(gp_Ax2d(gp_Pnt2d(0.0, 0.0), gp_Dir2d(1.0, 0.0)), 10.0, 5.0); + occ::handle anEllipse = new Geom2d_Ellipse(anElips); + + Geom2dLProp_CurAndInf2d anOld; + anOld.PerformCurExt(anEllipse); + ASSERT_TRUE(anOld.IsDone()); + + Geom2dProp_Curve aProp; + aProp.Initialize(anEllipse); + const Geom2dProp::CurveAnalysis aNew = aProp.FindCurvatureExtrema(); + ASSERT_TRUE(aNew.IsDone); + + compareExtrema(aNew, anOld, 1.0e-6); +} + +TEST(Geom2dProp_VsCurAndInf2dTest, Ellipse_NoInflections) +{ + gp_Elips2d anElips(gp_Ax2d(gp_Pnt2d(0.0, 0.0), gp_Dir2d(1.0, 0.0)), 10.0, 5.0); + occ::handle anEllipse = new Geom2d_Ellipse(anElips); + + Geom2dLProp_CurAndInf2d anOld; + anOld.PerformInf(anEllipse); + ASSERT_TRUE(anOld.IsDone()); + EXPECT_EQ(anOld.NbPoints(), 0); + + Geom2dProp_Curve aProp; + aProp.Initialize(anEllipse); + const Geom2dProp::CurveAnalysis aNew = aProp.FindInflections(); + ASSERT_TRUE(aNew.IsDone); + EXPECT_EQ(aNew.Points.Length(), 0); +} + +TEST(Geom2dProp_VsCurAndInf2dTest, Ellipse_HighEccentricity_Extrema) +{ + gp_Elips2d anElips(gp_Ax2d(gp_Pnt2d(0.0, 0.0), gp_Dir2d(1.0, 0.0)), 50.0, 2.0); + occ::handle anEllipse = new Geom2d_Ellipse(anElips); + + Geom2dLProp_CurAndInf2d anOld; + anOld.PerformCurExt(anEllipse); + ASSERT_TRUE(anOld.IsDone()); + + Geom2dProp_Curve aProp; + aProp.Initialize(anEllipse); + const Geom2dProp::CurveAnalysis aNew = aProp.FindCurvatureExtrema(); + ASSERT_TRUE(aNew.IsDone); + + compareExtrema(aNew, anOld, 1.0e-6); +} + +TEST(Geom2dProp_VsCurAndInf2dTest, Ellipse_FullPerform) +{ + gp_Elips2d anElips(gp_Ax2d(gp_Pnt2d(0.0, 0.0), gp_Dir2d(1.0, 0.0)), 8.0, 3.0); + occ::handle anEllipse = new Geom2d_Ellipse(anElips); + + Geom2dLProp_CurAndInf2d anOld; + anOld.Perform(anEllipse); + ASSERT_TRUE(anOld.IsDone()); + + Geom2dProp_Curve aProp; + aProp.Initialize(anEllipse); + const Geom2dProp::CurveAnalysis aNewExt = aProp.FindCurvatureExtrema(); + const Geom2dProp::CurveAnalysis aNewInfl = aProp.FindInflections(); + ASSERT_TRUE(aNewExt.IsDone); + ASSERT_TRUE(aNewInfl.IsDone); + + const int aNewTotal = aNewExt.Points.Length() + aNewInfl.Points.Length(); + EXPECT_EQ(aNewTotal, anOld.NbPoints()); +} + +// ============================================================================ +// Hyperbola - 1 extremum at vertex, no inflections +// ============================================================================ + +TEST(Geom2dProp_VsCurAndInf2dTest, Hyperbola_Extrema) +{ + gp_Hypr2d anHypr(gp_Ax2d(gp_Pnt2d(0.0, 0.0), gp_Dir2d(1.0, 0.0)), 6.0, 3.0); + occ::handle aHyperbola = new Geom2d_Hyperbola(anHypr); + + Geom2dLProp_CurAndInf2d anOld; + anOld.PerformCurExt(aHyperbola); + ASSERT_TRUE(anOld.IsDone()); + + Geom2dProp_Curve aProp; + aProp.Initialize(aHyperbola); + const Geom2dProp::CurveAnalysis aNew = aProp.FindCurvatureExtrema(); + ASSERT_TRUE(aNew.IsDone); + + compareExtrema(aNew, anOld, 1.0e-6); +} + +TEST(Geom2dProp_VsCurAndInf2dTest, Hyperbola_NoInflections) +{ + gp_Hypr2d anHypr(gp_Ax2d(gp_Pnt2d(0.0, 0.0), gp_Dir2d(1.0, 0.0)), 6.0, 3.0); + occ::handle aHyperbola = new Geom2d_Hyperbola(anHypr); + + Geom2dLProp_CurAndInf2d anOld; + anOld.PerformInf(aHyperbola); + ASSERT_TRUE(anOld.IsDone()); + EXPECT_EQ(anOld.NbPoints(), 0); + + Geom2dProp_Curve aProp; + aProp.Initialize(aHyperbola); + const Geom2dProp::CurveAnalysis aNew = aProp.FindInflections(); + ASSERT_TRUE(aNew.IsDone); + EXPECT_EQ(aNew.Points.Length(), 0); +} + +// ============================================================================ +// Parabola - 1 extremum at vertex, no inflections +// ============================================================================ + +TEST(Geom2dProp_VsCurAndInf2dTest, Parabola_Extrema) +{ + gp_Parab2d aParab(gp_Ax2d(gp_Pnt2d(0.0, 0.0), gp_Dir2d(1.0, 0.0)), 2.0); + occ::handle aParabola = new Geom2d_Parabola(aParab); + + Geom2dLProp_CurAndInf2d anOld; + anOld.PerformCurExt(aParabola); + ASSERT_TRUE(anOld.IsDone()); + + Geom2dProp_Curve aProp; + aProp.Initialize(aParabola); + const Geom2dProp::CurveAnalysis aNew = aProp.FindCurvatureExtrema(); + ASSERT_TRUE(aNew.IsDone); + + compareExtrema(aNew, anOld, 1.0e-6); +} + +TEST(Geom2dProp_VsCurAndInf2dTest, Parabola_NoInflections) +{ + gp_Parab2d aParab(gp_Ax2d(gp_Pnt2d(0.0, 0.0), gp_Dir2d(1.0, 0.0)), 2.0); + occ::handle aParabola = new Geom2d_Parabola(aParab); + + Geom2dLProp_CurAndInf2d anOld; + anOld.PerformInf(aParabola); + ASSERT_TRUE(anOld.IsDone()); + EXPECT_EQ(anOld.NbPoints(), 0); + + Geom2dProp_Curve aProp; + aProp.Initialize(aParabola); + const Geom2dProp::CurveAnalysis aNew = aProp.FindInflections(); + ASSERT_TRUE(aNew.IsDone); + EXPECT_EQ(aNew.Points.Length(), 0); +} + +TEST(Geom2dProp_VsCurAndInf2dTest, Parabola_FullPerform) +{ + gp_Parab2d aParab(gp_Ax2d(gp_Pnt2d(0.0, 0.0), gp_Dir2d(1.0, 0.0)), 5.0); + occ::handle aParabola = new Geom2d_Parabola(aParab); + + Geom2dLProp_CurAndInf2d anOld; + anOld.Perform(aParabola); + ASSERT_TRUE(anOld.IsDone()); + + Geom2dProp_Curve aProp; + aProp.Initialize(aParabola); + const Geom2dProp::CurveAnalysis aNewExt = aProp.FindCurvatureExtrema(); + const Geom2dProp::CurveAnalysis aNewInfl = aProp.FindInflections(); + ASSERT_TRUE(aNewExt.IsDone); + ASSERT_TRUE(aNewInfl.IsDone); + + const int aNewTotal = aNewExt.Points.Length() + aNewInfl.Points.Length(); + EXPECT_EQ(aNewTotal, anOld.NbPoints()); +} + +// ============================================================================ +// Bezier - numeric extrema and inflections +// ============================================================================ + +TEST(Geom2dProp_VsCurAndInf2dTest, Bezier_CubicS_Inflections) +{ + NCollection_Array1 aPoles(1, 4); + aPoles(1) = gp_Pnt2d(0.0, 0.0); + aPoles(2) = gp_Pnt2d(1.0, 2.0); + aPoles(3) = gp_Pnt2d(3.0, -1.0); + aPoles(4) = gp_Pnt2d(4.0, 1.0); + occ::handle aBezier = new Geom2d_BezierCurve(aPoles); + + Geom2dLProp_CurAndInf2d anOld; + anOld.PerformInf(aBezier); + ASSERT_TRUE(anOld.IsDone()); + + Geom2dProp_Curve aProp; + aProp.Initialize(aBezier); + const Geom2dProp::CurveAnalysis aNew = aProp.FindInflections(); + ASSERT_TRUE(aNew.IsDone); + + compareExtrema(aNew, anOld); +} + +TEST(Geom2dProp_VsCurAndInf2dTest, Bezier_CubicS_Extrema) +{ + NCollection_Array1 aPoles(1, 4); + aPoles(1) = gp_Pnt2d(0.0, 0.0); + aPoles(2) = gp_Pnt2d(1.0, 2.0); + aPoles(3) = gp_Pnt2d(3.0, -1.0); + aPoles(4) = gp_Pnt2d(4.0, 1.0); + occ::handle aBezier = new Geom2d_BezierCurve(aPoles); + + Geom2dLProp_CurAndInf2d anOld; + anOld.PerformCurExt(aBezier); + ASSERT_TRUE(anOld.IsDone()); + + Geom2dProp_Curve aProp; + aProp.Initialize(aBezier); + const Geom2dProp::CurveAnalysis aNew = aProp.FindCurvatureExtrema(); + ASSERT_TRUE(aNew.IsDone); + + compareExtrema(aNew, anOld); +} + +TEST(Geom2dProp_VsCurAndInf2dTest, Bezier_CubicS_FullPerform) +{ + NCollection_Array1 aPoles(1, 4); + aPoles(1) = gp_Pnt2d(0.0, 0.0); + aPoles(2) = gp_Pnt2d(1.0, 2.0); + aPoles(3) = gp_Pnt2d(3.0, -2.0); + aPoles(4) = gp_Pnt2d(4.0, 0.0); + occ::handle aBezier = new Geom2d_BezierCurve(aPoles); + + Geom2dLProp_CurAndInf2d anOld; + anOld.Perform(aBezier); + ASSERT_TRUE(anOld.IsDone()); + + Geom2dProp_Curve aProp; + aProp.Initialize(aBezier); + const Geom2dProp::CurveAnalysis aNewExt = aProp.FindCurvatureExtrema(); + const Geom2dProp::CurveAnalysis aNewInfl = aProp.FindInflections(); + ASSERT_TRUE(aNewExt.IsDone); + ASSERT_TRUE(aNewInfl.IsDone); + + const int aNewTotal = aNewExt.Points.Length() + aNewInfl.Points.Length(); + EXPECT_EQ(aNewTotal, anOld.NbPoints()); +} + +TEST(Geom2dProp_VsCurAndInf2dTest, Bezier_Quadratic_NoInflections) +{ + NCollection_Array1 aPoles(1, 3); + aPoles(1) = gp_Pnt2d(0.0, 0.0); + aPoles(2) = gp_Pnt2d(2.0, 4.0); + aPoles(3) = gp_Pnt2d(4.0, 0.0); + occ::handle aBezier = new Geom2d_BezierCurve(aPoles); + + Geom2dLProp_CurAndInf2d anOld; + anOld.PerformInf(aBezier); + ASSERT_TRUE(anOld.IsDone()); + + Geom2dProp_Curve aProp; + aProp.Initialize(aBezier); + const Geom2dProp::CurveAnalysis aNew = aProp.FindInflections(); + ASSERT_TRUE(aNew.IsDone); + + EXPECT_EQ(aNew.Points.Length(), anOld.NbPoints()); +} + +TEST(Geom2dProp_VsCurAndInf2dTest, Bezier_HighDegree_FullPerform) +{ + NCollection_Array1 aPoles(1, 6); + aPoles(1) = gp_Pnt2d(0.0, 0.0); + aPoles(2) = gp_Pnt2d(1.0, 3.0); + aPoles(3) = gp_Pnt2d(2.0, -1.0); + aPoles(4) = gp_Pnt2d(3.0, 2.0); + aPoles(5) = gp_Pnt2d(4.0, -2.0); + aPoles(6) = gp_Pnt2d(5.0, 1.0); + occ::handle aBezier = new Geom2d_BezierCurve(aPoles); + + Geom2dLProp_CurAndInf2d anOld; + anOld.Perform(aBezier); + ASSERT_TRUE(anOld.IsDone()); + + Geom2dProp_Curve aProp; + aProp.Initialize(aBezier); + const Geom2dProp::CurveAnalysis aNewExt = aProp.FindCurvatureExtrema(); + const Geom2dProp::CurveAnalysis aNewInfl = aProp.FindInflections(); + ASSERT_TRUE(aNewExt.IsDone); + ASSERT_TRUE(aNewInfl.IsDone); + + const int aNewTotal = aNewExt.Points.Length() + aNewInfl.Points.Length(); + EXPECT_EQ(aNewTotal, anOld.NbPoints()); +} + +// ============================================================================ +// BSpline - numeric with C3 interval subdivision +// ============================================================================ + +TEST(Geom2dProp_VsCurAndInf2dTest, BSpline_Degree4_FullPerform) +{ + NCollection_Array1 aPoles(1, 5); + aPoles(1) = gp_Pnt2d(0.0, 0.0); + aPoles(2) = gp_Pnt2d(1.0, 3.0); + aPoles(3) = gp_Pnt2d(2.0, -1.0); + aPoles(4) = gp_Pnt2d(3.0, 2.0); + aPoles(5) = gp_Pnt2d(4.0, 0.0); + + NCollection_Array1 aKnots(1, 2); + aKnots(1) = 0.0; + aKnots(2) = 1.0; + + NCollection_Array1 aMults(1, 2); + aMults(1) = 5; + aMults(2) = 5; + + occ::handle aBSpline = new Geom2d_BSplineCurve(aPoles, aKnots, aMults, 4); + + Geom2dLProp_CurAndInf2d anOld; + anOld.Perform(aBSpline); + ASSERT_TRUE(anOld.IsDone()); + + Geom2dProp_Curve aProp; + aProp.Initialize(aBSpline); + const Geom2dProp::CurveAnalysis aNewExt = aProp.FindCurvatureExtrema(); + const Geom2dProp::CurveAnalysis aNewInfl = aProp.FindInflections(); + ASSERT_TRUE(aNewExt.IsDone); + ASSERT_TRUE(aNewInfl.IsDone); + + const int aNewTotal = aNewExt.Points.Length() + aNewInfl.Points.Length(); + EXPECT_EQ(aNewTotal, anOld.NbPoints()); +} + +TEST(Geom2dProp_VsCurAndInf2dTest, BSpline_Cubic_Extrema) +{ + NCollection_Array1 aPoles(1, 6); + aPoles(1) = gp_Pnt2d(0.0, 0.0); + aPoles(2) = gp_Pnt2d(1.0, 3.0); + aPoles(3) = gp_Pnt2d(2.0, 1.0); + aPoles(4) = gp_Pnt2d(3.0, 4.0); + aPoles(5) = gp_Pnt2d(4.0, 2.0); + aPoles(6) = gp_Pnt2d(5.0, 0.0); + + NCollection_Array1 aKnots(1, 4); + aKnots(1) = 0.0; + aKnots(2) = 0.33; + aKnots(3) = 0.66; + aKnots(4) = 1.0; + + NCollection_Array1 aMults(1, 4); + aMults(1) = 4; + aMults(2) = 1; + aMults(3) = 1; + aMults(4) = 4; + + occ::handle aBSpline = new Geom2d_BSplineCurve(aPoles, aKnots, aMults, 3); + + Geom2dLProp_CurAndInf2d anOld; + anOld.PerformCurExt(aBSpline); + ASSERT_TRUE(anOld.IsDone()); + + Geom2dProp_Curve aProp; + aProp.Initialize(aBSpline); + const Geom2dProp::CurveAnalysis aNew = aProp.FindCurvatureExtrema(); + ASSERT_TRUE(aNew.IsDone); + + compareExtrema(aNew, anOld); +} + +TEST(Geom2dProp_VsCurAndInf2dTest, BSpline_Cubic_Inflections) +{ + NCollection_Array1 aPoles(1, 6); + aPoles(1) = gp_Pnt2d(0.0, 0.0); + aPoles(2) = gp_Pnt2d(1.0, 3.0); + aPoles(3) = gp_Pnt2d(2.0, 1.0); + aPoles(4) = gp_Pnt2d(3.0, 4.0); + aPoles(5) = gp_Pnt2d(4.0, 2.0); + aPoles(6) = gp_Pnt2d(5.0, 0.0); + + NCollection_Array1 aKnots(1, 4); + aKnots(1) = 0.0; + aKnots(2) = 0.33; + aKnots(3) = 0.66; + aKnots(4) = 1.0; + + NCollection_Array1 aMults(1, 4); + aMults(1) = 4; + aMults(2) = 1; + aMults(3) = 1; + aMults(4) = 4; + + occ::handle aBSpline = new Geom2d_BSplineCurve(aPoles, aKnots, aMults, 3); + + Geom2dLProp_CurAndInf2d anOld; + anOld.PerformInf(aBSpline); + ASSERT_TRUE(anOld.IsDone()); + + Geom2dProp_Curve aProp; + aProp.Initialize(aBSpline); + const Geom2dProp::CurveAnalysis aNew = aProp.FindInflections(); + ASSERT_TRUE(aNew.IsDone); + + compareExtrema(aNew, anOld); +} + +TEST(Geom2dProp_VsCurAndInf2dTest, BSpline_LowContinuity_FullPerform) +{ + NCollection_Array1 aPoles(1, 5); + aPoles(1) = gp_Pnt2d(0.0, 0.0); + aPoles(2) = gp_Pnt2d(1.0, 3.0); + aPoles(3) = gp_Pnt2d(2.0, 1.0); + aPoles(4) = gp_Pnt2d(3.0, 3.0); + aPoles(5) = gp_Pnt2d(4.0, 0.0); + + NCollection_Array1 aKnots(1, 4); + aKnots(1) = 0.0; + aKnots(2) = 0.33; + aKnots(3) = 0.66; + aKnots(4) = 1.0; + + NCollection_Array1 aMults(1, 4); + aMults(1) = 3; + aMults(2) = 1; + aMults(3) = 1; + aMults(4) = 3; + + occ::handle aBSpline = new Geom2d_BSplineCurve(aPoles, aKnots, aMults, 2); + + Geom2dLProp_CurAndInf2d anOld; + anOld.Perform(aBSpline); + ASSERT_TRUE(anOld.IsDone()); + + Geom2dProp_Curve aProp; + aProp.Initialize(aBSpline); + const Geom2dProp::CurveAnalysis aNewExt = aProp.FindCurvatureExtrema(); + const Geom2dProp::CurveAnalysis aNewInfl = aProp.FindInflections(); + ASSERT_TRUE(aNewExt.IsDone); + ASSERT_TRUE(aNewInfl.IsDone); + + const int aNewTotal = aNewExt.Points.Length() + aNewInfl.Points.Length(); + + // The new C3-subdivision solver analyzes each smooth interval independently, + // finding inflection points near knots that the old global solver misses. + // Verify the new API finds at least as many points as the old. + EXPECT_GE(aNewTotal, anOld.NbPoints()); + + // Verify all old points are found by the new API. + for (int i = 1; i <= anOld.NbPoints(); ++i) + { + const double anOldParam = anOld.Parameter(i); + bool aFound = false; + for (int j = 0; j < aNewExt.Points.Length(); ++j) + { + if (std::abs(aNewExt.Points[j].Parameter - anOldParam) < 1.0e-3) + { + aFound = true; + break; + } + } + if (!aFound) + { + for (int j = 0; j < aNewInfl.Points.Length(); ++j) + { + if (std::abs(aNewInfl.Points[j].Parameter - anOldParam) < 1.0e-3) + { + aFound = true; + break; + } + } + } + EXPECT_TRUE(aFound) << "Old point at param=" << anOldParam << " not found in new results"; + } +} + +// ============================================================================ +// Trimmed curve - should work through unwrapping +// ============================================================================ + +TEST(Geom2dProp_VsCurAndInf2dTest, TrimmedEllipse_Extrema) +{ + gp_Elips2d anElips(gp_Ax2d(gp_Pnt2d(0.0, 0.0), gp_Dir2d(1.0, 0.0)), 10.0, 5.0); + occ::handle anEllipse = new Geom2d_Ellipse(anElips); + occ::handle aTrimmed = new Geom2d_TrimmedCurve(anEllipse, 0.0, M_PI); + + Geom2dLProp_CurAndInf2d anOld; + anOld.PerformCurExt(aTrimmed); + ASSERT_TRUE(anOld.IsDone()); + + Geom2dProp_Curve aProp; + aProp.Initialize(aTrimmed); + const Geom2dProp::CurveAnalysis aNew = aProp.FindCurvatureExtrema(); + ASSERT_TRUE(aNew.IsDone); + + compareExtrema(aNew, anOld, 1.0e-6); +} + +TEST(Geom2dProp_VsCurAndInf2dTest, TrimmedBezier_FullPerform) +{ + NCollection_Array1 aPoles(1, 4); + aPoles(1) = gp_Pnt2d(0.0, 0.0); + aPoles(2) = gp_Pnt2d(1.0, 3.0); + aPoles(3) = gp_Pnt2d(3.0, -1.0); + aPoles(4) = gp_Pnt2d(4.0, 1.0); + occ::handle aBezier = new Geom2d_BezierCurve(aPoles); + occ::handle aTrimmed = new Geom2d_TrimmedCurve(aBezier, 0.1, 0.9); + + Geom2dLProp_CurAndInf2d anOld; + anOld.Perform(aTrimmed); + ASSERT_TRUE(anOld.IsDone()); + + Geom2dProp_Curve aProp; + aProp.Initialize(aTrimmed); + const Geom2dProp::CurveAnalysis aNewExt = aProp.FindCurvatureExtrema(); + const Geom2dProp::CurveAnalysis aNewInfl = aProp.FindInflections(); + ASSERT_TRUE(aNewExt.IsDone); + ASSERT_TRUE(aNewInfl.IsDone); + + const int aNewTotal = aNewExt.Points.Length() + aNewInfl.Points.Length(); + EXPECT_EQ(aNewTotal, anOld.NbPoints()); +} + +// ============================================================================ +// Offset curve - numeric +// ============================================================================ + +TEST(Geom2dProp_VsCurAndInf2dTest, OffsetEllipse_Extrema) +{ + gp_Elips2d anElips(gp_Ax2d(gp_Pnt2d(0.0, 0.0), gp_Dir2d(1.0, 0.0)), 10.0, 5.0); + occ::handle anEllipse = new Geom2d_Ellipse(anElips); + occ::handle anOffset = new Geom2d_OffsetCurve(anEllipse, 1.0); + + Geom2dLProp_CurAndInf2d anOld; + anOld.PerformCurExt(anOffset); + ASSERT_TRUE(anOld.IsDone()); + + Geom2dProp_Curve aProp; + aProp.Initialize(anOffset); + const Geom2dProp::CurveAnalysis aNew = aProp.FindCurvatureExtrema(); + ASSERT_TRUE(aNew.IsDone); + + compareExtrema(aNew, anOld); +} + +TEST(Geom2dProp_VsCurAndInf2dTest, OffsetEllipse_Inflections) +{ + gp_Elips2d anElips(gp_Ax2d(gp_Pnt2d(0.0, 0.0), gp_Dir2d(1.0, 0.0)), 10.0, 5.0); + occ::handle anEllipse = new Geom2d_Ellipse(anElips); + occ::handle anOffset = new Geom2d_OffsetCurve(anEllipse, 1.0); + + Geom2dLProp_CurAndInf2d anOld; + anOld.PerformInf(anOffset); + ASSERT_TRUE(anOld.IsDone()); + + Geom2dProp_Curve aProp; + aProp.Initialize(anOffset); + const Geom2dProp::CurveAnalysis aNew = aProp.FindInflections(); + ASSERT_TRUE(aNew.IsDone); + + EXPECT_EQ(aNew.Points.Length(), anOld.NbPoints()); +} + +TEST(Geom2dProp_VsCurAndInf2dTest, OffsetCircle_NoExtrema) +{ + gp_Circ2d aCirc(gp_Ax2d(gp_Pnt2d(0.0, 0.0), gp_Dir2d(1.0, 0.0)), 5.0); + occ::handle aCircle = new Geom2d_Circle(aCirc); + occ::handle anOffset = new Geom2d_OffsetCurve(aCircle, 2.0); + + Geom2dLProp_CurAndInf2d anOld; + anOld.PerformCurExt(anOffset); + ASSERT_TRUE(anOld.IsDone()); + + Geom2dProp_Curve aProp; + aProp.Initialize(anOffset); + const Geom2dProp::CurveAnalysis aNew = aProp.FindCurvatureExtrema(); + ASSERT_TRUE(aNew.IsDone); + + EXPECT_EQ(aNew.Points.Length(), anOld.NbPoints()); +} diff --git a/src/ModelingData/TKG2d/Geom2dProp/FILES.cmake b/src/ModelingData/TKG2d/Geom2dProp/FILES.cmake new file mode 100644 index 0000000000..b8c17d4bc8 --- /dev/null +++ b/src/ModelingData/TKG2d/Geom2dProp/FILES.cmake @@ -0,0 +1,25 @@ +# Source files for Geom2dProp package +set(OCCT_Geom2dProp_FILES_LOCATION "${CMAKE_CURRENT_LIST_DIR}") + +set(OCCT_Geom2dProp_FILES + Geom2dProp.hxx + Geom2dProp.cxx + Geom2dProp_BezierCurve.hxx + Geom2dProp_BezierCurve.cxx + Geom2dProp_BSplineCurve.hxx + Geom2dProp_BSplineCurve.cxx + Geom2dProp_Circle.hxx + Geom2dProp_Curve.hxx + Geom2dProp_Curve.cxx + Geom2dProp_Ellipse.hxx + Geom2dProp_Ellipse.cxx + Geom2dProp_Hyperbola.hxx + Geom2dProp_Hyperbola.cxx + Geom2dProp_Line.hxx + Geom2dProp_OffsetCurve.hxx + Geom2dProp_OffsetCurve.cxx + Geom2dProp_OtherCurve.hxx + Geom2dProp_OtherCurve.cxx + Geom2dProp_Parabola.hxx + Geom2dProp_Parabola.cxx +) diff --git a/src/ModelingData/TKG2d/Geom2dProp/Geom2dProp.cxx b/src/ModelingData/TKG2d/Geom2dProp/Geom2dProp.cxx new file mode 100644 index 0000000000..9fdd1d8f12 --- /dev/null +++ b/src/ModelingData/TKG2d/Geom2dProp/Geom2dProp.cxx @@ -0,0 +1,124 @@ +// Copyright (c) 2025 OPEN CASCADE SAS +// +// This file is part of Open CASCADE Technology software library. +// +// This library is free software; you can redistribute it and/or modify it under +// the terms of the GNU Lesser General Public License version 2.1 as published +// by the Free Software Foundation, with special exception defined in the file +// OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT +// distribution for complete text of the license and disclaimer of any warranty. +// +// Alternatively, this file may be used under the terms of Open CASCADE +// commercial license or contractual agreement. + +#include + +#include + +//================================================================================================== + +Geom2dProp::TangentResult Geom2dProp::ComputeTangent(const gp_Vec2d& theD1, + const gp_Vec2d& theD2, + const gp_Vec2d& theD3, + const double theTol) +{ + const double aTol2 = theTol * theTol; + + // Try first derivative + if (theD1.SquareMagnitude() > aTol2) + { + return {gp_Dir2d(theD1), true}; + } + + // Try second derivative + if (theD2.SquareMagnitude() > aTol2) + { + return {gp_Dir2d(theD2), true}; + } + + // Try third derivative + if (theD3.SquareMagnitude() > aTol2) + { + return {gp_Dir2d(theD3), true}; + } + + return {{}, false}; +} + +//================================================================================================== + +Geom2dProp::CurvatureResult Geom2dProp::ComputeCurvature(const gp_Vec2d& theD1, + const gp_Vec2d& theD2, + const double theTol) +{ + const double aTol2 = theTol * theTol; + const double aDD1 = theD1.SquareMagnitude(); + + // If first derivative is null, curvature is infinite (singular point). + if (aDD1 <= aTol2) + { + return {0.0, true, true}; + } + + const double aDD2 = theD2.SquareMagnitude(); + + // If second derivative is null, curvature is zero. + if (aDD2 <= aTol2) + { + return {0.0, true, false}; + } + + // Cross magnitude squared: |D1 x D2|^2 + const double aN = theD1.CrossSquareMagnitude(theD2); + + // If D1 and D2 are collinear, curvature is zero. + const double aT = aN / aDD1 / aDD2; + if (aT <= aTol2) + { + return {0.0, true, false}; + } + + // Curvature = |D1 x D2| / |D1|^3 + const double aCurvature = std::sqrt(aN) / aDD1 / std::sqrt(aDD1); + return {aCurvature, true, false}; +} + +//================================================================================================== + +Geom2dProp::NormalResult Geom2dProp::ComputeNormal(const gp_Vec2d& theD1, + const gp_Vec2d& theD2, + const double theTol) +{ + // First compute curvature to check if normal is defined. + const CurvatureResult aCurvRes = ComputeCurvature(theD1, theD2, theTol); + if (!aCurvRes.IsDefined || aCurvRes.IsInfinite || std::abs(aCurvRes.Value) <= theTol) + { + return {{}, false}; + } + + // Normal = D2 * (D1.D1) - D1 * (D1.D2) + // This is equivalent to D1 x (D2 x D1) in 2D using the vector triple product identity. + const gp_Vec2d aNorm = theD2 * theD1.Dot(theD1) - theD1 * theD1.Dot(theD2); + return {gp_Dir2d(aNorm), true}; +} + +//================================================================================================== + +Geom2dProp::CentreResult Geom2dProp::ComputeCentreOfCurvature(const gp_Pnt2d& thePnt, + const gp_Vec2d& theD1, + const gp_Vec2d& theD2, + const double theTol) +{ + const CurvatureResult aCurvRes = ComputeCurvature(theD1, theD2, theTol); + if (!aCurvRes.IsDefined || aCurvRes.IsInfinite || std::abs(aCurvRes.Value) <= theTol) + { + return {{}, false}; + } + + // Normal vector (unnormalized) = D2 * (D1.D1) - D1 * (D1.D2) + gp_Vec2d aNorm = theD2 * theD1.Dot(theD1) - theD1 * theD1.Dot(theD2); + aNorm.Normalize(); + aNorm.Divide(aCurvRes.Value); + + return {thePnt.Translated(aNorm), true}; +} diff --git a/src/ModelingData/TKG2d/Geom2dProp/Geom2dProp.hxx b/src/ModelingData/TKG2d/Geom2dProp/Geom2dProp.hxx new file mode 100644 index 0000000000..da24c60156 --- /dev/null +++ b/src/ModelingData/TKG2d/Geom2dProp/Geom2dProp.hxx @@ -0,0 +1,129 @@ +// Copyright (c) 2025 OPEN CASCADE SAS +// +// This file is part of Open CASCADE Technology software library. +// +// This library is free software; you can redistribute it and/or modify it under +// the terms of the GNU Lesser General Public License version 2.1 as published +// by the Free Software Foundation, with special exception defined in the file +// OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT +// distribution for complete text of the license and disclaimer of any warranty. +// +// Alternatively, this file may be used under the terms of Open CASCADE +// commercial license or contractual agreement. + +#ifndef _Geom2dProp_HeaderFile +#define _Geom2dProp_HeaderFile + +#include +#include +#include +#include +#include + +//! @brief Namespace containing result structures and free functions for 2D curve +//! differential property computation. +//! +//! Provides lightweight result structures with explicit validity flags instead of +//! exception-based APIs, and geometry-agnostic free functions that compute local +//! differential properties from derivative vectors. +namespace Geom2dProp +{ + +//! Result of tangent direction computation. +struct TangentResult +{ + gp_Dir2d Direction; //!< Tangent direction (valid only when IsDefined is true) + bool IsDefined = false; //!< True if the tangent is well-defined +}; + +//! Result of curvature computation. +struct CurvatureResult +{ + double Value = 0.0; //!< Curvature value (valid only when IsDefined is true) + bool IsDefined = false; //!< True if curvature could be computed + bool IsInfinite = false; //!< True if first derivative is null (singular point) +}; + +//! Result of normal direction computation. +struct NormalResult +{ + gp_Dir2d Direction; //!< Normal direction (valid only when IsDefined is true) + bool IsDefined = false; //!< True if the normal is well-defined +}; + +//! Result of centre of curvature computation. +struct CentreResult +{ + gp_Pnt2d Centre; //!< Centre of curvature (valid only when IsDefined is true) + bool IsDefined = false; //!< True if the centre is well-defined +}; + +//! Type of a special curve point (curvature extremum or inflection). +enum class CIType +{ + Inflection, //!< Inflection point (curvature changes sign) + MinCurvature, //!< Local minimum of the radius of curvature (maximum of |curvature|) + MaxCurvature //!< Local maximum of the radius of curvature (minimum of |curvature|) +}; + +//! A special point on a curve with its parameter and type. +struct CurveSpecialPoint +{ + double Parameter = 0.0; //!< Curve parameter + CIType Type = CIType::Inflection; //!< Point type +}; + +//! Result of global curve analysis (curvature extrema and inflection points). +struct CurveAnalysis +{ + NCollection_DynamicArray Points; //!< Special points sorted by parameter + bool IsDone = false; //!< True if analysis completed +}; + +//! Compute tangent direction from derivative vectors. +//! Tries D1 first; if D1 magnitude^2 <= theTol^2, tries D2, then D3. +//! @param[in] theD1 first derivative vector +//! @param[in] theD2 second derivative vector +//! @param[in] theD3 third derivative vector +//! @param[in] theTol linear tolerance for zero-vector detection +//! @return tangent result with validity flag +Standard_EXPORT TangentResult ComputeTangent(const gp_Vec2d& theD1, + const gp_Vec2d& theD2, + const gp_Vec2d& theD3, + double theTol); + +//! Compute curvature from first and second derivative vectors. +//! Curvature = |D1 x D2| / |D1|^3 +//! @param[in] theD1 first derivative vector +//! @param[in] theD2 second derivative vector +//! @param[in] theTol linear tolerance for zero-vector detection +//! @return curvature result with validity and infinity flags +Standard_EXPORT CurvatureResult ComputeCurvature(const gp_Vec2d& theD1, + const gp_Vec2d& theD2, + double theTol); + +//! Compute normal direction from first and second derivative vectors. +//! Normal = D1 x (D2 x D1) (normalized), perpendicular to tangent pointing toward center. +//! @param[in] theD1 first derivative vector +//! @param[in] theD2 second derivative vector +//! @param[in] theTol linear tolerance for zero-vector detection +//! @return normal result with validity flag +Standard_EXPORT NormalResult ComputeNormal(const gp_Vec2d& theD1, + const gp_Vec2d& theD2, + double theTol); + +//! Compute centre of curvature from point and derivative vectors. +//! Centre = Point + Normal / Curvature +//! @param[in] thePnt point on the curve +//! @param[in] theD1 first derivative vector +//! @param[in] theD2 second derivative vector +//! @param[in] theTol linear tolerance for zero-vector detection +//! @return centre result with validity flag +Standard_EXPORT CentreResult ComputeCentreOfCurvature(const gp_Pnt2d& thePnt, + const gp_Vec2d& theD1, + const gp_Vec2d& theD2, + double theTol); + +} // namespace Geom2dProp + +#endif // _Geom2dProp_HeaderFile diff --git a/src/ModelingData/TKG2d/Geom2dProp/Geom2dProp_BSplineCurve.cxx b/src/ModelingData/TKG2d/Geom2dProp/Geom2dProp_BSplineCurve.cxx new file mode 100644 index 0000000000..0f8c94e12a --- /dev/null +++ b/src/ModelingData/TKG2d/Geom2dProp/Geom2dProp_BSplineCurve.cxx @@ -0,0 +1,429 @@ +// Copyright (c) 2025 OPEN CASCADE SAS +// +// This file is part of Open CASCADE Technology software library. +// +// This library is free software; you can redistribute it and/or modify it under +// the terms of the GNU Lesser General Public License version 2.1 as published +// by the Free Software Foundation, with special exception defined in the file +// OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT +// distribution for complete text of the license and disclaimer of any warranty. +// +// Alternatively, this file may be used under the terms of Open CASCADE +// commercial license or contractual agreement. + +#include + +#include +#include +#include +#include +#include +#include + +#include + +namespace +{ + +constexpr double THE_CURVATURE_DERIV_COEFF = 3.0; //!< Coefficient in d(KC)/dU formula +constexpr double THE_DIFF_STEP_DIVISOR = 100.0; //!< Divisor for numerical differentiation step +constexpr double THE_D2_MAGNITUDE_THRESHOLD = 1.0e-4; //!< Threshold for second derivative magnitude +constexpr double THE_EPSILON_SCALE = 1.0e-4; //!< Scale factor for epsilon relative to domain +constexpr int THE_EXTREMA_NB_SAMPLES = 100; //!< Number of samples for curvature extrema search +constexpr int THE_INFLECTION_NB_SAMPLES = 30; //!< Number of samples for inflection search +constexpr double THE_INFLECTION_TOLERANCE = 1.0e-6; //!< Tolerance for inflection point finding + +//! Function for finding curvature extrema: F = d(curvature)/dU = 0 +class FuncCurExt +{ +public: + FuncCurExt(const Geom2dAdaptor_Curve* theCurve, const double theTol) + : myCurve(theCurve), + myEpsX(theTol) + { + } + + bool Value(const double X, double& F) + { + gp_Pnt2d aP; + gp_Vec2d aV1, aV2, aV3; + myCurve->D3(X, aP, aV1, aV2, aV3); + + const double aCPV1V2 = aV1.Crossed(aV2); + const double aCPV1V3 = aV1.Crossed(aV3); + const double aV1V2 = aV1.Dot(aV2); + const double aV1V1 = aV1.SquareMagnitude(); + const double aNV1 = std::sqrt(aV1V1); + const double aV13 = aV1V1 * aNV1; + const double aV15 = aV13 * aV1V1; + + if (aV15 < gp::Resolution()) + { + return false; + } + + F = aCPV1V3 / aV13 - THE_CURVATURE_DERIV_COEFF * aCPV1V2 * aV1V2 / aV15; + return true; + } + + bool Values(const double X, double& F, double& D) + { + double aDx = myEpsX / THE_DIFF_STEP_DIVISOR; + if (X + aDx > myCurve->LastParameter()) + { + aDx = -aDx; + } + + Value(X, F); + double aF2; + Value(X + aDx, aF2); + D = (aF2 - F) / aDx; + return true; + } + + bool IsMinKC(const double X) const + { + gp_Pnt2d aP; + gp_Vec2d aV1, aV2, aV3; + + myCurve->D3(X, aP, aV1, aV2, aV3); + const double aV1V1 = aV1.SquareMagnitude(); + const double aNV1 = std::sqrt(aV1V1); + const double aV13 = aV1V1 * aNV1; + + if (aV13 < gp::Resolution()) + { + return false; + } + const double aKC = aV1.Crossed(aV2) / aV13; + + double aDx = myEpsX; + if (X + aDx > myCurve->LastParameter()) + { + aDx = -aDx; + } + + myCurve->D3(X + aDx, aP, aV1, aV2, aV3); + const double aV1V1n = aV1.SquareMagnitude(); + const double aNV1n = std::sqrt(aV1V1n); + const double aV13n = aV1V1n * aNV1n; + + if (aV13n < gp::Resolution()) + { + return false; + } + const double aKP = aV1.Crossed(aV2) / aV13n; + + return std::abs(aKC) > std::abs(aKP); + } + +private: + const Geom2dAdaptor_Curve* myCurve; + double myEpsX; +}; + +//! Function for finding inflection points: F = (V1^V2) / (||V1|| * ||V2||) = 0 +class FuncCurNul +{ +public: + FuncCurNul(const Geom2dAdaptor_Curve* theCurve) + : myCurve(theCurve) + { + } + + bool Value(const double X, double& F) + { + double aD; + return Values(X, F, aD); + } + + bool Values(const double X, double& F, double& D) + { + gp_Pnt2d aP; + gp_Vec2d aV1, aV2, aV3; + myCurve->D3(X, aP, aV1, aV2, aV3); + + const double aCP1 = aV1.Crossed(aV2); + const double aCP2 = aV1.Crossed(aV3); + const double aV1V2 = aV1.Dot(aV2); + const double aV2V3 = aV2.Dot(aV3); + const double aNV1 = aV1.Magnitude(); + const double aNV2 = aV2.Magnitude(); + + F = 0.0; + D = 0.0; + + if (aNV2 < THE_D2_MAGNITUDE_THRESHOLD) + { + return true; + } + if (aNV1 * aNV2 < gp::Resolution()) + { + return false; + } + + F = aCP1 / (aNV1 * aNV2); + D = (aCP2 - aCP1 * aV1V2 / (aNV1 * aNV1) - aCP1 * aV2V3 / (aNV2 * aNV2)) / (aNV1 * aNV2); + return true; + } + +private: + const Geom2dAdaptor_Curve* myCurve; +}; + +//! Perform numeric curvature extrema finding on a curve interval. +void numericCurvatureExtrema(const Geom2dAdaptor_Curve* theCurve, + const double theUMin, + const double theUMax, + Geom2dProp::CurveAnalysis& theResult) +{ + const double aEpsH = THE_EPSILON_SCALE * (theUMax - theUMin); + + FuncCurExt aFunc(theCurve, aEpsH); + + MathRoot::MultipleConfig aConfig; + aConfig.NbSamples = THE_EXTREMA_NB_SAMPLES; + aConfig.XTolerance = aEpsH; + aConfig.FTolerance = aEpsH; + + MathRoot::MultipleResult aRoots = + MathRoot::FindAllRootsWithDerivative(aFunc, theUMin, theUMax, aConfig); + + if (aRoots.IsDone()) + { + for (int j = 0; j < aRoots.NbRoots(); ++j) + { + double aParam = aRoots[j]; + // Refine the solution. + MathUtils::Config aBrentCfg; + aBrentCfg.XTolerance = Precision::PConfusion(); + aBrentCfg.FTolerance = Precision::PConfusion(); + auto aBrent = MathRoot::Brent(aFunc, aParam - aEpsH, aParam + aEpsH, aBrentCfg); + if (aBrent.IsDone() && aBrent.Root.has_value()) + { + aParam = *aBrent.Root; + } + const bool aIsMin = aFunc.IsMinKC(aParam); + const Geom2dProp::CIType aType = + aIsMin ? Geom2dProp::CIType::MinCurvature : Geom2dProp::CIType::MaxCurvature; + theResult.Points.Append({aParam, aType}); + } + } + else + { + theResult.IsDone = false; + } +} + +//! Perform numeric inflection finding on a curve interval. +void numericInflections(const Geom2dAdaptor_Curve* theCurve, + const double theUMin, + const double theUMax, + Geom2dProp::CurveAnalysis& theResult) +{ + FuncCurNul aFunc(theCurve); + + MathRoot::MultipleConfig aConfig; + aConfig.NbSamples = THE_INFLECTION_NB_SAMPLES; + aConfig.XTolerance = THE_INFLECTION_TOLERANCE; + aConfig.FTolerance = THE_INFLECTION_TOLERANCE; + + MathRoot::MultipleResult aRoots = MathRoot::FindAllRoots(aFunc, theUMin, theUMax, aConfig); + + if (aRoots.IsDone()) + { + for (int j = 0; j < aRoots.NbRoots(); ++j) + { + theResult.Points.Append({aRoots[j], Geom2dProp::CIType::Inflection}); + } + } + else + { + theResult.IsDone = false; + } +} + +//! Remove duplicate points that may appear at shared interval boundaries. +//! Points are considered duplicates if their parameters are within theTol. +void removeDuplicatePoints(Geom2dProp::CurveAnalysis& theResult, const double theTol) +{ + const int aNbPts = theResult.Points.Size(); + if (aNbPts <= 1) + { + return; + } + + // Pre-check: detect if any duplicates exist before allocating. + bool aHasDuplicates = false; + for (int i = 1; i < aNbPts && !aHasDuplicates; ++i) + { + for (int j = 0; j < i; ++j) + { + if (std::abs(theResult.Points[i].Parameter - theResult.Points[j].Parameter) < theTol) + { + aHasDuplicates = true; + break; + } + } + } + if (!aHasDuplicates) + { + return; + } + + NCollection_DynamicArray aFiltered; + aFiltered.Append(theResult.Points[0]); + for (int i = 1; i < aNbPts; ++i) + { + bool aIsDuplicate = false; + for (int j = static_cast(aFiltered.Size()) - 1; j >= 0; --j) + { + if (std::abs(theResult.Points[i].Parameter - aFiltered[j].Parameter) < theTol) + { + aIsDuplicate = true; + break; + } + } + if (!aIsDuplicate) + { + aFiltered.Append(theResult.Points[i]); + } + } + + theResult.Points = std::move(aFiltered); +} + +} // namespace + +//================================================================================================== + +Geom2dProp::TangentResult Geom2dProp_BSplineCurve::Tangent(const double theParam, + const double theTol) const +{ + if (myAdaptor == nullptr) + { + return {{}, false}; + } + gp_Pnt2d aPnt; + gp_Vec2d aD1, aD2, aD3; + myAdaptor->D3(theParam, aPnt, aD1, aD2, aD3); + return Geom2dProp::ComputeTangent(aD1, aD2, aD3, theTol); +} + +//================================================================================================== + +Geom2dProp::CurvatureResult Geom2dProp_BSplineCurve::Curvature(const double theParam, + const double theTol) const +{ + if (myAdaptor == nullptr) + { + return {0.0, false, false}; + } + gp_Pnt2d aPnt; + gp_Vec2d aD1, aD2; + myAdaptor->D2(theParam, aPnt, aD1, aD2); + return Geom2dProp::ComputeCurvature(aD1, aD2, theTol); +} + +//================================================================================================== + +Geom2dProp::NormalResult Geom2dProp_BSplineCurve::Normal(const double theParam, + const double theTol) const +{ + if (myAdaptor == nullptr) + { + return {{}, false}; + } + gp_Pnt2d aPnt; + gp_Vec2d aD1, aD2; + myAdaptor->D2(theParam, aPnt, aD1, aD2); + return Geom2dProp::ComputeNormal(aD1, aD2, theTol); +} + +//================================================================================================== + +Geom2dProp::CentreResult Geom2dProp_BSplineCurve::CentreOfCurvature(const double theParam, + const double theTol) const +{ + if (myAdaptor == nullptr) + { + return {{}, false}; + } + gp_Pnt2d aPnt; + gp_Vec2d aD1, aD2; + myAdaptor->D2(theParam, aPnt, aD1, aD2); + return Geom2dProp::ComputeCentreOfCurvature(aPnt, aD1, aD2, theTol); +} + +//================================================================================================== + +Geom2dProp::CurveAnalysis Geom2dProp_BSplineCurve::FindCurvatureExtrema() const +{ + Geom2dProp::CurveAnalysis aResult; + aResult.IsDone = true; + + if (myAdaptor == nullptr) + { + aResult.IsDone = false; + return aResult; + } + + if (myAdaptor->Continuity() >= GeomAbs_C3) + { + numericCurvatureExtrema(myAdaptor, + myAdaptor->FirstParameter(), + myAdaptor->LastParameter(), + aResult); + } + else + { + // Subdivide into C3 intervals. + const int aNbInt = myAdaptor->NbIntervals(GeomAbs_C3); + NCollection_Array1 aParams(1, aNbInt + 1); + myAdaptor->Intervals(aParams, GeomAbs_C3); + for (int i = 1; i <= aNbInt; ++i) + { + numericCurvatureExtrema(myAdaptor, aParams(i), aParams(i + 1), aResult); + } + // Remove duplicate roots that may appear at shared interval boundaries. + const double aEpsH = + THE_EPSILON_SCALE * (myAdaptor->LastParameter() - myAdaptor->FirstParameter()); + removeDuplicatePoints(aResult, aEpsH); + } + + return aResult; +} + +//================================================================================================== + +Geom2dProp::CurveAnalysis Geom2dProp_BSplineCurve::FindInflections() const +{ + Geom2dProp::CurveAnalysis aResult; + aResult.IsDone = true; + + if (myAdaptor == nullptr) + { + aResult.IsDone = false; + return aResult; + } + + if (myAdaptor->Continuity() >= GeomAbs_C3) + { + numericInflections(myAdaptor, myAdaptor->FirstParameter(), myAdaptor->LastParameter(), aResult); + } + else + { + // Subdivide into C3 intervals. + const int aNbInt = myAdaptor->NbIntervals(GeomAbs_C3); + NCollection_Array1 aParams(1, aNbInt + 1); + myAdaptor->Intervals(aParams, GeomAbs_C3); + for (int i = 1; i <= aNbInt; ++i) + { + numericInflections(myAdaptor, aParams(i), aParams(i + 1), aResult); + } + // Remove duplicate roots that may appear at shared interval boundaries. + removeDuplicatePoints(aResult, THE_INFLECTION_TOLERANCE); + } + + return aResult; +} diff --git a/src/ModelingData/TKG2d/Geom2dProp/Geom2dProp_BSplineCurve.hxx b/src/ModelingData/TKG2d/Geom2dProp/Geom2dProp_BSplineCurve.hxx new file mode 100644 index 0000000000..56921b3d3b --- /dev/null +++ b/src/ModelingData/TKG2d/Geom2dProp/Geom2dProp_BSplineCurve.hxx @@ -0,0 +1,76 @@ +// Copyright (c) 2025 OPEN CASCADE SAS +// +// This file is part of Open CASCADE Technology software library. +// +// This library is free software; you can redistribute it and/or modify it under +// the terms of the GNU Lesser General Public License version 2.1 as published +// by the Free Software Foundation, with special exception defined in the file +// OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT +// distribution for complete text of the license and disclaimer of any warranty. +// +// Alternatively, this file may be used under the terms of Open CASCADE +// commercial license or contractual agreement. + +#ifndef _Geom2dProp_BSplineCurve_HeaderFile +#define _Geom2dProp_BSplineCurve_HeaderFile + +#include +#include +#include +#include + +//! @brief Local differential properties for a 2D B-spline curve. +//! +//! Uses numeric root-finding for curvature extrema and inflection points. +//! For B-splines with continuity less than C3, the parameter range is subdivided +//! into C3 intervals for more robust root-finding. +//! +//! @warning The caller must ensure that the adaptor pointer remains valid +//! for the entire lifetime of this object. This class does not manage +//! the adaptor's lifetime. +class Geom2dProp_BSplineCurve +{ +public: + DEFINE_STANDARD_ALLOC + + //! Constructor with adaptor pointer (non-owning). + //! @param theAdaptor the 2D curve adaptor (must wrap a B-spline curve, must not be null) + Geom2dProp_BSplineCurve(const Geom2dAdaptor_Curve* theAdaptor) + : myAdaptor(theAdaptor) + { + } + + //! Non-copyable and non-movable. + Geom2dProp_BSplineCurve(const Geom2dProp_BSplineCurve&) = delete; + Geom2dProp_BSplineCurve& operator=(const Geom2dProp_BSplineCurve&) = delete; + Geom2dProp_BSplineCurve(Geom2dProp_BSplineCurve&&) = delete; + Geom2dProp_BSplineCurve& operator=(Geom2dProp_BSplineCurve&&) = delete; + + //! Returns the adaptor pointer. + const Geom2dAdaptor_Curve* Adaptor() const { return myAdaptor; } + + //! Compute tangent at given parameter. + Standard_EXPORT Geom2dProp::TangentResult Tangent(double theParam, double theTol) const; + + //! Compute curvature at given parameter. + Standard_EXPORT Geom2dProp::CurvatureResult Curvature(double theParam, double theTol) const; + + //! Compute normal at given parameter. + Standard_EXPORT Geom2dProp::NormalResult Normal(double theParam, double theTol) const; + + //! Compute centre of curvature at given parameter. + Standard_EXPORT Geom2dProp::CentreResult CentreOfCurvature(double theParam, double theTol) const; + + //! Find curvature extrema using numeric root-finding. + //! For non-C3 B-splines, subdivides into C3 intervals. + Standard_EXPORT Geom2dProp::CurveAnalysis FindCurvatureExtrema() const; + + //! Find inflection points using numeric root-finding. + //! For non-C3 B-splines, subdivides into C3 intervals. + Standard_EXPORT Geom2dProp::CurveAnalysis FindInflections() const; + +private: + const Geom2dAdaptor_Curve* myAdaptor; +}; + +#endif // _Geom2dProp_BSplineCurve_HeaderFile diff --git a/src/ModelingData/TKG2d/Geom2dProp/Geom2dProp_BezierCurve.cxx b/src/ModelingData/TKG2d/Geom2dProp/Geom2dProp_BezierCurve.cxx new file mode 100644 index 0000000000..155e39cc20 --- /dev/null +++ b/src/ModelingData/TKG2d/Geom2dProp/Geom2dProp_BezierCurve.cxx @@ -0,0 +1,347 @@ +// Copyright (c) 2025 OPEN CASCADE SAS +// +// This file is part of Open CASCADE Technology software library. +// +// This library is free software; you can redistribute it and/or modify it under +// the terms of the GNU Lesser General Public License version 2.1 as published +// by the Free Software Foundation, with special exception defined in the file +// OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT +// distribution for complete text of the license and disclaimer of any warranty. +// +// Alternatively, this file may be used under the terms of Open CASCADE +// commercial license or contractual agreement. + +#include + +#include +#include +#include +#include + +#include + +namespace +{ + +constexpr double THE_CURVATURE_DERIV_COEFF = 3.0; //!< Coefficient in d(KC)/dU formula +constexpr double THE_DIFF_STEP_DIVISOR = 100.0; //!< Divisor for numerical differentiation step +constexpr double THE_D2_MAGNITUDE_THRESHOLD = 1.0e-4; //!< Threshold for second derivative magnitude +constexpr double THE_EPSILON_SCALE = 1.0e-4; //!< Scale factor for epsilon relative to domain +constexpr int THE_EXTREMA_NB_SAMPLES = 100; //!< Number of samples for curvature extrema search +constexpr int THE_INFLECTION_NB_SAMPLES = 30; //!< Number of samples for inflection search +constexpr double THE_INFLECTION_TOLERANCE = 1.0e-6; //!< Tolerance for inflection point finding + +//! Function for finding curvature extrema: F = d(curvature)/dU = 0 +//! KC = (V1^V2) / ||V1||^3 +//! F = d KC / dU +class FuncCurExt +{ +public: + FuncCurExt(const Geom2dAdaptor_Curve* theCurve, const double theTol) + : myCurve(theCurve), + myEpsX(theTol) + { + } + + bool Value(const double X, double& F) + { + gp_Pnt2d aP; + gp_Vec2d aV1, aV2, aV3; + myCurve->D3(X, aP, aV1, aV2, aV3); + + const double aCPV1V2 = aV1.Crossed(aV2); + const double aCPV1V3 = aV1.Crossed(aV3); + const double aV1V2 = aV1.Dot(aV2); + const double aV1V1 = aV1.SquareMagnitude(); + const double aNV1 = std::sqrt(aV1V1); + const double aV13 = aV1V1 * aNV1; + const double aV15 = aV13 * aV1V1; + + if (aV15 < gp::Resolution()) + { + return false; + } + + F = aCPV1V3 / aV13 - THE_CURVATURE_DERIV_COEFF * aCPV1V2 * aV1V2 / aV15; + return true; + } + + bool Values(const double X, double& F, double& D) + { + double aDx = myEpsX / THE_DIFF_STEP_DIVISOR; + if (X + aDx > myCurve->LastParameter()) + { + aDx = -aDx; + } + + Value(X, F); + double aF2; + Value(X + aDx, aF2); + D = (aF2 - F) / aDx; + return true; + } + + //! Test if parameter corresponds to a minimum of the radius of curvature + //! (maximum of |curvature|) by comparison with a neighboring point. + bool IsMinKC(const double X) const + { + gp_Pnt2d aP; + gp_Vec2d aV1, aV2, aV3; + + myCurve->D3(X, aP, aV1, aV2, aV3); + const double aCPV1V2 = aV1.Crossed(aV2); + const double aV1V1 = aV1.SquareMagnitude(); + const double aNV1 = std::sqrt(aV1V1); + const double aV13 = aV1V1 * aNV1; + + if (aV13 < gp::Resolution()) + { + return false; + } + const double aKC = aCPV1V2 / aV13; + + double aDx = myEpsX; + if (X + aDx > myCurve->LastParameter()) + { + aDx = -aDx; + } + + myCurve->D3(X + aDx, aP, aV1, aV2, aV3); + const double aCPV1V2n = aV1.Crossed(aV2); + const double aV1V1n = aV1.SquareMagnitude(); + const double aNV1n = std::sqrt(aV1V1n); + const double aV13n = aV1V1n * aNV1n; + + if (aV13n < gp::Resolution()) + { + return false; + } + const double aKP = aCPV1V2n / aV13n; + + return std::abs(aKC) > std::abs(aKP); + } + +private: + const Geom2dAdaptor_Curve* myCurve; + double myEpsX; +}; + +//! Function for finding inflection points: F = (V1^V2) / (||V1|| * ||V2||) = 0 +class FuncCurNul +{ +public: + FuncCurNul(const Geom2dAdaptor_Curve* theCurve) + : myCurve(theCurve) + { + } + + bool Value(const double X, double& F) + { + double aD; + return Values(X, F, aD); + } + + bool Values(const double X, double& F, double& D) + { + gp_Pnt2d aP; + gp_Vec2d aV1, aV2, aV3; + myCurve->D3(X, aP, aV1, aV2, aV3); + + const double aCP1 = aV1.Crossed(aV2); + const double aCP2 = aV1.Crossed(aV3); + const double aV1V2 = aV1.Dot(aV2); + const double aV2V3 = aV2.Dot(aV3); + const double aNV1 = aV1.Magnitude(); + const double aNV2 = aV2.Magnitude(); + + F = 0.0; + D = 0.0; + + if (aNV2 < THE_D2_MAGNITUDE_THRESHOLD) + { + return true; + } + if (aNV1 * aNV2 < gp::Resolution()) + { + return false; + } + + F = aCP1 / (aNV1 * aNV2); + D = (aCP2 - aCP1 * aV1V2 / (aNV1 * aNV1) - aCP1 * aV2V3 / (aNV2 * aNV2)) / (aNV1 * aNV2); + return true; + } + +private: + const Geom2dAdaptor_Curve* myCurve; +}; + +//! Perform numeric curvature extrema finding on a curve interval. +void numericCurvatureExtrema(const Geom2dAdaptor_Curve* theCurve, + const double theUMin, + const double theUMax, + Geom2dProp::CurveAnalysis& theResult) +{ + const double aEpsH = THE_EPSILON_SCALE * (theUMax - theUMin); + + FuncCurExt aFunc(theCurve, aEpsH); + + MathRoot::MultipleConfig aConfig; + aConfig.NbSamples = THE_EXTREMA_NB_SAMPLES; + aConfig.XTolerance = aEpsH; + aConfig.FTolerance = aEpsH; + + MathRoot::MultipleResult aRoots = + MathRoot::FindAllRootsWithDerivative(aFunc, theUMin, theUMax, aConfig); + + if (aRoots.IsDone()) + { + for (int j = 0; j < aRoots.NbRoots(); ++j) + { + double aParam = aRoots[j]; + // Refine the solution. + MathUtils::Config aBrentCfg; + aBrentCfg.XTolerance = Precision::PConfusion(); + aBrentCfg.FTolerance = Precision::PConfusion(); + auto aBrent = MathRoot::Brent(aFunc, aParam - aEpsH, aParam + aEpsH, aBrentCfg); + if (aBrent.IsDone() && aBrent.Root.has_value()) + { + aParam = *aBrent.Root; + } + const bool aIsMin = aFunc.IsMinKC(aParam); + const Geom2dProp::CIType aType = + aIsMin ? Geom2dProp::CIType::MinCurvature : Geom2dProp::CIType::MaxCurvature; + theResult.Points.Append({aParam, aType}); + } + } + else + { + theResult.IsDone = false; + } +} + +//! Perform numeric inflection finding on a curve interval. +void numericInflections(const Geom2dAdaptor_Curve* theCurve, + const double theUMin, + const double theUMax, + Geom2dProp::CurveAnalysis& theResult) +{ + FuncCurNul aFunc(theCurve); + + MathRoot::MultipleConfig aConfig; + aConfig.NbSamples = THE_INFLECTION_NB_SAMPLES; + aConfig.XTolerance = THE_INFLECTION_TOLERANCE; + aConfig.FTolerance = THE_INFLECTION_TOLERANCE; + + MathRoot::MultipleResult aRoots = MathRoot::FindAllRoots(aFunc, theUMin, theUMax, aConfig); + + if (aRoots.IsDone()) + { + for (int j = 0; j < aRoots.NbRoots(); ++j) + { + theResult.Points.Append({aRoots[j], Geom2dProp::CIType::Inflection}); + } + } + else + { + theResult.IsDone = false; + } +} + +} // namespace + +//================================================================================================== + +Geom2dProp::TangentResult Geom2dProp_BezierCurve::Tangent(const double theParam, + const double theTol) const +{ + if (myAdaptor == nullptr) + { + return {{}, false}; + } + gp_Pnt2d aPnt; + gp_Vec2d aD1, aD2, aD3; + myAdaptor->D3(theParam, aPnt, aD1, aD2, aD3); + return Geom2dProp::ComputeTangent(aD1, aD2, aD3, theTol); +} + +//================================================================================================== + +Geom2dProp::CurvatureResult Geom2dProp_BezierCurve::Curvature(const double theParam, + const double theTol) const +{ + if (myAdaptor == nullptr) + { + return {0.0, false, false}; + } + gp_Pnt2d aPnt; + gp_Vec2d aD1, aD2; + myAdaptor->D2(theParam, aPnt, aD1, aD2); + return Geom2dProp::ComputeCurvature(aD1, aD2, theTol); +} + +//================================================================================================== + +Geom2dProp::NormalResult Geom2dProp_BezierCurve::Normal(const double theParam, + const double theTol) const +{ + if (myAdaptor == nullptr) + { + return {{}, false}; + } + gp_Pnt2d aPnt; + gp_Vec2d aD1, aD2; + myAdaptor->D2(theParam, aPnt, aD1, aD2); + return Geom2dProp::ComputeNormal(aD1, aD2, theTol); +} + +//================================================================================================== + +Geom2dProp::CentreResult Geom2dProp_BezierCurve::CentreOfCurvature(const double theParam, + const double theTol) const +{ + if (myAdaptor == nullptr) + { + return {{}, false}; + } + gp_Pnt2d aPnt; + gp_Vec2d aD1, aD2; + myAdaptor->D2(theParam, aPnt, aD1, aD2); + return Geom2dProp::ComputeCentreOfCurvature(aPnt, aD1, aD2, theTol); +} + +//================================================================================================== + +Geom2dProp::CurveAnalysis Geom2dProp_BezierCurve::FindCurvatureExtrema() const +{ + Geom2dProp::CurveAnalysis aResult; + aResult.IsDone = true; + + if (myAdaptor == nullptr) + { + aResult.IsDone = false; + return aResult; + } + + numericCurvatureExtrema(myAdaptor, + myAdaptor->FirstParameter(), + myAdaptor->LastParameter(), + aResult); + return aResult; +} + +//================================================================================================== + +Geom2dProp::CurveAnalysis Geom2dProp_BezierCurve::FindInflections() const +{ + Geom2dProp::CurveAnalysis aResult; + aResult.IsDone = true; + + if (myAdaptor == nullptr) + { + aResult.IsDone = false; + return aResult; + } + + numericInflections(myAdaptor, myAdaptor->FirstParameter(), myAdaptor->LastParameter(), aResult); + return aResult; +} diff --git a/src/ModelingData/TKG2d/Geom2dProp/Geom2dProp_BezierCurve.hxx b/src/ModelingData/TKG2d/Geom2dProp/Geom2dProp_BezierCurve.hxx new file mode 100644 index 0000000000..8a3363c4d6 --- /dev/null +++ b/src/ModelingData/TKG2d/Geom2dProp/Geom2dProp_BezierCurve.hxx @@ -0,0 +1,72 @@ +// Copyright (c) 2025 OPEN CASCADE SAS +// +// This file is part of Open CASCADE Technology software library. +// +// This library is free software; you can redistribute it and/or modify it under +// the terms of the GNU Lesser General Public License version 2.1 as published +// by the Free Software Foundation, with special exception defined in the file +// OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT +// distribution for complete text of the license and disclaimer of any warranty. +// +// Alternatively, this file may be used under the terms of Open CASCADE +// commercial license or contractual agreement. + +#ifndef _Geom2dProp_BezierCurve_HeaderFile +#define _Geom2dProp_BezierCurve_HeaderFile + +#include +#include +#include +#include + +//! @brief Local differential properties for a 2D Bezier curve. +//! +//! Uses numeric root-finding for curvature extrema and inflection points. +//! +//! @warning The caller must ensure that the adaptor pointer remains valid +//! for the entire lifetime of this object. This class does not manage +//! the adaptor's lifetime. +class Geom2dProp_BezierCurve +{ +public: + DEFINE_STANDARD_ALLOC + + //! Constructor with adaptor pointer (non-owning). + //! @param theAdaptor the 2D curve adaptor (must wrap a Bezier curve, must not be null) + Geom2dProp_BezierCurve(const Geom2dAdaptor_Curve* theAdaptor) + : myAdaptor(theAdaptor) + { + } + + //! Non-copyable and non-movable. + Geom2dProp_BezierCurve(const Geom2dProp_BezierCurve&) = delete; + Geom2dProp_BezierCurve& operator=(const Geom2dProp_BezierCurve&) = delete; + Geom2dProp_BezierCurve(Geom2dProp_BezierCurve&&) = delete; + Geom2dProp_BezierCurve& operator=(Geom2dProp_BezierCurve&&) = delete; + + //! Returns the adaptor pointer. + const Geom2dAdaptor_Curve* Adaptor() const { return myAdaptor; } + + //! Compute tangent at given parameter. + Standard_EXPORT Geom2dProp::TangentResult Tangent(double theParam, double theTol) const; + + //! Compute curvature at given parameter. + Standard_EXPORT Geom2dProp::CurvatureResult Curvature(double theParam, double theTol) const; + + //! Compute normal at given parameter. + Standard_EXPORT Geom2dProp::NormalResult Normal(double theParam, double theTol) const; + + //! Compute centre of curvature at given parameter. + Standard_EXPORT Geom2dProp::CentreResult CentreOfCurvature(double theParam, double theTol) const; + + //! Find curvature extrema using numeric root-finding. + Standard_EXPORT Geom2dProp::CurveAnalysis FindCurvatureExtrema() const; + + //! Find inflection points using numeric root-finding. + Standard_EXPORT Geom2dProp::CurveAnalysis FindInflections() const; + +private: + const Geom2dAdaptor_Curve* myAdaptor; +}; + +#endif // _Geom2dProp_BezierCurve_HeaderFile diff --git a/src/ModelingData/TKG2d/Geom2dProp/Geom2dProp_Circle.hxx b/src/ModelingData/TKG2d/Geom2dProp/Geom2dProp_Circle.hxx new file mode 100644 index 0000000000..ebfd0e1006 --- /dev/null +++ b/src/ModelingData/TKG2d/Geom2dProp/Geom2dProp_Circle.hxx @@ -0,0 +1,116 @@ +// Copyright (c) 2025 OPEN CASCADE SAS +// +// This file is part of Open CASCADE Technology software library. +// +// This library is free software; you can redistribute it and/or modify it under +// the terms of the GNU Lesser General Public License version 2.1 as published +// by the Free Software Foundation, with special exception defined in the file +// OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT +// distribution for complete text of the license and disclaimer of any warranty. +// +// Alternatively, this file may be used under the terms of Open CASCADE +// commercial license or contractual agreement. + +#ifndef _Geom2dProp_Circle_HeaderFile +#define _Geom2dProp_Circle_HeaderFile + +#include +#include +#include + +//! @brief Local differential properties for a 2D circle. +//! +//! A circle has constant curvature = 1/R, well-defined tangent and normal +//! at every point, and no curvature extrema or inflection points. +//! +//! @warning The caller must ensure that the adaptor pointer remains valid +//! for the entire lifetime of this object. This class does not manage +//! the adaptor's lifetime. +class Geom2dProp_Circle +{ +public: + DEFINE_STANDARD_ALLOC + + //! Constructor with adaptor pointer (non-owning). + //! @param theAdaptor the 2D curve adaptor (must wrap a circle, must not be null) + Geom2dProp_Circle(const Geom2dAdaptor_Curve* theAdaptor) + : myAdaptor(theAdaptor) + { + } + + //! Non-copyable and non-movable. + Geom2dProp_Circle(const Geom2dProp_Circle&) = delete; + Geom2dProp_Circle& operator=(const Geom2dProp_Circle&) = delete; + Geom2dProp_Circle(Geom2dProp_Circle&&) = delete; + Geom2dProp_Circle& operator=(Geom2dProp_Circle&&) = delete; + + //! Returns the adaptor pointer. + const Geom2dAdaptor_Curve* Adaptor() const { return myAdaptor; } + + //! Compute tangent at given parameter. + //! @param[in] theParam curve parameter + //! @param[in] theTol tolerance (unused for circle) + //! @return tangent result (always defined) + Geom2dProp::TangentResult Tangent(double theParam, double theTol) const + { + (void)theTol; + gp_Pnt2d aPnt; + gp_Vec2d aD1; + myAdaptor->D1(theParam, aPnt, aD1); + return {gp_Dir2d(aD1), true}; + } + + //! Compute curvature at given parameter. + //! For a circle, curvature = 1/R (constant). + //! @param[in] theParam curve parameter (unused) + //! @param[in] theTol tolerance (unused) + //! @return curvature result (always defined, constant) + Geom2dProp::CurvatureResult Curvature(double theParam, double theTol) const + { + (void)theParam; + (void)theTol; + return {1.0 / myAdaptor->Circle().Radius(), true, false}; + } + + //! Compute normal at given parameter. + //! @param[in] theParam curve parameter + //! @param[in] theTol tolerance (unused) + //! @return normal result (always defined) + Geom2dProp::NormalResult Normal(double theParam, double theTol) const + { + (void)theTol; + gp_Pnt2d aPnt; + gp_Vec2d aD1, aD2; + myAdaptor->D2(theParam, aPnt, aD1, aD2); + // Normal = D2 * (D1.D1) - D1 * (D1.D2) + const gp_Vec2d aNorm = aD2 * aD1.Dot(aD1) - aD1 * aD1.Dot(aD2); + return {gp_Dir2d(aNorm), true}; + } + + //! Compute centre of curvature at given parameter. + //! For a circle, the centre of curvature is the geometric centre. + //! @param[in] theParam curve parameter (unused) + //! @param[in] theTol tolerance (unused) + //! @return centre result (always the circle centre) + Geom2dProp::CentreResult CentreOfCurvature(double theParam, double theTol) const + { + (void)theParam; + (void)theTol; + return {myAdaptor->Circle().Location(), true}; + } + + //! Find curvature extrema on the circle. + //! A circle has constant curvature, so no extrema. + //! @return empty analysis (always done) + Geom2dProp::CurveAnalysis FindCurvatureExtrema() const { return {{}, true}; } + + //! Find inflection points on the circle. + //! A circle has no inflection points. + //! @return empty analysis (always done) + Geom2dProp::CurveAnalysis FindInflections() const { return {{}, true}; } + +private: + const Geom2dAdaptor_Curve* myAdaptor; +}; + +#endif // _Geom2dProp_Circle_HeaderFile diff --git a/src/ModelingData/TKG2d/Geom2dProp/Geom2dProp_Curve.cxx b/src/ModelingData/TKG2d/Geom2dProp/Geom2dProp_Curve.cxx new file mode 100644 index 0000000000..a71d76f7a8 --- /dev/null +++ b/src/ModelingData/TKG2d/Geom2dProp/Geom2dProp_Curve.cxx @@ -0,0 +1,214 @@ +// Copyright (c) 2025 OPEN CASCADE SAS +// +// This file is part of Open CASCADE Technology software library. +// +// This library is free software; you can redistribute it and/or modify it under +// the terms of the GNU Lesser General Public License version 2.1 as published +// by the Free Software Foundation, with special exception defined in the file +// OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT +// distribution for complete text of the license and disclaimer of any warranty. +// +// Alternatively, this file may be used under the terms of Open CASCADE +// commercial license or contractual agreement. + +#include + +#include +#include + +//================================================================================================== + +void Geom2dProp_Curve::Initialize(const Adaptor2d_Curve2d& theCurve) +{ + if (theCurve.IsKind(STANDARD_TYPE(Geom2dAdaptor_Curve))) + { + const auto& aGeomAdaptor = static_cast(theCurve); + myAdaptor = new Geom2dAdaptor_Curve(aGeomAdaptor); + initFromAdaptor(); + return; + } + + // For non-Geom2dAdaptor, set uninitialized. + myAdaptor.Nullify(); + myCurveType = theCurve.GetType(); + myEvaluator.emplace(); +} + +//================================================================================================== + +void Geom2dProp_Curve::Initialize(const occ::handle& theCurve) +{ + if (theCurve.IsNull()) + { + myAdaptor.Nullify(); + myEvaluator.emplace(); + myCurveType = GeomAbs_OtherCurve; + return; + } + + myAdaptor = new Geom2dAdaptor_Curve(theCurve); + initFromAdaptor(); +} + +//================================================================================================== + +void Geom2dProp_Curve::initFromAdaptor() +{ + myCurveType = myAdaptor->GetType(); + const Geom2dAdaptor_Curve* aPtr = myAdaptor.get(); + + switch (myCurveType) + { + case GeomAbs_Line: + myEvaluator.emplace(aPtr); + break; + case GeomAbs_Circle: + myEvaluator.emplace(aPtr); + break; + case GeomAbs_Ellipse: + myEvaluator.emplace(aPtr); + break; + case GeomAbs_Hyperbola: + myEvaluator.emplace(aPtr); + break; + case GeomAbs_Parabola: + myEvaluator.emplace(aPtr); + break; + case GeomAbs_BezierCurve: + myEvaluator.emplace(aPtr); + break; + case GeomAbs_BSplineCurve: + myEvaluator.emplace(aPtr); + break; + case GeomAbs_OffsetCurve: + myEvaluator.emplace(aPtr); + break; + default: + myEvaluator.emplace(aPtr); + break; + } +} + +//================================================================================================== + +bool Geom2dProp_Curve::IsInitialized() const +{ + return !std::holds_alternative(myEvaluator); +} + +//================================================================================================== + +Geom2dProp::TangentResult Geom2dProp_Curve::Tangent(const double theParam, + const double theTol) const +{ + return std::visit( + [theParam, theTol](const auto& theEval) -> Geom2dProp::TangentResult { + using T = std::decay_t; + if constexpr (std::is_same_v) + { + return {{}, false}; + } + else + { + return theEval.Tangent(theParam, theTol); + } + }, + myEvaluator); +} + +//================================================================================================== + +Geom2dProp::CurvatureResult Geom2dProp_Curve::Curvature(const double theParam, + const double theTol) const +{ + return std::visit( + [theParam, theTol](const auto& theEval) -> Geom2dProp::CurvatureResult { + using T = std::decay_t; + if constexpr (std::is_same_v) + { + return {0.0, false, false}; + } + else + { + return theEval.Curvature(theParam, theTol); + } + }, + myEvaluator); +} + +//================================================================================================== + +Geom2dProp::NormalResult Geom2dProp_Curve::Normal(const double theParam, const double theTol) const +{ + return std::visit( + [theParam, theTol](const auto& theEval) -> Geom2dProp::NormalResult { + using T = std::decay_t; + if constexpr (std::is_same_v) + { + return {{}, false}; + } + else + { + return theEval.Normal(theParam, theTol); + } + }, + myEvaluator); +} + +//================================================================================================== + +Geom2dProp::CentreResult Geom2dProp_Curve::CentreOfCurvature(const double theParam, + const double theTol) const +{ + return std::visit( + [theParam, theTol](const auto& theEval) -> Geom2dProp::CentreResult { + using T = std::decay_t; + if constexpr (std::is_same_v) + { + return {{}, false}; + } + else + { + return theEval.CentreOfCurvature(theParam, theTol); + } + }, + myEvaluator); +} + +//================================================================================================== + +Geom2dProp::CurveAnalysis Geom2dProp_Curve::FindCurvatureExtrema() const +{ + return std::visit( + [](const auto& theEval) -> Geom2dProp::CurveAnalysis { + using T = std::decay_t; + if constexpr (std::is_same_v) + { + return {{}, false}; + } + else + { + return theEval.FindCurvatureExtrema(); + } + }, + myEvaluator); +} + +//================================================================================================== + +Geom2dProp::CurveAnalysis Geom2dProp_Curve::FindInflections() const +{ + return std::visit( + [](const auto& theEval) -> Geom2dProp::CurveAnalysis { + using T = std::decay_t; + if constexpr (std::is_same_v) + { + return {{}, false}; + } + else + { + return theEval.FindInflections(); + } + }, + myEvaluator); +} diff --git a/src/ModelingData/TKG2d/Geom2dProp/Geom2dProp_Curve.hxx b/src/ModelingData/TKG2d/Geom2dProp/Geom2dProp_Curve.hxx new file mode 100644 index 0000000000..e39974390c --- /dev/null +++ b/src/ModelingData/TKG2d/Geom2dProp/Geom2dProp_Curve.hxx @@ -0,0 +1,152 @@ +// Copyright (c) 2025 OPEN CASCADE SAS +// +// This file is part of Open CASCADE Technology software library. +// +// This library is free software; you can redistribute it and/or modify it under +// the terms of the GNU Lesser General Public License version 2.1 as published +// by the Free Software Foundation, with special exception defined in the file +// OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT +// distribution for complete text of the license and disclaimer of any warranty. +// +// Alternatively, this file may be used under the terms of Open CASCADE +// commercial license or contractual agreement. + +#ifndef _Geom2dProp_Curve_HeaderFile +#define _Geom2dProp_Curve_HeaderFile + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +//! @brief Unified local differential property evaluator for any 2D curve. +//! +//! Uses std::variant for compile-time type safety and zero heap allocation +//! for the evaluator itself. Automatically detects curve type from +//! Adaptor2d_Curve2d or Geom2d_Curve and dispatches to the appropriate +//! specialized evaluator. +//! +//! Supported curve types with optimized evaluation: +//! - Line: Trivial (zero curvature, constant tangent) +//! - Circle: Constant curvature 1/R +//! - Ellipse: Analytical curvature extrema at vertices +//! - Hyperbola: Analytical curvature extremum at vertex +//! - Parabola: Analytical curvature extremum at vertex +//! - BezierCurve: Numeric curvature extrema/inflection finding +//! - BSplineCurve: Numeric with C3 interval subdivision +//! - OffsetCurve: Numeric approach +//! - Other: Fallback using Geom2d_Curve virtual D1/D2/D3 +//! +//! Usage: +//! @code +//! Geom2dProp_Curve aProp; +//! aProp.Initialize(myGeom2dCurve); +//! Geom2dProp::CurvatureResult aCurv = aProp.Curvature(0.5, Precision::Confusion()); +//! if (aCurv.IsDefined) +//! { +//! double aValue = aCurv.Value; +//! } +//! @endcode +class Geom2dProp_Curve +{ +public: + DEFINE_STANDARD_ALLOC + + //! Variant type holding all possible 2D curve property evaluators. + using EvaluatorVariant = std::variant; + + //! Default constructor - uninitialized state. + Geom2dProp_Curve() + : myEvaluator(std::monostate{}), + myCurveType(GeomAbs_OtherCurve) + { + } + + //! Non-copyable and non-movable. + Geom2dProp_Curve(const Geom2dProp_Curve&) = delete; + Geom2dProp_Curve& operator=(const Geom2dProp_Curve&) = delete; + Geom2dProp_Curve(Geom2dProp_Curve&&) = delete; + Geom2dProp_Curve& operator=(Geom2dProp_Curve&&) = delete; + + //! Initialize from 2D adaptor reference (auto-detects curve type). + //! For Geom2dAdaptor_Curve, extracts underlying Geom2d_Curve for optimized evaluation. + //! @param[in] theCurve 2D curve adaptor reference + Standard_EXPORT void Initialize(const Adaptor2d_Curve2d& theCurve); + + //! Initialize from geometry handle (auto-detects curve type). + //! @param[in] theCurve 2D geometry to evaluate + Standard_EXPORT void Initialize(const occ::handle& theCurve); + + //! Returns true if properly initialized. + Standard_EXPORT bool IsInitialized() const; + + //! Returns the detected curve type. + GeomAbs_CurveType GetType() const { return myCurveType; } + + //! Compute tangent at given parameter. + //! @param[in] theParam curve parameter + //! @param[in] theTol linear tolerance + //! @return tangent result with validity flag + Standard_EXPORT Geom2dProp::TangentResult Tangent(double theParam, double theTol) const; + + //! Compute curvature at given parameter. + //! @param[in] theParam curve parameter + //! @param[in] theTol linear tolerance + //! @return curvature result with validity and infinity flags + Standard_EXPORT Geom2dProp::CurvatureResult Curvature(double theParam, double theTol) const; + + //! Compute normal at given parameter. + //! @param[in] theParam curve parameter + //! @param[in] theTol linear tolerance + //! @return normal result with validity flag + Standard_EXPORT Geom2dProp::NormalResult Normal(double theParam, double theTol) const; + + //! Compute centre of curvature at given parameter. + //! @param[in] theParam curve parameter + //! @param[in] theTol linear tolerance + //! @return centre result with validity flag + Standard_EXPORT Geom2dProp::CentreResult CentreOfCurvature(double theParam, double theTol) const; + + //! Find curvature extrema on the curve. + //! @return analysis result with special points sorted by parameter + Standard_EXPORT Geom2dProp::CurveAnalysis FindCurvatureExtrema() const; + + //! Find inflection points on the curve. + //! @return analysis result with inflection points sorted by parameter + Standard_EXPORT Geom2dProp::CurveAnalysis FindInflections() const; + +private: + //! Initialize from stored adaptor (dispatches to per-geometry evaluator). + //! Must be called after myAdaptor is set. Per-geometry evaluators receive + //! a non-owning pointer to myAdaptor; their lifetime is managed by the variant. + Standard_EXPORT void initFromAdaptor(); + + occ::handle myAdaptor; //!< Owns the adaptor (ensures lifetime). + EvaluatorVariant myEvaluator; //!< Per-geometry evaluator (non-owning pointer to myAdaptor). + GeomAbs_CurveType myCurveType; +}; + +#endif // _Geom2dProp_Curve_HeaderFile diff --git a/src/ModelingData/TKG2d/Geom2dProp/Geom2dProp_Ellipse.cxx b/src/ModelingData/TKG2d/Geom2dProp/Geom2dProp_Ellipse.cxx new file mode 100644 index 0000000000..26aa24dc98 --- /dev/null +++ b/src/ModelingData/TKG2d/Geom2dProp/Geom2dProp_Ellipse.cxx @@ -0,0 +1,122 @@ +// Copyright (c) 2025 OPEN CASCADE SAS +// +// This file is part of Open CASCADE Technology software library. +// +// This library is free software; you can redistribute it and/or modify it under +// the terms of the GNU Lesser General Public License version 2.1 as published +// by the Free Software Foundation, with special exception defined in the file +// OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT +// distribution for complete text of the license and disclaimer of any warranty. +// +// Alternatively, this file may be used under the terms of Open CASCADE +// commercial license or contractual agreement. + +#include + +#include + +#include + +namespace +{ +constexpr int THE_ELLIPSE_NB_EXTREMA = 4; //!< Number of curvature extrema on full ellipse +constexpr double THE_ELLIPSE_PERIOD = 2.0 * M_PI; //!< One full period of ellipse parameter +} // namespace + +//================================================================================================== + +Geom2dProp::TangentResult Geom2dProp_Ellipse::Tangent(const double theParam, + const double theTol) const +{ + if (myAdaptor == nullptr) + { + return {{}, false}; + } + gp_Pnt2d aPnt; + gp_Vec2d aD1, aD2, aD3; + myAdaptor->D3(theParam, aPnt, aD1, aD2, aD3); + return Geom2dProp::ComputeTangent(aD1, aD2, aD3, theTol); +} + +//================================================================================================== + +Geom2dProp::CurvatureResult Geom2dProp_Ellipse::Curvature(const double theParam, + const double theTol) const +{ + if (myAdaptor == nullptr) + { + return {0.0, false, false}; + } + gp_Pnt2d aPnt; + gp_Vec2d aD1, aD2; + myAdaptor->D2(theParam, aPnt, aD1, aD2); + return Geom2dProp::ComputeCurvature(aD1, aD2, theTol); +} + +//================================================================================================== + +Geom2dProp::NormalResult Geom2dProp_Ellipse::Normal(const double theParam, + const double theTol) const +{ + if (myAdaptor == nullptr) + { + return {{}, false}; + } + gp_Pnt2d aPnt; + gp_Vec2d aD1, aD2; + myAdaptor->D2(theParam, aPnt, aD1, aD2); + return Geom2dProp::ComputeNormal(aD1, aD2, theTol); +} + +//================================================================================================== + +Geom2dProp::CentreResult Geom2dProp_Ellipse::CentreOfCurvature(const double theParam, + const double theTol) const +{ + if (myAdaptor == nullptr) + { + return {{}, false}; + } + gp_Pnt2d aPnt; + gp_Vec2d aD1, aD2; + myAdaptor->D2(theParam, aPnt, aD1, aD2); + return Geom2dProp::ComputeCentreOfCurvature(aPnt, aD1, aD2, theTol); +} + +//================================================================================================== + +Geom2dProp::CurveAnalysis Geom2dProp_Ellipse::FindCurvatureExtrema() const +{ + Geom2dProp::CurveAnalysis aResult; + aResult.IsDone = true; + + if (myAdaptor == nullptr) + { + aResult.IsDone = false; + return aResult; + } + + const double aUFirst = myAdaptor->FirstParameter(); + const double aULast = myAdaptor->LastParameter(); + const double aUFPlus2PI = aUFirst + THE_ELLIPSE_PERIOD; + + // Ellipse curvature extrema at 0, PI/2, PI, 3*PI/2 + // At 0 and PI (major axis endpoints): min radius of curvature -> max |curvature| -> MinCurvature + // At PI/2 and 3*PI/2 (minor axis endpoints): max radius of curvature -> min |curvature| -> + // MaxCurvature + const double aCandidates[] = {0.0, M_PI / 2.0, M_PI, 3.0 * M_PI / 2.0}; + const bool aIsMin[] = {true, false, true, false}; + + for (int i = 0; i < THE_ELLIPSE_NB_EXTREMA; ++i) + { + const double aU = ElCLib::InPeriod(aCandidates[i], aUFirst, aUFPlus2PI); + if (aU >= aUFirst && aU <= aULast) + { + const Geom2dProp::CIType aType = + aIsMin[i] ? Geom2dProp::CIType::MinCurvature : Geom2dProp::CIType::MaxCurvature; + aResult.Points.Append({aU, aType}); + } + } + + return aResult; +} diff --git a/src/ModelingData/TKG2d/Geom2dProp/Geom2dProp_Ellipse.hxx b/src/ModelingData/TKG2d/Geom2dProp/Geom2dProp_Ellipse.hxx new file mode 100644 index 0000000000..1a5071f081 --- /dev/null +++ b/src/ModelingData/TKG2d/Geom2dProp/Geom2dProp_Ellipse.hxx @@ -0,0 +1,76 @@ +// Copyright (c) 2025 OPEN CASCADE SAS +// +// This file is part of Open CASCADE Technology software library. +// +// This library is free software; you can redistribute it and/or modify it under +// the terms of the GNU Lesser General Public License version 2.1 as published +// by the Free Software Foundation, with special exception defined in the file +// OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT +// distribution for complete text of the license and disclaimer of any warranty. +// +// Alternatively, this file may be used under the terms of Open CASCADE +// commercial license or contractual agreement. + +#ifndef _Geom2dProp_Ellipse_HeaderFile +#define _Geom2dProp_Ellipse_HeaderFile + +#include +#include +#include +#include + +//! @brief Local differential properties for a 2D ellipse. +//! +//! An ellipse has analytically known curvature extrema at the four vertices: +//! - Parameter 0 and PI: endpoints of major axis (min radius of curvature) +//! - Parameter PI/2 and 3*PI/2: endpoints of minor axis (max radius of curvature) +//! +//! @warning The caller must ensure that the adaptor pointer remains valid +//! for the entire lifetime of this object. This class does not manage +//! the adaptor's lifetime. +class Geom2dProp_Ellipse +{ +public: + DEFINE_STANDARD_ALLOC + + //! Constructor with adaptor pointer (non-owning). + //! @param theAdaptor the 2D curve adaptor (must wrap an ellipse, must not be null) + Geom2dProp_Ellipse(const Geom2dAdaptor_Curve* theAdaptor) + : myAdaptor(theAdaptor) + { + } + + //! Non-copyable and non-movable. + Geom2dProp_Ellipse(const Geom2dProp_Ellipse&) = delete; + Geom2dProp_Ellipse& operator=(const Geom2dProp_Ellipse&) = delete; + Geom2dProp_Ellipse(Geom2dProp_Ellipse&&) = delete; + Geom2dProp_Ellipse& operator=(Geom2dProp_Ellipse&&) = delete; + + //! Returns the adaptor pointer. + const Geom2dAdaptor_Curve* Adaptor() const { return myAdaptor; } + + //! Compute tangent at given parameter. + Standard_EXPORT Geom2dProp::TangentResult Tangent(double theParam, double theTol) const; + + //! Compute curvature at given parameter. + Standard_EXPORT Geom2dProp::CurvatureResult Curvature(double theParam, double theTol) const; + + //! Compute normal at given parameter. + Standard_EXPORT Geom2dProp::NormalResult Normal(double theParam, double theTol) const; + + //! Compute centre of curvature at given parameter. + Standard_EXPORT Geom2dProp::CentreResult CentreOfCurvature(double theParam, double theTol) const; + + //! Find curvature extrema on the ellipse. + //! Extrema occur analytically at 0, PI/2, PI, and 3*PI/2, filtered to [FirstParam, LastParam]. + Standard_EXPORT Geom2dProp::CurveAnalysis FindCurvatureExtrema() const; + + //! Find inflection points on the ellipse. + //! An ellipse has no inflection points. + Geom2dProp::CurveAnalysis FindInflections() const { return {{}, true}; } + +private: + const Geom2dAdaptor_Curve* myAdaptor; +}; + +#endif // _Geom2dProp_Ellipse_HeaderFile diff --git a/src/ModelingData/TKG2d/Geom2dProp/Geom2dProp_Hyperbola.cxx b/src/ModelingData/TKG2d/Geom2dProp/Geom2dProp_Hyperbola.cxx new file mode 100644 index 0000000000..3796e8a4ab --- /dev/null +++ b/src/ModelingData/TKG2d/Geom2dProp/Geom2dProp_Hyperbola.cxx @@ -0,0 +1,99 @@ +// Copyright (c) 2025 OPEN CASCADE SAS +// +// This file is part of Open CASCADE Technology software library. +// +// This library is free software; you can redistribute it and/or modify it under +// the terms of the GNU Lesser General Public License version 2.1 as published +// by the Free Software Foundation, with special exception defined in the file +// OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT +// distribution for complete text of the license and disclaimer of any warranty. +// +// Alternatively, this file may be used under the terms of Open CASCADE +// commercial license or contractual agreement. + +#include + +//================================================================================================== + +Geom2dProp::TangentResult Geom2dProp_Hyperbola::Tangent(const double theParam, + const double theTol) const +{ + if (myAdaptor == nullptr) + { + return {{}, false}; + } + gp_Pnt2d aPnt; + gp_Vec2d aD1, aD2, aD3; + myAdaptor->D3(theParam, aPnt, aD1, aD2, aD3); + return Geom2dProp::ComputeTangent(aD1, aD2, aD3, theTol); +} + +//================================================================================================== + +Geom2dProp::CurvatureResult Geom2dProp_Hyperbola::Curvature(const double theParam, + const double theTol) const +{ + if (myAdaptor == nullptr) + { + return {0.0, false, false}; + } + gp_Pnt2d aPnt; + gp_Vec2d aD1, aD2; + myAdaptor->D2(theParam, aPnt, aD1, aD2); + return Geom2dProp::ComputeCurvature(aD1, aD2, theTol); +} + +//================================================================================================== + +Geom2dProp::NormalResult Geom2dProp_Hyperbola::Normal(const double theParam, + const double theTol) const +{ + if (myAdaptor == nullptr) + { + return {{}, false}; + } + gp_Pnt2d aPnt; + gp_Vec2d aD1, aD2; + myAdaptor->D2(theParam, aPnt, aD1, aD2); + return Geom2dProp::ComputeNormal(aD1, aD2, theTol); +} + +//================================================================================================== + +Geom2dProp::CentreResult Geom2dProp_Hyperbola::CentreOfCurvature(const double theParam, + const double theTol) const +{ + if (myAdaptor == nullptr) + { + return {{}, false}; + } + gp_Pnt2d aPnt; + gp_Vec2d aD1, aD2; + myAdaptor->D2(theParam, aPnt, aD1, aD2); + return Geom2dProp::ComputeCentreOfCurvature(aPnt, aD1, aD2, theTol); +} + +//================================================================================================== + +Geom2dProp::CurveAnalysis Geom2dProp_Hyperbola::FindCurvatureExtrema() const +{ + Geom2dProp::CurveAnalysis aResult; + aResult.IsDone = true; + + if (myAdaptor == nullptr) + { + aResult.IsDone = false; + return aResult; + } + + const double aUFirst = myAdaptor->FirstParameter(); + const double aULast = myAdaptor->LastParameter(); + + // Hyperbola has maximum |curvature| at parameter 0 (vertex). + if (aUFirst <= 0.0 && aULast >= 0.0) + { + aResult.Points.Append({0.0, Geom2dProp::CIType::MinCurvature}); + } + + return aResult; +} diff --git a/src/ModelingData/TKG2d/Geom2dProp/Geom2dProp_Hyperbola.hxx b/src/ModelingData/TKG2d/Geom2dProp/Geom2dProp_Hyperbola.hxx new file mode 100644 index 0000000000..00d085ba0b --- /dev/null +++ b/src/ModelingData/TKG2d/Geom2dProp/Geom2dProp_Hyperbola.hxx @@ -0,0 +1,75 @@ +// Copyright (c) 2025 OPEN CASCADE SAS +// +// This file is part of Open CASCADE Technology software library. +// +// This library is free software; you can redistribute it and/or modify it under +// the terms of the GNU Lesser General Public License version 2.1 as published +// by the Free Software Foundation, with special exception defined in the file +// OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT +// distribution for complete text of the license and disclaimer of any warranty. +// +// Alternatively, this file may be used under the terms of Open CASCADE +// commercial license or contractual agreement. + +#ifndef _Geom2dProp_Hyperbola_HeaderFile +#define _Geom2dProp_Hyperbola_HeaderFile + +#include +#include +#include +#include + +//! @brief Local differential properties for a 2D hyperbola. +//! +//! A hyperbola has a single curvature extremum (maximum |curvature|) at parameter 0 +//! (the vertex). No inflection points exist. +//! +//! @warning The caller must ensure that the adaptor pointer remains valid +//! for the entire lifetime of this object. This class does not manage +//! the adaptor's lifetime. +class Geom2dProp_Hyperbola +{ +public: + DEFINE_STANDARD_ALLOC + + //! Constructor with adaptor pointer (non-owning). + //! @param theAdaptor the 2D curve adaptor (must wrap a hyperbola, must not be null) + Geom2dProp_Hyperbola(const Geom2dAdaptor_Curve* theAdaptor) + : myAdaptor(theAdaptor) + { + } + + //! Non-copyable and non-movable. + Geom2dProp_Hyperbola(const Geom2dProp_Hyperbola&) = delete; + Geom2dProp_Hyperbola& operator=(const Geom2dProp_Hyperbola&) = delete; + Geom2dProp_Hyperbola(Geom2dProp_Hyperbola&&) = delete; + Geom2dProp_Hyperbola& operator=(Geom2dProp_Hyperbola&&) = delete; + + //! Returns the adaptor pointer. + const Geom2dAdaptor_Curve* Adaptor() const { return myAdaptor; } + + //! Compute tangent at given parameter. + Standard_EXPORT Geom2dProp::TangentResult Tangent(double theParam, double theTol) const; + + //! Compute curvature at given parameter. + Standard_EXPORT Geom2dProp::CurvatureResult Curvature(double theParam, double theTol) const; + + //! Compute normal at given parameter. + Standard_EXPORT Geom2dProp::NormalResult Normal(double theParam, double theTol) const; + + //! Compute centre of curvature at given parameter. + Standard_EXPORT Geom2dProp::CentreResult CentreOfCurvature(double theParam, double theTol) const; + + //! Find curvature extrema on the hyperbola. + //! Single extremum at parameter 0 (the vertex), if within parameter range. + Standard_EXPORT Geom2dProp::CurveAnalysis FindCurvatureExtrema() const; + + //! Find inflection points on the hyperbola. + //! A hyperbola has no inflection points. + Geom2dProp::CurveAnalysis FindInflections() const { return {{}, true}; } + +private: + const Geom2dAdaptor_Curve* myAdaptor; +}; + +#endif // _Geom2dProp_Hyperbola_HeaderFile diff --git a/src/ModelingData/TKG2d/Geom2dProp/Geom2dProp_Line.hxx b/src/ModelingData/TKG2d/Geom2dProp/Geom2dProp_Line.hxx new file mode 100644 index 0000000000..2f9c7103e5 --- /dev/null +++ b/src/ModelingData/TKG2d/Geom2dProp/Geom2dProp_Line.hxx @@ -0,0 +1,112 @@ +// Copyright (c) 2025 OPEN CASCADE SAS +// +// This file is part of Open CASCADE Technology software library. +// +// This library is free software; you can redistribute it and/or modify it under +// the terms of the GNU Lesser General Public License version 2.1 as published +// by the Free Software Foundation, with special exception defined in the file +// OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT +// distribution for complete text of the license and disclaimer of any warranty. +// +// Alternatively, this file may be used under the terms of Open CASCADE +// commercial license or contractual agreement. + +#ifndef _Geom2dProp_Line_HeaderFile +#define _Geom2dProp_Line_HeaderFile + +#include +#include +#include + +//! @brief Local differential properties for a 2D line. +//! +//! A line has constant tangent, zero curvature, undefined normal and centre. +//! No curvature extrema or inflection points exist. +//! +//! @warning The caller must ensure that the adaptor pointer remains valid +//! for the entire lifetime of this object. This class does not manage +//! the adaptor's lifetime. +class Geom2dProp_Line +{ +public: + DEFINE_STANDARD_ALLOC + + //! Constructor with adaptor pointer (non-owning). + //! @param theAdaptor the 2D curve adaptor (must wrap a line, must not be null) + Geom2dProp_Line(const Geom2dAdaptor_Curve* theAdaptor) + : myAdaptor(theAdaptor) + { + } + + //! Non-copyable and non-movable. + Geom2dProp_Line(const Geom2dProp_Line&) = delete; + Geom2dProp_Line& operator=(const Geom2dProp_Line&) = delete; + Geom2dProp_Line(Geom2dProp_Line&&) = delete; + Geom2dProp_Line& operator=(Geom2dProp_Line&&) = delete; + + //! Returns the adaptor pointer. + const Geom2dAdaptor_Curve* Adaptor() const { return myAdaptor; } + + //! Compute tangent at given parameter. + //! For a line, the tangent is always the line direction. + //! @param[in] theParam curve parameter (unused) + //! @param[in] theTol tolerance (unused) + //! @return tangent result (always defined) + Geom2dProp::TangentResult Tangent(double theParam, double theTol) const + { + (void)theParam; + (void)theTol; + return {myAdaptor->Line().Direction(), true}; + } + + //! Compute curvature at given parameter. + //! For a line, curvature is always zero. + //! @param[in] theParam curve parameter (unused) + //! @param[in] theTol tolerance (unused) + //! @return curvature result (always zero) + Geom2dProp::CurvatureResult Curvature(double theParam, double theTol) const + { + (void)theParam; + (void)theTol; + return {0.0, true, false}; + } + + //! Compute normal at given parameter. + //! For a line, the normal is undefined (zero curvature). + //! @param[in] theParam curve parameter (unused) + //! @param[in] theTol tolerance (unused) + //! @return normal result (always undefined) + Geom2dProp::NormalResult Normal(double theParam, double theTol) const + { + (void)theParam; + (void)theTol; + return {{}, false}; + } + + //! Compute centre of curvature at given parameter. + //! For a line, the centre is undefined (zero curvature). + //! @param[in] theParam curve parameter (unused) + //! @param[in] theTol tolerance (unused) + //! @return centre result (always undefined) + Geom2dProp::CentreResult CentreOfCurvature(double theParam, double theTol) const + { + (void)theParam; + (void)theTol; + return {{}, false}; + } + + //! Find curvature extrema on the line. + //! A line has no curvature extrema. + //! @return empty analysis (always done) + Geom2dProp::CurveAnalysis FindCurvatureExtrema() const { return {{}, true}; } + + //! Find inflection points on the line. + //! A line has no inflection points. + //! @return empty analysis (always done) + Geom2dProp::CurveAnalysis FindInflections() const { return {{}, true}; } + +private: + const Geom2dAdaptor_Curve* myAdaptor; +}; + +#endif // _Geom2dProp_Line_HeaderFile diff --git a/src/ModelingData/TKG2d/Geom2dProp/Geom2dProp_OffsetCurve.cxx b/src/ModelingData/TKG2d/Geom2dProp/Geom2dProp_OffsetCurve.cxx new file mode 100644 index 0000000000..4f2f259bbe --- /dev/null +++ b/src/ModelingData/TKG2d/Geom2dProp/Geom2dProp_OffsetCurve.cxx @@ -0,0 +1,322 @@ +// Copyright (c) 2025 OPEN CASCADE SAS +// +// This file is part of Open CASCADE Technology software library. +// +// This library is free software; you can redistribute it and/or modify it under +// the terms of the GNU Lesser General Public License version 2.1 as published +// by the Free Software Foundation, with special exception defined in the file +// OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT +// distribution for complete text of the license and disclaimer of any warranty. +// +// Alternatively, this file may be used under the terms of Open CASCADE +// commercial license or contractual agreement. + +#include + +#include +#include +#include +#include + +#include + +namespace +{ + +constexpr double THE_CURVATURE_DERIV_COEFF = 3.0; //!< Coefficient in d(KC)/dU formula +constexpr double THE_DIFF_STEP_DIVISOR = 100.0; //!< Divisor for numerical differentiation step +constexpr double THE_D2_MAGNITUDE_THRESHOLD = 1.0e-4; //!< Threshold for second derivative magnitude +constexpr double THE_EPSILON_SCALE = 1.0e-4; //!< Scale factor for epsilon relative to domain +constexpr int THE_EXTREMA_NB_SAMPLES = 100; //!< Number of samples for curvature extrema search +constexpr int THE_INFLECTION_NB_SAMPLES = 30; //!< Number of samples for inflection search +constexpr double THE_INFLECTION_TOLERANCE = 1.0e-6; //!< Tolerance for inflection point finding + +//! Function for finding curvature extrema on offset curves. +class FuncCurExt +{ +public: + FuncCurExt(const Geom2dAdaptor_Curve* theCurve, const double theTol) + : myCurve(theCurve), + myEpsX(theTol) + { + } + + bool Value(const double X, double& F) + { + gp_Pnt2d aP; + gp_Vec2d aV1, aV2, aV3; + myCurve->D3(X, aP, aV1, aV2, aV3); + + const double aCPV1V2 = aV1.Crossed(aV2); + const double aCPV1V3 = aV1.Crossed(aV3); + const double aV1V2 = aV1.Dot(aV2); + const double aV1V1 = aV1.SquareMagnitude(); + const double aNV1 = std::sqrt(aV1V1); + const double aV13 = aV1V1 * aNV1; + const double aV15 = aV13 * aV1V1; + + if (aV15 < gp::Resolution()) + { + return false; + } + + F = aCPV1V3 / aV13 - THE_CURVATURE_DERIV_COEFF * aCPV1V2 * aV1V2 / aV15; + return true; + } + + bool Values(const double X, double& F, double& D) + { + double aDx = myEpsX / THE_DIFF_STEP_DIVISOR; + if (X + aDx > myCurve->LastParameter()) + { + aDx = -aDx; + } + + Value(X, F); + double aF2; + Value(X + aDx, aF2); + D = (aF2 - F) / aDx; + return true; + } + + bool IsMinKC(const double X) const + { + gp_Pnt2d aP; + gp_Vec2d aV1, aV2, aV3; + + myCurve->D3(X, aP, aV1, aV2, aV3); + const double aV1V1 = aV1.SquareMagnitude(); + const double aNV1 = std::sqrt(aV1V1); + const double aV13 = aV1V1 * aNV1; + if (aV13 < gp::Resolution()) + { + return false; + } + const double aKC = aV1.Crossed(aV2) / aV13; + + double aDx = myEpsX; + if (X + aDx > myCurve->LastParameter()) + { + aDx = -aDx; + } + + myCurve->D3(X + aDx, aP, aV1, aV2, aV3); + const double aV1V1n = aV1.SquareMagnitude(); + const double aNV1n = std::sqrt(aV1V1n); + const double aV13n = aV1V1n * aNV1n; + if (aV13n < gp::Resolution()) + { + return false; + } + const double aKP = aV1.Crossed(aV2) / aV13n; + + return std::abs(aKC) > std::abs(aKP); + } + +private: + const Geom2dAdaptor_Curve* myCurve; + double myEpsX; +}; + +//! Function for finding inflection points on offset curves. +class FuncCurNul +{ +public: + FuncCurNul(const Geom2dAdaptor_Curve* theCurve) + : myCurve(theCurve) + { + } + + bool Value(const double X, double& F) + { + double aD; + return Values(X, F, aD); + } + + bool Values(const double X, double& F, double& D) + { + gp_Pnt2d aP; + gp_Vec2d aV1, aV2, aV3; + myCurve->D3(X, aP, aV1, aV2, aV3); + + const double aCP1 = aV1.Crossed(aV2); + const double aCP2 = aV1.Crossed(aV3); + const double aV1V2 = aV1.Dot(aV2); + const double aV2V3 = aV2.Dot(aV3); + const double aNV1 = aV1.Magnitude(); + const double aNV2 = aV2.Magnitude(); + + F = 0.0; + D = 0.0; + + if (aNV2 < THE_D2_MAGNITUDE_THRESHOLD) + { + return true; + } + if (aNV1 * aNV2 < gp::Resolution()) + { + return false; + } + + F = aCP1 / (aNV1 * aNV2); + D = (aCP2 - aCP1 * aV1V2 / (aNV1 * aNV1) - aCP1 * aV2V3 / (aNV2 * aNV2)) / (aNV1 * aNV2); + return true; + } + +private: + const Geom2dAdaptor_Curve* myCurve; +}; + +} // namespace + +//================================================================================================== + +Geom2dProp::TangentResult Geom2dProp_OffsetCurve::Tangent(const double theParam, + const double theTol) const +{ + if (myAdaptor == nullptr) + { + return {{}, false}; + } + gp_Pnt2d aPnt; + gp_Vec2d aD1, aD2, aD3; + myAdaptor->D3(theParam, aPnt, aD1, aD2, aD3); + return Geom2dProp::ComputeTangent(aD1, aD2, aD3, theTol); +} + +//================================================================================================== + +Geom2dProp::CurvatureResult Geom2dProp_OffsetCurve::Curvature(const double theParam, + const double theTol) const +{ + if (myAdaptor == nullptr) + { + return {0.0, false, false}; + } + gp_Pnt2d aPnt; + gp_Vec2d aD1, aD2; + myAdaptor->D2(theParam, aPnt, aD1, aD2); + return Geom2dProp::ComputeCurvature(aD1, aD2, theTol); +} + +//================================================================================================== + +Geom2dProp::NormalResult Geom2dProp_OffsetCurve::Normal(const double theParam, + const double theTol) const +{ + if (myAdaptor == nullptr) + { + return {{}, false}; + } + gp_Pnt2d aPnt; + gp_Vec2d aD1, aD2; + myAdaptor->D2(theParam, aPnt, aD1, aD2); + return Geom2dProp::ComputeNormal(aD1, aD2, theTol); +} + +//================================================================================================== + +Geom2dProp::CentreResult Geom2dProp_OffsetCurve::CentreOfCurvature(const double theParam, + const double theTol) const +{ + if (myAdaptor == nullptr) + { + return {{}, false}; + } + gp_Pnt2d aPnt; + gp_Vec2d aD1, aD2; + myAdaptor->D2(theParam, aPnt, aD1, aD2); + return Geom2dProp::ComputeCentreOfCurvature(aPnt, aD1, aD2, theTol); +} + +//================================================================================================== + +Geom2dProp::CurveAnalysis Geom2dProp_OffsetCurve::FindCurvatureExtrema() const +{ + Geom2dProp::CurveAnalysis aResult; + aResult.IsDone = true; + + if (myAdaptor == nullptr) + { + aResult.IsDone = false; + return aResult; + } + + const double aUMin = myAdaptor->FirstParameter(); + const double aUMax = myAdaptor->LastParameter(); + const double aEpsH = THE_EPSILON_SCALE * (aUMax - aUMin); + + FuncCurExt aFunc(myAdaptor, aEpsH); + + MathRoot::MultipleConfig aConfig; + aConfig.NbSamples = THE_EXTREMA_NB_SAMPLES; + aConfig.XTolerance = aEpsH; + aConfig.FTolerance = aEpsH; + + MathRoot::MultipleResult aRoots = + MathRoot::FindAllRootsWithDerivative(aFunc, aUMin, aUMax, aConfig); + + if (aRoots.IsDone()) + { + for (int j = 0; j < aRoots.NbRoots(); ++j) + { + double aParam = aRoots[j]; + MathUtils::Config aBrentCfg; + aBrentCfg.XTolerance = Precision::PConfusion(); + aBrentCfg.FTolerance = Precision::PConfusion(); + auto aBrent = MathRoot::Brent(aFunc, aParam - aEpsH, aParam + aEpsH, aBrentCfg); + if (aBrent.IsDone() && aBrent.Root.has_value()) + { + aParam = *aBrent.Root; + } + const bool aIsMin = aFunc.IsMinKC(aParam); + const Geom2dProp::CIType aType = + aIsMin ? Geom2dProp::CIType::MinCurvature : Geom2dProp::CIType::MaxCurvature; + aResult.Points.Append({aParam, aType}); + } + } + else + { + aResult.IsDone = false; + } + + return aResult; +} + +//================================================================================================== + +Geom2dProp::CurveAnalysis Geom2dProp_OffsetCurve::FindInflections() const +{ + Geom2dProp::CurveAnalysis aResult; + aResult.IsDone = true; + + if (myAdaptor == nullptr) + { + aResult.IsDone = false; + return aResult; + } + + FuncCurNul aFunc(myAdaptor); + + MathRoot::MultipleConfig aConfig; + aConfig.NbSamples = THE_INFLECTION_NB_SAMPLES; + aConfig.XTolerance = THE_INFLECTION_TOLERANCE; + aConfig.FTolerance = THE_INFLECTION_TOLERANCE; + + MathRoot::MultipleResult aRoots = + MathRoot::FindAllRoots(aFunc, myAdaptor->FirstParameter(), myAdaptor->LastParameter(), aConfig); + + if (aRoots.IsDone()) + { + for (int j = 0; j < aRoots.NbRoots(); ++j) + { + aResult.Points.Append({aRoots[j], Geom2dProp::CIType::Inflection}); + } + } + else + { + aResult.IsDone = false; + } + + return aResult; +} diff --git a/src/ModelingData/TKG2d/Geom2dProp/Geom2dProp_OffsetCurve.hxx b/src/ModelingData/TKG2d/Geom2dProp/Geom2dProp_OffsetCurve.hxx new file mode 100644 index 0000000000..db3d9c6f80 --- /dev/null +++ b/src/ModelingData/TKG2d/Geom2dProp/Geom2dProp_OffsetCurve.hxx @@ -0,0 +1,73 @@ +// Copyright (c) 2025 OPEN CASCADE SAS +// +// This file is part of Open CASCADE Technology software library. +// +// This library is free software; you can redistribute it and/or modify it under +// the terms of the GNU Lesser General Public License version 2.1 as published +// by the Free Software Foundation, with special exception defined in the file +// OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT +// distribution for complete text of the license and disclaimer of any warranty. +// +// Alternatively, this file may be used under the terms of Open CASCADE +// commercial license or contractual agreement. + +#ifndef _Geom2dProp_OffsetCurve_HeaderFile +#define _Geom2dProp_OffsetCurve_HeaderFile + +#include +#include +#include +#include + +//! @brief Local differential properties for a 2D offset curve. +//! +//! Uses numeric root-finding for curvature extrema and inflection points. +//! Local properties are computed from the offset curve's own D1/D2/D3. +//! +//! @warning The caller must ensure that the adaptor pointer remains valid +//! for the entire lifetime of this object. This class does not manage +//! the adaptor's lifetime. +class Geom2dProp_OffsetCurve +{ +public: + DEFINE_STANDARD_ALLOC + + //! Constructor with adaptor pointer (non-owning). + //! @param theAdaptor the 2D curve adaptor (must wrap an offset curve, must not be null) + Geom2dProp_OffsetCurve(const Geom2dAdaptor_Curve* theAdaptor) + : myAdaptor(theAdaptor) + { + } + + //! Non-copyable and non-movable. + Geom2dProp_OffsetCurve(const Geom2dProp_OffsetCurve&) = delete; + Geom2dProp_OffsetCurve& operator=(const Geom2dProp_OffsetCurve&) = delete; + Geom2dProp_OffsetCurve(Geom2dProp_OffsetCurve&&) = delete; + Geom2dProp_OffsetCurve& operator=(Geom2dProp_OffsetCurve&&) = delete; + + //! Returns the adaptor pointer. + const Geom2dAdaptor_Curve* Adaptor() const { return myAdaptor; } + + //! Compute tangent at given parameter. + Standard_EXPORT Geom2dProp::TangentResult Tangent(double theParam, double theTol) const; + + //! Compute curvature at given parameter. + Standard_EXPORT Geom2dProp::CurvatureResult Curvature(double theParam, double theTol) const; + + //! Compute normal at given parameter. + Standard_EXPORT Geom2dProp::NormalResult Normal(double theParam, double theTol) const; + + //! Compute centre of curvature at given parameter. + Standard_EXPORT Geom2dProp::CentreResult CentreOfCurvature(double theParam, double theTol) const; + + //! Find curvature extrema using numeric root-finding. + Standard_EXPORT Geom2dProp::CurveAnalysis FindCurvatureExtrema() const; + + //! Find inflection points using numeric root-finding. + Standard_EXPORT Geom2dProp::CurveAnalysis FindInflections() const; + +private: + const Geom2dAdaptor_Curve* myAdaptor; +}; + +#endif // _Geom2dProp_OffsetCurve_HeaderFile diff --git a/src/ModelingData/TKG2d/Geom2dProp/Geom2dProp_OtherCurve.cxx b/src/ModelingData/TKG2d/Geom2dProp/Geom2dProp_OtherCurve.cxx new file mode 100644 index 0000000000..703d51a73d --- /dev/null +++ b/src/ModelingData/TKG2d/Geom2dProp/Geom2dProp_OtherCurve.cxx @@ -0,0 +1,322 @@ +// Copyright (c) 2025 OPEN CASCADE SAS +// +// This file is part of Open CASCADE Technology software library. +// +// This library is free software; you can redistribute it and/or modify it under +// the terms of the GNU Lesser General Public License version 2.1 as published +// by the Free Software Foundation, with special exception defined in the file +// OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT +// distribution for complete text of the license and disclaimer of any warranty. +// +// Alternatively, this file may be used under the terms of Open CASCADE +// commercial license or contractual agreement. + +#include + +#include +#include +#include +#include + +#include + +namespace +{ + +constexpr double THE_CURVATURE_DERIV_COEFF = 3.0; //!< Coefficient in d(KC)/dU formula +constexpr double THE_DIFF_STEP_DIVISOR = 100.0; //!< Divisor for numerical differentiation step +constexpr double THE_D2_MAGNITUDE_THRESHOLD = 1.0e-4; //!< Threshold for second derivative magnitude +constexpr double THE_EPSILON_SCALE = 1.0e-4; //!< Scale factor for epsilon relative to domain +constexpr int THE_EXTREMA_NB_SAMPLES = 100; //!< Number of samples for curvature extrema search +constexpr int THE_INFLECTION_NB_SAMPLES = 30; //!< Number of samples for inflection search +constexpr double THE_INFLECTION_TOLERANCE = 1.0e-6; //!< Tolerance for inflection point finding + +//! Function for finding curvature extrema. +class FuncCurExt +{ +public: + FuncCurExt(const Geom2dAdaptor_Curve* theCurve, const double theTol) + : myCurve(theCurve), + myEpsX(theTol) + { + } + + bool Value(const double X, double& F) + { + gp_Pnt2d aP; + gp_Vec2d aV1, aV2, aV3; + myCurve->D3(X, aP, aV1, aV2, aV3); + + const double aCPV1V2 = aV1.Crossed(aV2); + const double aCPV1V3 = aV1.Crossed(aV3); + const double aV1V2 = aV1.Dot(aV2); + const double aV1V1 = aV1.SquareMagnitude(); + const double aNV1 = std::sqrt(aV1V1); + const double aV13 = aV1V1 * aNV1; + const double aV15 = aV13 * aV1V1; + + if (aV15 < gp::Resolution()) + { + return false; + } + + F = aCPV1V3 / aV13 - THE_CURVATURE_DERIV_COEFF * aCPV1V2 * aV1V2 / aV15; + return true; + } + + bool Values(const double X, double& F, double& D) + { + double aDx = myEpsX / THE_DIFF_STEP_DIVISOR; + if (X + aDx > myCurve->LastParameter()) + { + aDx = -aDx; + } + + Value(X, F); + double aF2; + Value(X + aDx, aF2); + D = (aF2 - F) / aDx; + return true; + } + + bool IsMinKC(const double X) const + { + gp_Pnt2d aP; + gp_Vec2d aV1, aV2, aV3; + + myCurve->D3(X, aP, aV1, aV2, aV3); + const double aV1V1 = aV1.SquareMagnitude(); + const double aNV1 = std::sqrt(aV1V1); + const double aV13 = aV1V1 * aNV1; + if (aV13 < gp::Resolution()) + { + return false; + } + const double aKC = aV1.Crossed(aV2) / aV13; + + double aDx = myEpsX; + if (X + aDx > myCurve->LastParameter()) + { + aDx = -aDx; + } + + myCurve->D3(X + aDx, aP, aV1, aV2, aV3); + const double aV1V1n = aV1.SquareMagnitude(); + const double aNV1n = std::sqrt(aV1V1n); + const double aV13n = aV1V1n * aNV1n; + if (aV13n < gp::Resolution()) + { + return false; + } + const double aKP = aV1.Crossed(aV2) / aV13n; + + return std::abs(aKC) > std::abs(aKP); + } + +private: + const Geom2dAdaptor_Curve* myCurve; + double myEpsX; +}; + +//! Function for finding inflection points. +class FuncCurNul +{ +public: + FuncCurNul(const Geom2dAdaptor_Curve* theCurve) + : myCurve(theCurve) + { + } + + bool Value(const double X, double& F) + { + double aD; + return Values(X, F, aD); + } + + bool Values(const double X, double& F, double& D) + { + gp_Pnt2d aP; + gp_Vec2d aV1, aV2, aV3; + myCurve->D3(X, aP, aV1, aV2, aV3); + + const double aCP1 = aV1.Crossed(aV2); + const double aCP2 = aV1.Crossed(aV3); + const double aV1V2 = aV1.Dot(aV2); + const double aV2V3 = aV2.Dot(aV3); + const double aNV1 = aV1.Magnitude(); + const double aNV2 = aV2.Magnitude(); + + F = 0.0; + D = 0.0; + + if (aNV2 < THE_D2_MAGNITUDE_THRESHOLD) + { + return true; + } + if (aNV1 * aNV2 < gp::Resolution()) + { + return false; + } + + F = aCP1 / (aNV1 * aNV2); + D = (aCP2 - aCP1 * aV1V2 / (aNV1 * aNV1) - aCP1 * aV2V3 / (aNV2 * aNV2)) / (aNV1 * aNV2); + return true; + } + +private: + const Geom2dAdaptor_Curve* myCurve; +}; + +} // namespace + +//================================================================================================== + +Geom2dProp::TangentResult Geom2dProp_OtherCurve::Tangent(const double theParam, + const double theTol) const +{ + if (myAdaptor == nullptr) + { + return {{}, false}; + } + gp_Pnt2d aPnt; + gp_Vec2d aD1, aD2, aD3; + myAdaptor->D3(theParam, aPnt, aD1, aD2, aD3); + return Geom2dProp::ComputeTangent(aD1, aD2, aD3, theTol); +} + +//================================================================================================== + +Geom2dProp::CurvatureResult Geom2dProp_OtherCurve::Curvature(const double theParam, + const double theTol) const +{ + if (myAdaptor == nullptr) + { + return {0.0, false, false}; + } + gp_Pnt2d aPnt; + gp_Vec2d aD1, aD2; + myAdaptor->D2(theParam, aPnt, aD1, aD2); + return Geom2dProp::ComputeCurvature(aD1, aD2, theTol); +} + +//================================================================================================== + +Geom2dProp::NormalResult Geom2dProp_OtherCurve::Normal(const double theParam, + const double theTol) const +{ + if (myAdaptor == nullptr) + { + return {{}, false}; + } + gp_Pnt2d aPnt; + gp_Vec2d aD1, aD2; + myAdaptor->D2(theParam, aPnt, aD1, aD2); + return Geom2dProp::ComputeNormal(aD1, aD2, theTol); +} + +//================================================================================================== + +Geom2dProp::CentreResult Geom2dProp_OtherCurve::CentreOfCurvature(const double theParam, + const double theTol) const +{ + if (myAdaptor == nullptr) + { + return {{}, false}; + } + gp_Pnt2d aPnt; + gp_Vec2d aD1, aD2; + myAdaptor->D2(theParam, aPnt, aD1, aD2); + return Geom2dProp::ComputeCentreOfCurvature(aPnt, aD1, aD2, theTol); +} + +//================================================================================================== + +Geom2dProp::CurveAnalysis Geom2dProp_OtherCurve::FindCurvatureExtrema() const +{ + Geom2dProp::CurveAnalysis aResult; + aResult.IsDone = true; + + if (myAdaptor == nullptr) + { + aResult.IsDone = false; + return aResult; + } + + const double aUMin = myAdaptor->FirstParameter(); + const double aUMax = myAdaptor->LastParameter(); + const double aEpsH = THE_EPSILON_SCALE * (aUMax - aUMin); + + FuncCurExt aFunc(myAdaptor, aEpsH); + + MathRoot::MultipleConfig aConfig; + aConfig.NbSamples = THE_EXTREMA_NB_SAMPLES; + aConfig.XTolerance = aEpsH; + aConfig.FTolerance = aEpsH; + + MathRoot::MultipleResult aRoots = + MathRoot::FindAllRootsWithDerivative(aFunc, aUMin, aUMax, aConfig); + + if (aRoots.IsDone()) + { + for (int j = 0; j < aRoots.NbRoots(); ++j) + { + double aParam = aRoots[j]; + MathUtils::Config aBrentCfg; + aBrentCfg.XTolerance = Precision::PConfusion(); + aBrentCfg.FTolerance = Precision::PConfusion(); + auto aBrent = MathRoot::Brent(aFunc, aParam - aEpsH, aParam + aEpsH, aBrentCfg); + if (aBrent.IsDone() && aBrent.Root.has_value()) + { + aParam = *aBrent.Root; + } + const bool aIsMin = aFunc.IsMinKC(aParam); + const Geom2dProp::CIType aType = + aIsMin ? Geom2dProp::CIType::MinCurvature : Geom2dProp::CIType::MaxCurvature; + aResult.Points.Append({aParam, aType}); + } + } + else + { + aResult.IsDone = false; + } + + return aResult; +} + +//================================================================================================== + +Geom2dProp::CurveAnalysis Geom2dProp_OtherCurve::FindInflections() const +{ + Geom2dProp::CurveAnalysis aResult; + aResult.IsDone = true; + + if (myAdaptor == nullptr) + { + aResult.IsDone = false; + return aResult; + } + + FuncCurNul aFunc(myAdaptor); + + MathRoot::MultipleConfig aConfig; + aConfig.NbSamples = THE_INFLECTION_NB_SAMPLES; + aConfig.XTolerance = THE_INFLECTION_TOLERANCE; + aConfig.FTolerance = THE_INFLECTION_TOLERANCE; + + MathRoot::MultipleResult aRoots = + MathRoot::FindAllRoots(aFunc, myAdaptor->FirstParameter(), myAdaptor->LastParameter(), aConfig); + + if (aRoots.IsDone()) + { + for (int j = 0; j < aRoots.NbRoots(); ++j) + { + aResult.Points.Append({aRoots[j], Geom2dProp::CIType::Inflection}); + } + } + else + { + aResult.IsDone = false; + } + + return aResult; +} diff --git a/src/ModelingData/TKG2d/Geom2dProp/Geom2dProp_OtherCurve.hxx b/src/ModelingData/TKG2d/Geom2dProp/Geom2dProp_OtherCurve.hxx new file mode 100644 index 0000000000..eb35620fae --- /dev/null +++ b/src/ModelingData/TKG2d/Geom2dProp/Geom2dProp_OtherCurve.hxx @@ -0,0 +1,73 @@ +// Copyright (c) 2025 OPEN CASCADE SAS +// +// This file is part of Open CASCADE Technology software library. +// +// This library is free software; you can redistribute it and/or modify it under +// the terms of the GNU Lesser General Public License version 2.1 as published +// by the Free Software Foundation, with special exception defined in the file +// OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT +// distribution for complete text of the license and disclaimer of any warranty. +// +// Alternatively, this file may be used under the terms of Open CASCADE +// commercial license or contractual agreement. + +#ifndef _Geom2dProp_OtherCurve_HeaderFile +#define _Geom2dProp_OtherCurve_HeaderFile + +#include +#include +#include +#include + +//! @brief Fallback local differential properties for any 2D curve type. +//! +//! Uses adaptor D1/D2/D3 methods for property computation +//! and numeric root-finding for curvature extrema and inflection points. +//! +//! @warning The caller must ensure that the adaptor pointer remains valid +//! for the entire lifetime of this object. This class does not manage +//! the adaptor's lifetime. +class Geom2dProp_OtherCurve +{ +public: + DEFINE_STANDARD_ALLOC + + //! Constructor with adaptor pointer (non-owning). + //! @param theAdaptor the 2D curve adaptor (must not be null) + Geom2dProp_OtherCurve(const Geom2dAdaptor_Curve* theAdaptor) + : myAdaptor(theAdaptor) + { + } + + //! Non-copyable and non-movable. + Geom2dProp_OtherCurve(const Geom2dProp_OtherCurve&) = delete; + Geom2dProp_OtherCurve& operator=(const Geom2dProp_OtherCurve&) = delete; + Geom2dProp_OtherCurve(Geom2dProp_OtherCurve&&) = delete; + Geom2dProp_OtherCurve& operator=(Geom2dProp_OtherCurve&&) = delete; + + //! Returns the adaptor pointer. + const Geom2dAdaptor_Curve* Adaptor() const { return myAdaptor; } + + //! Compute tangent at given parameter. + Standard_EXPORT Geom2dProp::TangentResult Tangent(double theParam, double theTol) const; + + //! Compute curvature at given parameter. + Standard_EXPORT Geom2dProp::CurvatureResult Curvature(double theParam, double theTol) const; + + //! Compute normal at given parameter. + Standard_EXPORT Geom2dProp::NormalResult Normal(double theParam, double theTol) const; + + //! Compute centre of curvature at given parameter. + Standard_EXPORT Geom2dProp::CentreResult CentreOfCurvature(double theParam, double theTol) const; + + //! Find curvature extrema using numeric root-finding. + Standard_EXPORT Geom2dProp::CurveAnalysis FindCurvatureExtrema() const; + + //! Find inflection points using numeric root-finding. + Standard_EXPORT Geom2dProp::CurveAnalysis FindInflections() const; + +private: + const Geom2dAdaptor_Curve* myAdaptor; +}; + +#endif // _Geom2dProp_OtherCurve_HeaderFile diff --git a/src/ModelingData/TKG2d/Geom2dProp/Geom2dProp_Parabola.cxx b/src/ModelingData/TKG2d/Geom2dProp/Geom2dProp_Parabola.cxx new file mode 100644 index 0000000000..9d35ea4d4b --- /dev/null +++ b/src/ModelingData/TKG2d/Geom2dProp/Geom2dProp_Parabola.cxx @@ -0,0 +1,99 @@ +// Copyright (c) 2025 OPEN CASCADE SAS +// +// This file is part of Open CASCADE Technology software library. +// +// This library is free software; you can redistribute it and/or modify it under +// the terms of the GNU Lesser General Public License version 2.1 as published +// by the Free Software Foundation, with special exception defined in the file +// OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT +// distribution for complete text of the license and disclaimer of any warranty. +// +// Alternatively, this file may be used under the terms of Open CASCADE +// commercial license or contractual agreement. + +#include + +//================================================================================================== + +Geom2dProp::TangentResult Geom2dProp_Parabola::Tangent(const double theParam, + const double theTol) const +{ + if (myAdaptor == nullptr) + { + return {{}, false}; + } + gp_Pnt2d aPnt; + gp_Vec2d aD1, aD2, aD3; + myAdaptor->D3(theParam, aPnt, aD1, aD2, aD3); + return Geom2dProp::ComputeTangent(aD1, aD2, aD3, theTol); +} + +//================================================================================================== + +Geom2dProp::CurvatureResult Geom2dProp_Parabola::Curvature(const double theParam, + const double theTol) const +{ + if (myAdaptor == nullptr) + { + return {0.0, false, false}; + } + gp_Pnt2d aPnt; + gp_Vec2d aD1, aD2; + myAdaptor->D2(theParam, aPnt, aD1, aD2); + return Geom2dProp::ComputeCurvature(aD1, aD2, theTol); +} + +//================================================================================================== + +Geom2dProp::NormalResult Geom2dProp_Parabola::Normal(const double theParam, + const double theTol) const +{ + if (myAdaptor == nullptr) + { + return {{}, false}; + } + gp_Pnt2d aPnt; + gp_Vec2d aD1, aD2; + myAdaptor->D2(theParam, aPnt, aD1, aD2); + return Geom2dProp::ComputeNormal(aD1, aD2, theTol); +} + +//================================================================================================== + +Geom2dProp::CentreResult Geom2dProp_Parabola::CentreOfCurvature(const double theParam, + const double theTol) const +{ + if (myAdaptor == nullptr) + { + return {{}, false}; + } + gp_Pnt2d aPnt; + gp_Vec2d aD1, aD2; + myAdaptor->D2(theParam, aPnt, aD1, aD2); + return Geom2dProp::ComputeCentreOfCurvature(aPnt, aD1, aD2, theTol); +} + +//================================================================================================== + +Geom2dProp::CurveAnalysis Geom2dProp_Parabola::FindCurvatureExtrema() const +{ + Geom2dProp::CurveAnalysis aResult; + aResult.IsDone = true; + + if (myAdaptor == nullptr) + { + aResult.IsDone = false; + return aResult; + } + + const double aUFirst = myAdaptor->FirstParameter(); + const double aULast = myAdaptor->LastParameter(); + + // Parabola has maximum |curvature| at parameter 0 (vertex). + if (aUFirst <= 0.0 && aULast >= 0.0) + { + aResult.Points.Append({0.0, Geom2dProp::CIType::MinCurvature}); + } + + return aResult; +} diff --git a/src/ModelingData/TKG2d/Geom2dProp/Geom2dProp_Parabola.hxx b/src/ModelingData/TKG2d/Geom2dProp/Geom2dProp_Parabola.hxx new file mode 100644 index 0000000000..8a1378ffd5 --- /dev/null +++ b/src/ModelingData/TKG2d/Geom2dProp/Geom2dProp_Parabola.hxx @@ -0,0 +1,75 @@ +// Copyright (c) 2025 OPEN CASCADE SAS +// +// This file is part of Open CASCADE Technology software library. +// +// This library is free software; you can redistribute it and/or modify it under +// the terms of the GNU Lesser General Public License version 2.1 as published +// by the Free Software Foundation, with special exception defined in the file +// OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT +// distribution for complete text of the license and disclaimer of any warranty. +// +// Alternatively, this file may be used under the terms of Open CASCADE +// commercial license or contractual agreement. + +#ifndef _Geom2dProp_Parabola_HeaderFile +#define _Geom2dProp_Parabola_HeaderFile + +#include +#include +#include +#include + +//! @brief Local differential properties for a 2D parabola. +//! +//! A parabola has a single curvature extremum (maximum |curvature|) at parameter 0 +//! (the vertex). No inflection points exist. +//! +//! @warning The caller must ensure that the adaptor pointer remains valid +//! for the entire lifetime of this object. This class does not manage +//! the adaptor's lifetime. +class Geom2dProp_Parabola +{ +public: + DEFINE_STANDARD_ALLOC + + //! Constructor with adaptor pointer (non-owning). + //! @param theAdaptor the 2D curve adaptor (must wrap a parabola, must not be null) + Geom2dProp_Parabola(const Geom2dAdaptor_Curve* theAdaptor) + : myAdaptor(theAdaptor) + { + } + + //! Non-copyable and non-movable. + Geom2dProp_Parabola(const Geom2dProp_Parabola&) = delete; + Geom2dProp_Parabola& operator=(const Geom2dProp_Parabola&) = delete; + Geom2dProp_Parabola(Geom2dProp_Parabola&&) = delete; + Geom2dProp_Parabola& operator=(Geom2dProp_Parabola&&) = delete; + + //! Returns the adaptor pointer. + const Geom2dAdaptor_Curve* Adaptor() const { return myAdaptor; } + + //! Compute tangent at given parameter. + Standard_EXPORT Geom2dProp::TangentResult Tangent(double theParam, double theTol) const; + + //! Compute curvature at given parameter. + Standard_EXPORT Geom2dProp::CurvatureResult Curvature(double theParam, double theTol) const; + + //! Compute normal at given parameter. + Standard_EXPORT Geom2dProp::NormalResult Normal(double theParam, double theTol) const; + + //! Compute centre of curvature at given parameter. + Standard_EXPORT Geom2dProp::CentreResult CentreOfCurvature(double theParam, double theTol) const; + + //! Find curvature extrema on the parabola. + //! Single extremum at parameter 0 (the vertex), if within parameter range. + Standard_EXPORT Geom2dProp::CurveAnalysis FindCurvatureExtrema() const; + + //! Find inflection points on the parabola. + //! A parabola has no inflection points. + Geom2dProp::CurveAnalysis FindInflections() const { return {{}, true}; } + +private: + const Geom2dAdaptor_Curve* myAdaptor; +}; + +#endif // _Geom2dProp_Parabola_HeaderFile diff --git a/src/ModelingData/TKG2d/PACKAGES.cmake b/src/ModelingData/TKG2d/PACKAGES.cmake index f22570041f..731d38f1c2 100644 --- a/src/ModelingData/TKG2d/PACKAGES.cmake +++ b/src/ModelingData/TKG2d/PACKAGES.cmake @@ -8,4 +8,5 @@ set(OCCT_TKG2d_LIST_OF_PACKAGES Geom2dHash Geom2dGridEval Geom2dEval + Geom2dProp )