diff --git a/src/FoundationClasses/TKMath/GTests/FILES.cmake b/src/FoundationClasses/TKMath/GTests/FILES.cmake index 6ec8fc1ae2..490df12353 100644 --- a/src/FoundationClasses/TKMath/GTests/FILES.cmake +++ b/src/FoundationClasses/TKMath/GTests/FILES.cmake @@ -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 diff --git a/src/FoundationClasses/TKMath/GTests/MathRoot_Trig_Test.cxx b/src/FoundationClasses/TKMath/GTests/MathRoot_Trig_Test.cxx new file mode 100644 index 0000000000..04d8cfa284 --- /dev/null +++ b/src/FoundationClasses/TKMath/GTests/MathRoot_Trig_Test.cxx @@ -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 + +#include + +#include + +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(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); +} diff --git a/src/FoundationClasses/TKMath/MathRoot/MathRoot_Trig.hxx b/src/FoundationClasses/TKMath/MathRoot/MathRoot_Trig.hxx index 0e93fafecd..9563b91ab0 100644 --- a/src/FoundationClasses/TKMath/MathRoot/MathRoot_Trig.hxx +++ b/src/FoundationClasses/TKMath/MathRoot/MathRoot_Trig.hxx @@ -17,7 +17,8 @@ #include #include #include -#include +#include +#include #include #include @@ -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) diff --git a/src/FoundationClasses/TKMath/MathSys/FILES.cmake b/src/FoundationClasses/TKMath/MathSys/FILES.cmake index d5ec8faf59..c67864d744 100644 --- a/src/FoundationClasses/TKMath/MathSys/FILES.cmake +++ b/src/FoundationClasses/TKMath/MathSys/FILES.cmake @@ -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 ) diff --git a/src/FoundationClasses/TKMath/MathSys/MathSys_Newton2D.hxx b/src/FoundationClasses/TKMath/MathSys/MathSys_Newton2D.hxx new file mode 100644 index 0000000000..956056f3a7 --- /dev/null +++ b/src/FoundationClasses/TKMath/MathSys/MathSys_Newton2D.hxx @@ -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 +#include + +#include +#include + +//! @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 +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 +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 +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 diff --git a/src/FoundationClasses/TKMath/MathSys/README.md b/src/FoundationClasses/TKMath/MathSys/README.md index 27ca350554..ebdc81d441 100644 --- a/src/FoundationClasses/TKMath/MathSys/README.md +++ b/src/FoundationClasses/TKMath/MathSys/README.md @@ -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, diff --git a/src/FoundationClasses/TKMath/MathUtils/FILES.cmake b/src/FoundationClasses/TKMath/MathUtils/FILES.cmake index 6fbb097a42..271ef08838 100644 --- a/src/FoundationClasses/TKMath/MathUtils/FILES.cmake +++ b/src/FoundationClasses/TKMath/MathUtils/FILES.cmake @@ -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 diff --git a/src/FoundationClasses/TKMath/MathUtils/MathUtils_Domain.hxx b/src/FoundationClasses/TKMath/MathUtils/MathUtils_Domain.hxx new file mode 100644 index 0000000000..42ed92e474 --- /dev/null +++ b/src/FoundationClasses/TKMath/MathUtils/MathUtils_Domain.hxx @@ -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 + +//! @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