From 582134feb311ab15513eb6879bc702bd1a2e6deb Mon Sep 17 00:00:00 2001
From: Karsten Rink <karsten.rink@ufz.de>
Date: Tue, 8 Apr 2014 15:52:34 +0200
Subject: [PATCH] optimised LayerMapping() method
---
Gui/DataView/MeshLayerEditDialog.cpp | 4 +-
Gui/DataView/MshLayerMapper.cpp | 165 +++++++++++++--------------
Gui/DataView/MshLayerMapper.h | 2 +-
MathLib/MathTools.cpp | 9 --
MathLib/MathTools.h | 2 -
MeshLib/Node.h | 2 +-
6 files changed, 83 insertions(+), 101 deletions(-)
diff --git a/Gui/DataView/MeshLayerEditDialog.cpp b/Gui/DataView/MeshLayerEditDialog.cpp
index 94ac5605218..3240b9bcdab 100644
--- a/Gui/DataView/MeshLayerEditDialog.cpp
+++ b/Gui/DataView/MeshLayerEditDialog.cpp
@@ -191,7 +191,7 @@ void MeshLayerEditDialog::accept()
const std::string imgPath ( this->_edits[0]->text().toStdString() );
const double noDataReplacementValue = this->_noDataReplacementEdit->text().toDouble();
if (!imgPath.empty())
- result = MshLayerMapper::LayerMapping(new_mesh, imgPath, nLayers, 0, noDataReplacementValue);
+ result = MshLayerMapper::LayerMapping(*new_mesh, imgPath, nLayers, 0, noDataReplacementValue);
}
else
{
@@ -214,7 +214,7 @@ void MeshLayerEditDialog::accept()
const double noDataReplacement = (i==0) ? 0.0 : -9999.0;
if (!imgPath.empty())
{
- result = MshLayerMapper::LayerMapping(new_mesh, imgPath, nLayers, i, noDataReplacement);
+ result = MshLayerMapper::LayerMapping(*new_mesh, imgPath, nLayers, i, noDataReplacement);
if (result==0) break;
}
}
diff --git a/Gui/DataView/MshLayerMapper.cpp b/Gui/DataView/MshLayerMapper.cpp
index 547a7a16a5a..f218fa7482b 100644
--- a/Gui/DataView/MshLayerMapper.cpp
+++ b/Gui/DataView/MshLayerMapper.cpp
@@ -111,109 +111,102 @@ MeshLib::Mesh* MshLayerMapper::CreateLayers(const MeshLib::Mesh &mesh, const std
return new MeshLib::Mesh("SubsurfaceMesh", new_nodes, new_elems);
}
-int MshLayerMapper::LayerMapping(MeshLib::Mesh* new_mesh, const std::string &rasterfile,
+int MshLayerMapper::LayerMapping(MeshLib::Mesh &new_mesh, const std::string &rasterfile,
const unsigned nLayers, const unsigned layer_id, double noDataReplacementValue = 0.0)
{
- if (new_mesh == nullptr)
+ if (nLayers < layer_id)
{
- ERR("MshLayerMapper::LayerMapping() - Passed Mesh is NULL.");
+ ERR("MshLayerMapper::LayerMapping() - Mesh has only %d Layers, cannot assign layer %d.", nLayers, layer_id);
return 0;
}
- if (nLayers >= layer_id)
+ const GeoLib::Raster *raster(GeoLib::Raster::getRasterFromASCFile(rasterfile));
+ if (! raster) {
+ ERR("MshLayerMapper::LayerMapping - could not read raster file %s", rasterfile.c_str());
+ return 0;
+ }
+ const double x0(raster->getOrigin()[0]);
+ const double y0(raster->getOrigin()[1]);
+ const double delta(raster->getRasterPixelDistance());
+ 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 size_t nNodes (new_mesh.getNNodes());
+ const size_t nNodesPerLayer (nNodes / (nLayers+1));
+
+ const size_t firstNode (layer_id * nNodesPerLayer);
+ const size_t lastNode (firstNode + nNodesPerLayer);
+
+ const double half_delta = 0.5*delta;
+ const std::vector<MeshLib::Node*> &nodes = new_mesh.getNodes();
+ for (unsigned i = firstNode; i < lastNode; ++i)
{
- GeoLib::Raster *raster(GeoLib::Raster::getRasterFromASCFile(rasterfile));
- if (! raster) {
- ERR("MshLayerMapper::LayerMapping - could not read raster file %s", rasterfile.c_str());
- return 0;
- }
- const double x0(raster->getOrigin()[0]);
- const double y0(raster->getOrigin()[1]);
- const double delta(raster->getRasterPixelDistance());
- 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 double* coords (nodes[i]->getCoords());
- 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 size_t nNodes = new_mesh->getNNodes();
- const size_t nNodesPerLayer = nNodes / (nLayers+1);
-
- const size_t firstNode = layer_id * nNodesPerLayer;
- const size_t lastNode = firstNode + nNodesPerLayer;
-
- const double half_delta = 0.5*delta;
- const std::vector<MeshLib::Node*> nodes = new_mesh->getNodes();
- for (unsigned i = firstNode; i < lastNode; ++i)
+ if (!isNodeOnRaster(*nodes[i], xDim, yDim))
{
- const double* coords = nodes[i]->getCoords();
+ // use either default value or elevation from layer above
+ const double new_elevation = (layer_id == 0) ? noDataReplacementValue : (*nodes[i-nNodesPerLayer])[2];
+ nodes[i]->updateCoordinates(coords[0], coords[1], noDataReplacementValue);
+ continue;
+ }
- if (!isNodeOnRaster(*nodes[i], xDim, yDim))
+ // 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())
{
- if (layer_id == 0) // use default value
- nodes[i]->updateCoordinates(coords[0], coords[1], noDataReplacementValue);
- else // use z-value from layer above
- nodes[i]->updateCoordinates(coords[0], coords[1], (*nodes[i-nNodesPerLayer])[2]);
- continue;
+ weight[j] = 0;
+ no_data_count++;
}
+ }
- // 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)));
-
- // deviation of mesh node from centre of raster cell ( in [-1:1) because it is normalised by delta/2 )
- const double xShift = (xPos-xIdx-half_delta)/half_delta;
- const double yShift = (yPos-yIdx-half_delta)/half_delta;
-
- const int xShiftIdx = static_cast<int>((xShift>=0) ? ceil(xShift) : floor(xShift));
- const int yShiftIdx = static_cast<int>((yShift>=0) ? ceil(yShift) : floor(yShift));
-
- // determining the neighbouring pixels that add weight to the interpolation
- const int x_nb[4] = {0, xShiftIdx, xShiftIdx, 0};
- const int y_nb[4] = {0, 0, yShiftIdx, yShiftIdx};
-
- double locZ[4];
- locZ[0] = elevation[yIdx*width + xIdx];
- if (fabs(locZ[0] - no_data) > std::numeric_limits<double>::epsilon())
- {
- for (unsigned j=1; j<4; ++j)
- {
- locZ[j] = elevation[(yIdx+y_nb[j])*width + (xIdx+x_nb[j])];
- if (fabs(locZ[j] - no_data) < std::numeric_limits<double>::epsilon())
- locZ[j]=locZ[0];
- }
-
- double ome[4];
- double xi = 1-fabs(xShift);
- double eta = 1-fabs(xShift);
- MathLib::MPhi2D(ome, xi, eta);
-
- double z(0.0);
- for(unsigned j = 0; j < 4; ++j)
- z += ome[j] * locZ[j];
-
- nodes[i]->updateCoordinates(coords[0], coords[1], z);
- }
- else
+ // 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
{
- if (layer_id == 0) // use default value
- nodes[i]->updateCoordinates(coords[0], coords[1], noDataReplacementValue);
- else // use z-value from layer above
- nodes[i]->updateCoordinates(coords[0], coords[1], (*nodes[i-nNodesPerLayer])[2]);
+ const double new_elevation = (layer_id == 0) ? noDataReplacementValue : (*nodes[i-nNodesPerLayer])[2];
+ nodes[i]->updateCoordinates(coords[0], coords[1], noDataReplacementValue);
+ continue;
}
+ const double norm = 4/(4-no_data_count);
+ for_each(weight.begin(), weight.end(), [&norm](double &val){val*=norm;});
}
- delete raster;
- return 1;
+ // new value
+ double z = MathLib::scalarProduct<double,4>(weight.data(), pix_val.data());
+ nodes[i]->updateCoordinates(coords[0], coords[1], z);
}
- else
- ERR("MshLayerMapper::LayerMapping() - Mesh has only %d Layers, cannot assign layer %d.", nLayers, layer_id);
- return 0;
+
+ delete raster;
+ return 1;
}
bool MshLayerMapper::isNodeOnRaster(const MeshLib::Node &node,
@@ -231,7 +224,7 @@ MeshLib::Mesh* MshLayerMapper::blendLayersWithSurface(MeshLib::Mesh* mesh, const
// construct surface mesh from DEM
const MathLib::Vector3 dir(0,0,1);
MeshLib::Mesh* dem = MeshLib::MeshSurfaceExtraction::getMeshSurface(*mesh, dir);
- MshLayerMapper::LayerMapping(dem, dem_raster, 0, 0);
+ MshLayerMapper::LayerMapping(*dem, dem_raster, 0, 0);
const std::vector<MeshLib::Node*> dem_nodes (dem->getNodes());
const std::vector<MeshLib::Node*> mdl_nodes (mesh->getNodes());
diff --git a/Gui/DataView/MshLayerMapper.h b/Gui/DataView/MshLayerMapper.h
index daacaf2ae5e..d108e41cdda 100644
--- a/Gui/DataView/MshLayerMapper.h
+++ b/Gui/DataView/MshLayerMapper.h
@@ -46,7 +46,7 @@ public:
* Maps the z-values of nodes in the designated layer of the given mesh according to the given raster.
* Note: This only results in a valid mesh if the layers don't intersect each other.
*/
- static int LayerMapping(MeshLib::Mesh* msh, const std::string &rasterfile,
+ static int LayerMapping(MeshLib::Mesh &mesh, const std::string &rasterfile,
const unsigned nLayers, const unsigned layer_id, double noDataReplacementValue);
/**
diff --git a/MathLib/MathTools.cpp b/MathLib/MathTools.cpp
index 7de70163cbe..651fa75369b 100644
--- a/MathLib/MathTools.cpp
+++ b/MathLib/MathTools.cpp
@@ -86,13 +86,4 @@ double calcTetrahedronVolume(const double* x1, const double* x2, const double* x
+ (x1[2] - x4[2]) * ((x2[0] - x4[0]) * (x3[1] - x4[1]) - (x2[1] - x4[1]) * (x3[0] - x4[0]))) / 6.0;
}
-void MPhi2D(double* vf, double r, double s)
-{
- vf[0] = (1.0 + r) * (1.0 + s);
- vf[1] = (1.0 - r) * (1.0 + s);
- vf[2] = (1.0 - r) * (1.0 - s);
- vf[3] = (1.0 + r) * (1.0 - s);
- for (unsigned i = 0; i < 4; i++)
- vf[i] *= 0.25;
-}
} // namespace
diff --git a/MathLib/MathTools.h b/MathLib/MathTools.h
index f205d25b598..585d141b14c 100644
--- a/MathLib/MathTools.h
+++ b/MathLib/MathTools.h
@@ -177,8 +177,6 @@ T fastpow (T base, std::size_t exp)
return result;
}
-/// 2D linear interpolation function (TODO adopted from geo_mathlib)
-void MPhi2D(double* vf, double r, double s);
} // namespace
#endif /* MATHTOOLS_H_ */
diff --git a/MeshLib/Node.h b/MeshLib/Node.h
index e6000a6cd5d..2669445a8cc 100644
--- a/MeshLib/Node.h
+++ b/MeshLib/Node.h
@@ -39,7 +39,7 @@ class Node : public GeoLib::PointWithID
{
/* friend functions: */
#ifdef OGS_BUILD_GUI
- friend int MshLayerMapper::LayerMapping(MeshLib::Mesh* msh, const std::string &rasterfile, const unsigned nLayers,
+ friend int MshLayerMapper::LayerMapping(MeshLib::Mesh &msh, const std::string &rasterfile, const unsigned nLayers,
const unsigned layer_id, double noDataReplacementValue);
friend MeshLib::Mesh* MshLayerMapper::blendLayersWithSurface(MeshLib::Mesh* mesh, const unsigned nLayers, const std::string &dem_raster);
#endif
--
GitLab