ENH: Define itk::Math constants from literals and C++20 std::numbers

Replace the vnl_math:: aliases with self-contained definitions: the
correctly-rounded double literals (OEIS-referenced) are the fallback,
and when C++20 <numbers> is available the std::numbers constants, or
constexpr arithmetic on them, are used instead.

This decouples itk::Math from vnl_math, whose constants are slated for
deprecation and eventual removal in upstream vxl. itk::Math no longer
depends on those aliases remaining available.

sqrt2pi has no std::numbers counterpart and std::sqrt is not constexpr
before C++26, so it stays a literal in both paths.

Author: hjmjohnson

Modules/Core/Common/include/itkMath.h

@@ -33,6 +33,11 @@
 #include <cmath>
 #include <limits>
 #include <type_traits>
+// _MSVC_LANG tracks /std: even when __cplusplus stays 199711L (no /Zc:__cplusplus),
+// matching the STL condition under which <limits> defines __cpp_lib_math_constants.
+#if (defined(__cplusplus) && __cplusplus >= 202002L) || (defined(_MSVC_LANG) && _MSVC_LANG >= 202002L)
+#  include <numbers>
+#endif
 #include "itkMathDetail.h"
 #include "itkConceptChecking.h"
 #include <vnl/vnl_math.h>
@@ -46,59 +51,67 @@
 
 namespace itk::Math
 {
-// These constants originate from VXL's vnl_math.h. They are exposed as
-// inline constexpr namespace-scope constants so they are usable in
-// constant expressions with a single shared definition across translation
-// units.
-
+// These mathematical constants originate from VXL's vnl_math.h. The literal
+// fallback values are the correctly-rounded double-precision constants (OEIS
+// references in vnl_math.h). When C++20 <numbers> is available, the standard
+// constants (or constexpr arithmetic on them) are used instead.
+#if defined(__cpp_lib_math_constants)
+#  define itkMATH_STD_OR_LITERAL(std_value, literal_value) (std_value)
+#else
+#  define itkMATH_STD_OR_LITERAL(std_value, literal_value) (literal_value)
+#endif
 
 /** \brief \f[e\f] The base of the natural logarithm or Euler's number */
-inline constexpr double e = vnl_math::e;
+inline constexpr double e = itkMATH_STD_OR_LITERAL(std::numbers::e, 2.71828182845904523536);
 /** \brief  \f[ \log_2 e \f] */
-inline constexpr double log2e = vnl_math::log2e;
+inline constexpr double log2e = itkMATH_STD_OR_LITERAL(std::numbers::log2e, 1.44269504088896340736);
 /** \brief \f[ \log_{10} e \f] */
-inline constexpr double log10e = vnl_math::log10e;
+inline constexpr double log10e = itkMATH_STD_OR_LITERAL(std::numbers::log10e, 0.43429448190325182765);
 /** \brief \f[ \log_e 2 \f] */
-inline constexpr double ln2 = vnl_math::ln2;
+inline constexpr double ln2 = itkMATH_STD_OR_LITERAL(std::numbers::ln2, 0.69314718055994530942);
 /** \brief \f[ \log_e 10 \f] */
-inline constexpr double ln10 = vnl_math::ln10;
+inline constexpr double ln10 = itkMATH_STD_OR_LITERAL(std::numbers::ln10, 2.30258509299404568402);
 /** \brief \f[ \pi \f]  */
-inline constexpr double pi = vnl_math::pi;
+inline constexpr double pi = itkMATH_STD_OR_LITERAL(std::numbers::pi, 3.14159265358979323846);
 /** \brief \f[ 2\pi \f]  */
-inline constexpr double twopi = vnl_math::twopi;
+inline constexpr double twopi = itkMATH_STD_OR_LITERAL(2 * std::numbers::pi, 6.28318530717958647693);
 /** \brief \f[ \frac{\pi}{2} \f]  */
-inline constexpr double pi_over_2 = vnl_math::pi_over_2;
+inline constexpr double pi_over_2 = itkMATH_STD_OR_LITERAL(std::numbers::pi / 2, 1.57079632679489661923);
 /** \brief \f[ \frac{\pi}{4} \f]  */
-inline constexpr double pi_over_4 = vnl_math::pi_over_4;
+inline constexpr double pi_over_4 = itkMATH_STD_OR_LITERAL(std::numbers::pi / 4, 0.78539816339744830962);
 /** \brief \f[ \frac{\pi}{180} \f]  */
-inline constexpr double pi_over_180 = vnl_math::pi_over_180;
+inline constexpr double pi_over_180 = itkMATH_STD_OR_LITERAL(std::numbers::pi / 180, 0.01745329251994329577);
 /** \brief \f[ \frac{1}{\pi} \f]  */
-inline constexpr double one_over_pi = vnl_math::one_over_pi;
+inline constexpr double one_over_pi = itkMATH_STD_OR_LITERAL(std::numbers::inv_pi, 0.31830988618379067154);
 /** \brief \f[ \frac{2}{\pi} \f]  */
-inline constexpr double two_over_pi = vnl_math::two_over_pi;
+inline constexpr double two_over_pi = itkMATH_STD_OR_LITERAL(2 * std::numbers::inv_pi, 0.63661977236758134308);
 /** \brief \f[ \frac{180}{\pi} \f]  */
-inline constexpr double deg_per_rad = vnl_math::deg_per_rad;
+inline constexpr double deg_per_rad = itkMATH_STD_OR_LITERAL(180 * std::numbers::inv_pi, 57.2957795130823208768);
 /** \brief \f[ \sqrt{2\pi} \f]  */
-inline constexpr double sqrt2pi = vnl_math::sqrt2pi;
+inline constexpr double sqrt2pi =
+  2.50662827463100050242; // no std::numbers constant; std::sqrt not constexpr until C++26
 /** \brief \f[ \frac{2}{\sqrt{\pi}} \f]  */
-inline constexpr double two_over_sqrtpi = vnl_math::two_over_sqrtpi;
+inline constexpr double two_over_sqrtpi = itkMATH_STD_OR_LITERAL(2 * std::numbers::inv_sqrtpi, 1.12837916709551257390);
 /** \brief \f[ \frac{1}{\sqrt{2\pi}} \f]  */
-inline constexpr double one_over_sqrt2pi = vnl_math::one_over_sqrt2pi;
+inline constexpr double one_over_sqrt2pi =
+  itkMATH_STD_OR_LITERAL(std::numbers::inv_sqrtpi / std::numbers::sqrt2, 0.39894228040143267794);
 /** \brief \f[ \sqrt{2} \f]  */
-inline constexpr double sqrt2 = vnl_math::sqrt2;
+inline constexpr double sqrt2 = itkMATH_STD_OR_LITERAL(std::numbers::sqrt2, 1.41421356237309504880);
 /** \brief \f[ \sqrt{ \frac{1}{2}} \f] */
-inline constexpr double sqrt1_2 = vnl_math::sqrt1_2;
+inline constexpr double sqrt1_2 = itkMATH_STD_OR_LITERAL(std::numbers::sqrt2 / 2, 0.70710678118654752440);
 /** \brief \f[ \sqrt{ \frac{1}{3}} \f] */
-inline constexpr double sqrt1_3 = vnl_math::sqrt1_3;
+inline constexpr double sqrt1_3 = itkMATH_STD_OR_LITERAL(std::numbers::inv_sqrt3, 0.57735026918962576451);
 /** \brief euler constant */
-inline constexpr double euler = vnl_math::euler;
+inline constexpr double euler = itkMATH_STD_OR_LITERAL(std::numbers::egamma, 0.57721566490153286061);
+
+#undef itkMATH_STD_OR_LITERAL
 
 //: IEEE double machine precision
-inline constexpr double eps = vnl_math::eps;
-inline constexpr double sqrteps = vnl_math::sqrteps;
+inline constexpr double eps = std::numeric_limits<double>::epsilon();
+inline constexpr double sqrteps = 1.490116119384766e-08;
 //: IEEE single machine precision
-inline constexpr float float_eps = vnl_math::float_eps;
-inline constexpr float float_sqrteps = vnl_math::float_sqrteps;
+inline constexpr float float_eps = std::numeric_limits<float>::epsilon();
+inline constexpr float float_sqrteps = 3.4526698300e-4F;
 
 /** A useful macro to generate a template floating point to integer
  *  conversion templated on the return type and using either the 32

Modules/Filtering/PhaseSymmetry/include/itkSinusoidSpatialFunction.hxx

@@ -19,7 +19,7 @@
 #define itkSinusoidSpatialFunction_hxx
 
 #include <cmath>
-#include "vnl/vnl_math.h"
+#include "itkMath.h"
 
 namespace itk
 {
@@ -45,7 +45,7 @@ SinusoidSpatialFunction<TOutput, VImageDimension, TInput>::Evaluate(const TInput
   {
     frequencyTerm += this->m_Frequency[ii] * position[ii];
   }
-  frequencyTerm *= 2.0 * vnl_math::pi;
+  frequencyTerm *= 2.0 * itk::Math::pi;
   frequencyTerm += this->m_PhaseOffset;
   const double value = std::cos(frequencyTerm);
   return static_cast<TOutput>(value);

Modules/IO/NIFTI/test/itkNiftiReadWriteDirectionTest.cxx

@@ -259,7 +259,7 @@ itkNiftiReadWriteDirectionTest(int argc, char * argv[])
   const double trace = rRelative[0][0] + rRelative[1][1] + rRelative[2][2];
 
   // Calculate the angle of rotation between the two matrices
-  const double angle = std::acos((trace - 1.0) / 2.0) * 180.0 / vnl_math::pi;
+  const double angle = std::acos((trace - 1.0) / 2.0) * 180.0 / itk::Math::pi;
 
   // The angle between the two matrices will depend on the amount of shear in the sform. Some test images
   // have relatively large shear to make sure they trigger the sform correction. In practice, permissive mode