lili/OCCT
1
0
Fork 0
mirror of https://github.com/Open-Cascade-SAS/OCCT.git synced 2026-09-09 22:28:13 +08:00
Code Issues Packages Projects Releases Wiki Activity

Modeling Data - Add Geom2dProp package for modern 2D curve differential properties (#1113)

Browse Source 
Replace archaic LProp/Geom2dLProp macro-based (.gxx) pattern with a modern C++17 std::variant-dispatched package for computing local differential properties of 2D curves: tangent, curvature, normal, centre of curvature, curvature extrema, and inflection points.

New package Geom2dProp (TKG2d) provides:
- Geom2dProp: Result structs (TangentResult, CurvatureResult, NormalResult, CentreResult, CurveAnalysis) and geometry-agnostic free functions for property computation from derivatives.
- Geom2dProp_Curve: Unified variant dispatcher that auto-detects curve type from Geom2d_Curve or Adaptor2d_Curve2d and delegates to specialized evaluators. Owns the Geom2dAdaptor_Curve handle and passes non-owning raw pointers to per-geometry classes.
- Per-geometry evaluators with optimized evaluation:
  - Line (header-only): zero curvature, constant tangent
  - Circle (header-only): constant curvature 1/R
  - Ellipse: analytical extrema at 0, PI/2, PI, 3PI/2
  - Hyperbola: analytical extremum at vertex
  - Parabola: analytical extremum at vertex
  - BezierCurve: numeric curvature extrema/inflection finding
  - BSplineCurve: numeric with C3 interval subdivision
  - OffsetCurve: numeric approach
  - OtherCurve: fallback via adaptor virtual D1/D2/D3

Key design decisions:
- Uses Geom2dAdaptor_Curve for optimized derivative evaluation (D0-DN) with BSpline span caching.
- Non-owning raw pointers in per-geometry classes; lifetime managed by the Geom2dProp_Curve dispatcher which holds the adaptor handle.
- Returns result structs with IsDefined flags instead of throwing exceptions for degenerate cases.
- 106 GTests covering all curve types, free functions, adaptor/geometry initialization, trimmed curves, cross-validation against LProp.
This commit is contained in:
Pasukhin Dmitry
2026-02-24 17:10:36 +00:00
committed by GitHub
parent 9ee9dc02a4
commit 3266a82318
29 changed files with 6401 additions and 448 deletions
Download Patch File Download Diff File Expand all files Collapse all files
src/FoundationClasses/TKMath/MathRoot/FILES.cmake
+1
View File
@@ -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
src/FoundationClasses/TKMath/MathRoot/MathRoot_Multiple.hxx
+31 -448
View File
@@ -14,16 +14,7 @@
#ifndef _MathRoot_Multiple_HeaderFile
#define _MathRoot_Multiple_HeaderFile
#include <MathUtils_Types.hxx>
#include <MathUtils_Config.hxx>
#include <MathUtils_Core.hxx>
#include <MathRoot_Brent.hxx>
#include <math_Vector.hxx>
#include <math_IntegerVector.hxx>
#include <NCollection_Vector.hxx>
#include <cmath>
#include <MathRoot_MultipleUtils.hxx>
//! @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<double> Roots; //!< Found roots (sorted)
NCollection_Vector<double> 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<Function> aSampleFn{theFunc, aSamples, theConfig.Offset};
MultipleGetValueFn aGetValue{aSamples};
MultipleBrentValueWrapper<Function> aWrapper{theFunc, theConfig.Offset};
MultipleGetRootValueFn<Function> 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<double> 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<Function> aSampleFn{theFunc, aFValues, aDFValues, theConfig.Offset};
MultipleGetValueFn aGetValue{aFValues};
MultipleBrentDerivWrapper<Function> aWrapper{theFunc, theConfig.Offset};
MultipleGetRootDerivFn<Function> aGetRootValue{theFunc};
MultipleTangentialHandler<Function> 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<double> 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.
src/FoundationClasses/TKMath/MathRoot/MathRoot_MultipleUtils.hxx
+443
View File
@@ -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 <MathUtils_Types.hxx>
#include <MathUtils_Config.hxx>
#include <MathRoot_Brent.hxx>
#include <math_Vector.hxx>
#include <NCollection_Vector.hxx>
#include <cmath>
//! @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<double> Roots; //!< Found roots (sorted)
NCollection_Vector<double> 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 <typename Function>
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 <typename Function>
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 <typename Function>
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 <typename Function>
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 <typename Function>
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 <typename Function>
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 <typename Function>
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 <typename SampleFn,
typename GetValueFn,
typename BrentWrapperT,
typename GetRootValueFn,
typename IntervalExtraFn>
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
src/ModelingData/TKG2d/GTests/FILES.cmake
+3
View File
@@ -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
src/ModelingData/TKG2d/GTests/Geom2dProp_Test.cxx
+1502
View File
File diff suppressed because it is too large Load Diff
src/ModelingData/TKG2d/GTests/Geom2dProp_VsCLProps2d_Test.cxx
+603
View File
@@ -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 <Geom2d_BezierCurve.hxx>
#include <Geom2d_BSplineCurve.hxx>
#include <Geom2d_Circle.hxx>
#include <Geom2d_Ellipse.hxx>
#include <Geom2d_Hyperbola.hxx>
#include <Geom2d_Line.hxx>
#include <Geom2d_OffsetCurve.hxx>
#include <Geom2d_Parabola.hxx>
#include <Geom2d_TrimmedCurve.hxx>
#include <Geom2dLProp_CLProps2d.hxx>
#include <Geom2dProp.hxx>
#include <Geom2dProp_Curve.hxx>
#include <gp_Ax2d.hxx>
#include <gp_Circ2d.hxx>
#include <gp_Dir2d.hxx>
#include <gp_Elips2d.hxx>
#include <gp_Hypr2d.hxx>
#include <gp_Lin2d.hxx>
#include <gp_Parab2d.hxx>
#include <gp_Pnt2d.hxx>
#include <NCollection_Array1.hxx>
#include <Precision.hxx>
#include <gtest/gtest.h>
#include <cmath>
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<Geom2d_Line> 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<Geom2d_Line> 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<Geom2d_Circle> 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<Geom2d_Circle> 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<Geom2d_Circle> 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<Geom2d_Ellipse> 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<Geom2d_Ellipse> 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<Geom2d_Ellipse> 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<Geom2d_Hyperbola> 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<Geom2d_Hyperbola> 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<Geom2d_Parabola> 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<Geom2d_Parabola> 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<Geom2d_Parabola> 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<gp_Pnt2d> 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<Geom2d_BezierCurve> 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<gp_Pnt2d> 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<Geom2d_BezierCurve> 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<gp_Pnt2d> 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<Geom2d_BezierCurve> 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<gp_Pnt2d> 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<double> aKnots(1, 3);
aKnots(1) = 0.0;
aKnots(2) = 0.5;
aKnots(3) = 1.0;
NCollection_Array1<int> aMults(1, 3);
aMults(1) = 3;
aMults(2) = 1;
aMults(3) = 3;
occ::handle<Geom2d_BSplineCurve> 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<gp_Pnt2d> 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<double> aKnots(1, 4);
aKnots(1) = 0.0;
aKnots(2) = 0.33;
aKnots(3) = 0.66;
aKnots(4) = 1.0;
NCollection_Array1<int> aMults(1, 4);
aMults(1) = 4;
aMults(2) = 1;
aMults(3) = 1;
aMults(4) = 4;
occ::handle<Geom2d_BSplineCurve> 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<gp_Pnt2d> 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<double> aKnots(1, 2);
aKnots(1) = 0.0;
aKnots(2) = 1.0;
NCollection_Array1<int> aMults(1, 2);
aMults(1) = 5;
aMults(2) = 5;
occ::handle<Geom2d_BSplineCurve> 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<Geom2d_Circle> aCircle = new Geom2d_Circle(aCirc);
occ::handle<Geom2d_OffsetCurve> 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<Geom2d_Ellipse> anEllipse = new Geom2d_Ellipse(anElips);
occ::handle<Geom2d_OffsetCurve> 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<Geom2d_Ellipse> anEllipse = new Geom2d_Ellipse(anElips);
occ::handle<Geom2d_TrimmedCurve> 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<gp_Pnt2d> 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<Geom2d_BezierCurve> aBezier = new Geom2d_BezierCurve(aPoles);
occ::handle<Geom2d_TrimmedCurve> 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);
}
}
src/ModelingData/TKG2d/GTests/Geom2dProp_VsLProp_Test.cxx
+685
View File
@@ -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 <Geom2d_BezierCurve.hxx>
#include <Geom2d_BSplineCurve.hxx>
#include <Geom2d_Circle.hxx>
#include <Geom2d_Ellipse.hxx>
#include <Geom2d_Hyperbola.hxx>
#include <Geom2d_OffsetCurve.hxx>
#include <Geom2d_Parabola.hxx>
#include <Geom2d_TrimmedCurve.hxx>
#include <Geom2dLProp_CurAndInf2d.hxx>
#include <Geom2dProp.hxx>
#include <Geom2dProp_Curve.hxx>
#include <gp_Ax2d.hxx>
#include <gp_Circ2d.hxx>
#include <gp_Dir2d.hxx>
#include <gp_Elips2d.hxx>
#include <gp_Hypr2d.hxx>
#include <gp_Parab2d.hxx>
#include <gp_Pnt2d.hxx>
#include <LProp_CIType.hxx>
#include <NCollection_Array1.hxx>
#include <gtest/gtest.h>
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<Geom2d_Circle> 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<Geom2d_Circle> 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<Geom2d_Ellipse> 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<Geom2d_Ellipse> 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<Geom2d_Ellipse> 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<Geom2d_Ellipse> 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<Geom2d_Hyperbola> 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<Geom2d_Hyperbola> 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<Geom2d_Parabola> 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<Geom2d_Parabola> 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<Geom2d_Parabola> 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<gp_Pnt2d> 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<Geom2d_BezierCurve> 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<gp_Pnt2d> 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<Geom2d_BezierCurve> 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<gp_Pnt2d> 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<Geom2d_BezierCurve> 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<gp_Pnt2d> 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<Geom2d_BezierCurve> 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<gp_Pnt2d> 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<Geom2d_BezierCurve> 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<gp_Pnt2d> 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<double> aKnots(1, 2);
aKnots(1) = 0.0;
aKnots(2) = 1.0;
NCollection_Array1<int> aMults(1, 2);
aMults(1) = 5;
aMults(2) = 5;
occ::handle<Geom2d_BSplineCurve> 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<gp_Pnt2d> 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<double> aKnots(1, 4);
aKnots(1) = 0.0;
aKnots(2) = 0.33;
aKnots(3) = 0.66;
aKnots(4) = 1.0;
NCollection_Array1<int> aMults(1, 4);
aMults(1) = 4;
aMults(2) = 1;
aMults(3) = 1;
aMults(4) = 4;
occ::handle<Geom2d_BSplineCurve> 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<gp_Pnt2d> 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<double> aKnots(1, 4);
aKnots(1) = 0.0;
aKnots(2) = 0.33;
aKnots(3) = 0.66;
aKnots(4) = 1.0;
NCollection_Array1<int> aMults(1, 4);
aMults(1) = 4;
aMults(2) = 1;
aMults(3) = 1;
aMults(4) = 4;
occ::handle<Geom2d_BSplineCurve> 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<gp_Pnt2d> 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<double> aKnots(1, 4);
aKnots(1) = 0.0;
aKnots(2) = 0.33;
aKnots(3) = 0.66;
aKnots(4) = 1.0;
NCollection_Array1<int> aMults(1, 4);
aMults(1) = 3;
aMults(2) = 1;
aMults(3) = 1;
aMults(4) = 3;
occ::handle<Geom2d_BSplineCurve> 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<Geom2d_Ellipse> anEllipse = new Geom2d_Ellipse(anElips);
occ::handle<Geom2d_TrimmedCurve> 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<gp_Pnt2d> 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<Geom2d_BezierCurve> aBezier = new Geom2d_BezierCurve(aPoles);
occ::handle<Geom2d_TrimmedCurve> 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<Geom2d_Ellipse> anEllipse = new Geom2d_Ellipse(anElips);
occ::handle<Geom2d_OffsetCurve> 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<Geom2d_Ellipse> anEllipse = new Geom2d_Ellipse(anElips);
occ::handle<Geom2d_OffsetCurve> 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<Geom2d_Circle> aCircle = new Geom2d_Circle(aCirc);
occ::handle<Geom2d_OffsetCurve> 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());
}
src/ModelingData/TKG2d/Geom2dProp/FILES.cmake
+25
View File
@@ -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
)
src/ModelingData/TKG2d/Geom2dProp/Geom2dProp.cxx
+124
View File
@@ -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 <Geom2dProp.hxx>
#include <cmath>
//==================================================================================================
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};
}
src/ModelingData/TKG2d/Geom2dProp/Geom2dProp.hxx
+129
View File
@@ -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 <gp_Dir2d.hxx>
#include <gp_Pnt2d.hxx>
#include <gp_Vec2d.hxx>
#include <NCollection_DynamicArray.hxx>
#include <Standard.hxx>
//! @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<CurveSpecialPoint> 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
src/ModelingData/TKG2d/Geom2dProp/Geom2dProp_BSplineCurve.cxx
+429
View File
@@ -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 <Geom2dProp_BSplineCurve.hxx>
#include <GeomAbs_Shape.hxx>
#include <gp.hxx>
#include <MathRoot_Brent.hxx>
#include <MathRoot_Multiple.hxx>
#include <NCollection_Array1.hxx>
#include <Precision.hxx>
#include <cmath>
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<Geom2dProp::CurveSpecialPoint> aFiltered;
aFiltered.Append(theResult.Points[0]);
for (int i = 1; i < aNbPts; ++i)
{
bool aIsDuplicate = false;
for (int j = static_cast<int>(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<double> 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<double> 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;
}
src/ModelingData/TKG2d/Geom2dProp/Geom2dProp_BSplineCurve.hxx
+76
View File
@@ -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 <Geom2dAdaptor_Curve.hxx>
#include <Geom2dProp.hxx>
#include <Standard.hxx>
#include <Standard_DefineAlloc.hxx>
//! @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
src/ModelingData/TKG2d/Geom2dProp/Geom2dProp_BezierCurve.cxx
+347
View File
@@ -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 <Geom2dProp_BezierCurve.hxx>
#include <gp.hxx>
#include <MathRoot_Brent.hxx>
#include <MathRoot_Multiple.hxx>
#include <Precision.hxx>
#include <cmath>
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;
}
src/ModelingData/TKG2d/Geom2dProp/Geom2dProp_BezierCurve.hxx
+72
View File
@@ -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 <Geom2dAdaptor_Curve.hxx>
#include <Geom2dProp.hxx>
#include <Standard.hxx>
#include <Standard_DefineAlloc.hxx>
//! @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
src/ModelingData/TKG2d/Geom2dProp/Geom2dProp_Circle.hxx
+116
View File
@@ -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 <Geom2dAdaptor_Curve.hxx>
#include <Geom2dProp.hxx>
#include <Standard_DefineAlloc.hxx>
//! @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
src/ModelingData/TKG2d/Geom2dProp/Geom2dProp_Curve.cxx
+214
View File
@@ -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 <Geom2dProp_Curve.hxx>
#include <Geom2dAdaptor_Curve.hxx>
#include <Geom2d_TrimmedCurve.hxx>
//==================================================================================================
void Geom2dProp_Curve::Initialize(const Adaptor2d_Curve2d& theCurve)
{
if (theCurve.IsKind(STANDARD_TYPE(Geom2dAdaptor_Curve)))
{
const auto& aGeomAdaptor = static_cast<const Geom2dAdaptor_Curve&>(theCurve);
myAdaptor = new Geom2dAdaptor_Curve(aGeomAdaptor);
initFromAdaptor();
return;
}
// For non-Geom2dAdaptor, set uninitialized.
myAdaptor.Nullify();
myCurveType = theCurve.GetType();
myEvaluator.emplace<std::monostate>();
}
//==================================================================================================
void Geom2dProp_Curve::Initialize(const occ::handle<Geom2d_Curve>& theCurve)
{
if (theCurve.IsNull())
{
myAdaptor.Nullify();
myEvaluator.emplace<std::monostate>();
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<Geom2dProp_Line>(aPtr);
break;
case GeomAbs_Circle:
myEvaluator.emplace<Geom2dProp_Circle>(aPtr);
break;
case GeomAbs_Ellipse:
myEvaluator.emplace<Geom2dProp_Ellipse>(aPtr);
break;
case GeomAbs_Hyperbola:
myEvaluator.emplace<Geom2dProp_Hyperbola>(aPtr);
break;
case GeomAbs_Parabola:
myEvaluator.emplace<Geom2dProp_Parabola>(aPtr);
break;
case GeomAbs_BezierCurve:
myEvaluator.emplace<Geom2dProp_BezierCurve>(aPtr);
break;
case GeomAbs_BSplineCurve:
myEvaluator.emplace<Geom2dProp_BSplineCurve>(aPtr);
break;
case GeomAbs_OffsetCurve:
myEvaluator.emplace<Geom2dProp_OffsetCurve>(aPtr);
break;
default:
myEvaluator.emplace<Geom2dProp_OtherCurve>(aPtr);
break;
}
}
//==================================================================================================
bool Geom2dProp_Curve::IsInitialized() const
{
return !std::holds_alternative<std::monostate>(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<decltype(theEval)>;
if constexpr (std::is_same_v<T, std::monostate>)
{
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<decltype(theEval)>;
if constexpr (std::is_same_v<T, std::monostate>)
{
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<decltype(theEval)>;
if constexpr (std::is_same_v<T, std::monostate>)
{
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<decltype(theEval)>;
if constexpr (std::is_same_v<T, std::monostate>)
{
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<decltype(theEval)>;
if constexpr (std::is_same_v<T, std::monostate>)
{
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<decltype(theEval)>;
if constexpr (std::is_same_v<T, std::monostate>)
{
return {{}, false};
}
else
{
return theEval.FindInflections();
}
},
myEvaluator);
}
src/ModelingData/TKG2d/Geom2dProp/Geom2dProp_Curve.hxx
+152
View File
@@ -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 <Adaptor2d_Curve2d.hxx>
#include <Geom2d_Curve.hxx>
#include <Geom2dAdaptor_Curve.hxx>
#include <GeomAbs_CurveType.hxx>
#include <Geom2dProp.hxx>
#include <Geom2dProp_BezierCurve.hxx>
#include <Geom2dProp_BSplineCurve.hxx>
#include <Geom2dProp_Circle.hxx>
#include <Geom2dProp_Ellipse.hxx>
#include <Geom2dProp_Hyperbola.hxx>
#include <Geom2dProp_Line.hxx>
#include <Geom2dProp_OffsetCurve.hxx>
#include <Geom2dProp_OtherCurve.hxx>
#include <Geom2dProp_Parabola.hxx>
#include <Standard.hxx>
#include <Standard_DefineAlloc.hxx>
#include <variant>
//! @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<std::monostate,
Geom2dProp_Line,
Geom2dProp_Circle,
Geom2dProp_Ellipse,
Geom2dProp_Hyperbola,
Geom2dProp_Parabola,
Geom2dProp_BezierCurve,
Geom2dProp_BSplineCurve,
Geom2dProp_OffsetCurve,
Geom2dProp_OtherCurve>;
//! 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<Geom2d_Curve>& 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<Geom2dAdaptor_Curve> myAdaptor; //!< Owns the adaptor (ensures lifetime).
EvaluatorVariant myEvaluator; //!< Per-geometry evaluator (non-owning pointer to myAdaptor).
GeomAbs_CurveType myCurveType;
};
#endif // _Geom2dProp_Curve_HeaderFile
src/ModelingData/TKG2d/Geom2dProp/Geom2dProp_Ellipse.cxx
+122
View File
@@ -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 <Geom2dProp_Ellipse.hxx>
#include <ElCLib.hxx>
#include <cmath>
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;
}
src/ModelingData/TKG2d/Geom2dProp/Geom2dProp_Ellipse.hxx
+76
View File
@@ -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 <Geom2dAdaptor_Curve.hxx>
#include <Geom2dProp.hxx>
#include <Standard.hxx>
#include <Standard_DefineAlloc.hxx>
//! @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
src/ModelingData/TKG2d/Geom2dProp/Geom2dProp_Hyperbola.cxx
+99
View File
@@ -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_Hyperbola.hxx>
//==================================================================================================
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;
}
src/ModelingData/TKG2d/Geom2dProp/Geom2dProp_Hyperbola.hxx
+75
View File
@@ -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 <Geom2dAdaptor_Curve.hxx>
#include <Geom2dProp.hxx>
#include <Standard.hxx>
#include <Standard_DefineAlloc.hxx>
//! @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
src/ModelingData/TKG2d/Geom2dProp/Geom2dProp_Line.hxx
+112
View File
@@ -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 <Geom2dAdaptor_Curve.hxx>
#include <Geom2dProp.hxx>
#include <Standard_DefineAlloc.hxx>
//! @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
src/ModelingData/TKG2d/Geom2dProp/Geom2dProp_OffsetCurve.cxx
+322
View File
@@ -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 <Geom2dProp_OffsetCurve.hxx>
#include <gp.hxx>
#include <MathRoot_Brent.hxx>
#include <MathRoot_Multiple.hxx>
#include <Precision.hxx>
#include <cmath>
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;
}
src/ModelingData/TKG2d/Geom2dProp/Geom2dProp_OffsetCurve.hxx
+73
View File
@@ -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 <Geom2dAdaptor_Curve.hxx>
#include <Geom2dProp.hxx>
#include <Standard.hxx>
#include <Standard_DefineAlloc.hxx>
//! @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
src/ModelingData/TKG2d/Geom2dProp/Geom2dProp_OtherCurve.cxx
+322
View File
@@ -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 <Geom2dProp_OtherCurve.hxx>
#include <gp.hxx>
#include <MathRoot_Brent.hxx>
#include <MathRoot_Multiple.hxx>
#include <Precision.hxx>
#include <cmath>
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;
}
src/ModelingData/TKG2d/Geom2dProp/Geom2dProp_OtherCurve.hxx
+73
View File
@@ -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 <Geom2dAdaptor_Curve.hxx>
#include <Geom2dProp.hxx>
#include <Standard.hxx>
#include <Standard_DefineAlloc.hxx>
//! @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
src/ModelingData/TKG2d/Geom2dProp/Geom2dProp_Parabola.cxx
+99
View File
@@ -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_Parabola.hxx>
//==================================================================================================
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;
}
src/ModelingData/TKG2d/Geom2dProp/Geom2dProp_Parabola.hxx
+75
View File
@@ -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 <Geom2dAdaptor_Curve.hxx>
#include <Geom2dProp.hxx>
#include <Standard.hxx>
#include <Standard_DefineAlloc.hxx>
//! @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
src/ModelingData/TKG2d/PACKAGES.cmake
+1
View File
@@ -8,4 +8,5 @@ set(OCCT_TKG2d_LIST_OF_PACKAGES
Geom2dHash
Geom2dGridEval
Geom2dEval
Geom2dProp
)