interpolators/detail/bilinear_uniform_detail.hpp

0001 // Copyright Nick Thompson, 2021
0002 // Use, modification and distribution are subject to the
0003 // Boost Software License, Version 1.0.
0004 // (See accompanying file LICENSE_1_0.txt
0005 // or copy at http://www.boost.org/LICENSE_1_0.txt)
0006
0007 #ifndef BOOST_MATH_INTERPOLATORS_BILINEAR_UNIFORM_DETAIL_HPP
0008 #define BOOST_MATH_INTERPOLATORS_BILINEAR_UNIFORM_DETAIL_HPP
0009
0010 #include <stdexcept>
0011 #include <iostream>
0012 #include <string>
0013 #include <limits>
0014 #include <cmath>
0015 #include <utility>
0016
0017 namespace boost::math::interpolators::detail {
0018
0019 template <class RandomAccessContainer>
0020 class bilinear_uniform_imp
0021 {
0022 public:
0023     using Real = typename RandomAccessContainer::value_type;
0024     using Z = typename RandomAccessContainer::size_type;
0025
0026     bilinear_uniform_imp(RandomAccessContainer && fieldData, Z rows, Z cols, Real dx = 1, Real dy = 1, Real x0 = 0, Real y0 = 0)
0027     {
0028         using std::to_string;
0029         if(fieldData.size() != rows*cols)
0030         {
0031             std::string err = std::string(__FILE__) + ":" + to_string(__LINE__)
0032                + " 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.";
0033             throw std::logic_error(err);
0034         }
0035         if (rows < 2) {
0036             throw std::logic_error("There must be at least two rows of data for bilinear interpolation to be well-defined.");
0037         }
0038         if (cols < 2) {
0039             throw std::logic_error("There must be at least two columns of data for bilinear interpolation to be well-defined.");
0040         }
0041
0042         fieldData_ = std::move(fieldData);
0043         rows_ = rows;
0044         cols_ = cols;
0045         x0_ = x0;
0046         y0_ = y0;
0047         dx_ = dx;
0048         dy_ = dy;
0049
0050         if (dx_ <= 0) {
0051             std::string err = std::string(__FILE__) + ":" + to_string(__LINE__) + " dx = " + to_string(dx) + ", but dx > 0 is required. Are the arguments out of order?";
0052             throw std::logic_error(err);
0053         }
0054         if (dy_ <= 0) {
0055             std::string err = std::string(__FILE__) + ":" + to_string(__LINE__) + " dy = " + to_string(dy) + ", but dy > 0 is required. Are the arguments out of order?";
0056             throw std::logic_error(err);
0057         }
0058     }
0059
0060     Real operator()(Real x, Real y) const
0061     {
0062         using std::floor;
0063         if (x > x0_ + (cols_ - 1)*dx_ || x < x0_) {
0064             std::cerr << __FILE__ << ":" << __LINE__ << ":" << __func__ << "\n";
0065             std::cerr << "Querying the bilinear_uniform interpolator at (x,y) = (" << x << ", " << y << ") is not allowed.\n";
0066             std::cerr << "x must lie in the interval [" << x0_ << ", " << x0_ + (cols_ -1)*dx_ << "]\n";
0067             return std::numeric_limits<Real>::quiet_NaN();
0068         }
0069         if (y > y0_ + (rows_ - 1)*dy_ || y < y0_) {
0070             std::cerr << __FILE__ << ":" << __LINE__ << ":" << __func__ << "\n";
0071             std::cerr << "Querying the bilinear_uniform interpolator at (x,y) = (" << x << ", " << y << ") is not allowed.\n";
0072             std::cerr << "y must lie in the interval [" << y0_ << ", " << y0_ + (rows_ -1)*dy_ << "]\n";
0073             return std::numeric_limits<Real>::quiet_NaN();
0074         }
0075
0076         Real s = (x - x0_)/dx_;
0077         Real s0 = floor(s);
0078         Real t = (y - y0_)/dy_;
0079         Real t0 = floor(t);
0080         auto xidx = static_cast<Z>(s0);
0081         auto yidx = static_cast<Z>(t0);
0082         Z idx = yidx*cols_  + xidx;
0083         Real alpha = s - s0;
0084         Real beta = t - t0;
0085
0086         Real fhi;
0087         // If alpha = 0, then we can segfault by reading fieldData_[idx+1]:
0088         if (alpha <= 2*s0*std::numeric_limits<Real>::epsilon())  {
0089             fhi = fieldData_[idx];
0090         } else {
0091             fhi = (1 - alpha)*fieldData_[idx] + alpha*fieldData_[idx + 1];
0092         }
0093
0094         // Again, we can get OOB access without this check.
0095         // This corresponds to interpolation over a line segment aligned with the axes.
0096         if (beta <= 2*t0*std::numeric_limits<Real>::epsilon()) {
0097             return fhi;
0098         }
0099
0100         auto bottom_left = fieldData_[idx + cols_];
0101         Real flo;
0102         if (alpha <= 2*s0*std::numeric_limits<Real>::epsilon()) {
0103             flo = bottom_left;
0104         }
0105         else {
0106             flo = (1 - alpha)*bottom_left + alpha*fieldData_[idx + cols_ + 1];
0107         }
0108         // Convex combination over vertical to get the value:
0109         return (1 - beta)*fhi + beta*flo;
0110     }
0111
0112     friend std::ostream& operator<<(std::ostream& out, bilinear_uniform_imp<RandomAccessContainer> const & bu) {
0113         out << "(x0, y0) = (" << bu.x0_ << ", " << bu.y0_ << "), (dx, dy) = (" << bu.dx_ << ", " << bu.dy_ << "), ";
0114         out << "(xf, yf) = (" << bu.x0_ + (bu.cols_ - 1)*bu.dx_ << ", " << bu.y0_ + (bu.rows_ - 1)*bu.dy_ << ")\n";
0115         for (Z j = 0; j < bu.rows_; ++j) {
0116             out << "{";
0117             for (Z i = 0; i < bu.cols_ - 1; ++i) {
0118                 out << bu.fieldData_[j*bu.cols_ + i] << ", ";
0119             }
0120             out << bu.fieldData_[j*bu.cols_ + bu.cols_ - 1] << "}\n";
0121         }
0122         return out;
0123     }
0124
0125 private:
0126     RandomAccessContainer fieldData_;
0127     Z rows_;
0128     Z cols_;
0129     Real x0_;
0130     Real y0_;
0131     Real dx_;
0132     Real dy_;
0133 };
0134
0135
0136 }
0137 #endif