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

Foundation Classes - Enhance MathRoot and MathSys with new utilities (#954)

Browse Source 
- Added comprehensive testing infrastructure for trigonometric root finding with multiple edge cases
- Introduced optimized 2D Newton solver with gradient descent fallback and bounds checking
- Created domain management utilities for 1D/2D parameter spaces
- Fixed critical bugs in trigonometric root finding algorithm
This commit is contained in:
Pasukhin Dmitry
2025-12-25 21:21:59 +00:00
committed by GitHub
parent fee21ceac4
commit e7d1148e8f
8 changed files with 1243 additions and 25 deletions
Download Patch File Download Diff File Expand all files Collapse all files

1
src/FoundationClasses/TKMath/GTests/FILES.cmake
View File

@@ -79,6 +79,7 @@ set(OCCT_TKMath_GTests_FILES
MathRoot_Test.cxx
MathRoot_Multiple_Test.cxx
MathRoot_Comparison_Test.cxx
MathRoot_Trig_Test.cxx
# MathInteg tests
MathInteg_Test.cxx
MathInteg_New_Test.cxx

571
src/FoundationClasses/TKMath/GTests/MathRoot_Trig_Test.cxx Normal file
View File

@@ -0,0 +1,571 @@
// 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 <gtest/gtest.h>
#include <MathRoot_Trig.hxx>
#include <cmath>
namespace
{
constexpr double THE_TOL = 1.0e-10;
constexpr double THE_PI = 3.14159265358979323846;
constexpr double THE_2PI = 2.0 * THE_PI;
//! Helper to verify that a root satisfies the equation
//! a*cos^2(x) + 2*b*cos(x)*sin(x) + c*cos(x) + d*sin(x) + e = 0
double evaluateEquation(double theA,
double theB,
double theC,
double theD,
double theE,
double theX)
{
double aCos = std::cos(theX);
double aSin = std::sin(theX);
return theA * aCos * aCos + 2.0 * theB * aCos * aSin + theC * aCos + theD * aSin + theE;
}
//! Helper to verify all roots satisfy the equation
void verifyRoots(const MathRoot::TrigResult& theResult,
double theA,
double theB,
double theC,
double theD,
double theE,
double theTol = 1.0e-10)
{
for (int i = 0; i < theResult.NbRoots; ++i)
{
double aVal = evaluateEquation(theA, theB, theC, theD, theE, theResult.Roots[i]);
EXPECT_NEAR(aVal, 0.0, theTol)
<< "Root " << i << " = " << theResult.Roots[i] << " gives f(x) = " << aVal;
}
}
} // namespace
// ============================================================================
// Basic linear equations: d*sin(x) + e = 0
// ============================================================================
TEST(MathRoot_TrigTest, LinearSin_ZeroConstant)
{
// sin(x) = 0 => x = 0, PI
MathRoot::TrigResult aResult = MathRoot::TrigonometricLinear(1.0, 0.0);
ASSERT_TRUE(aResult.IsDone());
EXPECT_EQ(aResult.NbRoots, 2);
verifyRoots(aResult, 0.0, 0.0, 0.0, 1.0, 0.0);
}
TEST(MathRoot_TrigTest, LinearSin_HalfValue)
{
// sin(x) = 0.5 => x = PI/6, 5*PI/6
MathRoot::TrigResult aResult = MathRoot::TrigonometricLinear(1.0, -0.5);
ASSERT_TRUE(aResult.IsDone());
EXPECT_EQ(aResult.NbRoots, 2);
verifyRoots(aResult, 0.0, 0.0, 0.0, 1.0, -0.5);
}
TEST(MathRoot_TrigTest, LinearSin_NegativeHalf)
{
// sin(x) = -0.5 => x = -PI/6 + 2PI, PI + PI/6
MathRoot::TrigResult aResult = MathRoot::TrigonometricLinear(1.0, 0.5);
ASSERT_TRUE(aResult.IsDone());
EXPECT_EQ(aResult.NbRoots, 2);
verifyRoots(aResult, 0.0, 0.0, 0.0, 1.0, 0.5);
}
TEST(MathRoot_TrigTest, LinearSin_NoSolution)
{
// sin(x) = 2 => no solution
MathRoot::TrigResult aResult = MathRoot::TrigonometricLinear(1.0, -2.0);
ASSERT_TRUE(aResult.IsDone());
EXPECT_EQ(aResult.NbRoots, 0);
}
TEST(MathRoot_TrigTest, LinearSin_BoundaryValue)
{
// sin(x) = 1 => x = PI/2
MathRoot::TrigResult aResult = MathRoot::TrigonometricLinear(1.0, -1.0);
ASSERT_TRUE(aResult.IsDone());
EXPECT_GE(aResult.NbRoots, 1);
verifyRoots(aResult, 0.0, 0.0, 0.0, 1.0, -1.0);
}
// ============================================================================
// Basic equations: c*cos(x) + e = 0
// ============================================================================
TEST(MathRoot_TrigTest, LinearCos_ZeroConstant)
{
// cos(x) = 0 => x = PI/2, 3*PI/2
MathRoot::TrigResult aResult = MathRoot::Trigonometric(0.0, 0.0, 1.0, 0.0, 0.0);
ASSERT_TRUE(aResult.IsDone());
EXPECT_EQ(aResult.NbRoots, 2);
verifyRoots(aResult, 0.0, 0.0, 1.0, 0.0, 0.0);
}
TEST(MathRoot_TrigTest, LinearCos_HalfValue)
{
// cos(x) = 0.5 => x = PI/3, 5*PI/3
MathRoot::TrigResult aResult = MathRoot::Trigonometric(0.0, 0.0, 1.0, 0.0, -0.5);
ASSERT_TRUE(aResult.IsDone());
EXPECT_EQ(aResult.NbRoots, 2);
verifyRoots(aResult, 0.0, 0.0, 1.0, 0.0, -0.5);
}
TEST(MathRoot_TrigTest, LinearCos_NoSolution)
{
// cos(x) = 2 => no solution
MathRoot::TrigResult aResult = MathRoot::Trigonometric(0.0, 0.0, 1.0, 0.0, -2.0);
ASSERT_TRUE(aResult.IsDone());
EXPECT_EQ(aResult.NbRoots, 0);
}
// ============================================================================
// Linear combination: c*cos(x) + d*sin(x) + e = 0
// ============================================================================
TEST(MathRoot_TrigTest, CDE_CosPlusSin)
{
// cos(x) + sin(x) = 0 => tan(x) = -1 => x = 3*PI/4, 7*PI/4
MathRoot::TrigResult aResult = MathRoot::TrigonometricCDE(1.0, 1.0, 0.0);
ASSERT_TRUE(aResult.IsDone());
EXPECT_EQ(aResult.NbRoots, 2);
verifyRoots(aResult, 0.0, 0.0, 1.0, 1.0, 0.0);
}
TEST(MathRoot_TrigTest, CDE_CosPlusSinEqualsSqrt2)
{
// cos(x) + sin(x) = sqrt(2) => x = PI/4 (double root)
// This is a challenging case where F(x) and F'(x) are both zero at the root
MathRoot::TrigResult aResult = MathRoot::TrigonometricCDE(1.0, 1.0, -std::sqrt(2.0));
ASSERT_TRUE(aResult.IsDone());
EXPECT_GE(aResult.NbRoots, 1);
// With Halley's method fallback, precision should be much better
verifyRoots(aResult, 0.0, 0.0, 1.0, 1.0, -std::sqrt(2.0), 1.0e-9);
}
TEST(MathRoot_TrigTest, CDE_CosMinusSin)
{
// cos(x) - sin(x) = 0 => tan(x) = 1 => x = PI/4, 5*PI/4
MathRoot::TrigResult aResult = MathRoot::TrigonometricCDE(1.0, -1.0, 0.0);
ASSERT_TRUE(aResult.IsDone());
EXPECT_EQ(aResult.NbRoots, 2);
verifyRoots(aResult, 0.0, 0.0, 1.0, -1.0, 0.0);
}
TEST(MathRoot_TrigTest, CDE_Harmonic)
{
// 3*cos(x) + 4*sin(x) = 5 => one solution at x = atan(4/3) ~= 0.927
// This is R*cos(x - phi) = 5 where R = 5
MathRoot::TrigResult aResult = MathRoot::TrigonometricCDE(3.0, 4.0, -5.0);
ASSERT_TRUE(aResult.IsDone());
EXPECT_GE(aResult.NbRoots, 1);
verifyRoots(aResult, 0.0, 0.0, 3.0, 4.0, -5.0);
}
TEST(MathRoot_TrigTest, CDE_NoSolution)
{
// cos(x) + sin(x) = 2 => no solution (max value is sqrt(2))
MathRoot::TrigResult aResult = MathRoot::TrigonometricCDE(1.0, 1.0, -2.0);
ASSERT_TRUE(aResult.IsDone());
EXPECT_EQ(aResult.NbRoots, 0);
}
// ============================================================================
// Quadratic in cos: a*cos^2(x) + c*cos(x) + e = 0
// ============================================================================
TEST(MathRoot_TrigTest, QuadraticCos_TwoSolutions)
{
// cos^2(x) - 1 = 0 => cos(x) = +/-1 => x = 0, PI
MathRoot::TrigResult aResult = MathRoot::Trigonometric(1.0, 0.0, 0.0, 0.0, -1.0);
ASSERT_TRUE(aResult.IsDone());
EXPECT_GE(aResult.NbRoots, 2);
verifyRoots(aResult, 1.0, 0.0, 0.0, 0.0, -1.0);
}
TEST(MathRoot_TrigTest, QuadraticCos_FourSolutions)
{
// cos^2(x) = 0.25 => cos(x) = +/-0.5 => 4 solutions
MathRoot::TrigResult aResult = MathRoot::Trigonometric(1.0, 0.0, 0.0, 0.0, -0.25);
ASSERT_TRUE(aResult.IsDone());
EXPECT_EQ(aResult.NbRoots, 4);
verifyRoots(aResult, 1.0, 0.0, 0.0, 0.0, -0.25);
}
// ============================================================================
// Mixed terms: a*cos^2(x) + 2*b*cos(x)*sin(x) = 0
// ============================================================================
TEST(MathRoot_TrigTest, CosSinProduct_AE_Zero)
{
// 2*cos(x)*sin(x) = 0 => sin(2x) = 0 => x = 0, PI/2, PI, 3*PI/2
MathRoot::TrigResult aResult = MathRoot::Trigonometric(0.0, 1.0, 0.0, 0.0, 0.0);
ASSERT_TRUE(aResult.IsDone());
EXPECT_EQ(aResult.NbRoots, 4);
verifyRoots(aResult, 0.0, 1.0, 0.0, 0.0, 0.0);
}
TEST(MathRoot_TrigTest, Mixed_SinTimesCos)
{
// cos(x)*sin(x) + sin(x) = sin(x)*(cos(x) + 1) = 0
// Solutions: sin(x) = 0 => x = 0, PI
// or: cos(x) = -1 => x = PI
MathRoot::TrigResult aResult = MathRoot::Trigonometric(0.0, 0.5, 0.0, 1.0, 0.0);
ASSERT_TRUE(aResult.IsDone());
EXPECT_GE(aResult.NbRoots, 2);
verifyRoots(aResult, 0.0, 0.5, 0.0, 1.0, 0.0);
}
// ============================================================================
// General quartic cases
// ============================================================================
TEST(MathRoot_TrigTest, General_FourRoots)
{
// A general equation with 4 roots
// cos^2(x) + cos(x) - 0.5 = 0
// cos(x) = (-1 +/- sqrt(3))/2
MathRoot::TrigResult aResult = MathRoot::Trigonometric(1.0, 0.0, 1.0, 0.0, -0.5);
ASSERT_TRUE(aResult.IsDone());
// This should have 2-4 roots
EXPECT_GE(aResult.NbRoots, 2);
verifyRoots(aResult, 1.0, 0.0, 1.0, 0.0, -0.5);
}
TEST(MathRoot_TrigTest, General_AllCoefficients)
{
// Test with all coefficients non-zero
MathRoot::TrigResult aResult = MathRoot::Trigonometric(1.0, 0.5, 0.3, 0.4, -0.2);
ASSERT_TRUE(aResult.IsDone());
// Just verify that any roots found satisfy the equation
verifyRoots(aResult, 1.0, 0.5, 0.3, 0.4, -0.2, 1.0e-8);
}
// ============================================================================
// Bound constraint tests
// ============================================================================
TEST(MathRoot_TrigTest, Bounds_FirstQuadrant)
{
// sin(x) = 0.5 in [0, PI/2]
MathRoot::TrigResult aResult = MathRoot::TrigonometricLinear(1.0, -0.5, 0.0, THE_PI / 2.0);
ASSERT_TRUE(aResult.IsDone());
EXPECT_EQ(aResult.NbRoots, 1);
EXPECT_GE(aResult.Roots[0], 0.0 - THE_TOL);
EXPECT_LE(aResult.Roots[0], THE_PI / 2.0 + THE_TOL);
verifyRoots(aResult, 0.0, 0.0, 0.0, 1.0, -0.5);
}
TEST(MathRoot_TrigTest, Bounds_SecondQuadrant)
{
// sin(x) = 0.5 in [PI/2, PI]
MathRoot::TrigResult aResult = MathRoot::TrigonometricLinear(1.0, -0.5, THE_PI / 2.0, THE_PI);
ASSERT_TRUE(aResult.IsDone());
EXPECT_EQ(aResult.NbRoots, 1);
EXPECT_GE(aResult.Roots[0], THE_PI / 2.0 - THE_TOL);
EXPECT_LE(aResult.Roots[0], THE_PI + THE_TOL);
verifyRoots(aResult, 0.0, 0.0, 0.0, 1.0, -0.5);
}
TEST(MathRoot_TrigTest, Bounds_NarrowRange)
{
// cos(x) = 0 in narrow range around PI/2
MathRoot::TrigResult aResult =
MathRoot::Trigonometric(0.0, 0.0, 1.0, 0.0, 0.0, THE_PI / 2.0 - 0.1, THE_PI / 2.0 + 0.1);
ASSERT_TRUE(aResult.IsDone());
EXPECT_EQ(aResult.NbRoots, 1);
EXPECT_NEAR(aResult.Roots[0], THE_PI / 2.0, THE_TOL);
}
TEST(MathRoot_TrigTest, Bounds_RootOutsideRange)
{
// sin(x) = 0.5, roots at PI/6 and 5*PI/6
// Looking in [PI, 2*PI] where these roots don't exist
MathRoot::TrigResult aResult = MathRoot::TrigonometricLinear(1.0, -0.5, THE_PI, THE_2PI);
ASSERT_TRUE(aResult.IsDone());
EXPECT_EQ(aResult.NbRoots, 0);
}
TEST(MathRoot_TrigTest, Bounds_FullCircle)
{
// sin(x) = 0 in [0, 2*PI]
MathRoot::TrigResult aResult = MathRoot::TrigonometricLinear(1.0, 0.0, 0.0, THE_2PI);
ASSERT_TRUE(aResult.IsDone());
EXPECT_EQ(aResult.NbRoots, 2);
}
TEST(MathRoot_TrigTest, Bounds_MultipleCircles)
{
// sin(x) = 0 in [0, 4*PI] - but we only look in one period
MathRoot::TrigResult aResult = MathRoot::TrigonometricLinear(1.0, 0.0, 0.0, 4.0 * THE_PI);
ASSERT_TRUE(aResult.IsDone());
// Should still only return 2 roots in one period
EXPECT_EQ(aResult.NbRoots, 2);
}
// ============================================================================
// Degenerate cases
// ============================================================================
TEST(MathRoot_TrigTest, Degenerate_AllZero)
{
// 0 = 0 => infinite solutions
MathRoot::TrigResult aResult = MathRoot::Trigonometric(0.0, 0.0, 0.0, 0.0, 0.0);
EXPECT_TRUE(aResult.InfiniteRoots);
}
TEST(MathRoot_TrigTest, Degenerate_ConstantNonZero)
{
// e = 0 where e != 0 => no solutions
MathRoot::TrigResult aResult = MathRoot::Trigonometric(0.0, 0.0, 0.0, 0.0, 1.0);
ASSERT_TRUE(aResult.IsDone());
EXPECT_EQ(aResult.NbRoots, 0);
}
TEST(MathRoot_TrigTest, Degenerate_NearZeroCoefficients)
{
// Very small coefficients
MathRoot::TrigResult aResult = MathRoot::Trigonometric(1.0e-14, 1.0e-14, 1.0, 0.0, -0.5);
ASSERT_TRUE(aResult.IsDone());
// Should behave like cos(x) = 0.5
EXPECT_EQ(aResult.NbRoots, 2);
verifyRoots(aResult, 1.0e-14, 1.0e-14, 1.0, 0.0, -0.5, 1.0e-6);
}
// ============================================================================
// Special case: PI as root
// ============================================================================
TEST(MathRoot_TrigTest, SpecialCase_PiRoot)
{
// cos^2(x) - cos(x) = 0 => cos(x)*(cos(x) - 1) = 0
// cos(x) = 0 => x = PI/2, 3*PI/2
// cos(x) = 1 => x = 0
// Note: A - C + E = 1 - 0 + 0 = 1 != 0, so PI is not a special root here
MathRoot::TrigResult aResult = MathRoot::Trigonometric(1.0, 0.0, -1.0, 0.0, 0.0);
ASSERT_TRUE(aResult.IsDone());
EXPECT_GE(aResult.NbRoots, 2);
verifyRoots(aResult, 1.0, 0.0, -1.0, 0.0, 0.0);
}
TEST(MathRoot_TrigTest, SpecialCase_PiRootCheck)
{
// A - C + E = 0 triggers special PI check
// cos^2(x) + cos(x) = 0 => cos(x)*(cos(x) + 1) = 0
// cos(x) = 0 => x = PI/2, 3*PI/2
// cos(x) = -1 => x = PI
// Here A=1, C=-1, E=0 => A - C + E = 1 - (-1) + 0 = 2 != 0... not the right case
// Need A - C + E = 0: try A=1, C=1, E=0
// cos^2(x) - cos(x) = 0: that's A=1, C=-1, E=0 => 1-(-1)+0=2
// Let's try: A=1, C=0, E=-1 => 1-0+(-1)=0
// cos^2(x) - 1 = 0 => cos(x) = +/-1 => x = 0, PI
MathRoot::TrigResult aResult = MathRoot::Trigonometric(1.0, 0.0, 0.0, 0.0, -1.0);
ASSERT_TRUE(aResult.IsDone());
// Check if PI is in the roots
bool aFoundPi = false;
for (int i = 0; i < aResult.NbRoots; ++i)
{
if (std::abs(aResult.Roots[i] - THE_PI) < THE_TOL)
{
aFoundPi = true;
break;
}
}
EXPECT_TRUE(aFoundPi) << "PI should be a root of cos^2(x) - 1 = 0";
verifyRoots(aResult, 1.0, 0.0, 0.0, 0.0, -1.0);
}
// ============================================================================
// Numerical stability tests
// ============================================================================
TEST(MathRoot_TrigTest, Stability_LargeCoefficients)
{
// Equation with large coefficients
MathRoot::TrigResult aResult = MathRoot::Trigonometric(1000.0, 500.0, 300.0, 400.0, -200.0);
ASSERT_TRUE(aResult.IsDone());
verifyRoots(aResult, 1000.0, 500.0, 300.0, 400.0, -200.0, 1.0e-6);
}
TEST(MathRoot_TrigTest, Stability_SmallCoefficients)
{
// Equation with small (but not negligible) coefficients
MathRoot::TrigResult aResult = MathRoot::Trigonometric(1.0e-6, 0.5e-6, 0.3e-6, 0.4e-6, -0.2e-6);
ASSERT_TRUE(aResult.IsDone());
// Verify roots if any found
if (aResult.NbRoots > 0)
{
verifyRoots(aResult, 1.0e-6, 0.5e-6, 0.3e-6, 0.4e-6, -0.2e-6, 1.0e-4);
}
}
TEST(MathRoot_TrigTest, Stability_MixedMagnitude)
{
// Coefficients of very different magnitudes
MathRoot::TrigResult aResult = MathRoot::Trigonometric(1.0, 0.0, 1.0e6, 0.0, -0.5e6);
ASSERT_TRUE(aResult.IsDone());
verifyRoots(aResult, 1.0, 0.0, 1.0e6, 0.0, -0.5e6, 1.0e-4);
}
// ============================================================================
// Root ordering and uniqueness tests
// ============================================================================
TEST(MathRoot_TrigTest, RootOrdering_Sorted)
{
// sin(x) = 0 => x = 0, PI (should be sorted)
MathRoot::TrigResult aResult = MathRoot::TrigonometricLinear(1.0, 0.0);
ASSERT_TRUE(aResult.IsDone());
for (int i = 1; i < aResult.NbRoots; ++i)
{
EXPECT_LE(aResult.Roots[i - 1], aResult.Roots[i]) << "Roots should be in ascending order";
}
}
TEST(MathRoot_TrigTest, RootUniqueness_NoDuplicates)
{
// cos^2(x) = 0.25 => 4 distinct roots
MathRoot::TrigResult aResult = MathRoot::Trigonometric(1.0, 0.0, 0.0, 0.0, -0.25);
ASSERT_TRUE(aResult.IsDone());
for (int i = 0; i < aResult.NbRoots; ++i)
{
for (int j = i + 1; j < aResult.NbRoots; ++j)
{
EXPECT_GT(std::abs(aResult.Roots[i] - aResult.Roots[j]), 1.0e-10)
<< "Roots " << i << " and " << j << " are duplicates";
}
}
}
// ============================================================================
// Wrapper function tests
// ============================================================================
TEST(MathRoot_TrigTest, Wrapper_TrigonometricLinear)
{
// Test TrigonometricLinear wrapper
MathRoot::TrigResult aResult1 = MathRoot::TrigonometricLinear(2.0, -1.0);
MathRoot::TrigResult aResult2 = MathRoot::Trigonometric(0.0, 0.0, 0.0, 2.0, -1.0);
EXPECT_EQ(aResult1.NbRoots, aResult2.NbRoots);
for (int i = 0; i < aResult1.NbRoots; ++i)
{
EXPECT_NEAR(aResult1.Roots[i], aResult2.Roots[i], THE_TOL);
}
}
TEST(MathRoot_TrigTest, Wrapper_TrigonometricCDE)
{
// Test TrigonometricCDE wrapper
MathRoot::TrigResult aResult1 = MathRoot::TrigonometricCDE(1.0, 2.0, -0.5);
MathRoot::TrigResult aResult2 = MathRoot::Trigonometric(0.0, 0.0, 1.0, 2.0, -0.5);
EXPECT_EQ(aResult1.NbRoots, aResult2.NbRoots);
for (int i = 0; i < aResult1.NbRoots; ++i)
{
EXPECT_NEAR(aResult1.Roots[i], aResult2.Roots[i], THE_TOL);
}
}
// ============================================================================
// Result structure tests
// ============================================================================
TEST(MathRoot_TrigTest, Result_BoolConversion)
{
MathRoot::TrigResult aResult = MathRoot::TrigonometricLinear(1.0, 0.0);
EXPECT_TRUE(static_cast<bool>(aResult));
EXPECT_TRUE(aResult.IsDone());
}
TEST(MathRoot_TrigTest, Result_InfiniteRoots)
{
MathRoot::TrigResult aResult = MathRoot::Trigonometric(0.0, 0.0, 0.0, 0.0, 0.0);
EXPECT_TRUE(aResult.InfiniteRoots);
// IsDone() may or may not be true for infinite roots case
}
// ============================================================================
// Ellipse-related equation tests (from ExtremaPC use case)
// ============================================================================
TEST(MathRoot_TrigTest, Ellipse_PointOnMajorAxis)
{
// Equation from ellipse extrema: (b^2-a^2)/2 * 2*cos*sin + a*X*sin = 0
// For ellipse with a=20, b=10 and point at X=30, Y=0:
// (100-400)/2 * 2*cos*sin + 20*30*sin = 0
// -150 * 2*cos*sin + 600*sin = 0
// sin * (-300*cos + 600) = 0
// => sin=0 or cos=2 (impossible)
// Roots at sin=0: x = 0, PI
double aB2MinusA2Over2 = (100.0 - 400.0) / 2.0; // -150
double aAX = 20.0 * 30.0; // 600
double aBY = 0.0;
MathRoot::TrigResult aResult =
MathRoot::Trigonometric(0.0, aB2MinusA2Over2, -aBY, aAX, 0.0, 0.0, THE_2PI);
ASSERT_TRUE(aResult.IsDone());
// Should find x=0 and x=PI as roots
EXPECT_GE(aResult.NbRoots, 2);
// Verify roots
verifyRoots(aResult, 0.0, aB2MinusA2Over2, -aBY, aAX, 0.0);
}
TEST(MathRoot_TrigTest, Ellipse_GeneralPoint)
{
// Equation from ellipse extrema with point not on axis
// For ellipse with a=20, b=10 and point at X=15, Y=8:
double aB2MinusA2Over2 = (100.0 - 400.0) / 2.0; // -150
double aAX = 20.0 * 15.0; // 300
double aBY = 10.0 * 8.0; // 80
MathRoot::TrigResult aResult =
MathRoot::Trigonometric(0.0, aB2MinusA2Over2, -aBY, aAX, 0.0, 0.0, THE_2PI);
ASSERT_TRUE(aResult.IsDone());
// Verify any roots found
verifyRoots(aResult, 0.0, aB2MinusA2Over2, -aBY, aAX, 0.0);
}
// ============================================================================
// Negative parameter range tests
// ============================================================================
TEST(MathRoot_TrigTest, NegativeBounds_Basic)
{
// sin(x) = 0 in [-PI, 0]
MathRoot::TrigResult aResult = MathRoot::TrigonometricLinear(1.0, 0.0, -THE_PI, 0.0);
ASSERT_TRUE(aResult.IsDone());
// Should find x = -PI and x = 0 (or just one depending on boundary handling)
EXPECT_GE(aResult.NbRoots, 1);
for (int i = 0; i < aResult.NbRoots; ++i)
{
EXPECT_GE(aResult.Roots[i], -THE_PI - THE_TOL);
EXPECT_LE(aResult.Roots[i], THE_TOL);
}
}
TEST(MathRoot_TrigTest, NegativeBounds_CosEquation)
{
// cos(x) = 0.5 in [-PI, PI]
MathRoot::TrigResult aResult = MathRoot::Trigonometric(0.0, 0.0, 1.0, 0.0, -0.5, -THE_PI, THE_PI);
ASSERT_TRUE(aResult.IsDone());
// cos(x) = 0.5 => x = +/-PI/3
EXPECT_GE(aResult.NbRoots, 2);
verifyRoots(aResult, 0.0, 0.0, 1.0, 0.0, -0.5);
}

100
src/FoundationClasses/TKMath/MathRoot/MathRoot_Trig.hxx
View File

@@ -17,7 +17,8 @@
#include <MathUtils_Types.hxx>
#include <MathUtils_Config.hxx>
#include <MathUtils_Core.hxx>
#include <MathPoly.hxx>
#include <MathPoly_Quadratic.hxx>
#include <MathPoly_Quartic.hxx>
#include <cmath>
#include <algorithm>
@@ -164,21 +165,40 @@ inline TrigResult Trigonometric(double theA,
return aResult;
}
aZer[0] = std::acos(aVal);
aZer[1] = -aZer[0];
aNZer = 2;
double aPrincipal = std::acos(aVal);
aZer[0] = aPrincipal; // acos gives [0, PI]
aZer[1] = -aPrincipal; // Negative angle
aNZer = 2;
// For each solution, find the representative in or near the given bounds
for (size_t i = 0; i < aNZer; ++i)
{
if (aZer[i] <= -theEps)
double aAngle = aZer[i];
// Shift angle to be near the lower bound
double aK = std::floor((theInfBound - aAngle) / THE_2PI);
aAngle += (aK + 1) * THE_2PI; // Start from a value >= theInfBound - 2*PI
// Check both this value and the next period
for (int aPeriod = 0; aPeriod < 2; ++aPeriod)
{
aZer[i] = THE_2PI - std::abs(aZer[i]);
}
aZer[i] += std::trunc(aMod) * THE_2PI;
double aX = aZer[i] - aMyBorneInf;
if (aX >= -aDelta_Eps && aX <= aDelta + aDelta_Eps)
{
aResult.Roots[aResult.NbRoots++] = aZer[i];
double aTestAngle = aAngle + aPeriod * THE_2PI;
if (aTestAngle >= theInfBound - aDelta_Eps && aTestAngle <= theSupBound + aDelta_Eps)
{
// Avoid duplicates
bool aDup = false;
for (int k = 0; k < aResult.NbRoots; ++k)
{
if (std::abs(aTestAngle - aResult.Roots[k]) < theEps)
{
aDup = true;
break;
}
}
if (!aDup && aResult.NbRoots < 4)
{
aResult.Roots[aResult.NbRoots++] = aTestAngle;
}
}
}
}
return aResult;
@@ -186,12 +206,12 @@ inline TrigResult Trigonometric(double theA,
else
{
// c*cos(x) + d*sin(x) + e = 0
// Using t = tan(x/2): (e-c) + 2dt + (e+c)t^2 = 0
// Using t = tan(x/2): (e-c)*t^2 + 2d*t + (e+c) = 0
double aAA = theE - theC;
double aBB = 2.0 * theD;
double aCC = theE + theC;
MathPoly::PolyResult aPoly = MathPoly::Quadratic(aCC, aBB, aAA);
MathPoly::PolyResult aPoly = MathPoly::Quadratic(aAA, aBB, aCC);
if (!aPoly.IsDone())
{
aResult.Status = aPoly.Status;
@@ -355,26 +375,56 @@ inline TrigResult Trigonometric(double theA,
double aX = aTeta - aMyBorneInf;
if (aX >= -aDelta_Eps && aX <= aDelta + aDelta_Eps)
{
// Newton refinement
// Newton refinement with Halley's method fallback for double roots
auto aRefineRoot = [&](double theX) -> double {
constexpr int THE_MAX_ITER = 10;
constexpr double THE_TOL = 1.0e-15;
constexpr int THE_MAX_ITER = 20;
constexpr double THE_TOL = 1.0e-14;
for (int anIter = 0; anIter < THE_MAX_ITER; ++anIter)
{
double aCos = std::cos(theX);
double aSin = std::sin(theX);
double aF =
theA * aCos * aCos + 2.0 * theB * aCos * aSin + theC * aCos + theD * aSin + theE;
double aDF = -2.0 * theA * aCos * aSin + 2.0 * theB * (aCos * aCos - aSin * aSin)
- theC * aSin + theD * aCos;
double aCos = std::cos(theX);
double aSin = std::sin(theX);
double aCos2 = aCos * aCos;
double aSin2 = aSin * aSin;
double aCS = aCos * aSin;
if (std::abs(aDF) < 1.0e-30)
double aF = theA * aCos2 + 2.0 * theB * aCS + theC * aCos + theD * aSin + theE;
double aDF = -2.0 * theA * aCS + 2.0 * theB * (aCos2 - aSin2) - theC * aSin + theD * aCos;
// Check if already converged
if (std::abs(aF) < 1.0e-15)
{
break;
}
double aDelta = aF / aDF;
double aDelta;
if (std::abs(aDF) < 1.0e-10 * (std::abs(aF) + 1.0))
{
// Near double root: use Halley's method for better convergence
// F'' = -2*a*(cos^2-sin^2) - 4*b*cos*sin - c*cos - d*sin
double aD2F =
-2.0 * theA * (aCos2 - aSin2) - 4.0 * theB * aCS - theC * aCos - theD * aSin;
double aDenom = 2.0 * aDF * aDF - aF * aD2F;
if (std::abs(aDenom) < 1.0e-30)
{
// Can't improve further
break;
}
aDelta = 2.0 * aF * aDF / aDenom;
}
else
{
// Standard Newton step
aDelta = aF / aDF;
}
// Limit step size to avoid overshooting
constexpr double THE_MAX_STEP = 0.5;
if (std::abs(aDelta) > THE_MAX_STEP)
{
aDelta = (aDelta > 0) ? THE_MAX_STEP : -THE_MAX_STEP;
}
theX -= aDelta;
if (std::abs(aDelta) < THE_TOL)

1
src/FoundationClasses/TKMath/MathSys/FILES.cmake
View File

@@ -3,6 +3,7 @@ set(OCCT_MathSys_FILES_LOCATION "${CMAKE_CURRENT_LIST_DIR}")
set(OCCT_MathSys_FILES
MathSys_Newton.hxx
MathSys_Newton2D.hxx
MathSys_LevenbergMarquardt.hxx
README.md
)

335
src/FoundationClasses/TKMath/MathSys/MathSys_Newton2D.hxx Normal file
View File

@@ -0,0 +1,335 @@
// 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 _MathSys_Newton2D_HeaderFile
#define _MathSys_Newton2D_HeaderFile
#include <MathUtils_Types.hxx>
#include <MathUtils_Config.hxx>
#include <cmath>
#include <algorithm>
//! @file MathSys_Newton2D.hxx
//! @brief Optimized 2D Newton-Raphson solver for systems of 2 equations in 2 unknowns.
//!
//! This solver is specifically optimized for 2D problems like:
//! - Point-surface extrema (find u,v where gradient is zero)
//! - Curve intersection (find t1,t2 where curves meet)
//! - Surface intersection curves
//!
//! Optimizations compared to general Newton:
//! - No math_Vector/math_Matrix allocation overhead
//! - Cramer's rule for 2x2 system (faster than LU decomposition)
//! - Squared norm comparisons (avoid sqrt in convergence check)
//! - Step limiting to prevent wild oscillations
//! - Gradient descent fallback for singular Jacobian
namespace MathSys
{
using namespace MathUtils;
//! Result of 2D Newton iteration.
struct Newton2DResult
{
Status Status = Status::NotConverged; //!< Computation status
double U = 0.0; //!< Solution U coordinate
double V = 0.0; //!< Solution V coordinate
size_t NbIter = 0; //!< Number of iterations performed
double FNorm = 0.0; //!< Final |F| norm
//! Returns true if computation succeeded.
bool IsDone() const { return Status == Status::OK; }
//! Conversion to bool for convenient checking.
explicit operator bool() const { return IsDone(); }
};
//! Optimized 2D Newton-Raphson solver with bounds.
//!
//! Solves the system [F1, F2] = [0, 0] using Newton-Raphson iteration
//! with Cramer's rule for the 2x2 linear system at each step.
//!
//! The function type must provide a method:
//! @code
//! bool ValueAndJacobian(double theU, double theV,
//! double& theF1, double& theF2,
//! double& theDF1dU, double& theDF1dV,
//! double& theDF2dU, double& theDF2dV) const;
//! @endcode
//!
//! @tparam Function type with ValueAndJacobian method
//! @param theFunc function to solve (provides F and Jacobian)
//! @param theU0 initial U guess
//! @param theV0 initial V guess
//! @param theUMin U lower bound
//! @param theUMax U upper bound
//! @param theVMin V lower bound
//! @param theVMax V upper bound
//! @param theTol tolerance for convergence (|F| < tol)
//! @param theMaxIter maximum iterations (default 20)
//! @return Newton2DResult containing solution if converged
template <typename Function>
Newton2DResult Newton2D(const Function& theFunc,
double theU0,
double theV0,
double theUMin,
double theUMax,
double theVMin,
double theVMax,
double theTol,
size_t theMaxIter = 20)
{
Newton2DResult aRes;
aRes.U = theU0;
aRes.V = theV0;
// Pre-compute max step limit once
const double aMaxStep = 0.5 * std::max(theUMax - theUMin, theVMax - theVMin);
const double aTolSq = theTol * theTol;
for (size_t i = 0; i < theMaxIter; ++i)
{
aRes.NbIter = i + 1;
double aF1, aF2, aDF1dU, aDF1dV, aDF2dU, aDF2dV;
if (!theFunc.ValueAndJacobian(aRes.U, aRes.V, aF1, aF2, aDF1dU, aDF1dV, aDF2dU, aDF2dV))
{
aRes.Status = Status::NumericalError;
return aRes;
}
// Check convergence using squared norm (avoid sqrt)
const double aFNormSq = aF1 * aF1 + aF2 * aF2;
if (aFNormSq < aTolSq)
{
aRes.FNorm = std::sqrt(aFNormSq);
aRes.Status = Status::OK;
return aRes;
}
// Solve 2x2 linear system: J * [dU, dV]^T = -[F1, F2]^T
// J = [[DF1dU, DF1dV], [DF2dU, DF2dV]]
const double aDet = aDF1dU * aDF2dV - aDF1dV * aDF2dU;
double aDU, aDV;
if (std::abs(aDet) < 1.0e-30)
{
// Singular Jacobian - try gradient descent step
const double aNormSq = aDF1dU * aDF1dU + aDF1dV * aDF1dV + aDF2dU * aDF2dU + aDF2dV * aDF2dV;
if (aNormSq < 1.0e-60)
{
aRes.FNorm = std::sqrt(aFNormSq);
aRes.Status = Status::Singular;
return aRes;
}
const double aStep = std::sqrt(aFNormSq / aNormSq) * 0.1;
// Use steepest descent direction: -J^T * F
aDU = -aStep * (aDF1dU * aF1 + aDF2dU * aF2);
aDV = -aStep * (aDF1dV * aF1 + aDF2dV * aF2);
}
else
{
// Cramer's rule (optimized: compute inverse determinant once)
const double aInvDet = 1.0 / aDet;
aDU = (-aF1 * aDF2dV + aF2 * aDF1dV) * aInvDet;
aDV = (-aF2 * aDF1dU + aF1 * aDF2dU) * aInvDet;
// Limit step size (optimized: avoid sqrt when possible)
const double aStepNormSq = aDU * aDU + aDV * aDV;
if (aStepNormSq > aMaxStep * aMaxStep)
{
const double aScale = aMaxStep / std::sqrt(aStepNormSq);
aDU *= aScale;
aDV *= aScale;
}
}
// Update and clamp to bounds
aRes.U = std::max(theUMin, std::min(theUMax, aRes.U + aDU));
aRes.V = std::max(theVMin, std::min(theVMax, aRes.V + aDV));
}
// Final convergence check (only at max iterations)
double aF1, aF2, aDF1dU, aDF1dV, aDF2dU, aDF2dV;
if (theFunc.ValueAndJacobian(aRes.U, aRes.V, aF1, aF2, aDF1dU, aDF1dV, aDF2dU, aDF2dV))
{
aRes.FNorm = std::sqrt(aF1 * aF1 + aF2 * aF2);
if (aRes.FNorm < theTol)
{
aRes.Status = Status::OK;
}
else
{
aRes.Status = Status::MaxIterations;
}
}
else
{
aRes.Status = Status::NumericalError;
}
return aRes;
}
//! Unbounded 2D Newton-Raphson solver.
//!
//! @tparam Function type with ValueAndJacobian method
//! @param theFunc function to solve
//! @param theU0 initial U guess
//! @param theV0 initial V guess
//! @param theTol tolerance for convergence
//! @param theMaxIter maximum iterations (default 20)
//! @return Newton2DResult containing solution if converged
template <typename Function>
Newton2DResult Newton2D(const Function& theFunc,
double theU0,
double theV0,
double theTol,
size_t theMaxIter = 20)
{
constexpr double aInf = 1.0e100;
return Newton2D(theFunc, theU0, theV0, -aInf, aInf, -aInf, aInf, theTol, theMaxIter);
}
//! Optimized 2D Newton solver for symmetric Jacobian.
//!
//! Many 2D problems have symmetric Jacobians (e.g., gradient of scalar function).
//! This variant is optimized for the symmetric case where J12 = J21.
//!
//! The function type must provide a method:
//! @code
//! bool ValueAndJacobian(double theU, double theV,
//! double& theF1, double& theF2,
//! double& theJ11, double& theJ12, double& theJ22) const;
//! @endcode
//! where J = [[J11, J12], [J12, J22]] (symmetric).
//!
//! @tparam Function type with symmetric ValueAndJacobian method
//! @param theFunc function to solve
//! @param theU0 initial U guess
//! @param theV0 initial V guess
//! @param theUMin U lower bound
//! @param theUMax U upper bound
//! @param theVMin V lower bound
//! @param theVMax V upper bound
//! @param theTol tolerance for convergence
//! @param theMaxIter maximum iterations (default 20)
//! @return Newton2DResult containing solution if converged
template <typename Function>
Newton2DResult Newton2DSymmetric(const Function& theFunc,
double theU0,
double theV0,
double theUMin,
double theUMax,
double theVMin,
double theVMax,
double theTol,
size_t theMaxIter = 20)
{
Newton2DResult aRes;
aRes.U = theU0;
aRes.V = theV0;
// Pre-compute max step limit once
const double aMaxStep = 0.5 * std::max(theUMax - theUMin, theVMax - theVMin);
const double aTolSq = theTol * theTol;
for (size_t i = 0; i < theMaxIter; ++i)
{
aRes.NbIter = i + 1;
double aF1, aF2, aJ11, aJ12, aJ22;
if (!theFunc.ValueAndJacobian(aRes.U, aRes.V, aF1, aF2, aJ11, aJ12, aJ22))
{
aRes.Status = Status::NumericalError;
return aRes;
}
// Check convergence using squared norm (avoid sqrt)
const double aFNormSq = aF1 * aF1 + aF2 * aF2;
if (aFNormSq < aTolSq)
{
aRes.FNorm = std::sqrt(aFNormSq);
aRes.Status = Status::OK;
return aRes;
}
// Solve 2x2 symmetric linear system: J * [dU, dV]^T = -[F1, F2]^T
// J = [[J11, J12], [J12, J22]] (symmetric)
const double aDet = aJ11 * aJ22 - aJ12 * aJ12;
double aDU, aDV;
if (std::abs(aDet) < 1.0e-30)
{
// Singular Jacobian - try gradient descent step
const double aNormSq = aJ11 * aJ11 + 2.0 * aJ12 * aJ12 + aJ22 * aJ22;
if (aNormSq < 1.0e-60)
{
aRes.FNorm = std::sqrt(aFNormSq);
aRes.Status = Status::Singular;
return aRes;
}
const double aStep = std::sqrt(aFNormSq / aNormSq) * 0.1;
aDU = -aStep * aF1;
aDV = -aStep * aF2;
}
else
{
// Cramer's rule for symmetric matrix
const double aInvDet = 1.0 / aDet;
aDU = (-aF1 * aJ22 + aF2 * aJ12) * aInvDet;
aDV = (-aF2 * aJ11 + aF1 * aJ12) * aInvDet;
// Limit step size
const double aStepNormSq = aDU * aDU + aDV * aDV;
if (aStepNormSq > aMaxStep * aMaxStep)
{
const double aScale = aMaxStep / std::sqrt(aStepNormSq);
aDU *= aScale;
aDV *= aScale;
}
}
// Update and clamp to bounds
aRes.U = std::max(theUMin, std::min(theUMax, aRes.U + aDU));
aRes.V = std::max(theVMin, std::min(theVMax, aRes.V + aDV));
}
// Final convergence check (only at max iterations)
double aF1, aF2, aJ11, aJ12, aJ22;
if (theFunc.ValueAndJacobian(aRes.U, aRes.V, aF1, aF2, aJ11, aJ12, aJ22))
{
aRes.FNorm = std::sqrt(aF1 * aF1 + aF2 * aF2);
if (aRes.FNorm < theTol)
{
aRes.Status = Status::OK;
}
else
{
aRes.Status = Status::MaxIterations;
}
}
else
{
aRes.Status = Status::NumericalError;
}
return aRes;
}
} // namespace MathSys
#endif // _MathSys_Newton2D_HeaderFile

31
src/FoundationClasses/TKMath/MathSys/README.md
View File

@@ -22,6 +22,37 @@ At each iteration, it solves the linear system J(X) * dX = -F(X), where J is the
- `MathSys::Newton()` - Standard Newton-Raphson iteration
- `MathSys::NewtonBounded()` - Newton method with box constraints
### Optimized 2D Newton Solver (`MathSys_Newton2D.hxx`)
A highly optimized Newton solver specifically designed for 2D problems (systems of 2 equations
in 2 unknowns). Common use cases include point-surface extrema, curve intersections, and
surface-surface intersection curves.
**Features:**
- Zero allocation overhead (no math_Vector/math_Matrix)
- Cramer's rule for 2x2 linear system (faster than LU decomposition)
- Squared norm comparisons (avoids sqrt in convergence check)
- Step limiting to prevent wild oscillations
- Gradient descent fallback for singular Jacobian
- Symmetric Jacobian variant for gradient-based problems
**Functions:**
- `MathSys::Newton2D()` - General 2D Newton solver
- `MathSys::Newton2DSymmetric()` - Optimized for symmetric Jacobians (e.g., gradient of scalar function)
**Usage Example (Point-Surface Extrema):**
```cpp
// Function with symmetric Jacobian: F = [dD/dU, dD/dV] where D = ||S(u,v) - P||^2
auto aResult = MathSys::Newton2DSymmetric(myDistFunc, u0, v0,
uMin, uMax, vMin, vMax,
tolerance, maxIter);
if (aResult.IsDone())
{
double u = aResult.U;
double v = aResult.V;
}
```
### Levenberg-Marquardt Algorithm (`MathSys_LevenbergMarquardt.hxx`)
The Levenberg-Marquardt algorithm is designed for nonlinear least squares problems,

1
src/FoundationClasses/TKMath/MathUtils/FILES.cmake
View File

@@ -5,6 +5,7 @@ set(OCCT_MathUtils_FILES
# Types and configuration
MathUtils_Types.hxx
MathUtils_Config.hxx
MathUtils_Domain.hxx
# Core utilities
MathUtils_Core.hxx
MathUtils_Convergence.hxx

228
src/FoundationClasses/TKMath/MathUtils/MathUtils_Domain.hxx Normal file
View File

@@ -0,0 +1,228 @@
// 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 _MathUtils_Domain_HeaderFile
#define _MathUtils_Domain_HeaderFile
#include <cmath>
//! @file MathUtils_Domain.hxx
//! @brief Parameter domain types for curves and surfaces.
//!
//! Provides lightweight value types for representing 1D and 2D parameter domains
//! with utility methods for bounds checking, clamping, and domain analysis.
namespace MathUtils
{
//! @brief 1D parameter domain for curves.
//!
//! Represents a parameter range [Min, Max] with utility methods for:
//! - Bounds checking (Contains)
//! - Parameter clamping (Clamp)
//! - Domain analysis (IsLarge, IsFullPeriod)
//!
//! @note This is a lightweight value type designed for efficiency.
//! All methods are inline and constexpr where possible.
struct Domain1D
{
double Min = 0.0; //!< Lower bound of the domain
double Max = 0.0; //!< Upper bound of the domain
//! Default constructor - creates empty domain [0, 0].
constexpr Domain1D() = default;
//! Construct from bounds.
//! @param theMin lower bound
//! @param theMax upper bound
constexpr Domain1D(double theMin, double theMax)
: Min(theMin),
Max(theMax)
{
}
//! Returns the length of the domain.
constexpr double Length() const { return Max - Min; }
//! Returns the midpoint of the domain.
constexpr double Mid() const { return (Min + Max) * 0.5; }
//! Check if value is within domain.
//! @param theU parameter value to check
//! @param theTol tolerance for boundary check
//! @return true if theU is in [Min - theTol, Max + theTol]
bool Contains(double theU, double theTol = 0.0) const
{
return theU >= Min - theTol && theU <= Max + theTol;
}
//! Clamp value to domain bounds.
//! @param theU parameter value to clamp
//! @return clamped value in [Min, Max]
double Clamp(double theU) const { return theU < Min ? Min : (theU > Max ? Max : theU); }
//! Check if domain is "large" (effectively unbounded for optimization).
//! Large domains allow skipping bounds checking for performance.
//! @param theThreshold size threshold (default 1000)
//! @return true if Length() > theThreshold
bool IsLarge(double theThreshold = 1000.0) const { return Length() > theThreshold; }
//! Check if domain covers a full periodic range.
//! @param thePeriod period of the parameter (e.g., 2*PI for angles)
//! @param theTol tolerance
//! @return true if domain covers at least one full period
bool IsFullPeriod(double thePeriod, double theTol = 1.0e-10) const
{
return Length() >= thePeriod - theTol;
}
//! Interpolate within domain.
//! @param theT interpolation parameter in [0, 1]
//! @return Min + theT * Length()
constexpr double Lerp(double theT) const { return Min + theT * (Max - Min); }
//! Normalize parameter to [0, 1] range.
//! @param theU parameter value
//! @return (theU - Min) / Length(), or 0.5 if Length() == 0
double Normalize(double theU) const
{
double aLen = Length();
return aLen > 0.0 ? (theU - Min) / aLen : 0.5;
}
//! Check if domain has finite bounds (not effectively infinite).
//! @param theInfLimit threshold for "infinity" (default 1e100)
//! @return true if both Min and Max are within finite range
bool IsFinite(double theInfLimit = 1.0e100) const
{
return Min > -theInfLimit && Max < theInfLimit;
}
};
//! @brief 2D parameter domain for surfaces.
//!
//! Represents a rectangular parameter domain [UMin, UMax] x [VMin, VMax]
//! with utility methods for:
//! - Bounds checking (Contains)
//! - Parameter clamping (Clamp)
//! - Domain analysis (IsLarge, IsFullPeriod)
//! - Access to U and V subdomains
//!
//! @note This is a lightweight value type designed for efficiency.
struct Domain2D
{
double UMin = 0.0; //!< Lower U bound
double UMax = 0.0; //!< Upper U bound
double VMin = 0.0; //!< Lower V bound
double VMax = 0.0; //!< Upper V bound
//! Default constructor - creates empty domain.
constexpr Domain2D() = default;
//! Construct from bounds.
//! @param theUMin lower U bound
//! @param theUMax upper U bound
//! @param theVMin lower V bound
//! @param theVMax upper V bound
constexpr Domain2D(double theUMin, double theUMax, double theVMin, double theVMax)
: UMin(theUMin),
UMax(theUMax),
VMin(theVMin),
VMax(theVMax)
{
}
//! Construct from two 1D domains.
//! @param theUDomain U parameter domain
//! @param theVDomain V parameter domain
constexpr Domain2D(const Domain1D& theUDomain, const Domain1D& theVDomain)
: UMin(theUDomain.Min),
UMax(theUDomain.Max),
VMin(theVDomain.Min),
VMax(theVDomain.Max)
{
}
//! Returns the U subdomain.
constexpr Domain1D U() const { return {UMin, UMax}; }
//! Returns the V subdomain.
constexpr Domain1D V() const { return {VMin, VMax}; }
//! Returns U length.
constexpr double ULength() const { return UMax - UMin; }
//! Returns V length.
constexpr double VLength() const { return VMax - VMin; }
//! Returns U midpoint.
constexpr double UMid() const { return (UMin + UMax) * 0.5; }
//! Returns V midpoint.
constexpr double VMid() const { return (VMin + VMax) * 0.5; }
//! Check if UV point is within domain.
//! @param theU U parameter
//! @param theV V parameter
//! @param theTol tolerance for boundary check
//! @return true if (theU, theV) is in domain
bool Contains(double theU, double theV, double theTol = 0.0) const
{
return theU >= UMin - theTol && theU <= UMax + theTol && theV >= VMin - theTol
&& theV <= VMax + theTol;
}
//! Clamp UV to domain bounds.
//! @param theU U parameter (modified in place)
//! @param theV V parameter (modified in place)
void Clamp(double& theU, double& theV) const
{
theU = theU < UMin ? UMin : (theU > UMax ? UMax : theU);
theV = theV < VMin ? VMin : (theV > VMax ? VMax : theV);
}
//! Check if both U and V domains are "large".
//! @param theThreshold size threshold (default 1000)
//! @return true if both U and V lengths exceed threshold
bool IsLarge(double theThreshold = 1000.0) const
{
return ULength() > theThreshold && VLength() > theThreshold;
}
//! Check if U domain covers a full period.
//! @param thePeriod period of U parameter
//! @param theTol tolerance
bool IsUFullPeriod(double thePeriod, double theTol = 1.0e-10) const
{
return ULength() >= thePeriod - theTol;
}
//! Check if V domain covers a full period.
//! @param thePeriod period of V parameter
//! @param theTol tolerance
bool IsVFullPeriod(double thePeriod, double theTol = 1.0e-10) const
{
return VLength() >= thePeriod - theTol;
}
//! Check if domain has finite bounds (not effectively infinite).
//! @param theInfLimit threshold for "infinity" (default 1e100)
bool IsFinite(double theInfLimit = 1.0e100) const
{
return UMin > -theInfLimit && UMax < theInfLimit && VMin > -theInfLimit && VMax < theInfLimit;
}
};
} // namespace MathUtils
#endif // _MathUtils_Domain_HeaderFile