Interpolate a uniform grid with a bilinear function. (#643)

* Interpolate a uniform grid with a bilinear function.

* Typo removal.

* Invalid syntax in Jamfile.

* Do domain verification before computation.

* Fix OOB access on print.

* pimpl the class so it can be shared between threads.

* Add google/benchmark file to measure the performance of the bilinear interpolation.

* Fix up docs.

* Remove non-ASCII characters from print statements. Add a float128 test.

* Improve the documentation of the bilinear uniform class.

* Remove float128 as it doesn't support to_string.

* Don't use decltype(fieldData.size()) as the indexer; that makes MSVC 14.2 choke. Use RandomAccessContainer::size_type.

* Use ADL for to_string for compatibility with multiprecision.

* Improve error message which rows*cols != fieldData.size().

doc/interpolators/bilinear_uniform.qbk

@@ -0,0 +1,94 @@
+[/
+Copyright (c) 2021 Nick Thompson
+Use, modification and distribution are subject to the
+Boost Software License, Version 1.0. (See accompanying file
+LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+]
+
+[section:bilinear_uniform Bilinear Uniform Interpolation]
+
+[heading Synopsis]
+
+``
+#include <boost/math/interpolators/bilinear_uniform.hpp>
+
+namespace boost::math::interpolators {
+
+template <class RandomAccessContainer>
+class bilinear_uniform
+{
+public:
+    using Real = typename RandomAccessContainer::value_type;
+
+    bilinear_uniform(RandomAccessContainer && fieldData, decltype(fieldData.size()) rows, decltype(fieldData.size()) cols, Real dx = 1, Real dy = 1, Real x0 = 0, Real y0 = 0)
+
+    Real operator()(Real x, Real y) const;
+};
+``
+
+
+[heading Bilinear Uniform Interpolation]
+
+The bilinear uniform interpolator takes a grid of data points specified by a linear index and interpolates between each segment using a bilinear function.
+Note that "bilinear" does not mean linear-it is the product of two linear functions.
+The interpolant is continuous and its evaluation is constant time.
+An example usage is as follows:
+
+    std::vector<double> v{0.1, 0.2, 0.3,
+                          0.4, 0.5, 0.5};
+    using boost::math::interpolators::bilinear_uniform;
+    int rows = 2;
+    int cols = 3;
+    double dx = 1;
+    double dy = 1;
+    auto bu = bilinear_uniform(std::move(v), rows, cols, dx, dy);
+    // evaluate at a point:
+    double z = bu(0.0, 0.0);
+
+Periodically, it is helpful to see what data the interpolator has.
+This can be achieved via
+
+    std::cout << ub << "\n";
+
+Note that the interpolator is pimpl'd, so that copying the class is cheap, and hence it can be shared between threads.
+(The call operator is threadsafe.)
+
+Note that the layout of the field data follows the convention found in laying out images: The first value is associated with (x0, y0), and the last value is associate with (x0 + (cols - 1)dx, y0 + (rows - 1)dy).
+This matches with how we think about laying out matrices in C order, but of course there is no canonical choice and conventions must be made.
+For example, it is traditional in image processing the associate the first field value with the center of the pixel (which would be called a cell-centered field in VTK).
+This interpolator is point-centered, in the sense that (x0,y0) is associated with value v[0], and (x0+dx, y0) associated with v[1].
+If you have cell-centered data at (0,0), then just pass (x0,y0) = (0.5, 0.5) to this interpolator.
+
+Note that this interpolator does not provide the option for a rotation.
+We regarded that as too expensive to handle in this class.
+Rotating the arguments must be performed by the user.
+
+
+[heading Complexity and Performance]
+
+The google/benchmark in `reporting/performance/bilinear_uniform_performance.cpp` demonstrates the cost of the call operator for this interpolator:
+
+```
+Run on (16 X 4300 MHz CPU s)
+CPU Caches:
+  L1 Data 32K (x8)
+  L1 Instruction 32K (x8)
+  L2 Unified 1024K (x8)
+  L3 Unified 11264K (x1)
+Load Average: 0.92, 0.64, 0.35
+--------------------------------------
+Benchmark                      Time
+--------------------------------------
+Bilinear<double>/64            13.6 ns
+Bilinear<double>/128           13.3 ns
+Bilinear<double>/256           13.9 ns
+Bilinear<double>/512           13.7 ns
+Bilinear<double>/1024          13.2 ns
+Bilinear<double>/2048          13.1 ns
+Bilinear<double>/4096          13.2 ns
+Bilinear<double>/8192          13.2 ns
+```
+
+
+[endsect]
+[/section:bilinear_uniform]

doc/math.qbk

@@ -740,6 +740,7 @@ and as a CD ISBN 0-9504833-2-X  978-0-9504833-2-0, Classification 519.2-dc22.
 [include interpolators/makima.qbk]
 [include interpolators/pchip.qbk]
 [include interpolators/quintic_hermite.qbk]
+[include interpolators/bilinear_uniform.qbk]
 [endmathpart]
 
 [mathpart quadrature Quadrature and Differentiation]

include/boost/math/interpolators/bilinear_uniform.hpp

@@ -0,0 +1,52 @@
+// Copyright Nick Thompson, 2021
+// Use, modification and distribution are subject to the
+// Boost Software License, Version 1.0.
+// (See accompanying file LICENSE_1_0.txt
+// or copy at http://www.boost.org/LICENSE_1_0.txt)
+
+// This implements bilinear interpolation on a uniform grid.
+// If dx and dy are both positive, then the (x,y) = (x0, y0) is associated with data index 0 (herein referred to as f[0])
+// The point (x0 + dx, y0) is associated with f[1], and (x0 + i*dx, y0) is associated with f[i],
+// i.e., we are assuming traditional C row major order.
+// The y coordinate increases *downward*, as is traditional in 2D computer graphics.
+// This is *not* how people generally think in numerical analysis (although it *is* how they lay out matrices).
+// Providing the capability of a grid rotation is too expensive and not ergonomic; you'll need to perform any rotations at the call level.
+
+// For clarity, the value f(x0 + i*dx, y0 + j*dy) must be stored in the f[j*cols + i] position.
+
+#ifndef BOOST_MATH_INTERPOLATORS_BILINEAR_UNIFORM_HPP
+#define BOOST_MATH_INTERPOLATORS_BILINEAR_UNIFORM_HPP
+#include <memory>
+#include <boost/math/interpolators/detail/bilinear_uniform_detail.hpp>
+
+namespace boost::math::interpolators {
+
+template <class RandomAccessContainer>
+class bilinear_uniform
+{
+public:
+    using Real = typename RandomAccessContainer::value_type;
+    using Z = typename RandomAccessContainer::size_type;
+
+    bilinear_uniform(RandomAccessContainer && fieldData, Z rows, Z cols, Real dx = 1, Real dy = 1, Real x0 = 0, Real y0 = 0)
+    : m_imp(std::make_shared<detail::bilinear_uniform_imp<RandomAccessContainer>>(std::move(fieldData), rows, cols, dx, dy, x0, y0))
+    {
+    }
+
+    Real operator()(Real x, Real y) const
+    {
+        return m_imp->operator()(x,y);
+    }
+
+
+    friend std::ostream& operator<<(std::ostream& out, bilinear_uniform<RandomAccessContainer> const & bu) {
+        out << *bu.m_imp;
+        return out;
+    }
+
+private:
+    std::shared_ptr<detail::bilinear_uniform_imp<RandomAccessContainer>> m_imp;
+};
+
+}
+#endif

include/boost/math/interpolators/detail/bilinear_uniform_detail.hpp

@@ -0,0 +1,133 @@
+// Copyright Nick Thompson, 2021
+// Use, modification and distribution are subject to the
+// Boost Software License, Version 1.0.
+// (See accompanying file LICENSE_1_0.txt
+// or copy at http://www.boost.org/LICENSE_1_0.txt)
+
+#ifndef BOOST_MATH_INTERPOLATORS_BILINEAR_UNIFORM_DETAIL_HPP
+#define BOOST_MATH_INTERPOLATORS_BILINEAR_UNIFORM_DETAIL_HPP
+#include <stdexcept>
+#include <iostream>
+#include <string>
+
+namespace boost::math::interpolators::detail {
+
+template <class RandomAccessContainer>
+class bilinear_uniform_imp
+{
+public:
+    using Real = typename RandomAccessContainer::value_type;
+    using Z = typename RandomAccessContainer::size_type;
+
+    bilinear_uniform_imp(RandomAccessContainer && fieldData, Z rows, Z cols, Real dx = 1, Real dy = 1, Real x0 = 0, Real y0 = 0)
+    {
+        using std::to_string;
+        if(fieldData.size() != rows*cols)
+        {
+            std::string err = std::string(__FILE__) + ":" + to_string(__LINE__)
+               + " The field data must have rows*cols elements. There are " + to_string(rows)  + " rows and " + to_string(cols) + " columns but " + to_string(fieldData.size()) + " elements in the field data.";
+            throw std::logic_error(err);
+        }
+        if (rows < 2) {
+            throw std::logic_error("There must be at least two rows of data for bilinear interpolation to be well-defined.");
+        }
+        if (cols < 2) {
+            throw std::logic_error("There must be at least two columns of data for bilinear interpolation to be well-defined.");
+        }
+
+        fieldData_ = std::move(fieldData);
+        rows_ = rows;
+        cols_ = cols;
+        x0_ = x0;
+        y0_ = y0;
+        dx_ = dx;
+        dy_ = dy;
+        using std::to_string;
+        if (dx_ <= 0) {
+            std::string err = std::string(__FILE__) + ":" + to_string(__LINE__) + " dx = " + to_string(dx) + ", but dx > 0 is required. Are the arguments out of order?";
+            throw std::logic_error(err);
+        }
+        if (dy_ <= 0) {
+            std::string err = std::string(__FILE__) + ":" + to_string(__LINE__) + " dy = " + to_string(dy) + ", but dy > 0 is required. Are the arguments out of order?";
+            throw std::logic_error(err);
+        }
+    }
+
+    Real operator()(Real x, Real y) const
+    {
+        using std::floor;
+        if (x > x0_ + (cols_ - 1)*dx_ || x < x0_) {
+            std::cerr << __FILE__ << ":" << __LINE__ << ":" << __func__ << "\n";
+            std::cerr << "Querying the bilinear_uniform interpolator at (x,y) = (" << x << ", " << y << ") is not allowed.\n";
+            std::cerr << "x must lie in the interval [" << x0_ << ", " << x0_ + (cols_ -1)*dx_ << "]\n";
+            return std::numeric_limits<Real>::quiet_NaN();
+        }
+        if (y > y0_ + (rows_ - 1)*dy_ || y < y0_) {
+            std::cerr << __FILE__ << ":" << __LINE__ << ":" << __func__ << "\n";
+            std::cerr << "Querying the bilinear_uniform interpolator at (x,y) = (" << x << ", " << y << ") is not allowed.\n";
+            std::cerr << "y must lie in the interval [" << y0_ << ", " << y0_ + (rows_ -1)*dy_ << "]\n";
+            return std::numeric_limits<Real>::quiet_NaN(); 
+        }
+
+        Real s = (x - x0_)/dx_;
+        Real s0 = floor(s);
+        Real t = (y - y0_)/dy_;
+        Real t0 = floor(t);
+        Z xidx = s0;
+        Z yidx = t0;
+        Z idx = yidx*cols_  + xidx;
+        Real alpha = s - s0;
+        Real beta = t - t0;
+
+        Real fhi;
+        // If alpha = 0, then we can segfault by reading fieldData_[idx+1]:
+        if (alpha <= 2*s0*std::numeric_limits<Real>::epsilon())  {
+            fhi = fieldData_[idx];
+        } else {
+            fhi = (1 - alpha)*fieldData_[idx] + alpha*fieldData_[idx + 1];
+        }
+
+        // Again, we can get OOB access without this check.
+        // This corresponds to interpolation over a line segment aligned with the axes.
+        if (beta <= 2*t0*std::numeric_limits<Real>::epsilon()) {
+            return fhi;
+        }
+
+        auto bottom_left = fieldData_[idx + cols_];
+        Real flo;
+        if (alpha <= 2*s0*std::numeric_limits<Real>::epsilon()) {
+            flo = bottom_left;
+        }
+        else {
+            flo = (1 - alpha)*bottom_left + alpha*fieldData_[idx + cols_ + 1];
+        }
+        // Convex combination over vertical to get the value:
+        return (1 - beta)*fhi + beta*flo;
+    }
+
+    friend std::ostream& operator<<(std::ostream& out, bilinear_uniform_imp<RandomAccessContainer> const & bu) {
+        out << "(x0, y0) = (" << bu.x0_ << ", " << bu.y0_ << "), (dx, dy) = (" << bu.dx_ << ", " << bu.dy_ << "), ";
+        out << "(xf, yf) = (" << bu.x0_ + (bu.cols_ - 1)*bu.dx_ << ", " << bu.y0_ + (bu.rows_ - 1)*bu.dy_ << ")\n";
+        for (Z j = 0; j < bu.rows_; ++j) {
+            out << "{";
+            for (Z i = 0; i < bu.cols_ - 1; ++i) {
+                out << bu.fieldData_[j*bu.cols_ + i] << ", ";
+            }
+            out << bu.fieldData_[j*bu.cols_ + bu.cols_ - 1] << "}\n";
+        }
+        return out;
+    }
+
+private:
+    RandomAccessContainer fieldData_;
+    Z rows_;
+    Z cols_;
+    Real x0_;
+    Real y0_;
+    Real dx_;
+    Real dy_;
+};
+
+
+}
+#endif

reporting/performance/bilinear_interpolation_performance.cpp

@@ -0,0 +1,40 @@
+//  (C) Copyright Nick Thompson 2020.
+//  Use, modification and distribution are subject to the
+//  Boost Software License, Version 1.0. (See accompanying file
+//  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+
+#include <random>
+#include <benchmark/benchmark.h>
+#include <boost/math/interpolators/bilinear_uniform.hpp>
+#include <boost/multiprecision/cpp_bin_float.hpp>
+
+using boost::math::interpolators::bilinear_uniform;
+
+template<class Real>
+void Bilinear(benchmark::State& state)
+{
+    std::random_device rd;
+    std::mt19937_64 mt(rd());
+    std::uniform_real_distribution<Real> unif(0, 10);
+
+    std::vector<double> v(state.range(0)*state.range(0), std::numeric_limits<Real>::quiet_NaN());
+
+    for (auto& x : v) {
+        x = unif(mt);
+    }
+
+    auto ub = bilinear_uniform<decltype(v)>(std::move(v), state.range(0), state.range(0));
+    Real x = static_cast<Real>(unif(mt));
+    Real y = static_cast<Real>(unif(mt));
+    for (auto _ : state)
+    {
+        benchmark::DoNotOptimize(ub(x, y));
+        x += std::numeric_limits<Real>::epsilon();
+        y += std::numeric_limits<Real>::epsilon();
+    }
+     state.SetComplexityN(state.range(0));
+}
+
+BENCHMARK_TEMPLATE(Bilinear, double)->RangeMultiplier(2)->Range(1 << 6, 1 << 20)->Complexity();
+
+BENCHMARK_MAIN();

test/Jamfile.v2

@@ -1162,6 +1162,7 @@ test-suite interpolators :
    [ run septic_hermite_test.cpp  : : :  [ requires cxx17_if_constexpr cxx17_std_apply ]  [ check-target-builds ../config//has_float128 "GCC libquadmath and __float128 support" : <linkflags>-lquadmath ] ]
    [ run quintic_hermite_test.cpp  : : :  [ requires cxx17_if_constexpr cxx17_std_apply ]  [ check-target-builds ../config//has_float128 "GCC libquadmath and __float128 support" : <linkflags>-lquadmath ] ]
    [ run cubic_hermite_test.cpp  : : :  [ requires cxx17_if_constexpr cxx17_std_apply ]  [ check-target-builds ../config//has_float128 "GCC libquadmath and __float128 support" : <linkflags>-lquadmath ] ]
+   [ run bilinear_uniform_test.cpp  : : : [ requires cxx17_if_constexpr cxx17_std_apply ]  ]
    [ run catmull_rom_test.cpp ../../test/build//boost_unit_test_framework : : : <define>TEST=1 [ requires cxx11_hdr_array cxx11_hdr_initializer_list ] : catmull_rom_test_1 ]
    [ run catmull_rom_test.cpp ../../test/build//boost_unit_test_framework : : : <define>TEST=2 [ requires cxx11_hdr_array cxx11_hdr_initializer_list ] : catmull_rom_test_2 ]
    [ run catmull_rom_test.cpp ../../test/build//boost_unit_test_framework : : : <define>TEST=3 [ requires cxx11_hdr_array cxx11_hdr_initializer_list ] : catmull_rom_test_3 ]

test/bilinear_uniform_test.cpp

@@ -0,0 +1,109 @@
+/*
+ * Copyright Nick Thompson, 2021
+ * Use, modification and distribution are subject to the
+ * Boost Software License, Version 1.0. (See accompanying file
+ * LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+ */
+
+#include "math_unit_test.hpp"
+#include <numeric>
+#include <random>
+#include <array>
+#include <boost/core/demangle.hpp>
+#include <boost/math/interpolators/bilinear_uniform.hpp>
+#ifdef BOOST_HAS_FLOAT128
+#include <boost/multiprecision/float128.hpp>
+using boost::multiprecision::float128;
+#endif
+
+using boost::math::interpolators::bilinear_uniform;
+
+template<class Real>
+void test_four_values()
+{
+    Real x0 = 0;
+    Real y0 = 0;
+    Real dx = 1;
+    Real dy = 1;
+    Real value = 1.5;
+    std::vector<Real> v(2*2, value);
+    auto v_copy = v;
+    auto ub = bilinear_uniform<decltype(v)>(std::move(v_copy), 2, 2, dx, dy, x0, y0);
+    for (Real x = x0; x <= x0 + dx; x += dx/8) {
+        for (Real y = y0; y <= y0 + dx; y += dy/8) {
+            CHECK_ULP_CLOSE(value, ub(x, y), 1);
+        }
+    }
+
+    // Now we test the unit square:
+    std::random_device rd;
+    std::uniform_real_distribution<Real> dis(1,2);
+
+    int i = 0;
+    while (i++ < 300) {
+        v[0] = dis(rd);
+        v[1] = dis(rd);
+        v[2] = dis(rd);
+        v[3] = dis(rd);
+
+        // See https://en.wikipedia.org/wiki/Bilinear_interpolation, section: Unit square
+        auto f = [&v](Real x, Real y) {
+            return v[0]*(1-x)*(1-y) + v[1]*x*(1-y) + v[2]*(1-x)*y + v[3]*x*y;
+        };
+
+        v_copy = v;
+        ub = bilinear_uniform<decltype(v_copy)>(std::move(v_copy), 2, 2, dx, dy, x0, y0);
+        for (Real x = x0; x <= x0 + dx; x += dx/16) {
+            for (Real y = y0; y <= y0 + dx; y += dy/16) {
+                CHECK_ULP_CLOSE(f(x,y), ub(x, y), 3);
+            }
+        }
+    }
+}
+
+template<typename Real>
+void test_linear()
+{
+    std::random_device rd;
+    std::uniform_real_distribution<Real> dis(1,2);
+    std::array<Real, 4> a{dis(rd), dis(rd), dis(rd), dis(rd)};
+    auto f = [&a](Real x, Real y) {
+        return a[0] + a[1]*x + a[2]*y + a[3]*x*y;
+    };
+
+    for (int rows = 2; rows < 20; ++rows) {
+        for (int cols = 2; cols < 20; ++cols) {
+            Real dx = dis(rd);
+            Real dy = dis(rd);
+            Real x0 = dis(rd);
+            Real y0 = dis(rd);
+            std::vector<Real> v(rows*cols, std::numeric_limits<Real>::quiet_NaN());
+            for (int i = 0; i < cols; ++i) {
+                for (int j = 0; j < rows; ++j) {
+                    v[j*cols + i] = f(x0 + i*dx, y0 + j*dy);
+                }
+            }
+            auto ub = bilinear_uniform<decltype(v)>(std::move(v), rows, cols, dx, dy, x0, y0);
+
+           for (Real x = x0; x < x0 + (cols-1)*dx; x += dx/8) {
+                for (Real y = y0; y < y0 + (rows-1)*dy; y += dy/8) {
+                    if (!CHECK_ULP_CLOSE(f(x,y), ub(x, y), 13)) {
+                        std::cerr << " f(" << x << ", " << y << ") = " << f(x,y) << "\n";
+                        std::cerr << "ub(" << x << ", " << y << ") = " << ub(x,y) << "\n";
+                    }
+                }
+            }
+        }
+    }
+}
+
+
+
+int main()
+{
+    test_four_values<float>();
+    test_four_values<double>();
+    test_four_values<long double>();
+    test_linear<double>();
+    return boost::math::test::report_errors();
+}