ENH: Native Simpson, deprecate vnl integral classes, retire netlib quadrature

Modules/ThirdParty/VNL/src/CMakeLists.txt

@@ -26,7 +26,6 @@ endforeach()
 
 foreach(
   exe
-  netlib_integral_test
   netlib_lbfgs_example
   netlib_lbfgsb_example
   netlib_slamch_test

Modules/ThirdParty/VNL/src/vxl/core/vnl/algo/tests/test_integral.cxx

@@ -1,3 +1,4 @@
+#define ITK_LEGACY_TEST
 // not used? #include <iostream>
 #include <cmath>
 #include <utility>

Modules/ThirdParty/VNL/src/vxl/core/vnl/algo/vnl_adaptsimpson_integral.h

@@ -1,5 +1,18 @@
 #ifndef VNL_ADAPTSIMPSON_INTEGRAL_H_
 #define VNL_ADAPTSIMPSON_INTEGRAL_H_
+
+#if __has_include(<itkConfigure.h>)
+#  include <itkConfigure.h>
+#  if defined(ITK_FUTURE_LEGACY_REMOVE)
+#    error "vnl_adaptsimpson_integral was removed; supply an adaptive Simpson rule directly."
+#  elif defined(ITK_LEGACY_REMOVE) && !defined(ITK_LEGACY_SILENT) && !defined(ITK_LEGACY_TEST)
+#    if defined(_MSC_VER)
+#      pragma message("vnl_adaptsimpson_integral is deprecated; supply an adaptive Simpson rule directly.")
+#    else
+#      warning "vnl_adaptsimpson_integral is deprecated; supply an adaptive Simpson rule directly."
+#    endif
+#  endif
+#endif
 //:
 // \file
 // \author Kongbin Kang at Brown

Modules/ThirdParty/VNL/src/vxl/core/vnl/algo/vnl_simpson_integral.cxx

@@ -1,5 +1,4 @@
 #include "vnl_simpson_integral.h"
-#include <vnl/algo/vnl_netlib.h>
 
 double
 vnl_simpson_integral::int_fnct_(double * x)
@@ -10,13 +9,25 @@ vnl_simpson_integral::int_fnct_(double * x)
 double
 vnl_simpson_integral::integral(vnl_integrant_fnct * f, double a, double b, long n)
 {
-
-  double res = 0;
-
   // set the function
   pfnct_ = f;
 
-  v3p_netlib_simpru_(&vnl_simpson_integral::int_fnct_, &a, &b, &n, &res);
+  // Composite Simpson rule (Mathews Algorithm 7.2): n intervals, 2n subintervals.
+  const double h = (b - a) / (2 * n);
+
+  double sumeven = 0.0;
+  for (long k = 1; k <= n - 1; ++k)
+  {
+    double x = a + h * 2 * k;
+    sumeven += int_fnct_(&x);
+  }
+
+  double sumodd = 0.0;
+  for (long k = 1; k <= n; ++k)
+  {
+    double x = a + h * (2 * k - 1);
+    sumodd += int_fnct_(&x);
+  }
 
-  return res;
+  return h * (int_fnct_(&a) + int_fnct_(&b) + 2 * sumeven + 4 * sumodd) / 3;
 }

Modules/ThirdParty/VNL/src/vxl/core/vnl/algo/vnl_simpson_integral.h

@@ -1,5 +1,18 @@
 #ifndef VNL_SIMPSON_INTEGRAL_H_
 #define VNL_SIMPSON_INTEGRAL_H_
+
+#if __has_include(<itkConfigure.h>)
+#  include <itkConfigure.h>
+#  if defined(ITK_FUTURE_LEGACY_REMOVE)
+#    error "vnl_simpson_integral was removed; supply a composite Simpson rule (Mathews Algorithm 7.2) directly."
+#  elif defined(ITK_LEGACY_REMOVE) && !defined(ITK_LEGACY_SILENT) && !defined(ITK_LEGACY_TEST)
+#    if defined(_MSC_VER)
+#      pragma message("vnl_simpson_integral is deprecated; supply a composite Simpson rule directly.")
+#    else
+#      warning "vnl_simpson_integral is deprecated; supply a composite Simpson rule directly."
+#    endif
+#  endif
+#endif
 //:
 // \file
 // \author Kongbin Kang at Brown

Modules/ThirdParty/VNL/src/vxl/v3p/netlib/CMakeLists.txt

@@ -257,10 +257,6 @@ set(V3P_NETLIB_linalg_SOURCES
      linalg/lsmrBase.cxx linalg/lsmrBase.h
      linalg/lsmrDense.cxx linalg/lsmrDense.h
    )
-set(V3P_NETLIB_mathews_SOURCES
-     mathews/simpson.c mathews/simpson.h
-     mathews/trapezod.c mathews/trapezod.h
-   )
 set(V3P_NETLIB_sparse_SOURCES
      sparse/spAllocate.c
      sparse/spBuild.c
@@ -285,7 +281,6 @@ set(v3p_netlib_sources
      ${V3P_NETLIB_minpack_SOURCES}
      ${V3P_NETLIB_opt_SOURCES}
      ${V3P_NETLIB_linalg_SOURCES}
-     ${V3P_NETLIB_mathews_SOURCES}
      ${V3P_NETLIB_sparse_SOURCES}
    )
 

Modules/ThirdParty/VNL/src/vxl/v3p/netlib/README.md

@@ -73,7 +73,7 @@ respecitve sources.
 5. Convert the sources to C using `f2c` and replace the `f2c.h` header
 inclusion with `v3p_netlib.h`:
     ```bash
-    for d in blas linpack temperton eispack laso arpack lapack/complex16 lapack/double lapack/single lapack/util napack minpack opt linalg toms datapac mathews; do
+    for d in blas linpack temperton eispack laso arpack lapack/complex16 lapack/double lapack/single lapack/util napack minpack opt linalg toms datapac; do
       for f in ${d}/*.f; do
         b=`echo "$f" | sed 's/.f$//'`
         if [ ! -f "${b}.c" ]; then
@@ -120,7 +120,7 @@ be left alone.
 with the mangled interface.  Include them all in `v3p_netlib_prototypes.h`.
     ```bash
     echo "/* Include prototype headers. */" > v3p_netlib_prototypes.h
-    for f in blas/*.P linpack/*.P temperton/*.P eispack/*.P laso/*.P arpack/*.P lapack/*/*.P napack/*.P minpack/*.P opt/*.P linalg/*.P toms/*.P datapac/*.P mathews/*.P; do
+    for f in blas/*.P linpack/*.P temperton/*.P eispack/*.P laso/*.P arpack/*.P lapack/*/*.P napack/*.P minpack/*.P opt/*.P linalg/*.P toms/*.P datapac/*.P; do
       b=`echo "$f" | sed 's/.P//'`
       if [ ! -f "${b}.h" ] ; then
         echo "Converting prototype $b"