Fixed a problem in PiecewiseLinearInterpolation.

If the point to interpolate the value on is not located in any interval we choose either the first or the last interval and extrapolate the data. The algorithm for the determination of the interval is corrected.

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;
 }

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;