From cfc686b53d79870edfabfdfd6ff67263f34c5faf Mon Sep 17 00:00:00 2001
From: Thomas Fischer <thomas.fischer@ufz.de>
Date: Thu, 17 Jan 2013 13:20:21 +0100
Subject: [PATCH] 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.
---
.../PiecewiseLinearInterpolation.cpp | 58 ++++++++++---------
.../PiecewiseLinearInterpolation.h | 11 ++--
2 files changed, 35 insertions(+), 34 deletions(-)
diff --git a/MathLib/InterpolationAlgorithms/PiecewiseLinearInterpolation.cpp b/MathLib/InterpolationAlgorithms/PiecewiseLinearInterpolation.cpp
index 8e75dcdbf0c..f21978de258 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 975f3582fc6..df3fb2e45ae 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;
--
GitLab