diff --git a/MathLib/InterpolationAlgorithms/PiecewiseLinearInterpolation.cpp b/MathLib/InterpolationAlgorithms/PiecewiseLinearInterpolation.cpp index 8e75dcdbf0cce305ebc1363dd331a23638df9819..f21978de258d105d5c0204c2eeff859a29e747b7 100644 --- a/MathLib/InterpolationAlgorithms/PiecewiseLinearInterpolation.cpp +++ b/MathLib/InterpolationAlgorithms/PiecewiseLinearInterpolation.cpp @@ -11,57 +11,59 @@ * http://www.opengeosys.org/project/license * */ +#include <cmath> +#include <limits> #include "PiecewiseLinearInterpolation.h" -#include <iostream> - namespace MathLib { PiecewiseLinearInterpolation::PiecewiseLinearInterpolation(const std::vector<double>& supporting_points, - const std::vector<double>& values_at_supp_pnts) + const std::vector<double>& values_at_supp_pnts) : _supporting_points (supporting_points), _values_at_supp_pnts (values_at_supp_pnts) {} -PiecewiseLinearInterpolation::PiecewiseLinearInterpolation(const std::vector<double>& supporting_points, - const std::vector<double>& values_at_supp_pnts, - const std::vector<double>& points_to_interpolate, - std::vector<double>& values_at_interpol_pnts) +PiecewiseLinearInterpolation::PiecewiseLinearInterpolation( + const std::vector<double>& supporting_points, + const std::vector<double>& + values_at_supp_pnts, + const std::vector<double>& + points_to_interpolate, + std::vector<double>& + values_at_interpol_pnts) : _supporting_points (supporting_points), _values_at_supp_pnts (values_at_supp_pnts) { -// std::cout << "PiecewiseLinearInterpolation::PiecewiseLinearInterpolation support_points, values_at_supp_pnts: " << std::endl; -// for (size_t k(0); k<supporting_points.size(); k++) { -// std::cout << supporting_points[k] << " " << values_at_supp_pnts[k] << std::endl; -// } -// std::cout << std::endl; values_at_interpol_pnts.clear(); - for (size_t k(0); k < points_to_interpolate.size(); k++) + for (std::size_t k(0); k < points_to_interpolate.size(); k++) values_at_interpol_pnts.push_back (this->getValue (points_to_interpolate[k])); } PiecewiseLinearInterpolation::~PiecewiseLinearInterpolation() {} -double PiecewiseLinearInterpolation::getValue ( double pnt_to_interpolate ) +double PiecewiseLinearInterpolation::getValue(double pnt_to_interpolate) { // search interval that has the point inside - size_t interval_idx (std::numeric_limits<size_t>::max()); - for (size_t k(1); - k < _supporting_points.size() && interval_idx == std::numeric_limits<size_t>::max(); - k++) - if (_supporting_points[k - 1] <= pnt_to_interpolate && pnt_to_interpolate <= - _supporting_points[k]) + std::size_t interval_idx(std::numeric_limits<std::size_t>::max()); + for (std::size_t k(1); k < _supporting_points.size() + && interval_idx == std::numeric_limits<std::size_t>::max(); k++) + if (_supporting_points[k - 1] <= pnt_to_interpolate + && pnt_to_interpolate <= _supporting_points[k]) interval_idx = k - 1; + if (interval_idx == std::numeric_limits<std::size_t>::max()) { + const double dist_first(fabs(_supporting_points[0] - pnt_to_interpolate)); + const double dist_last(fabs(_supporting_points[_supporting_points.size() - 1] - pnt_to_interpolate)); + if (dist_first < dist_last) + interval_idx = 0; + else + interval_idx = _supporting_points.size() - 2; + } + // compute linear interpolation polynom: y = m * x + n - long double m ( - (_values_at_supp_pnts[interval_idx + - 1] - - _values_at_supp_pnts[interval_idx]) / - (_supporting_points[interval_idx + 1] - _supporting_points[interval_idx])); -// double m ((_values_at_supp_pnts[interval_idx] - _values_at_supp_pnts[interval_idx+1]) / (_supporting_points[interval_idx] - _supporting_points[interval_idx+1])); -// double n (_values_at_supp_pnts[interval_idx+1] - m * _supporting_points[interval_idx+1]); - long double n (_values_at_supp_pnts[interval_idx] - m * _supporting_points[interval_idx]); + long double m((_values_at_supp_pnts[interval_idx + 1] - _values_at_supp_pnts[interval_idx]) + / (_supporting_points[interval_idx + 1] - _supporting_points[interval_idx])); + long double n(_values_at_supp_pnts[interval_idx] - m * _supporting_points[interval_idx]); return m * pnt_to_interpolate + n; } diff --git a/MathLib/InterpolationAlgorithms/PiecewiseLinearInterpolation.h b/MathLib/InterpolationAlgorithms/PiecewiseLinearInterpolation.h index 975f3582fc6a8ee705a75e8a3c48079013bb486e..df3fb2e45ae37f356024247b088bc821ecfada7b 100644 --- a/MathLib/InterpolationAlgorithms/PiecewiseLinearInterpolation.h +++ b/MathLib/InterpolationAlgorithms/PiecewiseLinearInterpolation.h @@ -15,7 +15,6 @@ #ifndef PIECEWISELINEARINTERPOLATION_H_ #define PIECEWISELINEARINTERPOLATION_H_ -#include <limits> #include <vector> namespace MathLib @@ -24,14 +23,14 @@ class PiecewiseLinearInterpolation { public: PiecewiseLinearInterpolation(const std::vector<double>& supporting_points, - const std::vector<double>& values_at_supp_pnts); + const std::vector<double>& values_at_supp_pnts); PiecewiseLinearInterpolation(const std::vector<double>& supporting_points, - const std::vector<double>& values_at_supp_pnts, - const std::vector<double>& points_to_interpolate, - std::vector<double>& values_at_interpol_pnts); + const std::vector<double>& values_at_supp_pnts, + const std::vector<double>& points_to_interpolate, + std::vector<double>& values_at_interpol_pnts); virtual ~PiecewiseLinearInterpolation(); - double getValue ( double pnt_to_interpolate ); + double getValue(double pnt_to_interpolate); private: const std::vector<double>& _supporting_points;