diff --git a/BaseLib/DateTools.cpp b/BaseLib/DateTools.cpp
index 6f2b6e19c4bb7f1cf95e3f6a4dbf88a419bed906..157f572672df301afe55d8d6d9cced7e9c1f6550 100644
--- a/BaseLib/DateTools.cpp
+++ b/BaseLib/DateTools.cpp
@@ -44,8 +44,8 @@ std::string int2date(int date)
{
if (date > 10000000 && date < 22000000)
{
- int y = static_cast<int>(floor(date / 10000.0));
- int m = static_cast<int>(floor((date - (y * 10000)) / 100.0));
+ int y = static_cast<int>(std::floor(date / 10000.0));
+ int m = static_cast<int>(std::floor((date - (y * 10000)) / 100.0));
int d = date - (y * 10000) - (m * 100);
std::stringstream ss;
if (d < 10)
@@ -68,9 +68,9 @@ std::string date2string(double ddate)
}
int rest (static_cast<int>(ddate));
- int y = static_cast<int>(floor(rest / 10000.0));
+ int y = static_cast<int>(std::floor(rest / 10000.0));
rest = rest % (y * 10000);
- int m = static_cast<int>(floor(rest / 100.0));
+ int m = static_cast<int>(std::floor(rest / 100.0));
if (m < 1 || m > 12)
WARN("date2String(): month not in [1:12].");
rest = rest % (m * 100);
diff --git a/FileIO/AsciiRasterInterface.cpp b/FileIO/AsciiRasterInterface.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..a7dca0c835b6209578d1c555ef6b3dc8bdc96feb
--- /dev/null
+++ b/FileIO/AsciiRasterInterface.cpp
@@ -0,0 +1,221 @@
+/**
+ * @file AsciiRasterInterface.cpp
+ * @author Karsten Rink
+ * @date 2014-09-10
+ * @brief Implementation of the AsciiRasterInterface class.
+ *
+ * @copyright
+ * Copyright (c) 2012-2014, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ */
+
+#include "AsciiRasterInterface.h"
+
+// ThirdParty/logog
+#include "logog/include/logog.hpp"
+
+#include "Raster.h"
+
+// BaseLib
+#include "FileTools.h"
+#include "StringTools.h"
+
+namespace FileIO {
+
+GeoLib::Raster* AsciiRasterInterface::readRaster(std::string const& fname)
+{
+ std::string ext (BaseLib::getFileExtension(fname));
+ std::transform(ext.begin(), ext.end(), ext.begin(), tolower);
+ if (ext.compare("asc") == 0)
+ return getRasterFromASCFile(fname);
+ if (ext.compare("grd") == 0)
+ return getRasterFromSurferFile(fname);
+ return nullptr;
+}
+
+GeoLib::Raster* AsciiRasterInterface::getRasterFromASCFile(std::string const& fname)
+{
+ std::ifstream in(fname.c_str());
+
+ if (!in.is_open()) {
+ WARN("Raster::getRasterFromASCFile(): Could not open file %s.", fname.c_str());
+ return nullptr;
+ }
+
+ // header information
+ std::size_t n_cols(0), n_rows(0);
+ double xllcorner(0.0), yllcorner(0.0), cell_size(0.0), no_data_val(-9999);
+
+ if (readASCHeader(in, n_cols, n_rows, xllcorner, yllcorner, cell_size, no_data_val)) {
+ double* values = new double[n_cols*n_rows];
+ std::string s;
+ // read the data into the double-array
+ for (std::size_t j(0); j < n_rows; ++j) {
+ const size_t idx ((n_rows - j - 1) * n_cols);
+ for (size_t i(0); i < n_cols; ++i) {
+ in >> s;
+ values[idx+i] = strtod(BaseLib::replaceString(",", ".", s).c_str(),0);
+
+ }
+ }
+ in.close();
+ GeoLib::Raster *raster(new GeoLib::Raster(n_cols, n_rows, xllcorner, yllcorner,
+ cell_size, values, values+n_cols*n_rows, no_data_val));
+ delete [] values;
+ return raster;
+ } else {
+ WARN("Raster::getRasterFromASCFile(): Could not read header of file %s", fname.c_str());
+ return nullptr;
+ }
+}
+
+bool AsciiRasterInterface::readASCHeader(std::ifstream &in, std::size_t &n_cols, std::size_t &n_rows,
+ double &xllcorner, double &yllcorner, double &cell_size, double &no_data_val)
+{
+ std::string tag, value;
+
+ in >> tag;
+ if (tag.compare("ncols") == 0) {
+ in >> value;
+ n_cols = atoi(value.c_str());
+ } else return false;
+
+ in >> tag;
+ if (tag.compare("nrows") == 0) {
+ in >> value;
+ n_rows = atoi(value.c_str());
+ } else return false;
+
+ in >> tag;
+ if (tag.compare("xllcorner") == 0) {
+ in >> value;
+ xllcorner = strtod(BaseLib::replaceString(",", ".", value).c_str(), 0);
+ } else return false;
+
+ in >> tag;
+ if (tag.compare("yllcorner") == 0) {
+ in >> value;
+ yllcorner = strtod(BaseLib::replaceString(",", ".", value).c_str(), 0);
+ } else return false;
+
+ in >> tag;
+ if (tag.compare("cellsize") == 0) {
+ in >> value;
+ cell_size = strtod(BaseLib::replaceString(",", ".", value).c_str(), 0);
+ } else return false;
+
+ in >> tag;
+ if (tag.compare("NODATA_value") == 0) {
+ in >> value;
+ no_data_val = strtod(BaseLib::replaceString(",", ".", value).c_str(), 0);
+ } else return false;
+
+ return true;
+}
+
+GeoLib::Raster* AsciiRasterInterface::getRasterFromSurferFile(std::string const& fname)
+{
+ std::ifstream in(fname.c_str());
+
+ if (!in.is_open()) {
+ ERR("Raster::getRasterFromSurferFile() - Could not open file %s", fname.c_str());
+ return nullptr;
+ }
+
+ // header information
+ std::size_t n_cols(0), n_rows(0);
+ double xllcorner(0.0), yllcorner(0.0), cell_size(0.0), min(0.0), max(0.0);
+
+ if (readSurferHeader(in, n_cols, n_rows, xllcorner, yllcorner, cell_size, min, max))
+ {
+ const double no_data_val (min-1);
+ double* values = new double[n_cols*n_rows];
+ std::string s;
+ // read the data into the double-array
+ for (size_t j(0); j < n_rows; ++j)
+ {
+ const size_t idx (j * n_cols);
+ for (size_t i(0); i < n_cols; ++i)
+ {
+ in >> s;
+ const double val (strtod(BaseLib::replaceString(",", ".", s).c_str(),0));
+ values[idx+i] = (val > max || val < min) ? no_data_val : val;
+ }
+ }
+ in.close();
+ GeoLib::Raster *raster(new GeoLib::Raster(n_cols, n_rows, xllcorner, yllcorner, cell_size,
+ values, values+n_cols*n_rows, no_data_val));
+ delete [] values;
+ return raster;
+ } else {
+ ERR("Raster::getRasterFromASCFile() - could not read header of file %s", fname.c_str());
+ return nullptr;
+ }
+}
+
+bool AsciiRasterInterface::readSurferHeader(std::ifstream &in, size_t &n_cols, std::size_t &n_rows,
+ double &xllcorner, double &yllcorner, double &cell_size, double &min, double &max)
+{
+ std::string tag;
+
+ in >> tag;
+
+ if (tag.compare("DSAA") != 0)
+ {
+ ERR("Error in readSurferHeader() - No Surfer file.");
+ return false;
+ }
+ else
+ {
+ in >> n_cols >> n_rows;
+ in >> min >> max;
+ xllcorner = min;
+ cell_size = (max-min)/static_cast<double>(n_cols);
+
+ in >> min >> max;
+ yllcorner = min;
+
+ if (ceil((max-min)/static_cast<double>(n_rows)) == ceil(cell_size))
+ cell_size = ceil(cell_size);
+ else
+ {
+ ERR("Error in readSurferHeader() - Anisotropic cellsize detected.");
+ return 0;
+ }
+ in >> min >> max;
+ }
+
+ return true;
+}
+
+void AsciiRasterInterface::writeRasterAsASC(GeoLib::Raster const& raster, std::string const& file_name)
+{
+ GeoLib::Point const& origin (raster.getOrigin());
+ unsigned const nCols (raster.getNCols());
+ unsigned const nRows (raster.getNRows());
+
+ // write header
+ std::ofstream out(file_name);
+ out << "ncols " << nCols << "\n";
+ out << "nrows " << nRows << "\n";
+ out << "xllcorner " << origin[0] << "\n";
+ out << "yllcorner " << origin[1] << "\n";
+ out << "cellsize " << raster.getRasterPixelSize() << "\n";
+ out << "NODATA_value " << raster.getNoDataValue() << "\n";
+
+ // write data
+ double const*const elevation(raster.begin());
+ for (unsigned row(0); row < nRows; ++row)
+ {
+ for (unsigned col(0); col < nCols; ++col)
+ {
+ out << elevation[(nRows-row-1) * nCols + col] << " ";
+ }
+ out << "\n";
+ }
+ out.close();
+}
+
+} // end namespace GeoLib
diff --git a/FileIO/AsciiRasterInterface.h b/FileIO/AsciiRasterInterface.h
new file mode 100644
index 0000000000000000000000000000000000000000..e7e654915700f5b695c8ac31762c3b75e5c1ef0c
--- /dev/null
+++ b/FileIO/AsciiRasterInterface.h
@@ -0,0 +1,57 @@
+/**
+ * @file AsciiRasterInterface.h
+ * @author Karsten Rink
+ * @date 2014-09-10
+ * @brief Definition of the AsciiRasterInterface class.
+ *
+ * @copyright
+ * Copyright (c) 2012-2014, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ */
+
+#ifndef ASCIIRASTERINTERFACE_H_
+#define ASCIIRASTERINTERFACE_H_
+
+#include <fstream>
+
+namespace GeoLib {
+ class Raster;
+}
+
+namespace FileIO {
+
+/**
+ * Interface for reading and writing a number of ASCII raster formats.
+ * Currently supported are reading and writing of Esri asc-files and
+ * reading of Surfer grd-files.
+ */
+class AsciiRasterInterface {
+public:
+ /// Reads raster file by detecting type based on extension and then calling the apropriate method
+ static GeoLib::Raster* readRaster(std::string const& fname);
+
+ /// Reads an ArcGis ASC raster file
+ static GeoLib::Raster* getRasterFromASCFile(std::string const& fname);
+
+ /// Reads a Surfer GRD raster file
+ static GeoLib::Raster* getRasterFromSurferFile(std::string const& fname);
+
+ /// Writes an Esri asc-file
+ static void writeRasterAsASC(GeoLib::Raster const& raster, std::string const& file_name);
+
+
+private:
+ /// Reads the header of a Esri asc-file.
+ static bool readASCHeader(std::ifstream &in, std::size_t &n_cols, std::size_t &n_rows,
+ double &xllcorner, double &yllcorner, double &cell_size, double &no_data_val);
+
+ /// Reads the header of a Surfer grd-file.
+ static bool readSurferHeader(std::ifstream &in, size_t &n_cols, std::size_t &n_rows,
+ double &xllcorner, double &yllcorner, double &cell_size, double &min, double &max);
+};
+
+}
+
+#endif /* ASCIIRASTERINTERFACE_H_ */
diff --git a/FileIO/CMakeLists.txt b/FileIO/CMakeLists.txt
index 1c09aec5d5e711fa03bf26654af1cf10510827ec..13fafee37191aa57ffeb9366dace79095f62e8ff 100644
--- a/FileIO/CMakeLists.txt
+++ b/FileIO/CMakeLists.txt
@@ -1,6 +1,8 @@
# Source files
# GET_SOURCE_FILES(SOURCES_FILEIO)
SET( SOURCES
+ AsciiRasterInterface.h
+ AsciiRasterInterface.cpp
GMSInterface.h
GMSInterface.cpp
GMSHInterface.h
diff --git a/FileIO/TetGenInterface.cpp b/FileIO/TetGenInterface.cpp
index 1b4228f452c91727c9cc880569669a85e8d3fb88..1297eaff26ab1ce803e4f32a0ff53c2cc9bdb8f4 100644
--- a/FileIO/TetGenInterface.cpp
+++ b/FileIO/TetGenInterface.cpp
@@ -621,7 +621,7 @@ void TetGenInterface::write3dElements(std::ofstream &out,
// get position where number of facets need to be written and figure out worst case of chars that are needed
const std::streamoff before_elems_pos (out.tellp());
- const unsigned n_spaces (static_cast<unsigned>(floor(log(nElements*8))) + 1);
+ const unsigned n_spaces (static_cast<unsigned>(std::floor(log(nElements*8))) + 1);
out << std::string(n_spaces, ' ') << "\n";
unsigned element_count(0);
diff --git a/GeoLib/Raster.cpp b/GeoLib/Raster.cpp
index 3879f675a8b1fd585df78360e871abd7ccb371f2..603f5578bdb19b823cba07be6fc04fc1c950dbc6 100644
--- a/GeoLib/Raster.cpp
+++ b/GeoLib/Raster.cpp
@@ -73,8 +73,8 @@ Raster* Raster::getRasterFromSurface(Surface const& sfc, double cell_size, doubl
Point const& ll (sfc.getAABB().getMinPoint());
Point const& ur (sfc.getAABB().getMaxPoint());
- const std::size_t n_cols = static_cast<size_t>(fabs(ur[0]-ll[0]) / cell_size)+1;
- const std::size_t n_rows = static_cast<size_t>(fabs(ur[1]-ll[1]) / cell_size)+1;
+ const std::size_t n_cols = static_cast<size_t>(std::fabs(ur[0]-ll[0]) / cell_size)+1;
+ const std::size_t n_rows = static_cast<size_t>(std::fabs(ur[1]-ll[1]) / cell_size)+1;
const size_t n_triangles (sfc.getNTriangles());
double *z_vals (new double[n_cols*n_rows]);
size_t k(0);
@@ -101,13 +101,13 @@ 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)))
{
- int cell_x = static_cast<int>(floor((pnt[0] - _ll_pnt[0])/_cell_size));
- int cell_y = static_cast<int>(floor((pnt[1] - _ll_pnt[1])/_cell_size));
+ int cell_x = static_cast<int>(std::floor((pnt[0] - _ll_pnt[0])/_cell_size));
+ int cell_y = static_cast<int>(std::floor((pnt[1] - _ll_pnt[1])/_cell_size));
// use raster boundary values if node is outside raster due to rounding errors or floating point arithmetic
cell_x = (cell_x < 0) ? 0 : ((cell_x > static_cast<int>(_n_cols)) ? (_n_cols-1) : cell_x);
@@ -119,189 +119,60 @@ double Raster::getValueAtPoint(const GeoLib::Point &pnt)
return _no_data_val;
}
-void Raster::writeRasterAsASC(std::ostream &os) const
+double Raster::interpolateValueAtPoint(GeoLib::Point const& pnt) const
{
- // write header
- os << "ncols " << _n_cols << "\n";
- os << "nrows " << _n_rows << "\n";
- os << "xllcorner " << _ll_pnt[0] << "\n";
- os << "yllcorner " << _ll_pnt[1] << "\n";
- os << "cellsize " << _cell_size << "\n";
- os << "NODATA_value " << _no_data_val << "\n";
-
- // write data
- for (unsigned row(0); row<_n_rows; row++) {
- for (unsigned col(0); col<_n_cols; col++) {
- os << _raster_data[(_n_rows-row-1)*_n_cols+col] << " ";
- }
- os << "\n";
- }
-}
-
-Raster* Raster::readRaster(std::string const& fname)
-{
- std::string ext (BaseLib::getFileExtension(fname));
- std::transform(ext.begin(), ext.end(), ext.begin(), ::tolower);
- if (ext.compare("asc") == 0)
- return getRasterFromASCFile(fname);
- if (ext.compare("grd") == 0)
- return getRasterFromSurferFile(fname);
- return nullptr;
-}
-
-Raster* Raster::getRasterFromASCFile(std::string const& fname)
-{
- std::ifstream in(fname.c_str());
-
- if (!in.is_open()) {
- WARN("Raster::getRasterFromASCFile(): Could not open file %s.", fname.c_str());
- return NULL;
- }
-
- // header information
- std::size_t n_cols(0), n_rows(0);
- double xllcorner(0.0), yllcorner(0.0), cell_size(0.0), no_data_val(-9999);
-
- if (readASCHeader(in, n_cols, n_rows, xllcorner, yllcorner, cell_size, no_data_val)) {
- double* values = new double[n_cols*n_rows];
- std::string s;
- // read the data into the double-array
- for (size_t j(0); j < n_rows; ++j) {
- const size_t idx ((n_rows - j - 1) * n_cols);
- for (size_t i(0); i < n_cols; ++i) {
- in >> s;
- values[idx+i] = strtod(BaseLib::replaceString(",", ".", s).c_str(),0);
-
- }
- }
- in.close();
- Raster *raster(new Raster(n_cols, n_rows, xllcorner, yllcorner,
- cell_size, values, values+n_cols*n_rows, no_data_val));
- delete [] values;
- return raster;
- } else {
- WARN("Raster::getRasterFromASCFile(): Could not read header of file %s", fname.c_str());
- return NULL;
- }
+ // 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());
}
-bool Raster::readASCHeader(std::ifstream &in, std::size_t &n_cols, std::size_t &n_rows,
- double &xllcorner, double &yllcorner, double &cell_size, double &no_data_val)
+bool Raster::isPntOnRaster(GeoLib::Point const& pnt) const
{
- std::string tag, value;
-
- in >> tag;
- if (tag.compare("ncols") == 0) {
- in >> value;
- n_cols = atoi(value.c_str());
- } else return false;
-
- in >> tag;
- if (tag.compare("nrows") == 0) {
- in >> value;
- n_rows = atoi(value.c_str());
- } else return false;
-
- in >> tag;
- if (tag.compare("xllcorner") == 0) {
- in >> value;
- xllcorner = strtod(BaseLib::replaceString(",", ".", value).c_str(), 0);
- } else return false;
-
- in >> tag;
- if (tag.compare("yllcorner") == 0) {
- in >> value;
- yllcorner = strtod(BaseLib::replaceString(",", ".", value).c_str(), 0);
- } else return false;
-
- in >> tag;
- if (tag.compare("cellsize") == 0) {
- in >> value;
- cell_size = strtod(BaseLib::replaceString(",", ".", value).c_str(), 0);
- } else return false;
-
- in >> tag;
- if (tag.compare("NODATA_value") == 0) {
- in >> value;
- no_data_val = strtod(BaseLib::replaceString(",", ".", value).c_str(), 0);
- } else return false;
-
- return true;
-}
-
-Raster* Raster::getRasterFromSurferFile(std::string const& fname)
-{
- std::ifstream in(fname.c_str());
-
- if (!in.is_open()) {
- ERR("Raster::getRasterFromSurferFile() - Could not open file %s", fname.c_str());
- return NULL;
- }
-
- // header information
- std::size_t n_cols(0), n_rows(0);
- double xllcorner(0.0), yllcorner(0.0), cell_size(0.0), min(0.0), max(0.0);
-
- if (readSurferHeader(in, n_cols, n_rows, xllcorner, yllcorner, cell_size, min, max))
- {
- const double no_data_val (min-1);
- double* values = new double[n_cols*n_rows];
- std::string s;
- // read the data into the double-array
- for (size_t j(0); j < n_rows; ++j)
- {
- const size_t idx (j * n_cols);
- for (size_t i(0); i < n_cols; ++i)
- {
- in >> s;
- const double val (strtod(BaseLib::replaceString(",", ".", s).c_str(),0));
- values[idx+i] = (val > max || val < min) ? no_data_val : val;
- }
- }
- in.close();
- Raster *raster(new Raster(n_cols, n_rows, xllcorner, yllcorner,
- cell_size, values, values+n_cols*n_rows, no_data_val));
- delete [] values;
- return raster;
- } else {
- ERR("Raster::getRasterFromASCFile() - could not read header of file %s", fname.c_str());
- return NULL;
- }
-}
-
-bool Raster::readSurferHeader(std::ifstream &in, size_t &n_cols, std::size_t &n_rows,
- double &xllcorner, double &yllcorner, double &cell_size, double &min, double &max)
-{
- std::string tag;
-
- in >> tag;
-
- if (tag.compare("DSAA") != 0)
- {
- ERR("Error in readSurferHeader() - No Surfer file.");
+ 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;
- }
- else
- {
- in >> n_cols >> n_rows;
- in >> min >> max;
- xllcorner = min;
- cell_size = (max-min)/static_cast<double>(n_cols);
-
- in >> min >> max;
- yllcorner = min;
-
- if (ceil((max-min)/static_cast<double>(n_rows)) == ceil(cell_size))
- cell_size = ceil(cell_size);
- else
- {
- ERR("Error in readSurferHeader() - Anisotropic cellsize detected.");
- return 0;
- }
- in >> min >> max;
- }
-
- return true;
+ return true;
}
} // end namespace GeoLib
diff --git a/GeoLib/Raster.h b/GeoLib/Raster.h
index 52a0c74286d6c2eb14541e9b820d20e65ce62ffe..1f4473ce0ef8754058b2d548c42314eccd73e5c1 100644
--- a/GeoLib/Raster.h
+++ b/GeoLib/Raster.h
@@ -70,7 +70,7 @@ public:
/**
* get the distance between raster pixels
*/
- double getRasterPixelDistance() const { return _cell_size; }
+ double getRasterPixelSize() const { return _cell_size; }
/**
* get the origin of lower left raster cell
@@ -99,39 +99,20 @@ 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;
- ~Raster();
+ /// 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;
- /**
- * Write meta data and raw raster data as asci file into the output stream.
- * @param os the output stream
- */
- void writeRasterAsASC(std::ostream &os) const;
-
- /// Reads raster file by detecting type based on extension and then calling the apropriate method
- static Raster* readRaster(std::string const& fname);
+ ~Raster();
/// Creates a Raster based on a GeoLib::Surface
static Raster* getRasterFromSurface(Surface const& sfc, double cell_size, double no_data_val = -9999);
- /// Reads an ArcGis ASC raster file
- static Raster* getRasterFromASCFile(std::string const& fname);
-
- /// Reads a Surfer GRD raster file
- static Raster* getRasterFromSurferFile(std::string const& fname);
private:
- static bool readASCHeader(std::ifstream &in, std::size_t &n_cols, std::size_t &n_rows,
- double &xllcorner, double &yllcorner, double &cell_size, double &no_data_val);
- /**
- * Reads the header of a Surfer grd-file.
- * \param header The ascHeader-object into which all the information will be written.
- * \param in FileInputStream used for reading the data.
- * \return True if the header could be read correctly, false otherwise.
- */
- static bool readSurferHeader(std::ifstream &in, size_t &n_cols, std::size_t &n_rows,
- double &xllcorner, double &yllcorner, double &cell_size, double &min, double &max);
-
void setCellSize(double cell_size);
void setNoDataVal (double no_data_val);
diff --git a/Gui/DataView/DiagramView/DiagramScene.cpp b/Gui/DataView/DiagramView/DiagramScene.cpp
index f6f60a49dfd188fd0f831d986f1803aa610bc978..c4a49553444ff270546b3edb8f84265e64952072 100644
--- a/Gui/DataView/DiagramView/DiagramScene.cpp
+++ b/Gui/DataView/DiagramView/DiagramScene.cpp
@@ -134,15 +134,15 @@ void DiagramScene::adjustAxis(float &min, float &max, int &numberOfTicks)
{
const int MinTicks = 4;
double grossStep = (max - min) / MinTicks;
- double step = pow(10.0, floor(log10(grossStep)));
+ double step = pow(10.0, std::floor(log10(grossStep)));
if (5 * step < grossStep)
step *= 5;
else if (2 * step < grossStep)
step *= 2;
- numberOfTicks = int(ceil(max / step) - floor(min / step));
+ numberOfTicks = int(ceil(max / step) - std::floor(min / step));
if (numberOfTicks < MinTicks)
numberOfTicks = MinTicks;
- min = floor(min / step) * step;
+ min = std::floor(min / step) * step;
max = ceil(max / step) * step;
}
diff --git a/Gui/DataView/DirectConditionGenerator.cpp b/Gui/DataView/DirectConditionGenerator.cpp
index abb16bb120dda8dfa93d32d5650fbfd13cbce596..89d63df52262a81c0c7cae3fdc7735c86e1baa1e 100644
--- a/Gui/DataView/DirectConditionGenerator.cpp
+++ b/Gui/DataView/DirectConditionGenerator.cpp
@@ -19,6 +19,8 @@
#include "DirectConditionGenerator.h"
+#include "AsciiRasterInterface.h"
+
#include "Raster.h"
#include "MeshSurfaceExtraction.h"
#include "PointWithID.h"
@@ -32,7 +34,7 @@ const std::vector< std::pair<size_t,double> >& DirectConditionGenerator::directT
{
if (_direct_values.empty())
{
- GeoLib::Raster* raster (GeoLib::Raster::readRaster(filename));
+ GeoLib::Raster* raster (FileIO::AsciiRasterInterface::readRaster(filename));
if (! raster) {
ERR("Error in DirectConditionGenerator::directToSurfaceNodes() - could not load raster file.");
return _direct_values;
@@ -62,7 +64,7 @@ const std::vector< std::pair<size_t,double> >& DirectConditionGenerator::directW
{
if (_direct_values.empty())
{
- GeoLib::Raster* raster (GeoLib::Raster::readRaster(filename));
+ GeoLib::Raster* raster (FileIO::AsciiRasterInterface::readRaster(filename));
if (!raster) {
ERR("Error in DirectConditionGenerator::directWithSurfaceIntegration() - could not load raster file.");
return _direct_values;
diff --git a/Gui/DataView/GeoMapper.cpp b/Gui/DataView/GeoMapper.cpp
index 0e58afedbe991ef557ca00c2f20590e4854ed709..9541cb730bdcb41fe430bcc96d0507135253e0f8 100644
--- a/Gui/DataView/GeoMapper.cpp
+++ b/Gui/DataView/GeoMapper.cpp
@@ -19,6 +19,8 @@
#include <numeric>
+#include "FileIO/AsciiRasterInterface.h"
+
#include "AABB.h"
#include "Mesh.h"
#include "Elements/Element.h"
@@ -44,7 +46,7 @@ GeoMapper::~GeoMapper()
void GeoMapper::mapOnDEM(const std::string &file_name)
{
- this->_raster = GeoLib::Raster::getRasterFromASCFile(file_name);
+ this->_raster = FileIO::AsciiRasterInterface::getRasterFromASCFile(file_name);
if (! _raster) {
ERR("GeoMapper::mapOnDEM(): failed to load %s", file_name.c_str());
return;
@@ -104,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
@@ -113,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();
@@ -127,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->getRasterPixelDistance());
- 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 (std::abs(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/Gui/DataView/MeshLayerEditDialog.cpp b/Gui/DataView/MeshLayerEditDialog.cpp
index 91650be0a9b5ed3a51506437c8efbf8211a75488..231275565f02e5cfb72321cca63a1a8057b9c506 100644
--- a/Gui/DataView/MeshLayerEditDialog.cpp
+++ b/Gui/DataView/MeshLayerEditDialog.cpp
@@ -106,12 +106,12 @@ void MeshLayerEditDialog::createWithRasters()
this->_layerBox = new QGroupBox(this);
this->_layerBox->setTitle(selectText);
- for (unsigned i = 0; i <= _n_layers+1; ++i)
+ for (unsigned i = 0; i <= _n_layers; ++i)
{
QString text("");
if (i==0) text = "Surface";
- else if (i>_n_layers) text = "Layer" + QString::number(_n_layers) + "-Bottom";
- else text="Layer" + QString::number(i) + "-Top";
+ else if (i == _n_layers) text = "Layer" + QString::number(_n_layers) + "-Bottom";
+ else text="Layer" + QString::number(i+1) + "-Top";
QLineEdit* edit (new QLineEdit(this));
QPushButton* button (new QPushButton("...", _layerBox));
@@ -178,38 +178,30 @@ void MeshLayerEditDialog::createMeshToolSelection()
MeshLib::Mesh* MeshLayerEditDialog::createPrismMesh()
{
- const unsigned nLayers = _layerEdit->text().toInt();
- std::vector<float> layer_thickness;
- for (unsigned i=0; i<nLayers; ++i)
- {
- // "100" is just a default size to have any value for extruding 2D elements.
- // The actual mapping based on a raster file will be performed later.
- const float thickness = (_use_rasters) ? 100 : (this->_edits[i]->text().toFloat());
- layer_thickness.push_back(thickness);
- }
-
- MeshLib::Mesh* new_mesh = MeshLayerMapper::CreateLayers(*_msh, layer_thickness);
-
- if (_use_rasters)
- {
- for (unsigned i=0; i<=nLayers; ++i)
- {
- const std::string imgPath ( this->_edits[i+1]->text().toStdString() );
- const double noDataReplacement = (i==0) ? 0.0 : -9999.0;
- if (!MeshLayerMapper::LayerMapping(*new_mesh, imgPath, nLayers, i, noDataReplacement))
- {
- delete new_mesh;
- return nullptr;
- }
- }
- if (this->_edits[0]->text().length()>0)
- {
- MeshLib::Mesh* final_mesh = MeshLayerMapper::blendLayersWithSurface(*new_mesh, nLayers, this->_edits[0]->text().toStdString());
- delete new_mesh;
- new_mesh = final_mesh;
- }
- }
- return new_mesh;
+ const unsigned nLayers = _layerEdit->text().toInt();
+ std::vector<float> layer_thickness;
+ for (unsigned i=0; i<nLayers; ++i)
+ {
+ // "100" is just a default size to have any value for extruding 2D elements.
+ // The actual mapping based on a raster file will be performed later.
+ const float thickness = (_use_rasters) ? 100 : (this->_edits[i]->text().toFloat());
+ layer_thickness.push_back(thickness);
+ }
+
+ MeshLayerMapper mapper;
+ MeshLib::Mesh* new_mesh (nullptr);
+
+ if (_use_rasters)
+ {
+ std::vector<std::string> raster_paths;
+ for (int i=nLayers; i>=0; --i)
+ raster_paths.push_back(this->_edits[i]->text().toStdString());
+ if (mapper.createLayers(*_msh, raster_paths))
+ new_mesh= mapper.getMesh("SubsurfaceMesh");
+ }
+ else
+ new_mesh = mapper.createStaticLayers(*_msh, layer_thickness);
+ return new_mesh;
}
MeshLib::Mesh* MeshLayerEditDialog::createTetMesh()
@@ -226,12 +218,12 @@ MeshLib::Mesh* MeshLayerEditDialog::createTetMesh()
if (_use_rasters)
{
std::vector<std::string> raster_paths(nLayers+1);
- for (unsigned i=0; i<=nLayers; ++i)
- raster_paths[i] = this->_edits[i+1]->text().toStdString();
+ for (int i=nLayers; i>=0; --i)
+ raster_paths[i] = this->_edits[i]->text().toStdString();
LayeredVolume lv;
- lv.createGeoVolumes(*_msh, raster_paths);
+ lv.createLayers(*_msh, raster_paths);
- tg_mesh = lv.getMesh();
+ tg_mesh = lv.getMesh("SubsurfaceMesh");
QString file_path("");
if (tg_mesh)
@@ -246,7 +238,8 @@ MeshLib::Mesh* MeshLayerEditDialog::createTetMesh()
std::vector<float> layer_thickness;
for (unsigned i=0; i<nLayers; ++i)
layer_thickness.push_back(this->_edits[i]->text().toFloat());
- tg_mesh = MeshLayerMapper::CreateLayers(*_msh, layer_thickness);
+ MeshLayerMapper const mapper;
+ tg_mesh = mapper.createStaticLayers(*_msh, layer_thickness);
std::vector<MeshLib::Node> tg_attr;
FileIO::TetGenInterface tetgen_interface;
tetgen_interface.writeTetGenSmesh(filename.toStdString(), *tg_mesh, tg_attr);
@@ -291,7 +284,8 @@ void MeshLayerEditDialog::accept()
new_mesh = new MeshLib::Mesh(*_msh);
const std::string imgPath ( this->_edits[0]->text().toStdString() );
const double noDataReplacementValue = this->_noDataReplacementEdit->text().toDouble();
- if (!MeshLayerMapper::LayerMapping(*new_mesh, imgPath, nLayers, 0, noDataReplacementValue))
+ MeshLayerMapper const mapper;
+ if (!mapper.layerMapping(*new_mesh, imgPath, noDataReplacementValue))
{
delete new_mesh;
return;
diff --git a/Gui/VtkVis/VtkColorLookupTable.cpp b/Gui/VtkVis/VtkColorLookupTable.cpp
index 93b25955192e4ac3b1a96e4b85eea34319302ace..425705500703a03c2e4720eb6b570756bd8ae773 100644
--- a/Gui/VtkVis/VtkColorLookupTable.cpp
+++ b/Gui/VtkVis/VtkColorLookupTable.cpp
@@ -70,7 +70,7 @@ void VtkColorLookupTable::Build()
for (std::map<double, unsigned char*>::const_iterator it = _dict.begin(); it != _dict.end(); ++it)
{
double val = (it->first < range[0]) ? range[0] : ((it->first > range[1]) ? range[1] : it->first);
- nextIndex = static_cast<size_t>( floor(val-range[0]) );
+ nextIndex = static_cast<size_t>( std::floor(val-range[0]) );
this->SetTableValue(nextIndex, it->second);
@@ -181,7 +181,7 @@ void VtkColorLookupTable::getColor(vtkIdType indx, unsigned char rgba[4]) const
? static_cast<vtkIdType>(this->TableRange[0])
: (indx >=this->TableRange[1] ? static_cast<vtkIdType>(this->TableRange[1]) - 1 : indx));
indx =
- static_cast<size_t>( floor( (indx - this->TableRange[0]) *
+ static_cast<size_t>( std::floor( (indx - this->TableRange[0]) *
(this->NumberOfColors / (this->TableRange[1] - this->TableRange[0])) ) );
unsigned char* _rgba;
diff --git a/Gui/VtkVis/VtkRaster.cpp b/Gui/VtkVis/VtkRaster.cpp
index 7cfa4db9ce9d204ffd55c5177ae2a9be60c85910..60a7fcbe2f178852c31c0109f89604534f51b9ce 100644
--- a/Gui/VtkVis/VtkRaster.cpp
+++ b/Gui/VtkVis/VtkRaster.cpp
@@ -27,6 +27,8 @@
#include "StringTools.h"
+#include "AsciiRasterInterface.h"
+
// GeoLib
#include "Raster.h"
@@ -54,16 +56,16 @@ vtkImageAlgorithm* VtkRaster::loadImage(const std::string &fileName,
GeoLib::Raster *raster(nullptr);
if (fileInfo.suffix().toLower() == "asc") {
- raster = GeoLib::Raster::getRasterFromASCFile(fileName);
+ raster = FileIO::AsciiRasterInterface::getRasterFromASCFile(fileName);
}
else if (fileInfo.suffix().toLower() == "grd")
{
- raster = GeoLib::Raster::getRasterFromSurferFile(fileName);
+ raster = FileIO::AsciiRasterInterface::getRasterFromSurferFile(fileName);
}
if (raster) {
x0 = raster->getOrigin()[0];
y0 = raster->getOrigin()[1];
- delta = raster->getRasterPixelDistance();
+ delta = raster->getRasterPixelSize();
double const*const data (raster->begin());
return VtkRaster::loadImageFromArray(data, x0, y0,
raster->getNCols(), raster->getNRows(), delta,
diff --git a/MeshLib/MeshGenerators/ConvertRasterToMesh.cpp b/MeshLib/MeshGenerators/ConvertRasterToMesh.cpp
index c4d98231c9df6b3e81380976ccc2378c954a0fed..5fcfbfab7e1665a3942933ca2a0e736cfde7551c 100644
--- a/MeshLib/MeshGenerators/ConvertRasterToMesh.cpp
+++ b/MeshLib/MeshGenerators/ConvertRasterToMesh.cpp
@@ -76,7 +76,7 @@ MeshLib::Mesh* ConvertRasterToMesh::constructMesh(const double* pix_vals, const
const size_t height = _raster.getNRows()+1;
const size_t width = _raster.getNCols()+1;
size_t node_idx_count(0);
- const double distance(_raster.getRasterPixelDistance());
+ const double distance(_raster.getRasterPixelSize());
const double x_offset(_raster.getOrigin()[0]); // - distance / 2.0);
const double y_offset(_raster.getOrigin()[1]); // - distance / 2.0);
diff --git a/MeshLib/MeshGenerators/LayeredMeshGenerator.cpp b/MeshLib/MeshGenerators/LayeredMeshGenerator.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..d16106a875b0711133a68f293881369d6c8df60b
--- /dev/null
+++ b/MeshLib/MeshGenerators/LayeredMeshGenerator.cpp
@@ -0,0 +1,82 @@
+/**
+ * \file LayeredMeshGenerator.cpp
+ * \author Karsten Rink
+ * \date 2014-09-18
+ * \brief Implementation of the SubsurfaceMapper class.
+ *
+ * \copyright
+ * Copyright (c) 2012-2014, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#include "LayeredMeshGenerator.h"
+
+#include <vector>
+#include <fstream>
+
+#include "FileIO/AsciiRasterInterface.h"
+
+#include "Raster.h"
+
+#include "Mesh.h"
+#include "Node.h"
+#include "Elements/Element.h"
+#include "MeshQuality/MeshValidation.h"
+
+LayeredMeshGenerator::LayeredMeshGenerator()
+: _elevation_epsilon(0.0001)
+{
+}
+
+bool LayeredMeshGenerator::createLayers(MeshLib::Mesh const& mesh, std::vector<std::string> const& raster_paths, double noDataReplacementValue)
+{
+ if (mesh.getDimension() != 2 || !allRastersExist(raster_paths))
+ return false;
+
+ std::vector<GeoLib::Raster const*> rasters;
+ rasters.reserve(raster_paths.size());
+ for (auto path = raster_paths.begin(); path != raster_paths.end(); ++path)
+ rasters.push_back(FileIO::AsciiRasterInterface::getRasterFromASCFile(*path));
+
+ bool result = createRasterLayers(mesh, rasters, noDataReplacementValue);
+ std::for_each(rasters.begin(), rasters.end(), [](GeoLib::Raster const*const raster){ delete raster; });
+ return result;
+}
+
+MeshLib::Mesh* LayeredMeshGenerator::getMesh(std::string const& mesh_name) const
+{
+ if (_nodes.empty() || _elements.empty())
+ return nullptr;
+
+ MeshLib::Mesh* result (new MeshLib::Mesh(mesh_name, _nodes, _elements));
+ MeshLib::MeshValidation::removeUnusedMeshNodes(*result);
+ return result;
+}
+
+double LayeredMeshGenerator::calcEpsilon(GeoLib::Raster const& high, GeoLib::Raster const& low)
+{
+ const double max (*std::max_element(high.begin(), high.end()));
+ const double min (*std::min_element( low.begin(), low.end()));
+ return ((max-min)*1e-07);
+}
+
+bool LayeredMeshGenerator::allRastersExist(std::vector<std::string> const& raster_paths) const
+{
+ for (auto raster = raster_paths.begin(); raster != raster_paths.end(); ++raster)
+ {
+ std::ifstream file_stream (*raster, std::ifstream::in);
+ if (!file_stream.good())
+ return false;
+ file_stream.close();
+ }
+ return true;
+}
+
+void LayeredMeshGenerator::cleanUpOnError()
+{
+ std::for_each(_nodes.begin(), _nodes.end(), [](MeshLib::Node *node) { delete node; });
+ std::for_each(_elements.begin(), _elements.end(), [](MeshLib::Element *elem) { delete elem; });
+}
diff --git a/MeshLib/MeshGenerators/LayeredMeshGenerator.h b/MeshLib/MeshGenerators/LayeredMeshGenerator.h
new file mode 100644
index 0000000000000000000000000000000000000000..9cd3f81599642af2a2099426aac00241ae637fb8
--- /dev/null
+++ b/MeshLib/MeshGenerators/LayeredMeshGenerator.h
@@ -0,0 +1,79 @@
+/**
+ * \file LayeredMeshGenerator.h
+ * \author Karsten Rink
+ * \date 2014-09-18
+ * \brief Definition of the SubsurfaceMapper class
+ *
+ * \copyright
+ * Copyright (c) 2012-2014, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#ifndef LAYEREDMESHGENERATOR_H
+#define LAYEREDMESHGENERATOR_H
+
+#include <string>
+#include <vector>
+
+namespace GeoLib {
+ class Raster;
+}
+
+namespace MeshLib {
+ class Mesh;
+ class Node;
+ class Element;
+}
+
+/**
+ * \brief Base class for creation of 3D subsurface meshes based on raster data
+ */
+class LayeredMeshGenerator
+{
+public:
+ /**
+ * Returns a subsurface representation of a region represented by a 2D mesh by reading raster files and calling the appropriate construction method.
+ * @param mesh The 2D surface mesh that is used as a basis for the subsurface mesh
+ * @param raster_paths Containing all the raster-file-names for the subsurface layers from bottom to top (starting with the bottom of the oldest layer and ending with the DEM)
+ * @param noDataReplacementValue Default z-coordinate if there are mesh nodes not located on the DEM raster (i.e. raster_paths[0])
+ * @result true if the subsurface representation has been created, false if there was an error
+ */
+ virtual bool createLayers(MeshLib::Mesh const& mesh, std::vector<std::string> const& raster_paths, double noDataReplacementValue = 0.0) final;
+
+ /**
+ * Constructs a subsurface representation based on a 2D mesh and a number of rasters representing subsurface layer boundaries.
+ * @param mesh The 2D surface mesh that is used as a basis for the subsurface mesh
+ * @param rasters Containing all the raster-data for the subsurface layers from bottom to top (starting with the bottom of the oldest layer and ending with the DEM)
+ * @param noDataReplacementValue Default z-coordinate if there are mesh nodes not located on the DEM raster (i.e. raster_paths[0])
+ * @result true if the subsurface representation has been created, false if there was an error
+ */
+ virtual bool createRasterLayers(MeshLib::Mesh const& mesh, std::vector<GeoLib::Raster const*> const& rasters, double noDataReplacementValue) = 0;
+
+ /// Returns a mesh of the subsurface representation
+ MeshLib::Mesh* getMesh(std::string const& mesh_name) const;
+
+protected:
+ LayeredMeshGenerator();
+ ~LayeredMeshGenerator() {}
+
+ /// Adds another layer to the subsurface mesh
+ virtual void addLayerToMesh(MeshLib::Mesh const& mesh_layer, unsigned layer_id, GeoLib::Raster const& raster) = 0;
+
+ /// Calculates a data-dependent epsilon value
+ double calcEpsilon(GeoLib::Raster const& high, GeoLib::Raster const& low);
+
+ /// Checks if all raster files actually exist
+ bool allRastersExist(std::vector<std::string> const& raster_paths) const;
+
+ /// Cleans up the already created objects in case of an error
+ void cleanUpOnError();
+
+ double _elevation_epsilon;
+ std::vector<MeshLib::Node*> _nodes;
+ std::vector<MeshLib::Element*> _elements;
+};
+
+#endif //LAYEREDMESHGENERATOR_H
diff --git a/MeshLib/MeshGenerators/LayeredVolume.cpp b/MeshLib/MeshGenerators/LayeredVolume.cpp
index 185df44928111dcbc79a98e9f9544448795a2b7b..3263dcb0343c30bc69accbfb9605e2c0e53c5409 100644
--- a/MeshLib/MeshGenerators/LayeredVolume.cpp
+++ b/MeshLib/MeshGenerators/LayeredVolume.cpp
@@ -14,47 +14,17 @@
#include "LayeredVolume.h"
-#include <fstream>
-#include <numeric>
-
#include "Vector3.h"
-#include "GEOObjects.h"
-#include "PointVec.h"
-#include "Mesh.h"
-#include "convertMeshToGeo.h"
-#include "Elements/Element.h"
+#include "Raster.h"
+
#include "Elements/Tri.h"
#include "Elements/Quad.h"
#include "MeshGenerators/MeshLayerMapper.h"
-#include "MeshQuality/MeshValidation.h"
#include "MeshEditing/ElementExtraction.h"
-const double LayeredVolume::_invalid_value = -9999;
-
-LayeredVolume::LayeredVolume()
-: _elevation_epsilon(0.0001), _mesh(nullptr)
-{
-}
-
-bool LayeredVolume::createGeoVolumes(const MeshLib::Mesh &mesh, const std::vector<std::string> &raster_paths, double noDataReplacementValue)
-{
- if (mesh.getDimension() != 2 || !allRastersExist(raster_paths))
- return false;
-
- std::vector<GeoLib::Raster const*> rasters;
- rasters.reserve(raster_paths.size());
- for (auto path = raster_paths.begin(); path != raster_paths.end(); ++path)
- rasters.push_back(GeoLib::Raster::getRasterFromASCFile(*path));
-
- const bool result = this->createGeoVolumes(mesh, rasters, noDataReplacementValue);
- std::for_each(rasters.begin(), rasters.end(), [](GeoLib::Raster const*const raster){ delete raster; });
- return result;
-}
-
-
-bool LayeredVolume::createGeoVolumes(const MeshLib::Mesh &mesh, const std::vector<GeoLib::Raster const*> &rasters, double noDataReplacementValue)
+bool LayeredVolume::createRasterLayers(const MeshLib::Mesh &mesh, const std::vector<GeoLib::Raster const*> &rasters, double noDataReplacementValue)
{
if (mesh.getDimension() != 2)
return false;
@@ -74,34 +44,33 @@ bool LayeredVolume::createGeoVolumes(const MeshLib::Mesh &mesh, const std::vecto
const std::size_t nRasters (rasters.size());
for (size_t i=0; i<nRasters; ++i)
{
- const double replacement_value = (i==0) ? noDataReplacementValue : _invalid_value;
- if (!MeshLayerMapper::LayerMapping(*mesh_layer, *rasters[i], 0, 0, replacement_value))
+ const double replacement_value = (i==(nRasters-1)) ? noDataReplacementValue : rasters[i]->getNoDataValue();
+ if (!MeshLayerMapper::layerMapping(*mesh_layer, *rasters[i], replacement_value))
{
this->cleanUpOnError();
return false;
}
- this->addLayerToMesh(*mesh_layer, i);
+ this->addLayerToMesh(*mesh_layer, i, *rasters[i]);
}
// close boundaries between layers
this->addLayerBoundaries(*mesh_layer, nRasters);
this->removeCongruentElements(nRasters, mesh_layer->getNElements());
delete mesh_layer;
- _mesh = new MeshLib::Mesh("BoundaryMesh", _nodes, _elements);
- MeshLib::MeshValidation::removeUnusedMeshNodes(*_mesh);
return true;
}
-void LayeredVolume::addLayerToMesh(const MeshLib::Mesh &mesh_layer, unsigned layer_id)
+void LayeredVolume::addLayerToMesh(const MeshLib::Mesh &mesh_layer, unsigned layer_id, GeoLib::Raster const& raster)
{
const std::vector<MeshLib::Node*> &layer_nodes (mesh_layer.getNodes());
const std::size_t nNodes (layer_nodes.size());
const std::size_t node_id_offset (_nodes.size());
const std::size_t last_layer_offset (node_id_offset-nNodes);
+ const double no_data_value (raster.getNoDataValue());
for (std::size_t i=0; i<nNodes; ++i)
{
if (layer_id > 0 &&
- ((*layer_nodes[i])[2] == _invalid_value ||
+ ((*layer_nodes[i])[2] == no_data_value ||
(*_nodes[last_layer_offset+i])[2]-(*layer_nodes[i])[2] < _elevation_epsilon))
_nodes.push_back(new MeshLib::Node(*_nodes[last_layer_offset+i]));
else
@@ -196,36 +165,3 @@ void LayeredVolume::removeCongruentElements(std::size_t nLayers, std::size_t nEl
auto elem_vec_end = std::remove(_elements.begin(), _elements.end(), nullptr);
_elements.erase(elem_vec_end, _elements.end());
}
-
-bool LayeredVolume::exportToGeometry(GeoLib::GEOObjects &geo_objects) const
-{
- if (_mesh == nullptr)
- return false;
- MeshLib::convertMeshToGeo(*_mesh, geo_objects, std::numeric_limits<double>::min());
- return true;
-}
-
-double LayeredVolume::calcEpsilon(const GeoLib::Raster &high, const GeoLib::Raster &low)
-{
- const double max (*std::max_element(high.begin(), high.end()));
- const double min (*std::min_element( low.begin(), low.end()));
- return ((max-min)*1e-07);
-}
-
-bool LayeredVolume::allRastersExist(const std::vector<std::string> &raster_paths) const
-{
- for (auto raster = raster_paths.begin(); raster != raster_paths.end(); ++raster)
- {
- std::ifstream file_stream (*raster, std::ifstream::in);
- if (!file_stream.good())
- return false;
- file_stream.close();
- }
- return true;
-}
-
-void LayeredVolume::cleanUpOnError()
-{
- std::for_each(_nodes.begin(), _nodes.end(), [](MeshLib::Node *node) { delete node; });
- std::for_each(_elements.begin(), _elements.end(), [](MeshLib::Element *elem) { delete elem; });
-}
diff --git a/MeshLib/MeshGenerators/LayeredVolume.h b/MeshLib/MeshGenerators/LayeredVolume.h
index ee5ecb011aca5e65426fbf5870784e44cdd69ff1..9d6a2560df47d3ecc6ef2892f126fb9e10b6d0bd 100644
--- a/MeshLib/MeshGenerators/LayeredVolume.h
+++ b/MeshLib/MeshGenerators/LayeredVolume.h
@@ -18,58 +18,38 @@
#include <string>
#include <vector>
-#include "Raster.h"
#include "Node.h"
+#include "LayeredMeshGenerator.h"
namespace GeoLib {
class GEOObjects;
class Surface;
}
-namespace MeshLib {
- class Mesh;
- class Element;
-}
-
/**
* \brief Creates a volume geometry from 2D mesh layers based on raster data.
*/
-class LayeredVolume
+class LayeredVolume : public LayeredMeshGenerator
{
public:
- LayeredVolume();
+ LayeredVolume() {}
~LayeredVolume() {}
/**
* Constructs a subsurface representation of a mesh using only 2D elements (i.e. layer boundaries are represented by surfaces)
* @param mesh The 2D surface mesh that is used as a basis for the subsurface mesh
- * @param rasters Containing all the raster-data for the subsurface layers from top to bottom (starting with the DEM)
+ * @param rasters Containing all the raster-data for the subsurface layers from bottom to top (starting with the bottom of the oldest layer and ending with the DEM)
* @param noDataReplacementValue Default z-coordinate if there are mesh nodes not located on the DEM raster (i.e. raster_paths[0])
* @result true if the subsurface representation has been created, false if there was an error
*/
- bool createGeoVolumes(const MeshLib::Mesh &mesh, const std::vector<GeoLib::Raster const*> &rasters, double noDataReplacementValue = 0.0);
-
- /**
- * Constructs a subsurface representation of a mesh using only 2D elements (i.e. layer boundaries are represented by surfaces)
- * @param mesh The 2D surface mesh that is used as a basis for the subsurface mesh
- * @param raster_paths Containing all the raster-file-names for the subsurface layers from top to bottom (starting with the DEM)
- * @param noDataReplacementValue Default z-coordinate if there are mesh nodes not located on the DEM raster (i.e. raster_paths[0])
- * @result true if the subsurface representation has been created, false if there was an error
- */
- bool createGeoVolumes(const MeshLib::Mesh &mesh, const std::vector<std::string> &raster_paths, double noDataReplacementValue = 0.0);
-
- /// Returns the subsurface representation that can be used as input for the TetGenInterface
- MeshLib::Mesh* getMesh() const { return _mesh; }
-
+ bool createRasterLayers(const MeshLib::Mesh &mesh, const std::vector<GeoLib::Raster const*> &rasters, double noDataReplacementValue = 0.0);
+
/// Returns the region attribute vector necessary for assigning region attributes via TetGen
std::vector<MeshLib::Node> getAttributePoints() { return _attribute_points; }
- /// Converts the subsurface representation to geometry objects and adds them to the geo storage object
- bool exportToGeometry(GeoLib::GEOObjects &geo_objects) const;
-
private:
/// Adds another layer to the subsurface mesh
- void addLayerToMesh(const MeshLib::Mesh &mesh_layer, unsigned layer_id);
+ void addLayerToMesh(const MeshLib::Mesh &mesh_layer, unsigned layer_id, GeoLib::Raster const& raster);
/// Creates boundary surfaces between the mapped layers to make the volumes watertight
void addLayerBoundaries(const MeshLib::Mesh &layer, std::size_t nLayers);
@@ -77,21 +57,7 @@ private:
/// Removes duplicate 2D elements (possible due to outcroppings)
void removeCongruentElements(std::size_t nLayers, std::size_t nElementsPerLayer);
- /// Calculates a data-dependent epsilon value
- double calcEpsilon(const GeoLib::Raster &high, const GeoLib::Raster &low);
-
- /// Checks if all raster files actually exist
- bool allRastersExist(const std::vector<std::string> &raster_paths) const;
-
- /// Cleans up the already created object in case of an error
- void cleanUpOnError();
-
- static const double _invalid_value;
- double _elevation_epsilon;
- std::vector<MeshLib::Node*> _nodes;
- std::vector<MeshLib::Element*> _elements;
std::vector<MeshLib::Node> _attribute_points;
- MeshLib::Mesh* _mesh;
};
#endif //LAYEREDVOLUME_H
diff --git a/MeshLib/MeshGenerators/MeshLayerMapper.cpp b/MeshLib/MeshGenerators/MeshLayerMapper.cpp
index b302c5e113b04401ba68cda68dafc10061551db2..c56b846d7089d364e21ead757fad6eb349914d5a 100644
--- a/MeshLib/MeshGenerators/MeshLayerMapper.cpp
+++ b/MeshLib/MeshGenerators/MeshLayerMapper.cpp
@@ -12,28 +12,32 @@
*
*/
-// stl
+#include "MeshLayerMapper.h"
+
#include <algorithm>
-// ThirdParty/logog
#include "logog/include/logog.hpp"
-#include "MeshLayerMapper.h"
-// GeoLib
+#include "FileIO/AsciiRasterInterface.h"
+
#include "Raster.h"
-#include "Mesh.h"
-#include "Node.h"
-#include "Elements/Element.h"
#include "Elements/Tet.h"
#include "Elements/Hex.h"
#include "Elements/Pyramid.h"
#include "Elements/Prism.h"
#include "MeshSurfaceExtraction.h"
+
#include "MathTools.h"
+const unsigned MeshLayerMapper::_pyramid_base[3][4] =
+{
+ {1, 3, 4, 2}, // Point 4 missing
+ {2, 4, 3, 0}, // Point 5 missing
+ {0, 3, 4, 1}, // Point 6 missing
+};
-MeshLib::Mesh* MeshLayerMapper::CreateLayers(const MeshLib::Mesh &mesh, const std::vector<float> &layer_thickness_vector)
+MeshLib::Mesh* MeshLayerMapper::createStaticLayers(MeshLib::Mesh const& mesh, std::vector<float> const& layer_thickness_vector, std::string const& mesh_name) const
{
std::vector<float> thickness;
for (std::size_t i=0; i<layer_thickness_vector.size(); ++i)
@@ -45,11 +49,11 @@ MeshLib::Mesh* MeshLayerMapper::CreateLayers(const MeshLib::Mesh &mesh, const st
const std::size_t nLayers(thickness.size());
if (nLayers < 1 || mesh.getDimension() != 2)
{
- ERR("MshLayerMapper::CreateLayers(): A 2D mesh with nLayers > 0 is required as input.");
+ ERR("MeshLayerMapper::createStaticLayers(): A 2D mesh with nLayers > 0 is required as input.");
return nullptr;
}
- const size_t nNodes = mesh.getNNodes();
+ const std::size_t nNodes = mesh.getNNodes();
// count number of 2d elements in the original mesh
const std::size_t nElems (std::count_if(mesh.getElements().begin(), mesh.getElements().end(),
[](MeshLib::Element const* elem) { return (elem->getDimension() == 2);}));
@@ -57,8 +61,7 @@ MeshLib::Mesh* MeshLayerMapper::CreateLayers(const MeshLib::Mesh &mesh, const st
const std::size_t nOrgElems (mesh.getNElements());
const std::vector<MeshLib::Node*> &nodes = mesh.getNodes();
const std::vector<MeshLib::Element*> &elems = mesh.getElements();
- std::vector<MeshLib::Node*> new_nodes;
- new_nodes.reserve(nNodes + (nLayers * nNodes));
+ std::vector<MeshLib::Node*> new_nodes(nNodes + (nLayers * nNodes));
std::vector<MeshLib::Element*> new_elems;
new_elems.reserve(nElems * nLayers);
double z_offset (0.0);
@@ -68,304 +71,182 @@ MeshLib::Mesh* MeshLayerMapper::CreateLayers(const MeshLib::Mesh &mesh, const st
// add nodes for new layer
unsigned node_offset (nNodes * layer_id);
if (layer_id > 0) z_offset += thickness[layer_id-1];
- for (unsigned i = 0; i < nNodes; ++i)
- {
- const double* coords = nodes[i]->getCoords();
- new_nodes.push_back (new MeshLib::Node(coords[0], coords[1], coords[2]-z_offset, node_offset+i));
- }
+
+ std::transform(nodes.cbegin(), nodes.cend(), new_nodes.begin() + node_offset,
+ [&z_offset](MeshLib::Node* node){ return new MeshLib::Node((*node)[0], (*node)[1], (*node)[2]-z_offset); });
// starting with 2nd layer create prism or hex elements connecting the last layer with the current one
- if (layer_id > 0)
+ if (layer_id == 0)
+ continue;
+
+ node_offset -= nNodes;
+ const unsigned mat_id (nLayers - layer_id);
+
+ for (unsigned i = 0; i < nOrgElems; ++i)
{
- node_offset -= nNodes;
- const unsigned mat_id (nLayers - layer_id);
+ const MeshLib::Element* sfc_elem( elems[i] );
+ if (sfc_elem->getDimension() < 2) // ignore line-elements
+ continue;
+
+ const unsigned nElemNodes(sfc_elem->getNNodes());
+ MeshLib::Node** e_nodes = new MeshLib::Node*[2*nElemNodes];
- for (unsigned i = 0; i < nOrgElems; ++i)
+ for (unsigned j=0; j<nElemNodes; ++j)
{
- const MeshLib::Element* sfc_elem( elems[i] );
- if (sfc_elem->getDimension() != 2)
- {
- WARN("MshLayerMapper::CreateLayers() - Method can only handle 2D mesh elements.");
- WARN("Skipping Element %d of type \"%s\".", i, MeshElemType2String(sfc_elem->getGeomType()).c_str());
- continue;
- }
-
- const unsigned nElemNodes(sfc_elem->getNNodes());
- MeshLib::Node** e_nodes = new MeshLib::Node*[2*nElemNodes];
-
- for (unsigned j=0; j<nElemNodes; ++j)
- {
- const unsigned node_id = sfc_elem->getNode(j)->getID() + node_offset;
- e_nodes[j] = new_nodes[node_id+nNodes];
- e_nodes[j+nElemNodes] = new_nodes[node_id];
- }
- if (sfc_elem->getGeomType() == MeshElemType::TRIANGLE) // extrude triangles to prism
- new_elems.push_back (new MeshLib::Prism(e_nodes, mat_id));
- else if (sfc_elem->getGeomType() == MeshElemType::QUAD) // extrude quads to hexes
- new_elems.push_back (new MeshLib::Hex(e_nodes, mat_id));
+ const unsigned node_id = sfc_elem->getNode(j)->getID() + node_offset;
+ e_nodes[j] = new_nodes[node_id+nNodes];
+ e_nodes[j+nElemNodes] = new_nodes[node_id];
}
+ if (sfc_elem->getGeomType() == MeshElemType::TRIANGLE) // extrude triangles to prism
+ new_elems.push_back (new MeshLib::Prism(e_nodes, mat_id));
+ else if (sfc_elem->getGeomType() == MeshElemType::QUAD) // extrude quads to hexes
+ new_elems.push_back (new MeshLib::Hex(e_nodes, mat_id));
}
}
- return new MeshLib::Mesh("SubsurfaceMesh", new_nodes, new_elems);
+ return new MeshLib::Mesh(mesh_name, new_nodes, new_elems);
}
-bool MeshLayerMapper::LayerMapping(MeshLib::Mesh &new_mesh, const std::string &rasterfile,
- const unsigned nLayers, const unsigned layer_id, double noDataReplacementValue = 0.0)
+bool MeshLayerMapper::createRasterLayers(MeshLib::Mesh const& mesh, std::vector<GeoLib::Raster const*> const& rasters, double noDataReplacementValue)
{
- const GeoLib::Raster *raster(GeoLib::Raster::getRasterFromASCFile(rasterfile));
+ const std::size_t nLayers(rasters.size());
+ if (nLayers < 1 || mesh.getDimension() != 2)
+ {
+ ERR("MeshLayerMapper::createRasterLayers(): A 2D mesh and at least two rasters required as input.");
+ return nullptr;
+ }
+
+ MeshLib::Mesh* dem_mesh (new MeshLib::Mesh(mesh));
+ if (layerMapping(*dem_mesh, *rasters.back(), 0))
+ {
+ std::size_t const nNodes = mesh.getNNodes();
+ _nodes.reserve(nLayers * nNodes);
+
+ // number of triangles in the original mesh
+ std::size_t const nElems (std::count_if(mesh.getElements().begin(), mesh.getElements().end(),
+ [](MeshLib::Element const* elem) { return (elem->getGeomType() == MeshElemType::TRIANGLE);}));
+ _elements.reserve(nElems * (nLayers-1));
+
+ for (std::size_t i=0; i<nLayers; ++i)
+ addLayerToMesh(*dem_mesh, i, *rasters[i]);
+
+ return true;
+ }
+ return false;
+}
+
+void MeshLayerMapper::addLayerToMesh(const MeshLib::Mesh &dem_mesh, unsigned layer_id, GeoLib::Raster const& raster)
+{
+ std::size_t const nNodes = dem_mesh.getNNodes();
+ std::vector<MeshLib::Node*> const& nodes = dem_mesh.getNodes();
+ int const last_layer_node_offset = (layer_id-1) * nNodes;
+ double const no_data_value (raster.getNoDataValue());
+
+ // add nodes for new layer
+ for (std::size_t i=0; i<nNodes; ++i)
+ {
+ // min of dem elevation and layer elevation
+ double const elevation = std::min(raster.interpolateValueAtPoint(*nodes[i]), (*nodes[i])[2]);
+
+ if ((layer_id > 0) &&
+ ((std::abs(elevation - no_data_value) < std::numeric_limits<double>::epsilon()) ||
+ (elevation - (*_nodes[last_layer_node_offset + i])[2] < std::numeric_limits<double>::epsilon())))
+ _nodes.push_back(new MeshLib::Node((*nodes[i])[0], (*nodes[i])[1], elevation, _nodes[last_layer_node_offset +i]->getID()));
+ else
+ _nodes.push_back(new MeshLib::Node((*nodes[i])[0], (*nodes[i])[1], elevation, (layer_id * nNodes) + i));
+ }
+
+ if (layer_id == 0)
+ return;
+
+ std::vector<MeshLib::Element*> const& elems = dem_mesh.getElements();
+ std::size_t const nElems (dem_mesh.getNElements());
+
+ for (std::size_t i=0; i<nElems; ++i)
+ {
+ MeshLib::Element* elem (elems[i]);
+ if (elem->getGeomType() != MeshElemType::TRIANGLE)
+ continue;
+ unsigned node_counter(3), missing_idx(0);
+ std::array<MeshLib::Node*, 6> new_elem_nodes;
+ for (unsigned j=0; j<3; ++j)
+ {
+ new_elem_nodes[j] = _nodes[last_layer_node_offset + elem->getNodeIndex(j)];
+ new_elem_nodes[node_counter] = (_nodes[last_layer_node_offset + elem->getNodeIndex(j) + nNodes]);
+ if (new_elem_nodes[j]->getID() != new_elem_nodes[node_counter]->getID())
+ node_counter++;
+ else
+ missing_idx = j;
+ }
+
+ switch (node_counter)
+ {
+ case 6:
+ _elements.push_back(new MeshLib::Prism(new_elem_nodes, layer_id));
+ break;
+ case 5:
+ std::array<MeshLib::Node*, 5> pyramid_nodes;
+ pyramid_nodes[0] = new_elem_nodes[_pyramid_base[missing_idx][0]];
+ pyramid_nodes[1] = new_elem_nodes[_pyramid_base[missing_idx][1]];
+ pyramid_nodes[2] = new_elem_nodes[_pyramid_base[missing_idx][2]];
+ pyramid_nodes[3] = new_elem_nodes[_pyramid_base[missing_idx][3]];
+ pyramid_nodes[4] = new_elem_nodes[missing_idx];
+ _elements.push_back(new MeshLib::Pyramid(pyramid_nodes, layer_id));
+ break;
+ case 4:
+ std::array<MeshLib::Node*, 4> tet_nodes;
+ std::copy(new_elem_nodes.begin(), new_elem_nodes.begin() + node_counter, tet_nodes.begin());
+ _elements.push_back(new MeshLib::Tet(tet_nodes, layer_id));
+ break;
+ default:
+ continue;
+ }
+ }
+}
+
+bool MeshLayerMapper::layerMapping(MeshLib::Mesh &new_mesh, std::string const& rasterfile, double noDataReplacementValue = 0.0)
+{
+ const GeoLib::Raster *raster(FileIO::AsciiRasterInterface::getRasterFromASCFile(rasterfile));
if (! raster) {
- ERR("MshLayerMapper::LayerMapping - could not read raster file %s", rasterfile.c_str());
+ ERR("MshLayerMapper::layerMapping - could not read raster file %s", rasterfile.c_str());
return false;
}
- const bool result = LayerMapping(new_mesh, *raster, nLayers, layer_id, noDataReplacementValue);
+ const bool result = layerMapping(new_mesh, *raster, noDataReplacementValue);
delete raster;
return result;
}
-bool MeshLayerMapper::LayerMapping(MeshLib::Mesh &new_mesh, const GeoLib::Raster &raster,
- const unsigned nLayers, const unsigned layer_id, double noDataReplacementValue = 0.0)
+bool MeshLayerMapper::layerMapping(MeshLib::Mesh &new_mesh, GeoLib::Raster const& raster, double noDataReplacementValue = 0.0)
{
- if (nLayers < layer_id)
+ if (new_mesh.getDimension() != 2)
{
- ERR("MshLayerMapper::LayerMapping() - Mesh has only %d Layers, cannot assign layer %d.", nLayers, layer_id);
+ ERR("MshLayerMapper::layerMapping - requires 2D mesh");
return false;
}
const double x0(raster.getOrigin()[0]);
const double y0(raster.getOrigin()[1]);
- const double delta(raster.getRasterPixelDistance());
+ 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));
-
- 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)
+ for (unsigned i = 0; i < nNodes; ++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 (std::abs(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)
-{
- 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)
-{
- // construct surface mesh from DEM
- const MathLib::Vector3 dir(0,0,1);
- MeshLib::Mesh* dem = MeshLib::MeshSurfaceExtraction::getMeshSurface(mesh, dir);
- MeshLayerMapper::LayerMapping(*dem, dem_raster, 0, 0);
- const std::vector<MeshLib::Node*> &dem_nodes (dem->getNodes());
-
- const std::vector<MeshLib::Node*> &mdl_nodes (mesh.getNodes());
- const size_t nNodes (mesh.getNNodes());
- const size_t nNodesPerLayer = nNodes / (nLayers+1);
- std::vector<bool> is_surface_node(nNodes, false);
- std::vector<bool> nodes_below_surface(nNodes, false);
-
- // check if bottom layer nodes are below DEM
- const unsigned bottom_firstNode = nLayers * nNodesPerLayer;
- const unsigned bottom_lastNode = bottom_firstNode + nNodesPerLayer;
- for (unsigned i = bottom_firstNode; i < bottom_lastNode; ++i)
- {
- nodes_below_surface[i]=true;
- const double* coords = mdl_nodes[i]->getCoords();
- const double* dem_coords = dem_nodes[i-bottom_firstNode]->getCoords();
- if (coords[2] >= dem_coords[2])
- {
- WARN("Node %d (in bottom-layer) is above surface node %d. (%f, %f)", i, (i-bottom_firstNode), coords[2], dem_coords[2]);
- is_surface_node[i] = true;
- }
- }
-
- // for all other layers:
- // if node < dem-node: do nothing
- // if node > dem-node:
- // if first node above surface: map to dem and mark as surface node
- // else remove node (i.e. don't copy it)
- for (int layer_id=nLayers-1; layer_id>=0; --layer_id)
- {
- const size_t firstNode = layer_id * nNodesPerLayer;
- const size_t lastNode = firstNode + nNodesPerLayer;
-
- for(unsigned i = firstNode; i < lastNode; ++i)
- {
- if (is_surface_node[i+nNodesPerLayer])
- is_surface_node[i]=true;
- else
- {
- nodes_below_surface[i]=true;
- MeshLib::Node* node (mdl_nodes[i]);
- const double* coords = node->getCoords();
- const double* dem_coords = dem_nodes[i-firstNode]->getCoords();
- if (coords[2] > dem_coords[2])
- {
- node->updateCoordinates(dem_coords[0], dem_coords[1], dem_coords[2]);
- is_surface_node[i] = true;
- }
- }
- }
- }
- delete dem; // no longer needed
-
- // copy valid nodes to new node vector
- std::vector<MeshLib::Node*> new_nodes;
- std::vector<int> node_index_map(nNodes, -1);
- size_t node_count(0);
- for (unsigned j=0; j<nNodes; ++j)
- if (nodes_below_surface[j])
- {
- new_nodes.push_back(new MeshLib::Node(mdl_nodes[j]->getCoords(), mdl_nodes[j]->getID()));
- node_index_map[j]=node_count++;
- }
-
- // copy elements (old elements need to have at least 4 nodes remaining and form a 3d element
- const std::vector<MeshLib::Element*> &mdl_elements (mesh.getElements());
- const size_t nElems (mesh.getNElements());
- std::vector<MeshLib::Element*> new_elements;
- for (unsigned j=0; j<nElems; ++j)
- {
- const MeshLib::Element* elem = mdl_elements[j];
-
- size_t count(0);
- for (unsigned i=0; i<6; ++i) // check top surface of prism
- if (nodes_below_surface[elem->getNode(i)->getID()]) ++count;
-
- if (count==6) // copy prism elements if all six nodes are valid
- {
- MeshLib::Node** e_nodes = new MeshLib::Node*[count];
- for (unsigned i=0; i<6; ++i)
- e_nodes[i] = new_nodes[node_index_map[elem->getNode(i)->getID()]];
-
- MeshLib::Element* prism (new MeshLib::Prism(e_nodes, elem->getValue()));
- new_elements.push_back(prism);
- }
- else if (count==5) // change the current element to two tetrahedra if only five nodes are valid
- {
- MeshLib::Node** e_nodes = new MeshLib::Node*[count];
- unsigned top_idx(6);
- for (unsigned i=3; i<6; ++i) // find node that has been cut
- if (!nodes_below_surface[elem->getNode(i)->getID()])
- top_idx = i-3;
-
- // construct pyramid element based on missing node
- unsigned idx1 ((top_idx+1)%3);
- unsigned idx2 ((top_idx+2)%3);
- e_nodes[0] = new_nodes[node_index_map[elem->getNode(idx1)->getID()]];
- e_nodes[1] = new_nodes[node_index_map[elem->getNode(idx1+3)->getID()]];
- e_nodes[2] = new_nodes[node_index_map[elem->getNode(idx2+3)->getID()]];
- e_nodes[3] = new_nodes[node_index_map[elem->getNode(idx2)->getID()]];
- e_nodes[4] = new_nodes[node_index_map[elem->getNode(top_idx)->getID()]];
-
- MeshLib::Element* pyr (new MeshLib::Pyramid(e_nodes, elem->getValue()));
- new_elements.push_back(pyr);
- }
- else if (count==4) // change the current element to a tetrahedron if only four nodes are valid
- {
- MeshLib::Node** e_nodes = new MeshLib::Node*[count];
- for (unsigned i=0; i<3; ++i) // first three nodes are the bottom-face
- {
- unsigned idx (elem->getNode(i)->getID());
- if (nodes_below_surface[idx])
- e_nodes[i] = new_nodes[node_index_map[idx]];
- else
- e_nodes[i] = NULL;
- }
-
- if (e_nodes[0] && e_nodes[1] && e_nodes[2]) //make sure that the 4 remaining nodes don't form a quad
- {
- for (unsigned i=3; i<6; ++i) // last node
- {
- unsigned idx (elem->getNode(i)->getID());
- if (nodes_below_surface[idx])
- {
- e_nodes[3] = new_nodes[node_index_map[idx]];
- break;
- }
- }
-
- MeshLib::Element* tet (new MeshLib::Tet(e_nodes, elem->getValue()));
- new_elements.push_back(tet);
- }
- else delete e_nodes;
- }
- // else remove element, if less than four nodes are valid
- }
- return new MeshLib::Mesh("SubsurfaceMesh", new_nodes, new_elements);
-}
-
-
-
diff --git a/MeshLib/MeshGenerators/MeshLayerMapper.h b/MeshLib/MeshGenerators/MeshLayerMapper.h
index 33853f4045f3535b4dfa57fd02b7a0c4d6c970f7..cb17e764adb665f6439d3bc6a7a409b9c0775865 100644
--- a/MeshLib/MeshGenerators/MeshLayerMapper.h
+++ b/MeshLib/MeshGenerators/MeshLayerMapper.h
@@ -15,61 +15,55 @@
#ifndef MESHLAYERMAPPER_H
#define MESHLAYERMAPPER_H
-#include <string>
-#include "GeoLib/Raster.h"
-
-class QImage;
-
-namespace MeshLib {
- class Mesh;
- class Node;
-}
+#include "LayeredMeshGenerator.h"
/**
- * \brief Manipulating and adding layers to an existing mesh
+ * \brief Manipulating and adding prism element layers to an existing 2D mesh
*/
-class MeshLayerMapper
+class MeshLayerMapper : public LayeredMeshGenerator
{
public:
- MeshLayerMapper() {}
- ~MeshLayerMapper() {}
+ MeshLayerMapper() {}
+ ~MeshLayerMapper() {}
- /**
- * Based on a triangle-or quad mesh this method creates a 3D mesh with with a given number of prism- or hex-layers
- * \param mesh The triangle/quad mesh that is the basis for the new prism/hex mesh
- * \param nLayers The number of layers of prism/hex elements that will be extruded from the triangle/quad elements of the original mesh
- * \param thickness The thickness of each of these newly added layers
- * \return A mesh with the requested number of layers of prism/hex elements
- */
- static MeshLib::Mesh* CreateLayers(const MeshLib::Mesh &mesh, const std::vector<float> &layer_thickness_vector);
+ /**
+ * Based on a 2D triangle-or quad mesh this method creates a 3D mesh with a given number of prism- or hex-layers
+ * \param mesh The triangle/quad mesh that is the basis for the new prism/hex mesh
+ * \param nLayers The number of layers of prism/hex elements that will be extruded from the triangle/quad elements of the original mesh
+ * \param thickness The thickness of each of these newly added layers
+ * \param mesh_name The name of the newly created mesh
+ * \return A mesh with the requested number of layers of prism/hex elements
+ */
+ MeshLib::Mesh* createStaticLayers(MeshLib::Mesh const& mesh, std::vector<float> const& layer_thickness_vector, std::string const& mesh_name = "SubsurfaceMesh") const;
- /**
- * 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 bool LayerMapping(MeshLib::Mesh &mesh, const std::string &rasterfile,
- const unsigned nLayers, const unsigned layer_id, double noDataReplacementValue);
+ /**
+ * Based on a 2D triangle mesh this method creates a 3D mesh with a given number of prism-layers.
+ * Note: While this method would technically also work with quad meshes, this is discouraged as quad elements will most likely not
+ * be coplanar after the mapping process which result in invaled mesh elements.
+ * \param mesh The 2D triangle mesh that is the basis for the new 3D prism mesh
+ * \param rasters Containing all the raster-data for the subsurface layers from bottom to top (starting with the bottom of the oldest layer and ending with the DEM)
+ * \param noDataReplacementValue Default z-coordinate if there are mesh nodes not located on the DEM raster (i.e. raster_paths[0])
+ * \return A mesh with the requested number of layers of prism elements (also including Tet- & Pyramid-elements in case of degenerated prisms)
+ */
+ bool createRasterLayers(MeshLib::Mesh const& mesh, std::vector<GeoLib::Raster const*> const& rasters, double noDataReplacementValue = 0.0);
- /**
- * 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 bool LayerMapping(MeshLib::Mesh &mesh, const GeoLib::Raster &raster,
- const unsigned nLayers, const unsigned layer_id, double noDataReplacementValue);
+ /**
+ * Maps the elevation of nodes of a given 2D mesh according to the raster specified by the file path.
+ * At locations wher no information is given, node elevation is set to noDataReplacementValue.
+ */
+ static bool layerMapping(MeshLib::Mesh &mesh, const std::string &rasterfile, double noDataReplacementValue);
- /**
- * Blends a mesh with the surface given by dem_raster. Nodes and elements above the surface are either removed or adapted to fit the surface.
- * Note: It is unlikely but possible that the new nodes vector contains (very few) nodes that are not part of any element. This problem is
- * remedied at the end of method upon creating the actual mesh from the new node- and element-vector as the mesh-constructor checks for such
- * nodes and removes them. This note is just to call this issue to attention in case this methods is changed.
- */
- static MeshLib::Mesh* blendLayersWithSurface(MeshLib::Mesh &mesh, const unsigned nLayers, const std::string &dem_raster);
+ /**
+ * Maps the elevation of nodes of a given 2D mesh according to the raster. At locations wher no
+ * information is given, node elevation is set to noDataReplacementValue.
+ */
+ static bool layerMapping(MeshLib::Mesh &mesh, const GeoLib::Raster &raster, double noDataReplacementValue);
private:
- /// Checks if the given mesh is within the dimensions given by xDim and yDim.
- static bool isNodeOnRaster(const MeshLib::Node &node,
- const std::pair<double, double> &xDim,
- const std::pair<double, double> &yDim);
+ /// Adds another layer to a subsurface mesh
+ void addLayerToMesh(const MeshLib::Mesh &mesh_layer, unsigned layer_id, GeoLib::Raster const& raster);
+
+ static const unsigned _pyramid_base[3][4];
};
#endif //MESHLAYERMAPPER_H
diff --git a/MeshLib/Node.h b/MeshLib/Node.h
index ae0a21165a679a906c5e2e8f06f8757d5f13bc19..557d26bbdccc7614eebfc80fd28b0e4243c5829f 100644
--- a/MeshLib/Node.h
+++ b/MeshLib/Node.h
@@ -36,9 +36,7 @@ class Element;
class Node : public GeoLib::PointWithID
{
/* friend functions: */
- friend bool MeshLayerMapper::LayerMapping(MeshLib::Mesh &mesh, const GeoLib::Raster &raster, const unsigned nLayers,
- const unsigned layer_id, double noDataReplacementValue);
- friend MeshLib::Mesh* MeshLayerMapper::blendLayersWithSurface(MeshLib::Mesh &mesh, const unsigned nLayers, const std::string &dem_raster);
+ friend bool MeshLayerMapper::layerMapping(MeshLib::Mesh &mesh, const GeoLib::Raster &raster, double noDataReplacementValue);
friend MeshLib::Mesh* MeshSurfaceExtraction::getMeshSurface(const MeshLib::Mesh &mesh, const MathLib::Vector3 &dir, bool keep3dMeshIds);
/* friend classes: */
diff --git a/Utils/SimpleMeshCreation/createMeshElemPropertiesFromASCRaster.cpp b/Utils/SimpleMeshCreation/createMeshElemPropertiesFromASCRaster.cpp
index 557d506b6cb6ee5b86ac08d1e07174c036db6412..22b24ad0e1cf50a74eed133f236fe9d5937a9ea6 100644
--- a/Utils/SimpleMeshCreation/createMeshElemPropertiesFromASCRaster.cpp
+++ b/Utils/SimpleMeshCreation/createMeshElemPropertiesFromASCRaster.cpp
@@ -29,6 +29,7 @@
#include "FileTools.h"
#include "MeshIO.h"
#include "readMeshFromFile.h"
+#include "AsciiRasterInterface.h"
// GeoLib
#include "Raster.h"
@@ -122,9 +123,7 @@ int main (int argc, char* argv[])
raster_arg.getValue()));
new_raster_fname += "-" + std::to_string(raster->getNRows()) + "x" +
std::to_string(raster->getNCols()) + ".asc";
- std::ofstream out(new_raster_fname);
- raster->writeRasterAsASC(out);
- out.close();
+ FileIO::AsciiRasterInterface::writeRasterAsASC(raster, new_raster_fname);
}
}