From 769f4f690d38b72b8d398ef2b4c6ca65f617481c Mon Sep 17 00:00:00 2001 From: Nick Date: Tue, 22 Jun 2021 10:38:00 -0400 Subject: [PATCH] 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 | 94 +++++++++++++ doc/math.qbk | 1 + .../math/interpolators/bilinear_uniform.hpp | 52 +++++++ .../detail/bilinear_uniform_detail.hpp | 133 ++++++++++++++++++ .../bilinear_interpolation_performance.cpp | 40 ++++++ test/Jamfile.v2 | 1 + test/bilinear_uniform_test.cpp | 109 ++++++++++++++ 7 files changed, 430 insertions(+) create mode 100644 doc/interpolators/bilinear_uniform.qbk create mode 100644 include/boost/math/interpolators/bilinear_uniform.hpp create mode 100644 include/boost/math/interpolators/detail/bilinear_uniform_detail.hpp create mode 100644 reporting/performance/bilinear_interpolation_performance.cpp create mode 100644 test/bilinear_uniform_test.cpp diff --git a/doc/interpolators/bilinear_uniform.qbk b/doc/interpolators/bilinear_uniform.qbk new file mode 100644 index 000000000..0fff2fab2 --- /dev/null +++ b/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 + +namespace boost::math::interpolators { + +template +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 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/64 13.6 ns +Bilinear/128 13.3 ns +Bilinear/256 13.9 ns +Bilinear/512 13.7 ns +Bilinear/1024 13.2 ns +Bilinear/2048 13.1 ns +Bilinear/4096 13.2 ns +Bilinear/8192 13.2 ns +``` + + +[endsect] +[/section:bilinear_uniform] diff --git a/doc/math.qbk b/doc/math.qbk index 3d8677b99..d5f1525e0 100644 --- a/doc/math.qbk +++ b/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] diff --git a/include/boost/math/interpolators/bilinear_uniform.hpp b/include/boost/math/interpolators/bilinear_uniform.hpp new file mode 100644 index 000000000..3d3b438f5 --- /dev/null +++ b/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 +#include + +namespace boost::math::interpolators { + +template +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>(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 const & bu) { + out << *bu.m_imp; + return out; + } + +private: + std::shared_ptr> m_imp; +}; + +} +#endif diff --git a/include/boost/math/interpolators/detail/bilinear_uniform_detail.hpp b/include/boost/math/interpolators/detail/bilinear_uniform_detail.hpp new file mode 100644 index 000000000..c5114272f --- /dev/null +++ b/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 +#include +#include + +namespace boost::math::interpolators::detail { + +template +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::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::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::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::epsilon()) { + return fhi; + } + + auto bottom_left = fieldData_[idx + cols_]; + Real flo; + if (alpha <= 2*s0*std::numeric_limits::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 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 diff --git a/reporting/performance/bilinear_interpolation_performance.cpp b/reporting/performance/bilinear_interpolation_performance.cpp new file mode 100644 index 000000000..180697d1b --- /dev/null +++ b/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 +#include +#include +#include + +using boost::math::interpolators::bilinear_uniform; + +template +void Bilinear(benchmark::State& state) +{ + std::random_device rd; + std::mt19937_64 mt(rd()); + std::uniform_real_distribution unif(0, 10); + + std::vector v(state.range(0)*state.range(0), std::numeric_limits::quiet_NaN()); + + for (auto& x : v) { + x = unif(mt); + } + + auto ub = bilinear_uniform(std::move(v), state.range(0), state.range(0)); + Real x = static_cast(unif(mt)); + Real y = static_cast(unif(mt)); + for (auto _ : state) + { + benchmark::DoNotOptimize(ub(x, y)); + x += std::numeric_limits::epsilon(); + y += std::numeric_limits::epsilon(); + } + state.SetComplexityN(state.range(0)); +} + +BENCHMARK_TEMPLATE(Bilinear, double)->RangeMultiplier(2)->Range(1 << 6, 1 << 20)->Complexity(); + +BENCHMARK_MAIN(); diff --git a/test/Jamfile.v2 b/test/Jamfile.v2 index f3a3c2501..f6d2c4aea 100644 --- a/test/Jamfile.v2 +++ b/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" : -lquadmath ] ] [ run quintic_hermite_test.cpp : : : [ requires cxx17_if_constexpr cxx17_std_apply ] [ check-target-builds ../config//has_float128 "GCC libquadmath and __float128 support" : -lquadmath ] ] [ run cubic_hermite_test.cpp : : : [ requires cxx17_if_constexpr cxx17_std_apply ] [ check-target-builds ../config//has_float128 "GCC libquadmath and __float128 support" : -lquadmath ] ] + [ run bilinear_uniform_test.cpp : : : [ requires cxx17_if_constexpr cxx17_std_apply ] ] [ run catmull_rom_test.cpp ../../test/build//boost_unit_test_framework : : : 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 : : : 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 : : : TEST=3 [ requires cxx11_hdr_array cxx11_hdr_initializer_list ] : catmull_rom_test_3 ] diff --git a/test/bilinear_uniform_test.cpp b/test/bilinear_uniform_test.cpp new file mode 100644 index 000000000..fb1d3d502 --- /dev/null +++ b/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 +#include +#include +#include +#include +#ifdef BOOST_HAS_FLOAT128 +#include +using boost::multiprecision::float128; +#endif + +using boost::math::interpolators::bilinear_uniform; + +template +void test_four_values() +{ + Real x0 = 0; + Real y0 = 0; + Real dx = 1; + Real dy = 1; + Real value = 1.5; + std::vector v(2*2, value); + auto v_copy = v; + auto ub = bilinear_uniform(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 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(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 +void test_linear() +{ + std::random_device rd; + std::uniform_real_distribution dis(1,2); + std::array 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 v(rows*cols, std::numeric_limits::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(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(); + test_four_values(); + test_four_values(); + test_linear(); + return boost::math::test::report_errors(); +}