From cfb78529a670c1cdf159d07a1f216f33b0911fa1 Mon Sep 17 00:00:00 2001
From: rahv <karsten.rink@ufz.de>
Date: Wed, 3 Jun 2015 11:05:39 +0200
Subject: [PATCH] fixes for issues #677 and #678
---
GeoLib/Raster.cpp | 105 +++++++++++++++++++++++++---------------------
1 file changed, 57 insertions(+), 48 deletions(-)
diff --git a/GeoLib/Raster.cpp b/GeoLib/Raster.cpp
index 7b93f3edeea..6aa37aec30d 100644
--- a/GeoLib/Raster.cpp
+++ b/GeoLib/Raster.cpp
@@ -123,58 +123,67 @@ double Raster::getValueAtPoint(const MathLib::Point3d &pnt) const
double Raster::interpolateValueAtPoint(MathLib::Point3d const& pnt) const
{
- // position in raster
- double const xPos ((pnt[0] - _ll_pnt[0]) / _cell_size);
- double const yPos ((pnt[1] - _ll_pnt[1]) / _cell_size);
- // raster cell index
- std::size_t const xIdx (static_cast<size_t>(std::floor(xPos)));
- std::size_t const yIdx (static_cast<size_t>(std::floor(yPos)));
-
- // weights for bilinear interpolation
- double const half_delta = 0.5*_cell_size;
- double const xShift = std::fabs(xPos-(xIdx+half_delta)) / _cell_size;
- double const yShift = std::fabs(yPos-(yIdx+half_delta)) / _cell_size;
- std::array<double,4> weight = {{ (1-xShift)*(1-xShift), xShift*(1-yShift), xShift*yShift, (1-xShift)*yShift }};
-
- // neightbors to include in interpolation
- int const xShiftIdx = (xPos-xIdx-half_delta>=0) ? 1 : -1;
- int const yShiftIdx = (yPos-yIdx-half_delta>=0) ? 1 : -1;
- std::array<int,4> const x_nb = {{ 0, xShiftIdx, xShiftIdx, 0 }};
- std::array<int,4> const y_nb = {{ 0, 0, yShiftIdx, yShiftIdx }};
-
- // get pixel values
- std::array<double,4> pix_val;
- unsigned no_data_count (0);
- for (unsigned j=0; j<4; ++j)
- {
- pix_val[j] = _raster_data[(yIdx + y_nb[j])*_n_cols + (xIdx + x_nb[j])];
- if (std::fabs(pix_val[j] - _no_data_val) < std::numeric_limits<double>::epsilon())
- {
- weight[j] = 0;
- no_data_count++;
- }
- }
-
- // adjust weights if necessary
- if (no_data_count > 0)
- {
- if (no_data_count == 4) // if there is absolutely no data just use the default value
- return _no_data_val;
-
- const double norm = (double)(4)/(4-no_data_count);
- std::for_each(weight.begin(), weight.end(), [&norm](double &val){val*=norm;});
- }
-
- // new value
- return MathLib::scalarProduct<double,4>(weight.data(), pix_val.data());
-}
+ // position in raster
+ double const xPos ((pnt[0] - _ll_pnt[0]) / _cell_size);
+ double const yPos ((pnt[1] - _ll_pnt[1]) / _cell_size);
+ // raster cell index
+ std::size_t const xIdx (static_cast<size_t>(std::floor(xPos)));
+ std::size_t const yIdx (static_cast<size_t>(std::floor(yPos)));
+
+ // weights for bilinear interpolation
+ double const half_delta = 0.5*_cell_size;
+ double const xShift = std::fabs(xPos-(xIdx+half_delta)) / _cell_size;
+ double const yShift = std::fabs(yPos-(yIdx+half_delta)) / _cell_size;
+ std::array<double,4> weight = {{ (1-xShift)*(1-yShift), xShift*(1-yShift), xShift*yShift, (1-xShift)*yShift }};
+
+ // neightbors to include in interpolation
+ int const xShiftIdx = (xPos-xIdx-half_delta>=0) ? 1 : -1;
+ int const yShiftIdx = (yPos-yIdx-half_delta>=0) ? 1 : -1;
+ std::array<int,4> const x_nb = {{ 0, xShiftIdx, xShiftIdx, 0 }};
+ std::array<int,4> const y_nb = {{ 0, 0, yShiftIdx, yShiftIdx }};
+
+ // get pixel values
+ std::array<double,4> pix_val;
+ unsigned no_data_count (0);
+ for (unsigned j=0; j<4; ++j)
+ {
+ // check if neighbour pixel is still on the raster, otherwise substitute a no data value
+ if ( (xIdx + x_nb[j]) < 0 ||
+ (yIdx + y_nb[j]) < 0 ||
+ (xIdx + x_nb[j]) > _ll_pnt[0]+(_n_cols*_cell_size) ||
+ (yIdx + y_nb[j]) > _ll_pnt[1]+(_n_rows*_cell_size) )
+ pix_val[j] = _no_data_val;
+ else
+ pix_val[j] = _raster_data[(yIdx + y_nb[j])*_n_cols + (xIdx + x_nb[j])];
+
+ // remove no data values
+ if (std::fabs(pix_val[j] - _no_data_val) < std::numeric_limits<double>::epsilon())
+ {
+ weight[j] = 0;
+ no_data_count++;
+ }
+ }
+
+ // adjust weights if necessary
+ if (no_data_count > 0)
+ {
+ if (no_data_count == 4) // if there is absolutely no data just use the default value
+ return _no_data_val;
+
+ const double norm = 1.0 / (weight[0]+weight[1]+weight[2]+weight[3]);
+ std::for_each(weight.begin(), weight.end(), [&norm](double &val){val*=norm;});
+ }
+
+ // new value
+ return MathLib::scalarProduct<double,4>(weight.data(), pix_val.data());
+ }
bool Raster::isPntOnRaster(MathLib::Point3d const& pnt) const
{
- if ((pnt[0]<_ll_pnt[0]) || (pnt[0]>_ll_pnt[0]+(_n_cols*_cell_size)) ||
- (pnt[1]<_ll_pnt[1]) || (pnt[1]>_ll_pnt[1]+(_n_rows*_cell_size)))
+ if ((pnt[0]<_ll_pnt[0]) || (pnt[0]>_ll_pnt[0]+(_n_cols*_cell_size)) ||
+ (pnt[1]<_ll_pnt[1]) || (pnt[1]>_ll_pnt[1]+(_n_rows*_cell_size)))
return false;
- return true;
+ return true;
}
} // end namespace GeoLib
--
GitLab