diff --git a/GeoLib/Raster.cpp b/GeoLib/Raster.cpp
index 6616360280d01ccea674f0a9d924f1912f60320d..e3bd7d347e4e397549b25e9ff69464c0d52c78e2 100644
--- a/GeoLib/Raster.cpp
+++ b/GeoLib/Raster.cpp
@@ -101,7 +101,7 @@ Raster* Raster::getRasterFromSurface(Surface const& sfc, double cell_size, doubl
return new Raster(n_cols, n_rows, ll[0], ll[1], cell_size, z_vals, z_vals+n_cols*n_rows ,-9999);
}
-double Raster::getValueAtPoint(const GeoLib::Point &pnt)
+double Raster::getValueAtPoint(const GeoLib::Point &pnt) const
{
if (pnt[0]>=_ll_pnt[0] && pnt[0]<(_ll_pnt[0]+(_cell_size*_n_cols)) &&
pnt[1]>=_ll_pnt[1] && pnt[1]<(_ll_pnt[1]+(_cell_size*_n_rows)))
@@ -119,4 +119,60 @@ double Raster::getValueAtPoint(const GeoLib::Point &pnt)
return _no_data_val;
}
+double Raster::interpolateValueAtPoint(GeoLib::Point 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>(floor(xPos)));
+ std::size_t const yIdx (static_cast<size_t>(floor(yPos)));
+
+ // weights for bilinear interpolation
+ double const half_delta = 0.5*_cell_size;
+ double const xShift = fabs(xPos-(xIdx+half_delta)) / _cell_size;
+ double const yShift = 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 (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());
+}
+
+bool Raster::isPntOnRaster(GeoLib::Point 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)))
+ return false;
+ return true;
+}
+
} // end namespace GeoLib
diff --git a/GeoLib/Raster.h b/GeoLib/Raster.h
index 78406fedc7be37520d1b57e8acd24990d87e70bf..1f4473ce0ef8754058b2d548c42314eccd73e5c1 100644
--- a/GeoLib/Raster.h
+++ b/GeoLib/Raster.h
@@ -99,7 +99,13 @@ public:
/**
* Returns the raster value at the position of the given point.
*/
- double getValueAtPoint(const GeoLib::Point &pnt);
+ double getValueAtPoint(const GeoLib::Point &pnt) const;
+
+ /// Interpolates the elevation of the given point based on the 8-neighbourhood of the raster cell it is located on
+ double interpolateValueAtPoint(const GeoLib::Point &pnt) const;
+
+ /// Checks if the given point is located within the (x,y)-extension of the raster.
+ bool isPntOnRaster(GeoLib::Point const& node) const;
~Raster();
diff --git a/Gui/DataView/GeoMapper.cpp b/Gui/DataView/GeoMapper.cpp
index 658b01cea0e3a9d0253d238563ef4b3aa17e5ded..923838f44259e88549815d5e292c35bfc807e2d7 100644
--- a/Gui/DataView/GeoMapper.cpp
+++ b/Gui/DataView/GeoMapper.cpp
@@ -106,7 +106,7 @@ void GeoMapper::mapData()
{
GeoLib::Point* pnt ((*points)[j]);
(*pnt)[2] = (_grid) ? this->getMeshElevation((*pnt)[0],(*pnt)[1], min_val, max_val)
- : this->getDemElevation((*pnt)[0],(*pnt)[1]);
+ : this->getDemElevation(*pnt);
}
}
else
@@ -115,7 +115,7 @@ void GeoMapper::mapData()
{
GeoLib::Point* pnt ((*points)[j]);
double offset = (_grid) ? (this->getMeshElevation((*pnt)[0],(*pnt)[1], min_val, max_val) - (*pnt)[2])
- : (this->getDemElevation((*pnt)[0],(*pnt)[1]) - (*pnt)[2]);
+ : this->getDemElevation(*pnt);
GeoLib::StationBorehole* borehole = static_cast<GeoLib::StationBorehole*>(pnt);
const std::vector<GeoLib::Point*> layers = borehole->getProfile();
@@ -129,21 +129,12 @@ void GeoMapper::mapData()
}
}
-float GeoMapper::getDemElevation(double x, double y) const
+float GeoMapper::getDemElevation(GeoLib::Point const& pnt) const
{
- const double origin_x(_raster->getOrigin()[0]);
- const double origin_y(_raster->getOrigin()[1]);
- const double cellsize(_raster->getRasterPixelSize());
- const std::size_t width(_raster->getNCols());
- const std::size_t height(_raster->getNRows());
-
- if ((x<origin_x) || (x>origin_x+(width*cellsize)) || (y<origin_y) || (y>origin_y+(height*cellsize)))
- return 0;
-
- const unsigned x_index = static_cast<unsigned>((x-origin_x)/cellsize);
- const unsigned y_index = static_cast<unsigned>((y-origin_y)/cellsize);
-
- return static_cast<float>(*(_raster->begin() + (y_index*width+x_index)));
+ double const elevation (_raster->getValueAtPoint(pnt));
+ if (elevation-_raster->getNoDataValue() < std::numeric_limits<double>::epsilon())
+ return 0.0;
+ return static_cast<float>(elevation);
}
double GeoMapper::getMeshElevation(double x, double y, double min_val, double max_val) const
diff --git a/Gui/DataView/GeoMapper.h b/Gui/DataView/GeoMapper.h
index 6ae5b11e22def013105409d602458a3920e35319..db3160d1b6dc3739a2f5804c06dae904756d2b50 100644
--- a/Gui/DataView/GeoMapper.h
+++ b/Gui/DataView/GeoMapper.h
@@ -57,7 +57,7 @@ private:
double getMeshElevation(double x, double y, double min_val, double max_val) const;
// Returns the elevation at Point (x,y) based on a raster
- float getDemElevation(double x, double y) const;
+ float getDemElevation(GeoLib::Point const& pnt) const;
// Calculates the intersection of two lines embedded in the xy-plane
GeoLib::Point* calcIntersection(GeoLib::Point const*const p1, GeoLib::Point const*const p2, GeoLib::Point const*const q1, GeoLib::Point const*const q2) const;
diff --git a/MeshLib/MeshGenerators/MeshLayerMapper.cpp b/MeshLib/MeshGenerators/MeshLayerMapper.cpp
index 649acf81f714a914f0bbb6be03db3663511f4a09..38ea970050241e62ac84d8fc2cb1bd06897bf8c1 100644
--- a/MeshLib/MeshGenerators/MeshLayerMapper.cpp
+++ b/MeshLib/MeshGenerators/MeshLayerMapper.cpp
@@ -34,7 +34,6 @@
#include "Elements/Prism.h"
#include "MeshSurfaceExtraction.h"
-
#include "MathTools.h"
MeshLib::Mesh* MeshLayerMapper::createStaticLayers(MeshLib::Mesh const& mesh, std::vector<float> const& layer_thickness_vector, std::string const& mesh_name) const
@@ -132,12 +131,10 @@ bool MeshLayerMapper::layerMapping(MeshLib::Mesh &new_mesh, const GeoLib::Raster
const double y0(raster.getOrigin()[1]);
const double delta(raster.getRasterPixelSize());
const double no_data(raster.getNoDataValue());
- const std::size_t width(raster.getNCols());
- const std::size_t height(raster.getNRows());
double const*const elevation(raster.begin());
- const std::pair<double, double> xDim(x0, x0 + width * delta); // extension in x-dimension
- const std::pair<double, double> yDim(y0, y0 + height * delta); // extension in y-dimension
+ const std::pair<double, double> xDim(x0, x0 + raster.getNCols() * delta); // extension in x-dimension
+ const std::pair<double, double> yDim(y0, y0 + raster.getNRows() * delta); // extension in y-dimension
const size_t nNodes (new_mesh.getNNodes());
const size_t nNodesPerLayer (nNodes / (nLayers+1));
@@ -149,76 +146,22 @@ bool MeshLayerMapper::layerMapping(MeshLib::Mesh &new_mesh, const GeoLib::Raster
const std::vector<MeshLib::Node*> &nodes = new_mesh.getNodes();
for (unsigned i = firstNode; i < lastNode; ++i)
{
- const double* coords (nodes[i]->getCoords());
-
- if (!isNodeOnRaster(*nodes[i], xDim, yDim))
+ if (!raster.isPntOnRaster(*nodes[i]))
{
// use either default value or elevation from layer above
- nodes[i]->updateCoordinates(coords[0], coords[1], noDataReplacementValue);
+ nodes[i]->updateCoordinates((*nodes[i])[0], (*nodes[i])[1], noDataReplacementValue);
continue;
}
- // position in raster
- const double xPos ((coords[0] - xDim.first) / delta);
- const double yPos ((coords[1] - yDim.first) / delta);
- // raster cell index
- const size_t xIdx (static_cast<size_t>(floor(xPos)));
- const size_t yIdx (static_cast<size_t>(floor(yPos)));
-
- // weights for bilinear interpolation
- const double xShift = fabs(xPos-(xIdx+half_delta))/delta;
- const double yShift = fabs(yPos-(yIdx+half_delta))/delta;
- std::array<double,4> weight = {{ (1-xShift)*(1-xShift), xShift*(1-yShift), xShift*yShift, (1-xShift)*yShift }};
-
- // neightbors to include in interpolation
- const int xShiftIdx = (xPos-xIdx-half_delta>=0) ? 1 : -1;
- const int yShiftIdx = (yPos-yIdx-half_delta>=0) ? 1 : -1;
- const std::array<int,4> x_nb = {{ 0, xShiftIdx, xShiftIdx, 0 }};
- const std::array<int,4> 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] = elevation[(yIdx + y_nb[j])*width + (xIdx + x_nb[j])];
- if (fabs(pix_val[j] - no_data) < 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
- {
- nodes[i]->updateCoordinates(coords[0], coords[1], noDataReplacementValue);
- continue;
- }
- const double norm = (double)(4)/(4-no_data_count);
- std::for_each(weight.begin(), weight.end(), [&norm](double &val){val*=norm;});
- }
-
- // new value
- double z = MathLib::scalarProduct<double,4>(weight.data(), pix_val.data());
- nodes[i]->updateCoordinates(coords[0], coords[1], z);
+ double elevation (raster.interpolateValueAtPoint(*nodes[i]));
+ if (elevation - no_data < std::numeric_limits<double>::epsilon())
+ elevation = noDataReplacementValue;
+ nodes[i]->updateCoordinates((*nodes[i])[0], (*nodes[i])[1], elevation);
}
return true;
}
-bool MeshLayerMapper::isNodeOnRaster(const MeshLib::Node &node,
- const std::pair<double, double> &xDim,
- const std::pair<double, double> &yDim) const
-{
- if (node[0] < xDim.first || node[0] > xDim.second || node[1] < yDim.first || node[1] > yDim.second)
- return false;
-
- return true;
-}
-
MeshLib::Mesh* MeshLayerMapper::blendLayersWithSurface(MeshLib::Mesh &mesh, const unsigned nLayers, const std::string &dem_raster) const
{
// construct surface mesh from DEM
diff --git a/MeshLib/MeshGenerators/MeshLayerMapper.h b/MeshLib/MeshGenerators/MeshLayerMapper.h
index 543896d0650edcc192407493acb25b5adcd93c72..aab59cb115cc497e47a99f010a4806e933d242ec 100644
--- a/MeshLib/MeshGenerators/MeshLayerMapper.h
+++ b/MeshLib/MeshGenerators/MeshLayerMapper.h
@@ -67,9 +67,7 @@ public:
MeshLib::Mesh* blendLayersWithSurface(MeshLib::Mesh &mesh, const unsigned nLayers, const std::string &dem_raster) const;
private:
- /// Checks if the given mesh is within the dimensions given by xDim and yDim.
- bool isNodeOnRaster(const MeshLib::Node &node,
- const std::pair<double, double> &xDim, const std::pair<double, double> &yDim) const;
+
};
#endif //MESHLAYERMAPPER_H